Off-Grid DOA Estimation Using Alternating Block Coordinate Descent in Compressed Sensing

Abstract

This paper presents a novel off-grid direction of arrival (DOA) estimation method to achieve the superior performance in compressed sensing (CS), in which DOA estimation problem is cast as a sparse reconstruction. By minimizing the mixed k-l norm, the proposed method can reconstruct the sparse source and estimate grid error caused by mismatch. An iterative process that minimizes the mixed k-l norm alternately over two sparse vectors is employed so that the nonconvex problem is solved by alternating convex optimization. In order to yield the better reconstruction properties, the block sparse source is exploited for off-grid DOA estimation. A block selection criterion is engaged to reduce the computational complexity. In addition, the proposed method is proved to have the global convergence. Simulation results show that the proposed method has the superior performance in comparisons to existing methods.

1. Introduction

Direction of arrival (DOA) estimation has played an important role in many fields, such as radar, medical imaging and array signal processing [1,2]. In the last decades, many classical methods have been developed, among which multiple signal classification (MUSIC) [3] and estimation of signal parameter via rotational invariance technique. (ESPRIT) [4] are the most popular and have high resolution for DOA estimation. However, these methods are very sensitive to the number of snapshots, signal to noise ratio (SNR) and the correlation between sources. Small number of snapshots, low SNR and high correlation or coherent sources can all make the performance of these methods degrade severely. More recently, the emerging field of compressed sensing (CS) [5,6] has attracted enormous attention and it can reconstruct the sparse source using nonadaptive linear projection measurement obtained by the measurement matrix that satisfies the restricted isometry property (RIP) [7,8]. Support denotes the set that contains the indices of the nonzero elements in the sparse source. Once the support is determined, the sparse source can be reconstructed.

Due to the fact that sources are intrinsically sparse in the spatial domain, the DOA estimation problem can be regarded as a sparse reconstruction in the framework of CS. The CS-based estimation methods have much better estimation performance than conventional estimation methods. Malioutov et al. [9] firstly adopted the sparse signal reconstruction (SSR) perspective for DOA estimation and utilized the singular value decomposition (SVD) of the data matrix to propose $l_1$-SVD method. In [10], CS-MUSIC was proposed by revisiting the link between CS and MUSIC. This method identifies the parts of support using CS, after which the remaining parts are estimated by a novel generalized MUSIC criterion. Xu et al. [11] utilized the Capon spectrum to design a weighted $l_1$-norm penalty in order to further enforce the sparsity and approximate the original $l_0$-norm for DOA estimation. Wei et al. [12] proposed a modified greedy block coordinate descent (R-GBCD) method and the corresponding version with weight refinement (R-GBCD+) to improve the estimation performance.

The key to guarantee the performance of conventional CS-based estimation methods is that all true DOAs are exactly located on the grid. However, when true DOAs are not on the grid set, the performance may severely degrade due to the grid error caused by mismatch, which is defined as the distance from the true direction to the nearest grid. In order to address this issue, Zhu et al. [13] proposed the sparse total least square (STIS) to perform the off-grid DOA estimation, in which perturbations of the model are assume to be Gaussian. In [14], Yang et al. introduced the Bayesian theory in off-grid DOA estimation and proposed an off-grid sparse Bayesian inference based on the singular value decomposition (OGSBI-SVD). Liang et al. [15] proposed an off-grid synchronous approach based on distributed compressed sensing to obtain larger array aperture. Zhang et al. [16] formulated a novel model based on the sampling covariance matrix and solved the off-grid DOA estimation problem by the block sparse Bayesian method even if the number of sources are unknown.

In this paper, a novel alternating block coordinate descent method called ABCD is proposed for off-grid DOA estimation in CS. The proposed method solves the mixed k-l norm minimization problem to reconstruct the sparse source and estimate the grid error. Since joint estimation will lead to a nonconvex optimization problem, the proposed method adopts an iterative process that minimizes the mixed k-l norm alternately over two sparse vectors. Instead of conventional sparse source, the block sparse source is exploited to achieve better reconstruction properties. The block is updated by the proposed block selection criterion, which can improve efficiency of the proposed method. In addition, we give a detailed derivation process of proving the global convergence of the proposed method. Simulation results illustrate the superior performance of the proposed method as compared with existing methods.

The rest of the paper is organized as follows. An off-grid DOA estimation model is formulated in Section 2 and Section 3 introduces the proposed method in detail. The global convergence of the proposed method is proved in Section 4 and Section 5 shows the performance of the proposed method. Conclusions are provided in Section 6.

2. Problem Formulation

Consider K far-field narrowband sources $s_1(t),s_2(t),\cdots ,s_K(t)$ impinging on the uniform linear array (ULA) consisting of M omnidirectional sensors with inter-sensor spacing $\stackrel{⌢}{d}$. Assume that each source $s_k(t)$ is located at different direction ${\mathrm{\theta }}_k$ with the power ${\text{σ}}_k$, $k=1,2,\cdots ,K$. At time instant t, the received source by ULA can be expressed as

$$x(t)={\displaystyle \sum _{k=1}^{K}a({\mathrm{\theta }}_k)s_k(t)}+n(t) t=1,2\cdots T$$

(1)

where $a({\mathrm{\theta }}_k)\in {ℂ}^{M\times 1}$ and $n(t)\in {ℂ}^{M\times 1}$ denote the steering vector and noise vector, respectively. Since the first point of ULA is set as the origin of sensor array, the mth element of $a({\mathrm{\theta }}_k)$ is written as $e^{-j(\text{2π}/\text{λ})\widehat{d}(m-1)\cos {\mathrm{\theta }}_k}$ with the wavelength of source $\text{λ}$.

