# Detection via simultaneous trajectory estimation and long time integration

###### Abstract

In this work, we consider the detection of manoeuvring small objects with radars. Such objects induce low signal to noise ratio (SNR) reflections in the received signal. We consider both co-located and separated transmitter/receiver pairs, i.e., mono-static and bi-static configurations, respectively, as well as multi-static settings involving both types. We propose a detection approach which is capable of coherently integrating these reflections within a coherent processing interval (CPI) in all these configurations and continuing integration for an arbitrarily long time across consecutive CPIs. We estimate the complex value of the reflection coefficients for integration while simultaneously estimating the object trajectory. Compounded with this is the estimation of the unknown time reference shift of the separated transmitters necessary for coherent processing. Detection is made by using the resulting integration value in a Neyman-Pearson test against a constant false alarm rate threshold. We demonstrate the efficacy of our approach in a simulation example with a very low SNR object which cannot be detected with conventional techniques.

## I Introduction

The detection of manoeuvring and small objects with radars is a challenging task [1] and is a highly desirable capability in surveillance applications [2]. Radars emit modulated pulses towards a surveillance region and collect reflected versions of the transmitted waveforms from objects in this area. Small objects induce low signal-to-noise ratio (SNR) signals at the radar receiver. The decision on object presence is made by testing the hypothesis that the received signal contains reflections against the noise only signal hypothesis after the front-end input is filtered with a system response matching the probing waveform, which is known as the matched filter (MF) [3].

In order to detect low SNR objects, many such pulse returns (i.e., multiple measurements) need to be considered as each reflection is at a level similar to the noise background. The sufficient statistics of multiple pulse returns are found by summing the associated reflection coefficients across them, which is referred to as pulse integration [3, Chp.8]. This process is applied on the sampled outputs of the MF stage. These samples correspond to, in effect, measurements corresponding to resolution bins in an equally divided range space. In conventional processing, beam-forming and Doppler processing with these samples are used to further segment the bearing and Doppler space into resolution bins and find the corresponding measurements. Conventional methods for integration over time such as coherent and non-coherent integration integrate pulse returns in the same range-bearing and Doppler bins across time. When objects manoeuvre, however, these reflections follow a trajectory across these bins, and, these methods fail to collect evidence on object existence for a long time due to not taking into account this trajectory. On the other hand, longer integration time provides higher probability of detection for a given false alarm rate, in principle.

One possible solution to providing long time integration for manoeuvring objects is to design filters with long impulse responses that match multiple pulse returns along a selection of possible trajectories (see, e.g., [4]&[5]). The number of filters required in this approach easily becomes impractically excessive with increasing integration time. An alternative approach is to employ a dynamic programming perspective and use a regular probing pulse MF to integrate its outputs along a trajectory estimated simultaneously which corresponds, in a sense, to on-line adaptive synthesis of long time MFs.

Trajectory estimation using the outputs of a pulse MF is often referred to as track-before-detect (see, for example, [6, 7]). The sample that corresponds to the true object kinematic state (i.e., location and velocity) is a complex value that is a sum of the reflection coefficient and background noise [8]. Most track-before-detect algorithms, on the other hand, use the modulus of the MF within models which describe the statistics of the modulus of the MF output. These models are averaged and hence cannot fully exploit the information captured by the measurements. For example, it is well known that the detection performance of these methods can be improved by also taking into account the phase of the data samples [9], in addition to the modulus.

The best achievable detection performance is obtained by coherent processing [3], in which one needs to estimate the complex reflection coefficient from the complex values of the MF outputs, the latter of which are processed by the aforementioned algorithms. This corresponds to using a non-averaged model in which the reflection coefficient is a random variable that remains the same during what is known as a coherent processing interval (CPI), and, is generated randomly for consecutive CPIs [8]. This is challenging partly because estimation of this quantity with a reasonable accuracy requires more samples than one can collect at the pulse-width sampling rate in a coherent processing interval (CPI) [10]. For example, in [11], coherent processing and integration within a CPI is performed with a very high sampling rate that yields a large number of samples in the pulse interval.

