Tracking Multiple Marine Ships via Multiple Sensors with Unknown Backgrounds

In multitarget tracking, knowledge of the backgrounds plays a crucial role in the accuracy of the tracker. Clutter and detection probability are the two essential background parameters which are usually assumed to be known constants although they are, in fact, unknown and time varying. Incorrect values of these parameters lead to a degraded or biased performance of the tracking algorithms. This paper proposes a method for online tracking multiple targets using multiple sensors which jointly adapts to the unknown clutter rate and the probability of detection. An effective filter is developed from parallel estimation of these parameters and then feeding them into the state-of-the-art generalized labeled multi-Bernoulli filter. Provided that the fluctuation of these unknown backgrounds is slowly-varying in comparison to the rate of measurement-update data, the validity of the proposed method is demonstrated via numerical study using multistatic Doppler data.


Introduction
In a multitarget scenario, the targets set cardinality and their dynamic states randomly vary with time. The objective of tracking multiple targets is to estimate the number of targets and their trajectories using the data collected from sensor(s) in a joint manner [1][2][3][4]. Currently, there are three major paradigms for this field of study, namely Joint Probability Data Association (JPDA) [1], Multiple Hypotheses Tracking (MHT) [2] and Random Finite Set (RFS) [3,4]. While the first two formers involve modifying single target tracking filters to accommodate the problem of multitarget tracking, the latter applies estimation theory focusing on Bayesian optimality and provide a top-down formulation for solving the multitarget estimation problem [3,4].
Using RFS leads to the development of a series of multitarget estimation algorithms. Several RFS-based filters has been proposed in both the literature and practical applications, such as the Probability Hypothesis Density (PHD) [5], Cardinalized PHD (CPHD) [6,7], and the multi-Bernoulli filters [8]. While these filters and their extensions can give good estimates of the current target states, they do not produce target trajectories without using heuristics [9,10]. A theoretically rigorous and systematic consideration of the multitarget trajectory estimation based on RFS approach was proposed in [11]. This work also derives an exact closed-form solution to the multitarget tracking problem, known as Generalized Labeled multi-Bernoulli (GLMB) filter. This filter can estimate not only the number of the targets but also their trajectories, simultaneously [12]. It has been applied to several problems as tracking with merged measurements [13], track-before-detect [14,15], extended targets [16], cell biology [17,18], sensor scheduling [19], spawning targets [20], distributed data fusion [21], field robotics [22,23] and computer vision [24]. The GLMB filter for multitarget tracking with two sensors has been developed in [25,26]. An efficient implementation of the GLMB filter based on Gibbs sampling whose complexity depends linearly on the total number of measurements and quadratically on the number of hypothesized targets has been presented in [27]. This method has been extended to the multi-scan GLMB filter [28] and the multi-sensor GLMB filter [9].
In the multitarget tracking problem, clutter and detection profile are notable sources of uncertainty [29]. Clutter is the set of false measurements that do not originate from any true target and detection profile models the ability of the sensor to detect targets. Knowledge of these parameters are essential in Bayesian multitarget estimation. Mismatches in parameters of clutter and detection models lead to poor performance of filtering algorithms. While these parameters are unknown and randomly time-varying, they are normally assumed to be known in advance. This assumption is unrealistic in most practical applications and these parameters need to be estimated from training data or manually tuned [29].
Since the adaptability of the tracker to these unknown parameters are important in practice, several RFS filters have been proposed in the literature to perform multitarget tracking with mismatches in clutter and detection probability. Some of the proposed methods that accommodate the unknown clutter rate are given in [30][31][32][33]. A filter which bootstraps the clutter estimator of [29] into the CPHD filter [6] has been proposed in [34]. Several approaches for dealing with unknown detection probability have been presented in the literature, such as [29,35,36]. However, none of these filters can output target tracks. While the GLMB filter can output tracks, and has been applied to several problems without prior knowledge of clutter rate, as in [37][38][39], it is still computationally expensive. A low computational cost bootstrapping method using GLMB filter has been given in [40] for multisensor multitarget tracking with unknown detection probability.
Multisensor multitarget tracking with jointly unknown clutter rate and detection profile is far more complicated than those with a single unknown parameter. The use of multiple sensors leads to multidimensional ranked assignment problem which is the main hurdle in the implementation of the GLMB filter [9]. Furthermore, exploiting background information from training data for the multitarget estimation at each time frame is insufficient due to the time-varying nature of the two mentioned unknown parameters.
This work is aimed to contribute an efficient method for multitarget tracking that not only produces target trajectories but also estimates the jointly unknown clutter rate and detection profile online with low computational cost. By using a simple combination of the two well-known filters, the CPHD and GLMB filters, this method is not only fast in estimating the unknown parameters but also producing trajectories of the targets. Specifically, these two mentioned unknown parameters would be estimated separately by using the λ−CPHD and p D −CPHD filters before feeding to the GLMB filter for the purpose of tracking trajectories. The preliminary results of this work are reported in [40]. Particularly, in [40], the unknown detection probability is treated by the p D −CPHD filter before boostraped into the GLMB filter with known clutter rate. The soundness and effectiveness of the proposed solution are demonstrated in Section 4 via a multiple marine ships tracking application.
The remainder of this work is presented as follows. The backgrounds on GLMB filtering will be given in Section 2. The proposed bootstrapping method will be introduced in Section 3 followed by numerical studies in Section 4. Some concluding marks in Section 5.

