I Introduction
Maneuvering target tracking is a fundamental task in sensorbased applications, such as radar, sonar, and navigation [LIU2016183, 6507656, 1413764]. Based on the Bayesian framework, target tracking is usually solved via a recursive update of the posterior probability density function (pdf) of the target states. The application of Bayesian filters is based on dynamical motion models and sensor measurement models [1413764, 5730505]. However, there may be significant motion model uncertainty when targets undergo unknown or mixed maneuvering behaviors [1561886, JLYang], such that the evolution of the target state is too complex to be approximated by a reasonable number of mathematical models. The uncertainty in motion behavior can also be due to various parameters in generative models not being known a priori or if they are timevarying. The tracking performance of Bayesian filters would degrade and even become unacceptable if incorrect models or parameters are applied. This paper focuses on the issue of learning uncertain motion models with higher flexibility and incorporating the learned models into the Bayesian framework for target tracking.
Ia StateoftheArt in Bayesian Filters
MM methods are commonly used to deal with the model uncertainty problem [1561886]. The basic idea of multiple model (MM
) methods is to use a bank of elemental filters with different motion models to capture the mixed motion behaviors. Then, the overall estimate is generated based on the results achieved by each elemental filter
[1561886]. For the cases when the target’s behavior is a timevarying or doublystochastic process, [JLYang] and [LiuZX] incorporated adaptive parameter estimation methods into the Bayesian filter framework to deal with the uncertainty of unknown parameters, where the unknown parameter is estimated by approximating its distribution with particles and corresponding weights. MM methods may be limited in performance if the assumed trajectories are not modeled accurately enough to capture all aspects of the motion behaviors, and the estimation of unknown parameters results in increased computation complexity.IB StateoftheArt in Machine Learning Approaches
As an alternative to traditional Bayesian methods, machine learningbased tracking algorithms have been proposed recently
[8099886, 7088657, 8096087, 9011258, Batch_nonlinear, 9185002, 9272174]. In [8099886], a quadruplet convolutional neural network (CNN) based algorithm was designed to learn object association for multiobject tracking via 2D image processing. An extended target tracking method using a Gaussian process (GP) was proposed in [7088657] and [8096087] to track the irregularly shaped object whose moving trajectory follows a fixed motion model. These papers use a GP for shape modeling but not for capturing the motion behavior. A long shortterm memory (LSTM) neural network was used in
[9011258] to perform the prediction step of target localization. The proposed filter with LSTM prediction can improve the tracking accuracy by using computationally efficient lowdimensional state spaces but only considers the singletarget tracking (STT) scenario. In [9011258], the continuousdiscrete estimation of the target’s trajectory is viewed as a onedimensional sparse GP, with time as the independent variable and measurements acquired at discrete times. Recursive GPbased approaches for target tracking and smoothing were developed in [Batch_nonlinear, 9185002], with online training and parameter learning. Only STT with a linear observation model (i.e. axis and axis coordinates) is considered. Furthermore, it was assumed that the motion model could be decomposed into separate models for the X and Y axes. In [9272174], the motion behavior of a single target is modeled as a Markov process which switches between the constant velocity (CV) and coordinated turn (CT) models with the specified transition probabilities. The complex mixed model is learned as joint
GPs for position and velocity prediction. Then the trained GPs are integrated into the PF for tracking the maneuvering targets in realtime.The learningbased tracking algorithms proposed in [Batch_nonlinear, 9185002, 9272174] were for motion models based on Cartesian positions which are problematic because positions are not translationally invariant, i.e. the motion model needs to be constantly updated as the target moves. This paper proposes learning a naturally shift invariant motion (NSIM) model. Specifically, the change in the target’s Cartesian positions on the axis and axis in one instant (i.e. Cartesian velocities), instead of the Cartesian positions, is learned as multiindependent GPs.
IC Proposed Method and Contributions
This work follows a learningtotrack framework for tracking problems with uncertain motion models and proposes a new learning method with which the motion behaviors of targets are learned as GPs based on a training data set. Then we show how the learned motion models are integrated into the PF and belief propagation (BP) to realize target tracking. The three main contributions are:

We propose a new GPbased dynamical motion model learning method. The NSIM behavior of a target is learned using Gaussian process regression (GPR) based on the Cartesian velocities. The learned GPs are translationally invariant and can be applied to track objects within different surveillance regions from the surveillance region of the training data.

