Degree-of-Freedom Strengthened Cascade Array for DOD-DOA Estimation in MIMO Array Systems

Abstract

In spatial spectrum estimation, difference co-array can provide extra degrees-of-freedom (DOFs) for promoting parameter identifiability and parameter estimation accuracy. For the sake of acquiring as more DOFs as possible with a given number of physical sensors, we herein design a novel sensor array geometry named cascade array. This structure is generated by systematically connecting a uniform linear array (ULA) and a non-uniform linear array, and can provide more DOFs than some exist array structures but less than the upper-bound indicated by minimum redundant array (MRA). We further apply this cascade array into multiple input multiple output (MIMO) array systems, and propose a novel joint direction of departure (DOD) and direction of arrival (DOA) estimation algorithm, which is based on a reduced-dimensional weighted subspace fitting technique. The algorithm is angle auto-paired and computationally efficient. Theoretical analysis and numerical simulations prove the advantages and effectiveness of the proposed array structure and the related algorithm.

1. Introduction

In the past few decades, research on the theory and algorithms of estimating direction of arrival (DOA) in array processing usually focused on the overdetermined system model with uniform linear array (ULA) configuration. A basic piece of knowledge is that the number of sources that can be resolved by the MUSIC algorithm is $N-1$ if adopting a N-element ULA [1,2]. However, many scenarios imply an underdetermined system model in which the number of unknown sources may be actually more than the number of physical sensors. To solve such a challenging problem, most techniques are based on creating the effect of a equivalent array with many more virtual elements than actual array through suitable spatial-domain, time-domain or frequency-domain sampling strategies.

The earliest pioneering studies, by exploiting the possible element arrangements, such as space tapered arrays and fractal-based quasi random array [3,4,5], are designed for yielding reasonably low sidelobe level; and, recently, the random array has been integrated into a hardware prototype to achieve sub-Nyquist sampling in spatial and spectral domains [6]. One of the most popular virtual arrays is that of synthetic aperture radar (SAR), where the larger artificial aperture is created by means of the motion of radar antennas [7]. In wide-band array signal processing, all sub-band signals are usually aligned into one single frequency to antagonize the coherence of signals, which can be viewed as a special virtual array that takes advantage of the frequency diversity [8,9]. In recent years, one of the most popular strategies that can provide the ability of resolving more sources than sensors is the so-called co-array. For example, in the work of [10], the quasi-stationary signals are considered to construct a virtual ULA for improving the parameter identifiability. Basically, the concept of co-array can be categorized into sum co-array and difference co-array. The first one usually appears in some active sensing systems such as co-located MIMO radar [11], whereas the second one appears in the passive one such as traditional direction finding system or bistatic MIMO radar [12,13,14,15,16,17,18,19].

Generally speaking, for an N-element array, both the sum co-array and difference co-array can create up to $\mathcal{O}(N^2)$ elements, depending on the geometry of the physical sensor array. Theoretically, minimum redundancy array (MRA) [20] is the optimum linear geometry, which can create the maximum possible number of elements in difference co-array without any missing elements in between, but retains some unwilling repeated elements. However, the MRA with a given N has no closed-form expression either for its structure or the achievable degrees-of-freedom (DOFs), and hence its geometry is not analytically tractable. By combining two or more ULAs with increasing inter-sensor spacing, Ref. [13] designs a new array structure named a nested array, which is capable of providing a dramatic increase in DOFs and hence can resolve more sources than the actual number of sensors. Later, the coprime sensing array [14] and its generalized versions [15] are designed for efficient sampling of multidimensional signals, which can be directly applied to the aforementioned underdetermined model.

Although the above array geometries/structures all can provide $\mathcal{O}(N^2)$ DOFs, there still exist quantitative differences. Furthermore, such difference determines the capability of resolving more sources than sensors. By comparison, in the sense of difference co-array with a given number of physical sensors, MRA sets an upper-bound of DOFs; the DOFs provided by nested array [13] is greater than the one of the coprime arrays [14], and both configurations are inferior to the coprime array with displaced subarrays (CADiS) [15]. However, the CADiS cannot stay consecutive in the whole co-array because it has holes. Actually, promoting the potential DOFs closely to the upper-bound is still an open challenge problem, the kernel of which relies on a more effective design for sensors’ positions.

In addition, when we apply the concept of difference co-array into the multiple input multiple output (MIMO) array systems such as MIMO radar or some wireless communication systems, the other key point is the parameter estimation algorithm with respect to direction-of-arrival (DOA) or carrier frequency. For example, Refs. [13,14,21,22,23,24] adopt multiple signal classification (MUSIC) spectrum searching to estimate DOA information, and Refs. [25,26,27,28,29] discuss algorithms with super resolution from the perspective of sparsity-based recovery. Unfortunately, however, the above algorithms or tensor-based algorithms [19], Refs. [30,31] will bring huge computational complexity in practice if more than one-dimensional parameters are considered. There exist non-spectrum-searching choices, e.g., the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [32] and MUSIC-rooting algorithm, which can decrease the calculation amount to some extent, but they are usually at the cost of reducing angle estimation accuracy.

In this paper, motivated by the need for efficiently solving the aforementioned problems, a series of work is developed to improve the exist weaknesses. The main contributions are three-fold:

• A new array structure named cascade array is designed, which is essentially composed of one ULA and one non-uniform linear array. Through theoretical analysis, the difference co-array of the designed optimal cascade array is hole-free and can provide more DOFs than some state-of-the-art sensor array structures. That is to say, it manifests a strengthened resolving capability.
• We then apply the cascade array into bistatic MIMO array systems to achieve joint direction-of-departure (DOD) and DOA estimation for multiple targets localization. Furthermore, by parameterizing the orthogonal projector onto the null space of the equivalent joint steering matrix, a novel algorithm based on a reduced-dimensional weighted subspace fitting technique is proposed.
• The proposed algorithm transforms a two-dimensional estimation problem into a one-dimensional one. To do so, the DOD information can be acquired by MODE rooting [33,34], and the auto-paired DOA information can be extracted from the estimated receiving array manifold. The proposed algorithm avoids exhausted spectrum searching or tensor decomposition [19,34,35], thus it is computationally efficient.

The rest of this paper is organized as follows: in Section 2, we introduce the cascade array and its optimization. The application in MIMO array systems and related two-dimensional difference co-array are presented in Section 3. In Section 4, we provide an effective joint DOD and DOA estimation algorithm. Numerical simulations are shown in Section 5, followed by the conclusions in Section 6.

Notation: ${(·)}^{*}$, ${(·)}^T$, ${(·)}^H$, ${(·)}^{†}$ denote the complex conjugate, transpose, Hermitian transpose, pseudo-inverse, respectively. Symbol ⊗ represents a Kronecker product and symbol ⊙ is Khatri–Rao product, which is a column-wise Kronecker product. ${\mathbf{I}}_M$ is a $M\times M$ identity matrix and $\mathbf{0}$ symbolizes zero matrix. ${\mathbf{A}}^{(m)}$ is a submatrix of $\mathbf{A}$ formed by its last m rows. Operator $\angle [·]$ serves to get the phase. In addition, as a shorthand notation, the addition between a set and a scalar is defined as $\mathcal{S}\pm b=\{a\pm b|\forall a\in \mathcal{S}\}$ and the difference set between ${\mathcal{S}}_1$ and ${\mathcal{S}}_2$ is given by $\mathrm{Diff}({\mathcal{S}}_1,{\mathcal{S}}_2)=\{a-b|a\in {\mathcal{S}}_1,b\in {\mathcal{S}}_1\}$.