Background
Some fundamentals on multitarget state-space model, the CPHD filter, and GLMB filter will be summarized in this section. Following the convention in [11], single target states are denoted with lower-case letters (i.e., x) while upper-case letters denote multitarget states (i.e., X). The corresponding spaces are denoted by blackboard bold letters (X, L, Z, etc). The sequence of variable X i , X i+1 , ..., X j is abbreviated by X i:j . In this work the inner product f (x)g(x)dx is rewritten as f , g . Given a set S, the finite subsets of S is written as F (S), and 1 S (·) denotes the indicator function of S. For a finite set X, |X| represents its the number of elements, and the product ∏ x∈X f (x) for some real-valued function f is denoted by the multitarget exponential f X , with f ∅ = 1. Further, the generalized Kronecker-delta function δ Y whose arguments can be arbitrary sets, vectors, integers, etc., is defined as follows

Multitarget States
As mentioned in Section 1, algorithms using non-labeled RFS cannot produce trajectories without using heuristic techniques [10]. The Labeled RFS framework, introduced in [11,41], is a principled approach to produce target tracks. Moreover, it is the only method that can produce trajectories from the filtering density [10]. In the labeled RFS frame work, a labeled target at time k is represented by a kinematic target state vector x k in state space X and its unique label k in the (discrete) label space L, and hence x = (x, ) ∈ X × L. This unique label is characterized by two parameters: time of target birth τ and the index of individual targets born at the same time ρ, i.e., k = (τ, ρ) ∈ L [11]. Hence, formally, a trajectory of each target is a sequence of consecutive labeled states with the same label [11]. Note that the label space for all targets born up to time k is the disjoint union L k = L k−1 B k where B k is the label space for targets born at time k, and L k−1 is the label space of the targets born prior to time k. To distinguish the unlabeled states from labeled ones, the normal and bold letters (e.g., x, X, x, X) are used, respectively. Suppose that at time k, there are N targets with corresponding states x k,1 , . . . , x k,N , then the multitarget state can be represented as follows: Definition 1. [11] Let L : X × L → L be the projection L(x; ) = , and hence L(X) = {L (x) : x ∈ X} is the set of labels of X. A labeled RFS with space X and (discrete) label space L is an RFS on X × L such that each realization X has distinct labels, i.e., |L(X)| = |X|.

