Massive MIMO Radar for Target Detection

TL;DR

This paper investigates the use of Massive MIMO technology in radar systems, demonstrating that increasing the number of antennas enhances target detection performance even without prior clutter knowledge.

Contribution

It develops a robust Wald-type detection test for Massive MIMO radar and derives asymptotic detection performance expressions.

Findings

Detection performance improves with more antennas.

Performance guarantees hold regardless of clutter statistics.

Asymptotic analysis matches finite system results.

Abstract

Since the seminal paper by Marzetta from 2010, the Massive MIMO paradigm in communication systems has changed from being a theoretical scaled-up version of MIMO, with an infinite number of antennas, to a practical technology. Its key concepts have been adopted in the 5G new radio standard and base stations, where $64$ fully-digital transceivers have been commercially deployed. Motivated by these recent developments, this paper considers a co-located MIMO radar with $M_T$ transmitting and $M_R$ receiving antennas and explores the potential benefits of having a large number of virtual spatial antenna channels $N=M_TM_R$. Particularly, we focus on the target detection problem and develop a \textit{robust} Wald-type test that guarantees certain detection performance, regardless of the unknown statistical characterization of the clutter disturbance. Closed-form expressions for the probabilities of false alarm and detection are derived for the asymptotic regime $N\to \infty$. Numerical results are used to validate the asymptotic analysis in the finite system regime with different disturbance models. Our results imply that there always exists a sufficient number of antennas for which the performance requirements are satisfied, without any a-priori knowledge of the clutter statistics. This is referred to as the Massive MIMO regime of the radar system.

Equations

$$\begin{array}[l]{clc}H_{0}:&\mathbf{x}_{k}=\mathbf{c}_{k}&k=1,\ldots,K,\\ H_{1}:&\mathbf{x}_{k}=\bar{\alpha}_{k}\mathbf{v}_{k}+\mathbf{c}_{k}&k=1,\ldots,K.\end{array}$$

$$\Lambda (\mathbf{X})\stackrel{H_1}{\underset{H_0}{\gtrless }}\lambda$$

$$\mathrm{Pr}\left\{\Lambda (\mathbf{X})>\lambda |H_0\right\}={\int }_{\lambda }^{\infty }{p}_{\Lambda |H_0}(a|H_0)da=\overline{P_{FA}},$$

$$\mathbf{x}(t)=\bar{\alpha }{\mathbf{a}}_R(\bar{\mathbf{\varphi }}){\mathbf{a}}_T^T(\bar{\mathbf{\varphi }})\mathbf{s}(t-\bar{\tau })e^{\mathsf{j}\bar{\omega }t}+\mathbf{n}(t),t\in [0,T],$$

$$\begin{array}{rl}\mathbf{X}& (l,k)=\bar{\alpha }{\mathbf{a}}_R(\bar{\mathbf{\varphi }}){\mathbf{a}}_T(\bar{\mathbf{\varphi }}{)}^T\mathbf{WS}(l,k)+\mathbf{C}(l,k),\end{array}$$

$$\mathbf{S}(l,k)\triangleq {\int }_0^T{\mathbf{s}}_o(t-\bar{\tau }){\mathbf{s}}_o^H(t-l\Delta t)e^{-j(k\Delta \omega -\bar{\omega })t}dt$$

$$\mathbf{C}(l,k)={\int }_0^T\mathbf{n}(t){\mathbf{s}}_o^H(t-l\Delta t)e^{-jk\Delta \omega t}dt.$$

$${\mathbb{C}}^N\ni \mathbf{x}=\mathrm{vec}\left(\mathbf{X}\right)=\bar{\alpha }\mathbf{v}(\bar{\mathbf{\varphi }})+\mathbf{c},$$

$$\mathbf{v}(\bar{\mathbf{\varphi }})=\left({\mathbf{S}}^T\otimes {\mathbf{I}}_{M_R}\right)\left[{\mathbf{W}}^T{\mathbf{a}}_T(\bar{\mathbf{\varphi }})\otimes {\mathbf{a}}_R(\bar{\mathbf{\varphi }})\right]$$

$$\mathbf{\Gamma }\triangleq E\{\mathbf{c}{\mathbf{c}}^H\}$$

$$\times {\left[{\mathbf{s}}_o^{\ast }(t-l\Delta t)\otimes {\mathbf{I}}_{M_R}\right]}^H{e}^{-jk\Delta \omega (t-\tau )}dtd\tau$$

$$={\int }_0^T{\int }_0^T\left[{\mathbf{s}}_o^{\ast }(t-l\Delta t){\mathbf{s}}_o^T(t-l\Delta t)\otimes \mathbf{\Sigma }(t-\tau )\right]\times$$

$$\times e^{-jk\Delta \omega (t-\tau )}dtd\tau .$$

$$\mathbf{c}{=}_d\sqrt{\tau }\mathbf{m}$$

$$\begin{array}[l]{clc}H_{0}:&\mathbf{x}=\mathbf{c},\\ H_{1}:&\mathbf{x}=\bar{\alpha}\mathbf{v}+\mathbf{c},\end{array}$$

$$\hat{\alpha }=\underset{\alpha \in \mathbb{C}}{\mathrm{argmin}}{G}_N(\mathbf{x},\alpha ),\mathbf{x}\sim p_X,$$

$$\begin{array}{rl}{G}_N& (\mathbf{x},\alpha )\triangleq \sum _{n=1}^N|x_n-\alpha v_n{|}^2=(\mathbf{x}-\alpha \mathbf{v}{)}^H(\mathbf{x}-\alpha \mathbf{v})\\ & =||\mathbf{x}|{|}^2+||\mathbf{v}|{|}^2{\left|\alpha -\frac{{\mathbf{v}}^H\mathbf{x}}{||\mathbf{v}|{|}^2}\right|}^2-\frac{|{\mathbf{v}}^H\mathbf{x}{|}^2}{||\mathbf{v}|{|}^2}\end{array}$$

$$\hat{\alpha }=\frac{{\mathbf{v}}^H\mathbf{x}}{||\mathbf{v}|{|}^2},$$

$$\hat{\alpha }\stackrel{a.s.}{\underset{N\to \infty }{\to }}\bar{\alpha }$$

$$\sqrt{N}{\bar{B}}_N^{-1/2}{A}_N\times (\hat{\alpha }-\bar{\alpha })\underset{N\to \infty }{\sim }\mathcal{C}\mathcal{N}(0,1),$$

$$A_N$$

$${\bar{B}}_N$$

$${\hat{B}}_N$$

$$+2N^{-1}\sum _{m=1}^l\sum _{n=m+1}^N\mathrm{Re}\left\{v_n{v}_{n-m}^{\ast }{\hat{c}}_n{\hat{c}}_{n-m}^{\ast }\right\}$$

$${\hat{c}}_n=x_n-\hat{\alpha }{v}_n,\forall n$$

$${\hat{B}}_N=N^{-1}{\mathbf{v}}^H{\hat{\mathbf{\Gamma }}}_l\mathbf{v},$$

