Improved Spatial Differencing Scheme for 2-D DOA Estimation of Coherent Signals with Uniform Rectangular Arrays

Abstract

This paper proposes an improved spatial differencing (ISD) scheme for two-dimensional direction of arrival (2-D DOA) estimation of coherent signals with uniform rectangular arrays (URAs). We first divide the URA into a number of row rectangular subarrays. Then, by extracting all the data information of each subarray, we only perform difference-operation on the auto-correlations, while the cross-correlations are kept unchanged. Using the reconstructed submatrices, both the forward only ISD (FO-ISD) and forward backward ISD (FB-ISD) methods are developed under the proposed scheme. Compared with the existing spatial smoothing techniques, the proposed scheme can use more data information of the sample covariance matrix and also suppress the effect of additive noise more effectively. Simulation results show that both FO-ISD and FB-ISD can improve the estimation performance largely as compared to the others, in white or colored noise conditions.

1. Introduction

Two-dimensional direction of arrival (2-D DOA) (i.e., elevation and azimuth angles) estimation of multiple signals with different array geometries is an important problem in many practical applications such as radars and wireless communications. Various methods have been developed for solving this problem [1,2,3,4,5,6,7,8,9,10], such as the subspace-based methods [1,2,3,4,5,6], the sparse reconstruction methods [7,8,9], and the least-square approach [10].

However, for coherent signals [11,12,13], the traditional methods, such as the multiple signal classification (MUSIC) [14] and estimation of signal parameters via rotational invariance technique (ESPRIT) [15], suffer from performance degradation due to the rank-deficiency of the signal covariance matrix. The forward backward spatial smoothing (FBSS) technique [16] is very effective in removing the coherency but at the cost of few degrees of freedom (DOFs). Therefore, spatial smoothing techniques with URAs are developed by using virtual sensors to increase the DOFs, e.g., the unitary ESPRIT [17] and spatial smoothing MUSIC [18] algorithms. The 2-D spatial smoothing methods were also applied for DOA estimation [19] or joint DOA and direction of departure (DOD) estimation [20] with multi-input multi-output (MIMO) radar, where the transmission and reception diversity smoothing is derived by constructing a new covariance matrix with decorrelated signal subspace. Besides, compared with the higher-order cumulants methods [21,22], spatial smoothing methods also have a lower computational complexity.

Recently, spatial differencing techniques [23,24,25,26,27] have been introduced to suppress the effect of additive noise. Aiardi et al. [23] proposed a high-resolution DOA estimation method by performing the partial spatial differencing operation; the methods in [24,25] can suppress white noise or colored noise by using the difference between the first and backward subarrays or between the neighboring subarrays, respectively. Liu et al. [26] developed a generalized covariance differencing algorithm by using the difference between the FB smoothing matrix and its complex conjugation. In fact, these methods [23,24,25,26] cannot be directly applied for 2-D DOA estimation. Then, similar with the 2-D configuration, the method in [27] explored four kinds of smoothing techniques including the spatial difference smoothing (SDS), asymmetric SDS (A-SDS), transmit-receive diversity SDS (TRD-SDS), and asymmetric TRD-SDS (A-TRD-SDS) for coherent targets with MIMO radar.

Nevertheless, the aforementioned spatial differencing techniques haven’t fully explored the advantages of the URA. On the one hand, the involved spatial smoothing subarrays can only use part of the data information of the sample covariance matrix. On the other, the classical spatial differencing methods only focus on the suppression of additive noise and few ones consider the information loss caused by the difference-operation. Therefore, in this paper, we propose an improved spatial differencing (ISD) scheme for 2-D DOA estimation of coherent signals with URAs, where both the forward only ISD (FO-ISD) and forward backward ISD (FB-ISD) methods are developed. Simulation results show the usefulness of the proposed methods. For clarity, the main advantages are given as follows:

• Classic spatial differencing techniques only use the data information of overlapped smoothing subarrays, while FO-ISD and FB-ISD can extract all the data information of each row or column rectangle subarrays.
• Classic spatial differencing techniques perform difference-operation on the whole overlapping subarrays, while FO-ISD and FB-ISD calculate the differencing matrix for the auto-correlations and keep the cross-correlations unchanged. So SD-SMS has less information loss, resulting in a more effective noise suppression.
• FB-ISD can achieve a further improved performance than FO-ISD due to the increased number of smoothing submatrices.

The rest of the paper is listed as follows. We first introduce the basic signal model of URA for coherent signals in Section 2. In Section 3, we develop the FO-ISD and FB-ISD methods by using the row or column rectangular subarrays, where the Cramér-Rao bound (CRB) is also given. Simulation results are presented in Section 4.1 and we conclude this paper in Section 5.

In this paper, operators ${(·)}^T$, ${(·)}^{*}$, and ${(·)}^H$ represent transpose, conjugation, and conjugate transpose, respectively. ${\mathbf{I}}_N$ denotes an $N\times N$ identity matrix and ${\mathbf{J}}_M$ denotes an $M\times M$ exchange matrix with ones on its anti-diagonal and zeros elsewhere. ∘ and ⊕ represent the Khatri-Rao product and Hadamard product, respectively; $diag(·)$ and $blkdiag(·)$ denote the diagonal matrix or block diagonal matrix operator. $E[·]$ and $vec(·)$ denote expectation and vectorization, respectively.

2. System Model

