A Track-Before-Detect Strategy Based on Sparse Data Processing for Air Surveillance Radar Applications

Abstract

In this paper we consider the tracking problem of a moving target competing against noise and clutter in a surveillance radar scenario. For a single array-antenna multiple-target tracking system and according to the Track-Before-Detect paradigm, we present a novel approach based on a three-stage processing chain that involves the Sparse Learning via Iterative Minimization algorithm, the k-means clustering method and the ad hoc detector by exploiting the sparse nature of the operating scenario. Under the latter assumption, the detection strategy declares the presence of targets subsequently to the retrieval of their corresponding tracks performed by jointly processing the received echoes of multiple consecutive radar scans. Simulation results show that the proposed approach is able to provide good tracking and detection capabilities for different multiple target trajectories with low Signal-to-Interference-plus-Noise ratio and results in providing advantages when compared to a number of other reference Track-Before-Detect strategies based on sparse data processing techniques.

1. Introduction

In a multiple target tracking problem, the objective is to estimate an unknown and time varying number of targets as well as their trajectories from sensor data [1]. Traditional radar tracking algorithms are designed assuming that the sensor provides a set of point measurements at each scan, obtained by thresholding the output of a matched filter fed by a baseband version of the collected echoes [2,3,4]. This approach is effective in the limited case of Signal-to-Interference-plus-Noise ratio (SINR), while it experiences some limitations when the signal amplitude is weak compared to the background interference, and in general, in the detection of dim targets. Specifically, for dim targets the threshold must be low enough to allow for sufficient probability of target detection with the drawback that a low threshold leads to high volumes of false alarms.

An alternative approach, referred to as Track-Before-Detect (TBD), consists in supplying the processor with either thresholded or unthresholded data from several consecutive scans [5]. Unlike conventional techniques that declare the presence of a target at each scan of the radar, TBD-based procedures jointly process more consecutive scans and, exploiting the physically admissible target transitions, jointly declare the presence of a target and its track. Track-Before-Detect improves track accuracy and allows the tracker to effectively deal with low-observable targets, i.e., objects with a low SINR. The main difficulty in TBD techniques is that the measurements depend on the target state in a highly nonlinear way. Typically, the target state describes the kinematic evolution of the target and may also include its amplitude. A possible means to solve the nonlinear estimation problem is to resort to particle filtering [6]. An alternative is to discretize the target state space, since when the state is discrete, the linearity is irrelevant [7,8,9,10,11]. In this case, the problem with using a discrete state space is that it leads to high computation and memory resource requirements. Summarizing, if on one hand the TBD paradigm allows for a detection performance improvement, on the other hand, it leads high computational costs due to the huge amount of data collected at each scan. We refer the readers to [12] for a comparison of several TBD strategies to detect low amplitude targets.

The idea behind TBD has been developed in the past 20 years and its importance in low SINR radar application was first investigated in [13,14]. Since then, a number of TBD techniques and applications have been studied [15,16,17,18,19]. In recent years, security applications employing radar sensors for surveillance objectives are becoming increasingly important, so it is convenient to apply TBD techniques because targets with a low SINR appear [20,21,22,23,24]. In [21] the authors present a family of TBD procedures for early detection of moving targets from airborne radars exploiting a sectorization of the coverage area, the received echoes are jointly processed in the fast and slow time domains. In [22], TBD strategies for space-time adaptive processing radars are derived. The authors introduce the target and noise models in discrete-time form and a generalized likelihood ratio test strategy for target detection into two different scenarios. Results show that the proposed procedures might be a viable means to implement early detection and track initiation of weak moving targets. In [23], a number of new algorithms are proposed for adaptive detection and tracking based on spatial-time data. The possible spillover of target energy to adjacent range cells is taken into account at the design stage. The research proposed in [24] gives application of TBD processing against airborne radar surveillance data. The sensitivity of TBD processing to the choice of clutter model is examined through processing of real and simulated data containing non-coherent radar returns of a small maritime targets in sea clutter. The results confirm the good detection and tracking performance.

Leaving aside the TBD approach for the moment, we recall here that a radar system illuminates a given area and attempts to detect prospective targets present in its field of view by also estimating their locations and reflectivity. This task is accomplished by partitioning the region of interest into range, azimuth, and elevation bins. Remarkably, the prospective targets occupy only few bins as shown in Figure 1 and hence, if we neglect the interference, the 3D data matrix representing the observed scenario is sparse as most of its entries are approximately zero. Remarkably, it is important to highlight that considering multiple radar scans over time, the data sparsity increases even more. In [25], the authors exploit such features to devise radar processors based on compressed sensing techniques [26,27,28]. As a matter of fact, the exploitation of sparsity, originally proposed in statistics, signal processing, and machine learning communities, can improve the performance of radar systems under certain conditions and is therefore of considerable practical value.

Compressed sensing based methods have been applied to radar systems [29,30,31,32,33,34,35,36,37,38,39] to recover the target scene from a fewer number of measurements than traditional methods. In [29], it is demonstrated that compressed sensing can eliminate the need for matched filter at the receiver and has the potential to reduce the required sampling rate. Specifically, the authors present an adaptive clutter suppression method for airborne random pulse repetition interval radar by using prior knowledge of clutter boundary in Doppler spectrum. The research developed in [30] focuses on monostatic chaotic multiple-input-multiple-output radar systems and investigates the performance of sparsity-exploiting algorithms for the parameter estimation of targets with low-SINR. Remarkably, in the context of TBD, compressed sensing provides an effective means to address the case of multiple targets allowing to solve the nonlinear problem encountered in the traditional TBD algorithms. In [32] the authors propose an innovative, computational efficient TBD algorithm which exploits the sparse nature of measured data. They propose a two-step approach: the first step discards unreliable measurements, i.e., those with likelihood ratio below a preassigned thresholds; the second step, instead, uses the correlation among observations taken at different time instants to makes the final decision. In [33] a novel compressed sensing based TBD algorithm is devised. The latter reconstructs the whole radar scenario (direction of arrival (DOA)-Doppler plane) for each range gate at consecutive scans using an improved stagewise orthogonal matching pursuit (StOMP) algorithm, resulting in a three-dimensional range-DOA-Doppler space. The improved StOMP algorithm together with the temporal tracking, can effectively distinguish true targets from false targets and clutter based on information from multiple illuminations. In the context of synthetic aperture radar (SAR), refs. [34,35,36,37,38,39] present compressed sensing based data acquisition and imaging algorithms.

