Multiple extended target tracking by truncated JPDA in a clutter environment

Data association is a crucial part of target tracking systems with clutter measurements. In general, its complexity increases sharply with a number of targets and measurements. Recently, high ‐ resolution sensors have given rise to extended target tracking problems and more than one measurements can emerge from each target, making the association problems more complex. In this study, a tractable algorithm based on the Gaussian process measurement model and truncated joint probabilistic data association technique is proposed for multiple extended target tracking in the presence of measurement origin uncertainty. Based on the marginal association probabilities, the calculation amount is effectively reduced by truncating the association events with low probabilities in the shortest path problem. The effectiveness of the proposed algorithm is verified by the test of multiple extended targets tracking and compared with the linear ‐ time joint probabilistic data association as well as the algorithm on random finite sets. Simulation results show that the proposed algorithm can track multiple extended targets accurately, which is significant in high ‐ resolution radar tracking systems.


| INTRODUCTION
Early target tracking algorithms usually assume that each target generates at most one measurement in the vicinity of its real location, indicating the applications of the point target (PT) model. In PT model, only the kinematic states of targets are considered. Recently, with the emergence of high-resolution sensor technology, the PT model has no longer been suitable for modelling the target since multiple measurements have been generated according to both the kinematic states and shape of each target. The target that can simultaneously generate multiple measurements is denoted as an extended target (ET). The extended target tracking (ETT) represents an essential part in many transportation systems, including the driver assistance systems [1], crowd surveillance [2,3], autonomous driving systems [4] and maritime traffic control [5]. In these applications, both the kinematic and the extent (size and shape) states of the targets need to be estimated, thus, a suitable extended target model is required.
In recent years, various types of target models have been proposed to describe the ETs [6,7]. In Ref. [8] bicycle target was modelled as a thin stick determined by its orientation and length. An ellipse is a simple shape that can be parameterised by three factors: major axis, minor axis and orientation. In Ref. [9] symmetric positive definite random matrices (RMs) were used to describe the extent of targets, and the RMs were assumed to be inverse Wishart distributed. Assuming that multiple measurements can appear on the scaled version of target contours, elliptical targets can also be modelled with elliptic random hypersurface models (RHMs) [10]. As for nonelliptical or irregular targets, each target can be regarded as a combination of multiple ellipsoidal sub-targets, which leads to a multiple random matrices model [11,12]. An alternative way is to model an ET with a radial function that represents the distance from the centre of motion (COM) to the target boundary for each given angle. Using the star-convex RHMs [13], target radial function can be expressed by Fourier series expansion in the frequency domain. Another idea is to consider the target boundary as a curve determined by a number of control points. In Ref. [14], B-splines were used to model a target with an arbitrary extent, without the restrictions that targets should be multiple elliptical or star-convex.
The Gaussian process (GP) is a stochastic process with an infinite dimensional Gaussian distribution. The GP is capable of learning the unknown function by a training set, and it has been widely used in machine learning [15,16]. In Ref. [17], based on the assumption that radial function values in different angles have been multi-dimensional Gaussian distributed, a GP-based measurement model was applied to the single ETT to estimate the kinematic and extent states of an ET simultaneously. By combining the GP measurement model with the labelled multi-Bernoulli (LMB) filter [18], the multiple ETT can be realized. A generalized version of probabilistic data association (PDA) based on the GP was derived in Ref. [19] to achieve the single ETT with measurement origin uncertainty. In Ref. [20], the extended target Gaussian process probabilistic multi-hypothesis tracker (ET-GP-PMHT) was proposed to track multiple PTs and ETs seamlessly. Recently, the Gaussian processes based multiple model joint probabilistic data association (GP-MM-JPDA) filter [21] was used to track a maneuvering ship, where the shape of the ship's wake was modelled by the GP.
Data association algorithms are essential when measurements from targets as well as various types of clutter (generated by clouds, terrain) are received by sensors. The nearest neighbour (NN) is the simplest algorithm that selects only one measurement with the maximum likelihood and discards the others, but it does not perform well at a high false alarm rate or a low detection probability. Another data association method is to consider all the valid measurements with different weights [22]. The probabilistic data association filter was proposed by Bar-Shalom in Refs. [23,24] to track a single target in a clutter environment. In Ref. [25], the information about target existence was also taken into account, and the probabilities of target existence and association events were calculated.
In multiple target tracking (MTT), measurements from multiple targets and clutter are obtained by sensors. Thus, data association in the MTT is a challenging problem since it is much more complex than in a single target scenario. Random finite sets (RFSs)-based method is an innovative idea that can avoid data association problem and has been widely used in the MTT. Early works about RFSs include the probability hypothesis density (PHD) filter [26,27], cardinalized PHD (CPHD) filter [28], multitarget multi-Bernoulli (MeMBer) [29] and cardinality balanced MeMBer (CBMeMBer) filter [30]. Recently, following the theory of the conjugacy prior to [31] in the Bayesian framework, the generalized labelled multi-Bernoulli (GLMB) filter [32][33][34] and Poisson multi-Bernoulli mixture filter [35,36] have been proposed to track multiple (extended) targets. In the traditional data association algorithms, such as probabilistic multi-hypothesis tracking (PMHT) [37,38], data association is realized using measurements obtained by multiple scans. In other words, the association results depend on both current and historical measurements. Different from the PMHT, the joint probabilistic data association (JPDA) [39,40] methods calculate the probability of each association event according to the current measurements only and the estimation results represent the weighted sum of posteriors conditioned on each association event.
This study proposes an improved multiple star-convex ETT algorithm based on the GP measurement model in the presence of measurement origin uncertainty. In the existing JPDA algorithms, the gate techniques are used to exclude the false alarms and measurements that only lie in the validation gates are considered in the update of the target states. However, the computational cost is still too high due to the increase in the number of targets and measurements. In this study, in view of the Poisson spatial measurement model [41]-based linear-time calculation of marginal association probabilities [42], and a modified version of the data association technology called the truncated joint probabilistic data association (TJPDA) is derived by selecting a number of association events with maximum probabilities. Then, estimates conditioned on each selected association event are calculated, and the result represents the normalized weighted sum of these conditional estimates. That is , the proposed TJPDA only calculates the marginal association probabilities, whereas the existing JPDA needs to consider the probabilities of all the possible association events.
The article is organised as follows. In Section 2, the extended target model based on the Bayesian framework is discussed, and a brief introduction in the Gaussian process regression (GPR) and the single ET measurement model based on the GP is given. The main idea of the traditional JPDA and the calculation of GP-based marginal association probabilities, which leads to the derivation of a truncated joint probabilistic data association algorithm for multiple ETT, are presented in Section 3. The simulation results are evaluated in Section 4. Finally, the conclusions are summarised in Section 5.

