Non-Ellipsoidal Infrared Group/Extended Target Tracking Based on Poisson Multi-Bernoulli Mixture Filter and B-Spline

Abstract

This study provides a solution for multiple group/extended target tracking with an arbitrary shape. Many tracking approaches for extended/group targets have been proposed. However, these approaches make assumptions about the target shape, which have limitations in practical applications. To address this problem, in this work, an extended/group target tracking algorithm based on B-spline is proposed. Specifically, the extension of an extended or a group target was modeled as a spatial probability distribution characterized by the control points of a B-spline function that was then jointly propagated with the measurement rate model and kinematic component model over time using the Poisson multi-Bernoulli mixture (PMBM) filter framework. In addition, an amplitude-aided measurement partitioning approach is proposed to improve the accuracy caused by distance-based approaches. The simulation results demonstrate that the extension, shape and orientation of targets can be estimated better by the proposed algorithm, even if the shape changes. The tracking performance is also improved by about 10% and 13% compared to the other two algorithms.

1. Introduction

With the increase in human space activities, the number of space targets, including debris, abandoned spacecraft, missiles, and so forth is also increasing, thus affecting the safety of human space activities. This study focuses on tracking space targets based on space-based infrared sensors. Tracking problems related to space targets are part of the multi-target tracking problem. Multi-target tracking (MTT) is a hot topic both in military and civilian fields [1,2]. Most traditional tracking algorithms are based on the assumption that each target produces at most one measurement per time step limited by the resolution of the traditional sensors. However, with the rapid development of sensors, the resolution and accuracy have been greatly improved. Thus, the assumption no longer holds; that is, a target may generate multiple measurements at a given time. In this case, the target is preferably treated as an extended target [3] with a specific size, shape and orientation. When targets are close enough, they form an invisible group target, which acts as an extended target; therefore, the tracking of a group target can be regarded as the tracking of an extended target [4]. Generally, the tracking process of multiple group/extended targets contains three steps. Firstly, the target state is modeled. Secondly, the data received from sensors are processed, which involves partitioning the measurements, and data are associated with the measurements and the targets. Lastly, the iterative state parameters containing the predictions and updated processes of the parameters are outlined. The discussion that follows focuses on extended target tracking, which can also be applied to group targets.

There are two models that can be considered for extended targets: (1) a model of the number of measurements generated by each extended target and (2) a model of the spatial distribution of the targets. Generally, the number of measurements generated by the extended target is modeled using the Poisson distribution, whose rate parameter is modeled as a gamma distribution [5]. As for the target extension, it is critical for target identification in the tracking scenario and is the focus of the present study. It has been the subject of many studies.

When it comes to target extension, there are two assumptions in the literature: (1) the assumption that there are some general parameters that can describe the shape of the targets, such as an ellipse, rectangle [6], line [7], and so on, and (2) the assumption of an arbitrary shape. For the former, an ellipse is the most frequently used, and there are several approaches for ellipse modeling, such as random matrix (RM) [8] theory, which models the extension as a symmetric positive definite matrix. Generally, the gamma Gaussian inverse–Wishart (GGIW) distribution [5,9,10,11] is adapted. In these cases, the extended target is assumed to be elliptical. These methods are effective and have been applied to many scenarios [9]. However, their performance suffers when the shape is not elliptical. To estimate the non-elliptical extended target, [4,12] modeled the extension of the arbitrary shape using multiple ellipses (Em) and achieved the tracking of non-elliptical extended targets. However, a priori knowledge of the number of ellipses is needed. When the actual shape is unknown, it is challenging to choose an appropriate number of ellipses. To more accurately estimate the shape of the targets, random hypersurface models (RHMs) [13] were proposed by Baum. The shape constraint (i.e., elliptical constraint, star-convex constraint and level-set constraint) needs to be set to describe the target-extended state. That is, the extension can be obtained by estimating the parameters of the shape constraint. RHMs outperform RM models and can estimate extended targets with an irregular shape when using accurate preset shape constraints, but the computational complexity is higher. Many scholars have made great improvements based on these models. However, these methods all require assumptions about the shape of the target. In addition, some scholars have proposed other methods to estimate the shapes of extended targets. The extension–deformation approach [14] was proposed by Xiaorong to obtain the estimated shapes of targets. An extended target is considered to have a reference extension with control points on the boundary [14]. Then, the extension is estimated by moving some of the control points. The reference extension and the control points are a priori, which limits their application in practical scenarios due to the varying number of control points. Yulan Han applied an algorithm based on the level set and Gaussian surface fitting to extended target tracking [15]. Some researchers also use a B-spline curve to fit the shape [16,17]. A B-spline curve can be fitted to arbitrary shapes by adjusting the control points to estimate the shape of the extended target shape [17] assumed that the measurement points are generated using an edge, which is not applicable for some targets, especially group targets. In addition, non-uniform rational B-splines (NURBS) surfaces which can be considered as an extension of B-spline have been successfully applied to estimate a 3-dimensional (3D) target extension [18,19]. However, although NURBS is more flexible, it needs extra storage to define traditional curves and surfaces. In addition, weights have a significant impact on the shape estimation. Improper weights can cause shape distortion [20].

In addition, measurement partitioning is another important part of extended target tracking. Whether the measurement subset of each group target can be correctly divided at each specific moment determines the estimation accuracy of the group target state’s estimation accuracy and the tracking performance of the algorithm. Thus, distance partition, k-means ++ partition [21], prediction partition [22] and expectation maximization (EM) partition [21], DBSCAN [23] are proposed and successfully applied. However, the value of the distance threshold is difficult to determine. In addition, the distance information is insignificant in infrared images. For k-means ++, a priori knowledge of the selection and number of cluster centers is required. While the prediction information at the previous moment is inaccurate, such as the target maneuvering, the prediction partition and expectation maximization (EM) partition fail. For DBSCAN, the choice of two parameters is a challenge, and different parameter combinations have a significant effect on clustering. In addition, the performance degrades when the density of the sample set is not uniform. [23] proposed a grid-based DBSCAN algorithm which can deal with the non-equidistant sampling density and improve the performance. For infrared images, amplitude information is introduced to the partition measurements in this study.

Once models are defined, a multi-target filter is needed to implement the estimation of the target state, number and shape. To track multiple extended targets, random finite sets (RFSs) that model targets and measurements as random sets can provide an effective solution for data association and have received a great deal of attention in the literature. PMBM [24] and the Delta-generalized labeled multi-Bernoulli ($\delta$-GLMB) [25] are two well-established MTT conjugate priors that have better performance than others. Meanwhile, the simulation results show that PMBM filters outperform $\delta$-GLMB filters both in terms of their performance and computation cost [5,26]. Inspired by [16] and [5], in this study, a multiple extended target tracking algorithm based on B-spline and the PMBM filter is proposed. In addition, the amplitude-aided method is used for the measurement portioning part. Specifically, the target state is modeled in three parts. The first part denotes the number of measurements generated by an extended target, which is modeled using Poisson distribution, and the Poisson rate is modeled using gamma distribution. The second part is the kinematic state of the extended target center and is modeled using Gaussian distribution. The last part represents the extension of the extended targets and is modeled using a spatial probability distribution characterized by control points of a B-spline function. Then, the single target state is propagated by iterating parameters of models. Lastly, multi-target tracking is implemented under the PMBM filter.

