FISST Based Method for Multi-Target Tracking in the Image Plane of Optical Sensors

A finite set statistics (FISST)-based method is proposed for multi-target tracking in the image plane of optical sensors. The method involves using signal amplitude information in probability hypothesis density (PHD) filter which is derived from FISST to improve multi-target tracking performance. The amplitude of signals generated by the optical sensor is modeled first, from which the amplitude likelihood ratio between target and clutter is derived. An alternative approach is adopted for the situations where the signal noise ratio (SNR) of target is unknown. Then the PHD recursion equations incorporated with signal information are derived and the Gaussian mixture (GM) implementation of this filter is given. Simulation results demonstrate that the proposed method achieves significantly better performance than the generic PHD filter. Moreover, our method has much lower computational complexity in the scenario with high SNR and dense clutter.


Introduction
Optical sensors have been widely applied both in important military and civil areas due to their properties of long detection range, high concealment ability and large coverage area. Since there is usually a long distance between targets and sensor which generates a low SNR and dense clutter OPEN ACCESS scenario, multi-target tracking in the image plane of optical sensors is a very difficult problem. In multi-target tracking, the aim is to estimate the number of a set of targets and the state of each target from a set of measurements received. However, due to the variation of targets number with time in the field of view of the sensors and the existence of miss detection and dense clutters, multi-target tracking in image plane remains a challenging problem.
In addition to the location measurement, amplitude information has been proven to improve tracking performance [1][2][3][4][5]. One of the pioneering techniques was the approach proposed by Colegrove, Lerro and Bar-Shalom [1,2] where the probability data association (PDA) filter utilizing target amplitude was applied in the context of single-target tracking. The target amplitude has also been incorporated in the multiple hypothesis tracking (MHT) framework [3] and Viterbi data association scheme [4]. More recently, the significance of target amplitude has been explored for data association of closely spaced targets in [5]. Although significant progress has been made recently, the approaches mentioned above are based on data association technique, which requires expensive computational cost in most circumstances. Therefore, the traditional data association approach is not the optimal option for multi-target tracking in image plane of optical sensor in the scenarios with time-varying targets number, low SNR and dense clutters.
A good candidate is the emerging Bayesian approach in the framework of finite set statistics (FISST) proposed by Mahler [6]. FISST provides a set of mathematical tools that allows direct application of Bayesian inferencing to multi-target problems. PHD proposed by Mahler [7] is a computation tractable approximation for the optimal multi-object Bayes filter based on Finite Set Statistics (FISST). Operating on the single-target state space, the PHD filter avoids the combinatorial problem that arises from data association and thus leads to superior performance in comparison with traditional MHT algorithm in multi-target tracking [8,9]. These features render the PHD filter extremely attractive to many researchers. The primary PHD applications can be found in [10][11][12][13]. In [14], Clark et al. incorporated target amplitude into a PHD filter for the first time in multi-target tracking, which improves the tracking performance. However, their amplitude is modeled for radar and sonar application and the computational complexity analysis of their method is absent.
In this paper, we propose a FISST based method which uses signal amplitude information in a PHD filter for multi-target tracking in the image plane of optical sensors. Based on analyzing the imaging characteristics of the optical sensor, we model the signal amplitude and incorporate it into the PHD in the form of amplitude likelihood ratio. In the situation where the SNR of target is unknown, we then present an alternative method based on this amplitude model. Simulation results demonstrate a significant improvement of the proposed method in tracking performance. Furthermore, the computational complexity of our method in the scenarios with different clutter density and SNR is also discussed. This paper is organized as follows: Section 2 models the signal amplitude generated by the optical sensor. Section 3 demonstrates how to incorporate the amplitude information (for known and unknown SNR case) into a PHD filter and then the GM implementation for this filter is given. Section 4 presents the simulation results that validate the proposed method. The conclusions are given in Section 5.