As described in Figure 1, we consider K narrowband far-field coherent signals $s_k\left(t\right)$ ($k=1,2,\cdots ,K$) impinging on a URA ($M\times N$ sensors). We assume both x and y directions of the URA are separated by half a wavelength. Then the output can be expressed as [11]

$$\begin{array}{c}\hfill \mathbf{X}\left(t\right)=\sum _{k=1}^K{\mathbf{a}}_x({\alpha }_k,{\theta }_k){\mathbf{a}}_y^T({\alpha }_k,{\theta }_k)s_k\left(t\right)+\mathbf{Z}\left(t\right),\end{array}$$

(1)

where ${\theta }_k$ and ${\alpha }_k$ are the elevation and azimuth angles of the k-th signal; ${\mathbf{a}}_x({\alpha }_k,{\theta }_k)={\mathbf{a}}_x\left(u_k\right)={[1,{\mathrm{e}}^{-j\pi u_k},\dots ,{\mathrm{e}}^{-j\pi (M-1)u_k}]}^T$, ${\mathbf{a}}_y({\alpha }_k,{\theta }_k)={\mathbf{a}}_y\left(v_k\right)=[1,{\mathrm{e}}^{-j\pi v_k},\dots ,{\mathrm{e}}^{-j\pi (N-1)v_k}{]}^T$, $u_k=sin{\theta }_{k}cos{\alpha }_k$, $v_k=sin{\theta }_{k}sin{\alpha }_k$. The elements of $\mathbf{Z}\left(t\right)$ are temporally and spatially complex white Gaussian noises with zero-mean and variance ${\sigma }^2$. Then, vectorizing $\mathbf{X}\left(t\right)$ yields

$$\begin{array}{c}\hfill \mathbf{x}\left(t\right)=vec\left(\mathbf{X}\left(t\right)\right)=({\mathbf{A}}_x\circ {\mathbf{A}}_y)\mathbf{s}\left(t\right)+\mathbf{z}\left(t\right),\end{array}$$

(2)

where ${\mathbf{A}}_x={[{\mathbf{a}}_x\left(u_1\right),\cdots ,{\mathbf{a}}_x\left(u_K\right)]}^T$, ${\mathbf{A}}_y={[{\mathbf{a}}_y\left(v_1\right),\cdots ,{\mathbf{a}}_y\left(v_K\right)]}^T$, $\mathbf{s}\left(t\right)={[s_1\left(t\right),\cdots ,s_K\left(t\right)]}_{K\times 1}^T$, and $\mathbf{z}\left(t\right)$ = $vec\left(\mathbf{Z}\right(t\left)\right)$. With L snapshots, we calculate the sample covariance matrix as

$$\begin{array}{c}\hfill {\mathbf{R}}_0=E\left[\mathbf{x}\left(t\right){\mathbf{x}}^H\left(t\right)\right]=\frac{1}{L}\sum _{t=1}^L\mathbf{x}\left(t\right){\mathbf{x}}^H\left(t\right)=\mathbf{A}{\mathbf{R}}_s{\mathbf{A}}^H+{\sigma }^2{\mathbf{I}}_{MN}.\end{array}$$

(3)

where $t=1,2,\cdots ,L$, $\mathbf{A}={\mathbf{A}}_x\circ {\mathbf{A}}_y$, and ${\mathbf{R}}_s=E\left[\mathbf{s}\left(t\right){\mathbf{s}}^H\left(t\right)\right]$ represents the signal covariance matrix.

3. 2-D DOA Estimation with URA

In this section, we firstly review the classic spatial differencing technique in [26] and then we derive the FO-ISD and FB-ISD methods in detail.

3.1. Classic Spatial Differencing Technique

The main idea of the classic spatial differencing technique is to build $Q_x{Q}_y$ overlapping rectangular subarrays with size of $P_x\times P_y$, i.e., sliding windows, where $Q_x$ and $Q_y$ are the forward subarrays along the x and y directions, respectively, $Q_x=M-P_x+1$, $Q_y=N-P_y+1$. Then, we can get the $(q_x,q_y)$-th sliding window as

$$\begin{array}{c}\hfill {\mathbf{x}}_{q_x{q}_y}\left(t\right)={\mathbf{A}}_P{\mathsf{\Phi }}_x^{q_x-1}{\mathsf{\Phi }}_y^{q_y-1}\mathbf{s}\left(t\right)+{\mathbf{z}}_{q_x{q}_y}\left(t\right),\end{array}$$

(4)

where ${\mathbf{A}}_P={\mathbf{A}}_{P_x}\circ {\mathbf{A}}_{P_y}$, ${\mathbf{A}}_{P_x}$ is the submatrix of the array response matrix ${\mathbf{A}}_x$ consisting of the first row $P_x$ and the submatrix ${\mathbf{A}}_{P_y}$ of ${\mathbf{A}}_y$ consisting of the first row $P_y$; $q_x=1,\cdots ,Q_x$, $q_y=1,\cdots ,Q_y$, ${\mathsf{\Phi }}_x=diag$$[{\mathrm{e}}^{-j\pi u_1},\cdots ,{\mathrm{e}}^{-j\pi u_K}]$, ${\mathsf{\Phi }}_y=diag[{\mathrm{e}}^{-j\pi v_1},\cdots ,{\mathrm{e}}^{-j\pi v_K}]$, and ${\mathbf{z}}_{q_x{q}_y}\left(t\right)$ is the corresponding noise vector. Then the covariance submatrix of ${\mathbf{x}}_{q_x{q}_y}\left(t\right)$ can be given as

