Transmit Array Interpolation for DOA Estimation via Tensor Decomposition in 2D MIMO Radar

TL;DR

This paper introduces a 2D joint transmit array interpolation and tensor-based DOA estimation method for MIMO radar, improving accuracy and SNR by mapping the array to a desired configuration and compensating for interpolation bias.

Contribution

It proposes a novel tensor-based DOA estimation approach with a closed-form bias correction and CRB analysis for 2D MIMO radar systems.

Findings

Enhanced DOA estimation accuracy demonstrated in simulations

Bias correction improves estimation reliability

Performance approaches the derived CRB

Abstract

In this paper, we propose a two-dimensional (2D) joint transmit array interpolation and beamspace design for planar array mono-static multiple-input-multiple-output (MIMO) radar for direction-of-arrival (DOA) estimation via tensor modeling. Our underlying idea is to map the transmit array to a desired array and suppress the transmit power outside the spatial sector of interest. In doing so, the signal-tonoise ratio is improved at the receive array. Then, we fold the received data along each dimension into a tensorial structure and apply tensor-based methods to obtain DOA estimates. In addition, we derive a close-form expression for DOA estimation bias caused by interpolation errors and argue for using a specially designed look-up table to compensate the bias. The corresponding Cramer-Rao Bound (CRB) is also derived. Simulation results are provided to show the performance of the proposed method and compare its performance to CRB.

Tables

Table 1. TABLE I: Algorithm for Building Offline DOA Look-up Table

(i) Consider the transmit array interpolation model (26) and assume a single source located at $(\theta ,\varphi )\in (\Theta ,\Phi )$; map ($\theta ,\varphi$) to ($\mu ,\nu$) using $\mu =\sin \theta \cos \varphi$ and $\nu =\sin \theta \sin \varphi$.

(ii) Assuming that the observation model is noise free, build the corresponding noise-free tensor model of (29).

(iii) Apply the proposed DOA estimation methods to the noise-free tensor model to obtain ($\hat{\mu },\hat{\nu }$).

(iv) Repeat (i)–(iii) over a fine 2D grid within the 2D spatial sector of interest.

(v) Form two look-up tables mapping $\mu$ to $\hat{\mu }$ and $\nu$ to $\hat{\nu }$.

Equations

$$c\triangleq <\mathcal{A},\mathcal{D}>\stackrel{△}{=}\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\cdots \sum _{i_N=1}^{I_N}{d}_{i_1,i_2,\dots ,i_N}^{\ast }\cdot a_{i_1,i_2,\dots ,i_N}.$$

$$\Vert \mathcal{A}{\Vert }_{\text{H}}\triangleq \sqrt{<\mathcal{A},\mathcal{A}>}.$$

$$\mathcal{C}\triangleq \mathcal{A}{\times }_n\mathbf{D}$$

$$c_{i_1,i_2,\dots ,i_{n-1},j_n,i_{n+1},\dots ,i_N}=\sum _{i_n=1}^{I_n}{a}_{i_1,i_2,\dots ,i_N}\cdot d_{j,i_n}$$

$$\mathcal{C}\triangleq \mathcal{A}\circ \mathbf{D}\in {\mathbb{C}}^{I_1\times I_2\times \cdots \times I_N\times J_1\times J_2\times \cdots \times J_M}$$

$$c_{i_1,i_2,\dots ,i_N,j_1,j_2,\dots ,j_M}=a_{i_1,i_2,\dots ,i_N}\cdot d_{j_1,j_2,\dots ,j_M}.$$

$$\mathbf{a}(\theta ,\varphi )\triangleq [e^{-j2\pi {\mathbf{u}}^T(\theta ,\varphi ){\mathbf{p}}_1},\dots ,e^{-j2\pi {\mathbf{u}}^T(\theta ,\varphi ){\mathbf{p}}_M}{]}^T$$

$$\mathbf{b}(\theta ,\varphi )\triangleq [e^{-j2\pi {\mathbf{u}}^T(\theta ,\varphi ){\mathbf{p}}_1},\dots ,e^{-j2\pi {\mathbf{u}}^T(\theta ,\varphi ){\mathbf{p}}_N}{]}^T.$$

$$\zeta (t,\theta ,\varphi )={\mathbf{a}}^T(\theta ,\varphi )\mathbf{s}(t)=\sum _{m=1}^M{a}_m(\theta ,\varphi )s_m(t)$$

$$\mathbf{x}(t,q)=\sum _{k=1}^K{\beta }_k(q)({\mathbf{a}}^T({\theta }_k,{\varphi }_k)\mathbf{s}(t))\mathbf{b}({\theta }_k,{\varphi }_k)+\mathbf{n}(t,q)$$

$$\bm{y}_{m}(q)=\sum_{k=1}^{K}\big{(}a_{m}(\theta_{k},\phi_{k})\bm{b}(\theta_{k},\phi_{k})\big{)}\beta_{k}(q)+\bm{z}(q)$$

$$\mathbf{y}(q)=(\mathbf{A}(\theta ,\varphi )\odot \mathbf{B}(\theta ,\varphi ))\mathbf{\beta }(q)+\mathbf{z}(q)$$

$${\mathbf{E}}^H\mathbf{a}(\theta ,\varphi )=\tilde{\mathbf{a}}(\theta ,\varphi )$$

$$\tilde{\mathbf{a}}(\theta ,\varphi )=\mathbf{c}(\mu )\odot \mathbf{d}(\nu )$$

$$\mathbf{c}(\mu )\triangleq [1,e^{-j2\pi x_1\mu },\dots ,e^{-j2\pi x_{(M_1-1)}\mu }{]}^T$$

$$\mathbf{d}(\nu )\triangleq [1,e^{-j2\pi y_1\nu },\dots ,e^{-j2\pi y_{(M_2-1)}\nu }{]}^T.$$

$$\mathbf{A}(\Theta ,\Phi )\triangleq [\mathbf{a}({\theta }_1,{\varphi }_1),\dots ,\mathbf{a}({\theta }_{G_{\Theta }},{\varphi }_1),\dots ,\mathbf{a}({\theta }_{G_{\Theta }},{\varphi }_{G_{\varphi }})]$$

$$\tilde{\mathbf{A}}(\Theta ,\Phi )\triangleq [\tilde{\mathbf{a}}({\theta }_1,{\varphi }_1),\dots ,\tilde{\mathbf{a}}({\theta }_{G_{\Theta }},{\varphi }_1),\dots ,\tilde{\mathbf{a}}({\theta }_{G_{\Theta }},{\varphi }_{G_{\varphi }})].$$

$$\mathbf{E}={\left(\mathbf{A}(\Theta ,\Phi ){\mathbf{A}}^H(\Theta ,\Phi )\right)}^{-1}\mathbf{A}(\Theta ,\Phi ){\tilde{\mathbf{A}}}^H(\Theta ,\Phi ).$$

$$\underset{\mathbf{E}}{\min }\underset{{\theta }_g,{\varphi }_{g^{\prime }}}{\max }{\Vert {\mathbf{E}}^H\mathbf{a}({\theta }_g,{\varphi }_{g^{\prime }})-\tilde{\mathbf{a}}({\theta }_g,{\varphi }_{g^{\prime }})\Vert }_p$$

$${\theta }_g\in \Theta ,g=1,\dots ,G_{\theta },{\varphi }_{g^{\prime }}\in \Phi ,g^{\prime }=1,\dots ,G_{\varphi }$$

$$\mathrm{subject}\mathrm{to}{\Vert {\mathbf{E}}^H\mathbf{a}({\theta }_h,{\varphi }_{h^{\prime }})\Vert }_p\le \gamma ,$$

$${\theta }_h\in \bar{\Theta },h=1,\dots ,H_{\theta },{\varphi }_{h^{\prime }}\in \bar{\Phi },h^{\prime }=1,\dots ,H_{\varphi }$$

$$\underset{\mathbf{E}}{\min }\underset{{\theta }_h,{\varphi }_{h^{\prime }}}{\max }{\Vert {\mathbf{E}}^H\mathbf{a}({\theta }_h,{\varphi }_{h^{\prime }})\Vert }_p$$