Standard Multitarget Dynamic Model
Given a multitarget state X k at time k, each state (x k , k ) ∈ X k can either exist with probability P S,k+1|k (x k ) and evolve to a new state x k+1 at next time step k + 1 with probability density Bernoulli RFS of the surviving target with state x from time k to time k + 1 and B k+1 be the labeled multi-Bernoulli RFS of the new-born targets at time k + 1, then the multitarget state X k+1 is the union of the surviving targets and the new-born ones, Following the convention in [9], in this work, the set B k+1 is distributed according to the labeled multi-Bernoulli (LMB) density. Furthermore, for simplicity, the subscript k for the current time is omitted, and the next time step k + 1 is indicated by the subscript + .
Assuming that the appearance, disappearance, and movement of each target are independent of the others, the multitarget transition density (The Mahler's Finite Set Statistics (FISST) notion of density is used in this paper for consistency with the probability density [42]) is [11,41] (5) in which the distribution of new-born targets is given by where r B,+ ( ) is the birth probability of new target with new-born label , and p B,+ (·; ) is the distribution of its kinematic state [11]. The distribution of the survival targets is ).

Standard Multitarget Observation Model
Assuming that there are M sensors, each state (x, ) ∈ X can be either detected by sensor s, s = 1, . . . M with probability of detection P The likelihood function of a multitarget state X for sensor s is given as follows [9], is assigned to at most one target, then each target has a distinct label), Using the assumption that the sensors are conditionally independent (More concisely, the sensors do not interfere or influence each other while taking measurements or detections. The measurement noise, missed detections and clutter from each sensor in a multitarget scenario are, therefore, independent from the others), and let us define the following abbreviations then, the multi-sensor likelihood is written as Obviously, the form of the multi-sensor likelihood g (Z|X) in (17) and that of it its single-sensor counterpart in (9) are identical.

Multitarget Bayesian Recursion
Let π k−1 (·|Z 1:k−1 ) denotes the multitarget density of the multitarget state at time k − 1, where Z 1:k−1 = (Z 1 , . . . , Z k−1 ) is the set of all observation history up to time k − 1. For simplicity, we obmit the dependence on past measurements, i.e, we use π k−1 (·|Z k−1 ) instead of π k−1 (·|Z 1:k−1 ). The multitarget Bayes filter use the Chapman-Kolmogorov equation to predict the multitarget state to time k given posterior at time k − 1 as follows [3] where f k|k−1 (X k |X) is defined as the multitarget transition kernel from time k − 1 to time k, and the integral in Equation (18) is the set integral defined for any function f : F (X × L) → R, The multitarget state X k is partially observed at time k, and the RFS Z k is modeled by the multitarget likelihood function g k (Z k |X k ), thus the multitarget posterior at this time is given by Bayes rule:

GLMB Recursion with Bootstrapping Method
In this section, the generalized labeled multi-Bernoulli (GLMB) filter with its recursion is summarized. The proposed method for estimating unknown backgrounds before bootstrapping them into this filter is also introduced.

Definition 2.
A GLMB density is a labeled multitarget density given as follows [11] where the discrete space Ξ is the space of association map histories Θ 0:k Θ 0 × . . . × Θ k , each = (θ 1:k ) ∈ Ξ represents a history of the (multisensor) positive 1-1 map, the weight ω (J, ) and multitarget exponential Noting that, in Equation (21), while ω (J, ) (L (X)) is a function of only the labels of the multitarget state X, whereas p ( ) X depends on entire set X.
The cardinality distribution Pr (|X| = n), existence probability r ( ) and probability density p (x, ) of a track ∈ L are given as follows [11]: p (x, ) = 1

The GLMB Recursion
Since the GLMB filter is an exact closed-form multitarget Bayes filter under the standard multitarget dynamic and observation models [12], and the form of the likelihood function in a single sensor and multisensor cases are identical, the GLMB filter can be implemented via two separate steps (update and prediction) or the combined step (joint-predict-update process). In this work, for the convenience of proposed method, the two step GLMB recursion will be presented.

b. Prediction
Given the posterior multitarget density at current time is a GLMB filtering density with the form of (21), the predicted multitarget density at next time step is calculated under the standard multitarget dynamic model (4) as follows [11]: p (s, )