$$[\widehat{\bm{\Gamma}}_{l}]_{i,j}\triangleq\left\{\begin{array}[]{cc}\hat{c}_{i}\hat{c}_{j}^{*}&j-i\leq l\\ \hat{c}_{i}^{*}\hat{c}_{j}&i-j\leq l\\ 0&|i-j|>l\end{array}\right.$$

$${\hat{B}}_N-{\bar{B}}_N\stackrel{p}{\underset{N\to \infty }{\to }}0.$$

$${\Lambda }_{\mathsf{RW}}(\mathbf{x})=\frac{2N|\hat{\alpha }{|}^2}{A_N^{-2}{\hat{B}}_N}=\frac{2|{\mathbf{v}}^H\mathbf{x}{|}^2}{{\mathbf{v}}^H{\hat{\mathbf{\Gamma }}}_l\mathbf{v}},$$

$${\Lambda }_{\mathsf{RW}}(\mathbf{x}|H_0)\underset{N\to \infty }{\sim }{\chi }_2^2(0),$$

$${\Lambda }_{\mathsf{RW}}(\mathbf{x}|H_1)\underset{N\to \infty }{\sim }{\chi }_2^2\left(\varsigma \right),$$

$$\varsigma =\frac{2|\bar{\alpha }{|}^2{M}_R^2\Vert (\mathbf{WS}{)}^T{\mathbf{a}}_T(\bar{\mathbf{\varphi }}){\Vert }^4}{\mathrm{tr}\left(\mathbf{\Gamma }\left[(\mathbf{WS}{)}^T{\mathbf{a}}_T(\bar{\mathbf{\varphi }}){\mathbf{a}}_T^H(\bar{\mathbf{\varphi }})(\mathbf{WS}{)}^{\ast }\otimes {\mathbf{a}}_R(\bar{\mathbf{\varphi }}){\mathbf{a}}_R^H(\bar{\mathbf{\varphi }})\right]\right)}.$$

$$\varsigma =\frac{2|\bar{\alpha }{|}^2{M}_R\Vert (\mathbf{WS}{)}^T{\mathbf{a}}_T(\bar{\mathbf{\varphi }}){\Vert }^2}{{\iint }_0^T||{\mathbf{s}}_o(t-\bar{l}\Delta t)|{|}^2\mathrm{tr}\left[\mathbf{\Sigma }(t-\tau )\right]e^{-j\bar{k}\Delta \omega (t-\tau )}dtd\tau },$$

Full text

I Introduction

Consider a multiple antenna radar system characterized by $N$ spatial channels collecting $K$ temporal snapshots $\{\mathbf{x}_{k}\}_{k=1}^{K}\in\mathbb{C}^{N}$ from a specific resolution cell, defined in an absolute reference frame. The primary goal of any radar system is to discriminate between two alternative hypotheses: the presence ($H_{1}$) or absence ($H_{0}$) of the target, in the resolution cell under test. Among others, a common model for the signal of interest is $\bar{\alpha}_{k}\mathbf{v}_{k}$, where $\mathbf{v}_{k}\in\mathbb{C}^{N}$ is known at each time instant $k\in\{1,\ldots,K\}$ and $\bar{\alpha}_{k}\in\mathbb{C}$ is a deterministic, but unknown, scalar that may vary over $k$. Any measurement process involves a certain amount of disturbance. In radar signal processing, the disturbance is produced by two components, the clutter and white Gaussian measurement noise, and it is modelled as an additive random vector, say $\mathbf{c}_{k}$, whose statistics may vary over $k$. Formally, the detection problem can be recast as a composite binary hypothesis test (HT) [1]:

[TABLE]

To solve (1), a decision statistic $\Lambda(\mathbf{X})$ of the dataset $\mathbf{X}\triangleq[\mathbf{x}_{1},\ldots,\mathbf{x}_{K}]$ is needed and its value must be compared with a threshold:

[TABLE]

to discriminate between the null hypothesis $H_{0}$ and the alternative $H_{1}$. A common requirement in radar applications is that the probability of false alarm has to be maintained below a pre-assigned value, say $\overline{P_{FA}}$. Consequently, the threshold $\lambda$ should be chosen to satisfy the following integral equation:

[TABLE]

where $p_{\Lambda|H_{0}}$ is the probability density function (pdf) of $\Lambda(\mathbf{X})$ under the null hypothesis $H_{0}$.

I-A Motivation

Finding a solution to (3) is in general a challenge. The common way out relies upon some “ad-hoc” assumptions on the statistical model of the dataset $\mathbf{X}$. In order to clarify this point, let us have a closer look to the steps required to solve (3). Firstly, a closed-form expression for $p_{\Lambda|H_{0}}$ is needed. By definition, $p_{\Lambda|H_{0}}$ is a function of the chosen decision statistic $\Lambda(\mathbf{X})$ and of the joint pdf $p_{\bm{X}}(\mathbf{X})$ of $\mathbf{X}$. If all the $\bar{\alpha}_{k},\forall k$ are modelled as deterministic unknown scalars, $p_{\bm{X}}(\mathbf{X})$ is fully determined by the joint pdf $p_{\bm{C}}(\mathbf{C})$ of the disturbance $\mathbf{C}=[\mathbf{c}_{1},\ldots,\mathbf{c}_{K}]$. A first simplification comes from the assumption that the disturbance vectors $\{\mathbf{c}_{k}\}$ are independent and identically distributed (i.i.d.) random vectors such that $p_{\bm{C}}(\mathbf{C})=\prod\nolimits_{k=1}^{K}p_{C}(\mathbf{c}_{k})$ [2, 3]. This assumption is, however, not always valid in practice. A second simplification that is commonly adopted in the radar literature (see e.g. [4, 2, 3]) is to assume that the functional form of $p_{C}(\mathbf{c}_{k})\equiv p_{C}(\mathbf{c}_{k};\bm{\gamma})$ is perfectly known, up to a possible (finite-dimensional) deterministic nuisance vector parameter $\bm{\gamma}$; for example, the (vectorized) covariance matrix. In order to obtain a consistent estimate $\hat{\bm{\gamma}}$ of $\bm{\gamma}$, a secondary dataset111In radar terminology, a secondary dataset is a set of “signal-free” snapshots collected form resolution cells adjacent to the one under test and sharing the same statistical characterization. has to be exploited (see e.g. [5]). Note that the required $p_{\Lambda|H_{0}}$ is a function of $\hat{\bm{\gamma}}$ as well. A third simplifying assumption is that the signal parameters $\bar{\alpha}_{k}$ remain constant over $k$, i.e. $\bar{\alpha}_{k}\equiv\bar{\alpha},\forall k$ [4, 2, 3]. Under these three assumptions, a possible choice for the decision statistic is the generalized likelihood ratio (GLR) $\Lambda_{\mathsf{GLR}}(\mathbf{X})$ (see e.g. [6, Ch. 11] and references therein). However, a closed-form solution to (3) can be found only for a very limited class of disturbance models for which the Gaussianity assumption needs also to be imposed. An asymptotic approximation for the solution to (3) can be obtained by exploiting a well-known asymptotic property of the GLR. Under the hypothesis $H_{0}$ and for $K\to\infty$, the pdf of $\Lambda_{\mathsf{GLR}}(\mathbf{X})$ converges to the one of a central $\chi$-squared random variable with $2$ degrees of freedom, denoted as $\chi_{2}^{2}(0)$ [6, Ch. 11]. Hence, by using the properties of the $\chi$-squared distribution, it is immediate to verify that (3) is asymptotically satisfied by $\bar{\lambda}=-2\ln\overline{P_{FA}}$. This is a particularly simple result that has received a lot of attention in the literature. However, it relies on the four simplifying assumptions previously introduced and summarized as follows:

• A1

The disturbance vectors $\{\mathbf{c}_{k}\}_{k=1}^{K}$ are i.i.d. over the observation interval.

• A2

The pdf $p_{C}(\mathbf{c}_{k};\bm{\gamma})$ of the disturbance is perfectly known, up to a unknown nuisance parameter vector $\bm{\gamma}$.

• A3

The target complex amplitude $\bar{\alpha}_{k}$ is maintained constant over the observation interval, i.e. $\bar{\alpha}_{k}\equiv\bar{\alpha},\forall k$.

• A4

The number of temporal snapshots $K$ is assumed to be much larger that the spatial channels $N$.

Even if these assumptions make the (asymptotic) analysis of $\Lambda_{\mathsf{GLR}}(\mathbf{X})$ analytically tractable, they are seldom satisfied in practical applications.

I-B Contributions

This paper considers a co-located MIMO radar with $M_{T}$ transmitting and $M_{R}$ receiving antennas and aims at deriving a detector that satisfies pre-assigned performance requirements without relying on the four assumptions above. Inspired by the recent developments in Massive MIMO communications [7, 8, 9, 10], we aim at exploring the potential benefits of having a very large number of antennas. Particularly, we assume that a single time snapshot, i.e. $K=1$, is collected, and operate in the asymptotic regime where the number of virtual spatial antenna channels $N=M_{T}M_{R}$ grows unboundedly, i.e., $N\to\infty$. This makes the three assumptions A1, A3 and A4 no longer needed. Advances in robust and misspecified statistics ([11, 12, 13, 14, 15, 16] and [17, 18]) are used to dispose of the cumbersome and unrealistic assumption A2. By adopting a very general disturbance model taking into account the spatial correlation structure of the observed samples, we propose a robust Wald-type detector that is asymptotically distributed, when $N\to\infty$, as a $\chi$-squared random variable (under both $H_{0}$ and $H_{1}$) irrespective of the actual and unknown disturbance pdf $p_{C}(\mathbf{c})$. This asymptotic result is achieved without the need of any secondary dataset. Although the theoretical findings of this paper are valid for a very general disturbance model, numerical results are provided for two non-Gaussian, stable auto-regressive disturbance models of order $p=3$ and $6$. It turns out that a pre-assigned value of $\overline{P_{FA}}=10^{-4}$ is achieved for $N=M_{T}M_{R}\geq 10^{4}$ with both models. This number of virtual spatial antenna channels defines what we call the Massive MIMO regime of the radar system.

Compared to our previous paper [19], the main difference lies in the absence of any a priori knowledge of the disturbance model. In fact, in [19] the analysis was developed for an autoregressive model of order $1$, but with no a priori knowledge of its statistics.

I-C Relevant literature

The MIMO paradigm has been the subject of intensive research over the past 15 years in radar signal processing. Initially introduced in wireless communications as a new enabling technology, the MIMO framework has been recognized to have a great potential in boosting the capabilities of classical antenna array systems. Based on the array configurations used, MIMO radars can be classified into two main types. The first type uses widely separated antennas (so-called distributed MIMO) to capture the spatial diversity of the target’s radar cross section (RCS) [20]. The second type employs arrays of closely spaced antennas (so-called co-located MIMO) to coherently combine the probing signals in certain points of the search area [21]. Hybrid configurations are also possible.

While the advantages in terms of spatial resolution, parameter identifiability, direction-of-arrival estimation and interference mitigation have been largely investigated in the MIMO radar literature, the potential benefits that a large number of virtual spatial antenna channels can bring into the target detection problem in terms of robustness with respect to the generally unknown disturbance model have not been explored yet. Surprisingly, not only the highly desirable robustness property has been somehow disregarded but, as pointed out in [22], even the availability of reliable, non-trivial, disturbance models is scarce. Remarkable exceptions to the mainstream Gaussianity assumption have been recently discussed in [23] and in [24]. Particularly, in [23] the performance of the Adaptive Normalized Matched Filter (ANMF), exploiting robust estimators for the disturbance covariance matrix, has been investigated with non-Gaussian disturbance. Specifically, random matrix tools have been used to obtain asymptotic approximations of the probabilities of false alarm and misdetection of the ANMF for the regime in which both $N$ and $K$ go to infinity with a non-trivial ratio $N/K$. Similar random matrix tools have been adopted in [24] to derive some asymptotic (in random matrix regime) results about the direction-of-departure and direction-of-arrival estimation in a non-Gaussian disturbance setting. Again, the random matrix machinery has been exploited in [25] to investigate the asymptotic performance of a GLRT detector in cognitive radio applications. Specifically, the HT problem tackled in [25] is similar to the one in (1), but the steering vector $\mathbf{v}_{k}$ and $\bar{\alpha}_{k}$ are assumed unknown. In the same spirit of [25], [26] has recently investigated the possibility to derive eigenvalues-based detectors for HT problem of the form in (1) for spectrum sensing and sharing in cognitive radio. However, in both [25] and [26], the disturbance was assumed to be a simple white Gaussian process with distribution a priori known, up to its statistical power. The asymptotic analysis in [24, 25, 26] requires that both $N$ and $K$ grow unboundedly. This is different from this paper where the temporal dimension $K$ is kept fixed; specifically, we assume to collect a single snapshot vector.

I-D Outline and Notation

The reminder of the paper is organized as follows. In Section II, the system and signal models as well as the resulting HT problem are introduced, by focusing our attention on the disturbance model for co-located MIMO radars. Section III describes our main results. Specifically, the explicit form of the proposed robust Wald-type test is provided along with the theoretical derivation of its asymptotic distribution under both $H_{0}$ and $H_{1}$. Numerical validations and simulation results are provided in Section IV. Some concluding remarks and discussions are drawn in Section V.

Throughout this paper, italics indicates scalar quantities ($a$), lower case and upper case boldfaces indicate column vectors ($\mathbf{a}$) and matrices ($\mathbf{A}$), respectively. Each entry of a $N\times N$ matrix $\mathbf{A}$ is indicated as $a_{i,j}\triangleq[\mathbf{A}]_{i,j}$, while the $i$-th column vector of $\mathbf{A}$ is indicated as $\mathbf{a}_{i}$ such that $\mathbf{A}=[\mathbf{a}_{1},\ldots,\mathbf{a}_{N}]$. We use ∗, T, and H to indicate complex conjugation, transpose, and the Hermitian operators, respectively. For random variables or vectors, $=_{d}$ stands for “has the same distribution as”. Also, $\overset{a.s.}{\underset{N\rightarrow\infty}{\rightarrow}}$ indicates the almost sure (a.s.) convergence and $\overset{p}{\underset{N\rightarrow\infty}{\rightarrow}}$ indicates the convergence in $p$-probability. We call $Q_{1}(\cdot,\cdot)$ the Marcum $Q$ function of order 1. The symbol $\lfloor a\rfloor$ defines the nearest integer less than or equal to $a\in\mathbb{R}$.

II System Model and Problem Formulation

Consider a co-located MIMO radar system equipped with $M_{T}$ transmitting antennas and $M_{R}$ receiving antennas [21]. The transmitting array is characterized by the array manifold, also called steering vector, $\mathbf{a}_{T}(\bm{\phi})$, where $\bm{\phi}$ is the position vector defined in an absolute reference frame [27]. Similarly, the receiving array can be characterized by the steering vector $\mathbf{a}_{R}(\bm{\phi})$ since the positions of the antennas in the absolute reference frame are known; see Fig. 1.

II-A Signal model

Given a target with position vector $\bar{\bm{\phi}}$, the signal collected at the receiving array can be modelled as [28, 22]:

[TABLE]

where $\mathbf{x}(t)\in\mathbb{C}^{M_{R}}$ is the array output vector at time $t$, $\mathbf{s}(t)\in\mathbb{C}^{M_{T}}$ is the vector of transmitted signals, $\bar{\alpha}\in\mathbb{C}$ accounts for the radar cross section of the target and the two-way path loss, which is the same for each transmitter and receiver pair. This is generally verified in co-located MIMO radars [21]. The parameters $\bar{\tau}$ and $\bar{\omega}$ represent the actual time delay and Doppler shift, due to the target position and velocity. The complex, vector-valued, random process $\mathbf{n}(t)\in\mathbb{C}^{M_{R}}$ accounts for the disturbance. We assume that $\mathbf{s}(t)$ is obtained as a linear transformation of a set of nearly orthonormal signals $\mathbf{s}_{o}(t)\in\mathbb{C}^{M_{T}}$, i.e. $\mathbf{s}(t)=\mathbf{W}\mathbf{s}_{o}(t)$, where $\mathbf{W}=[\mathbf{w}_{1},\ldots,\mathbf{w}_{M_{T}}]^{T}\in\mathbb{C}^{M_{T}\times M_{T}}$ and $\mathbf{w}_{m}\in\mathbb{C}^{M}$ is the weighting vector of the transmit antenna element $m$ with power $||\mathbf{w}_{m}||^{2}$.

Let $l=1,\ldots,L$ and $k=1,\ldots,K$ be the indices characterizing each time $l\Delta t$ and frequency $k\Delta\omega$ samples, respectively. The output $\mathbf{X}(l,k)\in\mathbb{C}^{M_{R}\times M_{T}}$ of the linear filter matched to $\mathbf{s}_{o}(t)$ can be expressed as [28, 22]:

[TABLE]

where

[TABLE]

takes into account potential “straddling losses”, that are losses due to a not precise centering of the target in a range-Doppler gate or to a not exact orthogonality between waveforms, and

[TABLE]

After omitting the indexes $(l,k)$ for ease of notation, we may rewrite (5) in vectorial form as:

[TABLE]

where $N\triangleq M_{R}M_{T}$ and

[TABLE]

and $\mathbf{c}\triangleq\mathrm{vec}\left(\mathbf{C}\right)$. While $\mathbf{a}_{T}(\bm{\phi})$ and $\mathbf{a}_{R}(\bm{\phi})$ depend on the geometry of the transmitting and receiving arrays, respectively, the matrix $\mathbf{W}$ can be designed to shape arbitrarily the transmitting beam; see [29, Ch. 4] and references therein. For example, with $\mathbf{W}=\mathbf{I}_{M_{T}}$, the transmitted power is uniformly distributed over all possible directions. On the other hand, with $\mathbf{W}=\mathbf{a}_{T}(\bm{\phi})^{*}\mathbf{a}_{T}(\bm{\phi})^{T}$, it is fully directed towards the direction $\bm{\phi}$. Intermediate cases can be obtained [28].

We assume that $\mathbf{n}(t)$ is zero-mean and wide-sense stationary; that is, $E\{\mathbf{n}(t)\}=\mathbf{0}$, $\forall t$ and $E\{\mathbf{n}(t)\mathbf{n}(\tau)^{H}\}=\bm{\Sigma}(t-\tau)$. Hence, from (7) it easily follows that $E\{\mathbf{c}\}=\mathbf{0}$ and

[TABLE]

As seen, $\bm{\Gamma}$ is a function of $\bm{\Sigma}(t-\tau)$ (i.e., the covariance matrix of $\mathbf{n}(t)$) and $\mathbf{s}_{o}(t)$. In the literature (e.g., [30],[31, Ch. 4]), a simple model for $\mathbf{n}(t)$ is to assume that its samples are uncorrelated in both spatial (along the receiving array) and temporal (along $T$) domains. This implies that $\bm{\Sigma}(t-\tau)=\sigma^{2}\mathbf{I}_{M_{R}}\delta(t-\tau)$. If $\mathbf{s}_{o}(t)$ is a vector of orthonormal waveforms, it is thus immediate to verify that $\bm{\Gamma}$ in (II-A) reduces to $\bm{\Gamma}=\sigma^{2}\mathbf{I}_{N}$. Under the assumption of uncorrelated samples in the time domain only and perfect orthonormality of the transmitted waveforms, we have that $\bm{\Gamma}=\mathbf{I}_{M_{T}}\otimes\bm{\Sigma}_{R}$ where $\bm{\Sigma}_{R}$ denotes the receive spatial covariance matrix [32]. However, the above two conditions may not be satisfied in practice [33, 22]. This is why in this paper we do not make any a priori assumption on the structure of $\bm{\Gamma}$. We only assume that its $(i,j)$-th entry goes to zero at least polynomially fast as $|i-j|$ increases. This assumption will be discussed in the next section and formally introduced in Assumption 1.

II-B Disturbance model

As previously discussed, many simplified models have been proposed in the literature to statistically characterize the disturbance vector $\mathbf{c}$ at the output of the matched filter bank of a (co-located) MIMO radar system. We refer to [33, 22] for a comparison among various statistical models and for a discussion about the physical simplifying assumptions underlying them. Here, we simply note here that the two main hypotheses usually made about the statistical characterization of the disturbance vector are: i) $\mathbf{c}$ is temporally and spatially white, and ii) $\mathbf{c}$ is Gaussian-distributed. As we will show, advances in robust and mis-specified statistics allow us to drop these two strong assumptions in favor of much weaker conditions.