$${\theta }_h\in \bar{\Theta },n=1,\dots ,H_{\theta },{\varphi }_{h^{\prime }}\in \bar{\Phi },h^{\prime }=1,\dots ,H_{\varphi }$$

$$\mathrm{subject}\mathrm{to}{\Vert {\mathbf{E}}^H\mathbf{a}({\theta }_g,{\varphi }_{g^{\prime }})-\tilde{\mathbf{a}}({\theta }_g,{\varphi }_{g^{\prime }})\Vert }_p\le \Delta$$

$${\theta }_g\in \Theta ,g=1,\dots ,G_{\theta },{\varphi }_{g^{\prime }}\in \Phi ,g^{\prime }=1,\dots ,G_{\varphi }$$

$$\tilde{\zeta }(t,\theta ,\varphi )$$

$$={\mathbf{a}}^T(\theta ,\varphi ){\mathbf{E}}^{\ast }\tilde{\mathbf{s}}(t)=\sum _{\tilde{m}=1}^{\tilde{M}}{\mathbf{a}}^T(\theta ,\varphi ){\mathbf{e}}_{\tilde{m}}^{\ast }{s}_{\tilde{m}}(t)$$

$$\tilde{\mathbf{x}}(t,q)=\sum _{k=1}^K{\beta }_k(q)({\tilde{\mathbf{a}}}^T({\theta }_k,{\varphi }_k)\tilde{\mathbf{s}}(t))\mathbf{b}({\theta }_k,{\varphi }_k)+\mathbf{n}(t,q).$$

$${\tilde{\mathbf{y}}}_{\tilde{m}}(q)$$

$$\displaystyle=\sum_{k=1}^{K}\big{(}\tilde{\bm{a}}_{\tilde{m}}(\theta_{k},\phi_{k})\odot\bm{b}(\theta_{k},\phi_{k})\big{)}\beta_{k}(q)+\bm{z}(q)$$

$$=\sum _{k=1}^K{\mathbf{f}}_{\tilde{m}}({\theta }_k,{\varphi }_k){\beta }_k(q)+\mathbf{z}(q)$$

$$={\mathbf{F}}_{\tilde{m}}(\theta ,\varphi )\mathbf{\beta }(q)+\mathbf{z}(q)$$

$${\mathbf{F}}_{\tilde{m}}(\theta ,\varphi )\triangleq [{\mathbf{f}}_{\tilde{m}}({\theta }_1,{\varphi }_1),\dots ,{\mathbf{f}}_{\tilde{m}}({\theta }_K,{\varphi }_K)]$$

$$\tilde{\mathbf{y}}(q)=\mathbf{F}(\theta ,\varphi )\mathbf{\beta }(q)+\mathbf{z}(q)$$

$$\mathbf{F}(\theta ,\varphi )\triangleq [{\mathbf{F}}_1^T(\theta ,\varphi ),\dots ,{\mathbf{F}}_{\tilde{M}}^T(\theta ,\varphi ){]}^T\in {\mathbb{C}}^{\tilde{M}N\times K}$$

Full text

I Introduction

The multiple-input-multiple-output (MIMO) radar, including the MIMO radar with transmit beamspace design, has received considerable attention during the last decade [1]-[10]. Unlike phase-array radar, MIMO radar can simultaneously transmit several orthogonal waveforms through co-located antenna elements or widely separated antenna elements [4], [5]. Such waveform diversity brings many advantages to the MIMO radar over its phased-array radar counterpart among which are the possibility of flexible waveform design, improved parameter identifiability, higher spatial resolution, to name a few [1]-[8]. Moreover, MIMO radar has an inherent higher-dimensional structure that can be exploited to improve the estimation performance [11].

For a co-located MIMO radar, the transmit and receive antenna elements form particular array configurations, which provide a virtual array with increased number of antenna elements. This enables, for example, an improvement for direction-of-arrival (DOA) estimation performance. Many methods have been proposed to deal with the DOA estimation problem in MIMO radar, and most of them are transplanted from phased-array radar such as MUltiple SIgnal Classification (MUSIC), estimation of signal parameters via rotational invariance technique (ESPRIT) [9], maximum likelihood (ML) [12]. Another class of DOA estimation methods exploits the higher-dimensional structure inherent in the MIMO radar [11]. There are two main decomposition methods, parallel factor analysis (PARAFAC) [13], [14] and higher-order singular value decomposition (HOSVD) [15], [16]. PARAFAC methods require a more relaxed uniqueness condition and have shown better estimation performance than covariance matrix singular value decomposition (SVD)-based methods. HOSVD utilizes the tensor structure, i.e., it iteratively performs SVD to every dimension of the received data. In the related R-dimensional harmonic retrieval (HR) parameter estimation problem, it is shown that the tensor-based methods have better DOA estimation performance and resolution ability than the matrix-based ones [17], [18].

Both arrival and departure directions are often of interest, moreover, the transmit array or both the transmit and receive arrays are typically two-dimensional (2D). In the co-located MIMO radar, since the transmit waveforms are independently received by the receive array when transmitted from 2D array, the DOA estimation boils down to 2D (azimuth and elevation) parameter estimation. Several 2D estimation methods have been introduced in the context of MIMO radar [19], [20]. Typically, a regular array geometries such as uniform rectangular array (URA) [21] or L-shaped array [22], [23] are assumed or other properties such as the existence of orthogonally polarized components of linear arrays [24] are used to reduce the search complexity while designing the DOA estimation techniques.

The well known drawback of MIMO radar, which employs full waveform diversity, i.e., every transmit antenna element radiates orthogonal waveform, in the regime of parameter estimation is the reduction of the signal-to-noise ratio (SNR) gain comparing with phased-array radar at a fixed total transmit energy [7], [8]. In order to improve the SNR gain in the known spatial sector of interest, coherent transmit beamspace (TB) design techniques have been developed [8], [25]. These designs are called transmit beamspace MIMO radar, and they utilize the fact that the target spatial region is usually a priori known before solving the parameter estimation problem, thus the transmit beam could focus the transmit energy within the desired spatial region while suppressing the sidelobes outside of it. The TB MIMO radar takes the advantages of the phased-array and MIMO radars and enjoys both the waveform diversity and high SNR gain within the sector of interest. As a result, comparing to MIMO radar with full diversity, TB MIMO radar has better DOA estimation performance and lower Cramér-Rao bound (CRB), moreover, the number of waveforms for achieving the best perfomance can be found optimally [8].

Another merit of TB MIMO radar is the possibility to reduce the DOA estimation computational burden [26]-[30]. Specifically, for an arbitrary planar array MIMO radar, there are more than one parameter to be estimated such as elevation/azimuth, polarization, Doppler, etc. Thus, grid searching methods, e.g., spectral MUSIC, become very time consuming or infeasible. A nonuniform transmit array structure, which does not lead to Vandermonde structure for the steering vectors along any of the dimensions, prohibits the application of computationally efficient methods based on polynomial rooting instead of a grid searching. The TB design methods, however, allow to interpolate the transmit array into a virtual transmit array which has a desirable regular structure [26], [30]. The interpolation can be designed so that the signal subspace of the received snapshots will enjoy some desirable properties such as the rotational invariance property. In turn, it enables us to use polynomial rooting based parameter estimation techniques, which greatly reduces the computational burden even if the transmit array has an arbitrary planar structure. Note that different from the more traditional interpolation strategies applied used for receive array [31], [32], the interpolation in the context of TB MIMO radar needs to be performed as the transmit array. In addition to the above mentioned drawbacks, MIMO radar usually suffers from pulses storing condition. Comparing with the number of pulses, MIMO radar has a relative large dimension, which makes the estimation of the sample covariance matrix infeasible. Furthermore, performing a singular value decomposition (SVD) to a large scale covariance matrix demands a high computational burden.