To incorporate the CS theory with the DOA estimation, the entire angular space is divided into a fine grid $\tilde{\mathrm{\theta }}={[{\tilde{\mathrm{\theta }}}_1,{\tilde{\mathrm{\theta }}}_2\cdots {\tilde{\mathrm{\theta }}}_L]}^T$, where L$(L\gg K)$ denotes the number of the grid and ${[·]}^T$ denotes the transpose. Due to the fact that true directions $\mathrm{\theta }={[{\mathrm{\theta }}_1,{\mathrm{\theta }}_2\cdots {\mathrm{\theta }}_K]}^T$ are random in the entire angular space, ${\mathrm{\theta }}_k$ for some $k\in \{1,2\cdots K\}$ are likely to be not on the grid set. To reduce the grid error caused by mismatch, we formulate the off-grid model, which has a close relationship with the on-grid model. Let $\tilde{\mathrm{\theta }}$ satisfy the uniform distribution so that the grid interval $\text{τ}={\tilde{\mathrm{\theta }}}_{k+1}-{\tilde{\mathrm{\theta }}}_k\propto L^{-1}$. Without loss of generality, by assuming ${\mathrm{\theta }}_k\notin \{{\tilde{\mathrm{\theta }}}_1,{\tilde{\mathrm{\theta }}}_2\cdots {\tilde{\mathrm{\theta }}}_L\}$ and that ${\tilde{\mathrm{\theta }}}_{l_k}\in \{{\tilde{\mathrm{\theta }}}_1,{\tilde{\mathrm{\theta }}}_2\cdots {\tilde{\mathrm{\theta }}}_L\}$ is the nearest grid to ${\mathrm{\theta }}_k$, the steering vector $a({\mathrm{\theta }}_k)$ can be approximated as

$$a({\mathrm{\theta }}_k)\approx a({\tilde{\mathrm{\theta }}}_{l_k})+b({\tilde{\mathrm{\theta }}}_{l_k})({\mathrm{\theta }}_k-{\tilde{\mathrm{\theta }}}_{l_k})$$

(2)

where ${\text{ε}}_l={\mathrm{\theta }}_k-{\tilde{\mathrm{\theta }}}_{l_k}\in [-\frac{\text{τ}}{\text{2}},\frac{\text{τ}}{\text{2}}]$ is the grid error and $b({\tilde{\mathrm{\theta }}}_{l_k})$ is the partial derivative of $a({\tilde{\mathrm{\theta }}}_{l_k})$ with respective to ${\tilde{\mathrm{\theta }}}_{l_k}$. Then, if ${\tilde{\mathrm{\theta }}}_{l_1},{\tilde{\mathrm{\theta }}}_{l_2}\cdots {\tilde{\mathrm{\theta }}}_{l_K}$ are respectively nearest grids to the true directions ${\mathrm{\theta }}_1,{\mathrm{\theta }}_2\cdots {\mathrm{\theta }}_K$, we have