| PROBLEM FORMATION
This section presents the non-linear filtering technique that is used in this work and introduces the extended target model. In addition, the GP-based single extended target measurement model [17] in a two-dimensional space is presented.

| Extended target model
In target tracking, the Bayesian framework usually consists of two parts, the target dynamic model and the sensor measurement model, which are, respectively, expressed as follows: where f ð⋅; ⋅Þ and hð⋅; ⋅Þ represent the target state transition and measurement functions, w k−1 and v k denote the process and measurement noises, respectively. The state transition is processed by the Chapman-Kolmogorov equation as follows: where pðx k−1 jz 1:k−1 Þ denotes the posterior at time k, pðx k jz 1:k−1 Þ represents the predicted density and f kjk−1 ðx k jx k−1 Þ denotes the state transition density. The Bayesian formula is given by: where pðz k jx k Þ represents the measurement likelihood function, and pðx k jz 1:k Þ represents the posterior density at time k.
In this work, target states at time k are expressed as follows: where x k;c and x k;c 0 represent the position of the COM and its derivative, respectively, x k;φ denotes the target orientation, and its derivative is denoted as x k;φ 0 , x k;e describes the extent of target, which represents the radial function values for a series of given degrees uniformly distributed in the interval of ½0; 2π�. Assuming that the transitions of kinematic and extent states are independent of each other, the target dynamic model can be defined as follows: Where F ki and F e represent the transition matrices of kinematic and extent states, Q and R denote the process noise covariance values of the kinematic and extent states, and they are ,respectively, expressed as follows: where ⊗ denotes the Kronecker product, σ 2 p and σ 2 φ are the process noise variance values of position and orientation, respectively, M represents the dimension of the target extent state, α denotes a parameter for extent transition, kð⋅; ⋅Þ represents the covariance function of the GP and finally, T denotes the sampling time.
The Kalman filter is a special case of the recursive Bayesian formula, which is derived according to the minimum mean square error (MMSE) and under the assumptions that target dynamic and measurement models are linear Gaussian. In addition, the non-linear techniques, such as extended Kalman filter (EKF), unscented Kalman filter (UKF) and particle filter (PF), can be used to develop non-linear models. In this study, the EKF is applied to deal with non-linear problems.
In the EKF, non-linear models are linearised by the firstorder Taylor formula. The target dynamic model and sensor measurement model are, respectively, expressed as follows: The corresponding transition equations are given by: The measurement updating process of EKF is given as follows:

| Gaussian process regression
The Gaussian process is a stochastic process that can be seen as an infinite dimensional Gaussian distribution, which can be expressed as follows: where μð⋅Þ and kð⋅; ⋅Þ represent the mean and covariance functions, respectively, μð⋅Þ is usually set to a zero constant function. The purpose of the GPR is to learn an unknown LI ET AL.
-209 function using a training set. Consider the measurement model as follows: Suppose, there are a series of noisy measurements z ¼ ½z 1 ; :::; z m � T given by θ ¼ ½θ 1 ; :::; θ m � T , and that the aim is to learn the function values Then, the distribution of f p conditioned on z can be obtained by: it is a n-dimensional Gaussian distribution with mean Az and covariance P, where Covariance function is a very important factor in learning as it implies the assumptions on the unknown function f ð⋅Þ. The square exponential (SE) covariance function is the most common covariance function, and it is expressed as: The correlation coefficient between f ðθÞ and f ðθ 0 Þ increases as θ and θ 0 become closer to each other and vice versa. The SE covariance function can be modified to a periodic covariance function [17] as follows: where σ 2 f and σ 2 r represent the variances of the signal amplitude and Gaussian prior, respectively, n defines the period as follows T ¼ nπ. It should be noted that 2 is a common value of n in modelling the target extent because radial function for any starconvex shaped target has a period of 2π. By setting n to 1, the symmetric covariance function [17] that assumes that the target extent is central symmetric can be obtained. However, targets mostly tend to be axisymmetric, and the COMs are usually located on the symmetry axes, which leads to the axisymmetric function [18] that is expressed as: In Refs. [17,18,43], it has been demonstrated that more accurate estimates can be obtained by choosing appropriate covariance functions.

