A Robust TCPHD Filter for Multi-Sensor Multitarget Tracking Based on a Gaussian–Student’s t-Mixture Model

Abstract

To realize multitarget trajectory tracking under non-Gaussian heavy-tailed noise, we propose a Gaussian–Student t-mixture distribution-based trajectory cardinality probability hypothesis density filter (GSTM-TCPHD). We introduce the multi-sensor GSTM-TCPHD (MS-GSTM-TCPHD) filter to enhance tracking performance. Conventional cardinality probability hypothesis density (CPHD) filters typically assume Gaussian noise and struggle to accurately establish target trajectories when faced with heavy-tailed non-Gaussian distributions. Heavy-tailed noise leads to significant estimation errors and filter dispersion. Moreover, the exact trajectory of the target is crucial for tracking and prediction. Our proposed GSTM-TCPHD filter utilizes the GSTM distribution to model heavy-tailed noise, reducing modeling errors and generating a set of potential target trajectories. Since single sensors have a limited field of view and limited measurement information, we extend the filter to a multi-sensor scenario. To tackle the issue of data explosion from multiple sensors, we employed a greedy approximation method to assess measurements and introduced the MS-GSTM-TCPHD filter. The simulation results demonstrate that our proposed filter outperforms the CPHD/TCPHD filter and Student’s t-based TCPHD filter in terms of accurately estimating the trajectories of multiple targets during tracking while also achieving improved accuracy and shorter processing time.

1. Introduction

1.1. Literature Review and Motivation

Target tracking is widely used in intelligent driving, robot identification, and defense. The observation platform employs an assortment of sensors, including millimeter-wave radar, lidar, and camera, to acquire data; applies filters to preserve pertinent information; and ultimately estimates the target’s state information [1].

Conventional techniques for multitarget tracking rely on data association; these include nearest-neighbor association, multi-hypothesis tracking (MHT), and joint probabilistic data association [2]. However, these algorithms are susceptible to combinatorial explosion because of the data association concept. Random finite set (RFS) theory can simultaneously reduce correlation operations and has a solid theoretical demonstration inside the Bayesian framework. RFS filters commonly used today include probabilistic hypothesis density (PHD), cardinality PHD (CPHD) [3], multi-objective multi-Bernoulli (MeMBer), and generalized labeled multi-Bernoulli (GLMB) [4]. Given the GLMB filter’s numerous prediction and update processes, it demands a lot of processing power. The PHD filter uses fewer computing resources [5]. In recent years, trajectory cardinality probability hypothesis density (TCPHD) filters have emerged to address this shortcoming by enabling trajectory generation based on CPHD [6], which provides trajectory representation with superior estimation accuracy and computational complexity. The TCPHD filter addresses the issue of the standard CPHD filter’s inability to infer trajectories by utilizing the set of trajectories in the state variables instead of the set of targets. It employs IID-clustered multi-trajectory densities for the prediction and update steps, followed by KLD minimization after the update step. Compared with previous trajectory modeling methods used in CPHD filters, the TCPHD filter has a better theoretical foundation and trajectory estimation performance. It is also capable of updating the past state information of the trajectory [7]. The performance of filters dominated by Gaussian noise may rapidly deteriorate in the presence of outliers and non-Gaussian noise. This leads to poorer estimates due to model mismatches [8]. Therefore, the mainstream trend is to replace Gaussian models with robust estimators.

The Kalman filter is a technique that works under the assumption of Gaussian noise and tries to minimize root-mean-square errors. However, in practice, non-Gaussian noise is more common and frequently exhibits heavy-tailed characteristics that depart from the Gaussian distribution [9]. Non-Gaussian heavy-tailed noise is caused by the sensor’s performance becoming unstable under intense target motion [10]. When the target moves rapidly, the state displays heavy-tailed non-Gaussian noise [11]; changes in the detected signal’s reflection cause non-Gaussian heavy-tailed noise in the measurements [12]. Non-Gaussian heavy-tailed noise causes significant reductions in tracking performance and an increase in model mismatch, filter dispersion, and prior filter noise covariance modeling errors.

Research on target tracking with non-Gaussian noise can be divided into two categories: when the measurement noise is non-Gaussian and when both the process and measurement noise are non-Gaussian [13]. There are three mainstream methods for dealing with process noise, regardless of whether it is Gaussian or not: variational Bayes (VB) inference [10], the sequential Monte Carlo method [14], and Student’s t-distribution Kalman filtering [15]. The VB inference method is an adaptive Kalman filtering algorithm that approximates the target state via VB inference by introducing auxiliary noise parameters. For instance, the covariance of the observed noise is subjected to an inverse gamma distribution or an inverse Weiser distribution, and the variational Bayesian method is used to jointly estimate the parameters of the target. The sequential Monte Carlo method utilizes particle filtering to estimate measurement noise outliers. It employs weighted particle approximation to represent the posterior probability of the target. However, it requires a sufficient number of particles, resulting in a significant increase in computational complexity and low efficiency in high-dimensional space. The Kalman filtering method used by Student’s t-distribution model incorporates noise with heavy tail characteristics. Various nonlinear Kalman-filtering algorithms based on Student’s t-distribution have been proposed, depending on the integration operation [16]. Although Student’s t-distribution provides good tracking performance, it requires fixed-point iteration to calculate various parameters in variational Bayesian inference, which increases computational complexity. Other processing methods, such as Huber Kalman filtering and maximum correlation entropy criterion Kalman filtering, are also used by a small number of researchers [17,18,19,20]. In recent years, the GSTM distribution has maintained its characteristic of modeling outlier noise using Student’s t-distribution while also utilizing the Gaussian distribution to improve computational efficiency [20]. Further research is required to investigate the application of the GSTM distribution in target tracking with non-Gaussian noise.

Currently, there is limited research on outlier noise in RFS multitarget tracking. Two types of approaches are commonly used: Student’s t-distribution PHD filter and variational Bayesian inference. For the PHD filter method based on Student’s t-distribution [21], Liu Z proposed a PHD filter algorithm, modeling heavy-tailed noise with Student’s t-distribution, simultaneously expressing the same posterior probability density distribution as Student’s t-distribution alone. Later, the method was extended to nonlinear scenarios [22]. For the multitarget tracking method based on variational Bayesian approximation, Yang J and Li W proposed using inverse gamma distribution to represent the variance in measurement noise and expressed the posterior probability distribution as inverse gamma distribution, used the VB method to estimate parameters, and proposed the implementation of PHD [23]. Yang J L also proposed a CBMeMBer filter based on VB approximation [24]. Li W proposed modeling flicker noise with Student’s t-distribution and provided a mixed implementation of a Gaussian gamma distribution of the PHD filter [25]. The use of PHD filters in multitarget tracking has inherent drawbacks, such as the potential for large estimation errors and the inability to form trajectories. However, the TCPHD filter, which is based on CPHD, compensates for these drawbacks and has the added advantage of forming trajectories. This makes it a useful tool for multitarget tracking, particularly in the presence of outlier noise.