In [12], we demonstrated that this can be remedied using a phased array receiver structure. In particular, we introduced a simultaneous trajectory estimation and long time integration algorithm in which the integrated value is then tested against a constant false alarm rate (CFAR) threshold for declaring the existence or otherwise of an object in a Neyman-Pearson sense. In [13], we extend this approach for separated transmitter/receiver pairs, i.e., bi-static channels, with an unknown time reference shift. We recover the synchronisation term by diverting simultaneous beams towards the tested point of detection and the remote transmitter thereby relaxing the commonly used assumption that the remote transmitters and the local receiver are synchronised (see, e.g., [14, 15]).

In this work, we provide a complete exposition of our long time integration and trajectory estimation approach in mono-static and bi-static configurations as well as the multi-static case. In particular, we consider the system structure in Fig. 1 where there are multiple transmitters using mutually orthogonal waveforms. The receiver is a ULA and has the full knowledge of the transmission characteristics except the time reference shift of the separately located transmitters. The front-end signals at the receive elements are the superposition of noise, signals from direct channels, and, reflections from objects.

We consider a long time likelihood ratio test conditioned on a trajectory in a kinematic state space, reflection coefficients, and, synchronisation terms as unknown parameters. In order to estimate the kinematic quantities, we use a Markov state-space model in which the object state consists of location and velocity variables. The measurement model of this state space model involves the radar ambiguity function parametrised on the aforementioned reflection coefficients. These coefficients are estimated by using an expectation-maximisation algorithm [16] realising a maximum likelihood (ML) approach within Bayesian filtering recursions for state trajectory estimate. We show that this is an empirical Bayesian method [17] for realising the update stage of the filter. When these ML estimates are reasonably accurate, the empirical Bayes update is an accurate approximation to the otherwise intractable filtering update equations. For synchronisation, we employ a digital beam-forming technique to simultaneously divert beams towards both the test points of detection and the locations of the separately located transmitters in order to find the respective time reference shifts in the bi-static channels.

The resulting algorithm enables us to collect the entire evidence of object existence at the receiver by i) performing coherent integration in both mono-static and bi-static channels within a CPI, ii) non-coherently integrating across different (non-coherent) channels, e.g., local mono-static and remote bi-static channels, and, iii) continuing integration for an arbitrarily long interval that contains many CPIs. As a result, this approach enables us to detect manoeuvring and low SNR objects which cannot be detected using other techniques.

This article is organised as follows: Section II gives details of the problem scenario and introduces the mathematical statement of the problem. In Section III, we discuss trajectory estimation with the array measurements and detail the aforementioned empirical Bayes approach. In Section IV, we first introduce an expectation-maximisation (EM) algorithm for the ML estimation of the complex reflection coefficients. Then, we detail the ML estimation of the synchronisation term. We combine these estimators and specify the proposed detection scheme in Section V. The proposed detection algorithm is demonstrated in Section VI in comparison with a clairvoyant detector and a conventional scheme in a scenario with a manoeuvring and very low SNR object. Finally, we conclude in Section VII.

## Ii Problem statement

Let us consider the problem scenario in Fig. 1 with a ULA receiver (depicted by red dots), and, transmitters (depicted by triangles) one of which is co-located with the receiver forming a mono-static pair. The other transmitters are located elsewhere and form bi-static pairs with the receiver.

The receiver is comprised of elements spaced with a distance of which will be specified later in this section. Each element collects reflected versions of the transmitted waveforms emitted by both the co-located and the separately located transmitters thereby forming mono-static and bi-static pairs, respectively.

### Ii-a Spatio-temporal signal model

A detailed model for the signals induced at the receiver array by reflections from an object is as follows: We consider an interval of time in which each transmitter emits consecutive waveforms separated by a time length of after modulating with a common carrier that has an angular frequency of . The th transmitted signal is therefore given by

(1) |

where denotes the real part of its input complex argument and is known as the pulse repetition interval (PRI).

We assume that is an orthogonal set of waveforms of pulse duration and bandwidth , i.e.,

(2) | |||||

for , where is Kronecker’s delta function.

Use of such orthogonal transmit waveforms underlies the vision of multiple-input multiple-output (MIMO) radars [18, 19] a particular configuration of which is, hence, the system considered here. Design of orthogonal sets for MIMO sensing was investigated with various objectives such as maximisation of diversity [20] and waveform identifiability [21]. In this work, we consider a narrowband regime in which frequency division multiplexing can be used to achieve orthogonality in practice.