| Single extended target measurement model
Equations (20∼22) can be used to predict the values of the function f ð⋅Þ by conducting a number of noisy measurements. However, in order to make GPR applicable in the recursive Bayesian formula, it is needed to model the noisy measurements with a series of predictions, which leads to the single ET measurement model.
Since the radial function f ð⋅Þ represents the distance between the COM and the target boundary for a given angle, the jth measurement of a target at time k in two dimensions can be written as follows: where f ðθ L k;j Þ represents the predicted radial function value for a given angle θ L k;j ¼ θ G k;j − x k;φ conditioned on the priors, θ G k;j is calculated using the relative position between the COM and measurement z k;j , p k;j is the rotating vector expressed as ½cosðθ G k;j Þ; sinðθ G k;j Þ� T , e k;j represents the sensor measurement noise, which is a zero-mean two-dimensional Gaussian distribution with covariance R k;j and finally, s k;j denotes the scaling factor. Note that z k;j represents the contour measurement when s k:j equals to one, and when 0 ≤ s k;j < 1, it represents the surface measurement. The distribution of f ðθ L k;j Þ can be derived using the GPR as follows: By substituting (27∼29) into (26), the GP-based single ET measurement model can be obtained as follows: |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } þ e k;j |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where R is determined by sensor measurement noise, H k;j represents the measurement matrix for z k;j and e tr k;j denotes the true noise in the measurement model. The covariance of e tr k;j can be easily obtained based on the assumption that the GPR noise and sensor measurement noise are independent of each other, which is expressed as follows: when the scaling factor s k:j is assumed to be normally distributed as follows: The measurement function (30) is modified [17] to: An example of contour measurement z k;j , where s k:j ¼ 1, is shown in Figure 1. Two coordinate systems, the global and local coordinate systems, are constructed in the GP measurement model, such that the sensor is located at the origin of the global coordinate system, and the COM represents the origin of the local coordinate system, and the target orientation x k;φ decides the direction of the x-axis of the local coordinate system.

| GP-TJPDA ALGORITHM
The linear-time JPDA algorithm presented in Ref. [42] calculates the marginal association probabilities with the linear time complexity. These marginal probabilities only illustrate the association probability of each measurement to each target, which cannot be directly applied to the measurement update in the GP-based ETT algorithms. In this section, a GP-TJPDA algorithm for tracking multiple extended targets in a clutter environment using a number of association events with highest weights is derived. A diagram of the GP-TJPDA algorithm is shown in Figure 2. Firstly, the marginal association probabilities are calculated when receiving the current measurements. Then, the truncation process is used to obtain a series of association events with maximum probabilities and finally, the results can be obtained. It can be seen from Figure 2 that the main part of this algorithm is the truncation process. In addition, different from the traditional JPDA, in the GP-TJPDA, the measurement-to-target association event probabilities are calculated without using the gating techniques. The proposed algorithm follows several assumptions as below.
1. Each target has an extent and measurements from each target are randomly distributed along the target contour. The cardinality of measurements from each target is Poisson distributed. 2. Measurements from clutter are uniformly distributed in the surveillance region, and the number of clutter measurements follows the Poisson distribution. 3. All the measurements are conditionally independent of each other, when conditioned on the predicted target states.
The traditional JPDA algorithm and the calculation process of marginal association probabilities are presented in the following.

| Traditional JPDA
Define the set of measurements and target states at time k as follows: where M k and N k denote the numbers of measurements and targets, respectively. In traditional JPDA algorithm, before the arrival of the current measurements, a validation gate is built for each track. The jth measurement can be associated with the ith target only if it falls within the validation gate of the ith target, which is expressed as: measurement matrix of the ith target for the jth measurement, P kjk−1;i denotes the predicted state covariance matrix of the ith target, and finally, σ represents the gate threshold. The JPDA has the same estimation equations as the PDA, and by using the total probability rule, the posterior of the ith target can be approximately calculated bŷ where θ denotes a M k dimensional vector that represents the association event; θ j ¼ i when the jth measurement is associated with the ith target and θ j ¼ 0 when the jth measurement is considered as a clutter, Y k denotes the predicted target state Y k ¼ fx kjk−1;i ; P kjk−1;i g N k i¼1 ; x θ k;i and P θ k;i represent the posterior state and covariance conditioned on the association event θ, respectively; and r θ k represents the probability of θ. Note that in the traditional JPDA, θ cannot have repeated values greater than zero since each target is assumed to generate at most one measurement per time step.
The difference between the JPDA and PDA lies in the calculation of association event probabilities. In the PDA [24], all of the valid measurements lie in the same validation gate, and each association event assumes that at most one of the valid measurements is generated from one target and all of the measurements that are not associated with the target are considered to be Poisson-distributed clutter. However, in the JPDA, some targets can be close to each other (these are targets whose validation gates have intersections) and measurements can simultaneously lie in different validation gates, which assumes that one measurement can relate to different targets in different association events. The detailed process can be found in Ref. [39].