$$\begin{array}{c}\hfill {\mathbf{R}}_{q_x{q}_y}={\mathbf{A}}_P{\mathsf{\Phi }}_x^{q_x-1}{\mathsf{\Phi }}_y^{q_y-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-q_y}{\mathsf{\Phi }}_x^{1-q_x}{\mathbf{A}}_P^H+{\sigma }^2{\mathbf{I}}_{P_x{P}_y},\end{array}$$

(5)

Using the covariance submatrix in (5), we can build the $p_a$-th SDS (i.e., asymmetric SDS, A-SDS) matrix as [26]

$$\begin{array}{c}\hfill {\mathbf{R}}_{p_a}^f=\frac{1}{Q_x{Q}_y}\sum _{q_x=1}^{Q_x}\sum _{q_y=1}^{Q_y}{\mathbf{R}}_{q_x{q}_y}-\left(\frac{1}{Q_x{Q}_y-p_a}\right)\sum _{q_x=1,q_x\notin {\mathbf{m}}_a}^{Q_x}\sum _{q_y=1,q_y\notin {\mathbf{n}}_a}^{Q_y}{\mathbf{R}}_{q_x{q}_y}^b,\end{array}$$

(6)

where ${\mathbf{m}}_a=[m_a^1,m_a^2,\cdots ,m_a^{p_a}]$, $m_a^i\in [1,2,\cdots ,Q_x]$, and $1\le i\le p_a$; ${\mathbf{n}}_a=[n_a^1,$$n_a^2,\cdots ,n_a^{f_a}]$, $n_a^j\in [1,2,\cdots ,Q_y]$, $1\le j\le f_a$; ${\mathbf{R}}_{q_x{q}_y}^b={\mathbf{J}}_{P_x{P}_y}{\mathbf{R}}_{q_x{q}_y}^{*}{\mathbf{J}}_{P_x{P}_y}$.

We can see that A-SDS exploits the asymmetric difference between the complete forward spatially smoothed matrix and incomplete backward spatially smoothed matrix, and the noise can be suppressed by the differencing matrices in (6). However, since the forward and backward smoothed matrices have similar data structures, the difference-operation will have great information loss and performance will decrease greatly in the low SNR condition.

3.2. Improved Spatial Differencing (ISD) Scheme

3.2.1. Analysis for Row Rectangular Subarrays

As described in Figure 2a, we divide the URA into $Q_x$ row rectangular subarrays along the x direction, where each one contains $Q_y$ sliding windows along the y direction. We take the first one as an example and write it as ${\mathbf{x}}_{1x}\left(t\right)={[{\mathbf{x}}_{1,1}^T,\cdots ,{\mathbf{x}}_{1,N}^T]}^T$, where ${\mathbf{x}}_{1,n}$ denotes the first row $P_x$ of ${\mathbf{x}}_n$, ${\mathbf{x}}_n={\mathbf{A}}_x{\mathsf{\Phi }}_y^{n-1}\mathbf{s}\left(t\right)+{\mathbf{z}}_n\left(t\right)$ is the received signal of the sensors in the n-th column of URA, $n=1,\cdots ,N$. Figure 2b describes the covariance matrix ${\mathbf{R}}_{1x}$ of ${\mathbf{x}}_{1x}\left(t\right)$. We can see that the data information within the blue box is constructed by the covariance submatrices ${\mathbf{R}}_{1q_y}\left(t\right)$ ($q_x=1$) in (5). Therefore, in each row rectangular subarray, the sliding windows (within the blue box) can only use part of the data information of the covariance matrix ${\mathbf{R}}_{mx}$, ($m=1,2,\cdots M$), where ${\mathbf{R}}_{mx}$ represents the covariance submatrix of ${\mathbf{x}}_{mx}\left(t\right)$. To fully use this data information and also decrease the information loss caused by difference-operation, the FO-ISD and FB-ISD methods are given as follows.

3.2.2. Forward only ISD (FO-ISD) Method

We continue to take the covariance matrix ${\mathbf{R}}_{1x}$ in Figure 2b as an example. Since the matrix ${\mathbf{R}}_{1x}$ is symmetric, we can only extract the submatrices below the diagonal ones from top to bottom. Each column submatrix of ${\mathbf{R}}_{1x}$ can be divided into some column submatrices block, and the information of n-th $(n=1,\cdots ,Q_y-1)$ column can be set as

$${\mathbf{R}}_{1,n}=\left\{\left[\begin{array}{c}{\mathbf{x}}_{1,n}\\ {\mathbf{x}}_{1,n+1}\\ ⋮\\ {\mathbf{x}}_{1,n+P_y-1}\end{array}\right]{\mathbf{x}}_{1,n}^H,\cdots ,\left[\begin{array}{c}{\mathbf{x}}_{1,Q_y}\\ {\mathbf{x}}_{1,Q_y+1}\\ ⋮\\ {\mathbf{x}}_{1,N}\end{array}\right]{\mathbf{x}}_{1,n}^H\right\}=\left[\begin{array}{c}{\mathbf{A}}_{P_x}\\ {\mathbf{A}}_{P_x}{\mathsf{\Phi }}_y\\ ⋮\\ {\mathbf{A}}_{P_x}{\mathsf{\Phi }}_y^{P_y-1}\end{array}\right]{\mathsf{\Pi }}_{1,n}+{\sigma }^2{\mathbf{C}}_n={\mathbf{A}}_P{\mathsf{\Pi }}_{1,n}+{\sigma }^2{\mathbf{C}}_n,$$