In this paper, we develop a 2D transmit array interpolation for mono-static MIMO radar with irregular planar transmit array. We also develop a search-free 2D DOA estimation methods based on tensor modeling. To achieve high SNR gain within a known spatial sector/region of interest, we design a transmit array interpolation matrix limiting the interpolation errors between the actual steering vector and the desired one to a certain level within the sector of interest, while minimizing the sidelobes outside the sector. We have first reported the transmit array interpolation approach in [25] for the case of mapping an irregular transmit array into an L-shaped array. Here, we generalize the virtual array structure to several practically appealing structures, which satisfy the translational invariance property. The transmit array interpolation optimization problems are also generalized to arbitrary norm formulation to allow for smaller interpolation errors and sidelobes. With proper transmit array interpolation matrix design, computationally efficient methods can be applied to the 2D TB MIMO radar DOA estimation with an arbitrary array configuration. It is worth noting that, in addition, we do not require the receive array calibration information. Then we propose a general tensor modeling for the 2D TB MIMO radar with transmit array interpolation that fully capitalizes on its multidimensional structure. Under the small number of pulses condition, we consider the deterministic MIMO radar model, i.e., direct data approach, which requires significantly less computational burden. By performing HOSVD to the received MIMO radar snapshots, we achieve the tensor-based signal subspace and then develop tensor-based methods to obtain DOA estimates. Since the interpolation errors are not negligible in our model, we analyze the DOA estimation bias caused by the interpolation errors under high SNR condition. In addition, we explain how to build an offline look-up table to decrease the bias. Furthermore, we derive the corresponding transmit array interpolation 2D TB MIMO radar deterministic CRB as a benchmark for verifying the performance of the proposed methods.

The paper is organized as follows. The basic notations and MIMO radar signal model are given in Section II. In Section III, we develop the 2D TB MIMO radar interpolation method and present the corresponding direct data TB MIMO radar model, while the tensor modeling of 2D TB MIMO radar as well as matrix and tensor based ESPRIT type elevation and azimuth estimation methods are developed in Section IV. Section V is devoted to CRB derivation, the analysis of interpolation errors, and the development of a lookup table-based bias compensation technique in DOA estimation. Section VII reports simulation results showing the performance of the proposed approach under a variety of conditions, and Section VIII summarizes our conclusions. Some technical derivations of a supporting nature are given in Appendices.

II Notation and Signal Model

II-A Notation

Scalar is denoted by italic letters $a$, column vector by lower-case bold letter $\bm{a}$, matrix by upper-case bold letter $\bm{A}$, tensor by calligraphic bold-face letter $\bm{\mathcal{A}}$, respectively. The $i$-element of a vector is denoted as $a_{i}$, the $(i_{1},i_{2})$-element of a matrix $\bm{A}$ is denoted as $a_{i_{1},i_{2}}$, and the $(i_{1},i_{2},\ldots,i_{N})$-element of an $N$-order tensor $\bm{\mathcal{A}}\in\mathbb{C}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ as $a_{i_{1},i_{2},\ldots,i_{N}}$. We use the notations $(\cdot)^{*},(\cdot)^{T},(\cdot)^{H},(\cdot)^{-1},\circ,\odot,\otimes,\ast,\dagger$, and $\ddagger$ for denoting complex conjugate, transpose, Hermitian transpose, matrix inverse, outer product, Khatri-Rao product, Kronecker product, Hadamard product, pseudo-inverse, and projection, respectively. In addition, $\|\cdot\|_{p}$ and $\|\cdot\|_{\rm F}$ stand for the $p$-norm of a vector and the Frobenius norm of a matrix, respectively, $\mathbb{E}\{\cdot\}$ denoted the mathematical expectation, $\bm{I}_{M}$ is an $M\times M$ identity matrix, and $\text{diag}(\bm{a})$ denotes a diagonal matrix that holds the entries of $\bm{a}$ on its diagonal.

Definition 1 (The scalar product): The scalar product of two tensors $\bm{\mathcal{A}}$, $\bm{\mathcal{D}}\in\mathbb{C}^{I_{1}\times I_{2}\times\ldots I_{N}}$ is given by

[TABLE]

Using this definition, the higher-order norm of a tensor is given by

[TABLE]

A matrix unfolding of a tensor $\bm{\mathcal{A}}$ along $I_{n}$th mode is denoted as $\bm{\mathcal{A}}_{(I_{n})}$ following the notation in, for example, [17]. We define the concatenation of two tensors along $I_{n}$th mode as $\bm{\mathcal{A}}\sqcup_{I_{n}}\bm{\mathcal{D}}.$ By fixing $I_{n}$th mode to a specific value $k$, we obtain a sub-tensor denoted by $\bm{\mathcal{A}}_{I_{n}=k}$. Respectively, the vectorization of an $N$-order tensor $\bm{\mathcal{A}}$ along $I_{n}$th mode is defined as $\text{vec}(\bm{\mathcal{A}})\triangleq\text{vec}(\bm{\mathcal{A}}_{(I_{n})}^{T})$.

Definition 2 (The $n$-mode tensor-matrix product): The $n$-mode product of a tensor $\bm{\mathcal{A}}\in\mathbb{C}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ and a matrix $\bm{D}\in\mathbb{C}^{J\times I_{n}}$ along $n$th mode is given by

[TABLE]

where $\bm{\mathcal{C}}\in\mathbb{C}^{I_{1}\times I_{2}\times\ldots\times I_{n-1}\times J\times I_{n+1}\times\ldots\times I_{N}}$.

Definition 3 (The outer product): The outer product of $N$-order tensor $\bm{\mathcal{A}}\in\mathbb{C}^{I_{1}\times I_{2}\times\ldots I_{N}}$ and $M$-order tensor $\bm{\mathcal{D}}\in\mathbb{C}^{J_{1}\times J_{2}\times\ldots J_{M}}$ is given by

[TABLE]

II-B MIMO Radar Signal Model

Consider a mono-static MIMO radar system equipped with $M$ co-located transmit antenna elements and $N$ co-located receive antenna elements. The transmit and receive arrays are assumed to be placed on a plane and have arbitrary geometries. The receive array consists of antenna elements randomly selected from the transmit array. The transmit antenna elements are assumed to be located at the position $\bm{p}_{m}\triangleq[x_{m},\,y_{m}]^{T},m=1,2,\ldots,M$. Then the $M\times 1$ steering vector of the transmit array can be expressed as

[TABLE]

where $\bm{u}(\theta,\phi)\triangleq[\sin\theta\cos\phi,\,\sin\theta\sin\phi]^{T}$ represents the propagation vector, and $\theta$, $\phi$ are the elevation and azimuth, respectively. Receive antenna elements are randomly chosen out of $M$ elements in the transmit array. Similarly, the steering vector of the receive array can be then expressed as

[TABLE]

Let $s_{m}(t)$ be the complex envelope of the $m$th transmit signal where $t$ represents the fast time, and then $\{s_{m}(t)\}_{m=1}^{M}$ be a set of $M$ waveforms. Each waveform $s_{m}(t)$ has unit energy, and all waveforms are assumed to be orthogonal to each other during one pulse, i.e., $\int_{T}s_{m}(t)s_{m^{\prime}}^{*}(t)dt=\delta(m-m^{\prime})$, where $T$ is the radar pulse duration, $\delta(\cdot)$ denotes the Dirac delta function, and $L$ is the number of samples per pulse period. The signal radiated towards a spatial region of interest is therefore given by

[TABLE]

where $\bm{s}(t)\triangleq[s_{1}(t),\ldots,s_{M}(t)]^{T}$ and ${a}_{m}(\theta,\phi)$ is the $m$th element of $\bm{a}(\theta,\phi)$.

We assume that radar cross section (RCS) coefficients obey Swerling II model [33]. It means that the RCS coefficients remain unchanged during one pulse while vary for different pulses. There are $K$ targets to be located in the spatial sector of interest. Therefore, the received MIMO observation vector can be expressed as

[TABLE]

where $q$ represent the slow time index, $\beta_{k}(q)$ is the RCS coefficient of $k$th target with variance $\sigma_{\beta}^{2}$, and $\bm{n}(t,q)$ is the noise vector modeled as complex spatial and temporal white Gaussian process. Using the orthogonality property of the transmit waveforms, the received data vector corresponding to the $m$th waveform after matched-filtering can be obtained as

[TABLE]

where $\bm{y}_{m}(q)\in\mathbb{C}^{N\times 1}$, $\bm{z}(q)$ is the noise vector after matched-filtering whose covariance matrix is given by $\sigma^{2}_{n}\bm{I}_{M}$. Hence, the whole receive vector, i.e., the vector that is obtained by stacking $\bm{y}_{m}(q),\,m=1,\ldots,M$, one under another, can be written as