In order to specify the received signal at the array elements, let us consider the geometry of the problem which is illustrated in Fig. 2 for transmitters. The receiver array measures the superposition of signals from different channels which are depicted by coloured lines. In particular, there are i) a local (mono-static) channel (red line), ii) a remote (bi-static) channel (green line), and, iii) a direct channel from the remote transmitter (green dashed line). The first two are reflection channels propagating the reflected waveforms from the object (black circle) towards the receiver array. These channels can be fully separated given the array data by exploiting the orthogonality of the waveforms over time and the capability of spatial filtering thereby diverting multiple beams towards arbitrary arrival angles, simultaneously. These points will become clear in the sequel.

Let us model the signals in the reflection channels. We assume that the reflectivity of the object remains coherent (i.e., unchanged) during the collection of reflections from the pulses in (1). Such a time interval is known as a coherent processing interval (CPI). Modelling of the direct channel signals is introduced later in Section IV-B.

The kinematic state of the reflector (depicted by a black dot) in the 2D Cartesian plane is given by , where is the location, is the velocity, and denotes vector transpose. The distance of to the receiver is related to pulse time of flights. The overall distance a pulse emitted by the th transmitter at and reaches the receiver at after getting reflected at is given by

(3) | |||||

where and denote the distance from the object to the th transmitter and to the receiver, respectively.

The corresponding time of flights are found as

(4) | |||||

where is the speed of light.

The velocity of the object induces an angular frequency shift on reflections which is known as the Doppler shift. This quantity is given by

(5) | |||||

where and are the angle of arrival (AoA) of the reflections to the receiver and the bearing angle of the object with respect to the th transmitter, respectively. These quantities are given by

(6) |

For narrowband reflections, the signals induced at the array elements are characterised by a spatial steering vector as a function of which is given by [22, Chp.2]

where is the separation between the array elements selected as half of the carrier wavelength, i.e., . Substituting this quantity together with in the equation above leads to

(7) |

The superposition of the reflections after demodulation at the receiver is given using (7) and (5) by

(8) | |||||

where is a complex coefficient modelling the reflectivity in the th channel, and, is the time of flight of a pulse given in (4). Here, is an unknown time shift modelling the time reference difference between the th transmitter and the receiver (i.e., a synchronisation term).

The reflections in the received signal are optimally searched by matched filtering [8], i.e., by convolving the input with inverted versions of the probing waveforms. In our scenario, this corresponds to a bank of filters, (see, e.g. [19, Chp.3]). Owing to the orthogonality (asserted by (2)), the channels in (8) will have been separated at the filter outputs^{1}^{1}1Perfect orthogonality of waveforms might not be achievable in practice, nevertheless, design of waveforms with a fairly small mutual cross-correlation has been a productive research area which is also discussed, for example, in [19, Chp.2].. The output of the th filter is given by

(9) | |||||

where denotes convolution and is the auto-correlation of the th waveform given by

(10) | |||||

This output is sampled with a period that equals to the pulse duration . Let us assume that is an integer multiple of , i.e., where . samples of this discrete time vector sequence is given by

where

(12) | |||||

(13) |

The term will be referred to as the temporal steering vector.

Next, this vector sequence is arranged as a cube by folding the two dimensional data array in (II-A) in lengths of samples. The th layer of the resulting cube corresponds to the samples collected between the th and the following pulse, i.e.,

This processing chain is illustrated in Fig. 3 together with the cube which is also known as the radar data cube [3]. The axes of this cube are array index, slow time and fast time. In the fast time axis, we have samples of the filter output, each of which is associated with a time delay of the reflected signal. These time delays correspond to time of flights which can easily be converted to range values using (4). As a result, array measurements from range bin is a slice along the slow time axis given by

For convenience regarding the notation in the rest of this article, we stack columns of and form a data vector. Before specifying this vector, let us combine the signal model in a single entity as a function of the reflector kinematic state which induces the signals and the range bin which is the measurement index:

(15) | |||||

where denotes the Kronecker product operator. Here, is related to the data vector through the associated time of flight and found by evaluating (3)–(6). The th column measurement vector for the hypotheses that a reflector object exist at and the null hypothesis are hence given by

(16) |