$$\{\begin{array}{ccc}{g}_l(t)=s_{l_k}(t),{\text{ε}}_l={\mathrm{\theta }}_k-{\tilde{\mathrm{\theta }}}_{l_k}& ;& l=l_k(k=1,2,\cdots K)\hfill \\ g_l(t)={\text{ε}}_l=0& ;& elsewhere\hfill \end{array}$$

(3)

Then, by imposing the approximation error on the noise, the received source $x(t)$ can be rewritten as the following sparse form

$$x(t)=[A+Bdiag(\text{ε})]g(t)+\stackrel{⌢}{n}(t) t=1,2,\cdots T$$

(4)

where $A=[a({\tilde{\mathrm{\theta }}}_1),a({\tilde{\mathrm{\theta }}}_2),\cdots ,a({\tilde{\mathrm{\theta }}}_L)]$ is the $M\times L$ array manifold matrix corresponding to all potential directions, which is defined as an overcomplete dictionary in CS, and $B=[b({\tilde{\mathrm{\theta }}}_1),b({\tilde{\mathrm{\theta }}}_2),\cdots ,b({\tilde{\mathrm{\theta }}}_L)]$. In addition, the matrix $\stackrel{⌢}{N}=[\stackrel{⌢}{n}(1),\stackrel{⌢}{n}(2),\cdots ,\stackrel{⌢}{n}(T)]$ is the noise matrix and $diag(\text{ε})$ is a diagonal matrix with $\text{ε}$ being the diagonal elements. Since $g(t)={[g_1(t),g_2(t),\cdots ,g_L(t)]}^T$ has K nonzero elements in L elements, it is a K-sparse vector, where K is referred to as sparsity. More specifically, $\text{ε}={[{\text{ε}}_1,{\text{ε}}_2,\cdots ,{\text{ε}}_L]}^T$ is also a K-sparse vector and has the same support as $g(t)$. It is evident that the on-grid sparse model is a special case by simply setting $\text{ε}=\text{0}$ in Equation (4). Since ${\{g(t)\}}_{t=1}^T$ are jointly K-sparse, the matrix $G=[g(1),g(2),\cdots ,g(T)]\in {ℂ}^{L\times T}$ has K nonzero rows and is called row K-sparse.

To solve the off-grid DOA estimation problem, we need to jointly estimate the support of sparse sources, $g(t),t=1,2,\cdots ,T$, and grid error $\text{ε}$ from the matrix Y which is given by

$$Y=[y(1) y(2)\cdots y(T)]=\text{Φ}AG+\text{Φ}{\rm B}diag(\text{ε})G+\tilde{N}$$

(5)

where $\tilde{N}=\text{Φ}\stackrel{⌢}{N}$ is the $N\times T$ noise matrix, $\text{Φ}$ is the $N\times M$ measurement matrix with $N<M$ and N is the number of nonadaptive linear projection measurement.

3. Off-DOA Estimation

In this section, an alternating block coordinate descent (ABCD) method and a block selection criterion are elaborated in the CS scenario. The proposed method not only has the advantages of conventional BCD [17] strategy, but also uses an iterative process that minimizes the mixed k-l norm alternately over two sparse vectors. Note that due to solving the minimization alternately, a tractable convex problem is obtained and the global convergence of ABCD can be easily determined, which is proved in the next section.

By applying the central limit theorem, the components $\stackrel{⌢}{n}(t)$, $t=1,2,\cdots ,T$, of $\stackrel{⌢}{N}$ are independently white Gaussian noise with zero mean and covariance ${\text{σ}}^2{I}_M$, where $I_M$ denotes an $M\times M$ identity matrix. Thus, the covariance matrix of $Y$ with the size $M\times M$ is expressed as

$$R_Y=E[y(t)y^H(t)]=\text{Φ}A{R}_G{A}^H{\text{Φ}}^H+\text{Φ}{\rm B}diag(\text{ε})R_{G}dia{g}^H(\text{ε}){{\rm B}}^H{\text{Φ}}^H+{\text{σ}}^2{\text{ΦΦ}}^H$$

(6)

where $R_G=E[s(t)s^H(t)]$ is a $L\times L$ covariance matrix of the sparse source and ${(·)}^H$ denotes the conjugate transpose. Since all potential directions are one to one corresponding to the powers ${\text{σ}}_1^2,{\text{σ}}_2^2,\cdots ,{\text{σ}}_L^2$ and we are interested in estimating DOAs, $R_G$ can be reduced to a diagonal matrix $R_G=diag(p)$, where $p={[{\text{σ}}_1^2,{\text{σ}}_2^2,\cdots ,{\text{σ}}_L^2]}^T$ is a K-sparse vector. Then, by denoting $\widehat{{\rm A}}=\text{Φ}{\rm A}=[\widehat{a}({\tilde{\mathrm{\theta }}}_1),\widehat{a}({\tilde{\mathrm{\theta }}}_2),\cdots ,\widehat{a}({\tilde{\theta }}_L)]$ and $\widehat{{\rm B}}=\text{Φ}{\rm B}=[\widehat{b}({\tilde{\mathrm{\theta }}}_1),\widehat{b}({\tilde{\mathrm{\theta }}}_2),\cdots ,\widehat{b}({\tilde{\mathrm{\theta }}}_L)]$, (6) can be further rewritten as

$$R_Y={\displaystyle \sum _{l=1}^L\widehat{a}({\tilde{\mathrm{\theta }}}_l){\widehat{a}}^H({\tilde{\mathrm{\theta }}}_l){\text{σ}}_l^2}+{\displaystyle \sum _{l=1}^L\widehat{b}({\tilde{\mathrm{\theta }}}_l){\widehat{b}}^H({\tilde{\mathrm{\theta }}}_l){\text{σ}}_l^2{\text{ε}}_l^2}+{\text{σ}}^2{\text{ΦΦ}}^H$$

(7)

Due to the vector form of $R_Y$ in Equation (7), the following measurement vector is given by

$$z=vec(R_Y)=Cp+Dq+{\text{σ}}^2{\text{ΦΦ}}^H$$

(8)

with $C=[{\widehat{a}}^H({\tilde{\mathrm{\theta }}}_1)\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_1), {\widehat{a}}^H({\tilde{\mathrm{\theta }}}_2)\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_2),\cdots , {\widehat{a}}^H({\tilde{\mathrm{\theta }}}_L)\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_L)]$ and $D=[{\widehat{b}}^H({\tilde{\mathrm{\theta }}}_1)\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_1), {\widehat{b}}^H({\tilde{\mathrm{\theta }}}_2)\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_2),\cdots ,{\widehat{b}}^H({\tilde{\mathrm{\theta }}}_L)\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_L)]$, where $\otimes$ and $vec(·)$ denote the Kronecker product and the stack operation by placing the columns of a matrix on the top of one another in order, respectively. Moreover, $q={[{\text{σ}}_1^2{\text{ε}}_1^2,{\text{σ}}_2^2{\text{ε}}_2^2,\cdots ,{\text{σ}}_L^2{\text{ε}}_L^2]}^T$ is also a K-sparse vector, which has the same support as $p$. In the conventional sparse source, the nonzero elements of the sparse vector $p$ or $q$ can appear anywhere in the vector. However, in this paper, our goal is to explicitly take block sparse source into account, i.e., the nonzero elements of the sparse vector $p$ or $q$ tend to cluster in blocks. The motivations to exploit block sparse source are the following two main reasons. As can be seen in [18], the first reason is that block sparse source has been applied in many applications, such as unions of subspaces and multiband sources [19]. Secondly, block sparse source has better reconstruction properties than sparse source in the conventional sense, which is proved in [20]. To exploit block sparse source, denote $p[h]={[{\text{σ}}_{d(h-1)+1}^2,{\text{σ}}_{d(h-1)+2}^2,\cdots ,{\text{σ}}_{dh}^2]}^T$ and $q[h]={[{\text{σ}}_{d(h-1)+1}^2{\text{ε}}_{d(h-1)+1}^2,{\text{σ}}_{d(h-1)+2}^2{\text{ε}}_{d(h-1)+2}^2,\cdots ,{\text{σ}}_{dh}^2{\text{ε}}_{dh}^2]}^T$ as the hth blocks of $p$ and $q$ with the length d, respectively, so that we have

