Tracking Multiple Targets from Multistatic Doppler Radar with Unknown Probability of Detection

The measurements from multistatic radar systems are typically subjected to complicated data association, noise corruption, missed detection, and false alarms. Moreover, most of the current multistatic Doppler radar-based approaches in multitarget tracking are based on the assumption of known detection probability. This assumption can lead to biased or even complete corruption of estimation results. This paper proposes a method for tracking multiple targets from multistatic Doppler radar with unknown detection probability. A closed form labeled multitarget Bayes filter was used to track unknown and time-varying targets with unknown probability of detection in the presence of clutter, misdetection, and association uncertainty. The efficiency of the proposed algorithm was illustrated via numerical simulation examples.


Introduction
Initiated from the 1930s from a very simple device for aircraft detection, the radar has been developed into complicated systems in both civilian applications and modern warfare for the purposes of prevention as well as interception strategies [1]. As implied from its name, RADAR (RAdio Detection And Ranging), the main purposes of this system are not only detecting targets or obstacles but also estimating several parameters like velocity, range, and bearing of these targets from electromagnetic signals [2]. Since target tracking for radar is subjected to clutter (caused by environment) and various distortions (due to the signal propagation), the estimation accuracy and detection probability are limited. An improvement on tracking performance for radar is the use of a multistatic radar system (MRS). This system is equipped with multiple transmitting and receiving antenna pairs, which are spatially distributed in a large region of surveillance comparing to the antenna sizes [3] to maximize the estimate accuracy and probability of detection.
The large geographical separation of transmitters and receivers is an essential feature of the MRS [4] and leads to a notable increase in potentially useful information. Since different transmitter-receiver pairs can detect targets at different aspects, each target is observed from multi-directional perspectives and its accurate recognition is significantly enhanced [1]. Other benefits of using MRS are the increase of resolution capability, as well as jamming and clutter resistances [5]. In addition, the spatial distribution of radars prevents the whole radar system from physically being destroyed from attacks. multistatic Doppler measurement with bootstrapped unknown detection profile is shown. Third, validation of the proposed solution is illustrated by numerical simulations on marine ships, followed by some concluding remarks.

Multistatic Doppler Measurement Model
Typically, a passive multistatic radar system is realized with either one receiver and multiple transmitters or a single transmitter combined with multiple radar receivers. Although the latter is more costly than the former, its has better observability [40]. An illustration of a multitarget tracking using a multistatic radar scenario with one transmitter and several receivers is shown in Figure 1. The time difference of arrival (TDOA) between the transmitter-ith receiver and transmitter-jth target-ith receiver is proportional to the different range: R Tj + R Rij − L i , where R Tj , R Rij , and L i are the distances of the transmitter-jth target, jth target -ith receiver, and transmitter-ith receiver, respectively. Several solutions for the problems of localization and target tracking via radar signals have been proposed in the literature by means of measuring TDOA or the angle of arrival θ Rij , or the Doppler shift of the received echo [4]. For instance, some methods for single target localization using target range and bearing, range and Doppler measurements and Doppler-only measurements are proposed in Reference [41][42][43], and [44], respectively. In addition, the combination of range and range rate measurements-as well as the combination of range, range rate, and elevation measurements-for multisensor multitarget tracking has been presented in Reference [40]. All the conventional methods using multistatic radar-based target tracking consider the probability of target detection to be known a priori, while it is typically unknown and time-varying [11]. Such an assumption leads to biased estimation or even complete corruption of filtering results since this parameter directly influences observability of the target. Recent proposed methods for tackling unknown detection probability, (see, [45][46][47] for instance), do not produce target tracks. Tracking multiple targets using a multistatic radar system is still an attractive field of study.
Taking advantage of the MRS by using a type of MRS, like the multistatic Doppler radar system, has been used for tracking multiple targets with high accuracy [5,48]. In this work, we considered a multistatic passive Doppler radar system consisting of one cooperative transmitter and two spatially stationary distributed receivers (see Figure 1). The Doppler-only measurement method to track targets based on the RFS approach will also be proposed for further use of the tracking algorithm. By using the Doppler effect, the velocity of the target is calculated from analyzing the pulses of radio signals which are emitted from the transmitter, bouncing off the target then reflecting the receiver [49]. The Doppler measurement of a target with the state x k at the sth receiver is given by: in which the target position p k = [µ k , λ k ] T and velocity ν k = μ k ,λ k T are measured in longitudinal and latitudinal coordinates.
r ] T is the sth-receiver position; T denotes the transpose operation; w k is zero-mean Gaussian noise with covariance Q, w k ∼ N (0; Q k ); and f t and c are the emitted signal frequency of the transmitter and the speed of light, respectively.
Since targets can move in different directions with nonlinear dynamics, and the information collected from Doppler radar is subjected to environmental noise, the measurement in Equation (1) is highly nonlinear.
Several methods have been proposed to tackle those aforementioned challenges using target range and azimuthal direction [41,42], Range-Doppler maps [43], and Doppler-only measurements [44]. However, they are either applied to detect a single target or estimate target positions without producing target tracks in which the probability of target detection is assumed to be known a priori. Hence, a proper method that can produce target tracks with unknown detection probability needs to be considered.

