1. Introduction
In the evolution of Multi-Target Tracking (MTT) technology, foundational methods centered with data association, including Probabilistic Data Association (PDA), Joint Probabilistic Data Association (JPDA), and Multiple Hypothesis Tracking (MHT), have been crucial [
1]. These methods decompose the multi-target tracking problem into simpler single-target tracking problems through data association. The Bayesian multi-target filter, which propagates the multi-target posterior density over time, offers an optimal Bayesian solution for multi-target tracking challenges. The theory of Random Finite Set (RFS) has been widely applied in fields such as radar tracking [
2], air traffic management [
3], and autonomous driving [
4], constructing optimal Bayesian filters. Representative algorithms include the Probability Hypothesis Density (PHD) filter, the Cardinalized PHD (CPHD) filter, and the Cardinalized Balanced Multi-Target Multi-Bernoulli (CBMeMBer) filter. These foundational theoretical filters, such as the PHD and CPHD filters, directly propagate the posterior first moment [
5]. However, they encounter challenges in tracking trajectories.
For modern high-resolution millimeter-wave radars and lidars, the traditional point target assumption is no longer applicable, as multiple measurements can be generated at each time step. Therefore, Extended Object Tracking (EOT) has become a mainstream challenge in applications [
6]. For elliptical extended targets, commonly used measurement models include the random matrix model [
7], set cluster process [
8], physics-based model [
9], and Poisson space models [
10]. Among these, the Poisson space model is the most widely used. In this model, target detection is modeled by a non-homogeneous Poisson Point Process (PPP). The PPP measurement likelihood has a simple factorization that avoids the explicit representation of assumptions about the association between measurements and targets [
11]. At each time step, a Poisson distribution of random measurements is generated, spatially distributed around the measurements of the target. To address the issue of tracking multiple targets, the RFS framework is used for modeling. The RFS framework offers a systematic approach to Multi-Target Tracking (MTT) by accommodating changes in the number and state of targets over time. Both targets and measurements are modeled as random finite sets, and the PPP extended target model has been utilized in several computable complex models. The ensemble approach within this framework allows for an optimal Bayesian representation of the birth and death of targets. According to the literature, the state-of-the-art algorithm for tracking an unknown number of elliptical targets is the Poisson Multi-Bernoulli Mixture (PMBM) filter. This method acknowledges that the number, size, shape, and kinematic states of targets are interrelated and influence each other. However, due to the numerous prior conditions and associations that need to be considered, such associations present computational challenges.
According to [
12,
13], the most advanced algorithm for tracking an unknown number of extended targets is the PMBM filter. For estimating target trajectories, MTT methods based on random vectors link the current state estimates of an object with its past states, or represent the emergence of new targets. In the RFS-based MTT methods, one effective approach is to assign a unique label to each trajectory, allowing each object to be identified as it changes over time [
14]. Another effective method is to calculate the multi-target posterior probability over the set of trajectories, which contains all information about multiple trajectories and does not require the use of labels [
15]. In multi-target tracking scenarios, the Trajectory Probability Poisson Multi-Bernoulli Mixture (TPMBM) filter and its approximate filter [
16], the Trajectory Probability Poisson Multi-Bernoulli (TPMB) filter, have emerged [
17].
For methods based on gating and sampling, the assumption that negligible weights can be truncated avoids the exhaustive data association often required for large datasets [
18]. Therefore, in scenarios with high uncertainty in data association, performance degradation can occur due to information loss. To prevent the loss of information associated with explicit enumeration and local association assumptions, the method of computing the multi-target marginal posterior for data association is used, which marginalizes the unknown data association uncertainties [
10]. Currently, the best-performing multi-EOT algorithm, as proposed in [
19], is the Sum-Product Algorithm (SPA) based on particle Belief Propagation (BP), which is implemented through a factor graph of the joint posterior of multi-target states and measurement association variables. Simulation results in [
20] demonstrate that the TPMB filter based on factor graphs outperforms the TPMBM filter that uses gating clustering and truncated sampling in terms of estimation errors and runtime. The BP algorithm recursively iterates to obtain global hypothesis data association information, transmitting this association information between factor nodes and variable nodes via a directed graph. In the process of approximating PMB with PMBM, Kullback–Leibler Divergence (KLD) minimization is used to marginalize the uncertainty of data hypotheses, thereby reducing the amount of data processing. This process has led to the development of a closed-form filter that initiates the recursive process from the joint posterior factor graph model of object sets and association variables, along with the message-passing equations.
This paper presents a model-data-driven message-passing algorithm that adaptively adjusts the motion models of targets in different dynamics through an Interactive Multiple Model (IMM) to adapt to the switching motion state at different time steps [
21]. In the message-passing process using factor graphs, considering the measurement association factor’s message as an approximation of the posterior probability requires computation for each message-passing variable, which is complex. The proposed algorithm reduces the complexity of data association and minimizes message passing by gating hypothetical surviving targets, thereby enhancing algorithm efficiency. In the filter, a new particle filtering method is introduced that constructs the importance density function for each particle’s state via the Adaptive Square Root Cubature Kalman Filter (ASCKF), optimizing the covariance matrices of process noise and measurement noise to improve the tracking performance of high-maneuver targets under noise and clutter influences [
22,
23]. This paper is an important extension of [
19,
22], including the following contributions:
In addressing the complexity of model transition probability calculations in the IMM algorithm, we utilize changes in the posterior probability between two consecutive moments to reflect the transition probabilities between models. On the other hand, the standard IMM algorithm considers process noise as a constant during model switching, affecting the accuracy of the filtering algorithm. The improved IMM algorithm adaptively adjusts the magnitude of process noise, reducing errors in model switching.
In addressing the complex message passing between variable nodes and factor nodes in factor graphs, we have improved the posterior association step of the factor graph. The message-passing process is streamlined through adaptive gating, eliminating the need to consider all inter-node relationships when computing association nodes, thereby reducing running time.
To address the issue of particle degradation observed in the methods using particles, we employ the SCKF to perform particle importance sampling and use the Gate-control Belief Propagation (GBP) for the improved IMM-TPMB-GBP filter, updating the particle filter based on these sampling results. Due to the decrease in tracking accuracy of adaptively modeled maneuvering targets under noise and clutter, we have developed the ASCKF. This involves adaptively estimating the covariance of process and measurement noise, ensuring that the algorithm’s performance does not depend on noise statistics and is more robust when these statistics are unknown.
The rest of the paper is organized as follows.
Section 2 introduces the background on state and trajectory parameters, the trajectory set, and TPMB.
Section 3 provides the BP cycle for TPMB marginal posterior probabilities and the principles of factor graphs, further detailing the adaptive gating improvement for factor graphs and the ASCKF mixed particle filtering method. The implementation of ASCKF for unknown noise is presented in
Section 4. Simulation results are provided in
Section 5, and the conclusion is given in
Section 6.
2. Background
In this section, we first introduce the state and trajectory parameters of targets, then model the multi-target extended model state parameters and measurements, and finally present the TPMB approximation of the PMBM posterior density.
A finite spatial set is D, which is denoted as , and the cardinality of its subset is . We use ⨄ to denote the union of disjoint sets, to denote the inner product , and the multi-target index is represented as . For certain real-valued functions, the product is conventionally represented. Additionally, and , respectively, denote the Kronecker delta centered at , and the Kronecker delta function. Here, in represents a subspace within the Cartesian product space, .
2.1. State and Trajectory Variables
At time step , the single-target state includes its kinematic state and any potentially identified target extension states. At time step , the set of target states is a random finite set, and . The measurements collected by sensors at time step consist of the jth single-target measurement produced at time step and . The sequence from the initial moment to the current time is denoted as .
A target’s motion trajectory is its finite sequence of states over continuous time steps, represented as . Here, is the initial time step of the trajectory, is the duration of the trajectory steps, and is the finite sequence representing the target state trajectory. The parameter variable belongs to the set . Thus, the single-target trajectory reaching a finite time step at moment belongs to space , with ⨄ representing the union of disjoint sets. The trajectory set up to the current moment can be represented as . The single trajectory function can be decomposed into , where is defined on and .
2.2. Multi-Target Models for Sets of Trajectories
Multi-target dynamic models: At time step , existing targets survive to the next time step with probability , or perish with probability . If a target survives, its Markov state transition function is . At each time step, newly detected targets are independent of other surviving targets and follow a PPP with a birth intensity of .
Multi-target measurement models: The measurement set for targets is divided into two parts: one consists of actual measurement sequences generated by extended single-target detection by sensors, modeled as PPPs with measurement likelihoods of produced by each and with a Poisson intensity of . The other consists of clutter measurements induced by measurement noise, modeled as a PPP with a Poisson intensity of , a Poisson rate of and clutter intensity of .
Dynamic model for sets of all trajectories: For the set of surviving trajectories
at time step
, each trajectory
has a survival probability
. If the trajectory continues to survive in the next time step, its transition probability density is
Newly initiated trajectory sets
independently follow a PPP with an intensity of
Given the complete set
of trajectories, both surviving and newly initiated, at time step
for each
,
, the survival probability is 100%, i.e.,
, and the trajectory’s transition density is
Thus, the current moment’s survival probability of a trajectory depends solely on the final state of the trajectory in the previous moment and is independent of other moments. The probability density of the measurement sequence produced by a single trajectory is .
2.3. Standard Bayesian Model for Sets of Trajectories
By extending an auxiliary variable
over the single-target trajectory space,
is used. When
represents that the target trajectory was not detected, it corresponds to a PPP; when
, it indicates that the target trajectory corresponds to the
ith Bernoulli component. After adding the auxiliary variable, we use
to represent the set of target trajectories with the auxiliary variable.
denotes the weight, ⨄ denotes the union of disjoint sets;
is the Poisson part probability density function, representing the trajectories of objects assumed to exist but undetected;
is the Poisson intensity for trajectories;
is the multi-Bernoulli part probability density function, representing the probability of latent trajectories being detected at least once by some time
; in the multi-Bernoulli probability density Equation (5),
indicates the number of Bernoulli components that exist up to some time
, with each Bernoulli component having
independent hypotheses;
represents the global hypothesis index. For the
ith Bernoulli component density, there are
local hypotheses and
. For a global hypothesis
, a local hypothesis
is assigned to each Bernoulli component, and
is the set of global hypothesis indices; that is,
.
where
represents the weight of the local hypothesis Bernoulli component, and
represents the weight of the global hypothesis Bernoulli component. Normalization ensures that the sum of the probabilities of the global Bernoulli components equals 1; that is,
.
2.4. PMBM Conjugate Prior
TPMB is a special filter that approximates the TPMBM posterior probability density. In this approximation, trajectory-oriented MBMs are merged into a single MB by introducing the parameter variable, thus handling the mixture of densities. The preferred method for approximating PMBM is minimizing the Kullback–Leibler divergence (KLD).
Given the multi-target probability density defined in Equation (5), after the inclusion of parameter variables, the multi-target probability density
can be represented as follows:
For the given surviving targets
, the undetected targets
are included, and their predicted multi-target probability density after incorporating the parameter variables is the following:
where the auxiliary variable
may represent multiple undetected trajectories, and it does not necessarily imply that
equals zero.
2.5. Trajectory PMB Approximation
According to reference [
16], the TPMBM approximates the TPMB method, and this paper also utilizes special multi-trajectory states and measurement models. According to Equation (10), we derive the TPMB prediction equation.
Given the multi-trajectory set prediction density
and measurement set
at time step
, from Equation (5) to Equation (8), it is known that the number of updated Bernoulli components
, and the updated Poisson intensity is as follows:
For each Bernoulli component predicted at the previous measurements
, the current moment’s new measurements generate
Bernoulli components. At the same time, each Bernoulli component predicted from the previous moment generates
new local hypotheses at the current moment. These newly generated local hypotheses may be false detection or updates for non-empty subsets of measurements
. After augmenting the model using auxiliary variables, it is necessary to define the global hypothesis set of generalized measurement model parameter variables. We give the augmented representation method of measurement
for
. For the Bernoulli hypothesis component
i, a set of measurement collections represented by a local hypothesis
is
:
where
.
For each Bernoulli component
in the predicted probability density
and the new measurements detected at the current moment generating
Bernoulli components, there are
local hypotheses corresponding to either a false detection or an update at the current moment. For the false detection hypothesis of the Bernoulli component
, it can be represented as
.
Consider that the measurement set is empty in Equations (12)–(15). We use
to represent the non-null measurement set
in Equations (16)–(18). For the predicted Bernoulli component
at time step
, it survives to the next moment and is associated with measurement
, which yields,
,
,
.
Then, we construct the set
to associate with the new Bernoulli component
, with
, and use
to represent the Bernoulli component of the measurement set
, where
,
,
. For the new Bernoulli components
generated by new measurements
at time step
, there are two single-target hypotheses for each Bernoulli component: one corresponds to a previously undetected and non-existent Bernoulli component, represented as
,
,
; the other corresponds to a Bernoulli component that existed in the previous moment.
The set represents the local hypotheses related to the Bernoulli component measurement indicators. After the Bayesian update, the updated density of PMBM is obtained through the aforementioned process. The established generalized measurement of the model’s trajectory posterior set, utilizing the posterior conjugacy of PMBM, provides convenience for using general nonlinear target measurement models.
3. Model-Driven IMM-TPMB Filter
The TPMBM obtains a computationally less intensive TPMB posterior density through KLD approximation. However, during the approximation process, the uncertainty in the marginal density of the PMB trajectory set increases, and the uncertainty of the global hypotheses has been marginalized. For methods based on gating and sampling, negligible hypotheses can be ignored by truncating weights to avoid exhaustive data association. Nevertheless, in scenarios with high data association uncertainty, performance can degrade due to information loss. Therefore, to prevent the loss of information from local association hypotheses caused by explicit enumeration, a data association method that computes the multi-target marginal posterior is used to marginalize the unknown data association uncertainty. The model established in this paper is based on the factor graph probability model in [
17,
18], which enables SPA to have good computational complexity and utilizes the accuracy and flexibility of an SCKF-based implementation.
We define the measurement-vector association factor for multi-target trajectories. Due to the uncertainty in the posterior marginal association probability for each measurement
, the indicator factor for the Bernoulli component generated by each measurement is defined as follows:
where
only when the measurement
is associated with the
ith potential target; when
, it indicates that the
jth measurement is a false alarm.
The joint posterior probability of the trajectory set and the measurement vector association factor can be factored into a factor graph. According to [
10], the joint posterior of the trajectory set at time step
and the measurement vector association model during the TPMB model update process is
The joint posterior expression consists of three parts: the first part represents the updated posterior probability density of the undetected trajectory set; the second part represents the joint posterior of the survival probability density of the Bernoulli component from the previous moment and the measurement association vector at this moment; the third part represents the joint posterior of the probability density function of the newly detected trajectories and the measurement association vector after receiving measurements.
According to Equation (24), the factor graph is depicted with circles representing factor nodes and squares representing variable nodes. For undetected trajectories, only an additional factor node and a variable node are used to represent them, independent of the factor graph’s message-passing process. For the surviving trajectories at time step , it is predicted that the next moment will have Bernoulli components, represented here by factor nodes. Among these factor nodes, an equal number of variable nodes are generated, which connect to the measurement vector association nodes as the posterior marginal probability density of this Bernoulli component. For the Bernoulli components generated by the measurements at time step , it is indicated that factor nodes are generated. In these factor nodes, each factor node produces a different number of posterior marginal probability density nodes.
In the factor graph message-passing model, as the starting variable nodes in the message-passing process,
and
correspond to the prior probability densities of undetected and surviving targets, respectively, while
represents the posterior probability density of undetected targets and newly detected targets.
In
Figure 1, the factor nodes
and
in the information propagation are represented as follows:
3.1. Implementation of an IMM-TPMB Filter Based on Two Consecutive Time Points
In multi-maneuvering target tracking scenarios, the motion states at different times may be described by multiple motion models, and the filtering results must consider the state outputs from various models. The IMM method is an effective solution for such problems. The IMM method uses multiple different models to model the target’s state, and each model could be the optimal one at the current moment.
Assuming the number of tracking models is
, the motion matching model at time step
is
, and at time step
it is
, with a model mixing probability of
, it can be represented as follows:
where the Markov transition probability from
to
at time step
for model
,
represents the probability of model
being the matching model at time step
, and
represents the normalization constant, which can be expressed as follows:
Using the initial state mean
of model
, the initial conditions for model
are calculated.
where
and
represent the mixed state matrix and mixed covariance input to the filter, respectively. The filtering algorithm is used for prediction and updating. After each time update, it is necessary to recompute the Markov matrix. However, the computation of the transition matrix is complex and susceptible to noise. To address this issue, this section adopts the posterior probability density of each model to calculate the transition probabilities
. As the transition probabilities are closely related to state changes, and the state changes are mainly determined by the posterior probability
of model
at time step
, the transition probabilities between models can be reflected using changes in the posterior probabilities between two consecutive moments.
Given the posterior probability density
of model
at time step
, the relative rate of change of probability between two adjacent moments
can be given by the following proportional form:
As indicated by Equation (33), when the probability density
of model
at time step
increases,
. This means that the posterior probability density from other models at
moments transitioning to the current model gradually increases; the converse is also true. Since the posterior probability density
changes over time,
can be represented by the model transition probability from the previous moment
and the relative rate of change of probability
. Assuming a normalization constant
, the update of the transition probability
can be represented as follows:
According to Equation (34),
is represented as the product of
and
, and there exists
, thus satisfying the Markov transition matrix conditions. Obviously, the above parameters can inherently match the maneuvering dynamics, including the CT motion and the motion switching. As a result, the maneuvering targets with different dynamics can be detected effectively.
On the other hand, the intensity of target maneuvers directly determines the magnitude of process noise. The standard IMM algorithm considers process noise as a constant, which affects the accuracy of the filtering algorithm. Moreover, the process noise changes with unknown noise influences. Therefore, it is necessary to adapt the IMM algorithm to adjust the process noise to accommodate scenarios where the target is maneuvering or non-maneuvering and the noise is unknown. When the target is maneuvering, increase the process noise
to compensate for modeling errors, resulting in smaller filtering estimation errors; when the target is non-maneuvering, reduce
to maintain tracking accuracy. Assuming
is the maximum value of the posterior probability density of each model at time step
, then the covariance matrix
of the process noise
at time step
can be represented as follows:
where
is the initial value of the process noise covariance matrix for model
. According to Equation (35), when the models described by
and
are the same,
matches the current target motion model, and
takes the minimum value; when
and
describe different models,
does not match the current target motion model, and
takes the maximum value. Considering the impact from the process noise, the adaptive adjustment of
ensures consistency with the varying intensity of target maneuvers.
3.2. Implementation of an IMM-TPMB Filter Based on Gate-Control
3.2.1. Predict
Undetected and newly detected targets: at time step
, the PPP prediction intensity for all undetected and newly detected targets within the sensor-monitored area is as follows:
where the probability density of the single-target particle belief during propagation is
, the particle weight indexed by
is
, and
only when
; otherwise, it is 0.
indicates the number of newly detected target particles, and the target component is in the form of particles
.
Surviving targets: for the detected surviving targets, the target state of the
ith MBM component at time step
is represented as follows:
In the case of undetected and newly detected targets, there is no need to consider the target model transition probability. represents the motion model transition function, which performs maneuvering target particle sampling through multi-model mixing prediction.
3.2.2. Update
In the BP iterative propagation algorithm, the optimization of the propagated information is achieved through P iterations, leading to the convergence of belief for all nodes. At this point, the label of each node aligns with the label of the node with the highest posterior probability density among the related nodes, representing the optimal label. In the update step, we improve the way of information transmission and reduce the amount of data association calculation by using the GBP. The specific improvement method is in step 3. The update step in model-driven scenarios includes six stages: particle initialization, importance sampling of particles, iterative information transfer, belief calculation, pruning and resampling, and target state estimation.
Particle initialization: The Poisson intensity is represented as a set of particles denoted by . Updating for new, missed, and surviving targets requires considering whether targets are alive at the current moment.
Particle importance sampling: Utilize the ASCKF for importance sampling of particles. After iterating the prediction and updating process of prior and posterior information at time step
, the particles at time step
follow the distribution function
Here, denotes a Gaussian distribution with mean and variance , as determined by the ASCKF.
Iterative information transfer: Measurement assessment requires calculating the information
transferred from each factor node
to variable node
in the vector measurement. Factor node
represents the process of information transfer between surviving trajectories and both sides.
In iterations, based on the factor graph formula, the running time of calculating particle weights during the information propagation between each variable node and factor node increases. As the number of targets increases, we also consider the propagation time between factor nodes and variable nodes. The proposed GBP is set to reduce the running time of false hypothesis associations. Therefore, the main computational cost is effectively saved by only achieving the related weights. As a result, the tracking accuracy is well ensured.
There is no need to consider the information propagated across all nodes; instead, only the information within the gating needs to be calculated, thus reducing the weight calculation. The likelihood function for target generation is represented as follows:
where
represents the Euclidean distance between the surviving target particles and the current measurement location, with
being the gating threshold for particle-measurement association. After one iteration, particle weights converge towards values close to the true state of the target. In subsequent iterations, the information
propagated from the variable node
to the factor node
represents the unnormalized Bernoulli component density.
When
,
, predict the next moment’s surviving target generation of particles and update with current measurements:
When
,
, predict the next moment’s missed and new target generation of particles and update with current measurements. The information
and particle weights are expressed as follows:
- 4.
Belief calculation: After multiple iterations of the above steps, the information and example weights converge to the true posterior probability density of the target. The belief for the target state in each MBM component is represented as a set of particles, and its calculation method is similar to information . Additionally, subsequent iterations still require calculating particle weights and the existence probability of Bernoulli to ensure the validity of trajectory belief, followed by weight normalization.
For new targets
in the detection area at time step
, when
, then the updated existence probability of the target is as follows:
For undetected potential targets
, when
, then the updated existence probability of the target is the following:
The updated particle weight is normalized as follows:
- 5.
Pruning and resampling: To reduce the number of global hypotheses, pruning and resampling of the MBM components are necessary. By approximating MBMs with low existence probabilities as PPPs, the algorithm increases computational speed without significantly impacting accuracy. The calculation of distortion caused by minimizing the KLD under the PPP approximation is as follows:
- 6.
Target state estimation: Select the MB components with existence probabilities above a set threshold to estimate the current target state at time step , from a set of weighted particle collections, represented as the outputted target trajectory . New target estimates are achieved using a Poisson point process, outputting a set of new target trajectories, denoted as trajectory , where is the birth time, and is the trajectory length.
4. IMM-TPMB Filter Driven by Model and Data Collaboration
Addressing the particle degeneration issue identified in [
20], this paper proposes a SCKF-based particle importance sampling method for IMM-TPMB-GBP. The improved IMM-TPMB-GBP particle filter employs the ASCKF to obtain state estimates and their covariance matrices at each moment, using these as the mean and variance in the original filter to establish a proposal distribution function. Innovatively, in scenarios where noise statistics are unknown, the ASCKF provides an adaptive noise ASCKF, enhancing the accuracy of particle weights. While particle filters typically face degeneration, the ASCKF performs moderately in systems with non-Gaussian noise estimation. Utilizing the advantages of ASCKF in estimating unknown noise and nonlinear estimation, this paper combines the improved ASCKF with IMM-TPMB-GBP, proposing a new state estimation method. Further, the ASCKF in this paper adaptively estimates the covariance of process noise and measurement noise in interacting multiple models, not relying on noise system metrics, thus offering robust performance.
We omit the auxiliary variables
of the target state and use the node
i to represent the highest posterior probability density. In a nonlinear model, the target state equation and the measurement equation for a motion model
are formulated as follows:
where
denotes the input matrix, assuming the position and velocity components are
and
, respectively,
represents the sampling interval,
can be expressed as
,
is process noise with zero mean and the covariance matrix
.
represents the measurement set at time step
,
is the measurement noise with zero mean and the covariance matrix
.
4.1. Predict
Initialize the state mean
and the error covariance
and decompose the Gaussian terms’ covariance in the target trajectory composed by Equation (44) through Cholesky decomposition
.
Use the third-order spherical-radial volume criteria to generate volume points
, where
. The volume points after propagation are yielded through the nonlinear state equation:
Calculate the predicted state mean and square root covariance from the volume points after propagation:
where
represents the triangulation operation. The volume points after propagation are yielded through the measurement equation:
Calculate the measurement mean and the square root covariance matrix:
Calculate the cross-covariance between the state and measurement:
Calculate the Kalman gain:
4.2. Update
Calculate the updated state mean and the square root of the state covariance:
During the model update, calculate the probability of the
model:
Calculate the conditional probability:
As known from Equations (61) and (66),
and
significantly influence the values of
and
. However, in the target tracking process, the noise data change over time. Any mismatch between the actual noise affecting the system and the noise assumed in the SCKF can degrade the performance of the SCKF, as it may also diverge. Therefore, it is necessary to precisely estimate the matrices
and
. The posterior density function can be calculated as follows:
where
is considered a constant dependent on prior information, and
does not involve an optimization problem; hence, Equation (69) is equivalent to the following:
where
can be calculated using the conditional probability multiplication theorem as follows:
where
is a constant,
represents the dimension of process noise, and
can be expressed as follows:
where
represents the dimension of measurement noise and
is a constant. By substituting Equations (71) and (72) into Equation (70),
can be calculated:
where
. According to Equation (73), calculate the partial derivatives of the function
with respect to
and
, where the extremum points are as follows:
According to the calculation process of the ASCKF,
and
can be represented as follows:
By solving the first-order partial derivatives of the function
with respect to
and
, a suboptimal estimate of the unknown covariance matrices of process and measurement noise can be obtained, which is then incorporated into the nonlinear filter ASCKF for recursive computation.
Table 1 gives the algorithm process. Combined with the algorithm proposed in
Section 3, this method does not depend on noise statistics and exhibits good robustness. Usually, the detection probability and clutter rate are two important factors that affect tracking performance during the process of target estimation. Therefore, we analyze measurement noise and process noise individually, and then obtain the modification of
and
in the ASCKF. Both the over-estimation and the under-estimation derived from the factors can be corrected. On the other hand, the IMM filtering framework ensures the motion switching between two continuing dynamics, which can distinguish the close targets. Especially, these targets can be further distinguished by using the proposed filter based on the novel filtering mechanism.
This algorithm replaces the covariance matrix recursion of the CKF with the square root form of the covariance matrix, ensuring the efficiency of the filtering process while avoiding numerical instabilities caused by rounding errors, thus enhancing the precision and stability of the filtering algorithm. Moreover, ensuring the positive definiteness of the covariance matrix enhances its numerical stability. The performance of the proposed method, compared to CKF and UKF, features smaller estimation errors and faster convergence rates, and it effectively tracks the motion state of maneuvering targets, enhancing the system’s accuracy.
5. Simulation Results and Discussions
In the numerical study, a typical 2-dimensional scenario is performed to evaluate both the reliability and efficiency of the proposed filter, which is driven by model and data collaboration for extended target tracking. The experimental environment was the following: IntelTM CoreTM i5, RAM 16 GB, WindowsTM 10, and MATLABTM R2022b.
5.1. Simulation Scenario
This algorithm primarily addresses the tracking of extended targets under the conditions of multi-model switching and unknown noise influences. In this section, we present the results of 100 Monte Carlo simulations conducted in the first scenarios. We tested the proposed IMM-TPMB-GBP-ASCKF filter in two target motion scenarios and compared it with the GGIW-TPMB filter from [
24], which used clustering and assignment, as well as with the TPMB-BP [
20], TPMB-GBP, and IMM-TPMB-GBP-CKF filters.
The linear programming (LP) metric was used to compare the performance of these filters [
25,
26]. Thus, we set different numbers of maneuvering targets in two scenarios, each under different motion models. The target motion models include Constant Velocity (CV) and Constant Turning (CT), with their corresponding motion equations as follows:
Each object consists of a two-dimensional position and velocity, which we represent using the target matrix . denotes the plane position coordinates, and is used to represent the planar position. We consider a monitoring area of , with a total motion duration of . For the following two scenarios, we set the sampling period and the mean clutter rate . The Poisson rate of PPP birth is 0.01, the clutter is distributed in the region. For the birth density, the velocity follows a Gaussian distribution with zero mean and covariance 100, where represents an identity matrix of size 2. For the GGIW implementations, the position is Gaussian distributed with zero mean and covariance 150. In the two experimental scenarios, different numbers of target trajectories are set to test the algorithm’s performance. Assume a detection probability of , a survival probability of , and a forgetting factor for the measurement rate state set to . The experimental results indicate that using 1000 particles in our filter achieves a good balance between the runtime and estimation performance.
5.2. Simulation 1
For Simulation 1, we set up a challenging experiment to demonstrate the tracking error of target trajectories and the precision of extended state tracking under a fixed clutter rate, 10, by four filters. The parameters of target trajectories are shown in
Table 2. Note that the targets in three different motion states appear at different locations within the detection area at the 10th s, 20th s, and 30th s. With a total of 100 steps run, each trajectory initialized from the starting point to the midway step, with the remaining steps generated forwards and backwards.
In this scenario, we concentrate on estimating trajectories and compare the estimated extended states against the actual values. Among the various performance metrics for trajectory estimation, we employ the Linear Programming (LP) and the Gaussian Wasserstein Distance (GWD) metric
[
27]. We set the cut-off distance to 20, establish the order
, and set the track switch cost to 2 in order to apply these metrics at each time step.
The target trajectories are shown in
Figure 2. We know that the available measurement information contains both real targets and random clutter. From this figure, we observe three trajectories. Clearly, the three targets have nonlinear motion states, performing CV and CT motions in stages within the detection area and the model has different angular velocities during CT motion. During three stages, the respective motion states are designated as CV, CT, and CV in turn, where the turning velocity is given before the scenario. The clutter obeys Poisson distribution, and the measurement noise and covariance are unknown to the filter.
In
Figure 3a, it is evident that the filters proposed in this paper offer excellent extended state estimation and stable tracking performance for elliptical extended targets. In
Figure 3b, we note that the cardinality statistics change as time progresses. Obviously, all algorithms perform well in the potential estimation. The TPMB-BP, TPMB-GBP, IMM-TPMB-GBP-CKF algorithms, and the algorithm discussed in this paper all closely approximate the true number of targets, while the classical GGIW-TPMB algorithm underestimates the number of highly maneuverable targets; significant deviations in the potential estimation by the GGIW-TPMB algorithm from the true values are evident at the 30th s, 40th s, and 65th s. The TPMB-BP and TPMB-GBP algorithms have a faster computation speed compared to IMM-TPMB-BP-CKF, and both algorithms essentially provide similar estimates of target potential. Notably, the standard filters are affected by multiple measurements from random clutter, resulting in some positional estimation biases. The filters proposed in this paper effectively eliminate the impact of unknown noise on tracking performance, providing accurate estimates.
The track error of the trajectory in Scenario 1 under the LP metric is presented in
Table 3. For the extended target, the cost of state error dominates the overall LP estimation error. When the targets are in the same clutter, the measurement errors increase. For the most advanced TPMB-BP filter, we improved its correlation process and obtained the Gate-control BP association method that reduces the LP measurement error with adaptive gating. Among the five filters, the IMM-TPMB-GBP-ASCKF filter, with its superior noise adaptability, achieves the lowest LP metric. All five filters perform well when tracking targets at adequate intervals. However, as the targets move closer to each other, they become indistinguishable. Thanks to its adaptive noise estimation, the improved ASCKF enables the algorithm proposed in this study to achieve the lowest target state cost and false detection cost among the four filters. In such scenarios, the cost of losing targets predominates in the LP estimation error. Moreover, in situations where highly maneuvering targets frequently change their motion states, the improved algorithm allows the trajectory to achieve reduced switch cost during transitions.
Undoubtedly, the TPMB-BP filter is the current optimal filter, which avoids explicit assumption enumeration in data association compared to traditional filters, improves data association accuracy, and reduces the running time of the data association process. Our work is to further reduce the running time of the data association process and improve the anti-interference ability against noise based on the TPMB-BP filter, and its performance is demonstrated in Simulation 2.
5.3. Simulation 2
For Scenario 2, to further validate the performance of the proposed algorithm, we designed a set of comparative experiments under the low clutter rate of 10 and the high clutter rate of 30, respectively. The parameters of target trajectories are shown in
Table 4.
Figure 4 shows that all five filters can converge to the actual number of targets. It was observed that as the clutter density increases, the state estimation error of each algorithm also increases. The reason for the increase in state estimation error in standard filters is due to the generation of random clutter by the filter and environment while tracking maneuvering targets. In contrast, the other four filters exhibit better robustness to the state transitions of targets, resulting in more accurate state estimations. More importantly, this paper employs a method using interactive gating with GBP data association and an ASCKF that adapts the covariance matrix for unknown noise, addressing inaccuracies in state estimation under high clutter and unknown noise conditions. However, the cardinality estimated by standard filters is the least stable. During low clutter phases, false alarms occur due to the filter’s inability to adaptively model estimates for non-linear targets. Since the targets maintain non-maneuverability, standard filters mix random clutter and unknown noise, leading to inaccurate state estimation. Conversely, at the 55th s mark, there is a missed detection because two of the targets switch from CT to CV motion models, and the filters fail to immediately adapt to the motion changes. The IMM-TPMB-GBP-CKF, TPMB-GBP, and TPMB-BP algorithms exhibit lower state estimation errors and show stable performance throughout the detection period, though they also face the issues of increased estimation error with higher clutter rates. In contrast, the algorithm proposed in this work maintains low state estimation errors under both low and high clutter rates. It is less affected by increases in clutter and motion model transitions and controls the overestimation of target states under high clutter conditions.
The GOSPA distances for the five filters are shown in
Figure 5. As expected, throughout the detection period, the algorithm proposed in this paper consistently shows the lowest GOSPA distance. At the 10th s, 20th s, 45th s, and 65th s, the peaks in clutter density from the standard filter are caused by false alarms, mode transitions, and missed targets, involving both cardinality and location errors. The experiments reveal that due to the inability of the filters to switch motion models during the target’s motion cycle, significant tracking mismatches occur, leading to increased state estimation errors. Since the targets include CV and CT models, which appear with equal probability at certain times, the proposed filter automatically adjusts the gating thresholds based on the target state model, reducing the number of incorrect associations and association duration, thereby enhancing the accuracy of target measurement associations and correcting unstable cardinality estimates. Considering the actual tracking reliability of the filter proposed in this paper under different clutter rates, four types of filters were operated under two different clutter rates. This figure displays a comparison chart of the GOSPA errors of the five filters under these two clutter rates. Notably, as the clutter rate increases, the filter proposed in this paper provides stable cardinality estimates. In situations with noise and unknown covariance, the algorithm proposed demonstrates stable estimation performance. Moreover, in the ASCKF, it also performs well with noise having smaller covariance. This means that the filter proposed in this paper does not exhibit significant changes under various clutter environments.
Finally, due to the specificity of the GOSPA evaluation algorithm, it can be differentiated into two parts: the number of missed targets and the number of false targets. The missed and false target state estimates are compared under different clutter rates.
Figure 6 shows the missed target state estimation under different clutter rates. It can be seen that all four filters experience a sharp increase in the number of missed targets at the 10th s, as targets begin to appear, due to the classical filter’s lack of a multimodal assumption and the increasing number of targets under the influence of clutter. However, the remaining three filters, having a multimodal assumption, adjust their matching models with changes in the motion state of the models, thus significantly reducing missed targets compared to the classical filter. Compared with TPMB-GBP and TPMB-BP, the improved TPMB-GBP reduced the GOSPA error. Due to the lack of a mixture of adaptive multi-model algorithm for high maneuvering targets, the performance difference between the two algorithms is relatively small.
Table 5 shows that TBP-GBP has shorter running time under different clutter and particles. The IMM-TPMB-GBP-CKF, TPMB-GBP, and TPMB-BP algorithms show an increase in missed counts under different clutter levels, but the IMM-TPMB-GBP-CKF algorithm reduces the number of misses in high clutter rates due to adaptive gating adjustments and also reduces operation time. The algorithm proposed in this paper, with the inclusion of ASCKF, reduces target location errors and clutter interference, and the number of missed targets does not significantly change with increasing clutter.
The false target state estimation, shown in
Figure 7a, shows a significant increase at the 10th s as the number of emerging targets grows. It can be observed that at the 20th s, as the target motion models switch, the number of false targets reaches its peak, with the algorithm proposed in this paper having the lowest false target state estimation among the five filters. Studies show that, compared to the four filters analyzed alongside the proposed algorithm, this algorithm exhibits the lowest GOSPA error values and adapts to non-linear target estimation under the influence of different clutter rates.
In order to verify the average running time of the proposed filter in this article under different particle and clutter conditions, we set up a set of comparative experiments to verify its effectiveness. Due to the fact that the GGIW-TPMB filter is a Gaussian implementation method, we only compare the average running time of the other four filters under different particles. Note that TPMB-BP is a particle implementation based on the confidence propagation data association, which results in a large computational workload due to multiple iterations of calculating example weights. On the basis of the same tracking effect, the proposed algorithm reduces the computation of particle likelihood in the target measurement matching process by integrating the TPMB-GBP algorithm for comparison, thereby reducing the overall time. With the increase in trajectory, iteration count, and clutter rate, the time consumption of the IMM-TPMB-GBP-ASCKF filter is only slightly higher than that of the IMM-TPMB-GBP-CKF filter, which is reasonable due to the additional noise estimation process of the proposed filter. As the noise increases, the running time of the four filters significantly increases. Due to the improved ASCKF filter’s good adaptability to noise, it can effectively filter out clutter, and the time consumption does not change significantly compared to CKF. Additionally, under high clutter density conditions, the IMM-TPMB-GBP-ASCKF filter demonstrates higher accuracy. The results indicate that the proposed IMM-TPMB-GBP-ASCKF filter maintains robust tracking performance for highly maneuverable targets and under the impact of unknown noise, while incurring a moderate computational cost.
6. Conclusions
This paper introduces a new robust TPMB filter designed for conditions with uncertain process and measurement noise. By enhancing the adaptive multi-model switching method and performing process noise estimation, we provide robust recursive formulations for TPMB. In terms of model-driven aspects, we have improved the IMM by increasing the matching degree of target motion mode switching through model probability transformation at two adjacent moments. From a data-driven perspective, the GBP in the factor graph for the association of measurement vectors allows the model to adaptively adjust the gating threshold based on the association vector, effectively reducing the computational complexity of data association. Additionally, we present an improved method of importance sampling for IMM-TPMB’s ASCKF particles, optimizing the particle degradation issue by estimating the process and measurement noise in the context of target positioning. Simulation studies show that, compared to the GGIW-TPMB, TPMB-BP, TPMB-GBP, and IMM-TPMB-GBP-CKF, the IMM-TPMB-GBP-ASCKF achieves the best overall performance in terms of filtering efficiency and computation time across different target numbers, motion states, and clutter rates.
Future work includes conducting more extensive tests on different clutter models. Although the algorithm proposed in this paper mitigates particle degradation and running time to some extent, the running time continues to increase with the number of targets, and measurement estimation errors arise due to closely spaced measurements. For scenarios with a high number of measurements and closely spaced measurements, an improved measurement tracking gating method can be employed for processing.