2. Cascade Array

In this section, we will introduce the detailed structure of our designed cascade array, and then make a DOF comparison with other types of array configuration. For convenience, we first introduce the concept of difference co-array and its closed relationship with the Khatri–Rao product of two steering vectors [10] and [13], which plays an important role in the follow-up content.

2.1. Difference Co-Array

Given a sensor position set $\mathcal{P}=\{p_i,1\le i\le N\}$, the difference co-array with respect to set $\mathcal{P}$ can be created by

$${\mathcal{P}}_d=\{\overline{p}|\overline{p}=p_i-p_k,p_i\in \mathcal{P},p_k\in \mathcal{P}\},$$

(1)

which means that each element in ${\mathcal{P}}_d$ represents a virtual sensor.

It is worth mentioning that there exist some repeated values and/or ‘holes’ in this new set, which is completely determined by the array geometry. The number of DOFs with respect to a real set $\mathcal{P}$ is the cardinality of its difference co-array set. If defining $|{\mathcal{P}}_d|$ be the cardinality of set ${\mathcal{P}}_d$, i.e., the number of distinct elements, as we know, $|{\mathcal{P}}_d|=2N-1$ for a ULA, whereas that for MRA, nested array, coprime array, CADiS, the coprime array with compressed inter-element spacing (CACIS) [15] is $\mathcal{O}(N^2)$.

To establish a relationship between difference co-array and sensor array manifold $\mathbf{C}(\mathbf{\theta })\in {\mathbb{C}}^{N\times K}$, let us consider the Khatri–Rao product, $\tilde{\mathbf{C}}(\mathbf{\theta })={\mathbf{C}}^{*}(\mathbf{\theta })\odot \mathbf{C}(\mathbf{\theta })$. If defining p-th column of $\mathbf{C}(\mathbf{\theta })$ be $\mathbf{c}({\theta }_p)$, then correspondingly the p-th column of $\tilde{\mathbf{C}}(\mathbf{\theta })$ is given by

$$\tilde{\mathbf{c}}({\theta }_p)={\mathbf{c}}^{*}({\theta }_p)\otimes \mathbf{c}({\theta }_p)={[c_1^{*}{\mathbf{c}}^{\mathrm{T}}({\theta }_p),c_2^{*}{\mathbf{c}}^{\mathrm{T}}({\theta }_p),\cdots ,c_N^{*}{\mathbf{c}}^{\mathrm{T}}({\theta }_p)]}^{\mathrm{T}},$$

(2)

where $c_i$ is the i-th element of vector $\mathbf{c}({\theta }_p)$. It can be seen that there is a direct connection between $\tilde{\mathbf{c}}({\theta }_p)$ and ${\mathcal{P}}_d$ since the elements of $\tilde{\mathbf{c}}({\theta }_p)$ are given by

$${[\tilde{\mathbf{c}}]}_l=e^{j\pi (p_i-p_k)sin\theta },\hspace{1em}l=(m-1)N+n.$$

(3)

Obviously, each row of $\tilde{\mathbf{c}}({\theta }_p)$ has a one-to-one correspondence with the element in ${\mathcal{P}}_d$, and therefore $\tilde{\mathbf{C}}(\mathbf{\theta })$ implies the manifold of a virtual array whose sensors are located at the positions indicated by set ${\mathcal{P}}_d$.

A more effective array geometry means that it should yield as many distinct elements as possible in concrete natural numbers, which characterizes larger DOFs or a strong capability of resolving more sources than sensors; and simultaneously it should have no vacancy or holes, which characterizes a virtual ULA structure and also a free usage of traditional angle estimation algorithms. Both requirements will also be a criterion in the following designing and comparison.

2.2. Cascade Array Design

The kernel of designing array geometry is to determine and optimize the array elements’ positions. As we know, the nested array or coprime array is still rooted in uniform linear geometry with different element-spacing. Taking a two-level nested array [13], for example, as the array structure is shown in Figure 1, the second level ULA has a spacing that is determined by the number of the first level ULA. The biggest advantages of such structure are of easily designing and of easily figuring out the specific DOFs. However, the pursuit of more DOFs cannot be confined by that uniform feature, e.g., MRA, which is a typical example that can achieve the largest DOFs through non-uniform features. The cascade array in this research essentially borrows the non-uniform feature and meanwhile inherits the feature of nested array.

Assuming the total number of sensors is $N=N_1+N_2$ and a basic spacing unit d is set as half-wavelength, i.e., $d=\lambda /2$, where $\lambda$ is the carrier wavelength. Our cascade array can be partitioned into two parts. Just like the nested array, the first level array is still a $N_1$-element ULA with one spacing unit, but, on the contrary, the second level array is a $N_2$-element non-uniform linear array.

Compared with nested array, there are mainly two modifications. On the one hand, we enlarge the spacing between two arrays up to $(N_1+1)$ spacing unit, the purpose of which is to strengthen the DOF, but it will incur holes; therefore, on the other hand, we further shorten the spacing of the last two sensors in the second level array, which can mend that defect. In a word, the designed cascade array is described by a precise definition, i.e.,

Definition 1.

Assume $N_1$ and $N_2$ are integers satisfying $N_1\ge 1$ and $N_2\ge 2$. Second-order cascade array is specified by the integers set $\mathcal{L}$, defined by

$$\mathcal{L}={\mathcal{L}}^1\cup {\mathcal{L}}^2\cup {\mathcal{L}}^3,$$

(4)

where each subset is given as below:

$$\begin{array}{ccc}\hfill {\mathcal{L}}^1& =& \{0,1,2,\cdots ,N_1-1\},\hfill \\ \hfill {\mathcal{L}}^2& =& \{0,(N_1+1),\cdots ,(N_2-2)(N_1+1)\}+2N_1,\hfill \\ \hfill {\mathcal{L}}^3& =& \{(N_2-2)(N_1+1)+N_1\}+2N_1,\hfill \end{array}$$

And its difference co-array set is

$${\mathcal{P}}_d^{\mathcal{L}}=\mathrm{Diff}(\mathcal{L},\mathcal{L})=\{\overline{p}|\overline{p}=p_i-p_k,p_i\in \mathcal{L},p_k\in \mathcal{L}\}.$$

(5)

In set $\mathcal{L}$, ${\mathcal{L}}^1$ denotes the ULA and ${\mathcal{L}}^2\cup {\mathcal{L}}^3$ denotes the non-uniform one. The typical configurations for different array geometry with total $N=7$ sensors is shown in Figure 1. In this example, we make a graphical illustration for the designed cascade array, ULA, nested array [13], super nested array [30], CACIS [15], coprime array [14], nested CADiS (it is an optimum solution of CADiS under the requirement for highest number of consecutive lags) [15] and MR array [20]. We can see that the designed cascade array can generate more DOFs than other six array geometries and attain the same performance as MR array.

To completely exploit the DOFs of the cascade array, we summarize the properties of ${\mathcal{P}}_d^{\mathcal{L}}$ in the following theorem.

Theorem 1.

The second-order cascade array in Definition 1 is a kind of restricted array [13], i.e, its difference co-array is hole-free; and it can provide DOFs with the number of $2N_2(N_1+1)+2N_1-3$.

Proof. 