Labeled RFS
The key point in target tracking is using information collected from sensors to jointly estimate the number of targets and their states, as well as the target trajectories. However, both the number of the targets and their states in a multitarget system randomly vary with time, thus it is difficult to follow the target trajectories. Obviously, using vectors to represent multitarget state is insufficient, since the targets in a multitarget state are unordered and can be changed over time. A discussion about vector and finite set representations of multitarget state, given in Reference [16], has shown that a finite set representation is the most appropriate from an estimation viewpoint. The most appropriate model for multitarget state is, therefore, an RFS [10].
Individual targets in a multitarget state can distinguished by their labels [10]. Indeed, the key concept here is the assignment of a uniquely identifying track label to each kinematic target state x (i.e., x = (x, )). Moreover, this assignment, x = (x, ), can occur only once in the finite subset X [9], and the labeled finite subset is now denoted as X. By using this concept, the so-called labeled RFS proposed by Vo and Vo [18,19], the problem of the indistinguished targets can be solved [11]. Definition 1. [18] 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|.
It should be mentioned that the unlabeled version of a labeled RFS is obtained by simply discarding the labels. Unlabeled multitarget filtering formulation does not have the mechanism to address target tracks and, as a result, heuristic techniques are needed to produce target tracks [9].
In this paper, the unlabeled states are denoted by normal-faced letters to distinguish from labeled ones, which are denoted by the bold-faced letters (i.e., x, y, X, Y and x, y, X, Y), in which the single target state and multitarget state are denoted by lower-case letters and upper-case letters, respectively.
The spaces corresponding to variables are symbolized by blackboard bold letters (e.g., X, Z, L, etc.). The inclusion function 1 S (·) and the Kronecker delta function δ S (·) for a set S are given to support arbitrary arguments X (e.g., sets, vectors, and integers) as follows [19]: The class of finite subsets of S is represented by F (S), and the inner product is denoted by: The sequence of variables X i , X i+1 , . . . , X j is abbreviated as X i:j , and the cardinality of a finite set X is denoted by |X|.

Bayesian Multitarget Recursion
In classical Bayesian recursion, two assumptions are given, as follows [50]: (i) The hidden states follow a first-order Markov process on the state space according to a transition density f k|k−1 (x k |x k−1 ), and (ii) the observations are conditionally independent of the given states and are characterized by a likelihood g k (z k |x k ). By incorporating prior knowledge and observational evidence of a target, the Bayesian recursion has been formulated for the problem of single-target single-measurement systems.
Suppose that, at time k, there are n k targets i.e., the individual states are x k,1 , . . . , x k,n k , and the number of measurements is m k (i.e., the individual measurements are z k,1 , . . . , z k,m k ), then the multitarget state and multitarget observation are [19]: By using Bayes theorem, it could be seen that, conditioned on the measurement history Z 0:k = (Z 0 , . . . , Z k ), all the information on the set of target trajectories X 0:k = (X 0 , . . . , X k ) can be captured by the multitarget posterior density, which is given recursively for k ≥ 1 as follows: in which g k and f k|k−1 are the multitarget likelihood function and multitarget transition density to time k, respectively. By using g k and f k|k−1 , the underlying models for detection and false alarms and those for target motions, births, and deaths are encapsulated [19].
Assuming that, at the previous time step k − 1, the multitarget state is distributed according to a multitarget density π k−1 (·|Z 1:k−1 ). At time k, the measurement Z k is the superposition of detected targets and false measurements and is modeled by g k (·|·). The multitarget prediction and the multitarget posterior to time k is given by the Chapman-Kolmogorov equation and Bayes rule, respectively: Noting that the integral above is a set integral: defined for any function f : F (X × L) → R.