where models the noise background of the th channel and is a complex Gaussian random variable with zero mean and covariance of .

### Ii-B Problem definition

We would like to perform a hypothesis test based on the measurement model in (16). These measurements are complex numbers and we are interested in the evaluation of the sufficient statistics for the two hypothesis. Detection/processing using complex measurements are often referred to as coherent detection/processing and conventionally the input is the same resolution bin over multiple pulse returns [3]. Therefore, in order for this operation to maintain coherence, the target position should not be changed.

In order to extend coherent processing to the case of manoeuvring objects and remote transmitters, we introduce the mathematical statement of the problem as evaluation of a likelihood-ratio i) using complex versions of measurements (as opposed to, for example, using only their moduli) for all reflection channels, and, ii) for a time window of CPIs given an object trajectory where is the object kinematic state at the th CPI. This likelihood ratio will then be tested against a threshold in a Neyman-Pearson sense [23, Chp.3]. The detector we consider hence takes the form

(17) |

where are the data cubes for channel over . Here, and are reflectivity and synchronisation vectors across the channels, respectively, defined by

In order to carry out the test in (17), the trajectory needs to be estimated. This is also referred to as tracking and is the subject of Section III along with estimation of the reflection coefficients . Algorithmic strategies for estimating the synchronisation term are introduced in Section IV. These results are combined in Section V and threshold selection is detailed in order to evaluate the detection test in (17).

### Ii-C Sufficient statistics for the likelihood ratio

The likelihood ratio on the left hand side of (17) factorises over as the noise samples for different CPIs are also independent. Each time term also factorises over channel likelihood ratios as the related parameters are independent, i.e.,

(18) |

These measurements also satisfy a locality property in that the number of range bins which are associated with is limited by the support of in (10) which is of duration . Let us define the (range) extend of an object as

(19) |

where

(20) |

with denoting the nearest integer function, and gives the time of flight in the th channel associated with the object state . This range bin has the highest signal-to-noise ratio (in the th channel) given that as a time auto-correlation function typically vanishes towards tails.

As a result, the likelihood ratio in (18) further decomposes into factors over range bins as

(21) |

The numerator terms in (21) can easily be found using the distribution of the noise in the signal model in (16) as

(22) | |||||

The denominator in (21) regarding the noise only hypothesis is nothing but the noise density evaluated at . Therefore, the instantaneous likelihood ratio in (21) after substituting from (22) and the noise distribution is found as

(23) | |||||

where is the Hermitian of its argument, takes the real part of its complex argument, denotes conjugate and denotes modulus of a complex variable, respectively.

The likelihood ratio evaluation given in (23) is advantageous in that only a linear operation needs to be performed on the measurements which is in the form of a whitening transform with the inverse noise covariance followed by an inner product with the signal model. Because the signal model involves the spatial steering vector in (7), this inner product effectively performs beam-forming on the measurements filtering out contributions of other objects at the same range. Note that a second filtering is with respect to the Doppler as the temporal steering vector (13) is also in the signal model.

## Iii Simultaneous tracking and reflection coefficient estimation

In this section, we consider estimation of the object trajectory using coherent pulse returns during a CPI. Object trajectories are modelled as random vector sequences generated by a Markov state space model [24], i.e.,

(24) |

where the Markov transition density is selected as

(25) |

where is the time interval between two consecutive pulse train transmissions (or, the illumination period), models constant velocity motion, and is the covariance matrix specifying the level of the process noise modelling unknown manoeuvres.

The initial distribution is selected as a uniform distribution over the range-bearing interval for the detection test. These intervals often correspond to radar specific resolution bins. Let us denote the corresponding bounded set in the state space by , and a uniform distribution on by . Therefore,

(26) |

Sequential estimation of as data cubes arrive is performed by using Bayesian recursive filtering [24]. Suppose we have the given distribution of the state variable at the time step based on all the measurements collected up to and including CPI , i.e., . In order to update this prior information with the measurement at the th CPI, first, the Chapman-Kolmogorov equation is realised and a prediction density is found as

(27) |

where the first term inside the integral is the Markov transition given by (25).

The update stage of the filtering is given by multiplying this prediction and the measurement likelihood together with marginalising out all other variables, i.e.,

(28) | |||||