Amplitude Likelihood Ratio
In an optical multi-target tracking system, the original images taken by the sensor are usually processed by background suppression before being further used for target tracking. Assuming the noise is addictive, we define the SNR d of the residual image as [15][16][17]: (1) where s is the mean signal value of target, σ is the standard deviation of the residual image. Each measurement from the image consists of the two-dimensional position vector z in the image plane and the corresponding amplitude a ≥ 0, that is, a measurement vector has the form . For simplicity, we assume that the value of the target signal has no spreading in the image plane. We assume that the noise is Gaussian, then the probability densities of the amplitude of the false alarms and the target p 0 (a) and p 1 (a|d) can be written as [18]: This leads to the probabilities of false alarm and detection and with a detection threshold τ, is the probability error function. According to Equations (4) and (5), we have: In target tracking application we usually choose as a fixed value. Given , the threshold τ can be calculated via the inverse form of Equation (4) with the parameter σ, and the probability of detection can be calculated via Equations (5) or (6) for the target with SNR = d. Table 1 lists the values of for different SNR values d and specified with σ setting to be a normalized value, i.e., σ = 1.   When the target SNR d is known, we can get our amplitude likelihood functions for the false alarm and the target as: The amplitude likelihood ratio given a threshold τ is defined as [19]: (9) We use the notation c a (a) and ρ a (a|d) to denote the case where τ = 0, then we have: From Equations (5) and (9) we can see that the calculations of and rely on a specified known target SNR, however, this requirement cannot be satisfied in most practical tracking systems. We adopt an alternative approach to circumvent this issue next.

Method for Unknown SNR
When the SNR of target is unknown, one straightforward approach would be to estimate the unknown parameter d from the measurement amplitudes a. However, this approach requires a large number of measurements from the target to achieve an accurate estimate of d. Furthermore, due to the unknown association between measurements and targets in multi-target environment with clutter, the approaches of estimating d usually fail. Similar to the idea introduced in [14], we adopt an alternative approach where we do not attempt to estimate d at all. Instead we marginalize out the parameter d over the range of possible values and find a probability of detection for and a likelihood ratio for that is not conditional on d.
Consider that we always have some prior information about the targets being tracked, so we assume that p(d), defined on the possible SNR values [d 1 ,d 2 ], gives the expected probability distribution of SNR values. Since the amplitude distribution in Equations (2) and (3) are symmetric thus have no biases in high or low SNR targets, the reasonable choice for p(d) is the uniform distribution U[d 1 ,d 2 ]. Then we define the probability of detection and amplitude likelihood ratio where SNR is unknown as: From Equations (5) and (12) we have: (14) Note that the over the marginalized region [d 1 , d 2 ] can be computed with numerical integration offline since it does not need to be computed at each iteration. The computation of as in Equation (13) can be simplified and will be presented in Section 3.2.

PHD Filter with Signal Amplitude Information
Suppose that at time k there are N k target states , each taking values in a state space ; and measurements (detections) , each taking values in the observation space . In PHD filter, a multi-object state and a multi-object observation are represented by RFS: (15) where F(X) and F(Z) are the finite subsets of X and Z, respectively. The state of each target contains the position and velocity in the image plane, while the measurement z is defined in Section 2. We assume that each target follows a linear Gaussian dynamical model and the sensor has a linear Gaussian measurement model, i.e., where denotes a Gaussian density with mean m and covariance P, F k-1 is the state transition matrix, Q k-1 is the process noise covariance, H k is the observation matrix, and R k is the observation noise covariance.

The PHD Recursion with Amplitude Information
We abbreviate the PHD filter incorporated with amplitude information as AI-PHD filter. Next we derive the prediction and update equations of AI-PHD filter based on the amplitude likelihood ratio given by Section 2. For simplicity, we do not consider target spawning in this paper.
Step 1. Prediction: The prediction equation of AI-PHD filter is the same as generic PHD filter since their state vector and state transition matrix are the same, i.e., (18) where is the birth term for new targets, is the probability of target survival, is the transition density and is the PHD at time k − 1.
Step 2. Update: The update equation is changed when incorporated with amplitude information.
Analogized to the update equation of generic PHD filter in [7], we have the update equation of our AI-PHD filter as (19) where is the pseudo-likelihood function as where λ and V are the clutter density and area of image plane of optical sensor respectively. Assuming the amplitude a k is independent with target state x k , we can rewrite and as (22)   (23) where L z (x) is the measurement location likelihood function and c(z) is the probability density of the false alarm spatial distribution in the image plane. We assume that the targets are within the surveillance region of sensor, the probability of detection for a given threshold τ is then only dependent on d Substituting Equations (5), (8) and (22)(23)(24) into Equation (20) we have the pseudo-likelihood function of AI-PHD as (25) Equations (18) and (25) compose the recursion of AI-PHD filter. The probability of detection and amplitude likelihood ratio are replaced by and respectively for the unknown SNR case.
We can simplify the computation of by noting the fact that is calculated combined with in Equation (25). From Equations (8) and (9) we have Hence, instead of computing by Equation (13) we can compute the expression directly using the method introduced in Section 2.2, i.e., where is the standard normal distribution function which can be computed easily. Consequently, our approach incorporates the amplitude information into PHD filter with only a minor additional computational load. We show the consistency of our AI-PHD filter with the generic PHD filter. If the SNR of target is set as d = 0, from Equations (7-9) we have . This is the condition under which our AI-PHD filter degenerates to the generic PHD filter.

