A comparison of detection performance for several Track-Before-Detect algorithms

A typical sensor data processing sequence uses a detection algorithm prior to tracking to extract point-measurements from the observed sensor data. Track-before-detect (TkBD) is a paradigm which combines target detection and estimation by removing the detection algorithm and supplying the sensor data directly to the tracker. Various different approaches exist for tackling the TkBD problem. This paper compares the ability of several different approaches to detect low amplitude targets. The following algorithms are considered in this comparison: Bayesian estimation over a discrete grid, Dynamic Programming, Particle Filtering methods, and the Histogram Probabilistic Multi-Hypothesis Tracker. Algorithms are compared on the basis of detection performance and computation resource requirements.


I. INTRODUCTION
Traditional tracking algorithms are designed assuming that the sensor provides a set of point measurements at each scan. The tracking algorithm links measurements across time and estimates parameters of interest. However, a practical sensor may provide a data image, where each pixel corresponds to the received power in a particular spatial location (e.g. range bins and azimuth beams). In this case, the common approach is to apply a threshold to the data and to treat those cells that exceed the threshold as point measurements, perhaps using interpolation methods to improve accuracy. This is acceptable if the signal to noise ratio (SNR) is high. For low SNR targets the threshold must be low to allow sufficient probability of target detection. A low threshold also gives a high rate of false detections which cause the tracker to form false tracks. An alternative approach, referred to as track-before-detect (TkBD), is to supply the tracker with all of the sensor data without applying a threshold. This improves track accuracy and allows the tracker to follow low SNR targets.
The main difficulty in the TkBD problem is that the measurement, which is the whole sensor image, is a highly nonlinear function of the target state. Typically, the target state describes the kinematic evolution of the target, and may also include its amplitude. However, the sensor provides a map of received scatterer power, which may have a relatively high intensity response in the location corresponding to the target. One way to deal with this non-linearity is to discretise the state space. When the state is discrete, then linearity is irrelevant, © Commonwealth of Australia 2008. and estimation techniques such as the hidden Markov model (Baum-Welsh) filter or smoother [12] and the Viterbi algorithm [30] can be applied. Several approaches for TkBD have been developed using this method [1], [6], [19], [25], [29]. The problem with using a discrete state space is that it leads to high computation and memory resource requirements.
An alternative to discretising the state is to use a particle filter to solve the non-linear estimation problem [4], [20], [23]. The particle filter uses Monte Carlo techniques to solve the estimation integrals that are analytically intractable. Particle filtering has been used by a number of authors for TkBD, e.g. [5], [11], [21]. It is a numerical approximation technique that uses randomly placed samples instead of fixed samples as is the case for a discretised state space. Particle Filtering may be able to achieve similar estimation performance for lower cost by using less sampling points than would be required for a discrete grid.
Another alternative approach is the Histogram Probabilistic Multi-Hypothesis Tracker, H-PMHT [26], [28]. A key difference between H-PMHT and the algorithms above is that it uses a parametric representation of the target pdf rather than a numerical one. This reduces the computation load of the algorithm significantly. H-PMHT assumes the superposition of power from the scattering sources and probabilistically associates the received power in each sensor pixel with the target and clutter models. After the association phase it can exploit the estimation algorithms for point measurement tracking. H-PMHT is naturally a multi-target algorithm.
Rather than using the whole sensor image, Maximum Likelihood Probabilistic Data Association (ML-PDA) reduces the threshold to a low level and then applies a grid-based state model for estimation [2], [3], [13], [14]. The association of the high number of measurements is handled using PDA. An alternative version using PMHT for data association has also been used [31]. This algorithm will not be considered in this paper, which instead focuses on algorithms that use the sensor image directly.
In addition to estimating the target state, the TkBD algorithm needs to detect the presence or absence of targets. A simple method for this is to extend the state space to include a null state corresponding to the case that there is no target, e.g. [5], [25], [27]. In this case, a target is detected when any state other than the null state has the highest probability. A closely related concept is to use a separate Markov chain for The evolution of the state is modelled by the linear stochastic process where Vk is a noise process and the transition matrix is given by where qs is the power spectral density of the acceleration noise in the spatial dimensions and qi is the power spectral density of the noise in the rate of change of target return intensity. (6) if the target is present, or