To formally characterize the class of random processes to which the results of this paper apply, we need to introduce the concepts of uniform and strong mixing random sequences [34, 35, 36, 37]. Roughly speaking, the uniform and strong mixing properties characterize the dependence between two random variables extracted from a discrete process separated by $m$ lags. Without any claim of completeness or measure-theoretic rigor, the results of this paper apply to any random process that satisfies a restriction on the speed of decay of its auto-correlation function. Specifically, we limit ourselves to the following class of random processes:

Assumption 1**.**

Let $\{c_{n}:\forall n\}$ be a stationary discrete and circular complex-valued process [38] representing the true, and generally unknown, disturbance. Then, we assume that its autocorrelation function satisfies $r_{C}[m]\triangleq E\{c_{n}c^{*}_{n-m}\}=O(|m|^{-\gamma})$, $m\in\mathbb{Z}$, $\gamma>\varrho/(\varrho-1)$, $\varrho>1$.222Given a real-valued function $f(x)$ and a positive real-valued function $g(x)$, $f(x)=O(g(x))$ if and only if there exists a positive real number $a$ and a real number $x_{0}$ such that $|f(x)|\leq ag(x),\;\forall x\geq x_{0}$.

Assumption 1 implies that the volume of the MIMO radar must increase with the number of transmitting ($M_{T}$) and receiving ($M_{R}$) antennas. This means that the results of this paper do not apply to “space-constrained” array topologies. That said, Assumption 1 is very general and allows to account for most practical disturbance models. Indeed, any (not necessary Gaussian) stable second-order stationary (SOS) $\mathsf{ARMA}(p,q)$, and consequently any stable SOS $\mathsf{AR}(p)$ [39], satisfies Assumption 1, since the auto-correlation function of any stable SOS ARMA decays exponentially. The generality of the ARMA model is because it can approximate, for $p$ and $q$ sufficiently large, the second-order statistics of any complex discrete random processes having a continuous Power Spectral Density (PSD) [40, Ch. 3]. Moreover, a non necessarily Gaussian $\mathsf{ARMA}(p,q)$ is able to model the “spikiness” of heavy-tailed data as well. Another disturbance model of practical interest satisfying Assumption 1 is the Compound-Gaussian (CG) model [41, 42]. Indeed, recall that any CG-distributed random vector $\mathbf{c}$ admits a representation:

[TABLE]

for some real-valued positive random variable $\tau$, called texture, independent of the zero-mean, $N$-dimensional, circular, complex Gaussian random vector, called speckle, $\mathbf{m}\sim\mathcal{CN}(\mathbf{0},\bm{\Gamma})$, where $\bm{\Gamma}$ is its scatter matrix. Under a condition on the structure of $\bm{\Gamma}$, it easily follows that the $N$ entries of $\mathbf{m}$ can be interpreted as $N$ random variables extracted from a circular, Gaussian, SOS $\mathsf{ARMA}(p,q)$ process $\{m_{n}:\forall n\}$, with $p,q<N$.

We conclude by noticing that Assumption 1 is more general than the one adopted in our previous work [19]. In fact, the asymptotic results derived in [19] are obtained by assuming an AR disturbance model,333We considered only the $\mathsf{AR}(1)$ case, but the theory could be readily extended to a general $\mathsf{AR}(p)$ as discussed in [43]. driven by innovations with possibly unknown pdf. Unlike [19], this paper does not require any a priori information on the specific disturbance model; as stated in Assumption 1, only the polynomial decay of its auto-correlation function is needed.

II-C The hypothesis test problem

Based on the assumptions discussed above, the HT problem for target detection in (1) can be expressed as:

[TABLE]