Multitarget State Model
For the purpose of RFS-based multitarget tracking, the RFS formulation of the standard multitarget dynamic model which captures all targets, including appearance, disappearance, and evolution ones over the time, will be investigated. Let X k−1 represent the set of all individual target states at time k − 1. Each target state x k−1 can survive and evolve to a new state x k at time k with survival probability P S,k (x k−1 ) or disappear between the time duration of k − 1 and k with probability 1 − P S,k (x k−1 ), not mention to some new-appear targets at time k (see Figure 2). The complete multitarget state X k at time k is the superposition of survivals and births [18]: in which S k|k (x k−1 ) and B k are the labeled RFS of existing target state from time k − 1 to time k and that of the new born states at time k, respectively. The labeled Bernoulli RFS S k|k−1 (x k−1 ) at current time k is generated by a given state where f S (S|X) is the distribution function of the surviving target set S the next time, and: The set B k of new-born states is distributed according to [19]: where ω B and p B are given parameters of the multitarget birth density f B, defined on space X × B. Whenever the set Y contains an element y with L (y) / ∈ B, f B = 0. The birth model (Equation (6)) includes both labeled Poisson and labeled multi-Bernoulli [18].
The multitarget transition density is given by [19]:

Multitarget Observation Model
Given a target state x k = (x, ), the sth receiver can detect this state with probability p (s) D (x, ) and generates a measurement z at time k is the superposition of the detected targets and Poisson clutter with intensity function κ [10], meaning that the information collected by radar includes false alarms. The probability of target detection, therefore, has a significant effect on the performance of the tracker. This conclusion is indeed true to RFS-based recursive Bayesian multitarget filtering.
Assuming that condition on X, detections are independent of each other and clutter, the multitarget likelihood function of sensor s is given as follows [18,19]: where M (s) is the number of measurements from sensor s. The map θ (s) specifies which object generated which detection from sensor s, i.e., target generates detection z θ( ) ∈ Z (s) with undetected targets assigned to 0. θ (s) is 1-1 on : θ (s) ( ) > 0 . The positive 1-1 property means that each distinct label is only mapped to a distinct positive value. Therefore, this property ensures that any detection in Z (s) is assigned to, at most, one target.
Followed by Reference [52], the multisensor likelihood is given by: where: N is the sensor observation set. It could be seen that the form of the multisensor likelihood function and that of the single sensor likelihood function are identical. Basically, a GLMB density can be rewritten in the following form [18,19]: where each ξ = (θ 1:k ) ∈ Ξ represents each history of multisensor association maps, each ω (I,ξ) is non-negative satisfying ∑ ξ∈Ξ ∑ I⊆L ω (I,ξ) = 1, and each p (ξ) ( ) is a probability density on X.
The cardinality distribution, the existence probability, and probability density of track ∈ L are, respectively, p (x, ) = 1 Given the GLMB density (Equation (11)), an intuitive multi-object estimator is the multi-Bernoulli estimator [52]. Given a prescribed threshold existence probability, this estimator determines the set of labels L ⊆ L which have higher existence probabilities the threshold, then the states of the objects will be estimated by using mode/mean estimates from the densities p(·, ), ∈ L.