[TABLE]

where $\bm{A}(\theta,\phi)\triangleq[\bm{a}(\theta_{1},\phi_{1}),\ldots,\bm{a}(\theta_{K},\phi_{K})]$ is the transmit steering matrix, $\bm{B}(\theta,\phi)\triangleq[\bm{b}(\theta_{1},\phi_{1}),\ldots,\bm{b}(\theta_{K},\phi_{K})]$ is the receive steering matrix, and $\bm{\beta}(q)\triangleq[\beta_{1}(q),\ldots,\beta_{K}(q)]^{T}$ is the vector of RCS coefficients during $q$th pulse.

III 2D Transmit Array Interpolation

In many real world applications, we only need to consider a particular spatial sector instead of the whole space. The incident target spatial sector is thus assumed hereafter to be a prior knowledge. Then we can propose a 2D transmit array interpolation design aiming to focus the transmit energy within a pre-designed spatial sector while suppressing the sidelobes outside the sector of interest. In addition, it is needed to interpolate the actual steering vector of an arbitrary array into a virtual one which has a desired structure. In doing so, during the DOA estimation, we should be able to benefit from the use of the virtual array so that computationally efficient techniques can be applied and/or the higher-dimensional structure inherent in the 2D TB MIMO radar model can be exploited. Note that in this way, the initial waveform design and transmit array interpolation problems are decoupled, and the initial waveforms are multiplied by the desired interpolation matrix before being radiated from the antenna elements of the actual transmit array.

III-A Desirable 2D Virtual Transmit Array Structures

Let $\bm{E}\stackrel{{\scriptstyle\triangle}}{{=}}[\bm{e}_{1},\ldots,\bm{e}_{\tilde{M}}]$ denote a transmit beamspace and interpolation matrix, such that

[TABLE]

i.e., $\bm{a}(\theta,\phi)$ is mapped to $\bm{\tilde{a}}(\theta,\phi)$ that has a desired structure. A design technique for such interpolation matrix $\bm{E}$ can be developed in the same manner as the 2D transmit beamspace techniques [26], while enforcing additionally the condition (11). Here we discuss several desired structures of $\bm{\tilde{a}}(\theta,\phi)$.

If we map the transmit array steering vector to a coordinate product, i.e., Khatri-Rao product, of two arrays with $M_{1}$ and $M_{2}$ elements, respectively, the steering vector $\bm{\tilde{a}}(\theta,\phi)$, designed with enforcing (11), can be expressed as

[TABLE]

where

[TABLE]

Here, $\bm{c}(\mu)$ and $\bm{d}(\nu)$ represent the steering vectors of two virtual arrays, respectively. The parameters $\mu$ and $\nu$ in the steering vectors of the virtual arrays (12)–(14) are generalized direction cosine functions. In fact, the virtual arrays can have different structures depending on how we set the generalized parameters $\mu$ and $\nu$.

If we set $\mu=\sin\theta$ and $\nu=\sin\phi$, the virtual array (12) can be described as a Khatri-Rao product of two perpendicular uniform linear arrays (ULAs) placed along X-axis and Y-axis. In this case, the steering vectors of the vertical and horizontal ULAs carry the elevation and azimuth, respectively. The desired inter-element spacing has to be at least half-wavelength to avoid the ambiguity in angle estimation.

If we set $\mu=\sin\theta\cos\phi$ and $\nu=\sin\theta\sin\phi$, then the corresponding virtual array represents the URA.

If we map the transmit array into an L-shaped array, we have $\bm{\tilde{a}}(\theta,\phi)=[\bm{c}^{T}(\mu),\ \bm{d}^{T}(\nu)]^{T}$. Recalling (13) and (14), the virtual array in this case can be described as an L-shaped array placed on X-axis and Y-axis, which contains $M_{1}+M_{2}-1$ elements, i.e., one element is common for both the vertical and horizontal arrays.

Finally, the number of virtual antenna elements $\tilde{M}\leq M$ depends on which virtual array we design. For virtual two perpendicular/cross ULAs, we have $\tilde{M}=M_{1}+M_{2}$; for virtual URA, we have $\tilde{M}=M_{1}M_{2}$; and $\tilde{M}=M_{1}+M_{2}-1$ for virtual L-shaped array.

III-B Transmit Array Interpolation By Minimax Optimization

The spatial region of interest is the available prior knowledge Let $\theta_{g}\in\Theta,g=1,\ldots,G_{\theta}$ be the chosen angular grid which properly approximates the desired elevation sector $\Theta$ by a finite number $G_{\theta}$ of calibration points. Similarly, let $\phi_{g}\in\Phi,g=1,\ldots,G_{\phi}$ be the chosen angular grid which properly approximates the desired azimuth sector $\Phi$ by a finite number $G_{\phi}$ calibration points. The interpolation matrix ${\bm{E}}$ can be calculated as the least squares (LS) solution to (11). For this, we define the $M\times G_{\theta}\cdot G_{\phi}$ and the $\tilde{M}\times G_{\theta}\cdot G_{\phi}$ matrices ${\bm{A}}$ and $\tilde{\bm{A}}$, respectively, as follows

[TABLE]

[TABLE]

There are more calibration points than the antenna elements in the actual arrays ($G_{\theta}\cdot G_{\phi}\geq M$), thus $\bm{A}(\Theta,\Phi)$ and $\tilde{\bm{A}}(\Theta,\Phi)$ have full rank. The LS solution to (11) is given as

[TABLE]

The LS solution (17) provides the optimal approximation performance if we only concerned with minimizing the interpolation errors between the actual array and a desired one. However, for matrix $\bm{E}$ design in TB MIMO radar, LS solution (17) does not enable controlling the sidelobe levels of the transmit beampattern. Therefore, the resulting sidelobe levels can be even higher than the in-sector levels. This may result in wasting most of the transmit power in the out-of-sector areas which can lead to severe performance degradation.

To enable controlling sidelobe levels, we propose to use the minimax criterion to minimize the difference between the steering vectors of the actual array and the desired virtual array within the sector of interest, while keeping the sidelobe levels outside the spatial sector of interest bounded by a given maximum level. Therefore, the interpolation matrix design problem can be formulated as the minimax optimization problem as follows

[TABLE]

where $\bar{\Theta}$ and $\bar{\Phi}$ are the out-of-sector regions in the elevation and azimuth domains, respectively, $\theta_{h}\in\bar{\Theta},\ h=1,\ldots,H_{\theta}$ and $\phi_{h^{\prime}}\in\bar{\Phi},\ h^{\prime}=1,\ldots,H_{\phi}$ are angular grids which properly approximate the out-of-sector regions $\bar{\Theta}$ and $\bar{\Phi}$, respectively, $\gamma$ is a positive number of user choice used to upper-bound the worst-case sidelobe level, and any combination of $l_{1}$-, $l_{2}$-, or $l_{\infty}$-norms can be used in the objective and constraints of (18)–(19). It is worth noting that for $l_{2}$-norm, the optimization problem (18)–(19) becames a convex quadratically constrained quadratic programming (QCQP) problem. As such, it can be solved by interior-point methods [34] implemented, for example, in SeDuMi solver [35], which is used in the well known CVX MATLAB package. However, QCQP problems in general have heavy computational burden. Moreover and more significant, flat sidelobes over the elevation and azimuth cannot be achieved by (18)–(19) in the case of $l_{2}$-norm. Instead, if $l_{1}$- and $l_{\infty}$-norms are used in (18)–(19), the interpolation error for every value on the grid can be bounded, and sidelobes can be better controlled. In addition, for $l_{1}$- and $l_{\infty}$-norms, the optimization problem (18)–(19) can be cast as a linear programming (LP) problem, which can be efficiently solved by many existing LP solvers, including SeDuMi solver.

Alternatively, it is possible to minimize the worst-case out-of-sector sidelobe level while upper-bounding the norm of the difference between the interpolated array steering vector and the desired one. This can be formulated as the following optimization problem

[TABLE]