where $\mathbf{c}\in\mathbb{C}^{N}$ is the disturbance vector whose entries are sampled from a complex random process $\{c_{n}:\forall n\}$ satisfying Assumption 1. Note that, in practical radar scenarios, (12) needs to be solved for any radar resolution cell of interest. Specifically, let $\{\bm{\phi}_{i};i=1,\ldots,Q\}$ be the set of position vectors pointing at $Q$ radar resolution cells, defined in an absolute reference frame. Then, the presence (or absence) of a target has to be tested for all the $Q$ resolution cells. Consequently, (12) has to be solved for any different vector $\mathbf{v}(\bm{\phi}_{i})$, whose explicit form is defined in (9). Notice also that a single-snapshot, i.e., $K=1$, is used in (12).

To discriminate between $H_{0}$ and $H_{1}$ in the composite HT problem (12), a test statistic is needed. Classical model-based test statistics such as the Generalized Likelihood Ratio (GLR) test and the Wald test (see e.g. [44]) cannot be used since the functional form of the disturbance pdf $p_{C}$ and, consequently, of the data pdf $p_{X}$ are unknown. We propose the following approach. Since no a priori information on the functional form of $p_{X}$ is available, let us choose the simplest estimator, say $\hat{\alpha}$, for the signal parameter $\bar{\alpha}$. Then, building upon the asymptotic statistics of $\hat{\alpha}$ as $N\rightarrow\infty$, we define a Wald-type test to solve the HT problem in (12). In fact, unlike the GLR statistic that requires the explicit functional form of the data pdf, a Wald-type statistic only needs an asymptotically normal, $\sqrt{N}$-consistent estimator of $\bar{\alpha}$ and a consistent estimate of its error covariance. As shown next, this key fact allows us to derive a robust test statistic for Massive MIMO radar configurations.

III Main Results

The main result of this section is the derivation of a robust Wald-type test for (12) with the valuable property to have, under $H_{0}$, an asymptotic distribution invariant w.r.t. $p_{X}$. The closed form expressions of its asymptotic distribution under both $H_{0}$ and $H_{1}$ will be provided. Since (12) is a composite HT problem (i.e., it depends on the unknown deterministic signal parameter $\bar{\alpha}$) a prerequisite for the implementation of a decision statistic is the derivation of an estimator $\hat{\alpha}$ of $\bar{\alpha}$.

III-A Estimation of $\bar{\alpha}$ under dependent data

By relying on the outcomes of [18, 17], in the following we show that an asymptotically normal, $\sqrt{N}$-consistent estimator of $\bar{\alpha}$ can be easily implemented under any general disturbance models satisfying Assumption 1.

A standard procedure is to recast an estimation problem into a relevant (possibly constrained) optimization problem. In the application at hand, an estimate of ${\alpha}$ can be obtained as:

[TABLE]

where $G_{N}(\cdot,\cdot)$ is a suitable objective function and $\mathbf{x}=[x_{1},\ldots,x_{N}]^{T}\in\mathbb{C}^{N}$ is the available dataset characterized by an unknown joint pdf $p_{X}$.

Given the measurement model (under $H_{1}$) in (12), the most natural choice for $G_{N}(\cdot,\cdot)$ is a Least-Square (LS) based objective function:

[TABLE]

whose minimum is reached when the second addend vanishes. This yields

[TABLE]

which is a well-known result in the radar signal processing literature addressing decision rules in Gaussian environment. It can also be noted that, under a misspecified (white) Gaussian assumption on the disturbance vector $\mathbf{c}$, the LS estimator in (15) coincides with the Mismatched Maximum Likelihood estimator [13],[18].

By specializing the general findings about misspecified estimation under dependent data proposed in [17] and [18], the asymptotic properties of the LS estimator in (15) can be obtained as shown in the following Theorem 1.

Theorem 1**.**

Under Assumption 1, the LS estimator $\hat{\alpha}$ in (15) is $\sqrt{N}$-consistent

[TABLE]

and asymptotically normal

[TABLE]

where

[TABLE]

and $\bm{\Gamma}\triangleq E_{p_{C}}\{\mathbf{c}\mathbf{c}^{H}\}$, with $p_{C}$ being the true (but generally unknown) disturbance pdf.

Proof:

All the required regularity conditions and a measure-theoretic rigorous proof (for the real-valued case) can be found in [17] and [18], while in Appendix A of this paper we provide the reader with an “easy-to-follow” but still insightful sketch of the proof in the complex case. Here, only the principal facts underlying the proof of Theorem 1 are discussed. To prove the consistency of $\hat{\alpha}$, we need an extension of the Law of Large Numbers (LLN) to uniform and strong mixing random sequences (see Assumption 1). This result is stated in Theorem 2.3 of [17]. With this suitable generalization of the LLN, the (strong) consistency of the LS estimator can be readily established as shown by Theorem 3.1 in [17]. Regarding the asymptotic normality, a generalization to uniform and strong mixing random sequences of the Central Limit Theorem (CLT) is required. This extension can be found in [17, Th. 2.4] for the scalar case, and in [45, Th. 2] for the multivariate case. The asymptotic normality of $\hat{\alpha}$ can be established by a direct application of this version of the CLT as shown in [17, Th. 3.2]. Finally, the extension of these results to the complex field can be obtained simply by exploiting the natural set isomorphism between $\mathbb{C}$ and $\mathbb{R}^{2}$ and by using the fact that only circular random sequences are considered. ∎

Remark 1**.**

Note that, since $\hat{\alpha}$ is obtained as a linear combination of circular observations ${x_{1},\ldots,x_{N}}$, its second-order statistics are fully characterized by the variance $E_{p_{X}}\{|\hat{\alpha}-\bar{\alpha}|^{2}\}$, while its pseudo-variance $E_{p_{X}}\{(\hat{\alpha}-\bar{\alpha})^{2}\}$ is nil [38].

While $A_{N}$ in (18) is only a function of the known norm of $\mathbf{v}$, $\bar{B}_{N}$ in (19) requires to compute the expectation w.r.t. the true, but unknown, disturbance pdf $p_{C}$. It can be shown that, under a uniform (or strong) mixing condition for the disturbance process $\{c_{n}:\forall n\}$, a consistent estimator of $\bar{B}_{N}$ is [17]:

[TABLE]

where $l$ is the so-called truncation lag [17],

[TABLE]

and $\hat{\alpha}$ is given in (15). The estimate $\hat{B}_{N}$ in (III-A) can be rewritten in a more compact form as:

[TABLE]

where

[TABLE]

for $1\leq i,j\leq N$. The consistency of $\hat{B}_{N}$ is stated in the next theorem (see [17, Th. 3.5]).

Theorem 2**.**

Under Assumption 1, if $l\rightarrow\infty$ as $N\rightarrow\infty$ such that $l=o(N^{1/3})$ then 444Given a real-valued function $f(x)$ and a strictly positive real-valued function $g(x)$, $f(x)=o(g(x))$ if for every positive real number $a$, there exists a real number $x_{0}$ such that $|f(x)|\leq ag(x),\;\forall x\geq x_{0}$.:

[TABLE]

Proof:

The proof is given in [46, Th. 6.20]. ∎

Theorem 2 provides us with a useful criterion to choose the truncation lag $l$. In fact, to ensure the consistency of $\hat{B}_{N}$, $l$ has to grow with $N$, but more slowly than $N^{1/3}$.

III-B A robust Wald-type test

Given the asymptotic characterization of the LS estimator presented in Theorem 1, the Wald-type statistic can be set up as:

[TABLE]

where the entries of $\widehat{\mathbf{\Gamma}}_{l}$ are given in (23). The similarity between the proposed $\Lambda_{\mathsf{RW}}$ statistic and the Adaptive Matched Filter (AMF) proposed in [47] is evident. However, the following comments are in order.

First, in [47] a set of homogeneous secondary snapshots, i.e. a set of “signal-free” data collected from radar resolution cells adjacent to the cell under test or at different time instants, is used to estimate the disturbance covariance matrix. This approach, however, requires that the disturbance statistics remain constant over all the considered resolution cells or time instants. Unlike the AMF, the decision statistic $\Lambda_{\mathsf{RW}}$ in (25) does not need any secondary data since it is able to extract all the required information from the single snapshot collected in the cell under test. Secondly, the AMF in [47] is derived under the Gaussian assumption for the disturbance vector $\mathbf{c}$. Conversely, $\Lambda_{\mathsf{RW}}$ in (25) can handle all the disturbance distributions that satisfy Assumption 1, including the Gaussian one. Third, as shown in Theorem 3 below, $\Lambda_{\mathsf{RW}}$ in (25) is shown to have the Constant False Alarm Rate (CFAR) property as $N\to\infty$ for all the disturbance distributions satisfying Assumption 1 and without the need of any secondary data. This is a great advantage w.r.t. the AMF that is CFAR only if the disturbance is Gaussian-distributed and a set of homogeneous secondary data is available.

The asymptotic property of $\Lambda_{\mathsf{RW}}(\mathbf{x})$ can be stated as follows.

Theorem 3**.**

If Assumption 1 hold true, then

[TABLE]

where $\varsigma\triangleq 2|\bar{\alpha}|^{2}\frac{||\mathbf{v}||^{4}}{\mathbf{v}^{H}\bm{\Gamma}\mathbf{v}}$.

Proof:

The proof follows from Theorem 1 and a known result about circular Gaussian random variables [48]. In particular, if $a\sim\mathcal{CN}(\mu_{a},\sigma_{a}^{2})$, then $2|a-\mu_{a}|^{2}/\sigma_{a}^{2}\sim\chi_{2}^{2}(0)$. Details are given in Appendix B. ∎

The above theorem shows that the pdfs of (25) under $H_{0}$ and $H_{1}$ converge to $\chi$-squared pdfs with $2$ degrees of freedom when the number of virtual spatial antenna channels $N=M_{T}M_{R}$ goes to infinity. This means that (3) is asymptotically satisfied by $\bar{\lambda}=-2\ln\overline{P_{FA}}$. This is valid for any pre-assigned $\overline{P_{FA}}$ and, more importantly, for any disturbance process $\{c_{n}:\forall n\}$ satisfying Assumption 1. In other words, we could say that $\Lambda_{\mathsf{RW}}(\mathbf{x})$ achieves the CFAR property w.r.t. all the disturbance distributions satisfying Assumption 1 when a sufficiently large number of transmitting and receiving antennas is used. Numerical results will be used to show that such Massive MIMO regime is achieved for a feasible large number of antennas. Some considerations on the asymptotic distribution of $\Lambda_{\mathsf{RW}}(\mathbf{x})$ under $H_{1}$ can also be done. In particular, in order to make explicit the dependence of $\varsigma$ from $\mathbf{a}_{T}(\bar{\bm{\phi}})$, $\mathbf{a}_{R}(\bar{\bm{\phi}})$ and $\mathbf{W}$, one can substitute the definition of $\mathbf{v}$ given in (9) into $\varsigma$ to obtain:

[TABLE]

Remark 2**.**

Further manipulations to the expression of $\varsigma$ in (28) are allowed if the model in (II-A) is adopted for the covariance matrix $\bm{\Gamma}$. Specifically, by substituting (II-A) in (28), we get:

[TABLE]

where $\bar{l}$ and $\bar{\omega}$ defines the range-Doppler gate under test. Moreover, if $\mathbf{s}_{o}(t)$ is a vector of perfectly orthonormal waveforms, i.e. $\mathbf{S}=\mathbf{I}_{M_{T}}$, and if $\bm{\Sigma}(t-\tau)=\sigma^{2}\mathbf{I}_{M_{R}}\delta(t-\tau)$, (29) can be further simplified as:

[TABLE]

where $P(\bm{\phi})\triangleq\mathbf{a}_{T}^{H}(\bm{\phi})\mathbf{W}^{*}\mathbf{W}^{T}\mathbf{a}_{T}(\bm{\phi})$ is the beam pattern of the transmitting array as a function of the position vector $\bm{\phi}$. The same result can be found in [30], where a LR test, perfectly matched to white Gaussian disturbance, was exploited as a detector.

The expression in (30) suggests us an interesting fact. Under the specific and simplistic scenario of perfectly orthonormal waveforms and (spatial and temporal) white Gaussian disturbance, a perfectly matched LR-based detector has the same (asymptotic) detection performance of the robust Wald-type test proposed here in (25). In other words, under the simplistic scenario mentioned above, $\Lambda_{\mathsf{RW}}$ in (25) does not present any loss in detection w.r.t. the “optimal” LR test even if $\Lambda_{\mathsf{RW}}$ does not require any a priori knowledge about the Gaussianity of the disturbance, while the LR test does. It is also worth mentioning that, in a more involved and realistic scenario, an LR-based test is practically unfeasible due to the lack of a priori information on the functional form of the disturbance pdf. Moreover, even if such disturbance pdf were available, the derivation of LR statistics would likely be met with the impossibility to derive its closed form expression. On the contrary, the closed form expression and the asymptotic detection performance of $\Lambda_{\mathsf{RW}}$ in (25) remain unchanged under any known or unknown disturbance process satisfying Assumption 1.

The following corollary is found.

Corollary 1**.**

If Assumption 1 holds true, then the detection probability of (25) is such that

[TABLE]

where $Q_{1}(\cdot,\cdot)$ is the Marcum Q function of order 1 [49] and $\varsigma$ is still given by $\varsigma=2|\bar{\alpha}|^{2}\frac{||\mathbf{v}||^{4}}{\mathbf{v}^{H}\bm{\Gamma}\mathbf{v}}$.