(7)

where ${\mathbf{x}}_{1,n}\left(t\right)={\mathbf{A}}_{P_x}{\mathsf{\Phi }}_y^{n-1}\mathbf{s}\left(t\right)+{\mathbf{z}}_{1,n}\left(t\right)$, ${\mathbf{z}}_{1,n}\left(t\right)$ is the corresponding noise vector with $\mathbf{z}\left(t\right)$; ${\mathsf{\Pi }}_{1,n}=\left\{{\mathsf{\Phi }}_y^{n-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-n}{\mathbf{A}}_{P_x}^H,{\mathsf{\Phi }}_y^n{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-n}{\mathbf{A}}_{P_x}^H,\right.\cdots ,\left.{\mathsf{\Phi }}_y^{Q_y-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-n}{\mathbf{A}}_{P_x}^H\right\}$, and ${\mathbf{C}}_n$$=[{\mathbf{c}}_1,{\mathbf{0}}_{P_x{P}_y\times P_x},\cdots ,{\mathbf{0}}_{P_x{P}_y\times P_x}]$, and ${\mathbf{c}}_1=$${[{\mathbf{I}}_{P_x},{\mathbf{0}}_{P_x},\cdots ,{\mathbf{0}}_{P_x}]}^T$. From (9), We observe that only the first submatrix contains auto-correlations, while the others are constructed from the cross-correlations. So we can perform the difference-operation on the first one, i.e.,

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbf{x}}_{1,1}\\ {\mathbf{x}}_{1,2}\\ ⋮\\ {\mathbf{x}}_{1,P_y}\end{array}\right]{\mathbf{x}}_{1,1}^H-{\mathbf{J}}_{P_x{P}_y}{\left\{\left[\begin{array}{c}{\mathbf{x}}_{1,n}\\ {\mathbf{x}}_{1,n+1}\\ ⋮\\ {\mathbf{x}}_{1,n+P_y-1}\end{array}\right]{\mathbf{x}}_{1,n}^H\right\}}^{*}{\mathbf{J}}_{P_x}={\mathbf{A}}_P\left({\mathbf{R}}_s-{\mathsf{\Theta }}^{*}{\mathsf{\Phi }}_y^{1-n}{\mathbf{R}}_s^{*}{\mathsf{\Phi }}_y^{n-1}{\mathsf{\Phi }}_x^{P_x-1}\right){\mathbf{A}}_{P_x}^H,\end{array}$$

(8)

where ${\mathbf{J}}_{P_x{P}_y}{\mathbf{A}}_P^{*}={\mathbf{A}}_p{\mathsf{\Theta }}^{*}$, $\mathsf{\Theta }=diag\{{\mathrm{e}}^{-j\pi d_1},\cdots ,{\mathrm{e}}^{-j\pi d_K}\}$, and $d_k=(P_x-1)u_k+(P_y-1)v_k$; ${\mathbf{J}}_{P_x}{\mathbf{A}}_{P_x}^{*}={\mathbf{A}}_{P_x}{\mathsf{\Phi }}_x^{P_x-1}$. Using the differencing matrix in (8) to replace the first submatrix in (7), we have

$$\begin{array}{c}\hfill {\overline{\mathbf{R}}}_{1,n}={\mathbf{A}}_P{\overline{\mathsf{\Pi }}}_{1,n}={\mathbf{A}}_P{\overline{\mathbf{H}}}_{1,n}diag\left({\mathbf{A}}_{P_x}^H,\cdots ,{\mathbf{A}}_{P_x}^H\right).\end{array}$$

(9)

where ${\overline{\mathbf{H}}}_{1,n}=\left\{\left({\mathbf{R}}_s-{\mathsf{\Theta }}^{*}{\mathsf{\Phi }}_y^{1-n}{\mathbf{R}}_s^{*}{\mathsf{\Phi }}_y^{n-1}{\mathsf{\Phi }}_x^{P_x-1}\right),\cdots ,{\mathsf{\Phi }}_y^{Q_y-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-n}\right\}$, and the number of ${\mathbf{A}}_{P_x}^H$ is $Q_y$ − n+1. We see that the matrix ${\overline{\mathbf{R}}}_{1,n}$ can suppress the effect of noise by performing the difference-operation on the auto-correlations.

Then, the remaining $P_y$ columns of ${\mathbf{R}}_{1x}$ (i.e., the last sliding window) can be expressed as

$$\begin{array}{c}\hfill {\mathbf{R}}_{1,Q_y}={\mathbf{A}}_P{\mathsf{\Phi }}_y^{Q_y-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-Q_y}{\mathbf{A}}_P^H+{\sigma }^2{\mathbf{I}}_{P_x{P}_y},\end{array}$$

(10)

Likewise, performing the difference-operation on the matrix ${\mathbf{R}}_{1,Q_y}$, we have

$$\begin{array}{c}\hfill {\mathbf{R}}_{1,Q_y}^d={\mathbf{R}}_{1,1}-{\mathbf{J}}_{P_x{P}_y}{\mathbf{R}}_{1,Q_y}^{*}{\mathbf{J}}_{P_x{P}_y}={\mathbf{A}}_P{\mathsf{\Pi }}_{1,Q_y}={\mathbf{A}}_P{\mathbf{H}}_{1,Q_y}{\mathbf{A}}_P^H,\end{array}$$