B. Measurement Model
The measurement at each time is a 2-dimensional image consisting of a cells in the x-dimension and (3 cells in the ydimension. An example of the data used for simulation in this paper is shown in Figure 1. The measurement in each pixel of the image at time k, Zki,j) , is assumed to be the magnitude of a windowed complex sinusoid in Gaussian noise, as in [21]. Thus the pixel value will be Ricean distributed, if there is a target present, or Rayleigh distributed if there is no target [24]. The measurement pdf is The noise process is different for the discrete and continuous valued state models.
1) Discrete Valued State: Let X denote the set of all possible states. Assume that the states are uniformly sampled so that Xk = [~xq -,¥-r~yS~t h r' (4) for some integers q,r ,8 and t. The algorithms which use the discrete state do not estimate the intensity, but rely on a marginalised likelihood which is described in section II-C. The process noise, Vk, must also belong to X. To reduce the computation overhead, the probability mass function (pmt) of Vk is chosen so that p(Vk) == 0 for all Vk outside a tight region centred on the origin. The implementation for this paper restricts Vk to a single step in anyone dimension (the pmf is a 81 element matrix).
2) Continuous Valued State: In this case, the noise process is the usual Gaussian random variable with covariance Q given by [

II. PROBLEM DEFINITION
As in [20,Ch. 11], consider a sensor that collects a sequence of two-dimensional images. When present, a target moves in the plane according to a known statistical process. The algorithms use two different kinds of target model: the Bayesian and Viterbi algorithms use a discrete valued state space, whereas the particle and H-PMHT algorithms use a continuous valued state space. The true target state used to generate data for simulation analysis follows the latter model.

A. Target Model
For simplicity of notation, assume a discrete time model, with a fixed sampling period, T. The target state at time k, Xk, consists of position and velocity in the plane, and the intensity of the returned signal, i.e. the presence or absence of a target as originally introduced for PDA in [7]. This approach has been used for the particle filter [21] and a generalised version was applied to PMHT [8].
Although there are numerous algorithms for solving the TkBD problem, there is currently no TkBD benchmark, and existing comparisons between the competing algorithms are limited. The purpose of this paper is to compare a number of existing TkBD algorithms and to investigate their performance in terms of detection capability, estimation error, and required computation resource. [5] compared the RMS error of a particle based TkBD algorithm with a grid based Baum-Welsh algorithm. However, that comparison used a single initial target speed (with almost constant velocity) and a single amplitude. A preliminary comparison of the particle filter and H-PMHT algorithms was presented in [9]. This paper extends that comparison by including a broader set of algorithms, by using a more realistic measurement model ( [9] used Gaussian measurement noise), and by adding diversity in the target behaviour. The paper is a summary of results already published by the authors in [10].
This paper compares the performance of four TkBD algorithms on a radar-like simulation problem as a function of target speed and target amplitude. The target is assumed to be well approximated by a point scatterer, and its contribution to the received sensor image is via a known point-spreadfunction. Although the specific point-spread-function used here is the response of Fourier Transform windows, problems with extended targets (where the sensor resolution is relatively high) could easily be explored by instead using an appropriate target template, such as in [5].
The algorithms compared are the optimal Bayesian estimator for a discrete state space, detailed in [25], a Viterbi algorithm, much like that of [1], the particle filter of Rutten et al [21], and the H-PMHT [26]. The first two algorithms represent Maximum A Posteriori and Maximum Likelihood estimation over a fixed grid, the particle filter is a random sampling numerical approximation, and H-PMHT is a parametric approach.  The target peak SNR quantifies the height of the peak of the target point spread function relative to the noise floor, and represents a measure of how easy it is to detect the target. The point spread function is assumed to be normalised such that the contribution to cell i, j is I k when the target is located exactly on the sample point for the cell. Thus the peak SNR in dB is given by 20 log {Ik/(J'2}.

C. Likelihood Ratio
In many cases it is more convenient to deal with the likelihood ratio of the data, rather than the measurement pdf. For the measurement model described above, the likelihood ratio for cell (i, j) is Since the pixels are assumed to be conditionally independent, the likelihood of the whole image is simply the product over the pixels If a prior distribution is assumed for the target intensity,

p(I k ) , then an intensity independent marginal likelihood is given by
This integral can be approximated by a summation.

A. Bayesian Estimator
The posterior pdf of the target state can be recursively determined using the well known Bayesian relationship p(xkIZk) ex p(zklxk) Jp(xklxk-t}p(Xk-lIZk-t}dxk-l- (11) The Bayesian estimator in this paper is a direct approximation to (11) based on a discretisation of the state space. Choose a uniformly spaced set of states, X (which is not necessarily related to the discrete measurement function). (11) can then be approximated by where K is a normalising constant. The approximation is exact in the limit as X approaches~4. The first term in (12) is the intensity independent marginal likelihood, defined by (10).
The discrete state space is augmented with a null state, 0, to indicate the possibility that there is no target. Denote the probability of target death as Pd, and the probability of target birth as Pi; Then the evolution probability in (12) is given by where Vk == Xk -FXk-l and IX I is the number of discrete states in X . The parameters P b and Pd control the detection performance and can be tuned to optimise detection performance. The selection of the state space, X, is a tradeoff between estimation accuracy, which improves with finer resolution, and computation requirement, which increases with IX I. The process noise pmf also affects estimation accuracy, as well as providing some capacity to handle model mismatch between the assumed target model and the true target motion. The algorithm is initialised with p(xo == 0) == 1 and p(xo) == 0 V Xo i= 0.
Once the pdf of the state has been evaluated, a state estimate can be obtained by selecting the state with the highest probability. In the event that this state is the null state, then the algorithm reports that there is no target. To account for the case where the pdf has a peak that is spread over several grid cells, the implementation used in this paper finds the highest probability non-null state and accumulates the probability in the adjacent cells. If the accumulated probability is higher than the null state probability, then a detection is reported.

B. Dynamic Programming
The Bayesian algorithm above is a MAP estimator. It recursively defines the probability of the target occupying a particular location by the superposition of all of the possible paths to that position. An alternative is to use a Maximum (23)   (22) Likelihood (ML) estimator. Rather than accumulate the probability from alternate paths, an ML estimator selects the single best path. An ML algorithm for discrete states is the Viterbi algorithm, which has been applied to TkBD in [1]. The Viterbi algorithm is a batch processor that finds the most likely sequence of states. One advantage of this is that it always produces a estimate consistent with the dynamic model: if the transition model gives probability zero to the transition from state 1 to state 3, then the Viterbi estimate will never contain a transition from 1 to 3. In contrast, a MAP estimate may contain such a transition. The next algorithm considered is the Viterbi algorithm. The main difference between the algorithm used here and that of Barniv is the extension of the state space to include the possibility that there is no target.
The joint posterior probability of the sequence of states Xo ... Xk is given by The Viterbi algorithm is a recursive scheme for maximising the joint pdf above which has linear complexity in time (unlike the size of the joint state space, which is IXl k + 1 ) . Let Ck(Xk) denote the Viterbi cost metric, which is a normalised measure of the log-likelihood of the most likely sequence leading into state xi; As for the Bayesian algorithm, Xk == 0 denotes the case that there is no target. The previous state in the most likely sequence leading into state Xk is denoted (}k-1(xj ).
The algorithm proceeds as follows. The un-normalised cost of the null state, cZ, is given by which is used to define the normalised cost for all states Ck(Xk) == log l(Ztlxt) -cZ + max {Ck-1(Xk-1) + logp (xklxk-1)}. (16) Xk-l The most likely previous state is given by 3) The estimated state sequence is found by backtracking (18) As for the Bayesian algorithm above, the discrete state grid, Pd and P b are tuning parameters. The span of the state space includes a null state, so the algorithm reports that no target is present if the estimated state in (18) is the null state.

c. Particle Filter
The particle filter used in this paper is based on the algorithm derived in detail in [21]. Like the grid methods above, the state space is augmented with a null state to allow for automatic track initiation. This algorithm uses terminology similar to that used for target detection with the Probabilistic Data Association Filter [7], [18]. That is, a binary existence variable, Ek, is defined such that Ek == 1 --+ Xk =1= 0 and E k == 0 --+ Xk == 0. The algorithm makes a direct approximation of the target state posterior p(Xk IE k == 1, Zk) and the existence probability p(E k IZk). [22] demonstrated that this model is significantly more efficient than extending the state vector with a binary existence state.
The algorithm proceeds by constructing two sets of par- The unnormalised birth particle weights are calculated using the likelihood ratio (9) and proposal density (19) 2) Create a set of N; continuing particles using the system dynamics (2) as the proposal function, with weights -(e)i _ 1 r( I (e)i) (21) w k -N e J..., Zk x k 3) The mixing probabilities are calculated using sums of unnormalised weights Mb +M e Pk == - (24) 5) The particle weights are normalised and similarly for yk~).
The H-PMHT algorithm described above updates existing tracks, but does not provide a means for initiating new tracks or terminating old tracks. A typical two-stage approach based on the method in [28] is used for this function. The tracker maintains two sets of tracks: established tracks, that the tracker is confident correspond to targets, and tentative tracks, that the tracker is not confident in. The established tracks are updated first, and they vet the sensor data before it is presented to the tentative tracks. Similarly, the tentative tracks vet the data before it is passed to a new tentative track initiation stage. Established tracks are terminated when the estimated intensity drops below a threshold of -10dB for two consecutive scans. Tentative tracks are terminated when the estimated intensity drops below OdB for two consecutive scans. Tentative tracks are promoted to established tracks if the estimated intensity is greater than OdB for more than three scans. algorithm which alternates between data association and state estimation. The state and mixing proportion estimates at the pth iteration are denoted x~l and fl{~respectively. The algorithm is summarised as follows: 2) Calculate cell probabilities p(i,j) and p(i,j) The algorithms described so far are general numerical approximation techniques applied to the TkBD problem. The final algorithm is an approach specifically developed for TkBD. H-PMHT is derived by interpreting the sensor image as a histogram of observations of an underlying random process. The received energy in each cell is quantised, and the resulting integer is treated as a count of the number of measurements that fell within that cell. The sum over all of the cells is the total number of measurements taken. The probability mass function for these discrete measurements is modelled as a multinomial distribution where the probability mass for each cell is the superposition of target and noise contributions. The H-PMHT data association process probabilistically assigns each individual quantum to the target and noise models. For each model, the individual quanta and their assignment weights are combined to form a single synthetic measurement and measurement covariance. These are then used by a pointmeasurement based estimator. For the special case of linear Gaussian statistics, the synthetic measurement is formed using a weighted arithmetic mean and a Kalman Filter can be used as the estimator. The quantisation is an artificial process, and is removed by taking the limit as the quantisation step size approaches zero.
The H-PMHT measurement model is slightly different to the model in section II-B. Whereas section II-B explicitly represents the target amplitude, H-PMHT uses a relative power representation. In the H-PMHT model, the mean cell value is given by The two sets of particles can then be combined into one large set t, t -c, 6) Resample from N b + N; down to N; particles. Thus after completing the above steps the particles, {x~Ii == 1 ... N c } , with uniform weights, approximate the posterior target state density at time k and t, is an estimate of the probability of target existence.
The algorithm declares a target detected when the existence probability, i.e, 1 -p(Xk == 0), is above a tunable threshold.
The state estimate is then found by taking the mean of the state vectors for each particle. (34) The tracks vet the sensor data following the method proposed in [28]. This is done by scaling each cell based on the tracker model The result of the vetting process is a sensor image that is suppressed in the location of the existing tracks, but unchanged in other regions. New tentative tracks are formed by finding peaks in the vetted data. When there are peaks within a threshold distance in two consecutive scans, a new tentative track is formed.