$$p={[p^T[1],p^T[2],\cdots ,p^T[H]]}^T$$

(9)

$$q={[q^T[1],q^T[2],\cdots ,q^T[H]]}^T$$

(10)

where $L=dH$. Following the similar manner, the matrices $C$ and $D$ can be respectively viewed as a concatenation of block matrices $C[h]$ and $D[h]$ of the size $N^2\times d$, i.e.,

$$C=[C[1],C[2],\cdots ,C[H]]$$

(11)

$$D=[D[1],D[2],\cdots ,D[H]]$$

(12)

where $C[h]=[{\widehat{a}}^H({\tilde{\mathrm{\theta }}}_{d(h-1)+1})\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_{d(h-1)+1}),{\widehat{a}}^H({\tilde{\mathrm{\theta }}}_{d(h-1)+2})\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_{d(h-1)+2}),\cdots ,{\widehat{a}}^H({\tilde{\mathrm{\theta }}}_{dh})\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_{dh})]$ and $D[h]=[{\widehat{b}}^H({\tilde{\mathrm{\theta }}}_{d(h-1)+1})\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_{d(h-1)+1}),{\widehat{b}}^H({\tilde{\mathrm{\theta }}}_{d(h-1)+2})\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_{d(h-1)+2}),\cdots ,{\widehat{b}}^H({\tilde{\mathrm{\theta }}}_{dh})\otimes \widehat{b}({\tilde{\mathrm{\theta }}}_{dh})]$. Obviously, the conventional sparse source is a special case of block sparse source by simply setting $d=1$. If ${\Vert p[h]\Vert }_2,h=1,2,\cdots ,H,$ has at most K nonzero elements, the vector $p$ is referred to as block K-sparse, where ${\Vert \cdot \Vert }_2$ denotes the Euclidean norm for vectors. In contrast to conventional sparsity, K is called block-sparsity.

Since $p$ has the same support as $q$, $p$ and $q$ are jointly sparse. Thus, a mixed k-l norm minimization problem [21] is utilized to jointly reconstruct $p$ and $q$. Given a $L\times 1$ block sparse vector $p$, the mixed k-l norm of $p$ is defined as

$${\Vert p\Vert }_{k,l}^l={\displaystyle \sum _{h=1}^H{({\Vert p[h]\Vert }_k)}^l}$$

(13)

Combining the definition given by Equation (13), this mixed k-l norm minimization problem is formulated as

$$\underset{p,q}{\min }({\Vert p\Vert }_{k,l}^l+\text{β}{\Vert q\Vert }_{k,l}^l) s.t.z=Cp+Dq+{\text{σ}}^{\text{2}}{\text{ΦΦ}}^H$$

(14)

where $0\le l\le 1$. It is worth mentioning that an important class of methods for solving the constrained optimization problem is to form the auxiliary function. By introducing the Lagrange multiplier method, the Lagrange function with respect to Equation (14) is given by

$$\underset{p,q}{\min }\frac{1}{2}{\Vert z-Cp-Dq\Vert }_2^2+{\text{β}}_p{\displaystyle \sum _{h=1}^H{({\Vert p[h]\Vert }_k)}^l}+{\text{β}}_q{\displaystyle \sum _{h=1}^H{({\Vert q[h]\Vert }_k)}^l}$$

(15)

where ${\text{β}}_p$ and ${\text{β}}_q$ are regularized parameters. As can be seen in Equation (15), the minimization problem with respect to $p$ and $q$ is nonconvex so as to make DOA estimation intractable. But note that if we fix one of two sparse vectors, that is, either $p$ or $q$, the minimization problem in Equation (15) turns out to be convex with respect to the other sparse vector alone. Thus, $p$ and $q$ can be reconstructed by alternately solving the following two minimization problems under the condition that the other sparse vector is fixed.

$$\underset{p}{\min }\{\frac{1}{2}{\Vert z-Cp-Dq\Vert }_2^2+{\text{β}}_p{\displaystyle \sum _{h=1}^H{({\Vert p[h]\Vert }_k)}^l\}}$$

(16)

$$\underset{q}{\min }\{\frac{1}{2}{\Vert z-Cp-Dq\Vert }_2^2+{\text{β}}_q{\displaystyle \sum _{h=1}^H{({\Vert q[h]\Vert }_k)}^l}\}$$

(17)

It is clear that Equations (16) and (17) have the same structure and can be solved in a similar manner. In the following, we only need to find the optimal solution for Equation (16). Then, the minimization problem in Equation (17) can be handled in the same way.

The objective function in Equation (16) can be expressed as

$$F(p)=f(p)+v(p)$$

(18)

where $f(p)=\frac{1}{2}{\Vert z-Cp-Dq\Vert }_2^2$ and $v(p)={\text{β}}_p{\displaystyle \sum _{h=1}^H{({\Vert p[h]\Vert }_k)}^l}$. Assume that $p^{(k)}$ is obtained at the kth iteration. Then, by the quadratic approximation of $f(p)$ at the fixed point $p^{(k)}$, $F(p)$ can be approximated as

$$F(p)\approx f(p^{(k)})+{(p-p^{(k)})}^H\nabla f(p^{(k)})+\frac{{\displaystyle \sum _{h=1}^H{\Vert p[h]-p^{(k)}[h]\Vert }_2^2}}{2\Sigma }+v(p)$$

(19)