With the above remarks in mind, in this paper we propose a TBD strategy in the context of air surveillance radar that exploits the sparsity of measured data. The main technical contribution is represented by the development of a multiple target detection architecture comprising three processing stages using a monostatic array-antenna. More precisely, the first processing stage relies on the Sparse Learning via Iterative Minimization (SLIM) [40,41,42]. This choice is dictated by the fact that SLIM exhibits a good trade off between low computational cost and reconstruction performance, to estimate targets informations, i.e., position and speed, in the scene under consideration. More importantly, SLIM allows to avoid the nonlinear problem in the search space retrieval. The next step in the receiving chain is the formation of target tracks. To this end, the second stage is aimed at clustering the estimates provided by SLIM with each cluster representing a target trajectory. This stage connects the target information over time to create a trajectory. In the third and last stage, an ad hoc detector is derived resorting to a likelihood ratio test (LRT) that relies on the estimates provided by the clustering stage. The performance of the entire processing chain is assessed through simulated data accounting for different scenarios of practical interest. In addition, we compare the proposed technique with two compressed-sensing-based TBD algorithms exploiting the Orthogonal Matching Pursuit (OMP) [43] and the Compressive Sampling Matching Pursuit (CoSaMP) [44] algorithms, which represent the main benchmarks for compressive sensing recovery algorithms. Finally, the numerical analyses show the superiority of the proposed approach over the considered competitors.

The remainder of the paper is organized as follows. In Section 2, we formulate the TBD highlighting its sparse nature, then we provide a full description of the receiver architecture and the data processing framework. In Section 3 we assess the performance of the proposed architecture by means of numerical examples, while in Section 4 we analyse in detail the obtained results. Conclusions and future developments are provided in Section 5.

Notations

Vectors and matrices are denoted by boldface lower-case and upper-case letters, respectively. The $(k,l)$th entry of a generic matrix $\mathit{A}$ is denoted using $\mathit{A}(k,l)$ and, similarly, the lth entry of a generic vector $\mathit{a}$ with $\mathit{a}\left(l\right)$. Symbol ${(·)}^T$ denotes transpose operation, ${(·)}^{†}$ denotes complex conjugate transpose operation, $| · |$ denotes cardinality for a set while modulus for complex number, ⊗ denotes the Kronecker product, vec$(·)$ refers to the vectorization operation (e.g., stacking the columns of a matrix on top of each others), diag$(·)$ denotes the diagonalization operation, det$(·)$ and tr$(·)$ denote the determinant and the trace of a square matrix, respectively. $\mathit{I}$ and $\mathbf{0}$ represent the identity matrix and the null vector or matrix of suitable dimensions, respectively. As to the numerical sets, $\mathbb{R}$ is the set of real numbers, ${\mathbb{R}}^{N\times M}$ is the set of $(N\times M)$-dimensional real matrices (or vectors if $M=1$), $\mathbb{C}$ is the set of complex numbers, and ${\mathbb{C}}^{N\times M}$ is the set of $(N\times M)$-dimensional complex matrices (or vectors if $M=1$). Finally, the imaginary unit is denoted by j.

2. Material and Methods

2.1. Problem Formulation

The considered scenario involves a transmitter and a receiver, at the same location, equipped with an array of sensors, one or more point-like targets at a certain distance from the array in its far zone, and a narrowband signal that travels the round-trip between the radar and the target [45]. Specifically, we assume an uniform linear array of $N_s$ identical sensors with inter-element spacing h. The array illuminates the surveillance area by transmitting M coherent pulse trains, each consisting of $N_p$ rectangular pulses, before deciding whether or not a target is present. We denote by PRI the Pulse Repetition Interval, PRF $=1/$PRI the Pulse Repetition Frequency and SRT the Scan Repetition Time. The signal collected by the ith sensor is amplified, filtered, and down converted. Then, as customary, a discrete form for the signal received at the ith sensor is obtained by properly sampling the output of a filter matched to the transmitted pulse. For a fixed angle of view, the radar collects $N=N_s\times N_p$ samples from the rth range cell at mth scan; we denote by ${\mathit{z}}_{m,r}$ the received N-dimensional complex vector of noisy returns obtained by stacking up the $N_p$ vectors collected by the $N_s$ antenna sensors.

Let us assume that an unknown number L of point-like targets are present in the considered scenario at each scan with $L\in \{1,\cdots ,N_T\}$. Moreover, the number of range bins under test is $N_R\gg L$. In order to formulate the detection problem, we need to define $(r_{l_m},d_{l_m})$ for $l=1,\dots ,L$ and $m=1,\dots ,M$, that is the position of the lth target within the range-Doppler grid at the mth scan. The problem of detecting L multiple targets can be formulated in terms of the following multiple hypothesis test (1):