Gaussian Mixture Implementation
An analytic solution to the PHD filter can be found under linear assumptions on the system and observation equations with Gaussian process and observation noises as described in Equations (16) and (17) [20]. In this case, both the prediction and update equations of PHD are represented by a mixture of Gaussians where the means and covariances are updated with Kalman filter and the weights of the Gaussian components are found using the PHD filter equations. We use Gaussian Mixture implementation of our filter for its simplicity in calculation and convenience of target state extraction comparing to the Sequential Monte Carlo (SMC) method [21].
We assume that the survival probability is state independent, i.e., and the detection probabilities are and for the known SNR and the unknown SNR case respectively. The intensity of the target birth RFS are Gaussian mixture of the form (28) where , are given model parameters that determine the shape of the birth intensity.
We assume a uniform location distribution of clutter in the measurement space, so that the clutter location likelihood is not dependent on the state or the measurement. Hence the clutter location distribution is constant over the measurement space and equals to the reciprocal of the area of the image plane of optical sensor, i.e., c(z) = 1/V. Next we give the prediction and update equations of Gaussian mixture implementation of our AI-PHD filter.

Prediction:
The posterior intensity at time k − 1 is a Gaussian mixture of the form (29) where J k-1 is the number of Gaussian terms with the weights , means and covariances .
Then the prediction intensity is still a Gaussian mixture as where the birth intensity is given by Equation (28) and the means and covariances are computed with the Kalman filter prediction.
; , where is the predicted measurement and is the innovation covariance.
In the PHD update equation, Equation (32), and the weight update Equation (33), we replace the probability of detection and term by in Equation (14) and in Equation (27) when the SNR is unknown.

Simulation
In this Section, by setting up multi-target tracking simulation in the image plane of optical sensor, we examine the performance and computational complexity of our method for known and unknown SNR cases and benchmark them with generic PHD filter with different combinations of probability of false alarm and SNR value.

Simulation Scene and OSPA Metric
Consider a scenario with an unknown and time varying number of targets in clutter in the image region [−300, 300] × [2,000, 2,600] (pixel). Up to N k = 6 targets are generated in this region with the random birth and dieing time instants. Figure 1 shows the true trajectories of each target. All targets in each simulation had the same mean SNR (this is not necessary by the algorithm but simplifies the presentation of results).
Each target has survival probability p S,k = 0.99 and follows the linear Gaussian dynamics in Equation (16) where Δ = 1 s is the sampling period, and σ v = 0.5(pixel/s 2 ) is the standard deviation of process noise. The location measurement follows the observation model in Equation with To mitigate the exponential growth of mixture components, at each time step the number of Gaussian components is capped to a maximum of J max = 100 components, whilst pruning is performed with a weight threshold of T = 10 −5 , and merging is performed with a threshold of U = 2.
We adopt the optimal subpattern assignment (OSPA) metric [22] for the purpose of multi-target performance evaluation since it jointly captures errors in the target state and target number estimates. Given two arbitrary finite sets where n ∏ is the set of permutations on {1,···n}, d (c) (x,y) = min(c,d(x,y)), p is the order that penalizes error of individual element estimates, c is the cut-off parameter that penalizes error of cardinality estimate. We chose p = 2 and c = 30(pixel) in our simulation. Note that the chosen value of cut-off parameter c is significantly larger than the typical measurement noise but significantly smaller than the maximal distance between targets, thus maintaining a balance between the cardinality and localization components of the OSPA error [22].