Proof:

By definition

[TABLE]

Then, (31) follows directly from (27) and the properties of the non-central $\chi^{2}$ distribution [49]. ∎

Since $Q_{1}(\cdot,\cdot)$ is monotonic in its first argument, Corollary 1 states that the $P_{D}$ of the RW test in (25) goes to 1 as $N\to\infty$. Moreover, it shows that $P_{D}$ depends on the true, but unknown, covariance matrix $\mathbf{\Gamma}$ of the disturbance vector $\mathbf{c}$, the radar geometry, and the waveform matrix $\mathbf{W}$ (through the vector $\mathbf{v}$).

IV Numerical Analysis

Numerical results are now used to validate the theoretical findings on the asymptotic properties of $\Lambda_{\mathsf{RW}}$ as stated in Theorem 3. We consider a uniform linear array at the transmitter and receiver, and a single target located in the far-field. We assume that $\mathbf{W}=\mathbf{I}_{M_{T}}$ and the transmitted waveforms are orthogonal, i.e., $\mathbf{S}=\mathbf{I}_{M_{T}}$. Following [50], we choose the radar geometry that maximizes the parameter identifiability. This is achieved by using a receiving array characterized by $M_{R}$ antenna elements with an inter-element spacing of $d$ and a transmitting array whose $M_{T}$ elements are spaced by $M_{R}d$. This implies that

[TABLE]

[TABLE]

where the spatial frequency $\nu$ is

[TABLE]

where $f_{0}$ is the carrier frequency of the transmitted signal, $c$ is the speed of light and $g(\cdot)$ is a known function of the position vector $\bm{\phi}$. By substituting (33) and (34) into (9), the vector $\mathbf{v}(\bm{\phi})\in\mathbb{C}^{N}$ can be expressed as, for $i=1,\ldots,N=M_{R}M_{T}$

[TABLE]

which represents the steering vector of an equivalent phased-array with $N$ elements [50].

Numerical results are obtained by averaging over $10^{6}$ Monte Carlo simulations. Moreover, the truncation lag in Theorem 2 is chosen as $l=\lfloor N^{1/4}\rfloor$.

IV-A Disturbance models

Two different models are considered for the disturbance. In Scenario 1, the disturbance vector $\mathbf{c}$ is generated according to an underlying circular, SOS $\mathsf{AR}(p)$

[TABLE]

with $p=3$, driven by i.i.d., $t$-distributed innovations $w_{n}$ whose pdf $p_{w}$ is [14]:

[TABLE]

where $\lambda\in(1,\infty)$ and $\eta=\lambda/(\sigma_{w}^{2}(\lambda-1))$ are the shape and scale parameters. Specifically, $\lambda$ controls the tails of $p_{w}$. If $\lambda$ is close to 1, then $p_{w}$ is heavy-tailed and highly non-Gaussian. On the other hand, if $\lambda\rightarrow\infty$, then $p_{w}$ collapses to the Gaussian distribution. In our simulations, we set $\lambda=2$ and $\sigma_{w}^{2}=1$. The $\mathsf{AR}(3)$ coefficient vector is $\bar{\bm{\rho}}=[0.5e^{\mathsf{j}2\pi 0},0.3e^{-\mathsf{j}2\pi 0.1},0.4e^{\mathsf{j}2\pi 0.01}]^{T}$. The normalized PSD can be expressed as can be expressed as:

[TABLE]

and is shown in Fig. 2. As seen, Scenario 1 is characterized by a disturbance whose power is mostly concentrated around the spatial frequency $\nu=0$.

To prove the robustness of the proposed $\Lambda_{\mathsf{RW}}$ w.r.t. more general disturbance models, in Scenario 2 we increase the order of the AR process generating the disturbance vector $\mathbf{c}$. Particularly, we consider a circular, SOS $\mathsf{AR}(6)$ driven (as before) by i.i.d., $t$-distributed innovations $w_{n}$ and characterized by the following coefficient vector $\bar{\bm{\rho}}=[0.5e^{-\mathsf{j}2\pi 0.4},0.6e^{-\mathsf{j}2\pi 0.2},0.7e^{\mathsf{j}2\pi 0},0.4e^{\mathsf{j}2\pi 0.1},0.5e^{\mathsf{j}2\pi 0.3},$ $0.6e^{\mathsf{j}2\pi 0.35}]^{T}$. The normalized PSD is reported in Fig. 3 and shows that, differently from Scenario 1, the disturbance power is spread over the whole range of $\nu$. Moreover, it presents more than a single peak.

In both scenarios, the disturbance process is normalized in order to have $\sigma^{2}=r_{C}[0]=1$. Under hypothesis $H_{1}$, the signal-to-noise ratio is simply defined as $\mathrm{SNR}\triangleq 10\log_{10}(|\bar{\alpha}|^{2}/\sigma^{2})$.

IV-B Performance Analysis

Since the disturbance PSD in Figs. 2 and 3 is not constant w.r.t the spatial frequency $\nu$, the performance of $\Lambda_{\mathsf{RW}}(\mathbf{x})$ will be evaluated for three different values of $\nu$ corresponding to different disturbance power density levels (i.e. three different target DOAs): $\nu_{1}=-0.2$, $\nu_{1}=0$ and $\nu_{1}=0.2$. Figs. 4 and 5 plot the $P_{FA}$ of the proposed robust Wald test in (25) for the three considered values of $\nu$. As we can see, in both scenarios $\Lambda_{{\mathsf{RW}}}$ in (25) achieves the nominal value of $10^{-4}$ for $N\geq 10^{4}$. Moreover, in this Massive MIMO regime, i.e. for $N\geq 10^{4}$, the progress of the $P_{FA}$ curves for different scenarios and for different values of spatial frequency are almost identical. This provides a numerical validation of the theoretical result provided in (26) of Theorem 3.

Fig. 6 considers Scenario 1 and illustrates the estimated and the closed-form expression of $P_{D}$ given in Corollary 1 for three distinct SNR values. Particularly, we assume $\mathrm{SNR}=-20,-10$ and $-5$ dB. As seen, the $P_{D}$ tends to $1$ as the number of virtual spatial antenna channels $N$ increases. Specifically, for $\mathrm{SNR}\geq-20$ dB the $P_{D}$ approaches $1$ for $N\geq 10^{4}$. From Fig. 6, it is also immediate to verify that the $P_{D}$ estimated through Monte-Carlo runs is in perfect agreement with the theoretical one provided in Corollary 1. We conclude by noticing that similar numerical results have been obtained for the CG disturbance model discussed in Sub-section II-B. Since they are perfectly in line with the numerical analysis reported above, we decided to not include them here due to lack of space.

V Concluding remarks and discussions

The detection problem in co-located MIMO radar systems was analysed in this paper. A robust Wald-type detector was proposed and its asymptotic performance investigated when the virtual spatial antenna channels $N=M_{T}M_{R}$ goes to infinity. Specifically, the CFAR property of the proposed detector for the asymptotic regime $N\to\infty$ and the wide family of disturbance processes satisfying Assumption 1 was mathematically proved and validated through numerical simulations. The purpose of analysing the asymptotic performance when $N\to\infty$ is not that we advocate the deployment of radars with a nearly infinite number of virtual antennas. The importance of asymptotics is instead what it tells us about practical systems with a finite number of antennas. Indeed, our main results imply that we can always satisfy performance requirements by deploying sufficiently many virtual antennas $N$, without any a priori knowledge of the disturbance statistics. Our numerical results showed that a pre-assigned value of $\overline{P_{FA}}=10^{-4}$ can be achieved with $N=M_{T}M_{R}\geq 10^{4}$ with non-Gaussian, stable autoregressive disturbance models of order $p=3$ and $6$. This defines the so-called Massive MIMO regime of the radar.