where $\Delta$ is a positive number of user choice used to control the deviation of the interpolated array from the desired one, which is also called interpolation errors tolerance. However, the interpolation matrix design based on solving the above optimization problem will introduce interpolation errors to the virtual transmit array, which will cause a bias in the DOA estimation. The relationship between the DOA estimation performance and the interpolation errors will be discussed in Section VI. Finally, note that the transmit array interpolation design does not involve waveform design. In fact, the problems of initial orthogonal waveforms design and transmit array interpolation are decoupled. Then the transmission process can be presented so that the initial orthogonal waveforms are first multiplied by the transmit beamspace and interpolation matrix $\bm{E}$ and then the resulting waveforms are radiated from the antenna element of the actual transmit array. In turn, the reception process remains standard, that is, matched filtering to the initial orthogonal waveforms.

III-C 2D TB MIMO Radar Signal Model

Similar to (7), the signal radiated in the 2D TB MIMO radar after transmit array interpolation (11) is given by

[TABLE]

where $\bm{e}_{\tilde{m}}$ is the $\tilde{m}$th column of the transmit beamspace and interpolation matrix $\bm{E}$, $\tilde{\bm{s}}(t)\triangleq[s_{1}(t),\ldots,s_{\tilde{M}}(t)]^{T}$, and $\{s_{\tilde{m}}(t)\}_{\tilde{m}=1}^{\tilde{M}}$ is a set of $\tilde{M}$ initial orthogonal waveforms, which can be for example a subset of the set $\{s_{m}(t)\}_{m=1}^{M}$ of a larger number of $M$ waveforms in (7).

Then the receive array observations can be expressed as

[TABLE]

Thus, after matched-filtering to the ${\tilde{m}}$th waveform, the ${\tilde{m}}$th $N\times 1$ received vector can be expressed as

[TABLE]

where $\bm{f}_{\tilde{m}}(\theta_{k},\phi_{k})\triangleq\tilde{\bm{a}}_{\tilde{m}}(\theta_{k},\phi_{k})\odot\bm{b}(\theta_{k},\phi_{k})$ and

[TABLE]

is the desired steering matrix. Hence, the whole received vector, i.e., the vector obtained by stacking $\tilde{\bm{y}}_{m}(q)$, $m=1,\ldots,M_{1}M_{2}$, one under another, can be expressed as

[TABLE]

where $\tilde{\bm{y}}(q)\in\mathbb{C}^{\tilde{M}N\times 1}$, and

[TABLE]

denotes the array steering matrix consisting of individual steering vectors. It is worth noting that (26) does not fully capture the multidimensional structure of the data hidden in the 2D TB MIMO radar model when several pulses are considered. Therefore, we develop in the next section a tensor modeling for 2D TB MIMO radar capturing several pulses.

IV Tensor Formulation Based DOA Estimation

IV-A 2D TB MIMO Radar Tensor Modeling

In this section, we show how to fold the 2D TB MIMO radar data into a tensor form suitable for applying tensor-based ESPRIT method [25] for the elevation and azimuth direction finding. Here, we operate directly on the received data using the so-called direct data approach.

Consider $Q$ pulses, then the received 2D TB MIMO radar data matrix is given as

[TABLE]

where $\tilde{\bm{Y}}\in\mathbb{C}^{\tilde{M}N\times Q}$, $\bm{P}\triangleq[\bm{\beta}(1),\bm{\beta}(2),\ldots,\bm{\beta}(Q)]$ contains $K$ vectors of targets’ RCS coefficients for $Q$ pulses, and $\bm{Z}$ represents the noise, which is assumed to be Gaussian with zero mean.

Recalling (II-A), the corresponding 2D TB MIMO radar tensor model can be expressed as

[TABLE]

where steering tensor $\bm{\mathcal{A}}$ is composed by stacking $K$ targets’ steering tensor $\bm{\mathcal{A}}_{k}$ together, $\bm{\mathcal{Z}}$ stands of the noise samples the same as that in (28), and $R$ represents the $R$th mode tensor-matrix product. The relation between the tensor and its matrix unfolding is given by

[TABLE]

It is worth noting that (29) represents a general tensor model for the 2D TB MIMO radar with transmit array interpolation. Some practically important special cases, which we will use throughout the paper, are the virtual URA and L-shaped arrays. For the case of virtual L-shaped array, $\bm{\mathcal{Y}}$ becomes a third-order tensor, and its dimensions represent virtual transmit array, receive array, and the temporal pulses, respectively. While for the case of virtual URA, $\bm{\mathcal{Y}}$ becomes a fourth-order tensor with two dimensions carrying elevation and other two dimentions carrying azimuth information. Also, the value of $R$ depends on which type of virtual transmit array we design. For the case of virtual L-shaped array, $R=3$, and for the case of virtual URA, $R=4$.

Thus, we have obtained the tensor form representation of the 2D TB MIMO radar data. This form allows to use some computational efficient methods, e.g., ESPRIT, for elevation and azimuth estimation. Comparing with searching-based DOA estimation methods, such method allows to reduce significantly the computational burden, especially when the model is of large size.

IV-B Steering Vector and Steering Tensor

For the case of virtual transmit L-shaped array in the 2D TB MIMO radar model with transmit array interpolation, the received data $\bm{\mathcal{Y}}\in\mathbb{C}^{(M_{1}+M_{2}-1)}$ is a third-order tensor. The derivations in this case are similar to that of the case of virtual URA. In addition, the translational invariance property holds for both the virtual transmit L-shaped array and URA in the same way. Therefore, we discuss in details hereafter the case when the transmit array is mapped into a virtual URA, while some details for the case of L-shaped array can be found in [25].

The relationship between a rank-one third-order tensor and its vector expression is given as follows

[TABLE]

The steering tensor of $k$th target can be expressed as

[TABLE]

where $\bm{\mathcal{A}}_{k}\in\mathbb{C}^{M_{1}\times M_{2}\times N}$. Substituting (32) into (31), we obtain the following relationship between the steering vector and the steering tensor

[TABLE]

It is clear that the array steering tensor obeys the translational invariance property with respect to $\bm{c}(\theta_{k})$ or $\bm{d}(\phi_{k})$, which makes tensor-based ESPRIT methods applicable for this scheme.

IV-C Tensor-Based Signal Subspace

Consider the case of virtual transmit URA as the most general example. In this case, the received tensor has four modes, $R=4$. The main motivation behind expressing the 2D TB MIMO radar model in the tensor form is that it allows for a more accurate elevation and azimuth estimation through tensor decomposition. The HOSVD of the tensor (29) is given by

[TABLE]

where $\bm{\mathcal{S}}\in\mathbb{C}^{M_{1}\times M_{2}\times N\times Q}$ represents the core tensor, and $\hat{\bm{U}}_{r},r=1,\ldots,4$ are unitary matrices. The tensor-based signal subspace can be then obtained by truncating noisy observation tensor $\bm{\mathcal{Y}}$ as

[TABLE]

where $\bm{\hat{\mathcal{S}}}^{[s]}\in\mathbb{C}^{K\times K}$ and $\hat{\bm{U}}_{r}^{[s]}\in\mathbb{C}^{M_{r}\times K}$ are the core tensor and $r$th mode signal subspace estimate, respectively. Then we can give the relationship between the signal subspace of noisy observation matrix $\tilde{\bm{Y}}$ and the matrix form of the HOSVD-based subspace. The SVD of $\tilde{\bm{Y}}$ can be expressed as

[TABLE]

where $\hat{\bm{U}}_{s}$ and $\hat{\bm{U}}_{n}$ are the matrices of the signal and noise subspaces estimates, respectively, while $\hat{\bm{\varSigma}}_{s}$ and $\hat{\bm{\varSigma}}_{n}$ are the corresponding diagonal matrices of singular values. Indeed, HOSVD is actually performing SVD along different unfoldings of a tensor. Thus, the SVD of the $r$th mode unfolding of the data tensor $\bm{\mathcal{Y}}$ can be expressed as

[TABLE]