The strict and detailed derivations are shown in Appendix A. ☐

According to the above theorem, we can further reach a new corollary with respect to the problem that optimizes the distribution of sensors in the configuration of cascade array by finding $N_1$ and $N_2$ to maximize the total DOFs under a fixed number of sensors.

Corollary 1.

Given N physical sensors, e.g., $N=N_1+N_2$, the optimal choices of $N_1$ for the 1st level array and $N_2$ for the 2nd level array in the designed cascade array configuration are verified as

N	Optimal ${\mathit{N}}_{\mathbf{1}}$ and ${\mathit{N}}_{\mathbf{2}}$	The Number of DOFs
Even	$N_1=N_2=\frac{N}{2}$	$\frac{N^2}{2}+2N-3$
Odd	$N_1=\frac{N-1}{2},N_2=\frac{N+1}{2}$	$\frac{N^2-1}{2}+2N-3$

Proof. 

The solution is easily obtained by utilizing arithmetic mean-geometric mean (AM-GM) inequality, and we skip it for simplicity. ☐

The detailed comparison with respect to the DOFs, parameter identifiability and estimation accuracy of the designed cascade array will be introduced in Section 5.

3. MIMO Array Systems with Cascade Array

In this section, we will apply the aforementioned results into MIMO array systems to achieve higher accuracy of parameters estimation, and mainly focus on the problem of targets localization, i.e., the joint direction of departure (DOD) and direction of arrival (DOA) estimation.

3.1. Data Model

We herein consider a MIMO array systems with cascade array configuration at both transmitting and receiving end, in which M and N antenna elements are arranged, respectively. The whole system is frequency synchronization or fully calibrated [36]. It is also assumed that there are K targets, and the output baseband signal of the matched filters at the receiving end can be written as [18,19]

$$\mathbf{y}(t)=[\mathbf{b}({\varphi }_1)\otimes \mathbf{a}({\theta }_1),\cdots ,\mathbf{b}({\varphi }_K)\otimes \mathbf{a}({\theta }_K)]\mathbf{h}(t)+\mathbf{w}(t),$$

(6)

where the transmitting steering vector $\mathbf{a}({\theta }_k)$ and the receiving steering vector $\mathbf{b}({\varphi }_k)$, for $k=1,2,\cdots ,K$, are assumed to be unchanged during a coherent processing interval (CPI), and ${\theta }_k,{\varphi }_k$ are the DOD and DOA of the k-th target, respectively. We herein adopt a Swerling I target model, so the signal vector $\mathbf{h}(t)={[{\gamma }_1(t),\cdots ,{\gamma }_K(t)]}^{\mathrm{T}}$ relies on the Doppler frequency $f_{dk}$ and the radar cross section (RCS) coefficient ${\beta }_k$, i.e., then ${\gamma }_k(t)={\beta }_k{e}^{j2\pi f_{dk}(t-1)}$. Note that $\mathbf{a}({\theta }_k)={[e^{j\pi l_{1}sin{\theta }_k},\cdots ,e^{j\pi l_{M}sin{\theta }_k}]}^{\mathrm{T}}$ and $\mathbf{b}({\varphi }_k)={[e^{j\pi p_{1}sin{\varphi }_k},\cdots ,e^{j\pi p_{N}sin{\varphi }_k}]}^{\mathrm{T}}$ are mutilated Vandermonde vectors, where ${\{l_m\}}_{m=1}^M$ and ${\{p_n\}}_{N=1}^N$ is the antennas’ position defined by the cascade array configuration. $\mathbf{w}(t)$ is the additive zero-mean Gaussian noise with covariance ${\sigma }_n^2$.

Defining $\mathbf{A}=[\mathbf{a}({\theta }_1),\cdots ,\mathbf{a}({\theta }_K)]\in {\mathbb{C}}^{M\times K}$, $\mathbf{B}=[\mathbf{b}({\varphi }_1),\cdots ,\mathbf{b}({\varphi }_K)]\in {\mathbb{C}}^{N\times K}$, and after collecting Q consecutive pulses, Label (6) can be rewritten by the following compact form:

$$\mathbf{Y}=(\mathbf{B}\odot \mathbf{A})\mathbf{S}+\mathbf{W},$$

(7)

where $\mathbf{Y}=[\mathbf{y}(1),\mathbf{y}(2),\cdots ,\mathbf{y}(Q)]$, $\mathbf{S}=[\mathbf{h}(1),\mathbf{h}(2),\cdots ,\mathbf{h}(Q)]$ and $\mathbf{W}=[\mathbf{w}(1),\mathbf{w}(2),\cdots ,\mathbf{w}(Q)]$.

3.2. Two-Dimensional Difference Co-Array

Assuming the Doppler frequencies of all K targets are all different and also well separated. The covariance matrix of the received data $\mathbf{y}(t)$ is given by

$$\mathbf{R}=\mathbb{E}[\mathbf{y}(t){\mathbf{y}}^{\mathrm{H}}(t)]=(\mathbf{B}\odot \mathbf{A})\mathbf{\Lambda }{(\mathbf{B}\odot \mathbf{A})}^{\mathrm{H}}+{\sigma }_n^2{\mathbf{I}}_{MN}.$$

(8)

Due to the K target reflected signals in ${\{{\gamma }_k(t)\}}_{k=1}^K$ are sinusoidal signals with different frequencies and attenuation coefficients, they are temporally uncorrelated. If sampling them by a sufficient large snapshot rate (at least twice the maximum Doppler frequency), then the signal autocorrelation matrix $\mathbf{\Lambda }$ is diagonal, and $\mathbf{\Lambda }=\mathrm{diag}[\mathbf{p}]$ with $\mathbf{p}={[{\beta }_1^2,\cdots ,{\beta }_K^2]}^{\mathrm{T}}$.

We vectorize $\mathbf{R}$ to get the following vector:

$$\mathbf{r}=vec(\mathbf{R})=[{(\mathbf{B}\odot \mathbf{A})}^{*}\odot (\mathbf{B}\odot \mathbf{A})]\mathbf{p}+{\sigma }_n^2\mathbf{1},$$

(9)

where $vec(·)$ denotes the vectorizing operator. $\mathbf{1}={[{\mathbf{e}}_1^{\mathrm{T}},\cdots ,{\mathbf{e}}_{MN}^{\mathrm{T}}]}^{\mathrm{T}}$ with ${\mathbf{e}}_i$ being a vector of all zeros except a 1 at the i-th position.

Compared with one-dimensional scenario, it is a little more complicated when achieving difference co-array under the two-dimensional scenario, such as Label (9). Hence, we introduce the following lemma:

Lemma 1.

Consider a auxiliary matrix defined by $\mathbf{\Pi }={\mathbf{I}}_N\otimes \mathbf{\Gamma }\otimes {\mathbf{I}}_M$, where

$$\mathbf{\Gamma }={\left[\begin{array}{c}\mathbf{\Delta }\\ {\mathbf{0}}_{M\times 1}& \mathbf{\Delta }\\ ⋮& ⋮& \ddots \\ {\mathbf{0}}_{M\times 1}& {\mathbf{0}}_{M\times 1}& \dots & \mathbf{\Delta }\end{array}\right]}^{MN\times MN}and\mathbf{\Delta }={\left[\begin{array}{c}1\\ {\mathbf{0}}_{1\times N}& 1\\ ⋮& ⋮& \ddots \\ {\mathbf{0}}_{1\times N}& {\mathbf{0}}_{1\times N}& \dots & 1\end{array}\right]}^{M\times [N(M-1)+1];}$$