In Massive MIMO communications, linear combining, and precoding schemes can entirely eliminate the interference as the number of antennas grows unboundedly even with imperfect knowledge of propagation channels. We showed that a large-scale MIMO radar yields a target detector, which is robust to the unknown disturbance statistics. We foresee that further breakthroughs can be obtained by extending the Massive MIMO concept to other radar problems [10]. Clearly, by using very large arrays for waveform design, one can radically improve the spatial diversity gain and spatial resolution for target detection, parameter estimation, and interference rejection. However, the lesson learned from the last decade of research in communications is that Massive MIMO is not merely a system with many antennas but rather a paradigm shift with regards to the modelling, operation, theory, and implementation of MIMO systems. Our vision is that this paradigm can be applied also to radars, and open new research directions.

Appendix A Proof of Theorem 1

In this Appendix, the main ideas behind the proof of Theorem 1 are discussed. Specifically, we show the $\sqrt{N}$-consistency and the asymptotic normality of the mismatched LS estimator $\hat{\alpha}$ in (15) under Assumption 1.

A-A Consistency

Let us recall here the explicit expression of the LS objective function in (14): $G_{N}(\mathbf{x},\alpha)\triangleq\sum\nolimits_{n=1}^{N}|x_{n}-\alpha v_{n}|^{2}$. Moreover, let $\bar{G}_{N}(\alpha)$ be the following measurable and continuous function:

[TABLE]

Then, under Assumption 1, from [17, Th. 2.3] and [51, Th. 1], we have that

[TABLE]

where $\Omega$ is a compact subset of $\mathbb{C}(\equiv\mathbb{R}^{2})$. For the LS estimator $\hat{\alpha}$ in (15), the following convergence property holds.

Theorem 4**.**

If $\bar{G}_{N}(\alpha)$ has a unique minimum at $\alpha_{0,N}\in\mathbb{C}$, then from (41):

[TABLE]

and

[TABLE]

where $\alpha_{0,N}\triangleq\underset{\alpha\in\mathbb{C}}{\mathrm{argmin}}\;\bar{G}_{N}(\alpha)$.

Proof:

See the proof of Theorem 2.4 in [52]. ∎

Note that $\alpha_{0,N}$ represents the counterpart of the pseudo-true parameter defined in [12] for the i.i.d. data case. The consistency of $\hat{\alpha}$ can finally be established by showing that the pseudo-true parameter $\alpha_{0,N}$ is equal to the true one $\bar{\alpha}$. To this end, from (40) we obtain:

[TABLE]

By definition of data process in (12), $E_{p_{X}}\left\{x_{n}-\bar{\alpha}v_{n}\right\}=0,\;\forall n$, then the minimum of $\bar{Q}_{N}(\alpha)$ is attained at $\alpha_{0,N}=\bar{\alpha}$. This proves the consistency of $\hat{\alpha}$ under Assumption 1.

A-B Asymptotic normality

The proof for the asymptotic normality of the LS estimator $\hat{\alpha}$ mimics the one provided in standard statistical textbooks (see e.g. [53, Ch. 9]) for the Maximum Likelihood estimator. Let us start by taking the Taylor’s expansion of the real-valued objective function $G_{N}(\mathbf{x},\alpha)$ around the complex parameter $\bar{\alpha}$ [54, Th. 3.3]:

[TABLE]

where we used the fact that:

[TABLE]

By using the Mean Value Theorem [53, Ch. 9], for some $\tilde{\alpha}$ such that $|\tilde{\alpha}-\bar{\alpha}|<|\hat{\alpha}-\bar{\alpha}|$ we have:

[TABLE]

where $\hat{\alpha}$ is defined in (13). Consequently, we obtain

[TABLE]

For simplicity, we define

[TABLE]

[TABLE]

and rewrite (48) as

[TABLE]

Through direct calculation, it easily follows that

[TABLE]

where $c_{n}=x_{n}-\bar{\alpha}v_{n}$,$\forall n$. By substituting (52) and (53) into (51), we get:

[TABLE]

The asymptotic normality of $\hat{\alpha}$ follows from the application of the Central Limit Theorem established for (strong) mixing processes in [17, Th. 2.4] for the real scalar case and in [45] for the real multivariate case. In particular, let us define

[TABLE]

Under Assumption 1, by exploiting [17, Th.3.2] and recalling that $\{c_{n}:\forall n\}$ is a circular process, (17) follows.

Appendix B Asymptotic distribution of $\Lambda_{\mathsf{RW}}$

To establish the asymptotic distribution of $\Lambda_{\mathsf{RW}}$, we follow the standard procedure discussed in [53, Ch. 9].

B-A Asymptotic distribution of $\Lambda_{\mathsf{RW}}$ under $H_{0}$

We start by defining

[TABLE]

as the true parameter under the null hypothesis. Under $H_{0}$, from Theorem 1, we have that:

[TABLE]

[TABLE]

The inverse operator and the principal square root are both continuous operators on $\mathbb{R}^{+}$, then their composition is continuous on $\mathbb{R}^{+}$. Then, from (25) and by using the Continuous Mapping Theorem [55, Theo. 2.7], it follows that:

[TABLE]

Let us rewrite the test in (25) as:

[TABLE]

From (58) and (59), by a direct application of the Slutsky’s Theorem and of the properties of the complex Gaussian distribution [48], we immediately obtain that:

[TABLE]

where $\chi_{2}^{2}(0)$ indicates a central $\chi$-squared random variable with two degrees of freedom.

B-B Asymptotic distribution of $\Lambda_{\mathsf{RW}}$ under local alternatives

Following [53, Ch. 9], we suppose that the alternative to $H_{0}$ is of the form:

[TABLE]

Note that, as stated in [53, Sec. 9.3], this choice is made only for mathematical purposes, and has no direct physical significance. Specifically, (62) allows us to approximate the power of the test (or, in radar terminology, the probability of detection) locally, i.e. in the neighbourhood of the null hypothesis in (56). By defining $\bar{\alpha}_{H_{1}}^{(N)}\triangleq d/\sqrt{N}$ as the true parameter under local alternatives, we clearly have

[TABLE]

If $B_{N}(\bar{\alpha})$ in (55) is continuous in a neighbourhood of $\bar{\alpha}$, then from (63) and (56) it follows that

[TABLE]

Under $H_{1}$ in (62), from Theorem 1 we have

[TABLE]

[TABLE]

As before, by using the Continuous Mapping Theorem [55, Th. 2.7] and the limiting results in (64), we have

[TABLE]

By recasting the test in (25) as in (60), from (66) and (67), by a direct applications of the Slutsky’s Theorem, we immediately obtain that:

[TABLE]

Note that, from (55) we have $B_{N}(0)=N^{-1}\mathbf{v}^{H}\bm{\Gamma}\mathbf{v}$. This result can be used to approximate the power of the test, i.e. the $P_{D}$. By setting $d\equiv\sqrt{N}\bar{\alpha}$, we have:

[TABLE]

By using the properties of the non-central $\chi$-squared distribution with $2$ degrees of freedom [49], a closed form expression of the asymptotic probability of detection can be eventually expressed as:

[TABLE]