Adaptive to Unknown Backgrounds
In practice, the both the clutter rate and detection profile are unknown and unpredictably vary with time. Prior knowledge of background models, therefore, are typically unavailable. Mismatch in background models results in degradation of tracker performance [4]. In this section, based on the suite of methods for tackling the unknown clutter rate and detection probability introduced in [4], the bootstrapping method will be proposed.
A technique that accommodates the jointly unknown clutter rate λ and the unknown probability of detection p D has been introduced in [29]. This technique considers clutter as an RFS of "generator targets" or "false targets", and incorporates the non-homogeneous and unknown detection probability into each target state. Each real target state x ∈ X is corresponded to an augmented state x a = (x, a), in which a ∈ X d = [0, 1] is the variable on the probability detecting x. The augmented multitarget state now can be described as follows X a = (x a,1 , . . . , x a,n ) = {(x 1 , a 1 ) , . . . , (x n , a n )} (39) Similarly, the augmented generator target state is x c = [x, a c ] withx ∈ X c be the generator target state, and a c ∈ X d = [0, 1]. The augmented generator multitarget state is a c1 ) , . . . , (x n , a cm )} (40) Then the probability of detection is replaced by a and a c , respectively.
Assuming that the false and true targets are statistically independent, then each of the augmented generator targets can be modeled for their characteristics as appearances, disappearances, and transitions, together with likelihood, detection and missed detection. The multitarget state is then a combination of (augmented) actual targets and clutter generators. Meaning that, the augmented hybrid space X h involving the multitarget state can be defined as follows [29] where " " denotes the disjoint union, and " × " denotes the Cartesian product. The multitarget state (4) and multitarget observation (8)) at time k now become the hybrid ones: with Z a and Z c be the augmented multitarget and augmented generator observations, respectively. The integral of a function f (h) : X (h) → R is given by [29] X (h) Here, the set integral (19) has been applied to both augmented multitarget state and augmented generator multitarget state terms, i.e, [4] f Noting that the the measurement likelihood is kept unchanged a (x, a) g (s) (x) (49) g (x c ) = g (s) While the method proposed in [29] results in good estimates of targets, it do not produce the trajectories of the targets. Moreover, although this method is a closed-form solution of the CPHD recursion with jointly unknown clutter rate and detection profile, it is proposed for single-sensor multiple targets estimation solely. In this paper, we propose a method of using the technique introduced in [29] to estimate the mentioned unknown parameters then bootstrapping them into the GLMB filter for tracking on-the-fly. The structure of the proposed method is given in Figure 1.

Implementation
Since after each filtering iteration, the number of components in the GLMB density grows at an exponential rate, the low weight terms should be truncated for tractability. In this work, we use Murty's ranked assignment algorithm to sample a given number of hypotheses of the multitarget density with the highest probability to be the correct ones. Then these components are propagated through the filtering recursion only. Although the use of Murty's algorithm leads to a cubic complexity in the product of the number of Doppler measurements, its implementation is reasonable because there are maximum 10 targets in this work.

Numerical Study
The advantages of multi-static Doppler radar such as lightweight, wide range of surveillance with high accuracy, and low power consumption lead to its broad applications in both civilian and military applications [43][44][45]. However,the number of the sensors in conjunction with the non-linear nature and and low observability of the Doppler type measurement leads to many numerical difficulties [44,45]. The use of multistatic Doppler-only measurements in a scenario of 10 receivers and one cooperative transmitter has been proposed in [46] and its extended version [47] for joint detection and tracking of one target.
This numerical study based on the model mentioned in [40] with 10 marine ships. Each ship at time k is represented by a 5 − D state vector x k in the surveillance of interest where p k = [µ k , λ k ] T and ν k = μ k ,λ k T denote the position and velocity in the longitude and latitude, respectively; α k is the course of the target, and T denotes the transpose operation. The target dynamic model can be given as follows: where and t is sample period, n k is a Gaussian noise vector of velocity and course noise components with zero-mean. Note that latitudinal and longitudinal measurements are in degrees ( • ), the distance, speed and time are given in nautical miles (M), knots (kn), and hours (h), respectively.