then, the following result will hold:

$$\mathbf{\Pi }[{(\mathbf{B}\odot \mathbf{A})}^{*}\odot (\mathbf{B}\odot \mathbf{A})]=({\mathbf{B}}^{*}\odot \mathbf{B})\odot ({\mathbf{A}}^{*}\odot \mathbf{A}).$$

(10)

The above conclusion utilizes some properties of the Khatri–Rao product of multiple matrices, and one can refer to [37] for more detail.

According to the lemma above, after left-multiplying the permutation matrix $\mathbf{\Pi }$ on Label (9), we can acquire a new constructed data with the following form:

$$\overline{\mathbf{r}}=\mathbf{\Pi }\mathbf{r}=[({\mathbf{B}}^{*}\odot \mathbf{B})\odot ({\mathbf{A}}^{*}\odot \mathbf{A})]\mathbf{p}+{\sigma }_n^2\mathbf{\Pi }\mathbf{1}.$$

(11)

The vector $\overline{\mathbf{r}}$ can be viewed as one-pulse baseband observationof a virtual MIMO radar with equivalent transmitting steering matrix $({\mathbf{A}}^{*}\odot \mathbf{A})$ and equivalent receiving steering matrix $({\mathbf{B}}^{*}\odot \mathbf{B})$ in a deterministic noise environment.

We define $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$ as the virtual transmitting and the virtual receiving steering matrices after deleting the redundant items (It can be achieved by deleting the corresponding row observations in $\overline{\mathbf{r}}$). Therefore, Label (11) can be further rewritten as:

$$\tilde{\mathbf{r}}=[\tilde{\mathbf{B}}\odot \tilde{\mathbf{A}}]\mathbf{p}+{\sigma }_n^2\mathbf{e}.$$

(12)

The k-th column of $\tilde{\mathbf{B}}$ and $\tilde{\mathbf{A}}$, respectively, is represented as:

$$\begin{array}{ccc}\hfill \tilde{\mathbf{b}}({\varphi }_k)& =& {[e^{-j\pi \overline{N}sin{\varphi }_k},\cdots ,1,\cdots ,e^{j\pi \overline{N}sin{\varphi }_k}]}^{\mathrm{T}},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill \tilde{\mathbf{a}}({\theta }_k)& =& {[e^{-j\pi \overline{M}sin{\theta }_k},\cdots ,1,\cdots ,e^{j\pi \overline{M}sin{\theta }_k}]}^{\mathrm{T}},\hfill \end{array}$$

(14)

where $\overline{M}=l_M,\overline{N}=p_N$. $\mathbf{e}$ is a zero column vector, except a 1 in the middle.

3.3. Identifiability Analysis

Now, we have to discuss the parameter identifiability for data model Label (12), i.e., the maximal number of targets that can be uniquely identified without noise because it is a prerequisite for a well-posed parameter estimation. Generally speaking, the identifiability for a given data model is acting like a theoretical upper bound. We herein only consider the angle estimation algorithms based on signal/noise subspace technique.

Actually, most signal subspace based algorithms for parameter estimating require no rank deficiency in the observation matrix; however, we can see that the data $\tilde{\mathbf{r}}$ is rank one. Therefore, we have to make a spatial smoothing before signal subspace decomposition.

Define the following $(\overline{M}+1)(\overline{N}+1)\times (\overline{M}+1)(\overline{N}+1)$ selection operator

$${\Xi }_{n,m}=[{\mathbf{0}}_{(\overline{N}+1)\times (\overline{N}+1-n)}{\mathbf{I}}_{(\overline{N}+1)}{\mathbf{0}}_{(\overline{N}+1)\times (n-1)}]\otimes [{\mathbf{0}}_{(\overline{M}+1)\times (\overline{M}+1-m)}{\mathbf{I}}_{\overline{M}+1}{\mathbf{0}}_{(\overline{M}+1)\times (m-1)}],$$

(15)

where $1\le n\le \overline{N}+1$, $1\le m\le \overline{M}+1$. If we stack the selected observation vector $\tilde{\mathbf{r}}$ as follows, we have

$$\mathcal{R}(\tilde{\mathbf{r}})=[{\Xi }_{1,1}\tilde{\mathbf{r}}\cdots {\Xi }_{1,\overline{M}+1}\tilde{\mathbf{r}}{\Xi }_{2,1}\tilde{\mathbf{r}}\cdots {\Xi }_{2,\overline{M}+1}\tilde{\mathbf{r}}\cdots {\Xi }_{\overline{N}+1,\overline{M}+1}\tilde{\mathbf{r}}].$$

(16)

Then, the above two-dimensional virtual spatial smoothing operation yields

$$\overline{\mathbf{Y}}=\mathcal{R}(\tilde{\mathbf{r}})=[\mathcal{B}\odot \mathcal{A}]\overline{\mathbf{S}}+{\sigma }_n^2{\mathbf{I}}_{(\overline{M}+1)(\overline{N}+1)}=\overline{\mathbf{A}}\overline{\mathbf{S}}+{\sigma }_n^2{\mathbf{I}}_{(\overline{M}+1)(\overline{N}+1),}$$

(17)

where the equivalent signal matrix is

$$\overline{\mathbf{S}}=\mathbf{\Lambda }{[\mathcal{B}\odot \mathcal{A}]}^{\mathrm{H}}.$$

(18)

For convenience, we let $\mathcal{B}={\tilde{\mathbf{B}}}^{(\overline{N}+1)}$, $\mathcal{A}={\tilde{\mathbf{A}}}^{(\overline{M}+1)}$, and their k-th column is denoted by $\overline{\mathbf{b}}({\varphi }_k)$ and $\overline{\mathbf{a}}({\theta }_k)$, respectively.

Theorem 2.

For a MIMO array systems with cascade array configuration, see Label (7), if we perform two-dimensional difference co-array by Lemma 1 to get Label (12) and further to get Label (17), the parameter set $({\theta }_k,{\varphi }_k)$, $k=1,2,\cdots ,K$, can be uniquely identified if

$$\{\begin{array}{cc}K\le (\overline{N}+1)(\overline{M}+1)-1,& for 2D-MUSIC,\\ K\le \overline{N}(\overline{M}+1),& for 2D-ESPRIT,\end{array}$$

(19)

where the parameters’ set is drawn from a continuous distribution with respect to the Lebesgue measure in ${\mathbb{L}}^{2K}$, $\mathbb{L}:=[-\pi /2,\pi /2]$.

Proof. 

Suppose the parameters are all drawn from a continuous distribution. To view in retrospect the model described in Label (17) and illuminating this by the almost surely full column rank of the Khatri–Rao product of two Vandermonde matrices [37,38,39], i.e., $rank[{\mathbf{C}}_1^{M_1\times K}\odot {\mathbf{C}}_2^{M_2\times K}]=min(M_1{M}_2,K)$, we know that, providing $(\overline{N}+1)(\overline{M}+1)\ge K$, $\overline{\mathbf{A}}$ and ${\overline{\mathbf{A}}}^H$ are all almost surely K rank. The spectrum searching algorithms, such as two-dimensional MUSIC algorithm: $\varphi ,\theta \in [-\frac{\pi }{2},\frac{\pi }{2}]$

$$\mathcal{P}(\varphi ,\theta )=\frac{1}{{[\overline{\mathbf{b}}(\varphi )\otimes \overline{\mathbf{a}}(\theta )]}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{0}}{\mathbf{U}}_{\mathbf{0}}^{\mathrm{H}}[\overline{\mathbf{b}}(\varphi )\otimes \overline{\mathbf{a}}(\theta )],}$$