We explore the potential of using the learned motion model within the PF [978374, article] for STT and PF based BP [910572, koller2009probabilistic, 6301723] for multitarget tracking (MTT) to demonstrate that the proposed algorithm is a plugandplay learning method in target tracking systems.

The proposed learning method is firstly tested in STT with different motion scenarios to demonstrate its capability to learn uncertain motion models with higher flexibility. Then, the learned motion model is plugged into MTT, which shows significant improvement in tracking performance compared with interacting multiple model (IMM)PF methods. Additionally, offline learning is performed only once and then applied for MTT. Hence, the online tracking time is not affected by the learning process.
Compared with the available learningbased tracking methods using target positions, the motivations and advantages of the proposed GPbased NSIM model include:

The information available for use in the tracking problems may be limited depending on how we envisage getting the training data. For example, if the data was taken from other previously tracked targets, it could be challenging to access speed or acceleration.

The existing methods learn the motion models based on target positions, which are not translationally invariant and need to update the model as the target moves constantly. Moreover, the tested targets’ surveillance regions may differ from the training data’s surveillance region. On the other hand, the motion of an object is usually controlled by the application of a force, e.g. aircraft control surfaces. As such, the motion is translationally invariant. Therefore, the proposed NSIM model based on Cartesian velocities will be translationally invariant and can reduce the requirement of training data size and algorithm complexity.

The learned GPbased NSIM model can provide confidence intervals for particle sampling in the prediction step of the
PF to provide a reliable estimate of model uncertainty.
ID Paper Organization and Notations
The remainder of this paper is organized as follows. The problem solved in this investigation is defined in Section II. Section III reviews the GPR and proposes the GPbased learning method of NSIM behavior. The proposed algorithm, which incorporates the learned GPs
into the target tracking, is developed in Section IV. Simulation results and performance comparisons are presented in Section V. Conclusions are provided in Section VI. This paper presents matrices with uppercase letters, vectors with boldface lowercase letters and scalars with lowercase letters. The notations are defined and summarized in Table I.
Notation  Meaning 
Target states of training data sequence  
Observations of training data sequence  
Target states of test data sequence  
Observations of test data sequence  
Training time index  
Test tracking time index  
Target index  
Cartesian positions at 2D plane  
Cartesian velocities at 2D plane  
Ground speed  
Velocity heading angle  
Alongtrack acceleration  
Crosstrack acceleration  
Covariance matrix of process noise  
Covariance matrix of measurement noise  
Training data set  
Motion function  
Measurement function  
axis velocity motion function for training  
axis velocity motion function for training  
axis velocity motion function for testing  
axis velocity motion function for testing 
Ii Problem Definition and Background
In this section, we will introduce the system model and conventional Bayesian approaches of target tracking.
Iia System Model
The dynamic statespace model (DSSM) for target tracking is given by the following motion and measurement models:
(1)  
(2) 
Here, represents the time index, , is the target index, . For STT scenarios, . The state motion model and the measurement model are represented as and , respectively, where , , and are dimensions of the state, process noise, observation, and measurement noise vectors, respectively. The process noise and the measurement noise are assumed to be additive white Gaussian noise (AWGN). At each time instant, the sensor node detects each target with detection probability . Owing to missed detections and clutter, the number of detected measurements at time may not be equal to . is defined as the number of measurements at time and is timevarying. The set of measurements detected at time is denoted as . The measurement index is represented as and . The number of clutter points is modeled as a Poisson point process with intensity . The corresponding clutter measurements are independent and identically distributed (iid) with pdf .
IiB Sequential Bayesian Filters for Target Tracking
Based on the Bayesian theorem, most target tracking methods are sequential Bayesian filters which consist of the prediction and the update steps. In the prediction step, the prior distribution of targets’ states at time is calculated via the Chapman–Kolmogorov equation as [8290605]
(3) 
In the update step, the joint posterior pdf of the targets’ states is calculated as [8290605]
(4) 
For STT in nonlinear systems, the PF uses random particles with associated weights to approximate the posterior pdf as
(5) 
where is the Dirac delta function, the number of particles used at a given time is represented as . The particles are drawn from an importance density which can be set as [978374].
For MTT, the association between targets and measurements is described by two data association (DA) vectors shown in Fig. 1 [6978890, 7889057], i.e. the targetoriented association vector and the measurementoriented association vector where is the number of measurements at time . The DA vectors are defined as [6302287]

if at time , the th target generates the th measurement.

if at time , the th target does not generate a measurement, i.e. missed detection occurs.