(11)

where ${\mathbf{H}}_{1,Q_y}=({\mathbf{R}}_s-{\mathsf{\Theta }}^{*}{\mathsf{\Phi }}_y^{1-Q_y}{\mathbf{R}}_s^{*}{\mathsf{\Phi }}_y^{Q_y-1}\mathsf{\Theta })$.

Combining (9) and (11), the new differencing matrix constructed by the data information below the diagonal line can be rewritten as

$$\begin{array}{cc}\hfill {\overline{\mathbf{R}}}_1& =[{\overline{\mathbf{R}}}_{1,1},{\overline{\mathbf{R}}}_{1,2},\cdots ,{\overline{\mathbf{R}}}_{1,Q_y-1},{\mathbf{R}}_{1,Q_y}^d]\hfill \\ \hfill & ={\mathbf{A}}_P\left\{{\overline{\mathbf{H}}}_{1,1},{\overline{\mathbf{H}}}_{1,2},\cdots ,{\overline{\mathbf{H}}}_{1,Q_y-1},{\mathbf{H}}_{1,Q_y}\right\}diag\left({\mathbf{A}}_{P_x}^H,\cdots ,{\mathbf{A}}_{P_x}^H,{\mathbf{A}}_P^H\right),\hfill \end{array}$$

(12)

where the number of ${\mathbf{A}}_{P_x}^H$ is $\left((Q_y+1)Q_y/2-1\right)$. Similar with the processing of the first row rectangular subarray, we can extract the information of the $q_x$-th subarray and form the corresponding differencing matrix ${\mathbf{R}}_{q_x}$. As a result, the FO-ISD matrix can be defined as

$$\begin{array}{c}\hfill {\mathbf{R}}^f=\frac{1}{Q_x}{\sum }_{q_x=1}^{Q_x}{\mathbf{R}}_{q_x}.\end{array}$$

(13)

where

$$\begin{array}{c}\hfill {\overline{\mathbf{R}}}_{q_x}={\mathbf{A}}_P\left\{{\overline{\mathbf{H}}}_{q_x,1},{\overline{\mathbf{H}}}_{q_x,2},\cdots ,{\overline{\mathbf{H}}}_{q_x,Q_y-1},{\mathbf{H}}_{q_x,Q_y}\right\}diag\left({\mathbf{A}}_{P_x}^H,\cdots ,{\mathbf{A}}_{P_x}^H,{\mathbf{A}}_P^H\right),\end{array}$$

(14)

and ${\overline{\mathbf{H}}}_{q_x,n}=\left\{\left({\mathbf{R}}_s-{\mathsf{\Phi }}_x^{1-q_x}{\mathsf{\Theta }}^{*}{\mathsf{\Phi }}_y^{1-n}{\mathbf{R}}_s^{*}{\mathsf{\Phi }}_y^{n-1}{\mathsf{\Phi }}_x^{P_x-1}{\mathsf{\Phi }}_x^{q_x-1}\right),\cdots ,{\mathsf{\Phi }}_x^{q_x-1}{\mathsf{\Phi }}_y^{Q_y-1}{\mathbf{R}}_s{\mathsf{\Phi }}_y^{1-n}{\mathsf{\Phi }}_x^{1-q_x}\right\}$, ${\mathbf{H}}_{q_x,Q_y}$$=({\mathbf{R}}_s-{\mathsf{\Phi }}_x^{1-q_x}{\mathsf{\Theta }}^{*}{\mathsf{\Phi }}_y^{1-Q_y}{\mathbf{R}}_s^{*}{\mathsf{\Phi }}_y^{Q_y-1}\mathsf{\Theta }{\mathsf{\Phi }}_x^{q_x-1})$

Based on the definition in (13), we can prove that the FO-ISD matrix ${\mathbf{R}}^f$ has the following property.

Theorem 1. 

Consider a URA consisting of $M\times N$ sensors, where both the x and y directions are separated by half a wavelength. By performing partial difference operation for each row rectangular subarray, we can form the FO-ISD matrix ${\mathbf{R}}^f$ as in (13). Then, if $P_x{P}_y>K$, $Q_x{Q}_y>K$, $Q_x,Q_y>1$, the rank of ${\mathbf{R}}^f$ is equal to the number of coherent signals.

Proof. 

See the Appendix A. ☐

3.2.3. Forward Backward ISD (FB-ISD) Method

In this part, using the FB processing, we develop the FB-ISD method as follows.

As in (14), we can calculate that [20]

$$\begin{array}{c}{\mathbf{J}}_{((Q_y+1)Q_y/2-1)P_x+P_x{P}_y}{\mathbf{A}}_1={\mathbf{J}}_{((Q_y+1)Q_y/2-1)P_x+P_x{P}_y}diag({\mathbf{A}}_{P_x}^{*},\cdots ,{\mathbf{A}}_{P_x}^{*},{\mathbf{A}}_P^{*})\hfill \\ =diag\left({\mathbf{J}}_{P_x}{\mathbf{A}}_{P_x}^{*},\cdots ,{\mathbf{J}}_{P_x}{\mathbf{A}}_{P_x}^{*},{\mathbf{J}}_{P_x{P}_y}{\mathbf{A}}_P^{*}\right)=diag\left({\mathbf{A}}_{P_x}{\mathsf{\Phi }}_x^{1-P_x},\cdots ,{\mathbf{A}}_{P_x}{\mathsf{\Phi }}_x^{1-P_x},{\mathbf{A}}_P{\mathsf{\Theta }}^{*}\right)\hfill \\ ={\mathbf{A}}_{1}diag\left({\mathsf{\Phi }}_x^{1-P_x},\cdots ,{\mathsf{\Phi }}_x^{1-P_x},{\mathsf{\Theta }}^{*}\right),\hfill \end{array}$$