Remark 1. Equation (52)
is resulted from the assumption that the surveillance region is not very far from the Equator.
The new births are assumed to be distributed with labeled multi-Bernoulli RFS distributions of where r (i) B is the i th common existence probability , and p (i)  Table 1 lists out the initial state of ten targets with random time of appearance and disappearance, and the average course isᾱ = 2π/180(rad/s).
The parameters of the dynamic model are given in Table 2. Consider the configuration of multiple Doppler sensors system including two spatially distributed receivers and one cooperative transmitter located as in Figure 2. Based on Doppler effect, this system can measure the speed of a target at a distance by calculating the altered frequency of the returned signals which originate from the emitting pulses of radio signals and being reflected to radar after reaching target [48]. The observation of a target state x k at the s th receiver using Doppler measurement is given by in which p k and ν k have been defined above the Equation (51), p t = [µ t , λ t ] T is the position of the transmitter, and p (s) r ] T is the location of the s th −receiver ; w k is zero-mean Gaussian noise, w k ∼ N (0, Q k ), with covariance Q k = diag [1Hz 2 ] ; and f t is the signal frequency emitted from the transmitter, and c is the light speed.
Since the targets are dynamic in different directions, the value of observation z  Table 3. It can be seen that, not only the state equation but also the measurement one are highly nonlinear. By using the proposed B−GLMB filter, the configuration of multiple marine ships tracking using multiple Doppler radars with ground truths and their tracking results are illustrated in Figure 2. For better visualization of multiple targets, each target is assigned to a distinct color. The results of longitudinal-latitudinal co-ordinate target trajectories are demonstrated in Figure 3.  For evaluating the effectiveness of the proposed method comparing to the fixed-GLMB filter and the Joint-CPHD, 100 Monte -Carlo run has been used, and the distance, location and cardinality errors are calculated via Optimal Sub-Pattern Assignment, OSPA, [49] and shown in Figure 4a. By using this metric, the distance between the set of true multitarget states and that of estimated target states is calculated at each time step. For measuring the error between two set of tracks, the use of OSPA is insufficient, and OSPA (2) is needed. The OSPA (2) errors [50] of the B-GLMB and fixed-GLMB filters are compared and plotted against time in Figure 4b, respectively. For the B−GLMB filter and joint-CPHD filter, the clutter rate fluctuates in the range of λ c = [28,70], and the detection probability changes from 0.75 to 0.98, i.e., p D = [0.75, 0.98]. The fixed-GLMB filter is used with fixed p D of 0.75 and 0.98 and fixed λ of 28 and 70, respectively. The window length used in OSPA (2) to obtain the differences between the true and estimate sets of trajectories in Figure 4b is set at w l = 10. Both the OSPA and OSPA (2) are used with cut-off parameter c 0 = 0 and p = 1.
Obviously from Figure 4a, the errors in distance and location between the set of true targets and the estimated ones using B−GLMB filter is the smallest values comparing to those of the fixed-GLMB and joint-CPHD filters. In addition, the errors in cardinality statistics using fixed-GLMB and B−GLMB are almost identical and better than error measured by joint-CPHD filter. The results of measuring errors between the set of true target tracks and that of the estimated tracks are given in Figure 4b. Once again, the effectiveness of the proposed method in reducing the errors in distances and locations of the target tracks is validated. The cardinality statistics for the B−GLMB filter, fixed-GLMB filter and joint-CPHD over 100 Monte Carlo run are given in Figure 5.

Conclusions
This paper presented an efficient solution to the problem of tracking an unknown and time-varying number of marine ships from multiple sensors with unknown clutter rate and probability of detection. Particularly, these two unknown parameters are parallel estimated based on the λ−CPHD and the p D −CPHD filters, then bootstrapped into the cutting-edge GLMB filter. By using the bootstrapping method, the proposed filter utilizes the advantages of the two former estimators in accommodating the unknown backgrounds and reduces the computational cost from tracking algorithm of the latter filter. The effectiveness and correctness of the proposed method are demonstrated in Section 4. From our best knowledge, this is the first principled online algorithm for tracking marine ships via multiple sensors with unknown backgrounds in Doppler measurements. For future work, one of the forcuses would be investigating the combination of multi-scan GLMB [28] and multisensor GLMB [9] filters for multisensor multitarget tracking.