(20)

where ${\mathbf{U}}_{\mathbf{0}}$ denotes the noise subspace of $\overline{\mathbf{Y}}$, and it usually must be at least one-dimension; therefore, $K\le (\overline{N}+1)(\overline{M}+1)-1$ can guarantee the parameter identifiability with probability one. In addition, in the ESPRIT algorithm [18], partitioning the signal subspace ${\mathbf{U}}_{\mathbf{s}}$ for calculating steering matrix has to require $K\le \overline{N}(\overline{M}+1)$. ☐

Remark 1.

In some very special cases, it has $rank(\tilde{\mathbf{A}})<K$ even though $K\le (\overline{N}+1)(\overline{M}+1)-1$; however, the theorem tells us that such cases are measure-zero events, that is to say, $\tilde{\mathbf{A}}$ is almost surely full column rank.

Previous results on the maximum upper bound of the parameter identifiability such as MUSIC-like or rotational invariance algorithms are $MN-1$ when the uniform array configuration is used. Due to the non-uniform feature in the cascade array, it usually satisfies $\overline{M}>M$ and $\overline{N}>N$ so that a much stronger identifiability is acquired. It is worth mentioning that, although the upper-bound of parameter identifiability is derived under the noise-less case, it still works under the limited snapshot number and low signal-to-noise ratio (SNR) case.

4. Joint DOD and DOA Estimation

In this section, we mainly focus on introducing a new computational efficiency joint DOD and DOA estimating algorithm, which is based on a subspace fitting technique.

4.1. Weighted Subspace Fitting

As we know, an asymptotically statistically efficient estimation under a large number of snapshots or high signal-to-noise ratio (SNR) can be obtained by minimizing a weighted subspace fitting [1]:

$$\mathcal{F}(\mathbf{\theta },\mathbf{\varphi })=\mathrm{tr}\left\{{\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }{\mathbf{U}}_{\mathbf{s}}\mathbf{W}{\mathbf{U}}_{\mathbf{s}}^{\mathrm{H}}\right\},$$

(21)