where $\nabla f(p^{(k)})=C^H\cdot [{\displaystyle \sum _{h=1}^H(C[h]p^{(k)}[h]+D[h]q[h])-z}]$ is the partial derivation of $f(p^{(k)})$ with respect to $p^{(k)}$. $\Sigma$ is set to be $\Sigma =\frac{{\Vert \Delta \Vert }_1}{L}$, where ${\Vert \cdot \Vert }_1$ denotes one norm for vectors, $\Delta ={[{\Delta }_1,{\Delta }_2,\cdots ,{\Delta }_L]}^T$ is a $L\times 1$ vector and ${\Delta }_l=\frac{1}{{\Vert {\widehat{a}}^H({\tilde{\mathrm{\theta }}}_l)\otimes \widehat{a}({\tilde{\mathrm{\theta }}}_l)\Vert }_2^2}$, $l=1,2,\cdots ,L$. Moreover, $\Delta$ can be shown to be a block vector, i.e.,

$$\Delta ={[{\Delta }^T[1],{\Delta }^T[2],\cdots ,{\Delta }^T[H]]}^T$$

(20)

where $\Delta [h]$ is the hth block of $\Delta$ with the length d. For the convenience of analysis, denote $J^{(k)}=diag(\Delta )\nabla f(p^{(k)})$, and we can also represent $J^{(k)}$ as a block vector, i.e.,

$$J^{(k)}=diag(\Delta )C^H\cdot [{\displaystyle \sum _{h=1}^H(C[h]p^{(k)}[h]+D[h]q[h])-z}]={[{(J^{(k)}[1])}^T,{(J^{(k)}[2])}^T,\cdots ,{(J^{(k)}[H])}^T]}^T$$

(21)

where $J^{(k)}[h]$ is the hth block of $J^{(k)}$ with the length d. At the kth iteration, $p^{(k)}$ is updated by minimizing $F(p)$ in Equation (19) so that the next iteration $p^{(k+1)}$ is given by

$$p^{(k+1)}=\underset{p}{\min \{}f(p^{(k)})+{(p-p^{(k)})}^H\nabla f(p^{(k)})+\frac{{\displaystyle \sum _{h=1}^H{\Vert p[h]-p^{(k)}[h]\Vert }_2^2}}{2\text{Σ}}+v(p)\}$$

(22)

By further derivation, Equation (22) can be simplified to

$$p^{(k+1)}=\underset{p}{\min \{}{\displaystyle \sum _{h=1}^H[\frac{1}{2\text{Σ}}{\Vert p[h]-u^{(k)}[h]\Vert }_2^2}+{\lambda }_p{\Vert p[h]\Vert }_k^l]\}$$

(23)

where $u^{(k)}[h]=p^{(k)}[h]-J^{(k)}[h]$. As one may note, the objective function in Equation (23) is separable, and thus H blocks of $p^{(k+1)}$ can be solved in a parallel manner. Although it is hard to solve H blocks, the classical BCD method has provided a critical inspiration, fortunately. The solution to the hth block of $p^{(k+1)}$ is given by a soft-thresholding operator [22]

$$p^{(k+1)}[h]=\frac{u^{(k)}[h]}{{\Vert u^{(k)}[h]\Vert }_2}({\Vert u^{(k)}[h]\Vert }_2-{\text{λ}}_p{\Vert \Delta [h]\Vert }_2)I({\Vert u^{(k)}[h]\Vert }_2\ge {\text{λ}}_p{\Vert \Delta [h]\Vert }_2), h=1,2,\cdots ,H$$

(24)

where $I(\cdot )$ denotes an indicator function. Instead of reconstructing $p^{(k+1)}$ directly, $p^{(k+1)}$ can be reconstructed by H blocks, which may be zero vectors during the iteration. It can be seen in Equation (24) that we just need to determine the relation of ${\Vert u^{(k)}[h]\Vert }_2$ and ${\text{λ}}_p{\Vert \Delta [h]\Vert }_2$ to judge whether $p^{(k+1)}[h]$ is a zero vector. If ${\Vert u^{(k)}[h]\Vert }_2$ is less than ${\text{λ}}_p{\Vert \Delta [h]\Vert }_2$, $p^{(k+1)}[h]$ must be a zero vector. Therefore, a considerable amount of computations can be avoided in the process of solving H blocks. Similarly, by utilizing the soft-thresholding operator, the solution to the hth block of $q^{(k+1)}$ is expressed as

$$q^{(k+1)}[h]=\frac{{\overline{u}}^{(k)}[h]}{{\Vert {\overline{u}}^{(k)}[h]\Vert }_2}({\Vert {\overline{u}}^{(k)}[h]\Vert }_2-{\text{λ}}_q{\Vert \overline{\Delta }[h]\Vert }_2)I({\Vert {\overline{u}}^{(k)}[h]\Vert }_2\ge {\text{λ}}_q{\Vert \overline{\Delta }[h]\Vert }_2), h=1,2,\cdots ,H$$

(25)

where $\overline{\Delta }[h]$ and ${\overline{J}}^{(k)}[h]$ are the hth blocks of $\overline{\Delta }$ and $\overline{J}$ with the length d, respectively, and ${\overline{u}}^{(k)}[h]=q^k[h]-{\overline{J}}^{(k)}[h]$. Following the similar derivation, we have

$${\overline{J}}^{(k)}=diag(\overline{\Delta })D^H\cdot [{\displaystyle \sum _{h=1}^H(C[h]p[h]+D[h]q^{(k)}[h])-z}]$$

(26)

where $\overline{\Delta }={[{\overline{\Delta }}_1,{\overline{\Delta }}_2,\cdots ,{\overline{\Delta }}_L]}^T$ is a $L\times 1$ vector and ${\overline{\Delta }}_i=\frac{1}{{\Vert {\widehat{b}}^H({\tilde{\theta }}_l)\otimes \widehat{b}({\tilde{\theta }}_l)\Vert }_2^2}$, $l=1,2,\cdots ,L$. Based on Equations (24) and (25), $p$ and $q$ can be reconstructed alternately until the following criterion is satisfied