where $\hat{\bm{U}}_{r}^{[s]}$, $\hat{\bm{U}}_{r}^{[n]}$, $\hat{\bm{\varSigma}}_{r}^{[s]}$, $\hat{\bm{\varSigma}}_{r}^{[n]}$ are the matrices composed of the signal and noise eigenvectors, and signal and noise eigenvalues, respectively, of the $r$th mode unfolding of $\bm{\mathcal{Y}}$. We have called $\hat{\bm{U}}_{r}^{[s]}$ previously as $r$th mode signal subspace estimate.

There is the following relationship between the steering tensor and the signal subspace

[TABLE]

where $\bm{T}$ denotes a nonsingular transform matrix.

The rotational invariance property holds for the tensor-based signal subspace, and it can be exploited for estimating both elevation and azimuth. If there are no interpolation errors in the mapping of the actual transmit array to a virtual one, we have

[TABLE]

where $\bm{J}_{1}^{(r)}$ and $\bm{J}_{2}^{(r)}$ stand for the selection matrices of each virtual array, and are given as

[TABLE]

for $r=1,2$. In (39) and (40), the rotational invariance property is satisfied in between the left and right-hand sides of the expressions.

Performing eigenvalue decomposition (EVD) of $\bm{\varTheta}_{r}$ and $\bm{\varPhi}_{r}$, we obtain

[TABLE]

where $\bm{T}_{1r}$ and $\bm{T}_{2r}$ are nonsingular (invertible) matrices since the targets are independent, $\bm{\varTheta}$ and $\bm{\varPhi}$ are the diagonal matrices carrying the elevation and azimuth, respectively, and given as

[TABLE]

IV-D HOSVD-Based DOA Estimation

The above developed 2D TB MIMO radar model has a similar form with the so-called $R$-dimensional ($R$-D) harmonic retrieval (HR) problem [17], [18]. The main difference is that the proposed model does not necessarily obey the Vandermonde structure at the receive array. However, it is not needed, since we can still solve the DOA estimation problem through taking advantage of the translational invariance property enforced at the transmit array while neglecting the structure of the receive array. If we exploit the whole signal subspace, the problems of finding $\bm{\varTheta}_{r}$ and $\bm{\varPhi}_{r}$ can be expressed in the form of the following two tensor-based LS problems

[TABLE]

These two tensor-based LS problems are equivalent to

[TABLE]

where

[TABLE]

for $i=1,2$ and $p=1,2$. Solving (49) and (50), yields,

[TABLE]

It has to be noticed however that since this method exploits the whole signal subspace when estimating elevation and azimuth, the DOA estimation errors can be enlarged if there exist interpolation errors in mapping the actual transmit array into a virtual one, especially in the high SNR region.

IV-E TEV-Based DOA Estimation

Another tensor decomposition method proposed in [18] and called as tensor approach based on eigenvectors (TEV) can be used as well for DOA estimation in our framework as an alternative to (53)–(54). In accordance with the TEV-based method, we divide the tensor-based signal subspace into several small parts along the targets’ dimension as

[TABLE]

For $k$th signal dimension, we have

[TABLE]

where $\bm{T}_{k}$ is a non-singular matrix. Recalling (32), we find out that $\bm{\mathcal{Q}}_{k}$ represents the signal subspace associated with $k$th signal dimension. Applying HOSVD to $\bm{\mathcal{Q}}_{k}$, yields,

[TABLE]

It can be seen that $\bm{q}_{k,r}$ spans the same subspace as the $r$th dimension of $\bm{\mathcal{A}}_{k}$. Using the linear prediction property [36], we can apply ESPRIT method to $\bm{q}_{k,1}$ and $\bm{q}_{k,2}$ when estimating $\mu_{k}$ and $\nu_{k}$ as

[TABLE]

The TEV-based method operates on each dimension of the signal subspace separately, which makes the proposed method more robust to interpolation errors because of the mapping of the actual transmit array into a virtual one. In addition, the auto-pairing has to be performed. It is technical, but otherwise straightforward, and it has been addressed, for example, in [37] as well as in a number of other papers.

V CRB and Influence of Interpolation Errors on DOA Estimation

V-A CRB

A useful statistical bound for evaluating the limiting DOA estimation performance in our framework is the deterministic CRB [38]. It can be derived for the considered 2D TB MIMO radar model with transmit array interpolation, i.e., model (26), in close-form. It is worth mentioning that the proposed 2D TB MIMO radar model with transmit array interpolation differs from the MIMO radar model given in (10) due to the fact that the interpolated transmit array is used in (26), while the original transmit array is used in (10). Thus, the CRB for the model (26) will be a function of the transmit interpolation matrix, and hence, the interpolation error will also be considered.

It is straightforward to see that the model (26) can be regarded as a special type of the deterministic $R$-D HR model as we discussed in Subsection IV.D. The corresponding closed-form CRB expression for the deterministic $R$-D HR model has been derived in [39]. In our case, the unknown parameter vector of interest is

[TABLE]

and the CRB for accuracy of the parameter vector (60) estimate is therefore can be written in our notations as

[TABLE]

where $\sigma_{n}^{2}$ is the noise variance and

[TABLE]

with $\bm{W}_{\rm te},\bm{W}_{\rm ta},\bm{W}_{\rm re},\bm{W}_{\rm ra}$ being the first order derivation of the steering matrix of the transmit array and receive array with respect to elevation and azimuth, respectively. Here we also use notations $\bm{D}$, $\bm{D}_{\rm e}$, $\bm{D}_{\rm a}$, $\bm{A}$, $\bm{F}$, $\bm{B}$ instead of, respectively, the full notations $\bm{D}(\theta,\phi)$, $\bm{D}_{\rm e}(\theta,\phi)$, $\bm{D}_{\rm a}(\theta,\phi)$, $\bm{A}(\theta,\phi)$, $\bm{F}(\theta,\phi)$, $\bm{B}(\theta,\phi)$ for brevity.

V-B DOA Estimation Bias

Obviously, the interpolation of the actual transmit array into a virtual one within a sector of interest by solving, for example, the minimax optimization problem (20)–(21) is not perfect. Here we investigate the influence of the interpolation errors to the DOA estimation accuracy.

The interpolation errors can be expressed as

[TABLE]

where $\Delta\bm{A}(\theta,\phi)$ is regarded as model errors and can be calculated once the transmit array interpolation matrix is designed. As a result, the interpolation errors will cause a bias in DOA estimations through a perturbation of the signal subspace. The DOA estimation performance degradation caused by imperfect interpolation can even dominate over the finite samples effect in high SNR region [40], [41].

Let us thus analyse the effect of the interpolation errors on DOA estimation from the subspace perturbation point of view which is suitable for the direct data approaches. For simplicity, we isolate model errors from finite sample effects by neglecting the observation noise, i.e., the SNR is assumed to be large enough. The SVD of the noiseless observation model is then given as

[TABLE]

where $\bm{U}_{s}$ and $\bm{V}_{s}$ represent, correspondingly, the left and right singular vectors associated with $K$ non-zero singular values, which form the diagonal of the diagonal submatrix $\bm{\Sigma}_{s}$, while $\bm{U}_{n}$ and $\bm{V}_{n}$ contain, respectively, the left and right singular vectors associated with the zero singular values.

We note that the matrix-based signal subspace is equivalent to the tensor-based signal subspace in the absence of noise [17]. Thus, we can use just the matrix-based signal subspace hereafter. Here we define the perturbations of direction cosine functions as

[TABLE]

where $\hat{\mu}_{k}$ and $\hat{\nu}_{k}$ denote $k$th estimated direction cosine functions, while $\mu_{k}$ and $\nu_{k}$ denote the $k$th true direction cosine functions without the effect due to the interpolation errors.

In [42], the perturbed DOA estimation performance analysis has been performed for the case of virtual ULA and 1D matrix-based ESPRIT method. Our case of the 2D TB MIMO radar model can be regarded as a multi-dimension extension of the model considered in [43]. Otherwise, such extension is straightforward, and using the results in [42], the expectation of mean-squared error (MSE) with the perturbations of direction cosine functions $\Delta\mu_{k}$ and $\Delta\nu_{k}$ for the proposed 2D TB MIMO radar model with array interpolation can be found as

[TABLE]