| Marginal association probabilities calculation
As presented in Ref. [41], in the association problems of ETs, Z k is assumed to consist of N k þ 1 conditionally independent inhomogeneous Poisson point processes (PPPs), where the clutter measurements are assumed to be uniformly distributed with a density in the detection region and the number of clutter is Poisson-distributed with a parameter λ k;0 . Measurements generated from the ith target are Gaussian-distributed and the number of measurements of the ith target is Poisson-distributed with parameter λ k;i , so the likelihood function of the jth measurement for the ith target can be expressed as: where G ij denotes the Gaussian distribution when i > 0 and uniform distribution when i ¼ 0. The clutter density ρ ≜ 1=V is decided by the detection volume V .
Define the notations for an association event θ as follows: •Z θ k;i -a set of measurements associated with the ith target.
•Z θ k;0 -a set of measurements considered as clutter. By using the Bayesian formula, the probabilities of an association event θ conditioned on the current measurements and predictions in (40) can be divided into two parts as follows: The cardinality of each measurement set T θ ¼ fn θ is known when an association event θ is given. Therefore, the former part in (43) can be expressed as: F I G U R E 2 A diagram of the GP-TJPDA, note that the number of association events considered in the filter is related to the number of targets, which is where pðθjT θ ; Y k Þ represents the association probability conditioned on the cardinalities of each measurement set, and it can be obtained by: Based on the PPP hypothesis, the probability of cardinalities is calculated as follows: where λ k;i denotes the measurement rate of the ith target at time k and λ k;0 represents the clutter rate. By substituting (45) and (46) into (44): Since measurements are assumed to be conditionally independent of each other, the latter part in (43) is given by: Substitute (47) and (48) into (43), the probability of association event θ is obtained by: Therefore, the marginal association probability can be obtained [42] as follows: As demonstrated in Ref. [42], the time complexity of marginal association probabilities calculation is OðN k M k Þ, and it has a linear relationship with both the number of targets and the number of measurements.

| Truncated JPDA algorithm
Equation (49) can be used to calculate the probability of each association event. However, there are ðN k þ 1Þ M k possible association events in a total of N k given targets and M k measurements, so it is almost impossible to enumerate all the possible events when there are a large number of targets and measurements. In Ref. [19], the traditional PDA was extended to a generalized version that is suitable for single ETT. However, this generalized PDA is not suitable to be directly modified to the JPDA and used in multiple ETT. Although the gating techniques can be used to reduce the number of unnecessary association events, the computation will become intractable when plenty of measurements lie within the valid gates, especially when multiple measurements lie in the intersections of multiple valid gates. Therefore, a truncation technique is required for computational convenience.

| Marginal association probability matrix
Assume the elements of matrix denote the marginal association probabilities, such that the element in row i þ 1 and column j of A represents the probability of the jth measurement associated with the ith target, so it holds that where the zeroth target is set to be clutter. Therefore, an association event θ can be expressed by selecting one element from each column in A (θ j ¼ i for the (i+1)th element of the jth column in A), and its probability is the product of M k selected numbers, which is expressed as follows:

| Truncation process
In the following, a truncation process that can obtain a number of association events with highest probabilities in multiple ETT, is proposed. The number of possible association events expressed as ðN k þ 1Þ M k can be reduced to ðN k þ 1Þ by the truncation process. Although Murty's algorithm was used to extract a series of most likely hypotheses in MHT [44], this algorithm cannot be applied to ETT because of the constraint that each target can generate at most one measurement. The main idea of the proposed truncation process is to obtain the approximate estimates without computing all the conditional LI ET AL.
-213 posteriors and the corresponding weights using (40) and (41), which can be achieved by solving a problem similar to the shortest path problem. The truncation problem is depicted by a directed graph that is shown in Figure 3.
The cost matrix C can be defined based on the marginal association probability matrix A as follows The distance from node l ij to node l i 0 j 0 in Figure 3 can be calculated by Each path from the start node S to the end node E corresponds to a unique association event θ, which has a total distance of P M k j¼1 Cðl θ j j Þ. Since cost c ij and marginal association probability a ij are negatively correlated, a path with a lower distance will correspond to an association event with a higher probability. Assume the cost matrix C consists of N k þ 1 rows and M k columns (that includes N k targets and M k sensor measurements). Then, the steps of the truncation process are as follows: Step 1. Choose the first two columns, Cð:; 1Þ and Cð:; 2Þ of C and compute a ðN k þ 1Þ � ðN k þ 1Þ matrix by Q ¼ Cð:; 1Þ⊕Cð:; 2Þ; set i ¼ 3.
Step 2. Select N k þ 1 minimum elements in Q, record their positions and build a new vector p consisting of these N k þ 1 selected numbers.
Step 3. Update matrix Q by Q ¼ p ⊕ Cð:; iÞ and set i ¼ i þ1. If i ¼ M k þ 1, go to Step 4; otherwise, go to Step 2.
Step 4. At this point, Q contains distances of at least N k þ 1 shortest paths. Use all the recorded positions in the above steps to obtain association events with the highest probabilities.
The computational cost mainly lies in the Step 2, which is OððN k þ 1Þ 2 logðN k þ 1ÞÞ. Since Step 2 will be repeated M k − 1 times, the calculation of the truncation process is OðM k N 2 k log N k Þ. In contrast, the calculation for version without truncation is OðN M k k Þ.
Note that in this study, ⊕ denotes an operation that generates a matrix consisting of the sums of the elements of two vectors in pairs, such as the following matrix: � In order to illustrate the difference between TJPDA and JPDA algorithms, an example is given below. In Figure 4, there are four targets and 20 measurements, where 12 measurements lie in the validation gates. It is noted that three measurements are in the intersections of gates from two targets, and other three measurements are in the intersection of gates from three different targets. For the JPDA algorithm, there are 2 6 � 3 3 � 4 3 association events need to be numerated, however in the proposed TJPDA algorithm, only five association events with highest probabilities has to be obtained without using the gating technique.

| GP-TJPDA filter
After the introduction of the EKF based on the GP measurement model and the truncation process in the measurement-to-target associations, the pseudo-code of the GP-TJPDA filter is given in Algorithm 1.

| SIMULATION RESULTS
The performance of the GP-TJPDA was verified by several simulations and the proposed algorithm was compared with the linear-time JPDA [45] and RHM-PHD [46] algorithms, respectively. All of these three algorithms were tested in tracking multiple ETs with the clutter measurements and missed detection. The performances were evaluated by three metrics: position error, shape error and run time. Optimal subpattern assignment (OSPA) metric [47] and intersection over union (iou) were used to evaluate the position and shape errors, respectively. Each simulation was performed 50 Monte Carlo runs by using Matlab R2011b on an Intel(R) Core(TM) i5-7500 3.40 GHz PC with 8 GB RAM.