Tracking based on single sensors has the problems of limited field of view, measurement information, and instability, so the data fusion technology used in a multi-sensor is crucial [26]. Mahler first proposed the iterated corrector PHD (IC-PHD) filter, which sequentially processes information from different sensors [27]. In 2009, Mahler introduced the concept of the General PHD filter (G-PHD) for multi-sensor systems, and in 2015, Nannuru derived a rigorous theoretical proof and proposed a General CPHD filter [28]. In 2010, Mahler proposed the product multi-sensor PHD (PM-PHD) filter and the product multi-sensor CPHD (PM-CPHD) filter [29].

In the algorithm for multi-sensor multitarget tracking based on random finite sets, measurement partitioning is a crucial step. The algorithm for greedy approximation obtains the optimal measurement set for each target in sensor order while retaining the best measurement combination. Compared with Gibbs sampling, which assumes a finite number of measurements with higher weights as the final measurement partition result, the greedy partitioning method aims to preserve the measurement set that contains missed items when obtaining the measurement set for each target [30]. This ensures that the measurement partitioning results are in line with reality and are not affected by small likelihood weights, resulting in better tracking results. The greedy approximation algorithm is significant in multi-sensor and multitarget tracking scenarios where there is low probability and high maneuverability.

Based on the above research, multi-sensor multitarget tracking under heavy-tailed noise has important research significance and various application scenarios. The three key challenges include (1) generating the trajectory of many targets; (2) modeling non-Gaussian heavy-tailed noise to engage in iterative filter tracking [31]; and (3) leveraging the measurement data from multiple sensors to enhance tracking performance.

1.2. Our Contributions

The multi-sensor multitarget tracking problem under heavy-tailed noise has significant scientific value and promising application opportunities, according to the paper mentioned above. Consequently, we proposed the GSTM-TCPHD and MS-GSTM-TCPHD filters. The following are some of our approach’s primary innovations:

In addition to introducing some related multitarget tracking contents based on a random finite set, this study proposes a new multitarget trajectory tracking method under non-Gaussian noise, which uses GSTM to model non-Gaussian noise and implements multitarget tracking based on the TCPHD filter. Simultaneously, a new multi-sensor multitarget tracking method is introduced based on the GSTM-TCPHD filter. This method generates trajectories under non-Gaussian noise, and the implementation of the MS-GSTM-TCPHD filter is provided. Then, the proposed filter is simulated and evaluated to verify its tracking performance.

1.3. Organization of the Article

This article is structured as follows: Section 2 provides background knowledge. Section 3 presents a theoretical derivation process for the GSTM-TCPHD filter. Section 4 shows the implementation of the filter. Section 5 introduces the implementation of the MS-GSTM-TCPHD filter. In Section 6, the proposed method is validated, and the results of the simulation run are demonstrated and analyzed. Section 7 is the conclusion of this article.

2. Background

2.1. Notation

Assume that the variable $X=\left(t,x^{1:l}\right)$ with $x\in R^{n_x}$ is represented as a single trajectory variable, where

$t$ denotes the starting time;

$n_x$ is the dimension of the target state, $x$;

$x^{1:l}=\left(x^1,\dots ,x^l\right)$ denotes a set of target state sequences of length $l$.

At time k, the length of the trajectory is $t+l-1$; let $I_k=\{(t,i):0\le t\le k,1\le i\le k-t+1\}$. The trajectory space is represented as

$${\mathbb{T}}_k={\uplus }_{(t,i)\in I_k}\left\{t\right\}\times R^{i{n}_x}$$

(1)

where $\uplus$ denotes the mutually exclusive set.

The binomial coefficients, permutation coefficients, and generalized Kronecker delta function are expressed as follows:

$$C_m^n=\frac{n!}{m!(n-m)!},\hspace{1em}{P}_j^n=\frac{n!}{(n-m)!},\hspace{1em}{\delta }_A[B]\triangleq \left\{\begin{array}{cc}1,\hfill & \mathrm{if} A=B\hfill \\ 0,\hfill & \mathrm{otherwise} \hfill \end{array}\right.$$

(2)

The inner product of two real-valued functions is expressed as follows:

$$〈a,b〉={\displaystyle \int a\left(x\right)b\left(x\right)dx}$$

(3)

A j-th order primitive symmetric function of a finite set, $\mathcal{X}$, of real numbers is expressed as follows:

$$e_j(\mathcal{X})={\displaystyle \sum _{\sigma \subseteq \mathcal{X},\left|\sigma \right|=j}{\displaystyle \prod _{\xi \in \sigma }\xi }}$$

(4)

If p is the observed probability distribution function and q is used to approximate it [32], then the Kullback–Leibler divergence (KLD) between the two distributions can be calculated as follows:

$$D_{KL}(p||q)={\displaystyle \sum _{i=1}^{N}p}\left(x_i\right)\cdot \left(\log p\left(x_i\right)-\log q\left(x_i\right)\right)={\displaystyle \sum _{i=1}^{N}p}\left(x_i\right)\cdot \log \frac{p\left(x_i\right)}{q\left(x_i\right)}$$

(5)

2.2. Set Integral

Suppose that, in the trajectory space, the single trajectory variable, $\pi \left(X\right)$, and the trajectory RFS have the following representation:

$${\displaystyle \underset{\mathbb{T}}{\int }}\pi \left(X\right)dX={\displaystyle \sum _{(t,i)\in I_{(k)}}{\displaystyle \int \pi }}\left(t,x^{1:i}\right)d{x}^{1:i}$$

(6)

$${\displaystyle \int \pi }(X)\delta X={\displaystyle \sum _{n=0}^{\infty }\frac{1}{n!}}{\displaystyle \int \pi }\left(\left\{X_1,\dots ,X_n\right\}\right)d{X}_{1:n}$$

(7)

where

$X_k=\left\{X_1,\dots ,X_{N^k}\right\}\in \mathcal{F}\left({\mathbb{T}}_k\right)$;

$\mathcal{F}\left({\mathbb{T}}_k\right)$ denotes the set of all finite subsets in the trajectory space.

The cardinality distribution and PHD of the multi-trajectory density, $\pi (X)$, are expressed as

$${\rho }_{\pi }(n)=\frac{1}{n!}{\displaystyle \int \pi }\left(\left\{X_1,\dots ,X_n\right\}\right)d{X}_{1:n}$$

(8)

$$D_{\pi }(X)={\displaystyle \int \pi }(\{X\}\cup X)\delta X$$

(9)

Assuming that there is a Poisson multi-trajectory density, $\pi (X)$, at time $k$, the target’s multitarget density, ${\pi }_k^{\tau }$, is Poisson, and its PHD can be expressed as follows:

$$D_{{\pi }_k^{\tau }}(y)={\displaystyle \sum _{t=1}^k{\displaystyle \int D_{\pi }}}\left(t,x^{1:k-t},y\right)d{x}^{1:k-t}$$

(10)

2.3. TPHD and TCPHD Filter

The TPHD filter propagates a Poisson multi-trajectory density, its multi-trajectory density is expressed as follows:

$$p_k\left(\left\{X_1,\dots ,X_{N^k}\right\}\right)=e^{-{\lambda }_k}{\lambda }_k^n{\displaystyle \prod _{j=1}^{N^k}{\dot{p}}_k}\left(X_j\right)$$

(11)

where ${\lambda }_k$ represents the expected number of trajectories.

The TCPHD filter propagates IID cluster Poisson multi-trajectory density; its multi-trajectory density is expressed as follows:

$$p_k\left(\left\{X_1,\dots ,X_{N^k}\right\}\right)={\rho }_k(n)n!{\displaystyle \prod _{j=1}^{N^k}{\dot{p}}_k}\left(X_j\right)$$

(12)

where

$\dot{p}\left(X_j\right)$ is the density of single trajectories;

$\rho (n)$ is the cardinality distribution.

When the number of sets is distributed according to a Poisson distribution, the IID cluster random finite set degenerates into a Poisson random finite set. The PHD of the multi-trajectory densities of the Poisson RFS and the IID cluster RFS are

$$D_k(X)={\lambda }_k{\dot{p}}_k(X) , D_k(X)={\dot{p}}_k(X){\displaystyle \sum _{n=0}^{\infty }n}{\rho }_k(n)$$

(13)

At time $k-1$, the density of multiple trajectories in the state space, $X_{k-1}=\left\{X_1,\dots ,X_{N^{k-1}}\right\}\in \mathcal{F}\left({\mathbb{T}}_{k-1}\right)$, and the density of multiple trajectories in the measurement space, $Z_{k-1}=\left\{z_1,\dots ,z_{M^{k-1}}\right\}\in \mathcal{F}\left({\mathbb{Z}}_{k-1}\right)$, result in a Bayesian recursive representation as follows [3]

$${\pi }_{k|k-1}\left(X_k\right)={\displaystyle \int f}\left(X_k|X_{k-1}\right){\pi }_{k-1}\left(X_{k-1}\right)\delta X_{k-1}$$

(14)

$${\pi }_k\left(X_k\right)=\frac{{\mathcal{l}}_k\left(Z_k|X_k\right){\pi }_{k|k-1}\left(X_k\right)}{{\displaystyle \int {\mathcal{l}}_k}\left(Z_k|X_k\right){\pi }_{k|k-1}\left(X_k\right)\delta X_k}$$

(15)

where

$f\left(X_k|X_{k-1}\right)$ denotes the transfer density of the trajectory;

${\pi }_{k|k-1}\left(X_k\right)$ denotes the predicted multi-trajectory density at time $k$;

${\mathcal{l}}_k\left(Z_k|X_k\right)$ denotes the quantile likelihood function.

2.4. Multi-Sensor Observation Model

At time $k$, assuming the existence of $s$ sensors, the set of measurements is denoted as $Z_k^{1:s}=\left\{Z_k^1,\dots ,Z_k^s\right\}$. Let $W\subseteq Z_{k+1}^{1:s},|W{|}_j\le 1,|W{|}_j=\left|\left\{z\in W:z\in Z_{k+1}^j\right\}\right|$, denoting that the subset W takes, at most, one measurement from the set of measurements of each sensor, and all belong to the same PHD component; $\mathsf{\Omega }$ denotes the set of all measure subsets of W. For a mutually exclusive subset of W, let $V=Z_{k+1}^{1:s}\backslash \left({\cup }_{i=1}^n{W}_i\right)$ represents the set of clutters generated by all sensor measurements [33]. The set of measurements can be divided into

$$P=\left\{W_1,W_2,\dots W_{\left|P\right|-1},V\right\}$$

(16)

$${\displaystyle \underset{i=1}{\overset{\left|P\right|-1}{\cup }}{W}_i}\cup V=Z_{k+1}^{1:s},W_i\cap W_j=\varnothing ,W_i\cap V=\varnothing ,W_i,W_j\in P,i\ne j$$

(17)

where

P denotes the set of divisions;

$\left|P\right|$ denotes the number of partitions, P;

$|P{|}_j$ is the measurement produced by the sensor, $j$, on the target and $|P{|}_j={\displaystyle \sum _{i=1}^{\left|P\right|-1}{\left|W_i\right|}_j}$.

The number of measurements generated in sensor j that are considered to be clutter is $m_j-|P{|}_j$, and $m_j$ is the number of measurements generated in the j-th sensor.

3. Proposed GSTM-TCPHD Filters

Compared with the CPHD filter, the TCPHD filter has a wider range of scenarios given its ability to track multiple targets and produce target trajectories; however, under non-Gaussian heavy-tailed noise, it lacks relevant theory [34]. In this section, we introduce the GSTM-TCPHD filter and provide a theoretical derivation based on the L-scan TPHD and TCPHD filters [35].

3.1. GSTM Distribution

While non-Gaussian models are occasionally seen, Gaussian models account for the majority of noise [36]. Heavy tails are a defining feature of Student’s t-distribution, better simulating the heavy tail phenomenon [37]. Assume that the random variable, $x$, follows Student’s t-distribution with a mean, $\mu$; a scale matrix, $\Sigma$; and degrees of freedom, $v$ [38], and denote its probability density function as $St\left(x;\mu ,\Sigma ,v\right)$. The GSTM distribution adaptively learns the mixture of Student’s t-distribution and Gaussian probabilities to adjust to time-varying state and measurement noise.

Assume that the GSTM distribution with mixing probability, $e$, is represented as follows:

$$p(x)=eN(x;\mu ,\Sigma )+(1-e)\mathrm{St}(x;\mu ,\Sigma ,\nu )$$

(18)

In this paper,$Gst\left(x;\mu ,\Sigma ,v,e\right)$ is used to describe the Gaussian–Student’s t-mixture distribution.

The heavy-tailed characteristic of Student’s t-distribution is modulated by a degrees-of-freedom parameter, and the thickness on either side of the probability curve is correlated with the degrees of freedom. As the degree of freedom tends toward a Gaussian distribution, the more critical the heavy tail, and the smaller the heavy tail, the larger the degree of freedom [39]. The GSTM distribution uses the mixture probability parameter to adaptively adjust the probability of Student’s t-distribution and the Gaussian distribution to match the true noise distribution, as shown in Figure 1. Figure 1 shows that there are differences between the three distributions near the center point. As the distance range increases, the difference between the three distributions narrows and then increases. This suggests that the GSTM distribution may perform better than the Gaussian distribution and Student’s t-distribution, regardless of whether the situation is normal or abnormal.

3.2. Implementation of GSTM-TCPHD Filter

At time $k-1$, the multitarget posterior probability hypothesis density,$D_{k-1}(x)$, and the posterior potential distribution, ${\rho }_{k-1}$, are known. The GSTM-TCPHD filter is predicted and updated as follows [40].

We make the following assumptions in the prediction step:

P1: Each target survives with probability $p_S\left(x\right)=p_S$; moves to a new state with a transition density of $g\left(·\right)$; or dies with probability $1-p_S$.

P2: The multitarget state at the next time step is the union of the surviving targets and new targets, which are born independently with a Poisson multitarget density, $\beta \left(·\right)$.

P3: The state transition is represented as $g\left(x^k|x^{k-1}\right)=Gst\left(x_k;F{x}_{k-1},Q,v_1,\tau \right)$.

At time $k-1$, multi-trajectory PHD and cardinality distribution are as follows:

$$D_{k|k-1}(X)=p_S\left(x^{i-1}\right)f\left(x^i|x^{i-1}\right)D_{k-1}\left(t,x^{1:i-1}\right)+{\beta }_k(X)$$

(19)

$${\rho }_{k|k-1}(n)={\displaystyle \sum _{j=0}^n{\rho }_{\beta ,k}}(n-j){\displaystyle \sum _{l=j}^{\infty }{C}_l^j}{\rho }_{k-1}(l)\frac{{〈p_S,D_{k-1}^{\tau }〉}^j{〈1-p_S,D_{k-1}^{\tau }〉}^{l-j}}{{〈1,D_{k-1}^{\tau }〉}^l}$$

(20)

where

$$D_{k|k-1}^{\tau }\left(x^i\right)={\displaystyle \sum _{t=1}^k{\displaystyle \int D_{k|k-1}}}\left(t,x^{1:k-t+1}\right)d{x}^{1:k-t}$$

(21)

Assuming that the predicted TPHD, $D_{k|k-1}(x)$, and cardinality distribution, $p_{k|k-1}$, are known, the following assumptions are made in the update step [41].

U1: The probability of detection for each target is $p_D$, and each detection produces a measurement with density $l\left(·|x\right)$, or the probability of missing a detection is $1-p_D\left(x\right)$.

U2: Measurement $z_k$ is the union of the target-generated measurements with the IID cluster density, $c\left(·\right)$.

U3: The quantile likelihood function is expressed as $l(z|x)=Gst\left(z;Hx,R,v_2,\pi \right)$.

U4: The multi-trajectory density, ${\omega }_k\left(·\right)$, is the IID cluster, RFS.

Given the predicted cardinality distribution, ${\rho }_{k|k-1}(n)$, and TPHD, $D_{k-1}\left(t,x^{1:i-1}\right)$, at moment $k$, the update step of the TCPHD filter is as follows:

$${\rho }_k(n)=\frac{{\mathrm{Y}}_k^0\left[D_{k|k-1}^{\tau };z_k\right](n){\rho }_{k|k-1}(n)}{〈{\mathrm{Y}}_k^0\left[D_{k|k-1}^{\tau };z_k\right],{\rho }_{k|k-1}〉},$$

(22)

$$\begin{array}{l}{D}_k(X)=D_{k|k-1}(X)q_D\left(x^i\right)\times \frac{〈{\mathrm{Y}}_k^1\left[D_{k|k-1}^{\tau };Z_k\right],{\rho }_{k|k-1}〉}{〈{\mathrm{Y}}_k^0\left[D_{k|k-1}^{\tau };Z_k\right],{\rho }_{k|k-1}〉}\\ \hspace{1em}\hspace{1em}\hspace{1em}+D_{k|k-1}(X)p_D\left(x^i\right)\times {\displaystyle \sum _{z\in Z_k}\frac{l_k\left(z|x^i\right)}{\overline{c}(z)}}\frac{〈{\mathrm{Y}}_k^1\left[D_{k|k-1}^{\tau };Z_k\backslash \left\{z\right\}\right],{\rho }_{k|k-1}〉}{〈{\mathrm{Y}}_k^0\left[D_{k|k-1}^{\tau };Z_k\right],{\rho }_{k|k-1}〉}\end{array}$$

(23)

where

$q_{D,k}$ denotes the probability missing;

$\overline{c}(z)$ is the clutter density.

$$\begin{array}{l}{\mathrm{Y}}^u\left[D_{k|k-1}^{\tau },Z_k\right](n)={\displaystyle \sum _{j=0}^{\min \left(M^k,n-u\right)}\left(M^k-j\right)}!{\rho }_c\left(M^k-j\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \hspace{1em} P_{j+u}^n\frac{{〈q_D,D_{k|k-1}^{\tau }〉}^{n-j-u}}{{〈1,D_{k|k-1}^{\tau }〉}^n}{e}_j\left(\mathrm{\Xi }\left(D_{k|k-1}^{\tau },Z_k\right)\right)\end{array}$$

(24)

$$\mathrm{\Xi }\left(D_{k|k-1}^{\tau },Z_k\right)=\left\{{\displaystyle \int p_{D,k}}\left(x^i\right)\frac{l_k\left(z|x^i\right)}{\overline{c}(z)}{D}_{k|k-1}^{\tau }\left(x^i\right)d{x}^i:z\in Z_k\right\}$$

(25)

3.3. VB Approximation for Posterior Probabilities

The GSTM state model’s parameters are more difficult to calculate given their interdependence. By breaking down the multivariate variables into products of individual variables, VB approximation is used to approximate the posterior probability density and provide an analytical solution to the posterior probability. Figure 2 displays the filter’s flowchart.

The VB method uses an approximation of simply computed joint density, $q\left(X,\theta \right)$, to approximate the true joint a posteriori probability density, $q\left(X,\theta |Z\right)$, while considering the parameters to be independent of each other, computed as $q\left(X,\theta \right)=q\left(X\right)q\left(\theta \right)$. The errors on both sides of the equation are judged by the KL distance.

VB approximation of the solution equation for the posterior probability is expressed as follows:

$$\log q(\phi )={\mathrm{E}}_{{\theta }^{\left(-\phi \right)}}[\log p(X,\theta ,Z)]+c_{\phi },\hspace{1em}\mathrm{s}.\mathrm{t}. {\displaystyle \int q}(\phi )d\phi =1$$

(26)

where

$\theta$ denotes the set consisting of state vectors and unknown parameters;

$\phi \in \theta$ is denoted as any element of set $\theta$;

${\theta }^{\left(-\phi \right)}$ denotes the elements remaining in set $\theta$ after the removal of element $\phi$;

$q\left(\phi \right)$ denotes the approximate a posteriori density of $\phi$;

$E_x$ denotes the expectation of variable $x$;

$c_{\phi }$ represents a constant independent of $\phi$.

The GSTM distribution models the noise with the parameter $\left\{{\xi }_k,t_k,{\tau }_k\right\}$ describing the state noise and the parameter $\left\{{\lambda }_k,y_k,{\pi }_k\right\}$ describing the measurement noise. Its joint probability density is as follows:

$$p(X,\theta ,Z)=p(Z)g\left(x^k|x^{k-1}\right)l(z|x)$$

(27)

The VB approximation to the parameter decomposition is expressed as follows:

$$q\left(X,\theta |Z\right)=q\left(x\right)q\left(\xi \right)q\left(t\right)q\left(\tau \right)q\left(\lambda \right)q\left(y\right)q\left(\pi \right)$$

(28)

4. GSTM-TCPHD Filter Implementation under Single Sensor

The GSTM distribution is used to model the noise. In cases where the state noise and measurement noise in the GM-TCPHD filter mixture are non-Gaussian heavy-tailed distributions, the VB approximation is solved to update the state to a Gaussian mixture [42].

4.1. Predict Process

At time $k-1$, the PHD representation of a priori multi-trajectory density, ${\pi }^{k-1}\left(·\right)$, is as follows

$$D_{{\pi }^{k-1}}(X)={\displaystyle \sum _{j=1}^{J^{k-1}}{w}_j^{k-1}}\mathcal{N}\left(X;t_j^{k-1},m_j^{k-1},P_j^{k-1}\right)$$

(29)

where

$t_j^{k-1}+i_j^{k-1}-1=k-1$;

$i_j^{k-1}=\dim \left(m_j^{k-1}\right)/n_x$.

The PHD of the birth density, ${\beta }^k$, is

$$D_{{\beta }^k}(X)={\displaystyle \sum _{j=1}^{J_{\beta }^k}{w}_{\beta ,j}^k}\mathcal{N}\left(X;k,m_{\beta ,j}^k,P_{\beta ,j}^k\right)$$

(30)

where

$J_{\beta }^k\in \mathrm{N}$ is the number of components;

$w_{\beta ,j}^k$ is the weight of the j-th component;

$m_{\beta ,j}^k,P_{\beta ,j}^k$ is the mean and covariance matrix.

The TPHD is predicted as follows:

$$D_{S,k|k-1}=D_{{\beta }^k}\left(X\right)+p_S{\displaystyle \sum _{j=1}^{J^{k-1}}{w}_j^{k-1}}N\left(x_{k|k-1};t_j^{k-1},m_{\omega ,j}^k,P_{\omega ,j}^k\right)$$

(31)

where

$${\dot{F}}_j=\left[0_{1,i_j^{k-1}-1},1\right]\otimes F$$

(32)

$$m_{\omega ,j}^k={\left[{\left(m_j^{k-1}\right)}^T,{\left({\dot{F}}_j{m}_j^{k-1}\right)}^T\right]}^T$$

(33)

$$P_{\omega ,j}^k=\left[\begin{array}{cc}{P}_j^{k-1}& P_j^{k-1}{\dot{F}}_j^T\\ {\dot{F}}_j{P}_j^{k-1}& {\dot{F}}_j{P}_j^{k-1}{\dot{F}}_j^T+Q\end{array}\right]$$

(34)

Cardinality is predicted as follows:

$${\rho }_{{\omega }^k}(m)={\displaystyle \sum _{j=0}^m{\rho }_{{\beta }^k}}(m-j)\times {\displaystyle \sum _{l=j}^{\infty }\left(\begin{array}{c}n\\ j\end{array}\right)}{\rho }_{{\pi }^{k-1}}(n){\left(1-p_S\right)}^{l-j}{p}_S^j$$

(35)

4.2. Update Process

Suppose that the TPHD representation of the a priori density, ${\omega }^k$, in the update step is as follows:

$$D_{{\omega }^k}(X)={\displaystyle \sum _{j=1}^{J_{\omega }^k}{w}_{\omega ,j}^k}\mathcal{N}\left(X;t_k^{\omega ,j},m_k^{\omega ,j},P_k^{\omega ,j}\right)$$

(36)

updated to

$$\begin{array}{l}{D}_{{\pi }^k}(X)=\frac{〈{\mathsf{\Psi }}^1\left[w_{\omega }^k,z^k\right],{\rho }_{{\omega }^k}〉}{〈{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right],{\rho }_{{\omega }^k}〉}\left(1-p_D\right)D_{{\omega }^k}(X)\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{z\in z^k}{\displaystyle \sum _{j=1}^{J_{\omega }^k}{w}_j}}(z)\mathcal{N}\left(X;t_{\omega ,j}^k,m_j^k(z),P_j^k\right)\end{array}$$

(37)

$${\rho }_{{\pi }^k}(n)=\frac{{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right](n){\rho }_{{\omega }^k}(n)}{〈{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right],{\rho }_{{\omega }^k}〉}$$

(38)

$$\begin{array}{l}{\mathsf{\Psi }}^u\left[D_{k|k-1}^{\tau },Z_k\right](n)={\displaystyle \sum _{j=0}^{\min \left(\left|z^k\right|,n-u\right)}\left(\left|z^k\right|-j\right)}!{\rho }_c\left(\left|z^k\right|-j\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \frac{{\left(1-p_D\right)}^{n-(j+u)}}{{〈1,w_{\omega }^k〉}^{j+u}}\frac{n!}{(n-j-u)!}\times e_j\left(\mathsf{\Lambda }\left(w_{\omega }^k,z^k\right)\right)\end{array}$$

(39)

$$\mathsf{\Lambda }\left(w_{\omega }^k,z^k\right)=\left\{\frac{p_D}{\overline{c}(z)}{\left(w_{\omega }^k\right)}^{T}q(z):z\in z^k\right\}$$

(40)

where

$\overline{c}(z)$ represents clutter density;

$w_j(z),q(z),m_j^k(z),P_j^k$ is the value after VB approximation.

The approximation of $q(\theta )$ is obtained via iterative solving using the VB method by setting immobile points. The posterior density, $q(\theta )$, is updated to $q^{i+1}(\theta )$ by calculating the expectation of the above equation using $q^i({\theta }_k^{\left(-\phi \right)})$ at the $i+1$ iteration. The joint PDF decomposition is expressed as follows:

$$\begin{array}{l}p\left({\theta }_k,Z_k\right)=p\left(Z_{k-1}\right){\left[N\left(x_k;m_k^{\omega ,j},P_k^{\omega ,j}/{\xi }_k\right)\right]}^{\left(1-t_k\right)}{\left[N\left(x_k;m_k^{\omega ,j},P_k^{\omega ,j}\right)\right]}^{t_k}\\ \hspace{1em}\hspace{1em}\hspace{1em} \times {\left[N\left(z_k;\dot{H}{x}_k,R_k/{\lambda }_k\right)\right]}^{\left(1-y_k\right)}{\left[N\left(z_k;{\dot{H}}_k{x}_k,R_k\right)\right]}^{y_k}\\ \hspace{1em}\hspace{1em}\hspace{1em} \times \mathrm{G}\left({\xi }_k;\frac{{\omega }_k}{2},\frac{{\omega }_k}{2}\right)\mathrm{G}\left({\lambda }_k;\frac{{\nu }_k}{2},\frac{{\nu }_k}{2}\right){\tau }_k^{t_k}{\left(1-{\tau }_k\right)}^{\left(1-t_k\right)}{\pi }_k^{y_k}{\left(1-{\pi }_k\right)}^{\left(1-y_k\right)}\\ \hspace{1em}\hspace{1em}\hspace{1em} \times \mathrm{Be}\left({\tau }_k;g_0,1-g_0\right)\mathrm{Be}\left({\pi }_k;e_0,1-e_0\right)\end{array}$$

(41)

where

$\dot{H}=\left[0_{1,i_{\omega ,j}^k-1},1\right]\otimes H$;

$G$ and $Be$ represent gamma and beta distributions;

$e_0$ and $g_0$ are priori mixed probabilities.

At each iteration, $\theta$ is updated sequentially. In order, update $q^{i+1}(x_k)$ for Gaussian distribution, update $q^{i+1}(t_k)$ and $q^{i+1}(y_k)$ for Bernoulli distribution, update $q^{i+1}({\xi }_k)$ and $q^{i+1}({\lambda }_k)$ for Gamma distribution, update $q^{i+1}({\tau }_k)$ and $q^{i+1}({\pi }_k)$ for Beta distribution, and finally determine the state vector and covariance matrix.

The noise covariance, ${\tilde{R}}_{k|k-1}^{(j)(n)}$, and estimation error covariance, ${\tilde{R}}_{k|k-1}^{(j)(n)}$, are corrected during the process:

$${\tilde{P}}_{k|k-1}^{(j)(n)}=\frac{P_k^{\omega ,j}}{{\mathrm{E}}^{(n)}\left[t_{1,k}\right]+\left(1-{\mathrm{E}}^{(n)}\left[t_{1,k}\right]\right){\mathrm{E}}^{(n)}\left[{\xi }_{1,k}\right]}$$

(42)

$${\tilde{R}}_{k|k-1}^{(j)(n)}=\frac{R_k}{{\mathrm{E}}^{(n)}\left[y_{2,k}\right]+\left(1-{\mathrm{E}}^{(n)}\left[y_{2,k}\right]\right){\mathrm{E}}^{(n)}\left[{\lambda }_{2,k}\right]}$$

(43)

The procedure for updating the potential distribution is the same as for GM-TCPHD, but the weight computation is altered, and the weights are updated by utilizing the noise matrix’s optimal value following VB approximation.

$$w_j(z)=\frac{p_D{w}_{\omega ,j}^k{q}_k^{(j)(n+1)}(z)〈{\mathsf{\Psi }}^1\left[w_{\omega }^k,z^k\backslash \left\{z\right\}\right],{\rho }_{{\omega }^k}〉}{\overline{c}(z)〈{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right],{\rho }_{{\omega }^k}〉}$$

(44)

$$q_k^{(j)(n+1)}\left(z_k\right)=\mathrm{N}\left(z_k;{\dot{H}}_k{m}_k^{\omega ,j},{\tilde{R}}_{k|k-1}^{(j)(n)}+{\dot{H}}_k{\tilde{P}}_{k|k-1}^{(j)(n)}{\dot{H}}_k^T\right)$$

(45)

$$K_k^{(j)(n+1)}={\tilde{P}}_{k|k-1}^{(j)(n)}{\dot{H}}_k^T{\left({\dot{H}}_k{\tilde{P}}_{k|k-1}^{(j)(n)}{\dot{H}}_k^T+{\tilde{R}}_{k|k-1}^{(j)(n)}\right)}^{-1}$$

(46)

$$m_k^{(j)(n+1)}=m_k^{\omega ,j}+K_k^{(j)(n+1)}\left(z_k-{\dot{H}}_k{m}_k^{\omega ,j}\right)$$

(47)

$$P_k^{(j)(n+1)}=\left[I-K_k^{(j)(n+1)}{\dot{H}}_k\right]{\tilde{P}}_{k|k-1}^{(j)(n)}$$

(48)

5. GSTM-TCPHD Filter Implementation with a Multi-Sensor

For multitarget tracking, we employ the GSTM-TCPHD filter for non-Gaussian heavy-tailed noise systems [43]. Nevertheless, the tracking field of view is constrained with a single sensor, and there is a lot of clutter information and little target measurement information [44]. Consequently, multi-sensor scenarios must be studied; however, as the number of sensors increases, so do the measurements, which can lead to a combinatorial explosion. The corresponding PHD components are updated following the subset of measurement selections, and the greedy algorithm is utilized for both path selection and subset division [45].

By enhancing the algorithm from the previous section, the MS-GSTM-TCPHD filter is suggested in this section, along with an implementation. Below are the MS-GSTM-TCPHD filter’s prediction and updating steps.

5.1. Prediction Process

Under assumptions P1–P4, MS-GSTM-TCPHD filter predictions are obtained, and the prediction process is the same as that of GSTM-TCPHD filters in single-sensor conditions.

At time $k-1$, the posteriori density is as follows:

$$D_{k-1}(X)={\displaystyle \sum _{j=1}^{J^{k-1}}{w}_j^{k-1}}\mathcal{N}\left(X;t_j^{k-1},m_j^{k-1},P_j^{k-1}\right)$$

(49)

The prediction is as follows:

$${\rho }_{k|k-1}(m)={\displaystyle \sum _{j=0}^m{\rho }_{{\beta }^k}}(m-j)\times {\displaystyle \sum _{l=j}^{\infty }\left(\begin{array}{c}n\\ j\end{array}\right)}{\rho }_{{\pi }^{k-1}}(n){\left(1-p_S\right)}^{l-j}{p}_S^j$$

(50)

$$D_{k|k-1}(X)=D_{{\beta }^k}(X)+p_S{\displaystyle \sum _{j=1}^{J^{k-1}}{w}_j^{k-1}}\mathcal{N}\left(X;t_j^{k-1},m_{\omega ,j}^k,P_{\omega ,j}^k\right)$$

(51)

5.2. Two-Step Greedy Approximation

We can use a two-step greedy approximation algorithm to solve the data association problem for multi-sensor data. The first step involves processing the subset of sensor measurements one at a time, calculating the scoring function for the current sensor measurement and the target, keeping the m measurements with the highest weights as the subset of measurements for the current target, and processing each subsequent subset of sensor measurements until the last sensor, as shown in Figure 3.

For a subset of multi-sensing measures, W, the likelihood function is expressed as follows:

$$f\left(W|x\right)={\displaystyle \prod _{(j,l)\in T_W}\frac{p_d^j(x)h_j\left(z_l^j|x\right)}{c_j\left(z_i^l\right)}}{\displaystyle \prod _{(j,l)\notin T_W}\left(1-p_d^j(x)\right)}$$

(52)

where

$T_W=\left\{(j,l):z_l^j\in W\right\}$;

$p_d^j(·)$ represents the detection probability of the j-th sensor;

$h_j$ denotes the measurement matrix of the j-th sensor.

The score function in the measurement subset, W, is expressed as follows:

$$dW=\frac{{\displaystyle \int w}(x)\left({\displaystyle \prod _{(i,l)\in T_W}{p}_d^i}(x)h_i\left(z_l^i|x\right)\right){\displaystyle \prod _{j:(j,*)\notin T_W}{q}_d^j}(x)dx}{{\displaystyle \prod _{(i,l)\in T_W}{c}_i}\left(z_l^i\right)}$$

(53)

where

$q_d^j(·)$ represents the detection probability of the j-th sensor;

$dW$ is the ratio of the likelihood function of the subset of target-generated measurements, $W$, to the likelihood function of the subset of clutter-generated measurements.

In order to avoid splitting the measurement into distinct components, the subset of measures that yield Gaussian components is re-partitioned in the second step. The set of optimal paths is chosen by keeping the top few partitioning hypotheses with the highest weights, and the set of disjoint measures obtained in the first step is chosen in the order of Gaussian components. Figure 4 depicts the second step as follows.

Given a partitioned set of measures, $W_i$, the weights are calculated as follows.

5.3. Update Process

Each partition, $P=\left\{w^1,\dots ,w^n\right\}$, is created via greedy approximation, and each Gaussian component is updated using a partition, P. The updated cardinality distribution and TPHD can be expressed as follows:

$${\rho }_{{\pi }^k}(n)=\frac{{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right](n){\rho }_{{\omega }^k}(n)}{〈{\mathsf{\Psi }}^0\left[w_{\omega }^k,z^k\right],{\rho }_{{\omega }^k}〉}$$

(54)

$$D_{{\pi }_S^k}(X)={\displaystyle \sum _{j=1}^{J^k}{w}_j^k}\mathcal{N}\left(X;t_j^k,m_j^k,P_j^k\right)$$

(55)

The updated optimal value of the noise covariance matrix for the state is used in the VB approximation of the posterior state, which is performed as demonstrated in the preceding section as part of the update process. The difference is that the update of the weights in it has been changed, and the formula is as follows:

$$w_j(z)=\frac{p_D{w}_{\omega ,j}^k{q}_k^{(j)(n+1)}(z)〈{\mathsf{\Psi }}^1\left[w_{\omega }^k,P_w^i\backslash \left\{w\right\}\right],{\rho }_{{\omega }^k}〉}{\stackrel{⌣}{c}(z)〈{\mathsf{\Psi }}^0\left[w_{\omega }^k,P_w^i\right],{\rho }_{{\omega }^k}〉}$$

(56)

6. Simulation

To compare the CPHD, TPHD, and TCPHD filters, the proposed GSTM-TCPHD filter is simulated. Table 1 contains the target parameters’ designs, while Table 2 contains the TPHD filter’s design. A linear 2D scene with multiple targets appearing, disappearing, omitted, and cluttered is observed. There are four birth regions at positions $\left[1000,1500\right]$, $\left[250,750\right]$, $\left[-250,1000\right]$, and $\left[-1500,250\right]$. There are eight targets in the scenario, each with different start and end times. The equations of motion of the target are the CV model; the CV model is a linear model, which makes it convenient for implementing target-tracking algorithms and simplifying calculations. Its state vector is represented as $x=\left[x,v_x,y,v_y\right]$, where $x, y$ and $v_x,v_y$ are the position and velocity of the target. Its state transfer equation is as follows:

$$x_k=F{x}_{k-1}+Gw$$

(57)

$$F=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],G=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right]$$