Adaptive to Unknown Detection Probability
Normally, the probability of detection is assumed to be known a priori, however, this assumption is impractical. Therefore, the parameter p (s) D in Equation (9) much be estimated simultaneously with the update steps. It has been illustrated in Reference [28] that, with a specially chosen state space, the GLMB filter can be applied to solve the problem of tracking targets under the unknown probability of detection [47].
In this work, a new method called BpD-GLMB for tracking multitargets was proposed to estimate the unknown probability of target detection and produce target tracks online. Specifically, the unknown probability of detection parameter p D was estimated separately, then bootstrapped into the state-of-the-art GLMB filter for target tracking (see Figure 3). Inspired by Reference [53,54], this paper proposes using the pD-cardinalized PHD filter [11] to separately estimate the probability of detection, and then bootstrapped this parameter into the update stage of GLMB filter.
The idea of accommodating the unknown and non-homogeneous detection probability by incorporating this parameter into the target state variable has been proposed in Reference [11], in which each usual kinematic state x is replaced by an augmented statex = (a, x), where 0 ≤ a ≤ 1 is the unknown target detection probability of x. Obviously, the augmented state encompasses both the original kinematic target state x and the corresponding unknown and state-dependent detection probability (represented by a). Consequently, the augmented multitarget state has the form: x 1 ) , . . . , (a n , x n )} , (22) where n is the number of the target state, and the corresponding set integral has the form [11]: Since the probability of detection is unknown a priori, the filter should estimate this parameter such that the measurements are well adapted with the underlying target/object model [47]. By using the augmented state, the augmented space is given by: X =X × X, whereX = [0, 1] and X = R n x denote the spaces of detection probability and target kinematics, respectively. The augmented detection probability and single target likelihood function are described as follows:p Here, it is assumed that regardless of undetectability, a target will generate the same measurement. For the purposes of using RFS filters to tackle unknown probability of detection, the simple substitution of variables has been proposed in Reference [11] as follows: While keeping the likelihood function unchanged, the target state x and detection probability p  (a, x) and a, respectively. Consequently, the integral ·dx is substituted by 1 0 ·dadx.

Implementation
Because the number of terms in a GLMB updated multitarget density grows exponentially with time, practically only the terms with largest weights should be retained to limit the exhaustive evaluation of all terms. This work minimizes the L 1 approximation error [19] and can be formulated as a ranked implementation multi-dimensional assignment problem.
Although the K-short path algorithm can be used to compute the best terms of the prediction, it is easy for tracking loss unless the number of predicted terms is large enough. An alternative algorithm to mitigate this problem is the use of unscented transformation [55]. By using this algorithm, the predictions for target births and survivals are performed separately and then combined afterward [19]. For the purpose of reducing the cost of computation in the update state, measurement gating and pruning are applied in the GLMB filter. For the single-sensor case, two techniques, called the Murty's ranked assignment algorithm and Gibbs sampling, have been proposed to perform the truncation without having to propagate all the components [19,22]. While the Murty's algorithm can be used to determine a given number of highest weighted components of the multi-object filtering density without exhaustively generating all possible mappings, the Gibbs sampler can generate the significant components of the multi-object filtering density for a large number of targets to be tracked. Since the problem of tracking with Doppler measurements in this work requires multiple sensors, both the Murty's algorithm and Gibbs sampler implementation should be considered. However, the implementation of the two sensor GLMB filter developed in Reference [56] using Murty's algorithm showed that it has a cubic complexity in the product of the number of measurements from the sensors. Based on the extension of the Gibbs sampler implementation to multiple sensors proposed in Reference [52], Gibbs sampling is chosen as the most appropriate implementation for this problem. As a proof of concept of how GLMB filter addresses multistatic Doppler measurements, the simpler "iterated corrector" implementation that applies single sensor updates once for each sensor, in turn, has been used. This strategy would yield the exact solution if all components of the multitarget filtering density are kept.

