A Bayesian track-before-detect procedure for passive radars

This article presents a Bayesian algorithm for detection and tracking of a target using the track-before-detect framework. This strategy enables to detect weak targets and to circumvent the data association problem originating from the detection stage of classical radar systems. We first establish a Bayesian recursion, which propagates the target state probability density function. Since raw measurements are generally related to the target state through a nonlinear observation function, this recursion does not admit a closed form expression. Therefore, in order to obtain a tractable formulation, we propose a Gaussian mixture approximation. Our targeted application is passive radar, with civilian broadcasters used as illuminators of opportunity. Numerical simulations show the ability of the proposed algorithm to detect and track a target at very low signal-to-noise ratios.


Introduction
Most currently available civilian and military radars use collocated transmit and receive antennas to send an electromagnetic signal and detect the signal reflected back by a potential target [1]. However, it has been known since the 1930s that the antennas used for transmission and reception can also be located at different positions [2]. Such a configuration, known as passive radar, has received considerable attention during the last two decades [2,3]. The main reason for this renewed interest is that the transmitted signal needs neither extra hardware, nor extra power by using commercial FM or TV broadcasters as illuminators of opportunity. Moreover, the detection of targets is covert, since a passive radar does not radiate any pulsed signal.
In conventional detection strategies, a threshold is applied on the raw data at a constant false alarm rate (CFAR) to declare the presence of a potential target [1]. This detection stage generates missed detections and false alarms due to the presence of clutter. The main difficulty with this approach is the fact that it is not known a priori whether a thresholded measurement originates from a target or from clutter. This issue, known as the *Correspondence: frederic.lehmann@it-sudparis.eu 2 Institut Mines-Telecom/Telecom Sudparis/UMR-CNRS 5157, 9 rue Charles Fourier, 91011 Evry Cedex, France Full list of author information is available at the end of the article data association problem, can be solved using the wellknown multiple hypotheses tracker (MHT) [4] or the joint probabilistic data association filter (JPDAF) [5]. However, for low signal-to-noise ratio (SNR) targets, the detection threshold must be lowered to allow a sufficient probability of detection, thus generating an excessive number of false alarms.
An alternative strategy, known as track-before-detect (TBD), uses unthresholded measurements [6]. Therefore, TBD methods are generally more computationally demanding, since all available raw data are processed. However, TBD methods enable the detection of weak targets, since the loss of information due to the detection threshold is removed. The approaches available in the literature rely mainly on batch or recursive processing. Methods based on batch processing [7,8] use dynamic programming on consecutive scans of measurements. These batch methods have essentially two drawbacks. Firstly, the target state-space is discretized, thus introducing quantization errors. Secondly, a detection delay must be tolerated, since a decision is usually taken only after processing the entire batch of consecutive scans. Batch methods not relying on discretized state-spaces include the ML-PDA [9] and Histogram PMHT [10] algorithms. Since the focus of this article is on TBD methods processing raw measurements, we will not consider ML-PDA, which processes thresholded measurements (with a low http://asp.eurasipjournals.com/content/2013/1/45 detection threshold). The Histogram PMHT algorithm is able to update existing tracks but its drawback is that an external track confirmation or termination mechanism is needed. Methods based on recursive processing [11][12][13] use Bayesian filtering on a continuous-valued target state-space. However, since the observation model is a nonlinear function of the target state, the required Bayesian recursion does not admit a closed form. Existing implementations of the Bayesian recursion use particle filtering, which has the drawback to be computationally demanding for high dimensional state-spaces [14]. In this article, we introduce a novel TBD algorithm based on a recursive Bayesian methodology. The proposed structure is inherited from classical radar detection theory, where the delay/Doppler space is divided into regularly spaced intervals. Unlike the computationally intensive particle filtering solution retained in [12], we use a Gaussian mixture approximation [15] with a single Gaussian per delay/Doppler bin to propagate the target state probability density function (pdf ) over time. The resulting algorithm has the following interpretation: the weight (resp. the mean) of a Gaussian represents the a posteriori probability that a target is present in the corresponding delay/Doppler bin (resp. the target state estimate given that a target is present in the corresponding delay/Doppler bin). At first, a Gaussian mixture approach, as initially introduced in [15], may seem impractical since the embedded Kalman filtering requires the inverse of matrices of size the length of the observation vector, which is typically very large in TBD. By fully exploiting the statistical independencies in the received signal, we will show how to design a tractable algorithm requiring the inversion of matrices of very small dimension.
The main technical contributions of this article are as follows: • the development of a passive radar system model, enabling recursive Bayesian TBD filtering to take full advantage of the statistical independencies at the matched filter output • the derivation of a Gaussian mixture implementation suitable for a global surveillance of the state-space, by allocating a Gaussian for each delay/Doppler bin • the introduction of an entropy-based target detection rule.
Throughout the article, bold letters indicate vectors and matrices, while I m denotes the m × m identity matrix and 0 n×m the n × m all-zero matrix. A diagonal matrix, whose diagonal entries are stored in vector a and whose offdiagonal entries are zero, is denoted by diag{a}. N (x; m, P) denotes a Gaussian distribution of the variable x, with mean m and covariance matrix P. sinc(.) denote the sinus cardinal function. The dot product of two vectors This article is organized as follows. First, Section 2 describes a system model for passive radar, suitable for recursive Bayesian TBD. In Section 3, we introduce our Bayesian recursion for TBD target detection and tracking, using a tractable Gaussian mixture implementation. Finally, in Section 4, the performances of the proposed algorithm are assessed through numerical simulations and compared with existing methods.

Signal model
An illuminator of opportunity sends a continuous signal of bandwidth B, whose complex baseband equivalent signal is denoted by s(t). At the surveillance antenna, the contribution of a moving target has the form [2] The time-dependent parameters A, φ and τ denote the amplitude, the phase and the propagation delay, respectively. In particular, if ν(t) denotes the Doppler frequency due to the target motion, the first order derivative of φ(t) is given by 2πν(t). For simplicity, the contribution of clutter and ambient noise is modeled as a zero-mean complex additive white Gaussian noise (AWGN) w(t), with variance σ 2 . Let x e , x r and x(t) denote the position of the emitter, surveillance antenna and target in a 3D cartesian coordinate system. Let v(t) denote the target velocity vector. Let f c be the carrier frequency and c the speed of light, then τ (t) and ν(t) can be expressed as [2] Remark 2.1. The contribution of the direct path and ground clutter in (1) can be neglected, using the methods suggested in [3], namely physical shielding, Doppler processing, high gain antennas, sidelobe cancellation, adaptive beamforming or adaptive filtering.

Matched filtering
We assume that the receiver has a reference channel [2] able to recover s(t) perfectly. Therefore, coherent integration can be performed by cross correlating the received signal with the transmitted signal s(t), shifted in delay. Let T denote the integration time. Assuming that T is sufficiently small, the signal parameters A, φ, and τ in (1) can be considered as constant during each integration window. http://asp.eurasipjournals.com/content/2013/1/45 During the k-th integration window, the output of the matched filter corresponding to a delay shift t is given by Injecting (1) into (2), we obtain where n k (t) is the noise term Using the change of variable u = θ − t − kT, (3) becomes Define the autocorrelation function (AF) as then (5) can be written as The noise term n k (t) is Gaussian distributed and has the following first and second-order statistics Now, sampling the matched filter output at the Nyquist frequency, i.e. at delay shifts of the form where t 0 is the delay associated with the direct path from the emitter to the surveillance antenna, we obtain the vector of noisy observations y k =[ y k (t 0 ), . . . , y k (t I )] T for the k-th integration window. We introduce the notation y 1:k to denote the collection of past and present observation vectors {y 1 , . . . , y k }.
Assumption 2.2. The signal s(t) is a noiselike waveform. Therefore, the autocorrelation function (6) is a assimilated to a thumbtack function [1], i.e. Figure 1 gives an illustration of an autocorrelation function satisfying assumption (2.2).
It follows from (8) and (9), that the elements of y k can be considered as independent Gaussian variables.

State-space representation
According to Section 2.2, the dynamics of a target at the k-th integration window can be represented by a continuous-valued vector where a k + jb k , τ k and ν k denote the target's complex amplitude, propagation delay and Doppler frequency, respectively. Using the dynamical model for the complex amplitude introduced in [16], we obtain Considering that the Doppler frequency is proportional to the first-order derivative of the delay and using a constant velocity model, the dynamics of the target, at the discrete time instant k, are described by Equations (10) and (11) can be written as a discrete-time process equation where the process noise u k ∼ N (0 4×1 , Q) accounts for unmodeled perturbations and is assumed independent of the observation noise.

Observation likelihood
Assuming that observation y k (t m ) originates from the Gaussian distributed background noise, according to (8) its likelihood can be written as Using the independence of the observations, a property obtained as a result of assumption (2.2), the likelihood of the observation vector y k , given that all components originate from the background noise is given by the following factorization Let us now consider an hypothesized target, whose propagation delay lies in the i-th delay bin where i ∈ {1, 2, . . . , I}. Again, using the independence of the observations the likelihood of the observation vector y k conditioned on x k can be factorized as where correspond to the observations affected (resp. unaffected) by the presence of an hypothesized target in the i-th delay bin. In Bayesian filtering, the conditional likelihood needs to be known only up to a proportionality factor (see Section 3). Therefore, dividing (14) by the constant (13), we obtain the more convenient likelihood ratio [11] These factorizations will later prove useful in reducing drastically the complexity of the proposed Bayesian TBD recursion (see Remark 3.2). From (7) and the noise statistics in (8), we have where the observation function associated to i-th delay bin has the form and the noise covariance matrix is R = σ 2 2T I 4 .

Bayesian recursion for TBD multitarget detection and tracking
Recursive Bayesian filtering consists in propagating the a posteriori pdf p(x k−1 |y 1:k−1 ) forward in time, so as to obtain p(x k |y 1:k ), by taking into account the new measurement y k at instant k. It is well-known that this is achieved by applying successively the following steps [17]: (1) Prediction (2) Correction Unfortunately, in our case the integral in (18) and the multiplication in (19) do not admit a closed form due to the nonlinearities in the dynamics (see (10)) and in the observation model (see (16)). Therefore, some form of approximation is needed. For the purpose of target detection, we assume no prior knowledge about the location of a target and not even prior knowledge of its existence. Thus we seek a Bayesian recursion able to perform a global surveillance of the entire state-space. Monte Carlo approaches like particle filtering [11][12][13] are not well suited for this propose. The reason is that the resampling step of particle filtering has a natural tendency to eliminate prematurely entire regions of the state-space (corresponding to low particle weights) [18]. This phenomenon prevents long enough coherent integration for low SNR targets to generate particles with significant weights, unless a prohibitive number of particles is employed. As a remedy, we propose a parametric approach, where the pdfs in (18) and (19) belong to known distribution families.

Choice of a distribution family
A usual choice is the Gaussian distribution family [19], which leads to the simple extended Kalman filter (EKF) [20] for the desired recursion (18) and (19). Obviously, this approach would fail here because inherent approximations due to the linearization of the process equation and http://asp.eurasipjournals.com/content/2013/1/45 observation function [20] are invalid on the entire statespace. Therefore a more careful choice of distribution family is needed.
We propose to partition the state-space in delay/Doppler bins of equal size. Let us consider discrete values of the Doppler frequency variable ν, of the form: where f 0 denotes the lowest Doppler value and ν the discretization step. The discretization of the delay in (9) and Doppler frequency in (20) defines an implicit partition of the delay/Doppler plane into bins, as illustrated by Figure 2. We define the i-th delay bin as The observation function can be locally linearized with respect to the delay variable τ inside each delay bin. Similarly, we set the value of ν so that the process equation can be locally linearized with respect to the Doppler variable ν inside each Doppler bin. ν is thus a parameter of choice depending on the radar application at hand.
We adopt the Gaussian mixture distribution family [15], with a single Gaussian per delay/Doppler bin of the form The mixture weight w i,j k can be interpreted as the probability that a target is present in bin (i, j) at instant k.
represents the target state pdf, given that a target is present in bin (i, j) at instant k.
The reason for this choice is that each component of the Gaussian mixture now verifies locally the linearization approximation of an EKF. Next, we show how the desired recursion (18) and (19) can be expressed in closed form, while preserving the form (21) for each time instant.

Initialization
Assuming no prior knowledge, the probability of target presence must be the same in each bin. Also, given that a target is present in bin (i, j), the target state pdf must account for the initial uncertainty over the entire bin extent. Therefore we choose where where σ 2 a is related to the dynamic range of the target amplitude.

Prediction
Assuming that the a posteriori target state pdf at instant k − 1 belongs to the Gaussian mixture distribution family (21), it can be written as Injecting (26) into (18), the predicted target state pdf becomes and F i,j k is the jacobian matrix of f (.) with respect to the state The demonstration is postponed to Appendix 1.

Correction
Injecting (27) into (19), we obtain where and H i,j k is the jacobian matrix of the observation function h i k (.) with respect to the state The demonstration is postponed to Appendix 2.  (14), has a huge impact on the complexity of the proposed algorithm. Indeed we see from (30), that the correction step in each delay/Doppler bin requires only a 4 × 4 matrix inversion. Instead, a straightforward application of the original Gaussian sum methodology in [15] (i.e. using all the elements of y k for the correction step in each delay/Doppler bin) requires a full-fledged (I + 1) × (I + 1) inversion per delay/Doppler bin, which makes it unusable in practice, even for moderate values of I. In fact, the idea of reducing the size of on-line matrix inversions for a single EKF, using an information filter implementation, has appeared previously in [21]. Here, we use a similar idea in the context of a Gaussian mixture filter using a bank of parallel EKFs.

Per bin mixture reduction
A Gaussian component at instant k − 1 is initially located by design inside the delay/Doppler bin (i, j), i.e. its mean vector However, due to the target dynamics (during the prediction step) or the observations (during the correction step), the updated mean vector x i,j k|k at instant k, is not guaranteed to remain inside the delay/Doppler bin (i, j). Therefore two situations may arise.
In the first situation, delay/Doppler bin (i, j) hosts several Gaussian components (i.e. the target state is now estimated by a Gaussian mixture) including either the Gaussian component originally located in bin (i, j) at instant k − 1 or Gaussian components crossing delay or Doppler bin boundaries between instant k − 1 and instant k. For obvious engineering reasons, we cannot allow the number of Gaussian components to grow exponentially with time. Thus at instant k, all the Gaussian components verifying (31) belong to the delay/Doppler bin (i, j) and are collapsed to a single weighted Gaussian component using moment matching (see [22, p. 210]).
In the second situation, bin (i, j) is empty (i.e. no Gaussian component verifies (31) at instant k). In order to ensure proper surveillance of the entire state-space during subsequent time instants, we assign to bin (i, j) the Gaussian component where is a small weight (fixed to 10 −5 ) and the parameters x i,j 0 , P i,j 0 have been defined in Section 3.2. Finally, the weights of the Gaussian components must be renormalized, so that they sum to one.

Target detection and state estimation
Define Z k as the random variable associated to the determination of the position of a target in the delay/Doppler grid of Figure 2, at instant k. Then, Z k is a discrete random variable taking values in the ensemble of bins {(i, j)}, with 1 ≤ i ≤ I and 1 ≤ j ≤ J. In Section 3.1, the mixture weight w i,j k has been defined as the probability that a target is present in bin (i, j) at instant k. Therefore, the probability mass function of Z k is and the average uncertainty about the location of a target in the delay/Doppler grid at instant k, is the entropy of Z k [23] expressed in bits. We know from information theory, that the average uncertainty H k is maximum when Z k is an equiprobable random variable, which according to (23) happens when k = 0. As more and more observations are processed, H k decreases with k when a target is present. We consider that a target has been located within one delay/Doppler bin, if the average uncertainty is strictly less than 1 bit (which corresponds to an equiprobable choice between two bins). If H k < 1, the delay/Doppler bin containing the detected target, (î,ĵ), is obtained by applying the maximum posterior mode (MPM) criterion. Note that at very low SNR, H k must be first order low-pass filtered before thresholding, in order to eliminate most false alarms. Then the target state estimatex k and covarianceP k is obtained by appling the minimum mean-square error (MMSE) criterion

Complexity evaluation
The complete target detection and state estimation procedure is summarized in Algorithm 1.

end if end for
It is well known that the complexity of one recursion of the EKF is O (N 3 x ) [14], where N x is the dimension of the target kinematic state. Neglecting the contribution of occasional per bin mixture reductions (see Section 3.5) and of the target detection and state estimation stage (see Section 3.6), the computational complexity of the proposed algorithm can be evaluated as O(N 3 x IJ) per scan.

Simulation results
We consider a digital radio broadcaster as illuminator of opportunity, sending a digital audio broadcasting (DAB) signal using transmission mode I [24]. The modulation http://asp. According to [25], the AF of the transmitted signal has the form χ k (t) = sinc(Bt). σ a is fixed to 100. This corresponds to a 40 dB SNR difference between the lowest and highest possible target SNR, typical of radar applications. Moreover, the size of the Doppler bins is set to ν = 12.54 Hz. This value was found by trial and error, by augmenting progressively the size of the Doppler bins, until the linearization of the process equation inside each Doppler bin leads to an unacceptable deterioration of the proposed method at low SNR. Besides, due to the limitations imposed on target velocities, the frequency shifts of interest are in the interval [ −400, 400] Hz, so we set f 0 = −400 and J = 64.
For all simulations, a high-speed constant velocity target, whose parameters are listed in Table 1, is considered.

Benchmark batch TBD algorithm
In order to assess the performances of the proposed method, we seek a benchmark algorithm having similar features in order to provide a fair comparison. Namely, the benchmark algorithm must be a TBD method, performing a global surveillance of the state-space (i.e. of all delay/Doppler bins at each scan) and able to detect automatically the presence/absence of a target in the field of view. The batch processor proposed in [8] is good candidate. Joint tracking and detection is achieved using a generalized likelihood ratio testing strategy (GLRT). In order to obtain a fair comparison, the delay/Doppler space is oversampled in such a way that the average running time per scan is approximately the same as for the proposed method, that is Consequently, the benchmark method reduces to a Viterbi algorithm (VTA), whose cost metric is based on the squared modulus of raw matched filter outputs (the reader is referred to [8] for details). Here the raw matched filter output corresponding to the k-th integration window, associated to delay t i and Doppler shift f j , is expressed as Coherent integration is performed over consecutive scans indexed by k = 1, . . . , M forming a batch, where M is a parameter of choice. In its original version [8], the benchmark method waits until the end of each batch before making a decision and performing the backtracking stage if a target is declared. Here we use a modified detection rule for the benchmark algorithm. At each of the M available scans, the cumulated metric of the best path in the VTA is compared to a threshold, corresponding to a probability of false alarm fixed to 10 −4 . With this modification, a target is declared as soon as one of the M thresholds is exceeded. However, we choose to begin the backtracking stage only at the end of the M scans even when the target is detected before, since revisiting the state history at the end may lead to a better path.
Neglecting the contribution of the backtracking stage, the computational complexity of the benchmark batch TBD algorithm can be evaluated as O(54IJ) per scan.

Performance comparison
The matched filtering integration time T is chosen small enough so that the received signal's phase in (2) (resp. the received signal's Doppler in (36)) remains approximately  constant during one integration window. Therefore the proposed TBD (resp. the benchmark TBD) algorithm uses matched filtering with integration time equal to one (resp. 32) OFDM symbol(s). We assume no prior knowledge about the existence of the target. Also, no prior knowledge about the birth and death instants of a target is available. Therefore, both algorithms are reinitialized every 0.2 s in order to detect a new target appearing in the radar field of view (or drop a disappearing target). Consequently for the benchmark TBD method, the VTA processes a batch of 0.2 s of received signal, that is M = 5 consecutive scans given that a scan becomes available every 32 OFDM symbols. The computation resources are measured as the CPU time per scan, obtained for a Matlab © implementation of both methods on a 3.16 GHz Intel Xeon machine. Note that from the results in Table 2, the computation resources consumed by both algorithms are approximately the same, thus ensuring a fair comparison between both methods. We first compare the proposed and benchmark TBD algorithms in terms of detection performance for the target parameters in Table 1. Detection performance is measured in terms of mean-time-to-detect (MTTD) and probability of detection, P d . The MTTD is the average time delay between the onset of a target and its actual detection. The results in Table 2 show that both algorithms have approximately the same MTTD. Also, Figure 3 illustrates the evolution of the detection probability versus time, beginning at the onset of the target. The proposed algorithm and the benchmark approach have comparable detection probabilities. If the target SNR is further lowered with respect the value in Table 1, the detection probability drops sharply for both methods. Such an SNR threshold phenomenon is typical of TBD radar detection [6].
We now compare both algorithms in terms of estimation accuracy. Let us first consider a single run of the proposed TBD algorithm. Figure 4 depicts the evolution of the entropy H k (see Equation (32)) over time (t). As expected, in the presence of a target (i.e. for each window of duration 0.2 s such that t ∈[ 0.4, 1.2]s), the target is successfully detected since the entropy drops below the detection threshold. Otherwise when the target is absent, the entropy remains bounded away from the detection threshold and no target detection is declared. Figures 5  and 6 show that the normalized bistatic delay (τ B) and Doppler (ν) are estimated with very good precision, but not before the target is actually detected. The benchmark TBD algorithm has the opposite behavior. Figures 7 and 8 show that only rough estimates of the normalized bistatic delay (τ B) and Doppler (ν) are produced. This is due to the inherent quantization of the state-space in delay and Doppler bins (see Equation (35)). However, thanks to the VTA backtracking, an estimate is made available for every scan, even the first one. These results are confirmed by Monte Carlo simulations. Figures 9 and 10 illustrate the root mean square errror (RMSE) of the normalized bistatic delay and Doppler for the proposed method. The dashed vertical lines correspond to the MTTD after the beginning of each batch of 0.2 s. We observe that between the beginning of each batch and the next dashed vertical line, the RMSE can be quite high. This can be explained by the contribution of the runs during the Monte Carlo simulations, for which the target has not yet been detected (i.e. the entropy has not yet dropped below the detection threshold). We observe that the RMSE of the normalized bistatic delay is steadily decreasing with time after the MTTD is reached and converges to a small value, namely 0.05. A similar behavior is observed for the RMSE of the Doppler shift, which converges to 0.3 Hz.
For the benchmark method, Figures 11 and 12 illustrate the RMSE of the normalized bistatic delay and Doppler, respectively. Due to the quantization of the state space, those quantities need not be decreasing with time. For instance, the normalized bistatic delay drifts away from the nearest grid point during the target existence (see Figure 7). Moreover, the RMSE of the normalized bistatic delay (resp. the RMSE of the Doppler) has a maximum value of 0.25 (resp. 2.5 Hz). As expected, these values correspond approximately to 50% of the grid size in (35). Therefore, the proposed TBD algorithm outperforms the benchmark method in terms of estimation accuracy. If the target SNR is increased, the estimation RMSE decreases (resp. remains constant) for the proposed method (resp. for the benchmark method). This feature may be valuable in counter-battery radars or for weapon fire control systems, which need to locate a target as precisely as possible.