(15)

where ${\mathbf{A}}_1=diag\left({\mathbf{A}}_{P_x}^{*},\cdots ,{\mathbf{A}}_{P_x}^{*},{\mathbf{A}}_P^{*}\right)$. Combining (14) and (15), we can get the backward spatial differencing matrix of ${\mathbf{R}}_{q_x}$ as

$$\begin{array}{cc}\hfill {\mathbf{R}}_{q_x}^b& ={\mathbf{J}}_{P_x{P}_y}{\mathbf{R}}_{q_x}^{*}{\mathbf{J}}_{((Q_y+1)Q_y/2-1)P_x+P_x{P}_y}\hfill \\ \hfill & ={\mathbf{A}}_P\mathsf{\Theta }{\left\{{\overline{\mathbf{H}}}_{q_x,1},{\overline{\mathbf{H}}}_{q_x,2},\cdots ,{\overline{\mathbf{H}}}_{q_x,Q_y-1},{\mathbf{H}}_{q_x,Q_y}\right\}}^{*}diag\left({\mathsf{\Phi }}_x^{P_x-1},\cdots ,{\mathsf{\Phi }}_x^{P_x-1},\mathsf{\Theta }\right){\mathbf{A}}_1^H.\hfill \end{array}$$

(16)

Then, the FB-ISD matrix can be defined as

$$\begin{array}{c}\hfill {\mathbf{R}}^{fb}=\frac{1}{2Q_x}{\sum }_{q_x=1}^{Q_x}\left({\mathbf{R}}_{q_x}+{\mathbf{J}}_{P_x{P}_y}{\mathbf{R}}_{q_x}^{*}{\mathbf{J}}_{((Q_y+1)Q_y/2-1)P_x+P_x{P}_y}\right).\end{array}$$

(17)

Remark 1. 

Comparing FO-ISD and FB-ISD, ${\mathbf{R}}_{q_x}$ can extract all the information below the diagonal line, while ${\mathbf{R}}_{q_x}^{*}$ includes the information above the diagonal line. Thus, FB-ISD use more data information than FO-ISD. Besides, due to the increased number of smoothing sub-matrices, FB-ISD can achieve a further improved performance than FO-ISD, including accuracy and resolution.

3.2.4. Summary of FO-ISD and FB-ISD Methods

From (13) and (17), both FO-ISD and FB-ISD are proposed using the partial spatial differencing process for each row rectangular subarray. In fact, as described in Section 3.2.1, to improve the information utilization of the sample covariance matrix, we can first divide the URA into $Q_y$ column rectangular subarrays. Then, using the ISD scheme in Section 3.2.2 and Section 3.2.3, we can develop the corresponding FO-ISD matrix ${\mathbf{R}}_c^f$ and FB-ISD matrix ${\mathbf{R}}_c^{fb}$. In this case, the final ISD matrices can be set as

$$\begin{array}{c}\hfill {\mathbf{R}}_0^f=\left\{{\mathbf{R}}^f,{\mathbf{R}}_c^f\right\},{\mathbf{R}}_0^{fb}=\left\{{\mathbf{R}}^{fb},{\mathbf{R}}_c^{fb}\right\}.\end{array}$$

(18)

Clearly the final ISD matrices are more effective than ${\mathbf{R}}^f$ and ${\mathbf{R}}^{fb}$, due to the use of more data information. Then, the summary of proposed methods can be described as Algorithm 1.

Algorithm 1: FO-ISD and FB-ISD for 2-D DOA estimation of coherent signals

Input   : x$\left(t\right)$ = As$\left(t\right)$ + z$\left(t\right)$, t = 1, 2, …, L       Output: Estimated 2-D DOAs $\left({\theta }_k,{\alpha }_k\right),k=1,2,\dots ,K$ 1  Calculate the sample covariance matrix of ${\mathbf{x}}_{q_{x}x}\left(t\right)$ as ${\widehat{\mathbf{R}}}_{q_{x}x}=1/L{\sum }_{t=1}^L{\mathbf{x}}_{q_{x}x}\left(t\right){\mathbf{x}}_{q_{x}x}^H\left(t\right)$; 2  Construct the $q_x$-th spatial differencing matrix by performing partial difference-operation on the matrix ${\widehat{\mathbf{R}}}_{q_{x}x}$ as ${\widehat{\mathbf{R}}}_{q_x}$; 3  Extract the information above diagonal submatrices as ${\overline{\mathbf{R}}}_1=[{\overline{\mathbf{R}}}_{1,N},{\overline{\mathbf{R}}}_{1,N-1},\cdots ,{\overline{\mathbf{R}}}_{1,N-P_y+1},{\overline{\mathbf{R}}}_{1,1}]$ and compute its backward smoothing matrix ${\overline{\mathbf{R}}}_1^b$; 4  Build the FO-ISD matrix as ${\widehat{\mathbf{R}}}^f=1/Q_x{\sum }_{q_x=1}^{Q_x}{\widehat{\mathbf{R}}}_{q_x}$; 5  Compute the FB-ISD matrix as ${\widehat{\mathbf{R}}}^{fb}=1/\left(2Q_x\right){\sum }_{q_x=1}^{Q_x}\left({\widehat{\mathbf{R}}}_{q_x}+{\mathbf{J}}_{P_x{P}_y}{\widehat{\mathbf{R}}}_{q_x}^{*}{\mathbf{J}}_{((Q_y+1)Q_y/2-1)P_x+P_x{P}_y}\right)$; 6  Similar with Step 1∼6, for the column rectangular subarrays, construct the column spatial differencing matrix as ${\widehat{\mathbf{R}}}_c^f$ and ${\widehat{\mathbf{R}}}_c^{fb}$; 7  Compute the final ISD matrices as ${\widehat{\mathbf{R}}}_0^f=\left\{{\widehat{\mathbf{R}}}^f,{\widehat{\mathbf{R}}}_c^f\right\},{\widehat{\mathbf{R}}}_0^{fb}=\left\{{\widehat{\mathbf{R}}}^{fb},{\widehat{\mathbf{R}}}_c^{fb}\right\}$; 8  Perform singular value decomposition (SVD) operation on the final ISD matrices and use 2-D ESPRIT algorithm for 2-D DOA estimation [28].