if at time , the th measurement is hypothesized to be associated with th target.

if at time , the th measurement is hypothesized not to be generated by a target, i.e. clutter.
The objective of MTT is to calculate the marginal posterior pdf by
(6) 
where the marginal association probability, i.e. , approximates based on the loopy BP algorithm, as given by Eqn. (22) in [6978890].
It should be noted that most of the existing STT/MTT filters are fully modeldriven and, therefore, of use when the DSSM is explicitly available. However, assuming that the state motion function is given as a simple deterministic model is restrictive. These filters are intractable for some applications where the motion model is very complicated. Instead, information about the state transitions can be inferred from examples of preceding states and current states pairs. Hence, we assume that the motion model is learned using the nonparametric Gaussian process regression (GPR) [books/lib/RasmussenW06]. The measurement model is deterministic by different sensor types with no correlation over time and is in a known parametric form, i.e. we know the acquisition model. In the next section, we will discuss the application of GPR to learn motion behavior.
Iii Gaussian Progress Regression for Motion Behavior Learning
GP models are widely used to perform Bayesian nonlinear regression. A
GPmodel takes a training data set and a kernel as the input and outputs a joint Gaussian distribution that characterizes the value of the function at some given points
[phdthesisgp]. One practical advantage of GP is having confidence intervals for predictions that provide a reliable estimate of uncertainty. In this section, we briefly review the standard theory of GPR and then introduce the proposed GPbased method for motion behavior learning.Iiia Gaussian Process Regression
The GP model for regression can be considered as a distribution over a function and is defined as [7088657]
(7) 
where and represent the input variables. The mean function and the covariance function of two variables and are defined as [8882261]
(8)  
(9) 
Let us consider a general GPR problem with noisy observations from an unknown function described as
(10) 
Here,
is the variance of the noise. The training data set is denoted as
, where are the inputs with the corresponding outputs . The purpose of GPR is to derive the latent distributions of the vector at the test inputs , conditioned on the test measurements and the training data set. Specifically, the joint distribution of the training measurements
and the function at one test point, i.e. , is assumed to be Gaussian and is given by(11) 
Here, , is the covariance of , and denotes the covariance matrix for the training input data [9272174]. The conditional distribution of the function is then derived from the joint density in (11) and written as the following standard forms:
(12a)  
(12b)  
(12c) 
where is the abbreviation for . The most widely used kernel function is the squared exponential (SE) covariance function [books/lib/RasmussenW06], which is also used in this paper,
(13) 
Here, the hyperparameters
and represent the prior variance of the signal amplitude and the kernel function length scale, respectively.IiiB Gpbased Motion Behaviors Learning — Limited Information of Training State
This paper proposes a new GPbased learning method. I.e., the naturally shift invariant motion (NSIM) behavior of a target is learned as GPs based on the Cartesian velocities. The Cartesian positions and are converted into the Cartesian velocities as and , respectively. The NSIM models, i.e. and , are learned as GPs, where . As most target maneuvers are coupled across coordinates, the learned models capture the crucial characteristics of the target behavior during maneuvers by using both Cartesian velocities as the input of the GPs. Specifically, the motion behavior of the targets in (1) is rewritten as
(14a)  
(14b)  
(14c)  
(14d) 
The proposed GPbased method for prediction shown in Fig. 2 involves offline training and online prediction steps. The training data is used offline to capture the statistical characteristics of the motion behavior of the targets, which are then transferred to online tracking for state prediction.

[leftmargin=1.5mm]

Offline training: To train the Cartesian velocity motion functions, we take training data sets and . Specifically, the training inputs for both GPs are the same as
while the output data are different and set as
Here, denotes the number of training samples. The Cartesian / velocity motion functions in the offline training process, i.e. and , are modeled as two GPs shown in (15) and (16).
(15) (16) Here, and are the input variables, the squared exponential (SE) defined in (13) are used to calculate and with the same hyperparameters and . Examples of the learned GPs are shown in Fig. 3, where the motion behavior of the target follows the gradual coordinated turns (GCT) model, defined in Section V. The turn rate of the target follows the left and coordinated right turns (/s for 10s). In Fig. 3, the and axes represent the input of the GP, i.e., Cartesian velocities, and the Z axis represents the prediction value. The training samples are shown in blue points. The orangeedge and green surfaces represent the mean value of the predictive Gaussian distribution for the axis and axis velocity, respectively, which are joint functions of current Cartesian velocities.