(58)

where

$w$ represents process noise;

$T$ represents the sampling interval of the sensor.

The initial process noise covariance matrix, $Q$, is as follows:

$$Q=diag\left(\left[G{G}^T\right]\right)$$

(59)

Table 1. Configuration of target parameters.

Parameter	Value
Degrees of freedom	10
State noise	5 m/s2
Measurement of mixing probabilities	0.7
State mixing probability	0.7
Measurement of heavy-tailed noise probability	0.10
State heavy-tailed noise probability	0.05
Distance error	15 m

Table 1. Configuration of target parameters.

Parameter	Value
Degrees of freedom	10
State noise	5 m/s2
Measurement of mixing probabilities	0.7
State mixing probability	0.7
Measurement of heavy-tailed noise probability	0.10
State heavy-tailed noise probability	0.05
Distance error	15 m

Table 2. TPHD filtering parameters.

Parameter	Value
Detection probability	0.9
Survival probability	0.98
Clutter	40
Tracking period	1 s
Pruning threshold	10 × 10−5
Maximum number of PHDs	100
Trajectory scanning steps	5

Table 2. TPHD filtering parameters.

Parameter	Value
Detection probability	0.9
Survival probability	0.98
Clutter	40
Tracking period	1 s
Pruning threshold	10 × 10−5
Maximum number of PHDs	100
Trajectory scanning steps	5