Filtering Results for Multi-Target Tracking
The effectiveness of our AI-PHD filter for multi-target tracking in image plane of optical sensor is verified through simulation. We assume a moderately cluttered scenario that the probabilities of false    Table 1). For the unknown SNR case, the SNR region is set as [2,10] and the probability of detection is replaced by D p τ which can be computed by Equation (12). Other parameters for the filter are given as in Section 4.1. The true trajectories and filter estimates are shown in x and y coordinates of image plane versus time for AI-PHD filter with known and unknown target SNRs in Figures 2 and 3 respectively (denoted as case1 and case2 accordingly).   From the estimates of AI-PHD filter shown in Figures 2 and 3, we see both for the known and unknown SNR case, the filters can eliminate dense clutter in the scenario with time-varying number of targets. All targets are detected immediately after birth and tracked accurately, which highlights the track initialization, maintenance and termination capabilities of our algorithm. The results also demonstrate superior performance of the algorithm in targets number and target states estimation. We evaluate this performance by Monte Carlo simulation next.

Monte Carlo Results and Analysis
Average OSPA and computation time cost per frame are used to evaluate the performance and computational complexity of our AI-PHD filter. To prove the improvements of the proposed method, we have benchmarked the results against a generic PHD filter which does not use the amplitude information. To make the assessment as fair as possible, the probability of detection D p τ for filters without amplitude information was chosen to be the same as AI-PHD filter of the known SNR case, as in Table 1. 50 Monte Carlo runs were carried out for each combination of FA p τ and d on computer (Intel quad core processor 2.66 GHz, 32-bit operating system, 4 Gbytes RAM) using Matlab (R2009b). Average OSPA for AI-PHD filter of known SNR and unknown SNR case and generic PHD filter are given in Table 2 where the results are divided by '/' accordingly. respectively. This improvement in performance is mainly due to two reasons: firstly, as d increases, the false alarm distribution will poorly represent the target counterpart and there is a big distinction in the target and false alarm distributions; secondly, as FA p τ increases, having more measurements aids the method using the amplitude, since we discard less useful information. Table 2 also shows that the performance of generic PHD filter without amplitude decreases rapidly as FA p τ increases since there are more measurements from false alarms which by no means could be identified from those from targets. In contrast, we see no deterioration in the performance of AI-PHD filter in this case. For the known SNR case especially, the performance increases consistently as FA p τ increases, which means our method works even better in a scenario with dense clutters.The comparison of computational complexity between AI-PHD filter and generic PHD filter without amplitude information is shown in average computation time per frame versus target SNR for different FA p τ in Figure 4. Since we can achieve similar complexity for unknown SNR case with that of known SNR by computing Equation (27) with some fast algorithms, only the result for the known SNR case is given. We see that for different given FA p τ and d, the results from two filters are close with the maximum difference being no more than 1 s. Figure 4(a) shows that in the scenario with low clutter density, AI-PHD filter has only a minor increase in average computation time over the generic PHD counterpart. Furthermore, AI-PHD filter performs an even low value in high clutter density scenario which is shown in Figure 4(b). In the case of = 1 × 10 −3 , d = 8 where this reduction is most obvious, the average computation time of AI-PHD filter is reduced by 53.7% over the generic PHD counterpart, which means the AI-PHD filter has even lower computational complexity than the PHD filter without amplitude information in scenarios with dense clutter and high SNR. The primary reason for this trend is that in these scenarios, the computation time cost is mainly decided by the multi-target state extraction step given the same number of targets and measurements. Incorporated with amplitude information, the update for the AI-PHD filter (see Equation (33)) gives heavier weights to the Gaussian items updated by the measurements from targets, thus updating the PHD with comparatively higher intensity near the real target positions and at the same time, suppressing the intensity of PHD near clutter positions (see Figure 5). Therefore, the updated Gaussian items can be prune and merged quickly and accurately.

Conclusions
In this paper, we have proposed a FISST-based method using signal amplitude information in a PHD filter for multi-target tracking applications in the image plane of optical sensors. We extend the measurement model to include the signal amplitude of observations and then incorporate this information into a PHD recursion in the form of an amplitude likelihood ratio. Based on the assumption that the amplitudes of the measurements from true targets are stronger than those from clutter, we show our method can significantly improve the performance over the one without amplitude information. Furthermore, simulation results also demonstrate that our method has much lower computational complexity in the scenario with high SNR and dense clutter, which makes sense for its practical implementation. Future work will involve improving the tracking performance for the targets with much lower SNR.