| Evaluation metrics
Given two subsets X ¼ fx 1 ; :::; x m g and Y ¼ fy 1 ; :::; y n g, the OSPA metric can be calculated by where.
where d(x,y) denotes the distance between x and y, and ∏ n represents a set of all permutations; the cut-off parameter c defines the punishment degree of the cardinality error in the OSPA metric, and p represents the order parameter. Note that the OSPA metric can be regarded as the root mean square error (RMSE) when X and Y have the same cardinality and p = 2. In this study, it is set as c = 20 and p = 2.
The shape estimation performance was evaluated by the iou. For two sets of regions R ¼ fR 1 ; :::; R m g, R ¼ fR 1 ; :::;R n g corresponding to X and Y, the IOU of these two sets is defined as: The IOU for two regions is given by where R i andR πðiÞ represent the regions of the true target and target estimation, respectively. In (57), πðiÞ implies that in the calculation of the IOU, it is required to match the region of each true target with that of the estimated target using the results of the optimal point assignment from OSPA metric.

| Star-convex targets
In the first scenario, two triangular targets moved with a combination of the constant velocity (CV) and constant turn rate and velocity (CTRV) models [48] in a surveillance region of ½-10; 90� � ½-30; 70� m 2 , where a sensor was located at the origin. Both two targets existed between time 1 s and 600 s and the trajectories are presented in Figure 5. The initial kinematic states were as follows: The initial covariance of kinematic state was P i;0 ¼ diað10 −2 ; 10 −2 ; 10 −5 ; 10 −5 ; 10 −5 ; 10 −10 Þ. The extent state of each target was initialised to an eighty-dimensional vector with the values of three, which represented a circle with a radius of 3 m.
As given by (6∼10), the CV (with ω ¼ 0) model was used in the prediction of the kinematic state, where the process noises were σ p ¼ 10 −2 and σ φ ¼ 10 −3 , and the parameter α ¼ 10 −4 . The covariance function proposed in Ref. [43] with parameters σ f ¼ σ r ¼ 2 was used. The target measurements were generated randomly along the target boundaries with the sensor measurement noise R ¼ diagð0:01; 0:01Þ, and clutter measurements were uniformly distributed in the surveillance region. The measurement rates of each target and clutter were set to 15 and 30, respectively.
The tracking results of the three targets at time 480 s and 520 s are presented in Figure 5, and the enlarged views are shown in Figure 6.
The average values of the OSPA distance and iou are depicted in Figure 7 and Figure 8, respectively. Although the linear-time JPDA slightly outperformed the GP-TJPDA algorithm in the position error, the IOU for the GP-TJPDA was much higher than the linear-time JPDA at all time steps.
As shown in Figures 7 and 8, both algorithms performed well when two targets met with each other. The iou of the GP-TJPDA had significant improvement because the proposed algorithm can track any star-convex target, whereas the lineartime JPDA is limited to elliptical targets. The average CPU times of the GP-TJPDA and linear-time JPDA are listed in Table 1.
Linear-time JPDA has an advantage in runtime because it uses only three parameters to model a target extent. Thus, the proposed algorithm can obtain more detailed target shapes compared with the linear-time JPDA.

| Elliptical targets
In the second experiment, three elliptical extended targets were tracked in a surveillance region of ½-300; 1300� � ½-300; 1100� m 2 . The target trajectories were the results of the CTRV model, which are depicted in Figure 9. All the targets had the same size with the major axis of 9 m and minor axis of 7 m. Each target was initialised to a circle with a radius of 4 m and the initial kinematic states were as follows: Unlike the first scenario, the CTRV model was used in the prediction. In addition, the periodic covariance function k 2π ðθ; θ 0 Þ given by (24) was used.
The detailed estimation results of the three targets at time 200 s are shown in Figure 10. The average value of the OSPA and iou metrics are depicted in Figure 11 and Figure 12, respectively. Please note that for the results of the RHM-PHD algorithm, the 'pulses' that can be observed in Figures 11 and  12 were the consequence of the incorrect estimation for target number caused by clutter measurements. Similar phenomenon did not appear in the GP-TJPDA and linear-time JPDA algorithms because the proposed GP-TJPDA is not a comprehensive algorithm that considers the newborn and dead targets like the PHD filter.