E. Algorithm Tuning
Each of the algorithms has a number of different parameters which need to be selected to ensure good performance. It is of interest to characterise how algorithm performance varies with these parameters. However, it is impractical to investigate these characteristics in this paper. For this paper, each algorithm has been tuned to give approximately the same false alarm performance, and effort was made to optimise the algorithm code for speed. A more detailed analysis of algorithm performance as a function of various parameters is currently being undertaken by the authors.

IV. DETECTION PERFORMANCE
The performance of the various algorithms was investigated by simulating a scenario with a single target. Since this study is concerned with detection performance, only straight line targets were considered. Various scenarios were constructed by selecting a particular target speed and intensity. Intensity values of 12 dB, 6 dB and 3 dB were considered with target speeds from 0.25 pixels per frame to 2 pixels per frame. For each scenario, one hundred Monte Carlo trials were performed.
The algorithms which sample the state space on a fixed grid may be affected by the position of the target relative to the grid, that is, whether the target is close to a grid point or mid way between them. In order to average over this potential variation, the initial target position for each Monte Carlo trial was randomly sampled from the range 2.5 ... 3 independently in X and Y. The target heading was also randomly sampled from 0 degrees (East) to 45 degrees (North East).
False track performance was quantified using a single realisation of a scenario with no target present and 2000 scans. The long duration was chosen to test whether the particle filter algorithm suffered from degeneracy.
Six metrics were used to measure performance. Overall detection probability was defined as the fraction of Monte Carlo trials for which the target was detected for at least one frame. Instantaneous detection probability was defined as the total fraction of frames for which the target was detected. RMS position error was averaged over those frames when the target was detected. False track count was the number of false tracks formed in the no-target scenario. False track length was the average number of frames for which these false tracks persist. Computation resource was the total CPU time required to evaluate all of the scenarios. This figure was recorded both in seconds, and as a ratio compared with the fastest algorithm.
The overall detection probability is shown in figure 2, the instantaneous detection probability in figure 3, and the RMS position error in figure 4. In each of the figures, the metric is plotted as a function of scenario number. The horizontal (scenario) axis has two labels: the target speed for the scenario is shown below the axis, and the SNR for the scenario is shown above the axis. Vertical dotted lines delineate the scenarios with a particular SNR value. Table I shows the false track count and the computation resource. For comparison, a Probabilistic Data Association filter with Amplitude Information (PDAF-AI) [15]- [17] was also run on the false track scenario. PDAF-AI uses point measurements and includes amplitude as a non-kinematic measurement feature. The PDAF-AI algorithm was run assuming a known target SNR and a detection threshold to give ninety percent probability of detection for that SNR. The false track performance at the SNR values of interest is shown in table I. The false track performance of the PDAF-AI is clearly unacceptable below 12 dB. For the 3 dB case, the rate of false tracks is lower because the false tracks persist for much longer. The other performance metrics were not considered for PDAF-AI since the false track performance was so poor.
All of the TkBD algorithms were able to easily detect targets at 12 dB, and they all detected every target. The H-PMHT had a small delay in detecting some of the high speed (2 pixels per frame) targets. These targets proved to be difficult for the H-PMHT algorithm because of the track initialisation method that it used. The algorithm formed tentative tracks and then rejected or accepted the tentative tracks based on the estimated power of the track. For this strategy to work, the track formation logic must reliably form tentative tracks close to the true target state. Two consecutive measurements were used to initialise the tentative tracks, and the approach gave poor speed initialisation. This degraded the H-PMHT's ability to detect high speed targets.
For the 6dB targets, in the centre region of figures 2 and 3, the numerical approximation techniques continued to detect almost all of the targets, but the H-PMHT only detected about half of the high speed targets. The Viterbi algorithm gave the best instantaneous detection result because it was allowed to back-track. The H-PMHT also gave high instantaneous detection for the lower speed cases because the tentative track hpmht • particle ----+bayes viterbi I history was included. At 3dB, all of the algorithms showed degraded detection performance and found less than half of the targets. Again, the H-PMHT performed poorly for high speed targets. It may be counterintuitive that the Viterbi algorithm performed generally worse than the other numerical approximations. However, the smoothing it performs does not improve overall detection performance. It improves state estimation when the target is detected and allows the estimator to back-track once the target has been detected, but it does not increase the overall number of targets that are found.
The false track performance of all of the algorithms is comparable, since this was a requirement of the algorithm tuning. Thus, the overall conclusion is that the numerical approximation techniques considered give similar detection performance. If batch processing is acceptable, then backtracking can provide some improvement in the percentage of time that a target is tracked. Note that for batch processing, the Bayesian filter used here could be replaced with a Baum-Welsh such as in [5]. The H-PMHT gave worse performance than the numerical approaches when the target had high speed. However, this may be alleviated with an improved initialisation scheme.
The RMS estimation error curves shown in figure 4 show an obvious trend. The error increases as the target speed increases and as the amplitude is reduced. For all cases, the H-PMHT error is significantly lower than particle filter error, which is in tum better than the grid approximations. The smoothing used in the Viterbi algorithm reduced the RMS error a little over the Bayesian filter. As may be expected, the error from the grid-based algorithms was approximately half the grid size.
The computation resource required by the different algorithms shows a more marked difference than the detection performance. As table I shows, the H-PMHT was more than an order of magnitude faster than the particle filter, and more than two orders of magnitude faster than the grid-based algorithms. Significant effort was spent in optimising all of the algorithms. For the H-PMHT, there was no specific computation bottleneck: both the association and filtering code took similar resource. In contrast, the numerical approximation algorithms were all limited by the likelihood calculations. These incurred the vast majority of the processing cost, even though they used external library code. The H-PMHT was much faster because it does not perform likelihood computations.
In summary, the H-PMHT gave slightly worse detection results than the other algorithms, but at a fraction of the computation cost. Given that H-PMHT is already a multi-target algorithm, whereas the others assume only a single target, the results here suggest that a robust initialisation scheme would make it the algorithm of choice.
The results presented in this paper have not addressed some important issues. In particular, the grid spacing for the Bayesian filter and the Viterbi algorithm was fixed at the sensor resolution. Whether this is the best choice was not thoroughly investigated. However, anecdotal trials demonstrated that doubling the grid resolution gave little improvement in   detection performance and marginal improvement in RMS estimation error for four times the computation cost. Given the already high cost for these algorithms, it was decided not to pursue finer resolution at this stage. The number of particles in the particle filter is also a parameter that may be varied. This was not explored in this comparison since earlier work has demonstrated that reducing the number of particles degrades performance too much [9].

V. CONCLUSION
The detection performance of four alternative Track Before Detect algorithms has been investigated for a range of target SNR values and speeds. For most of the scenarios the difference in detection performance was minor, except for targets with high speed, which were not tracked well by H-PMHT. Both of the grid-based algorithms were very costly in terms of computation resource, whereas the particle filter was significantly faster. The H-PMHT was by far the fastest, requiring two orders of magnitude less computation resource than the grid-based algorithms. The difference in computation cost is largely due to the cost of calculating the likelihood ratio of the data: H-PMHT does not use the likelihood ratio, and the particle filter calculates it over fewer points than the grid algorithms.
The comparison also considered the RMS position estimation error of the algorithms, for which H-PMHT was clearly the best. The two grid-based algorithms had an RMS error approximately double that of H-PMHT, and the particle filter RMS error was part-way between.
The scenarios considered used only a single target, whereas many practical situations required the detection of multiple targets. The H-PMHT is a multi-target algorithm, so is already capable of handling such a problem, but the other algorithms require extension to the multi-target problem. Future work is planned to consider multi-target detection and tracking.