Conclusions
We have presented a novel TBD algorithm for weak radar target detection. The proposed method, derived by applying Bayesian filtering on raw matched filter outputs, cannot be obtained analytically due to nonlinearities in the process and observation models. Therefore, an approximation in the form of a Gaussian mixture implementation is introduced, that reduces to a bank of interacting EKFs. The proposed method has two distinctive features. Firstly, comparing to existing Gaussian mixture filters, the exploitation of the independencies at the matched filter output reduces drastically the computational complexity of the EKFs. Secondly, by allocating a Gaussian component to each delay/Doppler bin, a global surveillance of the state-space is ensured for each scan. With focus on a passive radar application using digital audio broadcasters as illuminators of opportunity, we have shown that the proposed approach outperforms classical TBD strategies based on the VTA. Future work will tackle the problem of finding a Markov Chain Monte Carlo implementation of the proposed work with acceptable complexity. Indeed, a naive particle-based implementation is unable to explore the entire state-space without omission, unless a prohibitive number of particles is used. This phenomenon is related to the resampling step of particle filtering, which has a natural tendency to eliminate prematurely entire regions of the state-space before a weak target gets even a chance to emerge with sufficient weight. Finally, an extension to multi-target scenarios will also be considered.
where F i,j k is the jacobian matrix of f (.) with respect to the state It follows that locally around x Finally, we obtain (see for instance [26, p.  × p(y k (t i−1 ), y k (t i )|x k ) p 0 (y k (t i−1 ))p 0 (y k (t i )) .
A local linearization of the observation function (17) around the x i,j k|k−1 leads to Finally, we obtain (see for instance [26, p. 40-41])