where and are prior densities for the reflection coefficient and the synchronisation term, respectively.

The measurement likelihood in (28) is the product of the numerator terms in the likelihood ratio in (21) over the object’s range bins and channels for the time step , i.e.,

(29) | |||||

and is easily computed by evaluating complex Gaussian densities as discussed in Section II-C.

The marginalisation of the reflection coefficients and synchronisation terms, however, is not straightforward: First, one needs to select prior densities for these terms. One reasonable approach is to use a non-informative prior such as Jeffrey’s prior [25, Chp.5]. These priors are useful when they lead to tractable computations in (28) (see, e.g., [26]). In our problem setting, however, Jeffrey’s priors for the reflection coefficients and the synchronisation terms are constant, and, do not help in finding a tractable form for the full Bayesian update in (28).

In order to tackle this challenge, we use an empirical Bayes approach [17]. These methods approximate the integration in (28) by solving an optimisation problem for finding the likely values of the unknowns and evaluating the integrand at those values. In other words, (28) is rewritten as

(30) |

Here, the reflection coefficients and the synchronisation terms act as model parameters to be selected and the second term inside the integration is similar to a prior for them. Because this prior is conditioned on the measurements, more probability mass should be concentrating around the maximum likelihood (ML) estimates of these values. Let us select this density as

(31) |

where denotes assignment and is Dirac’s delta distribution. In other words, we select the model densities given the measurements as a Dirac’s delta distribution concentrated around their ML estimates and , respectively.

After substituting from the empirical priors in (31) into (30), one obtains the empirical Bayes update as

(32) |

where denotes approximate proportionality. The approximation accuracy is better when these ML estimates are obtained using informative likelihoods (as quantified by their Fisher information) and equivalently have small CRLBs.

We will detail ML estimation of the reflection coefficients and the synchronisation terms in Section IV. For the remaining part of this section, let us assume that these estimates are given.

For realising the recursive filtering, a sequential Monte Carlo (SMC) approach known as the particle filter is used [27]. In particular, we use a bootstrap filtering approach for estimating the object trajectory.

The prediction stage at the time step is realised by forming a regular grid of points over representing samples generated from the initial state distribution in (26). These points constitute an equally weighted set of particles. For , we will have found weighted samples, or, particles, representing the state posterior in the previous step. Let us denote this set by

where is the weight of the th sample. The prediction stage is then realised by sampling from the Markov transition as

(33) |

The weights of these samples in the particle set is given by

(34) |

in order for this set to represent the prediction density in (27).

In the update stage, the same sample set is used to represent the state posterior in (32), i.e.,

(35) |

where denotes assignment.

The weights of these samples need to be adjusted using the measurement likelihood (as per the importance sampling principle [28]), i.e.,

(36) | |||||

After finding the normalised weights in (36), we test degeneracy of the weighted particles. The degeneracy test is performed by finding the number of effective particles using

(37) |

and, comparing it with a threshold . When , we perform re-sampling (see, e.g., [27]) and continue filtering with a new, equally weighted sample set

where is output by the re-sampler.

Using the above particle filter, the object state at the th CPI is estimated by using the empirical weighted average

(38) |

where denotes the estimated object state .

A remarkable feature of the processing scheme driven by the Bayesian recursions above is that no fixed selection of the spatio-temporal steering vectors are used. The evaluation of the likelihood in the update stage in (36) specifies the steering vectors through (22) and (15) as a function of the state value . Because are generated by sequential processing of the data cubes over CPIs, the resulting set of spatio-temporal steering vectors adapt to the measurements. This is in stark contrast with conventional processing chains in which the bearing and Doppler space is sampled with equal size steps leading to a fixed set of steering vectors and corresponding resolution bins. Thus, a super-resolution effect is achieved when finding the object locations as demonstrated in Section VI.

## Iv Maximum Likelihood estimation of unknown parameters

In this section, we first introduce the ML estimator for the reflection coefficients. This estimator is an iterative algorithm realising Expectation Maximisation at each step of the recursive filtering detailed in Section III– in particular when evaluating the tracking update in (36)–and, is also central to long time integration detailed later in Section V. Second, we derive the ML synchronisation term estimator used together with the reflectivity estimator in the filter update in Section III.

### Iv-a ML estimation of the reflection coefficients