where $\bm{\alpha}_{kr}=\bm{p}_{kr}^{T}(\bm{\varOmega}_{1r}\bm{U}_{s})^{\dagger}(\bm{\varOmega}_{2r}/\lambda_{kr}-\bm{\varOmega}_{1r})\bm{U}_{n}\bm{U}_{n}^{H}$, $\bm{\beta}_{kr}=\tilde{\bm{A}}(\theta,\phi)^{\ddagger}\bm{U}_{s}\bm{q}_{kr}$, $\bm{p}_{kr}$ and $\bm{q}_{kr}$ are the left and right orthgonal eigenvectors of $\bm{\varTheta}$ and $\bm{\varPhi}$ associated with $k$th eigenvalue $\lambda_{k}$, respectively, and $r=1,2$. The expressions for $\bm{\varOmega}_{ir},i=1,2$ have been given in (51), and $\sigma_{\rm app}^{2}$ denotes the array interpolation (approximation) errors, which can be formally expressed at a point $(\theta,\phi)$ as

[TABLE]

Note that $\sigma_{\rm app}^{2}$ is different from the interpolation error tolerance $\Delta$ as we defined it in (21). In particular, $\sigma_{\rm app}^{2}$ represents the exact array interpolation error at $(\theta,\phi)$, while the interpolation error tolerance $\Delta$ represents the worst array interpolation condition.

Based on (70), (71), and (72), we can finally derive the expectation of MSEs of the elevation and azimuth estimation biases. The details of the derivations are given in Appendix A, and the resulting $\mathbb{E}\{(\Delta\theta_{k})^{2}\}_{\Delta\bm{A}(\theta,\phi)}$ and $\mathbb{E}\{(\Delta\phi_{k})^{2}\}_{\Delta\bm{A}(\theta,\phi)}$ are given as

[TABLE]

In addition, the relationships between the perturbation of direction cosine functions and DOA estimation biases are also given in Appendix B.

V-C Look-Up Table

As shown by (73) and (74), the interpolation errors cause undesirable DOA estimation bias. However, we fortunately can compensate the bias to some extent through building an off-line look-up table over the spatial sector of interest. The look-up table building procedure is summarized in Table I. Then this table can be used to map the DOAs estimated from noisy model to more precise values in order to obtain the corresponding accurate estimates.

VI Simulation Results

In this section, we evaluate the effectiveness of the proposed method in terms of simulations. Since the TB and interpolation matrix design can be done off-line, we compare computational burden of the matrix-based and tensor-based DOA estimation methods with virtual transmit URA. It is worth noting that the computational complexity for the $R$-D HR approach has been calculated in [17] and [18]. The major computational complexities of the DOA estimation methods are the SVD and HOSVD parts, which require the complexities of $k_{t}K\prod_{i=1}^{2}M_{i}Q$ and $4k_{t}K\prod_{i=1}^{2}M_{i}Q$ for direct data approach, respectively. Here $k_{t}$ is a constant that depends on a specific method used for performing SVD [44]. Therefore, the tensor-based DOA estimation method has a higher, by a factor of 4, computational burden than the matrix-based one. Also, the computational burden of the DOA estimation method used in the case of mapping to virtual transmit L-shaped array is smaller than that of the case of mapping to virtual transmit URA.

In all examples, the mono-static MIMO radar with 64 co-located transmit antenna elements and 8 receive antenna elements is assumed. The transmit array is an irregular planar array for which each transmit antenna element is placed near the regular grid of $8\times 8$ URA of size of $4\lambda$ for each column and row. To realize this structure, we firstly generate an URA of size $8\times 8$ of length $4\lambda$ for each column and row, and then add random displacements, which obey uniform distributions $[-\frac{\lambda}{4},\frac{\lambda}{4}]$, to the antenna elements’ positions. The receive array consists of the antenna elements which are randomly selected from the transmit array. The transmit interpolation technique expressed in terms of the optimization problem (20)–(21) is used to obtain virtual transmit array. In (20)–(21), the $l_{\infty}$-norm is used in the objective function and $l_{1}$-norm is used in the constraints. The problem (20)–(21) minimizes sidelobes on a grid outside of the sector of interest, while keeping satisfactory accuracy within the sector of interest for each point on a grid. Thus, the use of $l_{\infty}$-norm in the minimax objective function (20) further guarantees that the worst (highest) sidelobe level is minimized, while the use of $l_{1}$-norm in the constraints (21) provides the least physical displacement of the actual antenna elements when they are mapped to the virtual elements. It is just because for each point on a grid within the sector of interest the sparsest solution, i.e., the solution that requires the smallest physical displacement between the actual and virtual antenna elements, is preferred by using $l_{1}$-norm.

VI-A TB and Interpolation Matrix

In all examples of this subsection, the targets are assumed to be located within a known spatial sector. The spatial sector of interest is set to be $\Theta=[30^{\circ},40^{\circ}]$ and $\Phi=[65^{\circ},75^{\circ}]$. We also design a transition zone at each side of the elevation and azimuth domains of $20^{\circ}$ and $15^{\circ}$, respectively.

First, we show the transmit beampatterns obtained for the 2D TB MIMO radar designs with interpolation, and then investigate the interpolation errors within the desired spatial sector.

In Fig. LABEL:fig:URA, the resulting transmit beampatterns obtained for the 2D TB MIMO radar with interpolation are shown. The actual transmit array is mapped here to the virtual URA with 4 antenna elements in each row and column, and the distance between adjacent antenna elements in the virtual URA is set to half wavelength. The interpolation errors are $0.1$ and $0.01$ in Fig. LABEL:fig:URA(a) and Fig. LABEL:fig:URA(b), respectively. It can be observed that the sidelobes of Fig. LABEL:fig:URA(a) are much smaller than the sidelobes in the Fig. LABEL:fig:URA(b), while the main beam of Fig. LABEL:fig:URA(a) is sharper than that in Fig. LABEL:fig:URA(a). Both figures show that the transmit beam powers are distributed mainly and almost uniformly within the spatial sector of interest. The TB with smaller interpolation errors focus more power with the spatial sector of interest.

In one more example, we investigate the transmit beampattern when the actual transmit array is mapped to the desired transmit array with L-shaped and 4 virtual antenna elements in the row and the column of such virtual array. The interpolation errors are set as before to $0.1$ and $0.01$. Fig. LABEL:fig:L shows the corresponding beampatterns. It can be observed from the figure that by TB design we can suppress the sidelobes and focus the transmit power within the spatial sector of interest. However, the transmit power is not uniformly distributed over the spatial sector. This may lead to DOA estimation performance degradation as compared to the previously tested example where the actual array has been mapped to the virtual URA of a larger number of virtual elements. For smaller interpolation tolerance, the sidelobes become larger because of the tradeoff between the interpolation error and transmit beamforming performance. Indeed, a sufficiently high accuracy of transmit interpolation has to be guaranteed within the sector of interest. When interpolation errors are smaller within the sector of interest, the accuracy of transmit array interpolation is higher, but there are less degrees of freedom left for sidelobes control. Therefore, sidelobes go higher as it can be seen by comparing to each other subfigures (a) and (b) in Figs. LABEL:fig:URA and LABEL:fig:L.

In the last example of this subsection, we study the performance of interpolation errors versus elevation and azimuth within the spatial sector of interest. For this, we express the normalized interpolation error with respect to elevation and azimuth at $(\theta,\phi)$ as

[TABLE]

The transmit array with 64 colocated antenna elements is interpolated to a virtual URA of size $4\times 4$. The array interpolation tolerance is set as $\Delta=0.1$. Fig. 3 shows the graph for the normalized interpolation errors (75) as a function of elevation and azimuth. It can be seen from the figure that the interpolation error is constrained below a certain level, and it is distributed almost uniformly within the spatial sector of interest.

VI-B DOA Estimation Performance

We evaluate the DOA estimation performance in terms of the root MSE (RMSE), which is computed, for elevation as

[TABLE]

The RMSE expression for azimuth is the same as (76) where elevation is replaced with azimuth parameter. SNR is defined as $10\text{log}_{10}(\|\tilde{\bm{A}}\bm{B}\|_{\rm F}^{2}/\|\bm{Z}\|_{\rm F}^{2})$. Performance of the matrix-based and tensor-based ESPRIT and 2D spectral MUSIC methods for elevation and azimuth estimation comparing with deterministic CRB as well as probability of resolving two targets are investigated to examine the effectiveness of the proposed methods.