Numerical Studies
In the present work, the problem of tracking 10 nonlinear birth-and-death time-varying marine ships was considered under the missed-detection-and-clutter observations and unknown detection probability. Consider the scenario illustrated in Figure 1, where we adopted the generic constant-turn model to better accommodate the unpredictable maneuvering behaviors of targets. The target state at time k is modeled using a 5-D vector x k = [p x k ,ṗ x k , p y k ,ṗ y k , ψ k ] T , comprising of its x-coordinate, x-velocity, y-coordinate, y-velocity, and course. Hence, the marine ship dynamic model can be expressed as follows: where: . (27) In which, ∆ is the sampling period, and n k = [nẋ k , nẏ k , n ψ k ] T is a zero-mean Gaussian noise vector of velocities and course noise components. σẋ = σẏ = σ v is the standard deviation (std.) of the velocity noise, and σ ψ is the std. of the course noise. Note that latitude, longitude, and course are measured in degrees ( • ), while distance, velocity, and time are calculated in nautical miles (M), i.e., knots (kn), and hours (h), respectively.
Remark 1 : The target model with transition matrix F (ψ) given in Equation (27) is based on the assumption that the surveillance area is located far enough from the North and South Poles.
The target model parameters, as well as the birth parameters, are given in Table 1. The multitarget scenario is conducted in the surveillance region [10 • N 30 • N, 100 • E 125 • E] with a total of 10 targets, which are random in positions, velocities, and the number of appearance (as in Table 2). Target ground truths and tracking results are depicted in Figure 4. By using BpD-GLMB filter, the results of 2-D coordinates target tracks are shown in Figure 5. The birth process is assumed as labeled Multi-Bernoulli where the common existence probabilities ω   Ten targets were assumed to be distributed around the birth model with the closest and farthest latitudinal distances being 2.85 km and 10.73 km, and the corresponding values for longitudinal distances being 2.6 and 10.6 km, respectively. The absolute velocities of the targets were assumed to be varied from 2 to 32 kn (approximately 3.5 to 60 km/h). The assumptions on positions of the multistatic passive Doppler radar transmitter and receivers and transmit frequency f t are given in Table 3.   The measurement space for each receiver is 200 Hz], and the measurement noise w k is zero-mean Gaussian noise with the covariance Q k = diag([0.5 Hz 2 ; 0.5 Hz 2 ]). The value of unknown detection probability p D,k (x) is first estimated by the corresponding pD-CPHD filter, and then bootstrapped into the GLMB filter. Clutter follows a Poisson RFS with average rates c 1 and c 2 , as mentioned in Table 3. It can be seen that the tracker can precisely track the targets in latitude and longitude with respect to the time. There are some delays and missed detection in observations, which may be due to the distances from the estimated positions and the actual positions. The tracker shows the effectiveness of the tracking algorithm when the targets are merged or close together.
In this paper, the tracking errors of the filters are evaluated and compared using both the Optimal Sub-Pattern Assignment (OSPA) [57] and the OSPA−on−OSPA, or OSPA (2) [36] with cut-off parameter c 0 = 100 and p = 1. The OSPA metric is a distance between two sets of points that jointly account for the dissimilarity in the number of points and the values of the points in the respective sets. By using OSPA metric, the errors between the true and estimated multitarget states at each time step is calculated and shown in Figure 6. The OSPA−on−OSPA, or OSPA (2) , distance has a different interpretation to that of the traditional OSPA distance. The OSPA (2) metric used in both GLMB and BpD-GLMB to capture the errors between the true and estimated sets of tracks over a period of time with window length set at l = 20 is illustrated in Figure 7. In this paper, we assumed that the GLMB with known p D = 0.98 as the optimal filter. Further, a comparison of estimated cardinality between the GLMB and BpD-GLMB is shown in Figure 8. The results demonstrate the capability of our proposed BpD-GLMB filter, which outperformed the pD-CPHD filter, and its comparability to the GLMB filter with known p D in tracking accuracy under OSPA and OSPA (2) metrics.

Conclusions
An online multiple targets tracking filter with multistatic Doppler measurements was given in this work. The unknown probability of detection parameter, which was assumed to vary slowly compared to the data rate, was estimated and bootstrapped into the state-of-the-art GLMB filter to output target tracks. The simulation results demonstrate the effectiveness of the proposed method, the BpD-GLMB, which is comparable to the optimal-GLMB with a known detection probability while surpassing pD-CPHD significantly. Since the application only involves a small number of targets, the Murty's algorithm was applied for the prediction and update performance of the tracker. For significantly faster implementation with a larger number of targets scenario, the joint prediction and update introduced in Reference [52] using Gibbs sampling could be used. Evidence has been shown in Reference [36] that the GLMB filter can be used to track over one million targets per scan in real time based on a C ++ platform. For strategies of realistic large-scale tracking implementation, MATLAB code is obviously insufficient, and it should be converted into C ++ or Python. The problem of tracking mobile targets using movable sensors will be further investigated for upcoming consideration.
Author Contributions: C.-T.D. devised the project, the main conceptual ideas, and proof outline. He also worked out the technical details and performed the numerical calculations for the numerical study. He wrote the manuscript and responded to the reviewers' questions. H.V.N. helped in correcting the visualization of numerical study and consulted the first author in responding to reviewers' questions.
Funding: This research was funded by CIPRS/MOET Scholarship.