$$\Vert \frac{{\text{η}}^{(k+1)}-{\text{η}}^{(k)}}{{\text{η}}^{(k)}}\Vert \le \text{γ}$$

(27)

where ${\text{η}}^{(k)}=[p^{(k)};q^{(k)}]$ and $\text{γ}$ is the small tolerance.

To reduce the computational complexity, a block selection criterion is given. This criterion is of great importance in the whole ABCD method. By utilizing the block selection criterion, we can only update the block that is the closest to $z$, i.e.,

$$h_0=\arg \min {\Vert w[h]-z\Vert }_2, h=1,2,\cdots ,H$$

(28)

where $w[h]=C[h]p[h]+D[h]q[h]$. This means that $p^{(k+1)}[h_0]$ and $q^{(k+1)}[h_0]$ at the kth iteration are updated as Equations (24) and (25) while the remaining blocks keep unchanged. The purpose of utilizing this block selection criterion is to avoid the update of repetitive and unnecessary blocks and reduce the computational complexity. The major steps of reconstructing $p$ and $q$ by the ABCD method are given as follows:

Initialization: set $k=0,p^{(0)}=0,q^{(0)}=0,C\in {ℂ}^{N^2\times L},D\in {ℂ}^{N^2\times L},z=vec(R_Y)\in {ℂ}^{N^2\times 1}$ and $\text{γ}={10}^{-5}$.

(1)

Calculate $J^{(k)}$ and ${\overline{J}}^{(k)}$ in terms of Equations (21) and (26).

(2)

Due to H blocks of $J^{(k)}$ and ${\overline{J}}^{(k)}$, calculate $u^{(k)}[h]$ and ${\overline{u}}^{(k)}[h]$, $h=1,2,\cdots ,H$.

(3)

Calculate $p^{(k+1)}[h]$ and $q^{(k+1)}[h]$, $h=1,2,\cdots ,H$, in terms of (24) and (25).

(4)

Choose the block index $h_0$ according to (28). Then, $p^{(k+1)}$ and $q^{(k+1)}$ are respectively updated as ${[p^{(k)}[1],\cdots ,p^{(k)}[h_0-1],p^{(k+1)}[h_0],p^{(k)}[h_0+1],\cdots ;p^{(k)}[H]]}^T$ and $[q^{(k)}[1],\cdots ,q^{(k)}[h_0-1],$ $q^{(k+1)}[h_0],q^{(k)}[h_0+1],\cdots ,q^{(k)}[H]{]}^T$.

(5)

If $\Vert \frac{{\text{η}}^{(k+1)}-{\text{η}}^{(k)}}{{\text{η}}^{(k)}}\Vert \le \text{γ}$, stop the iteration. Otherwise, set $k=k+1$ and return to step (1).

4. Global Convergence of the ABCD Method

The global convergence of the ABCD method is proved in this section. By combining the existing convergence proof of the general BCD framework [23] with ABCD method, a detailed derivation process for proving the global convergence is shown as follows.

First, we introduce the general BCD framework. Note that $f(p)$ in Equation (18) is a continuous convex function and $v(p)$ in (18) is a non-smooth convex function. Given a fixed point $p^{(k)}$, $F(p)$ in Equation (18) can be approximated as the following form by exploiting the second order Taylor expansion of $f(p)$ in the general BCD framework.