The sensor converts its observations into a two-dimensional distance measurement between the target and the sensor. The measurement model is as follows:

$$z_j=H_{j}x+{\epsilon }_j$$

(60)

where:

${\epsilon }_j$ represents the measurement noise of the j-th sensor.

$H_j$ represents the measurement matrix of the j-th sensor.

$$H_j=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$$

(61)

The measurement noise covariance matrix is expressed as follows:

$$R=r{I}_2$$

(62)

where $r=25$. Heavy-tailed noise in the system is caused by reflections of signals and instability in the sensors. As a result, we incorporate heavy-tailed noise points into our measurements with a 10% probability and a value that is 100 times greater than the measurement’s original value. The maneuver incorporates a heavy-tailed noise point with a probability of 5% and a value equal to 50 times the initial value that was predicted.

$$w_k\sim \left\{\begin{array}{c}N\left(0,Q\right),\\ N\left(0,100Q\right),\end{array}\right.\begin{array}{c}w.p. 1-p_{mo}\\ w.p. p_{mo}\end{array}$$

(63)

$${\epsilon }_k\sim \left\{\begin{array}{c}N\left(0,\mathrm{R}\right),\\ N\left(0,50\mathrm{R}\right),\end{array}\right.\begin{array}{c}w.p. 1-p_{po}\\ w.p. p_{po}\end{array}$$

(64)

where w.p. is the size of the probability, denoting heavy-tailed noise with probabilities $p_{mo}$ and $p_{po}$. Its model is shown in Figure 5.

The optimal subpattern assignment (OSPA) distance was used to assess the tracking results, which is defined as

$${\overline{d}}_p^{\left(c\right)}\left(X,Y\right)=\left\{\begin{array}{l}0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}m=n=0\\ {\left(\frac{1}{n}\left(\underset{\pi \in {\displaystyle {\prod }_n}}{\min }{\displaystyle \sum _{i=1}^m{d}^{(c)}}{(x_i,y_{\pi \left(i\right)})}^p+c^p(n-m)\right)\right)}^{\frac{1}{p}}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}m\le n\\ {\overline{d}}_p^{\left(c\right)}\left(Y,X\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}m>n\end{array}\right.$$

(65)

where

$c$ is the truncation distance; the parameter is set to 100;

$p$ is the order, which is set to 1.

The accuracy of both target position and target number estimations are measured by the OSPA distance.

6.1. Single-Sensor Simulation

According to the scene description, the target moves with a constant velocity model. If outlier value noise is present, the target’s motion model will change. Figure 6 shows the result of the motion path. From the simulation of the target path, it can be seen that the target’s maneuvering has significant changes at certain times, indicating that the target has strong maneuvering at these times.

In the absence of heavy-tailed noise, the noise exactly follows a Gaussian distribution. Several filters were simulated against this background, and the results are shown in Figure 7 and Figure 8. Then, we only introduced heavy-tailed noise into the measurements. The results of 100 Monte Carlo simulations are shown in Figure 9 and Figure 10.

In the scenario where noise outliers are not considered, Figure 7 and Figure 8 show that traditional methods can achieve accurate tracking. The number of targets and the OSPA distance simulated by the two strong maneuvering anti-outlier algorithms are relatively large. This is because the two filters deliberately increase the covariance calculation during the filtering process to improve the anti-noise performance, resulting in an increase in errors.

When process noise outliers are not considered, multiple algorithms can achieve stable tracking. Figure 9 illustrates that several algorithms have relatively small errors in estimating the number of targets. The GSTM-TCPHD algorithm has a slightly larger estimate of the number of targets than other methods, particularly when the number of targets increases and there is target extinction. Figure 10 shows that several algorithms have similar performance in estimating the OSPA distance.

Comparing the simulation results without considering outlier noise and without considering process noise outliers, when measurement noise outliers are added, the tracking number error and OSPA distance of several algorithms increase. The OSPA distance is composed of a distance error and a cardinality error. Therefore, when the estimated number of targets is similar, the improvement in OSPA distance performance mainly comes from the increased accuracy of target position estimation. This suggests that the proposed GSTM-TCPHD algorithm has higher position estimation accuracy in this scenario.

Simulations were conducted after adding heavy-tailed noise to the state and measurement noise. The results of 100 MC simulation runs are shown in Figure 11 and Figure 12. The figures illustrate the relatively large error and gradual tendency toward divergence in the standard TPHD/TCPHD algorithm. The GSTM-TCPHD filter performs exceptionally well.

Figure 11 and Figure 12 demonstrate that the GSTM-TCPHD proposed in this study outperforms the remaining filters in terms of overall performance, exhibiting the smallest tracking number error and OSPA distance error, as well as the most stable target tracking, when considering process outlier noise and measurement outlier noise. The tracking effects of the remaining filters tend to diverge with increasing tracking time.

Among the better-performing filters, the ST-TCPHD filter has a runtime of 21.91 s and an average OSPA distance of 53.19, while the GSTM-TCPHD filter has a runtime of 41.09 s and an average OSPA distance of 32.39. When comparing the simulation results while excluding outlier noise, only considering measurement outlier noise, and considering both process and measurement outlier noise, the addition of outlier noise leads to an increase in the number of target tracking errors and OSPA distance and a sharp decrease in the tracking effect. The proposed filter performs well in all three scenarios, particularly after the addition of outlier noise. It combines the characteristics of a short running time and high tracking accuracy.

The performance of the filter was further evaluated to investigate the impact of different measurement noises with the same process noise and different process noises with the same measurement noise. Figure 13 displays the simulation results. Figure 13 shows that process outlier noise has a greater impact on target-tracking accuracy than measurement outlier noise. OSPA distance is more sensitive to process noise when the probability of outlier noise changes. Our proposed method exhibits the smallest increase in OSPA distance and more stable tracking when the probability of outlier noise changes.

6.2. Multi-Sensor Simulation

To explore the performance of the proposed MS-GSTM-TCPHD filter, six sensors are arranged in the scene. The sensor positions are $\left[800,300\right]$, $\left[700,300\right]$, $\left[600,300\right]$, $\left[-100,1000\right]$, $\left[0,1000\right]$. The sensor location and the simulated primary path are shown in Figure 14, where the green hollow circle is the sensor location, and the black line represents the trajectory map of the primary simulation. The red circle is the born region. Each sensor has a certain probability of appearing as a quantitative wild value, and there is a total of eight target appearances and extinctions in the scene.

The tracking performance of the proposed algorithm in a multi-sensor background was verified. Firstly, a simulation was conducted without considering outlier noise, and the results are shown in Figure 15 and Figure 16. Secondly, a simulation was conducted when considering measurement outlier noise, and the results are shown in Figure 17 and Figure 18.

Figure 15, Figure 16, Figure 17 and Figure 18 demonstrate that the methods employed achieve complete tracking and reduce target-tracking error as measurement outlier noise increases.

Table 3. OSPA and runtime with different numbers of sensors under Gaussian noise.

Filter	OSPA for MS	OSPA for SS	Time for MS	Time for SS
CPHD	3.68	11.84	2.27 s	0.68 s
TCPHD	4.03	12.42	4.37 s	0.29 s
ST-TCPHD	4.19	23.50	119.18 s	17.70 s
GSTM-TCPHD	3.88	28.89	27.81 s	24.11 s

shows a comparison of the tracking results of a single sensor and multi-sensors running CPHD filters, TCHHD filters, ST-TCPHD filters, and GSTM-TCPHD filters in terms of average OSPA distance and running time without considering outlier noise.

Table 4. OSPA and runtime with different numbers of sensors under outlier measurement noise.

Filter	OSPA for MS	OSPA for SS	Time for MS	Time for SS
CPHD	4.08	17.49	2.74 s	0.79 s
TCPHD	4.35	18.71	5.15 s	0.28 s
ST-TCPHD	4.21	25.54	152.07 s	12.15 s
GSTM-TCPHD	3.92	17.54	48.88 s	34.61 s

shows the same comparison without considering process outlier noise. The comparison of computational efficiency and accuracy demonstrates that our proposed method outperforms the others.

The results of 100 Monte Carlo simulations, considering the process outlier noise and measurement outlier noise, are shown in Figure 19 and Figure 20. When considering both process outlier noise and measurement outlier noise simultaneously, the tracking error of traditional Gaussian noise filters increases sharply. Even under multiple sensors, tracking accuracy cannot be guaranteed. Compared with the operating results of a single sensor against the same background, a filter that can effectively track significantly reduces the OSPA distance, highlighting the advantages of multiple sensors.

Finally, the effect of the number of sensors is shown in Figure 21. The effect of different outlier noise probabilities on the target-tracking results was further analyzed. The simulation results are shown in Figure 22.

It has been verified that, by increasing the number of sensors, the error continues to decrease, and after reaching six sensors, the OSPA distance gradually converges to the limiting value. In a multi-sensor scenario, the change in the probability of outlier noise has little effect on the tracking results, and the stability is greater than that of a single sensor. The average OSPA distance and running time of the four methods are shown in

Table 5. OSPA and runtime with different numbers of sensors under outlier noise.

Filter	OSPA for MS	OSPA for SS	Time for MS	Time for SS
CPHD	34.71	51.77	2.74 s	0.65 s
TCPHD	35.55	55.47	5.15 s	0.19 s
ST-TCPHD	7.76	54.18	221.07 s	14.09 s
GSTM-TCPHD	6.75	32.52	48.66 s	41.91 s

. Even if the probability of outlier noise changes, the increase in OSPA distance is not significant. The resulting graph and table demonstrate that our proposed method enhances computational efficiency while maintaining tracking accuracy.

7. Conclusions

In order to handle the multitarget tracking problem in the presence of non-Gaussian heavy-tailed noise, in this paper, we propose a GSTM-TCPHD filter for a single-sensor background a the MS-GSTM-TCPHD filter for a multi-sensor background. By utilizing trajectory sets to create multitarget trajectories, the TCPHD filter allows us to model non-Gaussian heavy-tailed noise with a GSTM distribution. The true posterior probability is then approximated using VB approximation. We derive the GSTM-TCPHD filter’s implementation procedure under heavy-tailed noise scenarios. Furthermore, we create the MS-GSTM-TCPHD filter, which extends the GSTM-TCPHD filter to multi-sensor systems. We use a greedy approximation technique to ensure that no measurements with low-likelihood weights are missed in order to address the problem of data explosion in multi-sensor scenarios. This leads to a workable implementation of the MS-GSTM-TCPHD filter for multiple sensors.

We conducted Monte Carlo simulations of the filters to compare OSPA distances and cardinality errors. The simulation results show that the error of the CPHD/TPHD/TCPHD filter gradually increases in the presence of heavy-tailed noise, eventually leading to divergence. In contrast, the GSTM-TCPHD filter achieves stable tracking with minimal cardinality errors. Furthermore, compared with the GSTM-TCPHD filter, the MS-GSTM-TCPHD filter further reduces OSPA distance and exhibits improved tracking accuracy in multi-sensor scenarios. Additionally, we also compared filters using Student’s t-distribution to model heavy-tailed noise. According to the results, our suggested filter operates more efficiently and requires less time to operate, all while exhibiting smaller OSPA distance and cardinality errors. To sum up, a greater variety of scenarios can be covered by our suggested filter.