This study focuses on tracking the multiple non-ellipsoidal targets based on the PMBM filter. There are three contributions of this work, which are as follows:

(1)

The B-spline is applied to model the extension of extended targets, thus solving the inaccurate modeling of targets with an arbitrary shape. In addition, the performance is also improved using this algorithm.

(2)

The amplitude information is introduced to partition the measurement, which can accurately partition the measurement set, especially when the targets are close.

(3)

The updated prediction and likelihood formulas of the algorithm based on the B-spline model are derived.

The remainder of the paper is organized as follows. In Section 2, the background of the proposed algorithm is provided. Section 3 provides the implementation of the proposed algorithm. Section 4 presents the numerical simulation to verify the proposed algorithm, and several methods are compared with the proposed algorithm. Section 5 provides the conclusion.

2. Background

This section provides the background on the implementation of the proposed algorithm.

2.1. PMBM Density

Once the models of kinematic state and extension and the number of measurements generated by an extended target are defined, the multiple extended target tracking algorithm can be proposed based on the models and the PMBM filter. The PMBM filter is a promising filter that models the undetected targets and the detected targets using a Poisson point process (PPP) and multi-Bernoulli mixture (MBM) process, respectively. Thus, the set of targets can be divided into two-joint subsets [27]:

$$X_k=X_k^u{\displaystyle \cup }{X}_k^d,$$

(1)

where $X_k^u$ denotes the undetected targets and $X_k^d$ is the detected targets at time $k$. The density can be expressed as

$$f_{k|k}\left(X_k|Z_k\right)={{\displaystyle \sum }}_{X_k^u\uplus X_k^d}{f}_{k|k}^P\left(X_k^u\right)f_{k|k}^{mbm}\left(X_k^d\right),$$

(2)

$$f_{k|k}^P\left(X_k^u\right)=e^{-\langle D^u;1\rangle }{{\displaystyle \prod }}_{x\in X_k^u}{D}^u=e^{-\mu \left(x\right)}{{\displaystyle \prod }}_{x\in X_k^u}\mu f\left(x\right),$$

(3)

$$f_{k|k}^{mbm}\left(X_k^d\right)={{\displaystyle \sum }}_{j\in \mathbb{J}}{\omega }_j{{\displaystyle \sum }}_{{\uplus }_{i\in {\mathbb{I}}^j}{X}^i=X_k^d}{{\displaystyle \prod }}_{i=1}^{\left|{\mathbb{I}}^j\right|}\begin{array}{c}{f}_{j,i}\left(X^i\right),\end{array}$$

(4)

where $f_{k|k}^P\left(X_k^u\right)$ denotes the Poisson density, $D^u$ denotes the intensity function of PPP, $\mu$ is the Poisson rate, $f\left(x\right)$ represents the spatial distribution, and $D^u=\mu f\left(x\right)$. $f^{mbm}\left(X^d\right)$ is a mixture of multi-Bernoulli denoting the potential targets that are detected at least once. $\mathbb{J}$ and ${\mathbb{I}}^j$are an index set for the multi-Bernoullis (MBs) in the MBM and an index set for the Bernoullis in the $jth$ MB, respectively. $\left|\mathbb{J}\right|$ is the number of MB components in the MBM, and $\left| {\mathbb{I}}^j\right|$ is the number of Bernoulli components in the $jth$ component. ${\omega }_j$ is the hypothesis weight, and ${\displaystyle \sum }{\omega }^{j,i}=1$, $f_{j,i}\left(X^i\right)$ denotes the $ith$ Bernoulli density in the $jth$ global hypothesis, which is defined by