The reflection coefficients associated with an object at state are unknown constants for the duration of a CPI and vary across consecutive CPIs due to the change of the object position, orientation etc. The likelihood of these reflectivities is found by multiplying the likelihood in (29) with the priors for the other parameters and marginalising them out. Let us use the empirical prior for the synchronisation term (see, e.g., (31)) obtained using the ML estimator detailed in the next section. The likelihood to be maximised for estimation is hence found as

(39) |

It is not straightforward to optimise this function due to the marginalisation involved. In ML problems involving such latent variables as the state variable , the expectation maximisation (EM) method offers an iterative and gradient-free solution [16]. In this method, starting from an initial parameter configuration , an expectation that replaces the original likelihood is maximised. For the problem at hand, these iterations are given for by

(40) | |||||

Let us focus on the computation of the expectation in (40) and its maximisation. The state density function underlying the expectation is a state posterior conditioned on the previously found value of the reflectivities, i.e.,

(41) |

where the density function on the right hand side is nothing but the predictive density of the Bayesian filtering recursions given in Section III. Thus, the samples generated in the prediction stage in (33) and (34) lead to an importance sampling [28] estimate of the expectation. Given , this importance sampling estimate is given by

(42) | |||||

(43) |

where denotes the estimate of the term proportional to in (40).

This approximation is a sum of terms quadratic in . This can easily be seen by substituting from (29) and (22) to (42). The resulting expression is given in (44) (see the top of the next page). After taking the first order partial derivative of this expression with respect to and setting it to zero, the the ML estimate of the th reflection channel is found in closed form given in (45) (see the top of this page).

Note that the ML estimator in (45) takes the inner product of the “whitened” measurements with the signal model given in (15) for each state particle . This operation effectively performs digital beam-forming towards the particle state in space, and, matches its approach speed through its Doppler frequency encoded in . As a results, the estimator will not be interferred by other objects unless they appear very close to the state value in terms of the achievable spatial and Doppler resolution.

After finding for reflection coefficients using (45), convergence is tested by comparing the norm of the difference between the parameter configurations found in consecutive time steps with a threshold, i.e., iterations are terminated at if

where denotes the complex Euclidean norm. A pseudo-code of these steps are given in Algorithm 1.

(44) | |||||

(45) |

### Iv-B Synchronisation of the local processor with remote transmitters

In this section, we detail the ML estimation of the unknown synchronisation term parametrising the time origin shift between the local receiver and the th separately located transmitter. Our approach exploits the fact the data cube for the th bi-static channel contains direct path signals from the transmitter that can be recovered by diverting a digital beam towards the transmitters spatial state simultaneously with other processing tasks on the data cube, e.g., those related to trajectory estimation and reflectivity estimation for other spatio-temporal points.

The direct path signal in the th channel can easily be modelled using the spatial and temporal steering vectors defined in Section II in (7) and (13), respectively. The state of the th transmitter is given by which is associated with the time-of-flight given in (4) as the time to receiver. The angle of arrival is denoted by which is computed using (6). Different from a reflection channel, the unknown reflectivity is replaced with a known pulse energy term. Thus, the CPI measurement vector at the th range bin obtained by sampling the th matched filter output is given by

(46) | |||||

where is the pulse energy, is the noise free signal model associated with the transmitter state , and, is an all ones vector^{2}^{2}2Note that differs from in (15) in that the latter uses the bi-static time-of-flight in both the temporal steering vector and the waveform auto-correlation delay, whereas, the former uses direct path time-of-flight. Because the transmitters are of zero Doppler frequency, the temporal steering vector reduces to an all ones vector scaled with ..

In the presence of reflectors, we will have received a superposition of this signal and reflections from different spatio-temporal states. In order to recover the direct path signal, a spatio-temporal steering vector that matches in (46) is used which is given by

Note that this filter is nothing but a (scaled) beam-forming vector diverting a beam towards and maps the measurement vector to a single complex value given by

(47) | |||||

Here, the noise term is the inner product of the beam-forming vector and the complex Gaussian measurement noise in (16), i.e.,

which itself is a random variable with a complex Gaussian distribution of mean zero and variance .

As a result, the likelihood to be maximised is

(48) |

where the expected value of the complex Gaussian distributions as a function of is given by