Online Prediction: In the following, we describe the estimation of the axis position . The estimation of the axis position is done similarly and so is not discussed. The estimated parameters at time include the Cartesian position and velocity estimates represented as . The predictive distribution of the tracking state is partitioned as
(17) Here, . The pdfs in (17) are further calculated in (18) – (19),
(18) (19) where the predictive mean and covariance are calculated in (20) and (21), respectively.
(20) (21) where and are the abbreviations for and .
Iv Incorporation of Learned Gpbased Nsim Model into Target Tracking
For nonlinear and nonGaussian motion and observation models, the pdf of the hidden states in (6) and (17) cannot be evaluated in closed form. Therefore, particlebased implementation algorithms are commonly used for target tracking. The proposed GPbased NSIM model learning method is plugged into the PF and PFbased belief propagation (BP) algorithms in this section. When the parametric motion model is unavailable, particles are drawn based on the learned GPs. The learned GPbased PFMP implementation for MTT is summarized in Algorithm I. The GPbased PF for STT is implemented the same way without executing Step 5, i.e. the loopy BP. Note that the data association for MTT is obtained by running an iterative BP algorithm following [6978890]. Therefore, the proposed algorithm is still valid when the target number is unknown and time varying.
Iva Draw Particles based on Learned GPs
The particle vectors for the th target at the preceding time slot are represented as
with the corresponding posterior weights , . The particle vector consists of the Cartesian position particles on the axis and axis, i.e. , and the Cartesian velocity particles on the axis and axis, i.e.
. The posterior probability of
th target’s positions and velocities at time is approximated as(22) 
The new particles at time are drawn based on the importance densities represented as . Specifically, the particles of the Cartesian velocities are sampled from the Gaussian distributions based on the learned GPs as
(23)  
(24) 
Here, is the training data. The mean and covariance are calculated based on (20) and (21), respectively. Next, the position particles and are sampled according to (25) and (26) as
(25)  
(26) 
Then, the prediction of the position state at time is approximated as,
(27) 
where .
IvB Data Association
In the PFbased BP, the marginal association probabilities, and , are calculated based on the loopy BP, i.e. updating the messages between factor nodes and , as shown in Fig. 1. Specifically, the two sets of messages are represented as and to denote the message passed from node to node and the message passed from the node to the node , respectively. The messages are updated as [6978890]
(28)  
(29) 
where the factor indicates the association probability calculated based on the likelihood ratio (LR) as
(30) 
where when no measurement is correctly associated to th target, and the predictive distribution is approximated as (27). The messages are passed between the latent association variables until convergence, i.e. (28) and (29) are repeatedly processed until the conditions of convergence is satisfied [DBLP:journals/jmlr/IhlerFW05]. Finally, the approximations of the marginal association probabilities are obtained as
(31)  
(32) 
IvC PFbased State Estimation
After these marginal association probabilities in (31) and (32) are obtained, the marginal posterior pdf in (6) is then approximated as
(33) 
where the updated weights are calculated by
(34)  
(35) 
Based on the key idea of the PF [978374, 1232326], the maximum a posteriori (MAP) estimation of hidden states for the th target at time , i.e. , are approximated as
(36) 
where represents the estimates of the th target’s position, and represents the estimates of th target’s Cartesian velocities between and .
Based on the above descriptions, the implementation of the proposed learned GP based PFBP filter for target tracking is summarized in Algorithm I.
V Simulation Results
This section introduces the dynamic statespace model (DSSM) used to generate simulation data and different tracking maneuvering scenarios in the first two subsections. Then, the performance comparisons for STT and MTT are presented. Specifically, the proposed learned GP method is applied for STT over different maneuvering scenarios. The tracking performance demonstrates the validity and robustness of the proposed training method to track the targets within different surveillance regions and motion behaviors. Then, the datadriven MTT algorithm is compared with the oracle PF with fully known system information and the interacting multiple modelparticle filter (IMMPF) methods that know the correct DA [1561886]. The numbers of the used particle are set to be the same for all compared filters.
Va Dynamic Statespace Model (DSSM)
This paper generates the simulated target trajectories and measurements following the standard curvilinearmotion kinematics model [1261132]. The state of th target at time in (1) and (2) is represented as , where are the target Cartesian position, and represent the ground speed and the velocity heading angle, respectively. The alongtrack and crosstrack accelerations in the horizontal plane are denoted as and , respectively. By setting different values of