$$\left\{\begin{array}{cc}\hfill H_0:& {\mathit{z}}_{m,r}={\mathit{n}}_{m,r},r\in {\mathsf{\Omega }}_R,m=1,\dots ,M,\hfill \\ \hfill H_1:& \left\{\begin{array}{cc}{\mathit{z}}_{m,r}={\alpha }_{m,r,d}({\mathit{s}}_d\otimes \mathit{a})+{\mathit{n}}_{m,r,d},\hfill & (r,d)=(r_{1_m},d_{1_m})\in {\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D,m=1,\dots ,M,\hfill \\ {\mathit{z}}_{m,r}={\mathit{n}}_{m,r},\hfill & (r,d)\in {\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D\backslash \left\{(r_{1_m},d_{1_m})\right\},m=1,\dots ,M,\hfill \end{array}\right.\hfill \\ \hfill H_2:& \left\{\begin{array}{cc}{\mathit{z}}_{m,r}={\alpha }_{m,r,d}({\mathit{s}}_d\otimes \mathit{a})+{\mathit{n}}_{m,r,d},\hfill & (r,d)\in \{(r_{l_m},d_{l_m}),l=1,2\}\in {\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D,m=1,\dots ,M,\hfill \\ {\mathit{z}}_{m,r}={\mathit{n}}_{m,r},\hfill & (r,d)\in ({\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D)\backslash \{(r_{l_m},d_{l_m}),l=1,2\},m=1,\dots ,M,\hfill \end{array}\right.\hfill \\ \hfill ⋮\\ \hfill H_L:& \left\{\begin{array}{cc}{\mathit{z}}_{m,r}={\alpha }_{m,r,d}({\mathit{s}}_d\otimes \mathit{a})+{\mathit{n}}_{m,r,d},\hfill & (r,d)\in \{(r_{l_m},d_{l_m}),l=1,\dots ,L\}\in {\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D,m=1,\dots ,M,\hfill \\ {\mathit{z}}_{m,r}={\mathit{n}}_{m,r},\hfill & (r,d)\in {\mathsf{\Omega }}_R\times {\mathsf{\Omega }}_D\backslash \{(r_{l_m},d_{l_m}),l=1,\dots ,L\},m=1,\dots ,M,\hfill \end{array}\right.\hfill \end{array}\right.$$

(1)

where:

• ${\mathsf{\Omega }}_R=\{1,\dots ,N_R\}$ is the set of the range bins, and ${\mathsf{\Omega }}_D=\{1,\dots ,N_D\}$ indexes a discrete set of Doppler values covering the unambiguous Doppler interval;
• ${\alpha }_{m,r,d}\in \mathbb{C}$ is a factor representative of the target response and channel effects;
• ${\mathit{n}}_{m,r}\in {\mathbb{C}}^{N\times 1}$ for $m=1,\dots ,M$ and $r=1,\dots ,N_R,$ denote the additive interference components in the rth range bin at the mth scan, and are independent and identically distributed (iid) complex normal random vectors with zero mean and unknown positive definite Interference Covariance Matrix (ICM) ${\mathit{M}}_m\in {\mathbb{C}}^{N\times N}$;
• ${\mathit{s}}_d$ is the temporal steering vector defined as:

  $${\mathit{s}}_d={[1e^{j2\pi f_{d}PRI}\cdots e^{j2\pi f_{d}PRI(N_p-1)}]}^T\in {\mathbb{C}}^{N_p\times 1},$$

  (2)

  with $f_d=2v_d/\lambda$ and $v_d$ the prospective target radial velocity, positive if the target is approaching the radar;
• $\mathit{a}$ is the spatial steering vector defined as:

  $$\mathit{a}={[1e^{j2\pi \frac{h}{\lambda }cos{\theta }_0}\cdots e^{j2\pi \frac{h}{\lambda }cos{\theta }_0(N_a-1)}]}^T\in {\mathbb{C}}^{N_a\times 1},$$

  (3)

  with ${\theta }_0$ the fixed polar angle of target respect to the away direction. Figure 2 shows the radar geometry problem.

Thus, for a specific mth radar scan, the received signal matrix ${\mathit{Z}}_m\in {\mathbb{C}}^{N\times N_r}$ for the $N_R$ range bins and $N_D$ Doppler points can be written as:

$${\mathit{Z}}_m=\left(\sum _{r=1}^{N_R}\sum _{d=1}^{N_D}{\alpha }_{m,r,d}({\mathit{s}}_d\otimes \mathit{a}){\mathit{q}}_r^T\right)+{\mathit{N}}_m,$$

(4)

where ${\alpha }_{m,r,d}\ne 0$ if a target exists in the rth range bin and at the dth Doppler point, ${\mathit{N}}_m=[{\mathit{n}}_1,\cdots ,{\mathit{n}}_{N_R}]$ and ${\mathit{q}}_r$ represents the shift vector used to describe reflected waveforms from different range bins, and it is defined by:

$${\mathit{q}}_r={\left[\underset{r\mathrm{th}\mathrm{position}}{\underbrace{00\cdots 01}}0\cdots 0\right]}^T\in {\mathbb{C}}^{N_R\times 1}.$$

(5)

In order to recast (4) in vector form, we define:

$${\mathit{z}}_m=\mathrm{vec}\left({\mathit{Z}}_m\right)=\mathit{A}{\mathit{\alpha }}_m+{\mathit{n}}_m,$$

(6)

where $\mathit{A}$ is the space-time steering matrix defined as:

$$\mathit{A}=\left[{\mathit{v}}_{1,1}{\mathit{v}}_{1,2}\cdots {\mathit{v}}_{N_R,N_D}\right]\in {\mathbb{C}}^{N{N}_R\times N_R{N}_D},$$

(7)

with:

$${\mathit{v}}_{r,d}=\mathrm{vec}\left(({\mathit{s}}_d\otimes \mathit{a}){\mathit{q}}_r^T\right)\in {\mathbb{C}}^{N{N}_R\times 1},$$

(8)

${\mathit{\alpha }}_m$ is signal vector:

$${\mathit{\alpha }}_m={\left[{\alpha }_{m,1,1}{\alpha }_{m,1,2}\cdots {\alpha }_{m,N_R,N_D}\right]}^T\in {\mathbb{C}}^{N_R{N}_D\times 1},$$

(9)

and ${\mathit{n}}_m=\mathrm{vec}\left({\mathit{N}}_m\right)$ is a $N{N}_R\times 1$ complex normal random vector with zero mean and covariance matrix $\mathit{I}\otimes {\mathit{M}}_m$.

As $L\ll N_R{N}_D$, it follows that vector ${\mathit{\alpha }}_m$ in the linear model in (6) is sparse, i.e., most of its elements are zeros. Moreover, it is clear that with a limited number of measurement, sparsity constraints is required since, without it, the linear system has infinite solutions.

Furthermore, following the compact vector notation defined in (6), the problem (1) can be reformulated as follows:

$$\left\{\begin{array}{cc}\hfill H_0:& {\mathit{z}}_m\sim \mathcal{CN}(\mathbf{0},\mathit{I}\otimes {\mathit{M}}_m),m=1,\dots ,M,\hfill \\ \hfill H_1:& {\mathit{z}}_m\sim \mathcal{CN}(\mathit{A}{\mathit{\alpha }}_m,\mathit{I}\otimes {\mathit{M}}_m),m=1,\dots ,M,\hfill \end{array}\right.$$

(10)

where ${\mathit{\alpha }}_m$ accounts for the unknown number of targets.

2.2. Proposed TBD Framework

The detection architecture is composed of three stages and its schematic overview is shown in Figure 3. Specifically, the first stage is aimed at estimating the sparse vector ${\mathit{\alpha }}_m$ exploiting the SLIM algorithm, the second stage is responsible for the track formation and is based on the k-means clustering method, and, finally, the third stage implements an ad hoc detector, obtained through the LRT approach.

2.2.1. Sparse Learning

As stated before, the main objective of the first stage is to estimate the sparse vector ${\mathit{\alpha }}_m$ defined in the linear model (6). To this end, we apply the SLIM framework [40,42] assuming that ${\mathit{\alpha }}_m$ is a random complex vector ruled by a sparsity promoting prior [46,47]. It follows that:

$${\mathit{z}}_m|{\mathit{\alpha }}_m\sim \mathcal{CN}(\mathit{A}{\mathit{\alpha }}_m,\mathit{I}\otimes {\mathit{M}}_m),$$

(11)

whereas the prior of ${\mathit{\alpha }}_m$ is:

$$f({\mathit{\alpha }}_m;q)\propto \prod _{i=1}^{N_R{N}_D}{e}^{-\frac{2}{q}(\mid {\alpha }_{m,i}{\mid }^q-1)},$$

(12)

where $q\in (0,1]$ is a parameter controlling the sparsity level. In fact, when $q=1$, the prior becomes Laplacian which has a finite peak at 0 while when $q\to 0$ it has an infinite peak at 0. Note, also, ${\mathit{M}}_m$ is not known and, hence, must be estimated from data. Specifically, we assume that a set of training data collected in the spatial proximity of the range under test is available and use them to obtain the following ICM estimate:

$${\hat{\mathit{M}}}_m={\displaystyle \frac{1}{K}}\sum _{k=1}^K{\mathit{z}}_{m,k}{\mathit{z}}_{m,k}^{†},$$

(13)

where ${\mathit{z}}_{m,1},\cdots ,{\mathit{z}}_{m,K}$ are the training data at the mth scan. Thus, we can estimate ${\alpha }_m$ as follows:

$$\underset{{\mathit{\alpha }}_m}{max}f({\mathit{z}}_m;\mathit{I}\otimes {\hat{\mathit{M}}}_m|{\mathit{\alpha }}_m)f({\mathit{\alpha }}_m;q),$$

(14)

where:

$$f({\mathit{z}}_m;\mathit{I}\otimes {\hat{\mathit{M}}}_m|{\mathit{\alpha }}_m)=\frac{1}{{\pi }^{N{N}_R}\det (\mathit{I}\otimes {\hat{\mathit{M}}}_m)}{e}^{[-{({\mathit{z}}_m-\mathit{A}{\mathit{\alpha }}_m)}^{†}{(\mathit{I}\otimes {\hat{\mathit{M}}}_m)}^{-1}({\mathit{z}}_m-\mathit{A}{\mathit{\alpha }}_m)]},$$

(15)

is the conditional probability density function (pdf) of ${\mathit{z}}_m$ given ${\mathit{\alpha }}_m$. By taking the negative logarithm it is equivalent to

$$\underset{{\mathit{\alpha }}_m}{min}\underset{g_q\left({\mathit{\alpha }}_m\right)}{\underbrace{\left\{{∥{\mathit{y}}_m-\mathit{V}{\mathit{\alpha }}_m∥}_2^2+\sum _{i=1}^{N_R{N}_D}\frac{2}{q}(\mid {\alpha }_{m,i}{\mid }^q-1)\right\}}},$$

(16)

where ${\mathit{y}}_m={\left(\mathit{I}\otimes {\hat{\mathit{M}}}_m\right)}^{-1/2}{\mathit{z}}_m$ and $\mathit{V}={\left(\mathit{I}\otimes {\hat{\mathit{M}}}_m\right)}^{-1/2}\mathit{A}$. Note that the the first addendum of $g_q\left({\mathit{\alpha }}_m\right)$ corresponds to a fitting term, whereas the second term promotes the sparsity.

Setting to zero the first derivative of $g_q\left({\mathit{\alpha }}_m\right)$ with respect to ${\mathit{\alpha }}_m$ leads to:

$$\frac{d}{dx}\left[g_q\left({\mathit{\alpha }}_m\right)\right]={\mathit{V}}^{†}\mathit{V}{\mathit{\alpha }}_m-{\mathit{V}}^{†}{\mathit{y}}_m+{\mathit{P}}_q^{-1}{\mathit{\alpha }}_m=0,$$

(17)

where ${\mathit{P}}_q=$ diag$\left({\mathit{p}}_q\right)$, with ${\mathit{p}}_q={[\mid {\alpha }_{m,1}{\mid }^{2-q},\mid {\alpha }_{m,2}{\mid }^{2-q},\cdots ,\mid {\alpha }_{m,N_r{N}_d}{\mid }^{2-q}]}^T$. Supposing that an initial estimate of ${\mathit{\alpha }}_m$ is available, it is possible to apply a cyclic optimization procedure as in [40,42], and the step at the nth iteration can be expressed as:

$${\mathit{\alpha }}_{m,q}^{\left(n\right)}={\mathit{P}}_q^{(n-1)}{\mathit{V}}^{†}{\left(\mathit{V}{\mathit{P}}_q^{(n-1)}{\mathit{V}}^{†}+\mathit{I}\right)}^{-1}{\mathit{y}}_m,$$

(18)

where ${\mathit{P}}_q^{(n-1)}=$ diag$\left({\mathit{p}}_q^{n-1}\right)$ comes from the $(n-1)$th iteration. The optimization procedure can terminate after a fixed number of iterations or when the following convergence criterion is satisfied:

$$\frac{{∥{\mathit{\alpha }}_{m,q}^{\left(n\right)}-{\mathit{\alpha }}_{m,q}^{(n-1)}∥}_2}{{∥{\mathit{\alpha }}_{m,q}^{\left(n\right)}∥}_2}<{\mathsf{\Delta }}_{SLIM},$$

(19)

with ${\mathsf{\Delta }}_{SLIM}$ a user defined suitable small positive number to obtain the desired estimation accuracy. Specifically, as described in [40] SLIM algorithms converges quickly and in most cases it shows no improvements after 15 iterations. For the SLIM initialization, we apply the same condition based on matched filter as defined in [40,42]:

$${\mathit{\alpha }}_{m,i}^{\left(0\right)}=\frac{{\mathit{v}}_i^{†}{\left(\mathit{I}\otimes {\hat{\mathit{M}}}_m\right)}^{-1}{\mathit{y}}_m}{{\mathit{v}}_i^{†}{\left(\mathit{I}\otimes {\hat{\mathit{M}}}_m\right)}^{-1}{\mathit{v}}_i},$$

(20)

for $i=1,2,\dots ,N_r{N}_d$ and ${\mathit{v}}_i$ is the ith column of $\mathit{A}$. To automate the estimation of q, so that SLIM can essentially become user parameter free, the Bayesian Information Criterion (BIC) [48], a model selection tool, is applied [40]:

$$BI{C}_q=2{∥{\mathit{y}}_m-\mathit{V}{\mathit{\alpha }}_m∥}_2+4h\left(q\right)\times ln\left(N{N}_R\right),$$

(21)

where $h\left(q\right)$ represents the number of peaks in the SLIM result for a particular value of q. The first term in (21) represents a fitting term for the estimate ${\widehat{\mathit{\alpha }}}_m$, the second term denotes a penalty term, which increases as the number of peaks in the estimate increases and hence is larger for values of q approaching 1. The coefficient 4 in (21) represents the number of unknowns target parameters: complex amplitude, range bin and Doppler shift. It is important to observe that from a practical point of view, after this first stage a data buffer is necessary to temporarily store all the estimated signal vector during the M scans so that it collects $\mathsf{\Gamma }=[{\widehat{\mathit{\alpha }}}_1,\cdots ,{\widehat{\mathit{\alpha }}}_M]\in {\mathbb{C}}^{N{N}_R\times M}$.

2.2.2. Track Formation

The second processing stage is devoted to the extraction of target trajectories from the estimated matrix $\mathsf{\Gamma }$. This task is accomplished by means of two consecutive steps. The first one concerns the estimation of the target number (i.e., L). This value is automatically computed supposing that, for a specific radar scan, a target cannot be located in more than one range-Doppler grid. Therefore, L is computed by searching column-wise the maximum number of non-zero element in $\mathsf{\Gamma }$. More precisely, let $Y_m$ the set of non-zero element for the estimated vector ${\widehat{\mathit{\alpha }}}_m$:

$$Y_m=\{{\widehat{\alpha }}_{m,r,d}\in {\widehat{\mathit{\alpha }}}_m:{\widehat{\alpha }}_{m,r,d}\ne 0,m=1,\dots ,M\},$$

(22)

therefore, the estimation of maximum target number is computed as:

$$L=\underset{m=1,\dots ,M}{max}\left(|Y_m|\right)$$

(23)

The target to track association is carried out by the application of a clustering algorithm to a dataset containing the information about target positions during radar scanning. Specifically, the k-means clustering algorithm, which can be considered one of the simplest and popular unsupervised machine learning algorithms, is applied [49,50,51]. This algorithm works iteratively, starting from an initial dataset and given the number of clusters, to produce dataset partitions representing the tracks. The iterative refinement clustering is based on a well-defined objective function, also called error-function. In this context, the k-means algorithm converges at least to local optimum (i.e., not necessarily the best possible outcome) [51].

For the case at hand, the number of clusters corresponds to the number of targets L, and the dataset, namely X, to be clustered is formed starting from $\mathsf{\Gamma }$ matrix, as follows:

$$X=\left\{\bigcup _{m,r,d}{\mathit{x}}_{m,r,d}\right\},$$

(24)

where ${\mathit{x}}_{m,r,d}={[m,r,d,{\widehat{\alpha }}_{m,r,d}]}^T$, for all ${\widehat{\alpha }}_{m,r,d}\in Y_m$, with $m=1,\dots ,M$. The k-means algorithm starts with initial estimates for the centroids, which in general can either be randomly generated or randomly selected from the dataset, then it iterates between two steps:

• data assignment step: each centroid defines one cluster, therefore each vector data is assigned to its nearest centroid. This association in based on the Euclidean metric which is defined as objective function to be locally minimized:

  $$\underset{{\mu }^{\left(k\right)}}{argmin}\sum _{l=1}^{\left|X\right|}\sum _{k=1}^L{∥\mathit{W}{\mathit{x}}_l-{\mu }^{\left(k\right)}∥}_2^2\Rightarrow {\mathit{x}}_l={\mathit{x}}_l^{\left(k\right)}:X^{\left(k\right)}=\left\{\bigcup _k{\mathit{x}}^{\left(k\right)}\right\},$$

  (25)

  where $\mathit{W}=$ diag$\left(\mathit{w}\right)$ is a weight matrix, ${\mu }^{\left(k\right)}$ represents the kth cluster centroid; ${\mathit{x}}_l^{\left(k\right)}$ stands for the assignment of the generic vector ${\mathit{x}}_l$ to k-cluster, and $X^{\left(k\right)}$ indicates the kth cluster;
• centroid update step: the centroids are recomputed by taking the mean of all vector data assigned to that cluster centroid:

  $${\mu }_k=\frac{1}{|X^{\left(k\right)}|}\sum _{{\mathit{x}}^{\left(k\right)}\in X^{\left(k\right)}}^{}{\mathit{x}}^{\left(k\right)}.$$

  (26)

The above steps are repeated until a stopping criteria is met (i.e., no data points change clusters, the sum of the distances is minimized, or some maximum number of iterations is reached). The output of interest of the k-means algorithm, that is also the output for this processing stage, is the final L-partitioning of the starting dataset: $X=\{X^{\left(1\right)},\cdots ,X^{\left(L\right)}\}$. Each partition contains the range-Doppler point coordinates, that combined in time through the radar scanning information they make the target track.

2.2.3. Ad Hoc Detector

The last stage performs the validation of the tracks formed in the previous stages by comparing a suitable statistics exploiting such tracks to a threshold that allows to control the probability of false track detection. To this end, we resort to the LRT where the unknown parameters are replaced by the previously obtained estimated. Specifically, such decision scheme has the following expression:

$$\frac{f_1({\mathit{Z}}_1,\cdots ,{\mathit{Z}}_M;\mathsf{\Gamma },{\hat{\mathsf{\Omega }}}_1^{\left(L\right)},\cdots ,{\hat{\mathsf{\Omega }}}_M^{\left(L\right)},{\hat{\mathit{M}}}_1,\cdots ,{\hat{\mathit{M}}}_M)}{f_0({\mathit{Z}}_1,\cdots ,{\mathit{Z}}_M;{\hat{\mathsf{\Omega }}}_1^{\left(L\right)},\cdots ,{\hat{\mathsf{\Omega }}}_M^{\left(L\right)},{\hat{\mathit{M}}}_1,\cdots ,{\hat{\mathit{M}}}_M)}\underset{H_0}{\stackrel{H_1}{\gtrless }}\eta ,$$

(27)

where:

$$\begin{array}{c}{f}_1({\mathit{Z}}_1,\cdots ,{\mathit{Z}}_M;\mathsf{\Gamma },{\hat{\mathsf{\Omega }}}_1^{\left(L\right)},\cdots ,{\hat{\mathsf{\Omega }}}_M^{\left(L\right)},{\hat{\mathit{M}}}_1,\cdots ,{\hat{\mathit{M}}}_M)={\displaystyle \prod _{m=1}^M}{\displaystyle \prod _{(r,d)\in {\hat{\mathsf{\Omega }}}_m^{\left(L\right)}}}f({\mathit{z}}_{m,r};{\widehat{\alpha }}_{m,r,d}({\mathit{s}}_d\otimes \mathit{a}),{\hat{\mathit{M}}}_{\mathit{m}}),\\ f_0({\mathit{Z}}_1,\cdots ,{\mathit{Z}}_M;{\hat{\mathsf{\Omega }}}_1^{\left(L\right)},\cdots ,{\hat{\mathsf{\Omega }}}_M^{\left(L\right)},{\hat{\mathit{M}}}_1,\cdots ,{\hat{\mathit{M}}}_M)={\displaystyle \prod _{m=1}^M}{\displaystyle \prod _{(r,d)\in {\hat{\mathsf{\Omega }}}_m^{\left(L\right)}}}f({\mathit{z}}_{m,r};\mathbf{0},{\hat{\mathit{M}}}_{\mathit{m}}),\end{array}$$

with ${\hat{\mathsf{\Omega }}}_m^{\left(L\right)}$ the set of integers indexing the range bins and Doppler points corresponding to the estimated target position at the mth scan, and

$$\begin{array}{c}f({\mathit{z}}_m;{\widehat{\alpha }}_{m,r,d}({\mathit{s}}_d\otimes \mathit{a}),{\hat{\mathit{M}}}_{\mathit{m}})=\frac{1}{{\pi }^N\det \left({\hat{\mathit{M}}}_m\right)}{e}^{\{-{[{\mathit{z}}_{m,r}-{\widehat{\alpha }}_{m,r,d}({\mathit{s}}_d\otimes \mathit{a})]}^{†}{\hat{\mathit{M}}}_m^{-1}[{\mathit{z}}_{m,r}-{\widehat{\alpha }}_{m,r,d}({\mathit{s}}_d\otimes \mathit{a})]\}},\\ f({\mathit{z}}_m;\mathbf{0},{\hat{\mathit{M}}}_{\mathit{m}})=\frac{1}{{\pi }^N\det \left({\hat{\mathit{M}}}_m\right)}{e}^{\{-{\mathit{z}}_{m,r}^{†}{\hat{\mathit{M}}}_m^{-1}{\mathit{z}}_{m,r}\}}.\end{array}$$

It is interesting to note that the LRT under the noise only assumption is readily obtained from the above formulations as well. Once the likelihoods of the states are evaluated, a state estimate can be obtained according to the ratio test. In the event that this state is the null state, then the detector reports that there are no targets, otherwise a detection is reported.

3. Results

Numerical simulation are performed following the specifications of the Lockheed Martin FPS-117 surveillance radar [52]. FPS-117 is a 3-D solid-state radar designed to provide long-range accurate aircraft identification and position data for air defence, navigational assistance and tactical control for both close air support and counter air operations [53]. All simulations refer to a radar system with $N_s=4$, $N_p=4$, carrier frequency $f_c=1.388$ GHz, PRF equal to 917 Hz and SRT $=0.68$ s calculated considering a mechanical rotation scanning with an horizontally array-antenna with a 3 dB beamwidth equal to $2.3^{\circ }$.

We consider three specific scenarios, with two moving targets, to analyse how the sparse learning and the track formation stages work. These scenarios involves $N_R=15$ range bins located in the far field of the radar with spatial resolution $R_{min}=30$ m. The Doppler interval considered, according to the unambiguous frequency Doppler condition $\mid f_d\mid =$ PRF$/2$, is taken from $-458.5$ Hz to $+458.5$ Hz, with $18.34$ Hz separation between Doppler points (thus $N_d=51$). The SINR is defined after pulse integration considering the actual steering vector at the mth scan:

$$\mathrm{SINR}=\mid {\alpha }_{m,r,d}{\mid }^2{\mathit{v}}_{r,d}^{†}{\mathit{M}}_m^{-1}{\mathit{v}}_{r,d}.$$

(28)

In each considered scenario, targets move according to a constant velocity dynamic model, that is positive if target is approaching the radar. Neglecting the edge effects between adjacent range cells, Q accounts for the maximum number of range bins transitions of the target between two successive scans:

$$Q=⌈{\displaystyle \frac{\mid v_l\mid SRT}{R_{min}}}⌉,$$

(29)

where $v_l$ for $l=1,\dots ,L$ is the target velocity. Thus, according to the frequency Doppler definition, the maximum target velocity is $\mid v_{l,max}\mid =49.5$ m/s so that it results $Q_{max}=2$.

Both SLIM and clustering are iterative algorithms, so it is essential to set some convergence criteria. For both algorithms, a convergence threshold is set equal to ${\mathsf{\Delta }}_{SLIM}={\mathsf{\Delta }}_{\mathit{k}-\mathit{means}}={10}^{-3}$, that means the estimation performed can be accepted because it doesn’t vary significantly between two consecutive iterations. Moreover, the SLIM maximum number of iterations is 25.

Since closed-form expression of probability of detection and tracking ($P_{dt}$), defined as the probability to detects the targets and to estimate their tracks correctly in range-Doppler domain, and the desired false track detection probability ($P_{ft}$) are not available, we resort to Monte Carlo counting technique. More precisely, in order to evaluate the threshold necessary to ensure a preassigned value $P_{ft}$, we resort to $100/P_{ft}$ and ${10}^4$ independent trials, respectively. Finally, we compute also the Root Mean Square Error (RMSE) defined as:

$$RMSE=\sqrt{\frac{1}{n}\sum _{m=1}^M{\displaystyle \sum _{(r,d)\in {\hat{\mathsf{\Omega }}}_m^{\left(L\right)}}}(\mathsf{\Delta }{r}_r^2+\mathsf{\Delta }{d}_d^2)},$$

(30)

where n is the number of the Monte Carlo trials, $\mathsf{\Delta }{r}_r$ and $\mathsf{\Delta }{d}_d$ are the range and Doppler error respectively between the true and the estimated range-Doppler target coordinates.

As natural competitors for the proposed TBD algorithm, we consider two compressed-sensing-based TBD algorithms obaitined by exploiting OMP and CoSaMP which are referred to in the following as OMP-based and CoSaMP-based TBD algorithms, respectively. This choice is dictated by the fact OMP and CoSaMP are considered the main benchmark in the context of compressed sensing [54].

For the sake of clarity,

Table 1. Main simulation parameters. In bold italics are indicated values adopted in Monte Carlo simulation for performance analysis as described in Section 4.

Parameter	Value
Number of antenna sensors ($N_s$)	4
Antennas inter-element spacing (h)	$\lambda /2$
Number of pulses ($N_p$)	4
Carrier frequency ($f_c$)	1.388 GHz
Pulse Repetition Frequency (PRF)	917 Hz
Number of range bins ($N_R$)	15, 12
Range spatial resolution ($R_{min}$)	30 m
Doppler interval	$\mid f_d\mid \le$ PRF$/2$
Number of Doppler points ($N_D$)	51, 25
Azimuth view angle (${\theta }_0$)	$\pi /2$
Number of targets (L)	2
Maximum radial target velocity	49.5 m/s
SLIM convergence threshold (${\mathsf{\Delta }}_{SLIM}$)	${10}^{-3}$
Maximum SLIM iterations	25
k-means convergence threshold (${\mathsf{\Delta }}_{\mathit{k}-\mathit{means}}$)	${10}^{-3}$
Number of radar scans (M)	5, 4
Scan Repetition Time (SRT)	0.68 s
Indipendent Monte Carlo trials (n)	1000

lists the simulation parameters used in the three test cases and performance analysis, and a detailed discussion can be found in Section 3 and Section 4 respectively.

3.1. Case Study 1

In the first test case, we consider the two separate targets approaching the radar with different and constant velocities. At the beginning of the radar scanning, the first target is located in the range bin $r_{1_1}=6$ (starting from a relative position of 15 m) with velocity $v_1=+25.76$ m/s while the second target in the range bin $r_{2_1}=12$ (starting from a relative position of 15 m) with velocity $v_2=+41.62$ m/s. Considering the carrier frequency value and according to the velocity-frequency Doppler relationship, the velocities of the targets correspond to the Doppler frequency $f_{d_1}=+238.42$ Hz (Doppler point $d_{1_1}=39$) and $f_{d_2}=+385.14$ Hz (Doppler point $d_{2_1}=47$) for the first and the second target respectively. Figure 4 shows the outputs of interest obtained for a single Monte Carlo trial at a fixed SINR $=9$ dB according to the definition given in (28). In detail, Figure 4a shows the output of the matched filter at the receiver, that is the input of the overall processing chain. In Figure 4b the output of the SLIM algorithm is presented. By this figure, it is evident that the estimated signal vector (i.e., ${\widehat{\mathit{\alpha }}}_{\mathit{m}}$) during the scanning is sparse, in fact most of its value are zeros. However, non-zero values are representative of the range-Doppler target coordinate and its amplitude. This concept is better highlighted in Figure 4c where the Doppler points are stacked over the correspondent range bin. Thanks to this particular representation, is evident in an intuitive way the target paths during the radar scanning, but it is important to remark that this kind of representation is useful for a visual understanding and it is not used in the next stage of data processing, the track formation, because it leads to a loss of informations in the Doppler domain. Figure 4d–f illustrate in different ways the output obtained by the track retrieval stage. In particular, Figure 4d,e show the range bins and the Doppler point estimated for the detected targets during the radar scanning compared to true target positions. In this specific simulation no errors occur in range bins estimation (Figure 4d), while from Figure 4e is possible to identify some errors in Doppler point estimations attributable to the SLIM algorithm. Finally, Figure 4f summarizes the range-Doppler bin estimation or each targets, highlighting the spatial targets’ movements during the radar scanning.

3.2. Case Study 2

In the second test case we consider again the two targets approaching the radar with an intersection between their trajectories. Specifically, at the beginning of the radar scanning, the first target is located in the range bin $r_{1_1}=10$ (starting from a relative position of 25 m) and the second one in the range bin $r_{2_1}=8$ (starting from a relative position of 15 m). First target velocity is $v_1=+47.56$ m/s (corresponding to frequency Doppler $f_{d_1}=+440.16$ Hz and Doppler point $d_{1_1}=50$), while for the second target it is $v_2=+7.92$ m/s (corresponding to the frequency Doppler $f_{d_2}=+73.36$ Hz and Doppler point $d_{2_1}=30$) so that during the M radar scans the first target overtakes the second one. This test case is quite interesting because, being the targets spatially close, it is possible to understand some limitations of the proposed processing chain. As in the previous test case, Figure 5 shows the specific results obtained for a fixed SINR $=9$ dB. In particular, in Figure 5d is shown the correct range bin trajectory for each detected target while from Figure 5e it is possible to notice some errors, due to the SLIM algorithm, in determining the correct Doppler position. Regarding range bins track estimation, it is also interesting to note that the first target has a double range bin transition between the second and the third scans. This test case is equivalent to which one where targets are spatially distant but with comparable velocities.

3.3. Case Study 3

For the third test case, we consider the two targets moving in opposite directions that means the first target moves approaching the radar while the second one moves away from it. In detail, at the beginning of the radar scanning, the first target is located in the range bin $r_{1_1}=10$ (starting from a relative position of 15 m) with velocity $v_1=+39.63$ m/s (corresponding to the frequency Doppler $f_{d_1}=+366.80$ Hz and Doppler point $d_{1_1}=46$) while the second target is in the range bin $r_{2_1}=6$ (starting from a relative position of 15 m) and moves with velocity $v_2=-25.76$ m/s (corresponding to the frequency Doppler $f_{d_2}=-238.42$ Hz and Doppler point $d_{2_1}=13$). According this configuration targets trajectories cross during the radar scans. Figure 6 shows the specific results for SINR $=9$ dB. In Figure 6d is shown the correct Range bins retrieval for both targets during radar scans while Figure 6e reveals the Doppler points estimation during scanning. It is possible to note a single error occurring in Doppler point retrieval only for the first target in the fist scan.

4. Discussion

In the previous Section, three different scenarios are presented to show the tracking performance of the proposed SLIM-based TBD method. Here, we focus on the performance of the overall proposed processing chain compared with other TBD strategies based on alternative compressive sensing algorithms. More precisely, we resort to the OMP and CoSaMP [43,44,54] the main benchmarks in the context of the compressed sensing. In what follows, we refer to these competitors as OMP-based TBD and CoSaMP-based TBD algorithm. The performance evaluation of the proposed processing chain is carried out resorting to Monte Carlo simulation with 1000 independent trials for each algorithm.

We take into account the scenario as defined in the first test case and, to avoid long simulation times, some parameter intervals are suitably reduced without losing in generality. Specifically, in this numerical simulations we consider $N_R=12$, $N_D=25$ and $M=4$.

Simulation results are shown in Figure 7. More precisely, in Figure 7a we plot the $P_{dt}$ versus SINR, while in Figure 7b we plot RMSE versus the SINR, for the three algorithm respectively. By inspecting Figure 7a, it can be seen that the proposed SLIM-based TBD algorithm obtains unity probability of track detection at the lowest SINR (15 dB) compared with the OMP-based TBD and CoSaMP-based TBD algorithms which is almost the same (21 dB). This means that TBD strategy based on SLIM guarantees a significant of gain of 6 dB respect to the matching pursuit approaches. This also shows that the proposed SLIM-based TBD algorithm is resilient to the measurement noise. It can be seen from Figure 7b that the proposed algorithm obtains the lowest RMSE compared with the other two algorithms. The proposed TBD strategy can perform perfect detection of targets with measurement noise when the SINR is above 15 dB.

Finally, we describe the software and the hardware to do used for the numerical results and performance analysis. Specifically, we use a workstation equipped with an Intel^® Xenon^® CPU E-5-1620 v2 @ 3.70 GHz and 16 GB RAM. MATLAB^® scripts are implemented and developed to simulate the overall processing chain.

5. Conclusions

In this paper we have proposed a novel approach to multiple targets TBD for air surveillance radar applications based on sparse data processing. In the first part of the paper we have presented the problem formulation, then we have proposed the relative processing chain. This is based on three stages: the sparse learning, the track formation, and the ad hoc detector. To test it, we have analysed initially three operating scenarios which differ for targets configuration in terms of relative positions and velocities. The adopted radar parameters refer to the Lockheed Martin FPS-117, a military radar system currently in use. The performance assessment has highlighted that the proposed detection architecture exhibits satisfactory performances and, more importantly, the SLIM-based detector guarantees significant gain in terms of detection probability and track estimation, compared with the natural competitors OMP-based and CoSaMP-based TBD algorithms. Moreover, simulation results show that the proposed processing chain can track and detect multiple targets successfully.

Regarding future developments, considerations could be, firstly, the extension of the present work to the range-azimuth-Doppler domain. Furthermore, particular configurations could be taken into consideration, such as targets with different scattering properties or that change their trajectory. Improvements in the processing chain mainly concern the clustering algorithm for the track formation. In fact, specific constraints could be adopted in the k-means algorithm application to achieve better clustering capabilities, or, in addition, it could be useful to investigate other different clustering or path finding algorithms, such as Viterbi.