Remark 2. 

AF-SDS and AFB-SDS only use the data information of overlapped smoothing subarrays, while the proposed scheme can extract all the data information of each row or column rectangle subarrays. Besides, both FO-ISD and FB-ISD only perform the difference-operation on the auto-correlations and the cross-correlations are kept unchanged. In this case, both the differencing submatrices and cross-correlations can be used to suppress the effect of additive noise. Therefore, the ISD scheme can achieve a performance improvement due to less information loss.

Remark 3. 

FO-ISD and FB-ISD almost have the same computational complexity, which includes the formation of ISD matrices, SVD operation, and eigenvalue decomposition (EVD) operation. To avoid the increase of computational complexity, we can form the ISD matrices by using the sample covariance matrix of received signals, the cost of which is about $L{M}^2{N}^2+{\left(P_x{P}_y\right)}^2\left(P_x{P}_y+Q_y\left(Q_y+1\right)/2\right)+2K^3$. Then, similar to the proposed methods, we can also get the covariance matrix of the sliding windows for the spatial smoothing technique [17] or A-SDS method [27] by the same way, the cost of which is about $L{M}^2{N}^2+{\left(P_x{P}_y\right)}^3+2K^3$. So the proposed methods can achieve performance improvement with slightly higher computations. To be clear, we show the runtime of relevant methods in Figure 3 . We can see that FB-ISD has a heavier runtime load than that of FBSS and A-ISD.

Remark 4. 

As described in Section 3.2, the proposed methods are developed for ideal sensors without any mutual coupling. Then, just like the methods in [29,30], the ISD scheme is also suitable for direction finding with unknown mutual coupling. However, when considering the array imperfections, the ISD scheme will fail due to the breakdown of array response matrices.

3.3. Cramér-Rao Bound (CRB)

As described in Section 2, according to [31], the CRB can be obtained as

$$\begin{array}{c}\hfill CRB=\frac{{\sigma }^2}{2L}{\left\{\mathrm{Re}\left[{\mathbf{D}}^H{\mathsf{\Pi }}_{\mathbf{A}}^{\perp }\mathbf{D}\oplus {\widehat{\mathbf{R}}}_s^T\right]\right\}}^{-1},\end{array}$$

(19)

where $\mathbf{D}=\left[\frac{\partial {\mathbf{a}}_1}{\partial {\theta }_1},\cdots ,\frac{\partial {\mathbf{a}}_K}{\partial {\theta }_K},\frac{\partial {\mathbf{a}}_1}{\partial {\varphi }_1},\cdots ,\frac{\partial {\mathbf{a}}_K}{\partial {\varphi }_K}\right]$, ${\widehat{\mathbf{R}}}_s=\left[\begin{array}{cc}{\mathbf{R}}_s& {\mathbf{R}}_s\\ {\mathbf{R}}_s& {\mathbf{R}}_s\end{array}\right]$, ${\mathsf{\Pi }}_{\mathbf{A}}^{\perp }={\mathbf{I}}_{MN}-\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{\mathbf{A}}^H$, ${\mathbf{a}}_k$ is the kth column of $\mathbf{A}$, k = 1, ⋯, K.

4. Simulation Results

We now evaluate the estimation performance of the ISD methods through in-depth numerical experiments. We assume the number of sensors is $M=N=9$. The wavelength of transmitted signals is set as 1m and the estimation performance is examined over 500 Monte Carlo trials. To evaluate the performance, the root-mean-square-error (RMSE) can be defined as

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{1}{2\overline{N}K}\sum _{i=1}^{\overline{N}}\sum _{k=1}^K\left[{\left({\alpha }_k-{\widehat{\alpha }}_k^{\left(i\right)}\right)}^2+{\left({\theta }_k-{\widehat{\theta }}_k^{\left(i\right)}\right)}^2\right]}.\end{array}$$

(20)

where $\overline{N}$ denotes the total independent trials and (${\alpha }_k$, ${\theta }_k$) and (${\widehat{\alpha }}_k$, ${\widehat{\theta }}_k$) represent the true and estimated 2-D DOAs of the k-th signals, respectively.