$$f_{j,i}\left(X^i\right)=\{\begin{array}{c}1-r_{j,i} X^i=\varnothing \\ r_{j,i}{f}_{j,i}\left(x\right) X^i=\left\{x\right\}\\ 0 \left|X^i\right|\ge 2\end{array}$$

(5)

with $r_{j,i}$ denoting the probability of existence of the Bernoulli component and $f_{j,i}\left(x\right)$ denoting the corresponding density.

2.1.1. Standard Extended Target Measurement Model

$Z_k$ denotes the set of measurements at time $k$ and includes clutter and targets. Clutter and targets are assumed to be independent. The clutter at time $k$ is modeled using PPP, the Poisson rate is $\lambda$, and the spatial distribution is expressed as $c\left(z\right)$; thus, the intensity of clutter PPP is $\kappa \left(z\right)=\lambda c\left(z\right)$. A target with state $x_k$ can be detected at time $k$ with the probability $p_D\left(x_k\right)$. Generally, the measurements generated by an extended target are modeled using PPP, and the intensity is expressed as $\gamma \left(x_k\right)\varphi (z_k|x_k)$ if a target is detected, where $\gamma \left(x_k\right)$ is the Poisson rate and $\varphi (z_k|x_k)$ is the Poisson density.

Assuming that the set of measurements at time $k$ is non-empty, that is, $\left|Z_k\right|>0$, the measurement likelihood $l_z\left(x_k\right)$ of the extended target can be denoted by the product of detection probability and PPP density.

$$l_z\left(x_k\right)=p_D\left(x_k\right)p\left(Z_k|x_k\right)=p_D\left(x_k\right)e^{-\gamma \left(x_k\right)}{{\displaystyle \prod }}_{z\in Z_k}\gamma \left(x_k\right)\varphi (z_k|x_k),$$

(6)

The Poisson probability that an extended target with state $x_k$ will generate at least one measurement at time $k$ is $1-e^{-\gamma \left(x_k\right)}$; thus, the effective detection probability of an extended target is $p_D\left(x_k\right)\left(1-e^{-\gamma \left(x_k\right)}\right)$. Accordingly, the probability of not being detected is

$$q_D\left(x_k\right)=1-p_D\left(x_k\right)+p_D\left(x_k\right)e^{-\gamma \left(x_k\right)},$$

(7)

In addition, the likelihood of an empty set of measurements is $l_{\varnothing }\left(x_k\right)=q_D\left(x_k\right)$.

2.1.2. Standard Extended Target Dynamic Model

The targets survive from time $k-1$ to time $k$ with the probability $p_s({\xi }_{k-1})$, which is called the survival probability and moves with a single target transition density $f_{k|k-1}(\xi |\zeta )$ of state $\xi$ given $\zeta$, where $\xi$ is an extended state of the extended target consisting of kinematic state and extension. The newly born target is modeled using a PPP with the intensity $D_k^b\left(\xi \right)$.

2.1.3. Amplitude-Aided Measurement Partitioning

The traditional algorithm used to partition the measurements is based on the distance. However, the distance threshold setting presents challenges. Setting the threshold too large will affect the estimation of the number of targets, and setting it too small will increase the complexity. Figure 1 gives the partitioning results, corresponding to three different distance thresholds. The left column, that is, (a), (c) and (e), is the original image, the right column, that is, (b), (d) and (f), is a partial enlargement of the image.

Extended targets are presented as multiple adjacent pixels in the infrared image plane with connectivity. Thus, the measurement partitioning of the infrared image can use the characteristics of connectivity and amplitude. As shown in Figure 2a, there are extended targets and clutter. Figure 2b is the enlarged local image which shows the details of the measurements. Figure 2c gives the amplitude of these clusters, which shows the difference between targets and clutter. Clearly, clutter occupies small areas and has a low amplitude. Thus, the extended target measurements can be obtained. In this study, all the measurements with an area smaller than $3 pixel\times 3 pixel$ are eliminated. Additionally, those with an amplitude less than $\tau =20$ are also eliminated. The partitioning steps are as follows.

Step 1: The Hoshen–Kopelman (H-K) [28] algorithm is used to detect the connected regions in the image plane. The results are shown in Figure 3a,b. Figure 3a is the partitioning results of the 36th frame; Figure 3b is the enlarged image of two close extended targets which shows that the blue ellipse contains multiple local maxima, which can be further divided.

Step 2: Detect the local maxima in the clusters generated by H-K algorithm. As shown in Figure 2c, when the targets are close, the amplitude is superimposed, and wave peaks are formed. Thus, if there are multiple local maxima, each local maximum and its neighborhood contains at least one extended target.

Step 3: The clusters containing multiple local maxima are further partitioned using the k-means algorithm to produce new clusters. The results of partitioning are shown in Figure 3c,d. Figure 3c presents the partitioning results of frame 36 based on H-K with k-means; Figure 3d is an enlarged image of two close extended targets (H-K with k-means).

2.1.4. PMBM Filter for Extended Target

The recursion of the PMBM filter for the extended target is presented, which includes the prediction process and updated process. The PMBM RFS is given in the Ref. [5] for multiple extended targets and in the Ref. [27] for multiple point targets. According to these studies, the density is propagated by $D_k^u$, ${\left\{{\omega }_k^j,{\left\{r_k^{j,i},f_k^{j,i}\right\}}_{i\in {\mathbb{I}}^j}\right\}}_{j\in \mathbb{J}}$.

A.

Prediction process

Given the posterior PMBM density at time $k-1$ is $D_{k-1}^u$, ${\left\{{\omega }_{k-1}^j,{\left\{r_{k-1}^{j,i},f_{k-1}^{j,i}\right\}}_{i\in {\mathbb{I}}^j}\right\}}_{j\in \mathbb{J}}$ and the dynamic model in Section 2.1.2, the predicted parameters at time $k$ can be expressed as

$$D_{k|k-1}^u=D_k^b+{\displaystyle \int }{D}_{k-1}^u\left(\xi \right)p_s{f}_{k|k-1}(\xi |\zeta )d\xi d\zeta ,$$

(8)

$$r_{k|k-1}^{j,i}=r_{k-1}^{j,i}{\displaystyle \int }{f}_{k-1}^{j,i}\left(\zeta \right)p_{s}d\zeta ,$$

(9)

$$f_{k|k-1}^{j,i}=\frac{{\displaystyle \int }{f}_{k|k-1}\left(\xi |\zeta \right)p_s{f}_{k-1}^{j,i}\left(\zeta \right)d\zeta }{{\displaystyle \int }{p}_s{f}_{k-1}^{j,i}\left(\zeta \right)d\zeta },$$

(10)

$${\omega }_{k|k-1}^j={\omega }_{k|k-1}^{j,i},$$

(11)

and the detailed proof process can be found in the Ref. [24].

B.

Updated process

Assuming that the predicted density is $D_{k|k-1}^u$, ${\left\{{\omega }_{k|k-1}^j,{\left\{r_{k|k-1}^{j,i},f_{k|k-1}^{j,i}\right\}}_{i\in {\mathbb{I}}_{k|k-1}^j}\right\}}_{j\in {\mathbb{J}}_{k|k-1}}$, the measurement model and measurement set is given as Section 2.1.2, and the updated PMBM density is also in PMBM form.

$$f_{k|k}\left(X_k|Z_k\right)={{\displaystyle \sum }}_{X_k^u\uplus X_k^d}{f}_{k|k}^P\left(X_k^u\right)f_{k|k}^{mbm}\left(X_k^d\right),$$

(12)

$$f_{k|k}^P\left(X_k^u\right)=e^{-\langle D^u;1\rangle }{{\displaystyle \prod }}_{x\in X_k^u}{D}^u,$$

(13)

$$f_{k|k}^{mbm}\left(X_k^d\right)={{\displaystyle \sum }}_{j\in \mathbb{J}}{\omega }_j{{\displaystyle \sum }}_{{\uplus }_{i\in {\mathbb{I}}^j}{X}^i=X_k^d}{{\displaystyle \prod }}_{i=1}^{\left|{\mathbb{I}}^j\right|}\begin{array}{c}{f}_{j,i}\left(X^i\right),\end{array}$$

(14)

where for the undetected targets,

$$D_{k|k}^u=q_D{D}_{k|k-1}^u\left(\zeta \right),$$

(15)

and for the detected targets (including the first detected targets, the missed detected targets and the targets detected in all the steps),

$$r_{k|k}^j=\{\begin{array}{cc}1\hfill & \hfill C{\displaystyle \cap }{\mathbb{I}}^j=\varnothing ,\left|C_c\right|>1\\ \frac{{\displaystyle \int }{D}_{k|k-1}^u\left(\zeta \right)l_{C_c}d\zeta }{{\kappa }^{C_c}+{\displaystyle \int }{D}_{k|k-1}^u\left(\zeta \right)l_{C_c}d\zeta }\hfill & \hfill C{\displaystyle \cap }{\mathbb{I}}^j=\varnothing ,\left|C_c\right|=1\\ \frac{r_{k|k-1}^{j,i}{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right)q_{D}d\zeta }{1-r_{k|k-1}^{j,i}+r_{k|k-1}^{j,i}{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right)q_{D}d\zeta }\hfill & \hfill C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|=\varnothing \\ 1\hfill & \hfill C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|\ne \varnothing \end{array},$$

(16)

$$f_{k|k}^{j,i}=\{\begin{array}{c}\frac{l_{C_c}{D}_{k|k-1}^u\left(\zeta \right)}{{\displaystyle \int }{D}_{k|k-1}^u\left(\zeta \right)l_{C_c}d\zeta } C{\displaystyle \cap }{\mathbb{I}}^j=\varnothing \\ \frac{q_D{f}_{k|k-1}^{j,i}\left(\zeta \right)}{{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right)q_{D}d\zeta } C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|=\varnothing \\ \frac{l_{C_c}{f}_{k|k-1}^{j,i}\left(\zeta \right)}{{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right)l_{C_c}d\zeta } C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|\ne \varnothing \end{array},$$

(17)

$${\omega }_{k|k,A}^j=\frac{{\omega }_{k|k-1}^j{{\displaystyle \prod }}_{c\in A}{L}_c}{{{\displaystyle \sum }}_{{{\displaystyle \cup }}_{j\in {\mathbb{J}}_{k|k-1}}}{{\displaystyle \sum }}_{A\in A^j}{\omega }_{k|k-1}^j{{\displaystyle \prod }}_{c\in A}{L}_c},$$

(18)

$$L_c=\{\begin{array}{c}{\kappa }^{C_c}+{\displaystyle \int }{D}_{k|k-1}^u\left(\zeta \right)l_{C_c}d\zeta C{\displaystyle \cap }{\mathbb{I}}^j=\varnothing ,\left|C_c\right|=1\\ {\displaystyle \int }{D}_{k|k-1}^u\left(\zeta \right)l_{C_c}d\zeta C{\displaystyle \cap }{\mathbb{I}}^j=\varnothing ,\left|C_c\right|>1\\ 1-r_{k|k-1}^{j,i}+r_{k|k-1}^{j,i}{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right) q_{D}d\zeta C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|=\varnothing \\ r_{k|k-1}^{j,i}{\displaystyle \int }{f}_{k|k-1}^{j,i}\left(\zeta \right) l_{C_c}d\zeta C{\displaystyle \cap }{\mathbb{I}}^j\ne \varnothing ,\left|C_c\right|\ne 1\end{array},$$

(19)

where $l_{C_c}$ is the predicted likelihood, and $C_c$ denotes the measurement cell which corresponds to the data association $c$. $\left|C_c\right|$ is the number of measurements in $C_c$. ${\omega }_{k|k,A}^j$ represents the hypothesis weight. ${\kappa }^{C_c}$ is the clutter intensity. $A^j$ is the space of all data associations $A$ for the $j$th predicted MB. $C\in A$ denotes non-empty cells, which contains measurements from a source (either clutter or a single target).

2.2. B-Spline

A B-spline is a piecewise polynomial function which can describe the shape of any curve by adjusting the location of the control points. More information on the B-spline can be found in the Refs. [29,30]. The B-spline curve with $l$ order can be described as

$$s\left(t\right)={{\displaystyle \sum }}_{i=1}^n{c}_i{N}_{i,l}\left(t\right),$$

(20)

where $c_i\in C=\left[c_1,c_2,\dots c_n\right]$ is a control point, $n$ denotes the total number of control points and $N_{i,l}\left(t\right)$ is the B-spline basis function, which is defined over a knot vector. The mathematical expression is as follows:

$$N_{i,l}\left(t\right)=\frac{t-t_i}{t_{i+l-1}-t_i}{N}_{i,l-1}\left(t\right)+\frac{t_{i+l}-t}{t_{i+l}-t_{i+1}}{N}_{i+1,l-1}\left(t\right),$$

(21)

$$N_{i,1}\left(t\right)=\{\begin{array}{c}1 t_i\le t\le t_{i+1}\\ 0 otherwise\end{array},$$

(22)

where the variables $t_i$ are knot elements.

Generally, the B-spline curve is open, as Figure 4a shows, and a close curve can be obtained by repeating the first three control points, as Figure 4b shows. The B-spline used in this study is a closed curve, and the order is 3, that is, $l=3$. The B-spline approach has been used in target tracking applications [16,17,31,32,33] in a continuous state space.

In this study, the control points are obtained as in the Ref. [16]. As shown in Figure 5 [34], suppose that ${\tilde{Z}}_k$ is a pseudo-measurement set. More information on the construction of the pseudo-measurement set can be found in the Ref. [16]. The mean of ${\tilde{Z}}_k$ is taken as the coordinates’ origin, and $\left[0,2\pi \right]$ is divided into $n$ equal angles. Thus, the pseudo-measurement set can be divided by these angles. Let $W_{k,i}$ be the measurement set in the area of the $i$th partition at time $k$, which can be calculated by

$$W_{k,i}=\left\{z_{k,j}|d_{L_{i,1}}\left(z_{k,j}\right)<d,C\left(z_{k,j}\right)=true\right\},$$

(23)

$$d_{L_{i,1}}\left(z_{k,j}\right)=\frac{\Vert b1z_{k,j}^{\left(1\right)}+b2z_{k,j}^{\left(2\right)}\Vert }{\sqrt{b{1}^2+b{2}^2}},$$

(24)

$$C\left(z_{k,j}\right)=\{\begin{array}{c}true L_{i,2}:a1z_{k,j}^{\left(1\right)}+a2z_{k,j}^{\left(2\right)}>0,{\theta }_i\in \left(0,\pi \right]\\ or L_{i,2}:a1z_{k,j}^{\left(1\right)}+a2z_{k,j}^{\left(2\right)}<0,{\theta }_i\in \left[\pi ,2\pi \right]\\ false otherwise\end{array},$$

(25)

where $z_{k,j}$ is the measurement in ${\tilde{Z}}_k$, $z_{k,j}^{\left(1\right)}$ and $z_{k,j}^{\left(2\right)}$, denoting the position of the x- and y-coordinates at time $k$. $d$ gives the width of a partition. $L_{i,1}$ (the red dotted line) is the line which goes through the origin of the coordinate and along the $ith$ partition angle direction; it is defined as $L_{i,1}:b1x+b2y=0$, where $\frac{b1}{b2}=-\tan \left({\theta }_i\right)$. $L_{i,2}$ (the blue dotted line) is the line which is perpendicular to line $L_{i,1}$ and also goes through the orign of the coordinate. $L_{i,2}$ is defined as $a1x+a2y=0$, where $\frac{a2}{a1}=\tan \left({\theta }_i\right)$.

Let ${\rho }_i$ and ${\theta }_i$ denote the polar radius and polar angles, respectively. ${\rho }_i$ can be obtained by

$${\rho }_i=\frac{2}{\left|W_{k,i}\right|}{{\displaystyle \sum }}_j^{\left|W_{k,i}\right|}\frac{\Vert a1z_{k,j}^{\left(1\right)}+a2z_{k,j}^{\left(2\right)}\Vert }{\sqrt{a{1}^2+a{2}^2}},$$

(26)

and thus, the matrix of control points can be expressed by $M_k={\left\{\left({\rho }_{k,i},{\theta }_{k,i}\right)\right\}}_{i=1}^n$, which can describe the shape of the targets. The matrix of control points is updated and predicted by implementing a one-dimensional Kalman filter. The specific equations can be found in the Ref. [16] (Equations (13)–(15)). The target shape can be estimated according to $M_k$. First, convert the control points to Descartian coordinates, then fit the shape using a B-spline curve.

$$P_{k,i}={\left[{\rho }_{k,i}·\cos \left({\theta }_{k,i}\right),{\rho }_{k,i}·\sin \left({\theta }_{k,i}\right)\right]}^T,$$

(27)

$$s\left(t\right)={\displaystyle \sum }_{i=1}^n{P}_{k,i}{N}_{i,l}\left(t\right),$$

(28)

where $P_{k,i}$ is the control point positon in Descatian coordinates, and let $C_k={\left\{P_{k,i}\right\}}_i^n$. In addition, to obtain the closed curve, add three points $P_{k,n+1}$, $P_{k,n+2}$, $P_{k,n+3}$, and make them equal to $P_{k,1}$, $P_{k,2}$ and $P_{k,3}$, respectively.

3. Proposed Algorithm

The proposed algorithm is presented in Figure 6, which can be divided into three steps: (1) amplitude-aided measurement partitioning; (2) shape modeling based on the control points of the B-spline; and (3) a multi-target tracking algorithm based on the PMBM filter. Step 1 can be implemented based on Section 2.1.3. Step 2 is performed in 2.2.

3.1. Single Extended Target

An extended target state ${\xi }_k$ at time $k$ is modeled as ${\xi }_k=\left({\gamma }_k,x_k,E_k\right)$, where ${\gamma }_k$ is the Poisson measurement rate parameter, $x_k$ is the kinematic state of the group center and $E_k$ denotes the shape state. As in the Ref. [5], the density of the rate parameter and kinematic state was modeled using the gamma distribution and Gaussian distribution, respectively. The shape state was modeled using a spatial probability distribution characterized by the control points of a B-spline. Thus, the target state distribution based on this model is a gamma-Gaussian-spline (GGS) and is denoted as

$$f_k\left({\xi }_k\right)=\mathcal{G}\left({\gamma }_k;{\alpha }_k,{\beta }_k\right)\mathcal{N}\left(x_k;m_k,P_k+{\mathsf{\Delta }}_k\right)\mathcal{S}\left(E_k;C_k\right),$$

(29)

where ${\alpha }_k$ and ${\beta }_k$ are parameters of gamma distribution, $m_k$ and $P_k$ denote the mean and covariance of Gaussian distribution, ${\mathsf{\Delta }}_k$ is the covariance of the shape and $C_k$ denotes control points.

The measurement likelihood for a single measurement $z_k$ at time $k$ is expressed as

$$\varphi \left(z_k|{\xi }_k\right)=\mathcal{N}\left(z_k;H_k{x}_k,X_k\right),$$

(30)

where $H_k$ is a known measurement model. The updated process and prediction process of the parameter are given in

Table 1. Updated GGS parameters of a single extended target.

Input ${\xi }_{k|k-1}$ and set of detection $Z_k$

${\xi }_{k|k}=\{\begin{array}{c}{\alpha }_{k|k}={\alpha }_{k|k-1}+\left|Z_k\right|\\ {\beta }_{k|k}={\beta }_{k|k-1}+1\\ m_{k|k}=m_{k|k-1}+K_{k|k-1}\left({\overline{z}}_k-H_k{m}_{k|k-1}\right)\\ P_{k|k}=P_{k|k-1}+K_{k|k-1}{H}_k{P}_{k|k-1}\\ C_{k|k}=\begin{array}{c}{C}_{k|k-1}+K_{k|k-1}\left(Z_k-H{C}_{k|k-1}\right) \end{array}\end{array}$

where

${\overline{z}}_k=\frac{1}{\left|Z_k\right|}{\displaystyle {\displaystyle \sum }_{z_k\in Z_k}}{z}_k$ $K_{k|k-1}=P_{k|k-1}{H}_k^T{\left(S_{k|k-1}\right)}^{-1}$ $S_{k|k-1}=H_k{P}_{k|k-1}{H}_k^T+\frac{{\mathsf{\Delta }}_k}{\left|Z_k\right|}$ ${\mathsf{\Delta }}_k=\frac{1}{\left|Z_k\right|-1}{\displaystyle {\displaystyle \sum }_{z_k\in Z_k}}\left(z_k-{\overline{z}}_k\right){\left(z_k-{\overline{z}}_k\right)}^T$

Likelihood:

${ℒ}_k=\frac{1}{\left|Z_k\right|!}\frac{\mathsf{\Gamma }\left({\alpha }_{k|k}\right){\left({\beta }_{k|k-1}\right)}^{{\alpha }_{k|k-1}}}{\mathsf{\Gamma }\left({\alpha }_{k|k-1}\right){\left({\beta }_{k|k}\right)}^{{\alpha }_{k|k}}}$

Output: ${\xi }_{k|k}$ and likelihood

and

Table 2. Predicted GGS parameters of a single extended target.

Input: ${\xi }_{k-1}$

${\xi }_{k|k-1}=\{\begin{array}{c}{\alpha }_{k|k-1}=\frac{{\alpha }_{k-1}}{\eta }\\ {\beta }_{k|k-1}=\frac{{\beta }_{k-1}}{\eta }\\ m_{k|k-1}=F_{k-1}{m}_{k-1}\\ P_{k|k-1}=F_{k-1}{P}_{k-1}{\left(F_{k-1}\right)}^T+Q\\ C_{k|k-1}=F_{k-1}{C}_{k-1}+{\omega }_{k-1}\end{array}$

Output: ${\xi }_{k|k-1}$

where $F_{k-1}={\nabla }_{x}f\left(x\right){|}_{x=m_{k-1}}$ is the transition matrix and $Q$ is the process noise.

, respectively.

3.2. GGS-PMBM Filter

The multiple extended target tracking algorithm is proposed based on the models presented above and the PMBM filter. The recursion is based on the following assumptions: (1) the survival probability is state-independent, that is, $p_{s,k}\left(\xi \right)=p_{s,k}$; (2) the detection probability is state-independent, that is, $p_{d,k}\left(\xi \right)=p_{d,k}$; and (3) the density of new targets at time $k$ is in GGS form.

The prediction process and updated process of the parameters of PPP components and Bernoulli components are given in the following:

Prediction step: Given the posterior GGS-PMBM density at time $k-1$ as $D_{k-1}^u$, ${\left\{{\omega }_{k-1}^j,{\left\{r_{k-1}^{j,i},f_{k-1}^{j,i}\right\}}_{i\in {\mathbb{I}}^j}\right\}}_{j\in \mathbb{J}}$, $D_{k-1}^u={\displaystyle \sum }_{j=1}^{N_k^u}{\omega }_{k-1}^{\left(u,j\right)}GGS\left(\xi ;{\zeta }_{k-1}^{\left(u,j\right)}\right)$, $f_{k-1}^{j,i}={\displaystyle \sum }_{j=1}^{N_k^u}{\omega }_{k-1}^{\left(u,j\right)}GGS\left(\xi ;{\zeta }_{k-1}^{\left(u,j\right)}\right)$.

Undetected targets: Poisson rate of predicted PPP:

$${\mu }_{k|k-1}^u={\mu }_k^b+{\mu }_{k-1}^u{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}{\omega }_{k-1}^{\left(u,j\right)}{p}_s,$$

(31)

The predicted spatial distribution is

$$\begin{gathered} f_{k|k-1}^u\left({\xi }_{k|k-1}\right)=\frac{{\mu }_k^b}{{\mu }_k^b+{\mu }_{k-1}^u{{\displaystyle \sum }}_{j=1}^{N_{k-1}^u}{\omega }_{k-1}^{\left(u,j\right)}{p}_s}{{\displaystyle \sum }}_{j=1}^{N_k^b}{\omega }_k^{\left(b,j\right)}GGS({\xi }_k;{\zeta }_k^{\left(b,j\right)}) \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{{\mu }_{k-1}^u}{{\mu }_k^b+{\mu }_{k-1}^u{{\displaystyle \sum }}_{j=1}^{N_{k-1}^u}{\omega }_{k-1}^{\left(u,j\right)}{p}_s}{{\displaystyle \sum }}_{j=1}^{N_k^u}{\omega }_{k-1}^{\left(u,j\right)}{p}_{s}GGS\left(\xi ;{\zeta }_{k|k-1}^{\left(u,j\right)}\right) \end{gathered}$$

(32)

where the predicted parameters ${\zeta }_{k|k-1}^{\left(u,j\right)}$ are computed as in Table 2. ${\zeta }_k^{\left(b,j\right)}$ is the parameter of new targets. ${\mu }_k^b$ and ${\mu }_{k-1}^u$ are the Poisson rate of birth targets and undetected targets, respectively. $N_k^b$ and $N_k^u$ are the number of birth targets and undetected targets, respectively.

Detected targets: Prediction weights and number of components are unchanged. The probability of existence and spatial distribution are

$$r_{k|k-1}^{j,i}=p_s{r}_{k-1}^{j,i}$$

(33)

$$f_{k|k-1}^{j,i}\left({\xi }_{k|k-1}\right)=GGS(\xi ;{\zeta }_{k|k-1}^{j,i})$$

where the parameters ${\zeta }_{k|k-1}^{j,i}$ can be obtained using Table 2.

Updated step: Suppose that the predicted parameters are given as $D_{k|k-1}^u$, ${\left\{{\omega }_{k|k-1}^j,{\left\{r_{k|k-1}^{j,i},f_{k|k-1}^{j,i}\right\}}_{i\in {\mathbb{I}}_{k|k-1}^j}\right\}}_{j\in {\mathbb{J}}_{k|k-1}}$, $D_{k|k-1}^u={\displaystyle \sum }_{j=1}^{N_k^u+N_k^b}{\omega }_{k|k-1}^{\left(u,j\right)}GGS\left(\xi ;{\zeta }_{k|k-1}^{\left(u,j\right)}\right)$, $f_{k|k-1}^{j,i}=GGS(\xi ;{\zeta }_{k|k-1}^{j,i})$.

The updated parameters for the detected targets and undetected targets are obtained using the following:

Undetected targets:

$${\mu }_{k|k}^u={\mu }_{k|k-1}^u{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}{q}_D^{u,j},$$

(34)

where $q_D^{u,j}$ is the probability that the target is not detected, and is defined as

$$q_D^{u,j}=1-p_D\left({\xi }_{k|k-1}^{\left(u,j\right)}\right)+p_D\left({\xi }_{k|k-1}^{\left(u,j\right)}\right){\left(\frac{{\beta }_{k|k-1}^{\left(u,j\right)}}{{\beta }_{k|k-1}^{\left(u,j\right)}+1}\right)}^{{\alpha }_{k|k-1}^{\left(u,j\right)}}.$$

(35)

The updated density is

$$\begin{array}{c}{f}_{k|k}^u\left(\xi \right)=\frac{{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}\left(1-p_D\left({\xi }_{k|k-1}^{\left(u,j\right)}\right)\right){\omega }_{k|k-1}^{\left(u,j\right)}}{{{\displaystyle \sum }}_{j^{\prime }=1}^{N_{k|k-1}^u}{q}_D^{u,j^{\prime }}{\omega }_{k|k-1}^{\left(u,j\right)}}GGS\left({\xi }_k;{\zeta }_{k|k}^{\left(u,j\right)}\right)+\frac{{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}{p}_D\left({\xi }_{k|k-1}^{\left(u,j\right)}\right){\left(\frac{{\beta }_{k|k-1}^{\left(u,j\right)}}{{\beta }_{k|k-1}^{\left(u,j\right)}+1}\right)}^{{\alpha }_{k|k-1}^{\left(u,j\right)}}{\omega }_{k|k-1}^{\left(u,j\right)}}{{{\displaystyle \sum }}_{j^{\prime }=1}^{N_{k|k-1}^u}{q}_D^{u,j^{\prime }}{\omega }_{k|k-1}^{\left(u,j\right)}}\hfill \\ \times \mathcal{G}\left({\gamma }_k;{\alpha }_{k|k-1}^{\left(u,j\right)},{\beta }_{k|k-1}^{\left(u,j\right)}+1\right)\times \mathcal{N}\left(x_k;m_{k|k-1}^{\left(u,j\right)},P_{k|k-1}^{\left(u,j\right)}\right)\times \mathcal{S}\left(E_k;C_{k|k-1}^{\left(u,j\right)}\right)\end{array}$$

(36)

where ${\omega }_{k|k-1}^{\left(u,j\right)}={\omega }_{k-1}^{\left(u,j\right)}{p}_s$. There are two parts to the density. The first part corresponds to the detection process modeled using $p_D$, which means missed detection. The second means that the Poisson random number of detections is zero.

Targets detected for the first time: A target is detected for the first time, and the set of detection is denoted by $D$. The existence probability and spatial distribution can be expressed as

$$r_D=\{\begin{array}{c}1 \left|D\right|>1\\ \frac{{ℒ}_D}{{\kappa }^D+{ℒ}_D} \left|D\right|=1\end{array}$$

(37)

$$f_D\left(\xi \right)=\frac{{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}{\omega }_{k|k-1}^{u,j}{ℒ}_k^{\left(u,j,D\right)}{p}_{D}GGS\left(\xi ;{\xi }_{k|k}^{\left(u,j,D\right)}\right)}{{{\displaystyle \sum }}_{j=1}^{N_{k|k-1}^u}{\omega }_{k|k-1}^{u,j}{ℒ}_k^{\left(u,j,D\right)}{p}_D},$$

(38)

where

$${ℒ}_D=\{\begin{array}{c}{\displaystyle \sum }_{i=1}^{N_{k|k-1}^u}{p}_D{\omega }_{k|k-1}^{u,j}{ℒ}_k^{\left(u,j,D\right)} \left|D\right|>1\\ {\displaystyle \sum }_{i=1}^{N_{k|k-1}^u}{p}_D{\omega }_{k|k-1}^{u,j}{ℒ}_k^{\left(u,j,D\right)}+{\kappa }^D \left|D\right|=1\end{array},$$

(39)

${\xi }_{k|k}^{\left(u,j,D\right)}$ and ${ℒ}_k^{\left(u,j,D\right)}$ are obtained using Table 1.

Targets detected in all steps: The Bernoulli component that exists all the time can be updated by a non-empty measurement set $V$. The probability of existence is

$$r_{k|k}^{j,i,V}=1,$$

(40)

and the spatial distribution can be expressed as

$$f_{k|k}^{j,i,V}\left(\xi \right)=GGS\left(\xi ;{\xi }_{k|k}^{j,i,V}\right),$$

(41)

$${ℒ}_V=r_{k|k-1}^{j,i}{ℒ}_k^{\mathrm{j},\mathrm{i},\mathrm{V}},$$

(42)

${\xi }_{k|k}^{j,i,V}$ and ${ℒ}_k^{\mathrm{j},\mathrm{i},\mathrm{V}}$ are computed in Table 1.

Targets that missed detection: If the measurement set used to update the $i$th Bernoulli in $j$th MB components is empty, that is, $V=\varnothing$, the existence probability is

$$r_{k|k}^{j,i,V}=\frac{r_{k|k-1}^{j,i}{q}_D^{j,i}}{1-r_{k|k-1}^{j,i}+r_{k|k-1}^{j,i}{q}_D^{j,i}},$$

(43)

where

$$q_D^{j,i}=1-p_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right)+p_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right){\left(\frac{{\beta }_{k|k-1}^{\left(j,i\right)}}{{\beta }_{k|k-1}^{\left(j,i\right)}+1}\right)}^{{\alpha }_{k|k-1}^{\left(j,i\right)}}.$$

(44)

The updated spatial distribution is

$$\begin{array}{c}{f}_{k|k}^{j,i,V}\left(\xi \right)=\frac{1-p_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right)}{q_D^{j,i}}GGS\left({\xi }_k;{\zeta }_{k|k-1}^{\left(j,i\right)}\right)+\frac{p_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right){\left(\frac{{\beta }_{k|k-1}^{\left(j,i\right)}}{{\beta }_{k|k-1}^{\left(u,j\right)}+1}\right)}^{{\alpha }_{k|k-1}^{\left(j,i\right)}}}{q_D^{j,i}}\times \mathcal{G}\left({\gamma }_k;{\alpha }_{k|k-1}^{\left(j,i\right)},{\beta }_{k|k-1}^{\left(j,i\right)}+1\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \times \mathcal{N}\left(x_k;m_{k|k-1}^{\left(j,i\right)},P_{k|k-1}^{\left(j,i\right)}\right)\times \mathcal{S}\left(E_k;C_{k|k-1}^{\left(j,i\right)}\right)\hfill \end{array}$$

(45)

$${ℒ}_{\varnothing }=1-r_{k|k-1}^{j,i}{p}_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right)+r_{k|k-1}^{j,i}{p}_D\left({\xi }_{k|k-1}^{\left(j,i\right)}\right){\left(\frac{{\beta }_{k|k-1}^{\left(j,i\right)}}{{\beta }_{k|k-1}^{\left(j,i\right)}+1}\right)}^{{\alpha }_{k|k-1}^{\left(j,i\right)}}$$

(46)

There are also two parts to this formula: the first part corresponds to the case in which the target is not detected, while the second part denotes that the target does not generate a measurement.

For multimodal density, mixture reduction [35,36] can be used to reduce this to a unimodal density.

4. Results

To demonstrate the performance of the proposed algorithm, two long-wave infrared scenarios were simulated using a Satellite Tool Kit (STK) to demonstrate the performance of the proposed algorithm. The image size was 512 $\times$ 512 pixels, and there were 224 frames in total. The first scenario contains four targets, and there is no cross-trajectory. There are two targets in the second scenario. The details can be found in the following.

4.1. Scenario 1 (No-Crossing Track)

Figure 7 shows the real trajectories of the extended targets. All four extended targets appear in the first frame and disappear in frame 122, frame 138, frame 224 and frame 224. The target state ${\xi }_k=\left({\gamma }_k,x_k,E_k\right)$, where $x_k=\left(p_x,v_x,p_y,v_y\right)$, concludes the position and velocity. The initial state in which targets are born is listed in

Table 3. The initial states of the targets.

Target	State	Survival Time (Frame)
1	$\left[268.5;0;216.5;0\right]$	$\left[1,224\right]$
2	$\left[284.5;0;210.5;0\right]$	$\left[1,224\right]$
3	$\left[260.5;0;334.5;0\right]$	$\left[1,122\right]$
4	$\left[332.5;0;304.5;0\right]$	$\left[1,138\right]$

. The initial control points are defined as $M_k={\left\{\left({\rho }_{0,i},{\theta }_{0,i}\right)\right\}}_{i=1}^{12}$, where ${\rho }_{0,i}=10$ and ${\theta }_{0,i}=\frac{i-1}{6}\pi$. The sampling time is $T=1$. The detection probability is $p_D=0.95$, and the survival probability is $p_s=0.99$. The transition matrix $F$ and the process noise covariance matrix $Q$ are defined as:

$$F=I_d\otimes \left(\begin{array}{c}1 T\\ 0 1\end{array}\right),Q=q{I}_2\otimes \left(\begin{array}{c}\frac{T^3}{3} \frac{T^2}{2}\\ \frac{T^2}{2} T\end{array}\right)$$

(47)

where $\otimes$ denotes the Kronecker product, $q=0.01$. The proposed algorithm is compared to the approach based on RM and the approach based on multiple ellipses, which are denoted by GGIW-PMBM and Em-PMBM, respectively. In addition, the prior of the multiple ellipses, including the number, size and direction of the ellipses, needs to be given for the Em-PMBM-AP filter. The initial set of multiple ellipses is shown in Figure 8.

To demonstrate the effect of the measurement partition method on tracking performance, three methods are compared in this study: the amplitude-aided method (AP), the distance-based method (DP) and grid-based DBSCAN (GBDBSCAN) algorithm for which the density criterion was set to be a comparison of the mean amplitude in the search area to an amplitude threshold [23]. The modified optimal sub-pattern (mOSPA) [4,17] assignment metric was employed to assess the performance. If we suppose that the true state of an extended target is ${\xi }_k^t=\left({\gamma }_k^t,x_k^t,X_k^t\right)$ and the estimated state is ${\xi }_k^e=\left({\gamma }_k^e,x_k^e,X_k^e\right)$, the mOSPA is defined as

$$d\left({\xi }_k^t,{\xi }_k^e\right)\triangleq \frac{w_{\gamma }}{c_{\gamma }}{d}_{t,e}^{\left(c_{\gamma }\right)}+\frac{w_x}{c_x}{d}_{t,e}^{\left(c_x\right)}+\frac{w_X}{c_X}{d}_{t,e}^{\left(c_X\right)}$$

(48)

where ${\xi }_k^t$ and ${\xi }_k^e$ are the true target state and the estimated target state, respectively, and

$$w_{\gamma }+w_x+w_X=1$$

(49)

$$\{\begin{array}{c}{d}_{t,e}^{\left(c_{\gamma }\right)}=min\left(c_{\gamma },\left|{\gamma }_k^t-{\gamma }_k^e\right|\right)\hfill \\ d_{t,e}^{\left(c_x\right)}=min\left(c_x,\Vert x_k^t-x_k^e\Vert {}_2\right)\hfill \\ d_{t,e}^{\left(c_X\right)}=min(c_X,\frac{1}{M}{\displaystyle \sum }_M{\left|r\left(X_k^t\right)-r\left(X_k^e\right)\right|}^2)\hfill \end{array}$$

(50)

where $c_{\gamma }$, $c_x$ and $c_X$ were chosen to satisfy the maximum expected error for the measurement rate, kinematic state and extension state, respectively. $r\left(·\right)$ is a radial function that maps an angle to the radius of an arbitrary shape from its centroid (from 0 to $2\pi$; $M$ is the number of points that $r\left(·\right)$ was evaluated at). The details can be found in the Ref. [17].

The effectiveness of the extension estimation of the proposed algorithm (GGS-PMBM-AP), the PMBM filter based on GGIW distribution with amplitude-aided partitioning (GGIW-PMBM-AP) and the PMBM filter based on multiple ellipses with amplitude-aided partitioning (Em-PMBM-AP) are demonstrated in Figure 9 and Figure 10, which present the images of frames 36 and 90, and have been plotted using green, blue and black lines, respectively. The truth shape is shown using a red line. It can be concluded from the results that the extension, shape and orientation of the targets can be estimated better by the GGS-PMBM-AP filter even if the shape changes. This is due to the application of the B-spline model, which can accurately fit the shape by changing the position of the control points.

Figure 11 and Figure 12 show the averaged results over 100 Monte Carlo (MC) runs, corresponding to the performance metrics on mOSPA error and the number of targets. As can be seen in Figure 11a, the averaged mOSPA over 224 steps is 1.235572, 1.386461 and 1.425068 of the GGS-PMBM-AP filter, Em-PMBM-AP filter and GGIW-PMBM-AP filter, respectively. The GGS-PMBM-AP filter outperforms the Em-PMBM-AP filter and GGIW-PMBM-AP filter by 10% and 13%, respectively. This is due to the fact that the proposed B-spline model can more accurately estimate the extension even if the shape changes. Additionally, the Em-PMBM-AP filter offers better tracking accuracy than the GGIW-PMGM-AP filter due to its application of multiple ellipses. Figure 11b shows that the measurement partition method also has an impact on tracking accuracy, where the GGS-PMBM-AP filter outperforms the GGS-PMBM-DP filter and GGS-PMBM-GBDBSCAN filter. Since the amplitude threshold of the GBDBSCAN method is difficult to determine, the performance of the latter two methods is similar. Combining Figure 11b and Figure 13, it can be seen that the larger the number error, the worse the performance. Figure 12 shows the averaged number and number errors of all three filters. This error shows the difference between the estimated and true number of multiple extended targets. It can be seen that all filters are able to estimate the number of targets with minimal error.

The filters were run separately on an AMD Core 3.20 GHz CPU PC with 16 GB RAM and MATLAB R2021b. The computational complexity is illustrated by comparing the cost time. Based on 100 Monte Carlo runs, the average computational times of the filters are shown in

Table 4. Average computational times.

Filter	GGS-PMBM-AP	Em-PMBM-AP	GGIW-PMBM-AP	GGS-PMBM-DP	GGS-PMBM-GBDBSCAN
Time	10.07 s	9.95 s	9.91 s	10.34 s	9.57s

. It can be seen that the computational complexity of GGS-PMBM-DP is slightly higher than the GGS-PMBM-AP filter, proving that the amplitude-aided measurement portioning method is more effective than the method based on distance. The GGS-PMBM-GBDBSCAN algorithm has the smallest running time, which is due to the faster running speed of the clustering algorithm. In addition, the GGS-PMBM-AP filter takes more time to tackle the multiple control points. Thus, the complexity of the GGS-PMBM-AP filter is the highest.

4.2. Scenario 2 (Crossing Tracks)

There are two extended targets, the trajectories of which are shown in Figure 14. In this scenario, tracks cross at step 102. The initial state is as in

Table 5. The initial states of targets.

Target	State	Survival Time (Frame)
1	$\left[100;0;40;0\right]$	$\left[1,398\right]$
2	$\left[130;0;215.5;0\right]$	$\left[1,235\right]$

. Other parameters are the same as Scenario 1. The target shape is simulated as in Figure 8b.

Figure 15 shows the averaged mOSPA over 100 MC runs. There is a peak when the tracks cross. The difference in performance remains consistent with Scenario 1, thus demonstrating the effectiveness of the proposed algorithm.

5. Discussion

An amplitude-aided PMBM filter based on B-spline has been proposed for tracking multiple extended targets. As observed in Section 4, it can be concluded that the proposed algorithm can estimate the shape of targets more accurately compared to the other two methods. This is mainly because the control points can be adjusted in real time according to the estimated parameters, which will help to fit the real shape.

In addition, although the original intention was to estimate the shape of the target more accurately, the proposed algorithm shows better performance on the kinematic state. This can be attributed to the inherent coupling relationship between shape and motion. Because the target is captured accurately, there is an obvious improvement in the estimation performance of the kinematic state. Conversely, an accurate estimation of the target position will help to describe the shape of the target. This also explains why an accurate estimation of the shape can also improve the performance of the kinematic state. Meanwhile, the introduction of the amplitude information further contributes to the improvement of performance due to the accurate partitioning of measurement data.

6. Conclusions

A new tracking algorithm of multiple infrared non-ellipsoidal extended targets, namely GGS-PMBM, was proposed in this study to estimate the shape of the extended target. This algorithm utilizes the PMBM framework in order to implement the tracking of multiple targets. For the state of the extended target, a B-spline was employed to fit the shape by iterating control points. In addition, the measurement partition was also improved using the infrared image characteristic. Specifically, the local maximum of amplitudes in the connected domain was used as the clustering center, and clustering was performed. To verify the algorithm, a scenario using STK was simulated, and a series of infrared images were obtained. The two different scenario simulation results showed that the proposed algorithm can accurately estimate the shape of targets, and outperformed the Em-PMBM filter and the GGIW-PMBM filter. This algorithm can also be applied to group target tracking with slight modifications. In future work, it would be worthwhile to verify these findings using real data. In addition, the splitting and merging of group targets and track management issues were not addressed in this study, which are of great importance in practical scenarios.