where $\mathbf{\theta }={[{\theta }_1,{\theta }_2,\cdots ,{\theta }_K]}^{\mathrm{T}}$ and $\mathbf{\varphi }={[{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_K]}^{\mathrm{T}}$; the diagonal weighted matrix

$$\mathbf{W}={({\mathbf{\Sigma }}_{\mathbf{s}}-{\widehat{\sigma }}_n^2\mathbf{I})}^2{\mathbf{\Sigma }}_{\mathbf{s}}^{-1}$$

(22)

with the estimated noise variance

$${\widehat{\sigma }}_n^2=\frac{1}{(\overline{M}+1)(\overline{N}+1)-K}\sum _{i=K+1}^{(\overline{M}+1)(\overline{N}+1)}{\lambda }_i$$

(23)

${\mathbf{U}}_{\mathbf{s}}$ is the signal subspace of $\overline{\mathbf{Y}}$, and its corresponding eigenvalues ${\{{\lambda }_i\}}_{i=1}^K$ forms ${\mathbf{\Sigma }}_{\mathbf{s}}$.

In addition, ${\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }$ in Label (21) stands for the orthogonal projector onto the null space of ${(\mathcal{B}\odot \mathcal{A})}^{\mathrm{H}}$, the expression of which is given by,

$${\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }=\mathbf{I}-(\mathcal{B}\odot \mathcal{A}){[{(\mathcal{B}\odot \mathcal{A})}^{\mathrm{H}}(\mathcal{B}\odot \mathcal{A})]}^{-1}{(\mathcal{B}\odot \mathcal{A})}^{\mathrm{H}}.$$

(24)

For minimizing function Label (21), two-dimensional searching in the whole angle-domain is a very direct approach, but it is of computational inefficiency. We herein adopt a MODE-like algorithm that makes use of polynomial rooting. The MODE algorithm was originally proposed in [40] to deal with one-dimensional parameters estimation. For our two-dimensional case, the biggest difficult relies on how to parameterize ${\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }$ because $\mathcal{B}$ and $\mathcal{A}$ are coupling. Therefore, we have to introduce an substitute.

To begin, we first try to parameterize the above projector by a coefficient vectors $\mathbf{b}={[b_0,b_1,\cdots ,b_K]}^{\mathrm{T}}$. These coefficients construct a polynomial with the following form:

$$\sum _{i=0}^K{b}_i{z}^{K-i}=b_0\prod _{i=1}^K(z-e^{j\pi sin{\varphi }_i}),b_0\ne 0.$$

(25)

If we introduce the following set

$$\mathbb{L}=\left\{\{b_i\}\mid \mathcal{C}(z)=\sum _{i=0}^K{b}_i{z}^{K-i}\ne 0\mathrm{for}\mid z\mid \ne 1\right\},$$

(26)

it can be seen that the mapping from $\{{\varphi }_i\}\in \mathbb{R}$ to $\{b_i\}\in \mathbb{L}$ is one-to-one providing we eliminate the non-uniqueness implied by the introduction of $b_0\ne 1$.

Let ${\mathbf{G}}_{\mathbf{b}}\in {\mathbb{C}}^{(\overline{N}+1)\times (\overline{N}+1-K)}$ be the following Toeplitz matrix:

$${\mathbf{G}}_{\mathbf{b}}^{\mathrm{H}}=\left[\begin{array}{cccccc}{b}_K& \cdots & b_1& b_0& \cdots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & 0\\ 0& \cdots & b_K& \cdots & b_1& b_0\end{array}\right].$$

(27)

It is observed that $\mathrm{rank}\{{\mathbf{G}}_{\mathbf{b}}\}=\overline{N}+1-K$, and

$${\mathbf{G}}_{\mathbf{b}}^{\mathrm{H}}\mathcal{B}=\mathbf{0}.$$

(28)

Based on the above relation, we can conclude the following theorem that can be utilized to estimate all the DOA information.

Theorem 3.

Let the columns of ${\mathbf{G}}_{\mathbf{b}}$ span the null space of ${\mathbf{B}}^{\mathrm{H}}$, and if defining $\mathbf{G}={\mathbf{G}}_{\mathbf{b}}\otimes {\mathbf{I}}_{\overline{M}+1}$, then $\mathit{span}\{\mathbf{G}\}\subset \mathit{span}\{{\mathbf{U}}_0\}$.

Proof. 

According to Label (28), the columns of $\mathbf{G}$ span a column space that can guarantee the following result:

$${\mathbf{G}}^{\mathrm{H}}(\mathcal{B}\odot \mathcal{A})=({\mathbf{G}}_{\mathbf{b}}^{\mathrm{H}}\mathcal{B})\odot \mathcal{A}=\mathbf{0}.$$

On the other hand, if considering the rank of matrix $\mathbf{G}$, it is shown that

$$\mathrm{rank}\{\mathbf{G}\}=\mathrm{rank}\{{\mathbf{G}}_{\mathbf{b}}\}\times \mathrm{rank}\{{\mathbf{I}}_M\}=(\overline{N}+1-K)(\overline{M}+1).$$

(29)

As we know, the rank of noise subspace ${\mathbf{U}}_0$ is $\mathrm{rank}\{{\mathbf{U}}_0\}=(\overline{M}+1)(\overline{N}+1)-K$. Obviously, if and only if $\overline{M}\ge 1$, then

$$\mathrm{rank}\{{\mathbf{U}}_0\}>\mathrm{rank}\{\mathbf{G}\}.$$

(30)

Such result directly demonstrates that the column space spanned by the columns of $\mathbf{G}$ is included in the null space spanned by ${\mathbf{U}}_0$. This completes the proof. ☐

Therefore, we can construct a projection matrix

$${\mathbf{P}}_{\mathbf{G}}=\mathbf{G}{({\mathbf{G}}^{\mathrm{H}}\mathbf{G})}^{-1}{\mathbf{G}}^{\mathrm{H}}={\mathbf{P}}_{{\mathbf{G}}_{\mathbf{b}}}\otimes {\mathbf{I}}_{\overline{M}+1},$$

(31)

where ${\mathbf{P}}_{{\mathbf{G}}_{\mathbf{b}}}={\mathbf{G}}_{\mathbf{b}}{({\mathbf{G}}_{\mathbf{b}}^{\mathrm{H}}{\mathbf{G}}_{\mathbf{b}})}^{-1}{\mathbf{G}}_{\mathbf{b}}^{\mathrm{H}}$ ; then, we substitute ${\mathbf{P}}_{\mathbf{G}}$ for ${\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }$ in Label (21), and, correspondingly, a new objective function that needs to be minimized is reformulated as

$$\mathcal{F}(\mathbf{b})=\mathrm{tr}\left\{({\mathbf{P}}_{{\mathbf{G}}_{\mathbf{b}}}\otimes {\mathbf{I}}_M){\mathbf{U}}_{\mathbf{s}}\mathbf{W}{\mathbf{U}}_{\mathbf{s}}^{\mathrm{H}}\right\}.$$

(32)

If comparing function $\mathcal{F}(\mathbf{\theta },\mathbf{\varphi })$ with function $\mathcal{F}(\mathbf{b})$, we can see that: on the one hand, the original two-dimensional optimization problem is transformed into an one-dimensional problem, that is to say, it greatly decreases the complexity; on the other hand, the projector matrix ${\mathbf{P}}_{\mathcal{B}\odot \mathcal{A}}^{\perp }$ is parameterized by a series of coefficients. However, such operation combined with the subsequent angle estimation algorithm inevitably incurs performance loss if comparing with the two-dimensional spectrum searching method for optimizing $\mathcal{F}(\mathbf{\theta },\mathbf{\varphi })$. The related analysis and comparison for such performance loss will be reserved for future work.

4.2. Angle Estimation

In order to implement the minimization of Label (32) more conveniently, we define

$$\overline{\mathbf{U}}={\mathbf{U}}_{\mathbf{s}}\mathbf{W}{\mathbf{U}}_{\mathbf{s}}^{\mathrm{H}}=[\begin{array}{ccc}{\overline{\mathbf{U}}}_{11}& \cdots & {\overline{\mathbf{U}}}_{1N}\\ ⋮& \ddots & ⋮\\ {\overline{\mathbf{U}}}_{N1}& \cdots & {\overline{\mathbf{U}}}_{NN}\end{array}],$$

(33)

where ${\overline{\mathbf{U}}}_{uv}$ is a $(\overline{M}+1)\times (\overline{M}+1)$ matrix, $u,v=1,\cdots ,\overline{N}+1$. Furthermore, we can get

$$\widehat{\mathbf{b}}=argmin\mathrm{tr}\left\{{\mathbf{P}}_{{\mathbf{G}}_{\mathbf{b}}}\left[\begin{array}{ccc}\mathrm{tr}({\overline{\mathbf{U}}}_{11})& \cdots & \mathrm{tr}({\overline{\mathbf{U}}}_{1N})\\ ⋮& \ddots & ⋮\\ \mathrm{tr}({\overline{\mathbf{U}}}_{N1})& \cdots & \mathrm{tr}({\overline{\mathbf{U}}}_{NN})\end{array}\right]\right\}.$$

(34)

In addition, we also exert such constraints on the unknown parameters ${\{b_i\}}_{i=1}^{\overline{N}+1}$, i.e., $b_i=b_{\overline{N}+1-i}^{*}$. The detailed discussion with respect to the above constraint and procedures for minimizing the above quadratic function can be found in [33,41].

Once we get $\widehat{\mathbf{b}}$, the angular phase of the roots of the estimated polynomial in Label (25) will give the DOA information of all targets, i.e., ${\{{\widehat{\varphi }}_i\}}_{i=1}^K$.

For DOD information, the following method is adopted. We convert the two-dimensional MUSIC algorithm, see Label (20), into two optimization problems [42,43]

$$\begin{array}{c}\underset{\varphi }{max}{\mathbf{e}}^{\mathrm{T}}{\mathbf{E}}^{-1}(\varphi )\mathbf{e},\hfill \end{array}$$

(35)

$$\begin{array}{c}\underset{\theta }{min}\overline{\mathbf{a}}{(\theta )}^{\mathrm{H}}\mathbf{E}(\varphi )\overline{\mathbf{a}}(\theta ),s.t.{\mathbf{e}}^{\mathrm{T}}\overline{\mathbf{a}}(\theta )=1,\hfill \end{array}$$

(36)

where $\mathbf{E}(\varphi )={[\overline{\mathbf{b}}(\varphi )\otimes {\mathbf{I}}_{\overline{M}+1}]}^{\mathrm{H}}{\mathbf{U}}_{\mathbf{0}}{{\mathbf{U}}_{\mathbf{0}}}^{\mathrm{H}}[\overline{\mathbf{b}}(\varphi )\otimes {\mathbf{I}}_{\overline{M}+1}]$ and $\mathbf{e}={[1,0,\cdots ,0]}^{\mathrm{T}}$.

Due to the fact that we have obtained DOA information ${\{{\varphi }_k\}}_{k=1}^K$ through coefficient vector $\widehat{\mathbf{b}}$, the auto-paired DOD information ${\theta }_k$ can be drawn from the estimated transmit steering vectors, i.e.,

$$\widehat{\overline{\mathbf{a}}({\theta }_k)}=\frac{\mathbf{E}{({\widehat{\varphi }}_k)}^{-1}\mathbf{e}}{{\mathbf{e}}^{\mathrm{T}}\mathbf{E}{({\widehat{\varphi }}_k)}^{-1}\mathbf{e}},$$

(37)

which can be achieved easily by utilizing Lagrange multiplier method to Label (36), and consequently

$${\widehat{\theta }}_k=arcsin\left\{\frac{1}{\pi \overline{M}}\sum _{\overline{m}=1}^{\overline{M}}\angle \frac{{\overline{\mathbf{a}}}_k[\overline{m}+1]}{{\overline{\mathbf{a}}}_k[\overline{m}]}\right\},$$

(38)

where ${\overline{\mathbf{a}}}_k[\overline{m}]$ denotes the $\overline{m}$-th element of $\widehat{\overline{\mathbf{a}}({\theta }_k)}$.

4.3. Computational Complexity

We now analyze the computational complexity of the above algorithm. By neglecting some trivial operations,

Table 1. Computational complexity analysis.

Operation	Dimension Size	Required Flops
Eigenvalue Decompositon of $\overline{\mathbf{Y}}$: ${\mathbf{U}}_{\mathbf{s}}$ and ${\mathbf{\Sigma }}_{\mathbf{s}}$	$MN\times Q$	$Q{M}^2{N}^2+{\tilde{M}}^3{\tilde{N}}^3$
Matrix $\overline{\mathbf{U}}$ Partition in Label (34)	$\tilde{M}\tilde{N}\times \tilde{M}\tilde{N}$	$\frac{1}{2}\tilde{M}{\tilde{N}}^2$
Minimization of $\mathcal{F}(\mathbf{b})$	$\tilde{N}\times \tilde{N}$	$\frac{1}{2}\zeta K(\tilde{N}-K)(20K^2+33K+17)$
$\overline{\mathbf{a}}(\theta )$ Estimation	$\tilde{M}\times \tilde{M}$	${\tilde{M}}^3$
Total flops		$\mathcal{O}\{{\tilde{M}}^3{\tilde{N}}^3\}$

lists the number of flops required in major steps. Here, the SVD is assumed to be performed by the Golub–Reinsch SVD algorithm [44]. For convenience, let $\tilde{M}=\overline{M}+1$ and $\tilde{N}=\overline{N}+1$.

In Table 1, parameter $\zeta$ denotes the iterative number for optimization of MODE algorithm, which is usually selected as 4 [33]. Compared with the proposed algorithm, tensor-based algorithms [34,35] have to pay out large computational cost for calculating high-order singular value decomposition; and they also have to afford two-times of computational burden because two quadratic functions need separately minimization in the procedures of performing MODE algorithm. Consequently, the angle pairing is inevitable. In addition, the proposed algorithm requires more flops than two-dimensional ESPRIT algorithm (it requires $Q{M}^2{N}^2+{\tilde{M}}^3{\tilde{N}}^3+3\overline{N}\tilde{M}{K}^2+26K^3+\tilde{M}\tilde{N}{K}^2$ flops) and much less than the two-dimensional MUSIC algorithm.

5. Numerical Simulation

In this section, we will fully demonstrate the effectiveness of the designed array configuration and the proposed angle estimation algorithm by Monte Carlo simulations. In the following examples, the Doppler frequency $f_{dk}$ is generated by $f_{dk}=(2\pi {\upsilon }_k{T}_p)/\lambda$, where ${\upsilon }_k$ is the target velocity, $T_p=5\times {10}^{-6}$ is the pulse duration in seconds, and $\lambda =3\times {10}^8/f_c$ with $f_c=3$ GHz. The average root mean square error (RMSE) vs. signal-noise-ratio (SNR) is used as our performance assessment, where the SNR is defined according to Lable (7):

$$\mathrm{SNR}=10{log}_{10}\frac{{\parallel (\mathbf{B}\odot \mathbf{A})\mathbf{H}\parallel }^2}{MN{\sigma }_n^2},$$

and the average root-mean-square error (RMSE) is defined as

$$(1/K)\sum _{k=1}^K\sqrt{E\left[{({\widehat{\theta }}_k-{\theta }_k)}^2+{({\widehat{\varphi }}_k-{\varphi }_k)}^2\right].}$$

5.1. DOF Comparison

We first compare different array configurations from the perspective of maximum consecutive DOFs. As the results shown in Figure 2, we can conclude that the designed cascade array can generate more DOFs than other exist arrays such as ULA, nested array [13], super nested array [16], CACIS and nested CADiS [15]; and reaches a closed performance to the DOFs upper bound of the MR array [20].

5.2. Identifiability Performance for Cascade Array

In order to sufficiently examine the identifiability performance of the designed cascade array under the underdetermined and normal scenarios, we present a one-dimensional DOA case and MIMO DOD-DOA case through the following examples.

1D underdetermined scenario: We arrange such a underdetermined scenario that nine far-field narrow-band uncorrelated signals impinge an array with only $N=5$ sensors from direction-of-arrival (DOA) of $\{-{60}^{\circ },-{45}^{\circ },-{30}^{\circ },-{15}^{\circ },0^{\circ },{15}^{\circ },{30}^{\circ },{45}^{\circ },{60}^{\circ }\}$. The sources are modeled as random Gaussian processes and all with equal power. The signal-to-noise ratio is 10 dB and the number of snapshots is 500. From Figure 3, we can see that all sources can be resolved sufficiently well by cascade array: $\{0,1,4,7,9\}$. However, other array geometries cannot provide enough DOFs to fulfill the requirement of such spectrum estimation, e.g., the nested array: $\{0,1,2,5,8\}$; the CACIS: $\{0,1,2\}$ and $\{0,3,6\}$; the ULA: $\{0,1,2,3,4\}$; the super nested array cannot be constructed due to its requirement for $n\ge 7$. In addition, the nested CADiS: $\{0,1,2,6,9\}$ can perform the same resolving ability as our cascade array under that underdetermined case, but its identifiability becomes smaller than the cascade array when sensor number n increases (see Figure 2).

2D underdetermined scenario: We consider an underdetermined MIMO array systems with $M=N=5$ sensors. For traditional MIMO array with ULA configuration, according to the almost surely full column rank of the Khatri–Rao product of two Vandermonde matrixes, the upper bound of the parameter identifiability is $MN-1=24$; however, the virtual MIMO array based on the cascade array is $(\overline{N}+1)(\overline{M}+1)-1=99$. We set 25 targets with reflected signals coming from two-dimensional angle-domain uniformly, the RCS coefficients of which are set as ${\{{\beta }_k\}}_{k=1}^{25}=1$. The SNR is 10 dB and the number of snapshots is 500. The MUSIC spectrum $\mathcal{P}(\varphi ,\theta )$ we used is given by Label (20). From Figure 4, we can see that, under such two-dimensional underdetermined scenario, the designed cascade array can distinguish all 30 targets successfully.

2D normal scenario: Furthermore, Figure 5 plots the two-dimension MUSIC spectrum under a normal scenario: $K=5$ targets with $\theta =\{{40}^{\circ },{35}^{\circ },{30}^{\circ },-{40}^{\circ },{65}^{\circ }\}$ and $\varphi =\{{20}^{\circ },{25}^{\circ },{30}^{\circ },{50}^{\circ },-{45}^{\circ }\}$, i.e., for three closely spaced targets and two targets widely spaced from the others. The other parameters are $\{{\beta }_k\}=\{1,1,1,1,1\}$, $Q=50$ and SNR $=8$ dB. In this simulation, we also consider the random array configuration. According to [28,29], the total aperture for both transmitting and receiving array is set as 10 half-wavelengths. To ensure the aperture length, we fixed the location of two elements at the extremities and placed the rest of the elements uniformly at random in between. In this single realization of random array, the transmitting array is $\{0,2.7581,3.1378,4.7333,10\}$ and the receiving array is $\{0,0.6562,2.8234,7.8127,10\}$. From the simulation results, one can see that: (1) for $M=N=5$ ULA configuration, the two-dimensional MUSIC spectrum cannot give a accurate localization for that three closely spaced targets because that three peaks merged, while the other targets are well located; (2) for $M=N=5$ random array configuration, the distinguishing performance goes better; (3) for the cascade array configuration, it significantly improves the spatial resolution not only in the closely spaced targets becoming distinguishable, but also in the sidelobe suppression.

5.3. RMSE Performance for Joint DOD and DOA Estimation

In Figure 6, we compare the performance of some existing algorithms, including the tensor MODE algorithm with $M=N=8$ uniform linear array configuration [34], 2D-ESPRIT algorithm with $M=N=7$ nested array configuration [13] (Note that the super nested array [16] has the same DOF as the nested array when the number of antenna is larger than 7, so we herein just consider the nested array configuration), 2D-ESPRIT algorithm [19] with $M=N=8$ coprime array configuration [14] (It consists of two ULAs with five sensors and $2\times 2$ sensors, respectively, and their first sensors are of coincidence) and the proposed algorithm with $M=N=7$ cascade array configuration. In addition, we also compare our earlier proposed method [18], which is based on the MR array. We consider two cases, i.e., $K=6$ widely separated targets: $(\theta ,\varphi )=\{(-{80}^{\circ },{70}^{\circ }),(-{60}^{\circ },{10}^{\circ }),(-{40}^{\circ },{50}^{\circ }),(-{20}^{\circ },-{30}^{\circ }),$$(0^{\circ },-{10}^{\circ }),({20}^{\circ },{30}^{\circ })\}$ with the RCS coefficients $\{{\beta }_k\}=\{1,0.9,0.8,0.7,0.6,0.5\}$ and closely separated targets: $(\theta ,\varphi )=\{(-{80}^{\circ },{70}^{\circ }),(-{75}^{\circ },{65}^{\circ }),(-{40}^{\circ },{50}^{\circ }),(-{35}^{\circ },{45}^{\circ }),$ $(0^{\circ },-{10}^{\circ }),(5^{\circ },-{15}^{\circ })\}$ with the RCS coefficients $\{{\beta }_k\}=\{1,0.85,0.7,0.55,0.4,0.25\}$. We fix the number of pulses $Q=200$ for both cases. For the average RMSE performance evaluation, the global bound—Ziv–Zakai bound (ZZB) [45,46] and the local bound, Cramér–Rao lower bound (CRLB) [47]—are all provided as benchmarks.

From the first simulation results in Figure 6, we can see that the designed cascade array with the proposed algorithm performs the same as the method based on the MR array [18] for the widely separated targets, and performs a little worse for the closely separated targets. However, it outperforms other types of array configurations. The major manifestation lies in the lower average RMSE of DOD-DOA estimation can be attained. Furthermore, in the second example, the performance gap between the nested array and cascade array when $M=N=8$ is larger than the one when $M=N=7$. Such advantage benefits from the strengthened DOF generated by the designed cascade array. According to Section 2, for transmitting array or receiving array, the DOF difference of the co-arrays between cascade array and nested array is 4 when $M=N=7$ (The cascade array provides 35 DOFs and the nested array provides 31 DOFs); then, it increases to 6 when $M=N=8$ (The cascade array provides 45 DOFs and the nested array provides 39 DOFs). On the other hand, such advantage cannot persist in the case with many more antennas. For instance, seeing the case of $M=N=9$ in the second example, the performance improvement will become smaller and smaller with the SNR increasing. The major reason behind this phenomenon is that there exists a model error during the procedures of achieving co-array, i.e., the matrix $\mathbf{\Lambda }$ of Label (8) is approximately diagonal in estimating the actual correlation matrix $\mathbf{R}$ so that it plays a dominant adverse role in the high SNR region. In the second example, undoubtedly, the method based on the MR array [18] is optimal because it has an upper-bound of DOFs. For the RMSE lower bound, in calculating ZZB, the prior distribution of the DOD and DOA ($sin\theta$ and $sin\varphi$) are uniform on the unit disc with covariance $\frac{1}{4}\mathbf{I}$ [46]. For the widely separated targets scenario in Figure 6a, the ZZB merges with the CRLB, but it becomes a loose bound under the closely separated targets scenario. The reason is that the ZZB depends on the prior distribution of parameters rather than the parameters so that both scenarios share the same ZZB; however, the CRLB varies with different scenarios. In addition, we do not plot the performance at much lower SNR. The reason comes from Figure 7, that is to say, the detection of targets number, which is a prerequisite for angle estimation, cannot guarantee being absolutely right when the SNR is very low.

It is worth mentioning that the above result does not mean the ZZB is inferior to the CRLB. As we know [45,46,48], the performance of angle estimation is characterized by the presence of distinct SNR regions. Besides the asymptotic region we consider in this simulation, there still exist the a priori performance region (characterized by small snapshots and/or lower SNR) and an additional ambiguity region (or transition region). The ZZB can provide a much tighter bound in all regions and even can accurately identify the SNR thresholds, whereas the CRLB is not qualified.

5.4. Other Related Performance

Targets number Estimation: In the aforementioned simulation examples, the targets number is assumed to be known as a prior condition. In Figure 7, we try to consider the targets number detection performance under several different sensor array configurations. In this simulation, the second order statistic of the eigenvalues (SORTE) algorithm [49] is adopted, which is an eigenvalue-based strategy. The all parameters setting is the same as the previous example, that is to say, the targets number $K=6$ that needs to be detected. The detection probability (DP), ${\mathcal{N}}_K/{\mathcal{N}}_{tol}$, versus SNR is depicted, in which the ${\mathcal{N}}_{tol}$ is the total trials and ${\mathcal{N}}_K$ is the number of times that K is successfully detected. From the simulation results, we can see that the designed cascade array configuration outperforms the nested array, coprime array and ULA with the same number of antennas.

Performance with different $M,N$ and different snapshot number: In the following, we first fix the total number of transmitting-receiving antennas, i.e., $M+N=14$, to examine which combination of $\{M,N\}$ can provide the optimal average RMSE performance. The angle setting is $(\theta ,\varphi )=\{({45}^{\circ },{50}^{\circ }),({35}^{\circ },{30}^{\circ }),({15}^{\circ },{10}^{\circ }),(-{15}^{\circ },-{10}^{\circ }),$ $(-{25}^{\circ },-{30}^{\circ })\}$. From Figure 8, we can see that the simulation results under $M=N=7$ have the least average RMSE. In addition, we also consider the angle estimation performance versus different number of pulses, as is shown in Figure 9, which demonstrates that the increase of pulse number can also improve the target’s localization accuracy. It is worth mentioning that a large pulse number will make matrix $\mathbf{\Lambda }$ much more diagonally dominant; consequently, it can diminish to some extent the model error.

6. Conclusions

In this paper, we mainly focus on two problems in the multiple input multiple output array systems: one is to design a new linear array structure, which can provide more degrees of freedom than some existing arrays from the perspective of different co-arrays; the other is to propose a new targets’ localization algorithm with respect to DOA and DOD, which contains only one-time polynomial rooting and two-dimension angle estimation auto-pairing. We not only analyze systematically the optimal configuration of the designed cascade array under a given number of antennas and the maximal parameter identifiability when applying it into the MIMO array systems, but also provide some simulation proofs from the angle estimation accuracy and the target’s number detection probability to verify the effectiveness and advantages of the designed array structure and the proposed targets’ localization algorithm.

In the future work, we will mainly focus on three aspects. First, more efficient angle estimation algorithms should be devised because there exists an obvious performance gap with the Cramér–Rao lower bound; second, the cascade array or other array geometries should be optimized or designed to generate more DOFs because there exists a DOF gap with the MR array; finally, the applications to other different scenarios with the proposed algorithm and array geometry are also worth considering.