4.1. Effectiveness Evaluation

In this experiment, we examine the effectiveness of the FB-ISD method for Gaussian white noise and colored noise, respectively. The colored noise is of a second-order autoregressive (AR) model with coefficients $\mathbf{a}=[1,-0.7,0.6]$[11,12,13,14,15,16]. The size of the subarrays are $P_x=P_y=6$ and the signals are located at $\mathbf{\alpha }=[{10}^{°},{20}^{°},{30}^{°},{40}^{°}]$, $\mathbf{\theta }=[{20}^{°},{30}^{°},{40}^{°},{50}^{°}]$. The number of snapshots is L = 300 and the signal to noise ratio (SNR) is 15 dB. Figure 4 shows the estimation results of FB-ISD with 200 Monte Carlo trials in white and colored noise conditions, respectively. As expected, all the 2-D DOAs can be estimated effectively and accurately. Besides, FB-ISD with white noise performs better than that of colored noise, especially for the first two signals.

Then, we compare the resolution ability by resolving two closely located signals in white and colored noise conditions, where the signals are located at $\mathbf{\alpha }=[{10}^{°},{12}^{°}]$, $\mathbf{\theta }=[{20}^{°},{22}^{°}]$, the SNR is 10 dB, and the number of snapshots is 200. Figure 5 describes the estimated DOAs of FB-ISD and A-SDS with 200 Monte Carlo trials. We can observe that FB-ISD can resolve the two signals successfully by using the difference-operation, while A-SDS fails to distinguish the signals.

4.2. RMSE Performance in the Case of White Noise

In this experiment, we examine the performance of the proposed methods versus SNR and the number of snapshots in the white noise condition. Here we assume the signal locations as $\mathbf{\alpha }=[{10}^{°},{20}^{°},{30}^{°}]$, $\mathbf{\theta }=[{20}^{°},{30}^{°},{40}^{°}]$, and we compare the proposed methods with other existing methods, including the forward only spatial smoothing (FOSS) method, FBSS method [16], and A-SDS method [27] ($p_a=3$). Moreover, the CRB is provided for comparison.

Performance versus SNR: Figure 6 shows the RMSE versus SNR in the white noise condition, where we assume $P_x=P_y=7$ and $L=300$. It is observed that, the performance of FB-ISD is better than that of other methods due to the full use of data information. Since the partial difference-operation has less information loss than that of A-SDS, the performance of FB-ISD is much better, especially in the low SNR condition. Then, FB-ISD performs a little better than FO-ISD due to the FB processing. We also see that the RMSE curve of FB-ISD is very close to the CRB. In addition, combining with Figure 3, we can conclude that the proposed methods have a better performance but also higher computations.

Performance versus the number of snapshots: Here, we evaluate the performance in terms of the number of snapshots, where the SNR is 0 dB and $P_x=P_y=7$. As shown in Figure 7, when the number of snapshots is small, FB-SMS still outperforms other methods and also matches closely to the CRB due to the use of more data information. Then, A-SDS performs much worse for the small number of snapshots, and the reason is that the information loss caused by difference-operation becomes greater with the decrease of snapshots. Besides, just like Figure 6, since spatial smoothing based methods suffer from aperture loss, all these methods cannot converge to the CRB in the limit for very high SNR condition.

4.3. RMSE Performance in Case of Colored Noise

In this experiment, we focus on the performance analysis of the proposed methods in a colored noise condition, where the performance for the white noise and colored noise conditions is also compared systematically.

Performance versus SNR: Figure 8 presents the RMSE curves in terms of SNR in the colored noise condition, where the parameters are the same as Figure 7 except for the SNR varying from 0 dB to 35 dB. From Figure 8, we can see that the performance of FB-ISD is better than those of methods in [16,27]. Then, in the colored noise condition, FB-ISD and FO-ISD perform much better than others due to the partial difference-operation. Compared with Figure 6, we can summarize that the proposed scheme can suppress the effect of colored noise more effectively, and the reason is that the colored noise covariance matrix has significant values for diagonal and non-diagonal elements. Then, compared with A-SDS, we can also reach the following conclusion: FB-ISD can improve the performance greatly by only by performing the difference-operation on auto-correlations, for both colored noise and white noise.

Performance versus the number of snapshots: Figure 9 shows the RMSE against the number of snapshots in the colored noise condition, where the parameters are the same as Figure 7. As described in Figure 9, the performance of FB-ISD is superior to the other methods. Then, compared to Figure 7, the major distinction is that the performance for these methods is not much affected by the number of snapshots in the colored noise condition. The reason is that the information loss caused by difference-operation can also change with the number of snapshots. Besides, the proposed methods still have better performance, especially for the small number of snapshots.

5. Conclusions

In this paper, we have proposed an ISD scheme for 2-D DOA estimation of coherent signals with URAs, including FO-ISD and FB-ISD methods. By extracting all the data information of each row or column rectangular subarrays, we only performed the difference-operation on the auto-correlations, while the cross-correlations were kept unchanged. In this case, the reconstructed submatrices can use more data information of the sample covariance matrix and also suppress the effect of additive noise more effectively. Then, both the FO-ISD and FB-ISD methods were developed using the spatial smoothing submatrices. Simulation results demonstrated that, compared with other recently spatial smoothing and spatial differencing techniques, the performance of proposed methods was superior in white or colored noise conditions, in terms of accuracy and resolution ability.