For all simulation examples in this subsection, the number of Monte Carlo trials for each experiment is set as 1000. Note that the number of pulses is smaller than the size of the 2D TB MIMO radar.

In the first example, we plot the CRB associated with the L-shaped array and URA to testify the performance of the two desired virtual transmit array structures. We set the interpolation error tolerance to $0.1$ for both virtual transmit array structures tested. The virtual L-shaped array has 4 elements along X-axis and 4 elements along Y-axis, while the size of the virtual URA is $4\times 4$. The spatial sector of interest is set as $\Theta=[30^{\circ},40^{\circ}]$ and $\Phi=[65^{\circ},75^{\circ}]$. It can be observed from Fig. 5 that the deterministic CRB on the DOA estimation for the case of virtual transmit URA is lower than that of the virtual L-shaped transmit array.

In the second example, we examine the influence of the interpolation errors $\sigma^{2}_{a}$ on the DOA estimation accuracy. Indeed, the matrix- and tensor-based ESPRIT algorithms are sensitive to interpolation errors since the interpolated transmit array satisfies the translational invariance property only approximatly. Here, the transmit array with 64 colocated antenna elements is interpolated to a $4\times 4$ virtual URA. For given $\Delta$, we obtain interpolation matrix through transmit array interpolation design, thus, the interpolation error $\sigma^{2}_{a}$ can be calculated using (72). The empirical DOA estimation error is compared to the one obtained based on the analytical expressions derived in Subsection V-B, and the results are shown in Fig. 5. The figure verifies that the empirical DOA estimation bias agrees with the analytically computed one. For example, if $\Delta=0.1$ in the TB and interpolation matrix design in (11), the corresponding interpolation error is $\sigma_{\rm app}^{2}=5\times 10^{-3}$. As a result, we have $\Delta\theta=0.1^{\circ}$ and $\Delta\phi=0.2^{\circ}$ as shown in Fig. 5.

By using the look-up table, the DOA estimation error can be reduced, especially in the high SNR region. Therefore, in our third example, we examine the elevation and azimuth RMSE performance for the proposed methods when the look-up table is also used for both single and multiple signal cases. The transmit array has 64 co-located antenna elements. The size of the virtual URA is $4\times 4$ with half wavelength displacement. The spatial sector of interest is set as $\Theta=[30^{\circ},\,40^{\circ}]$ and $\Phi=[65^{\circ},\,75^{\circ}]$. For the single target case, the DOA is set as $\phi=66^{\circ}$ and $\theta=34^{\circ}$, while for the multiple targets case, the DOAs are set as $\bm{\theta}=\{33^{\circ},\,39^{\circ}\}$, $\bm{\phi}=\{66^{\circ},\,71^{\circ}\}$. The corresponding RMSEs are plotted versus SNR in Figs. LABEL:fig:s1q6 and LABEL:fig:q8, respectively, for the matrix- and tensor-based ESPRIT and 2D spectral MUSIC methods. Their corresponding CRBs are also plotted in the figures. It can be seen in Fig. LABEL:fig:s1q6 that the matrix- and tensor-based ESPRIT methods have almost the same elevation and azimuth RMSE performance. In comparison, the 2D spectral MUSIC method has a better elevation and azimuth RMSE performance and it achieves CRB because it is independent on the interpolation errors. However, the 2D spectral MUSIC requires an exhaustive 2D spectral search over the whole spatial region, which has a significantly higher computational complexity than the matrix- and tensor-based ESPRIT. In the multiple signal case, it can be observed in Fig. LABEL:fig:q8(a) that the HOSVD-based ESPRIT and TEV-based method have slightly better elevation RMSE performance than that of the matrix-based ESPRIT method at SNR $<$10dB. Moreover, all three methods have almost the same RMSE performance above SNR $>$10dB. In addition, as can be seen in Fig. LABEL:fig:q8(b), all three methods have almost the same azimuth RMSE performance at SNR $<$10dB, TEV-based method has a much better azimuth RMSE performance than the other two methods when SNR $>$10dB, because it is least influenced by the interpolation error as we discussed in Subsection V-B and shown in Fig. 5. Moreover, the 2D spectral MUSIC method suffers from the small samples threshold effect. As a result, the 2D spectral MUSIC method shown poor elevation and azimuth RMSE performance, which gets closer to the one by the ESPRIT methods only for SNR$>$25dB and SNR$>$20dB for elevation and azimuth, respectively.

In the last experiment, we investigate the resolution ability of the 2D TB and interpolation MIMO radar with matrix- and tensor-based ESPRIT methods. The elevation and azimuth of two targets are set as $\bm{\theta}=\{36^{\circ},39^{\circ}\}$ and $\bm{\phi}=\{66^{\circ},69^{\circ}\}$, respectively. The transmit array has 64 co-located antenna elements. The size of the virtual URA is $4\times 4$. The spatial sector of interest is set as $\Theta=[30^{\circ},40^{\circ}]$, $\Phi=[65^{\circ},75^{\circ}]$. The targets are considered as resolvable if the following condition is satisfied for the elevation [45]

[TABLE]

where $\hat{\theta}_{p},p=1,2$ represent the estimated elevation parameters of two targets. This criterion is also suitable for testing the azimuth resolution ability. Fig. LABEL:fig:resolution shows that the resolution probability vesus SNR for matrix- and tensor-based ESPRIT methods. It can be seen that both methods achieve 100% correct resolution at the SNR of 25dB for elevation and at the SNR of 27dB for azimuth. Also, the tensor-based ESPRIT method has a higher resolution ability than that of the matrix-based method when SNR $<$25dB.

VII Conclusion

The problem of the 2D transmit array interpolation and beamspace design for mono-static MIMO radar with application to elevation and azimuth estimation has been addressed in this paper. The 2D transmit array interpolation and beamspace design problem has been formulated, for example, as the minimax convex optimization problem with constraints on array interpolation errors within a spatial sector of interest while minimizing the transmit power outside the sector. The desired structure of the virtual transmit array is also enforced. The virtual transmit array can be of any desirable shape, but practically most appealing shapes are the L-shaped array and URA because they allow to benefit from translational invariance property when estimating elevation and azimuth parameter at the receiver. By taking advantage of the high-dimensional structure inherent in the received data, the proposed MIMO radar model has been folded into a higher-order tensor. Then tensor-based ESPRIT methods with direct data approach have been devised. Furthermore, we have analysed the DOA estimation bias caused by transmit array interpolation errors and have built an offline look-up table aiming to decrease the DOA estimation bias. The corresponding deterministic CRB on DOA estimation in the proposed 2D TB and interpolation MIMO radar has been also provided. Simulation results have verified the effectiveness of the proposed methods.

[]

-A Interpolation Errors Analysis

We derive the MSE of $\mu_{k}$ as an example. In the absence of noise, the solution of (53) can be expressed with the matrix-based signal subspace as

[TABLE]

If there exists interpolation errors, we have

[TABLE]

Neglecting the cross item, $\Delta\bm{\Theta}$ can be calculated as

[TABLE]

where

[TABLE]

The first-order perturbation of $k$th eigenvector of $\bm{\Theta}$ is therefore given as

[TABLE]

Moreover, the relationship between $\Delta\mu_{k}$ and $\Delta\lambda_{k}$ can be obtained as

[TABLE]

where $c_{k}$ depends on the array structure. Based on the angle-root relationship [46], we have $c_{k}=1/\pi$ for URA.

Substituting (83) into (84), we finally obtain

[TABLE]

-B Conversion of Direction Cosine Functions to DOA Estimation Biases

The first-order partial derivatives of $\theta$ and $\phi$ with respect to $\mu$ and $\nu$ for URA are

[TABLE]

The multi-parameters Taylor series expansion formula is given as

[TABLE]

Based on (88), we can easily obtain the approximation of $\hat{\theta}$ and $\hat{\phi}$ around, respectively, $\theta$ and $\phi$ as

[TABLE]

Substituting (86) and (87) into (89) and (90), respectively, and neglecting the second-order and cross terms, we finally find the relationships between the perturbation of direction cosine functions and DOA estimation biases.