$$F(p)\approx f(p^{(k)})+{(p-p^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{(p-p^{(k)})}^{H}H(p-p^{(k)})+v(p)$$

(29)

where $H=\nabla f^2(p^{(k)})=C^{H}C$ is the $L\times L$ Hessian matrix of $f(p^{(k)})$ with respect to $p^{(k)}$. It is clear that minimizing Equation (29) involves finding the next iteration $p^{(k+1)}$. Assume that $\text{χ}=[1,2,\cdots ,H]$ is an H-dimensional index set and ${\text{χ}}_0$ is a subset of $\text{χ}$ consisting of at most K indexes obtained by Equation (28). Hence, the next iteration $p^{(k+1)}$ is represented as

$$p^{(k+1)}=\Lambda (p^{(k)},{\text{χ}}_0)=\arg \underset{p}{\min }\{F(p)/p^{(k)}[h]=0,h\in \text{χ},h\notin {\text{χ}}_0\}$$

(30)

Then, to prove the global convergence, we give the modified Armijo rule and modified Gauss-Southwell-r rule that are prerequisites to guarantee the global convergence. These two rules are described in the following.

(1)

Modified Armijo rule: If $0<{\text{α}}^{(k)}<1$, $0<\text{ξ}<1$ and $0<\text{κ}<1$, the following inequality holds:

$$F(p^{(k)}+{\text{α}}^{(k)}{\text{μ}}^{(k)})\le F(p^{(k)})+{\text{α}}^{(k)}\text{κ}{\Gamma }^{(k)}$$

(31)

where ${\text{μ}}^{(k)}=p-p^{(k)}$ and ${\Gamma }^{(k)}={({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\text{ξ}{({\text{μ}}^{(k)})}^{H}H{\text{μ}}^{(k)}+v(p)-v(p^{(k)})$.

(2)

Modified Gauss-Southwell-r rule: In the iteration, the index set ${\text{χ}}_0$ obtained by Equation (28) must satisfy

$${\Vert \Lambda (p^{(k)},{\chi }_0)\Vert }_2\le \text{ς}{\Vert \Lambda (p^{(k)},\text{χ})\Vert }_2$$

(32)

where $0<\text{ς}<1$.

It is well known that the problem of proving the global convergence is quite complex and intractable. However, fortunately, since Equations (16) and (17) are both convex and have both only one global point, we only need to show that Equations (16) and (17) satisfy the modified Armijo rule and modified Gauss-Southwell-r rule to prove the global convergence of the ABCD method. Furthermore, since Equation (16) has the same structure as Equation (17), it is enough to just prove that Equation (16) satisfies the modified Armijo rule and modified Gauss-Southwell-r rule. Regarding Equation (17), the derivation process is given in the same way.

To see the first, the following inequality holds according to Equation (30).

$$F(p)=F(p^{(k)}+{\text{μ}}^{(k)})\le F(p^{(k)}+{\text{α}}^{(k)}{\text{μ}}^{(k)})$$

(33)

By substituting Equation (33) into Equation (29), we have

$$\begin{array}{l}f(p^{(k)})+{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p^{(k)}+{\text{μ}}^{(k)})\le \\ f(p^{(k)})+{({\text{α}}^{(k)}{\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{α}}^{(k)}{\text{μ}}^{(k)})}^{H}H({\text{α}}^{(k)}{\text{μ}}^{(k)})+v(p^{(k)}+{\text{α}}^{(k)}{\text{μ}}^{(k)})\end{array}$$

(34)

Based on the fact that $v(p^{(k)}+{\text{α}}^{(k)}{\text{μ}}^{(k)})\le {\text{α}}^{(k)}v(p^{(k)}+{\text{μ}}^{(k)})+(1-{\text{α}}^{(k)})v(p^{(k)})$, Equation (34) can be further simplified to

$$(1-{\text{α}}^{(k)}){({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}[1-{({\text{α}}^{(k)})}^2]{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+(1-{\text{α}}^{(k)})[v(p^{(k)}+{\text{μ}}^{(k)})-v(p^{(k)})]\le 0$$

(35)

Since $0<{\text{α}}^{(k)}<1$, by dividing by $1-{\text{α}}^{(k)}$ and setting ${\text{α}}^{(k)}\to 1$, we have

$${({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p^{(k)}+{\text{μ}}^{(k)})-v(p^{(k)})\le 0$$

(36)

For $0<\text{ξ}<1$, it can be deduced from (36) that

$${\Gamma }^{(k)}={({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\text{ξ}{({\text{μ}}^{(k)})}^{H}H{\text{μ}}^{(k)}+v(p)-v(p^{(k)})\le (\text{ξ}-1){({\text{μ}}^{(k)})}^{H}H{\text{μ}}^{(k)}\le 0$$

(37)

Subsequently, by exploiting the convexity of $f(p)$, we obtain

$$f(p)=f(p^{(k)}+{\text{μ}}^{(k)})\le f(p^{(k)})+{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})$$

(38)

Since $F(p)=f(p)+v(p)$, it is nature to have

$$F(p)=F(p^{(k)}+{\text{μ}}^{(k)})\le F(p^{(k)})+{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p)-v(p^{(k)})$$

(39)

Following the fact in Equation (37), the following inequality holds

$$\begin{array}{l}F(p^{(k)})+{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p)-v(p^{(k)})\\ \le F(p^{(k)})+\text{κ}[{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p)-v(p^{(k)})]\end{array}$$

(40)

where $0<\text{κ}<1$. Combining Equations (39) and (40), we have

$$F(p^{(k)}+{\text{μ}}^{(k)})\le F(p^{(k)})+\text{κ}[{({\text{μ}}^{(k)})}^H\nabla f(p^{(k)})+\frac{1}{2}{({\text{μ}}^{(k)})}^{H}H({\text{μ}}^{(k)})+v(p)-v(p^{(k)})]$$

(41)

It is worth pointing out that Equation (41) is equal to Equation (31) with ${\text{α}}^{(k)}=1$ and $\text{ξ}=\frac{1}{2}$. Thus, it has been proved that Equation (16) satisfies the modified Armijo rule.

Secondly, in order to prove that Equation (16) satisfies the modified Gauss-Southwell-r rule, the following form is shown

$$h_0=\arg \min {\Vert w[h]-z\Vert }_2=\arg \min {\Vert \Lambda (p^{(k)},{\text{χ}}_i)\Vert }_2$$

(42)

where the cardinality of ${\text{χ}}_i$ is the same as that of ${\text{χ}}_0$. Without loss of generality, consider the worst case, i.e., the cardinality of ${\text{χ}}_i$ is the maximization K, and assume $H=KZ$ so that $\text{χ}$ is expressed as $\text{χ}=[{\text{χ}}_1，{\text{χ}}_2，\cdots ，{\text{χ}}_Z]$. Thus, we have

$${\Vert \Lambda (p^{(k)},{\text{χ}}_0)\Vert }_2\le {\Vert \Lambda (p^{(k)},{\text{χ}}_i)\Vert }_2$$

(43)

Based on the following equation

$${\Vert \Lambda (p^{(k)},\text{χ})\Vert }_2=\sqrt{{\displaystyle \sum _{z=1}^Z{\Vert \Lambda (p^{(k)},{\text{χ}}_z)\Vert }_2^2}}$$

(44)

it is easy to have

$${\Vert \Lambda (p^{(k)},{\text{χ}}_0)\Vert }_2\le \frac{1}{\sqrt{Z}}{\Vert \Lambda (p^{(k)},\text{χ})\Vert }_2$$

(45)

Since Equation (45) is equal to Equation (32) with $\text{ς}=\frac{1}{\sqrt{Z}}$, Equation (16) satisfies the modified Gauss-Southwell-r rule. Therefore, based on the above analysis, the global convergence of the ABCD method has been proved.

5. Simulation Results

This section presents several simulations to validate the superior performance of the proposed method as compared with R-GBCD+ and OGSBI-SVD. The angular space [−90°, 90°] is taken the grid with grid interval τ = 3° to perform three methods for off-grid targets. We set the length of block, the number of ULA sensors and spacing between adjacent sensors to be $d=3$, $M=9$ and $\stackrel{⌢}{d}=\frac{\text{λ}}{2}$, respectively. In the simulation, the root mean squared error (RMSE) and success rate of DOA estimation are two significant performance indexes. RMSE is defined as

$$RMSE=\sqrt{{\displaystyle \sum _{i=1}^{Mc}{\displaystyle \sum _{k=1}^K\frac{{({\overline{\mathrm{\theta }}}_{k,i}-{\mathrm{\theta }}_k)}^2}{K\cdot Mc}}}}$$

(46)

where $Mc$ is the number of Monte Carlo runs and ${\overline{\mathrm{\theta }}}_{k,i}$ is the estimate of ${\mathrm{\theta }}_k$ in the ith Monte Carlo run, and success rate is declared if the estimation error is within a certain small Euclidean distance of the true directions.

In the first simulation, we compare the spatial spectra of R-GBCD+, OGSBI-SVD and ABCD. Consider four far-field narrowband sources impinging on the ULA from [−30.4° −3.8° 10.1° 15.3°], where the latter two most closely spaced sources are coherent and the remaining sources are independent of other sources. Figure 1 presents the spatial spectra of R-GBCD+, OGSBI-SVD and ABCD with SNR 3 dB and number of snapshots 100. For the convenience of analysis, the spatial spectra are normalized. We can see from Figure 1 that the spatial spectra of three methods are able to detect four sources, but the spatial spectrum obtained by R-GBCD+ has obvious bias at the true directions and OGSBI-SVD can yield slight bias in the vicinity of the coherent sources. Note that ABCD has a nearly ideal spatial spectrum, and thus it outperforms R-GBCD+ and OGSBI-SVD in terms of the spatial spectrum.

Figure 1. Spatial spectra of R-GBCD+, OGSBI-SVD and ABCD, where the pink circles denote the true DOAs.

Figure 1. Spatial spectra of R-GBCD+, OGSBI-SVD and ABCD, where the pink circles denote the true DOAs.

The success rates of three methods vs. SNR and the number of snapshots are analyzed in the second simulation. The source mode is the same as the first simulation. Figure 2 shows the success rates of three methods vs. SNR with the fixed number of snapshots 100, whereas the success rates of three methods vs. number of snapshots are depicted with the fixed SNR 0 dB in Figure 3. The following facts can be acquired from Figure 2 and Figure 3 that three methods can estimate correctly for high SNR or large number of snapshots and ABCD has a higher success rate than the other two methods for low SNR or a small number of snapshots.

Figure 2. Success rates vs. SNR with the fixed number of snapshots 100.

Figure 2. Success rates vs. SNR with the fixed number of snapshots 100.

Figure 3. Success rates vs. number of snapshots with the fixed SNR 0 dB.

Figure 3. Success rates vs. number of snapshots with the fixed SNR 0 dB.

Figure 4. RMSE vs. SNR with the fixed number of snapshots 100.

Figure 4. RMSE vs. SNR with the fixed number of snapshots 100.

Figure 5. RMSE vs. number of snapshots with the fixed SNR 0 dB.

Figure 5. RMSE vs. number of snapshots with the fixed SNR 0 dB.

The third simulation considers the RMSE of three methods vs. SNR and the number of snapshots. All the conditions are the same as the second simulation. Figure 4 and Figure 5 show the RMSE of three methods vs. SNR and the number of snapshots, respectively. It is indicated in Figure 4 and Figure 5 that ABCD has the best estimation accuracy among all three methods. Moreover, the accuracy of three methods is gradually improving with SNR or the number of snapshots increasing.

Finally, we test the resolving ability by showing the relation between RMSE of three methods and angle separation of sources, which is illustrated in Figure 6 Consider two coherent sources impinging on the ULA from $30.7$ and $30.7+\mathrm{\Delta }\mathrm{\theta }$, where the step of $\mathrm{\Delta }\mathrm{\theta }$ is 1°. The SNR is 0 dB and the number of snapshots is 100. As can be seen from Figure 6, the performance of R-GBCD+ and OGSBI-SVD degrades severely as angle separation is 3°, while ABCD can still provide a precise estimation as long as angle separation is no less than 3°. The proposed ABCD is the most accurate method and has higher resolution than the other two methods.

Figure 6. RMSE vs. angle separation with the fixed SNR 0 dB and number of snapshots 100 for coherent sources.

Figure 6. RMSE vs. angle separation with the fixed SNR 0 dB and number of snapshots 100 for coherent sources.

6. Conclusions

In this paper, a novel ABCD method is proposed for off-grid DOA estimation in CS. The proposed method minimizes the mixed k-l norm to reconstruct the sparse source and estimate the grid error. In order to make the minimization problem tractable, an iterative process that minimizes the mixed k-l norm alternately over two sparse vectors is adopted. By reconstructing the block sparse source instead of conventional sparse source, the proposed method can achieve the better reconstruction properties. A block selection criterion is given to update the block so that the proposed method can reduce computational complexity. It is proved that the proposed method has the global convergence. Simulation results show that the proposed method has more notable performance advantages than R-GBCD+ and OGSBI-SVD in terms of spatial spectrum, RMSE and success rate.