A Multi-Model Polynomial-Based Tracking Method for Targets with Complex Maneuvering Patterns

Abstract

In the absence of a priori knowledge about target motion characteristics, the task of tracking complex maneuvering targets remains challenging. A multi-model polynomial-based complex target tracking method is presented to address this issue. Observation sequences of varying lengths are fitted by time polynomials of different orders, which are used to create a set of target motion models. Subsequently, the multi-model framework is employed to track maneuvering targets with uncertain motion characteristics. To verify the effectiveness of the suggested approach, three datasets were created with kinematic equation, the gazebo platform and real watercrafts. Based on the above three datasets, the proposed method is compared with classical multi-model methods and a deep learning method. Theoretical analysis and experimental results reveal that, in the lack of a priori information of target maneuvering features, the tracking error of the proposed method can be reduced by 12.5~30% compared with the traditional MM method. Moreover, the proposed method is able to overcome the problem of accuracy degradation caused by model misalignment and parameter tuning faced by the deep learning based methods.

1. Introduction

The investigation of accurately tracking targets with unknown maneuvering characteristics has emerged as a prominent area of sensor data processing. The primary components of this problem can be divided into two sub-problems. The first sub-problem involves estimating the parameters of the unknown system, and the estimated target motion parameters are then used to establish a set of target motion models. The second sub-problem focuses on dynamically estimating the target state based on the constructed models [1]. By acquiring a more accurate target motion model, whether it is linear or nonlinear, the current methods [2,3,4,5] can achieve precise target tracking.

The challenges of the maneuvering target tracking problem stem from the timely and accurate creation of the target motion model. Currently, the most common approaches are multi-model (MM)-based algorithms, such as the conventional Interacting Multiple Model (IMM) algorithm [6], the MIE-BLUE-IMM algorithm [7], and the Hybrid Grid Multiple Model (HGMM) algorithm [8]. These methods perform simultaneous filtering on numerous priori target motion models and then fuse the results of the multi-model filtering. Compared to traditional tracking methods that rely on maneuver detection and single-model adaptive maneuvering target tracking methods [9,10,11,12,13,14], these MM methods clearly show improved performance. However, it is important to note that these methods rely on a preset set of models that can better match the target’s motion characteristics. If the model set is inaccurate, the tracking performance will decrease. Neural networks have been employed in the field of target tracking due to their robust fitting capabilities. For instance, Back Propagation (BP) neural networks [15] and hybrid recurrent networks [16] have been utilized to learn and predict target trajectories. However, the end-to-end training approach of neural networks is not suitable for tracking complex motion model sets. Moreover, they necessitate a substantial amount of data during the training process. The DeepMTT algorithm [17] used a bidirectional Long Short Term Memory (LSTM) network to acquire knowledge of the discrepancies between the Unscented Kalman Filter (UKF) and the ground truth. The network is utilized to rectify the UKF as a secondary processing. Consequently, the approach enhanced the accuracy of the tracking process. However, the network’s input/output space was excessively large, making it challenging for the network training to reach the desired accuracy. In practical applications, it is necessary to adjust the data to make the features more similar to the training set, which is inconvenient for practical application.

In addition to deep learning methods, the current research trends are mainly data-driven methods and multi-sensor fusion methods.

On the data-driven side, [18] proposed a joint model and data-driven track segment association (JMDD-TSA) algorithm to improve the track continuity of manoeuvring targets. In order to suppress the peak error when sudden maneuvering happens for non-cooperative target tracking, [19] proposed a new data-driven maneuvering target tracking method aided with partial models. Ref. [20] proposed a hybrid-driven approach for tracking multiple highly maneuvering targets, and the time-varying constant velocity model is integrated into the Gaussian process (GP) of online learning to improve the performance of GP prediction.

In the context of multi-sensor techniques, [21] proposed a multi-sensor interactive greedy cardinalized probability hypothesis density (MS-IMM-Greedy-CPHD) filter. Ref. [22] addressed a new scheme to solve the registration problem in a distributed network. The acceleration of targets can be estimated simultaneously with higher precision. On the basis of IMM, [23] presented a robust plot-to-plot level fusion-based tracking approach, which combines data from Primary Radar (PR) and Identification of Friend or Foe (IFF) sensors. The multi-sensor multiple-model generalized labelled multi-Bernoulli filter (MS-MM-GLMB) is presented in [24], and the Gibbs-based implementation was substituted by belief-propagation-based method, which performs better in more challenging scenarios.

Time polynomials can be used to approximate the target’s state, and the complex maneuvering targets can be flexibly characterized by polynomials of different orders and coefficients [25,26]. In this case, the fixed target models is no longer required. Using temporal polynomials of varying orders in multi-model tracking algorithms can improve the dynamics and flexibility of model approximation. This paper presented a novel approach for targets tracking based on multi-model polynomials. Initially, a set of time polynomials with varying orders is formulated with constant order. The parameters of these polynomials are determined by fitting the most recent observation sequence with a specific length by the least-squares method. Subsequently, the set of polynomials serves as the target state to constitute the multi-model polynomials algorithm. Additionally, at each instant, an adaptive polynomial model set is generated based on the mentioned polynomial model set in a moment-matching way. The technique aims to provide a more accurate characterization of the model’s distribution at each moment. Finally, the process of parallel filtering is implemented based on the aforementioned model sets. The state estimations are derived by fusing the filtering result with likelihood probabilities.

It is worth mentioning that the methodology of this paper is different from the current category of methods in the literature called “polynomial Kalman filtering algorithms” [27,28,29]. The latter use polynomials to approximate the state of a nonlinear system, which is used to improve the accuracy of the nonlinear approximation.

Traditional multi-model methods, deep learning-based methods, and the proposed method show their respective strengths in different application scenarios.

Mainstream traditional multi-model algorithms [6,7,8] mainly regard the target maneuvering process as a first-order Markov process. That is, the target is considered to be transformed between a number of fixed-parameter maneuvering models, which comes with a series of problems.

Firstly, the priori model set of targets and the corresponding transfer probability matrix (TPM) are difficult to obtain accurately. When there is a mismatch between the pre-defined model set and the actual model, the performance is reduced. To solve this problem, adaptive model set design [7,8] was introduced. When the model set is difficult to cover the state space accurately, the adaptive model set is able to complement the fixed model set. In turn, the accuracy of the description of the target motion model is improved with a smaller number of preset models.

Secondly, when the target maneuvering pattern changes, the traditional multi-model approach re-converges the tracking accuracy by adaptively changing the multi-model weights. It makes the traditional multi-model methods maintain low tracking accuracy for a long time when the target has high frequency and multiple maneuvering. Part of the existing deep learning methods address this problem. For example, [17] considers that after a target maneuvering, a ‘lag’ residual is generated between the tracker and the target. The residual can be corrected by a neural network, which can improve the tracking accuracy.

Further, traditional multi-model algorithms require a known state noise covariance matrix and observation noise covariance matrix. When the above two parameter matrices are not reasonably estimated, it will make the algorithm accuracy degraded.

However, the main problem faced by deep learning methods is that the performance of the neural network greatly depends on the level of training and the match of the training data. If the test data are mismatched with the training data (e.g., factors such as noise, observation intervals, etc.), then the accuracy of the deep learning method will be significantly degraded.

The proposed method employs multi-model polynomials to estimate the target motion state. Under the condition of lacking target motion information, the method does not need to preset too many fixed-parameter motion models to obtain a relatively accurate estimation of the target motion. Moreover, the polynomial can continuously describe the target motion state when the target maneuvering pattern jumps. Hence, the tracking performance of targets with high maneuvering frequency can be improved. Compared with deep learning methods, the performance of this proposed method is not excessively dependent on parameter tuning and training.

In summary, the traditional multi-model approach is suitable for tracking targets with known maneuvering models and low maneuvering frequency. The deep learning method is suitable for tracking targets for which a large amount of observation data is available. The proposed method is suitable for tracking targets with unknown maneuvering models, high maneuvering frequency, and multiple maneuvering patterns.

Three maneuvering target trajectory datasets were constructed by the mass kinematics model, the Gazebo simulation platform and outdoor experiments, respectively. The method of this paper was compared with mainstream tracking algorithms on the three datasets. The results show that in the absence of a priori information such as target maneuver characteristics, the proposed method has smaller tracking root mean squared error (RMSE) compared with the current mainstream maneuvering target tracking algorithms, which can be reduced by 12.5~30%. The applicable scenarios of different mainstream algorithms were also analyzed and summarized.

2. Multi-Model Polynomials Tracking Algorithm

A description of the main symbols in this paper is shown in

Table 1. Description of main symbols.

$F(·)$	State transfer function	$W$	Matrix for fitting weights
$G(·)$	Control input function	$M^f$	Fixed model set
$v_k$	State noise	$f_s(n,l)$	Probability density function of mode parameters
$p_k$	Observation vectors	$F_s(n,l)$	Cumulative probability function of mode parameters
$H(·)$	Observation function	$r$	Coordinate transformation angle
$w_k$	Observation noise	$S$	mode probability space
$Q$	State noise covariance matrix	$\Theta$	Model parameter set
$R$	Observational noise covariance matrix	$F_m$	Cumulative probability of the model sample
$\beta (t)$	State vector functions for polynomials	$m$	Sample of models
$b_i$	Polynomial parameter vector	$\Sigma$	Sample variance of models
$t$	Time variable	$N(·)$	Normal distribution function
$P$	Observation Sequence Matrix	$\mu$	Probability value
$l$	Observation sequence length	$L$	Likelihood value
$n$	Polynomial order	$\left|•\right|$	Cardinal number
$B$	Polynomial parameter matrix	${\widehat{x}}_{k|k}$	State estimates
$T$	Time vector	$P_{k|k}$	Covariance matrix for state estimation
${\lambda }_{\mathrm{m}}$	Parameter of exponential distribution	$a$	Acceleration
$E(·)$	expectation operator

.

2.1. Kalman Filter for Single Polynomial Fitting

Let the equation of state of the target be:

$$x_{k+1}=F(x_k)+G(u_k)+v_k$$

(1)

where $x_k$ is the state vector of the system at the time $k$, $F$ is the state transfer function, $G$ is the controlling input function, $u_k$ is the control input at the time $k$, $v_k$ and $w_k$ are the state noise matrix at the time $k$.

The observation equation for the system is:

$$p_k=H(x_k)+w_k$$

(2)

where $p_k$ is the observation vector at time $k$, $H$ is the observation function, and $w_k$ is the observation noise. Consider $v_k$, $w_k$ as zero-mean white Gaussian sequences with covariance matrices $Q$ and $R$, respectively, which are satisfied:

$$\begin{array}{l}E(v_i{v_j}^T)=Q{\delta }_{ij}\\ E(w_i{w_j}^T)=R{\delta }_{ij}\end{array}$$

(3)

where ${\delta }_{ij}$ is the Kronecker Delta function.

Estimating the system input $G(u_k)$ in a timely and accurate manner is typically challenging when dealing with complex maneuver patterns and high maneuver frequency. As a result, the maneuvering target is directly represented as an uncertain dynamic system.

$$x_{k+1}=F_k(x_k)+v_k$$

(4)

The ability to track large-scale, high-frequency maneuvering targets is improved by making a relatively accurate approximation of $F_k$. The following polynomials are used to fit the equations of motion of the target to track the state of the maneuvering target. The following polynomial is used to fit the target motion state:

$$\beta (k)={\displaystyle \sum _{i=0}^n{b}_i{k}^i}$$

(5)

where $\beta (k)$ is the target state vector, $b_i$ is the polynomial parameter vector, $t$ is the time variable, and $n$ is the polynomial order. The target is constrained by the capacity of the dynamical system, and its acceleration is usually a continuous function of time, so its position and velocity vectors can also be described by continuously differentiable functions. This makes it reasonable to use polynomials to fit the target state function.

The coefficients of the polynomial fitting (1) are determined through estimation using the least squares criterion, utilizing an observation sequence of a specific length. Let $P\in {ℝ}^{3\times l}$ be the observation vector representing the coordinates of the target’s 3D position, and let $l$ be the length of this observation sequence. The target observation sequence is fitted using a polynomial of order $n$. The polynomial parameter matrix is denoted as $B\in {ℝ}^{3\times (n+1)}$, the time vector as $T\in {ℝ}^l$, and the fitting weight matrix as $W\in {ℝ}^{l\times l}$. They take the following forms, respectively. More fitting regression problems can be found in the literature [30,31]

$$P=\left[\begin{array}{c}{p}_{x1}\\ p_{\mathrm{y}1}\\ p_{\mathrm{z}1}\end{array}\begin{array}{c},p_{x2},\\ ,p_{\mathrm{y}2},\\ ,p_{\mathrm{z}2},\end{array}\begin{array}{c}\cdots ,\\ \cdots ,\\ \cdots ,\end{array}\begin{array}{c}{p}_{xl}\\ p_{\mathrm{y}l}\\ p_{\mathrm{z}l}\end{array}\right]$$

(6)

$$B=\left[b_1,b_2,\cdots ,b_n\right],b_i=\left[\begin{array}{c}{b}_{ix}\\ b_{iy}\\ b_{iz}\end{array}\right]$$

(7)

$$W=\left[\begin{array}{ccc}{w}_1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & w_l\end{array}\right]$$

(8)

$$T=\left[t_1,t_2\cdots t_l\right]$$

(9)

$T$ is expanded to a Vandermonde matrix $T_v\in {ℝ}^{(n+1)\times l}$:

$$T_v=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ t_1& t_2& \cdots & t_l\\ ⋮& ⋮& \ddots & ⋮\\ t_1^n& t_2^n& \cdots & t_l^n\end{array}\right]$$

(10)

The error of the fitted polynomial concerning the observed sequence is:

$$J={‖(B{T}_v-P)W‖}_F$$

(11)

Let the partial derivatives ${\displaystyle \frac{\partial (J)}{\partial (B)}}=0$, then:

$$B=PW{T_v}^T{(T_v{W}^2{T}_v)}^{-1}$$

(12)

Since $T_v{W}^2{T}_v$ is full rank when $w_i\ne 0$, it follows:

$${\displaystyle \frac{\partial (J^2)}{\partial (B)}}=2(B{T}_v-P)W^2{T_v}^T=0$$

(13)

$$B{T}_v{W}^2{T_v}^T=P{W}^2{T_v}^T$$

(14)

Equation (14) has a unique solution and (13) has the same solution as ${\displaystyle \frac{\partial (J)}{\partial (B)}}=0$. Let the polynomial fitted by the observed sequence be:

$${\beta }^{*}(k)={\displaystyle \sum _{i=0}^n{b_i}^{*}{k}^i}$$

(15)

According to (4) and (15), the next moment state of the target is estimated as:

$$x_{k+1}={\beta }^{*}(k+1)$$

(16)

(16) is brought into the Kalman filter formulation system to form a Kalman filter for single polynomial fitting.

Approximating all the states of a maneuvering target accurately with a single polynomial is challenging. The effectiveness of the approximation depends mostly on the polynomial order and the length of the fitted observation sequence. When the polynomial order is low, it becomes difficult to accurately approximate complex maneuvering targets, but the accuracy of the approximation is not significantly influenced by observation errors. Conversely, when the polynomial order is excessively high, it improves the capability to fit complex motion states, but it also becomes more susceptible to noise. In Figure 1, neglecting observation noise, the 5th- and 8th-order polynomials exhibit a better fitting effect for the trajectory compared to the 3rd-order polynomials. However, in Figure 2, when observation noise is taken into account, the high-order polynomials show a decrease in fitting effect due to their susceptibility to noise.

2.2. Polynomial Model Set Design

The parameters of the polynomial model include the polynomial order, denoted by $n$, and the length of the fitted sequence, denoted by $l$. The design of the polynomial model set is conducted in a two-dimensional integer domain and comprises two primary components: the fixed model set design and the adaptive model set design.

2.2.1. The Fixed Polynomial Model Set Design

The fixed model set $M^f$ is pre-designed based on the probability distribution space of the mode parameters, which is non-temporal-varying. Let the probability density function of the model parameters be $f_s(n,l)$ and its cumulative distribution function be $F_s(n,l)$. The model set is designed based on $f_s(n,l)$ and $F_s(n,l)$, such that the statistical properties of the model samples are consistent with the $f_s(n,l)$ and $F_s(n,l)$. Literature [32,33] proposes the design of model sets based on the probability distribution space of mode parameters and the assessment of model set quality under the condition that the distributions of mode parameters are not coupled. However, in polynomial model parameters, there is a coupling between model order and observation length. When noise is present, higher order polynomials are more sensitive to noise, and longer observation sequences can effectively suppress this property; for lower order polynomials, their fitting ability is limited, and it is not advisable to fit them with too long observation sequences, otherwise the accuracy of the approximation of the system is reduced at the current moment.

In summary, the higher-order polynomials fit long observation sequences, while the lower-order polynomials fit short observation sequences, and in this paper, we approximate that the two parameters satisfy the linear relationship:

$${\displaystyle \frac{E(nl)}{E(n)}}=r$$

(17)

Hence, before implementing the model set design, it is necessary to decouple the parameters using coordinate transformation. Let $n^{\prime },l^{\prime }$ represent the converted parameter. The relationship between the transformed parameter and the original parameter can be expressed as:

$$\left[\begin{array}{c}{n}^{\prime }\\ l^{\prime }\end{array}\right]=\left[\begin{array}{cc}\sin \psi & \cos \psi \\ \cos \psi & -\sin \psi \end{array}\right]\left[\begin{array}{c}n\\ l\end{array}\right]$$

(18)

where $\tan \psi =r$. The probability density functions of the mode parameters are satisfied:

$${f^{\prime }}_s(n^{\prime },l^{\prime })={f^{\prime }}_{s\mathrm{n}}(n^{\prime }){f^{\prime }}_{s\mathrm{l}}(l^{\prime })$$

(19)

where ${f^{\prime }}_{s\mathrm{n}}(n^{\prime }), {f^{\prime }}_{s\mathrm{l}}(l^{\prime })$ are the marginal probability density functions of each transformed parameter, respectively. The distribution of the model samples is subsequently optimized according to (19) to fit with the parameter space of the mode.

Hickernell [34] proposed using the F-centered discrepancy as a measure to assess the degree of similarity between the distribution of the model sample and the probability density of the mode parameters:

$$C{D}_{pF}(m,s)={\left[{\displaystyle \sum _{S\ne \varnothing }{\displaystyle \underset{S}{\int }{\left|F_m(\Theta )-F_s(\Theta )\right|}^p\mathrm{d}\Theta }}\right]}^{1/p}$$

(20)

where $\Theta =\{{n^{\prime }}_p,{l^{\prime }}_p\}$ is the set of mode parameters, $F(\Theta )$ denotes a point in the parameter space $S$, $F_s(\Theta )=Vol(J(F(\Theta ),a_{\Theta }))$ represents the volume of the region defined by $F(\Theta )$ and its nearest vertex. $F_m(\Theta )={\displaystyle \frac{\left|M^f\cap J(F(\Theta ),a_{\Theta })\right|}{\left|M^f\right|}}$ where $\left|•\right|$ denotes cardinal number of a set. The geometric representation of this is depicted in Figure 3.

To reduce the computer’s consumption of exponential operations while ensuring that the index function is smooth. We take $p=2$. At this point, (20) can be stated as follows:

$$\begin{array}{l}C{D}_{2F}(m,s)={\left({\displaystyle \frac{13}{12}}\right)}^{\left|\Theta \right|}-{\displaystyle \frac{2}{\left|M^f\right|}}{\displaystyle \sum _{k=1}^{\left|M^f\right|}{\displaystyle \prod _{i=1}^{\left|\Theta \right|}\left[1+\frac{1}{2}\left|F_{im}({\theta }_{ki})-\frac{1}{2}\right|-{\left|F_{im}({\theta }_{ki})-\frac{1}{2}\right|}^2\right]}}\\ +{\displaystyle \frac{1}{{\left|M^f\right|}^2}}{\displaystyle \sum _{k,l\in \left|M^f\right|}{\displaystyle \prod _{j=1}^{\left|\Theta \right|}\left[1+\frac{1}{2}\left|F_{jm}({\theta }_{kj})-\frac{1}{2}\right|+\frac{1}{2}\left|F_{jm}({\theta }_{lj})-\frac{1}{2}\right|-\frac{1}{2}\left|F_{jm}({\theta }_{kj})-F_{jm}({\theta }_{lj})\right|\right]}}\end{array}$$

(21)

The more closely the model sample set resembles the distributional characteristics of its parameters, the smaller its $C{D}_{2F}$ value is in the mode parameter space $S$.

Taking $C{D}_{2F}$ as the loss function for the optimization of the model set $M^f$, a non-gradient optimization algorithm (e.g., simulated annealing algorithm [35]) is used to optimize the model set, thus obtaining a model set with a smaller $C{D}_{2F}$ value. When the design is completed, the original model parameters are obtained through the inverse transformation of (18).

The pre-designed fixed model set $M^f$ may not be adequate to depict the target motion pattern. The tracking accuracy can be further enhanced by improving the model approximation accuracy if the model set can be adaptively built at each time based on $M^f$. The capacity of adaptive model sets to enhance the model approximation of multi-model algorithms based on the Kullback–Leibler distance has been shown and expanded upon in the literature [36,37].

2.2.2. The Adaptive Polynomial Model Set Design

Moment matching creates the adaptive model set $M_k^a$, which is based on the fixed model set $M^f$ and the adaptive model set of the preceding moment $M_{k-1}^a$ [29]. Let $\overline{m}$ and $\Sigma$ represent the first- and second-order moments of the estimated models based on $M^f$ and $M_{k-1}^a$. $M_k^a$ satisfies

$$\overline{m}={\displaystyle \sum _{i=1}^{r}P(m_k^{i\left|a\right.}|M_k^a)m_a^i}$$

(22)

$$\Sigma ={\displaystyle \sum _{i=1}^{r}P(m_k^{i\left|a\right.}|M_k^a)[{\Sigma }^i+(\overline{m}-m_k^{i\left|a\right.})}{(\overline{m}-m_k^{i\left|a\right.})}^T]$$

(23)

where $r=\left|M_k^a\right|$.

We refer to experiments in the literature [28,29,30] and find that they used the Gaussian distribution assumption, which we follow. The distribution area is divided into r regions with equal probability.

Then each region $S_i$ satisfies

$$\begin{array}{l}P(m_k^{i\left|a\right.}|M_k^a),i\ne j\\ P(m_k^{i\left|a\right.}|M_k^a)=P(m\in S_i)={\displaystyle \underset{S_i}{\int }p(m)}\mathrm{d}m\end{array}$$

(24)

The expectation of $m_a^i$ can then be calculated

$$m_a^i={\displaystyle \underset{S_i}{\int }m}p(m|m\in S_i)\mathrm{d}m={\displaystyle \frac{1}{P(m_k^{i\left|a\right.}|M_k^a)}}{\displaystyle \underset{S_i}{\int }m}p(m)\mathrm{d}m$$

(25)

Correspondingly, the estimation variance is

$$\begin{array}{c}{\Sigma }_a^i={\displaystyle \underset{S_i}{\int }(m}-\overline{m}){(m-\overline{m})}^{T}p(m|m\in S_i)\mathrm{d}m\\ \hspace{1em}\hspace{1em} ={\displaystyle \frac{1}{P(m_k^{i\left|a\right.}|M_k^a)}}{\displaystyle \underset{S_i}{\int }(m-\overline{m}){(m-\overline{m})}^T}p(m)\mathrm{d}m\end{array}$$

(26)

Consequently, the adaptive model set $M_k^a=\{m_a^1,m_a^2,\cdots ,m_a^r\}$ fulfils Equations (22) and (23). However, it is very difficult to divide regions with equal probability in two or higher dimensions of space Then, the mode vectors are projected onto a constant vector, which then degenerates into the one-dimensional Gaussian distribution $N({\alpha }^T\mathrm{m},{\alpha }^T\sum \alpha )$. The model is subsequently constructed by moment matching based on the one-dimensional Gaussian distribution. We considers $\alpha =(\sin \psi ,\cos \psi )$ as the constant vector, which is identified by our previous assumptions about the model’s correlation coefficient. Then, (24) is projected onto a single dimension.

$$P(m\in S_i)=P(d_{i-1}\le {\alpha }^{T}m\le d_i)$$

(27)

$$P(m_a^i|M_k^a)={\displaystyle \underset{d_{i-1}}{\overset{d_i}{\int }}N(m,{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )}\mathrm{d}m$$

(28)

$$m_a^i=\overline{m}+Y(\lambda -{\alpha }^T\overline{m})$$

(29)

$${\sum }^i=\sum -Y({\alpha }^T\sum \alpha -\Gamma )Y^T$$

(30)

where

$$\lambda ={\displaystyle \frac{{\alpha }^T\sum \alpha }{P(m_k^{i\left|a\right.}|M_k^a)}}\left[N(d_i,{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )-N(d_{i-1},{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )\right]+{\alpha }^T\overline{\mathrm{m}}$$

(31)

$$\begin{array}{c}\Gamma ={\displaystyle \frac{{\alpha }^T\sum \alpha }{P(m_k^{i\left|a\right.}|M_k^a)}}[(d_i+{\alpha }^T\overline{\mathrm{m}})N(d_i,{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )-\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} (d_{i-1}+{\alpha }^T\overline{\mathrm{m}})N(d_{i-1},{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )]+{({\alpha }^T\overline{\mathrm{m}})}^2+{\alpha }^T\sum \alpha +{\lambda }^2\end{array}$$

(32)

From (29)–(32), the set of adaptive models $M_k^a$ and their corresponding covariances ${\sum }^i(i=1,2\cdots r)$ can be calculated.

2.3. Multi-Model Polynomial Tracking Algorithm

The multi-model polynomial tracking algorithm is developed by the single polynomial Kalman filter described in Section 2.1 and the model set designed in Section 2.2. The state estimation ${\widehat{x}}_{k|k}$ of the target is a fusion of the estimation results under the fixed model sets $M^f$ and adaptive designed model set $M_k^a$. Let the respective estimations based on the two model sets be ${\widehat{x}}_{k|k}^f$, ${\widehat{x}}_{k|k}^a$. Then the estimation of the state is

$$\begin{array}{cc}\hfill E(x_k,M^f,M_k^a)& =E(x_k\left|M^f\right.)P(M^f\left|M^f,M_k^a,p_k\right.)+E(x_k\left|M_k^a\right.)P(M_k^a\left|M^f,M_k^a,p_k\right.)\hfill \\ \hfill {\widehat{x}}_{k|k}& ={\widehat{x}}_{k|k}^f{\mu }_k^f+{\widehat{x}}_{k|k}^a{\mu }_k^a\hfill \end{array}$$

(33)

where ${\mu }_k^f=P(M^f|M^f,M_k^a,p_k)$ and ${\mu }_k^a=P(M_k^a|M^f,M_k^a,p_k)$. Furthermore ${\widehat{x}}_{k|k}^f$ and ${\mu }_k^f$ can be computed by multi-model polynomial Kalman filtering with the fixed model set $M^f$, which can be referred to the IMM framework [6]. The state estimation process in [6] can be replaced by Equations (6)–(16). We therefore focus on the computation of ${\widehat{x}}_{k|k}^a$ and ${\mu }_k^a$.

${\widehat{x}}_{k|k}^a$ is estimated based on the adaptive model set. Therefore, ${\widehat{x}}_{k|k}^a$ is the conditional expectation of the target state under the adaptive model set:

$${\widehat{x}}_{k|k}^a=E(x_k|M_k^a)P(M_k^a\left|M_{k-1}^a,M^f,p_k\right.)$$

(34)

It should be noted that the adaptive model set for moment $k$ is generated based on the fixed model set and the adaptive model set for moment $k-1$.Therefore, (34) is rewritten as

$$\begin{array}{cc}\hfill {\widehat{x}}_{k|k}^a& ={\displaystyle \sum _{i=1}^{r}E(x_k|m_k^{i\left|a\right.},M_k^a,p_k)}P(m_k^{i\left|a\right.}|M_{k-1}^a,M^f,p_k)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^r{\widehat{x}}_{k|k}^{i\left|a\right.}{\mu }_k^{i|a}}\hfill \end{array}$$

(35)

Correspondingly, the estimated covariance is

$$P_{k|k}^a={\displaystyle \sum _{i=1}^r\left[P_{k|k}^{i|a}+({\widehat{x}}_{k|k}^{i|a}-{\widehat{x}}_{k|k}^a){({\widehat{x}}_{k|k}^{i|a}-{\widehat{x}}_{k|k}^a)}^T\right]{\mu }_k^{i|a}}$$

(36)

where ${\widehat{x}}_{k|k}^{i\left|a\right.}$ is calculated by substituting the polynomial model $m_k^{i\left|a\right.}$ into the Kalman filter formula system:

$${\widehat{x}}_{k|k}^{i\left|a\right.}=KF({\widehat{x}}_{k-1|k-1}^{i\left|a\right.},p_k,P_{k-1|k-1}^{i|a})$$

(37)

where $KF(·)$ is the polynomial Kalman filter formula system. Correspondingly, the conditional probability of ${\widehat{x}}_{k|k}^{i\left|a\right.}$ is calculated by Bayesian theory:

$$\begin{array}{cc}\hfill {\mu }_k^{i|a}& ={\displaystyle \frac{P(p_k|m_k^{i\left|a\right.},M_k^a,M^f,p_{k-1})}{P(p_k|M_k^a,p_{k-1})}}P(m_k^{i\left|a\right.}|M_k^a,M^f,p_{k-1})\hfill \\ \hfill & ={\displaystyle \frac{1}{c_a}}{L}_k^i{\mu }_{k|k-1}^{i|a}\hfill \end{array}$$

(38)

where $c_a={\displaystyle \sum _{i=1}^r{L}_k^i{\mu }_{k|k-1}^{i|a}}$, $L_k^i$ is the likelihood function of $m_a^i$:

$$L_k^i={\displaystyle \frac{1}{\sqrt{{(2\pi )}^N\det (S_k^i)}}}\left[\exp -{\displaystyle \frac{1}{2}}({v_k^i}^T{S}^{-1}{v}_k^i)\right]$$

(39)

$v_k^i$ and $S_k^i$ are the residuals and covariances obtained from polynomial Kalman filtering.

The model set probability ${\mu }_k^a$ is

$$\begin{array}{c}{\mu }_k^a=P(M_k^a|M^f,p_k)\\ \hspace{1em}\hspace{1em} ={\displaystyle \sum _{i=1}^{r}P(m_k^{i\left|a\right.}|M^f,p_k)}\end{array}$$

(40)

Which can be calculated by Bayesian theory:

$$\begin{array}{cc}\hfill {\mu }_k^a& ={\displaystyle \frac{{\displaystyle \sum _{i=1}^{r}P(p_k|m_k^{i\left|a\right.},M^f,M_{k-1}^a,p_{k-1})P(m_k^{i\left|a\right.}|M^f,M_k^a,p_{k-1})}}{P(p_k|M^f,M_k^a,p_{k-1})}}\hfill \\ \hfill & ={\displaystyle \frac{{\displaystyle \sum _{i=1}^r{L}_k^i{\mu }_{k|k-1}^i}}{P(p_k|M^f,M_k^a,p_{k-1})}}\hfill \end{array}$$

(41)

where ${\mu }_{k|k-1}^i=P(m_i^a|M^f,p_k)$ is the estimated probability of the adaptive model. Unlike fixed model sets $M^f$, adaptive model set $M_k^a$ changes over time. Therefore it is not possible to calculate ${\mu }_{k|k-1}^i$ by a constant probability transfer matrix (TPM). According to reference [7], it is assumed that there is no transfer probability between the model set $M^f$ and $M_k^a$. Additionally, ${\mu }_{k-1}^a$ and ${\mu }_{k|k-1}^a$ are assumed to be approximately equivalent. Then ${\mu }_{k|k-1}^i$ can be calculated by

$${\mu }_{k|k-1}^i={\mu }_{k|k-1}^{i|a}{\mu }_{k|k-1}^a={\mu }_{k|k-1}^{i|a}{\mu }_{k-1}^a$$

(42)

According to (40)~(42):

$${\mu }_k^a={\displaystyle \frac{{\displaystyle \sum _{i=1}^r{L}_k^i{\mu }_{k|k-1}^{i|a}{\mu }_{k-1}^a}}{P(p_k|M^f,M_k^a,p_{k-1})}}={\displaystyle \frac{c_a{\mu }_{k-1}^a}{c_f{\mu }_{k-1}^f+c_a{\mu }_{k-1}^a}}$$

(43)

where $c_f={\displaystyle \sum _{j=1}^{\left|M^f\right|}{L}_k^j{\mu }_{k|k-1}^{j|f}}$, $L_k^j$, ${\mu }_{k|k-1}^{j|f}$ are the likelihood and predicted probability of the elements in $M^f$, which can be obtained by the IMM framework [6].

The structure of the multi-model polynomial Kalman filter is shown in Figure 4 and Algorithm 1. Where MM-A denotes the polynomial multi-model filter based on $M_k^a$ and MM-F is the polynomial multi-model filter based on $M^f$, which has an algorithmic framework similar to IMM [6].

Algorithm 1. Pseudo-code of one cycle of multi-model polynomial tracking algorithm.
Cycle ofMulti-Model Polynomial Tracking Algorithm
Input:	${\mu }_{k-1}^{j\left|f\right.}$, ${\mu }_{k-1}^f$, ${\widehat{x}}_{k-1|k-1}^{j|f}$, $P_{k-1|k-1}^{j|f}$, $m^{i\left|f\right.}$, $S_{k-1}^{i|f}$, $v_{k-1}^{i|f}$, ${\mu }_{k-1}^a$, ${\mu }_{k|k-1}^{i|a}$, ${\widehat{x}}_{k-1|k-1}^{i|a}$, $P_{k-1|k-1}^{i|a}$, $S_{k-1}^{i|a}$, $v_{k-1}^{i|a}$$m_{k-1}^{i\left|a\right.}$, $P$, $\alpha$
1.	Multi-mode polynomial filtering with a fixed model set (for i = 1,…, $\left|M^f\right|$)
1.1	Calculate the fusion probability	${\mu }_{k-1}^{ij\left|f\right.}={\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}/{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}$ (where ${\pi }_{ij}$ is an element of the TPM)
1.2	Fusion state and covariance	${\overline{x}}_{k-1|k-1}^{i|f}={\displaystyle {\sum }_j}{\widehat{x}}_{k-1|k-1}^{j|f}{\mu }_{k-1}^{ji\left|f\right.}$, ${\overline{P}}_{k-1|k-1}^{i|f}={\displaystyle {\sum }_j}\left[P_{k-1|k-1}^{j|f}+\left({\widehat{x}}_{k-1|k-1}^{j|f}-{\overline{x}}_{k-1|k-1}^{i|f}\right){\left(·\right)}^T\right]{\mu }_{k-1}^{ij\left|f\right.}$
1.3	Approximating the target state with fitted polynomials	$B^i=PW{T_v}^T{(T_v{W}^2{T}_v)}^{-1}$, $m^{i\left|f\right.}={\displaystyle \sum _{j=0}^n{b_j}^{*}{k}^j}$
1.4	Kalman filtering based on fitted -polynomials	$\left[{\widehat{x}}_{k|k}^{i|f},P_{k|k}^{i|f},S_k^{i\left|f\right.},v_k^{i\left|f\right.}\right]=KF\left({\widehat{x}}_{k-1|k-1}^{i|f},P_{k-1|k-1}^{i|f},S_{k-1}^{i|f},v_{k-1}^{i|f},m^{i\left|f\right.}\right)$
1.5	Calculate the likelihood probability	$L_k^{i\left|f\right.}={\displaystyle \frac{1}{\sqrt{{(2\pi )}^N\det (S_k^{i|f})}}}\left[\exp -{\displaystyle \frac{1}{2}}\left({\left(v_k^{i|f}\right)}^T{\left(S^{i|f}\right)}^{-1}{v}_k^{i|f}\right)\right]$
1.6	Estimated fusion	${\widehat{x}}_{k|k}^f={\displaystyle \frac{{\displaystyle {\sum }_i{\widehat{x}}_{k|k}^{i|f}\left(L_k^{i\left|f\right.}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}/{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}\right)}}{{\displaystyle {\sum }_i}\left(L_k^{i\left|f\right.}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}/{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}\right)}}$ $P_{k|k}^f={\displaystyle \frac{{\displaystyle {\sum }_i\left(P_{k|k}^{i|f}+\left({\widehat{x}}_{k|k}^f-{\widehat{x}}_{k|k}^{i|f}\right){\left(·\right)}^T\right)\left(L_k^{i\left|f\right.}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}/{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}\right)}}{{\displaystyle {\sum }_i}\left(L_k^{i\left|f\right.}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}/{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}{\pi }_{ji}{\mu }_{k-1}^{j\left|f\right.}\right)}}$
2.	Adaptive model set design
2.1	Calculate the mean and covariance of the model set	$\overline{m}=\left({\displaystyle \sum {\mu }_{k-1}^{i|a}{m}_{k-1}^{i\left|a\right.}}+{\displaystyle \sum {\mu }_{k-1}^{i|f}{m}^{\left.i\right|f}}\right)$ $\Sigma ={\displaystyle {\sum }_i\left({\sum }_{k-1}^i+\left(\overline{m}-m_{k-1}^{i\left|a\right.}\right){\left(·\right)}^T\right){\mu }_{k-1}^{i|a}}+{\displaystyle {\sum }_i\left(\overline{m}-m^{i\left|f\right.}\right){\left(·\right)}^T{\mu }_{k-1}^{i|f}}$
2.2	Model design based on moment matching (for i = 1,…, $\left|M^a\right|$)	$P(m\in S_i)=P(d_{i-1}\le {\alpha }^{T}m\le d_i)$ $\lambda ={\displaystyle \frac{{\alpha }^T\sum \alpha }{P(m_a^i|M_k^a)}}\left[N(d_i,{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )-\right.\left.N(d_{i-1},{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )\right]+{\alpha }^T\overline{\mathrm{m}}$ $\begin{array}{c}\Gamma ={\displaystyle \frac{{\alpha }^T\sum \alpha }{P(m_a^i|M_k^a)}}[(d_i+{\alpha }^T\overline{\mathrm{m}})N(d_i,{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )-\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} (d_{i-1}+{\alpha }^T\overline{\mathrm{m}})N(d_{i-1},{\alpha }^T\overline{\mathrm{m}},{\alpha }^T\sum \alpha )]+{({\alpha }^T\overline{\mathrm{m}})}^2+{\alpha }^T\sum \alpha +{\lambda }^2\end{array}$ $m_k^{i\left|a\right.}=\overline{m}+Y(\lambda -{\alpha }^T\overline{m})$ ${\sum }_k^i=\sum -Y({\alpha }^T\sum \alpha -\Gamma )Y^T$
3.	Multi-mode polynomial filtering with a adaptive designed model set (for i = 1,…, $\left|M^a\right|$)
3.1	Approximating the target state with fitted polynomials	$B^i=PW{T_v}^T{(T_v{W}^2{T}_v)}^{-1}$, $m_k^{i\left|a\right.}={\displaystyle \sum _{j=0}^n{b_j}^{*}{k}^j}$
3.2	Kalman filtering based on fitted polynomials	$\left[{\widehat{x}}_{k|k}^{i|a},P_{k|k}^{i|a},S_k^{i\left|a\right.},v_k^{i\left|a\right.}\right]=KF\left({\widehat{x}}_{k-1|k-1}^{i|a},P_{k-1|k-1}^{i|a},S_{k-1}^{i|a},v_{k-1}^{i|a},m_k^{i\left|a\right.}\right)$
3.3	Calculate the likelihood probability	$L_k^{i\left|a\right.}={\displaystyle \frac{1}{\sqrt{{(2\pi )}^N\det (S_k^{i|a})}}}\left[\exp -{\displaystyle \frac{1}{2}}\left({\left(v_k^{i|a}\right)}^T{\left(S^{i|a}\right)}^{-1}{v}_k^{i|a}\right)\right]$
3.4	Estimated fusion	${\widehat{x}}_{k|k}^a={\displaystyle {\sum }_i{\widehat{x}}_{k|k}^{i\left|a\right.}{\mu }_k^{i|a}}$ $P_{k|k}^a={\displaystyle \sum _{i=1}^r\left[P_{k|k}^i+({\widehat{x}}_{k|k}^{i\left|a\right.}-{\widehat{x}}_{k|k}^a){({\widehat{x}}_{k|k}^{i\left|a\right.}-{\widehat{x}}_{k|k}^a)}^T\right]{\mu }_k^{i|a}}$ ${\mu }_k^{i|a}={\displaystyle \frac{L_k^i{\mu }_{k|k-1}^{i|a}}{{\displaystyle {\sum }_i{L}_k^i{\mu }_{k|k-1}^{i|a}}}}$
3.5	Probability normalization	${\mu }_k^a={\displaystyle \frac{{\displaystyle {\sum }_i{L}_k^i{\mu }_{k|k-1}^{i|a}}{\mu }_{k-1}^a}{{\displaystyle {\sum }_i{L}_k^j{\mu }_{k|k-1}^{j|f}{\mu }_{k-1}^f}+{\displaystyle {\sum }_i{L}_k^i{\mu }_{k|k-1}^{i|a}}{\mu }_{k-1}^a}}$ ${\mu }_k^f=1-{\mu }_k^a$
4	Filter result fusion	${\widehat{x}}_{k|k}={\mu }_k^a{\widehat{x}}_{k|k}^a+{\mu }_k^f{\widehat{x}}_{k|k}^f$ $P_{k|k}={\mu }_k^a\left[P_{k|k}^a+\left({\widehat{x}}_{k|k}-{\widehat{x}}_{k|k}^a\right)\left(·\right)\right]+{\mu }_k^f\left[P_{k|k}^f+\left({\widehat{x}}_{k|k}-{\widehat{x}}_{k|k}^f\right)\left(·\right)\right]$

Remark 1. 

We provide further explanation of the choice of polynomial order and the length of fitted observation sequence. Let the state trajectory of the target over a period of time $[t_i,t_{i+1}]$ be $x(t)$ ($x(t)$ is second-order differentiable). Assume that the target is sampled k times at $[t_i,t_{i+1}]$ . At each sampling point, the state of the target is $x(k)$, and the observation noise is $w(k)$, which is considered as Gaussian white noise. Then the observations satisfy

$$p(k)=x(k)+w(k)$$

(44)

When the noise is ignored, According to the Weierstrass first approximation theorem [38]:

For any $x(t)$ continuous in $[t_i,t_{i+1}]$, there exist polynomials $\beta$ and $\epsilon >0$, such that $\left|\beta -x\right|<\epsilon$. Where $\beta$ is of the form:

$$\begin{array}{cc}\hfill \beta (x,t)& ={\displaystyle \sum _{i=0}^{m}x(t_i+\frac{k(t_{i+1}-t_i)}{n})\left(\begin{array}{c}n\\ k\end{array}\right)}{x}^k{(1-k)}^{n-k}\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^n{b}_i{t}^i}\hfill \end{array}$$

(45)

Therefore, when n is sufficiently large, there exists an optimal polynomial ${\beta }^{*}(t)$ satisfying ${\displaystyle {\int }_{t_i}^{t_{i+1}}\left|{\beta }^{*}(t)-x(t)\right|\mathrm{dt}\to 0}$. As shown in Figure 1.

When system noise is taken into account, it is assumed that there exists a virtual continuous function $p(t)$:

$$p(k)=x(k)+w(k)$$

(46)

When n is sufficiently large, $\beta (t)$ satisfies ${\displaystyle {\int }_{t_i}^{t_{i+1}}\left|{\beta }^{*}(t)-p(t)\right|\mathrm{dt}\to 0}$. Then $\beta (t)$ will deviate from the true target state. As shown in Figure 2. Thus in the polynomial approximation process, we can model the optimization as follows:

$$\begin{array}{l}\underset{b_i,l}{\min } E\left({\left(\beta (k+1)-x(k+1)\right)}^2\left|p(k-l),\cdots ,p(k)\right.\right)\\ s.t.\\ \beta (k)={\displaystyle \sum _{i=0}^n{b}_i{k}^i}\\ {\displaystyle \frac{\partial \left({\displaystyle \sum _k{\left({\displaystyle \sum _{i=0}^n{b}_i{k}^i}-p(k)\right)}^2}\right)}{\partial b_i}}=0\\ p(k)=x(k)+w(k)\end{array}$$

(47)

This optimization problem can be solved by a multi-model framework. That is, the target motion state is fitted by presetting several polynomials with different observation lengths and orders. Finally, the optimal estimate of the target motion model is obtained based on likelihood probability fusion.

Table 2. Scenarios of application of the polynomial model.

	High Polynomial Order	Low Polynomial Order
Long fitted sequence	Low noise level	—
Short fitted sequence	—	High noise level

summarizes the scenarios for which polynomial order and fit length are applicable.

Remark 2. 

We determine the applicability of the proposed methods compared to traditional MM methods by quantitative methods.

For the estimation of state $x(t)$ and observation $p(k)$ at time $k+1$, the inference form of IMM is:

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}^{IMM}& =E\left(x_{k+1}\left|m_1^f,\cdots ,m_n^f,p(1),\right.\cdots p(n)\right)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^{n}E\left(x_{k+1}\left|m_i^f,p(1),\right.\cdots p(n)\right)P\left(m_i^f\left|m_1^f,\cdots ,m_n^f,p(1),\cdots ,p(n)\right.\right)}\hfill \end{array}$$

(48)

where $m_i^f$ is the i th model in the model set. The inference form of HGMM is:

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}^{HGMM}& =E\left(x_{k+1}\left|m_1^f,\cdots ,m_n^f,m_1^a,\cdots ,m_r^f,p(1),\right.\cdots p(n)\right)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^{n}E\left(x_{k+1}\left|m_i^f,p(1),\right.\cdots p(n)\right)P\left(m_i^f\left|m_1^f,\cdots ,m_n^f,m_1^a,\cdots ,m_r^f,p(1),\cdots ,p(n)\right.\right)}\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^{r}E\left(x_{k+1}\left|m_i^a,p(1),\right.\cdots p(n)\right)P\left(m_i^a\left|m_1^f,\cdots ,m_n^f,m_1^a,\cdots ,m_r^f,p(1),\cdots ,p(n)\right.\right)}\hfill \end{array}$$

(49)

where $m_i^a$ is the i th model in the adaptive designed model set. For the proposed method

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}^{PRO}& =E\left(x_{k+1}\left|{\beta }_1^f,\cdots ,{\beta }_n^f,{\beta }_1^a,\cdots ,{\beta }_r^f,p(1),\right.\cdots p(n)\right)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^{n}E\left(x_{k+1}\left|{\beta }_i^f,p(1),\right.\cdots p(n)\right)P\left({\beta }_i^f\left|{\beta }_1^f,\cdots ,{\beta }_n^f,{\beta }_1^a,\cdots ,{\beta }_r^a,p(1),\cdots ,p(n)\right.\right)}\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^{r}E\left(x_{k+1}\left|{\beta }_i^a,p(1),\right.\cdots p(n)\right)P\left({\beta }_i^a\left|{\beta }_1^f,\cdots ,{\beta }_n^f,{\beta }_1^a,\cdots ,{\beta }_r^a,p(1),\cdots ,p(n)\right.\right)}\hfill \end{array}$$

(50)

where ${\beta }_1^f$ and ${\beta }_i^a$ are fixed and adaptive polynomial models, respectively.

The three methods have a formally similar estimated variance calculation

$$\begin{array}{l}{P}_k={\displaystyle \sum _i\left(P_k^i+\left({\widehat{x}}_k^i-{\widehat{x}}_k\right){\left({\widehat{x}}_k^i-{\widehat{x}}_k\right)}^T\right)}P\left(m_i\left|M,p(1),\cdots ,p(n)\right.\right)\\ {\displaystyle \sum _{i}P\left(m_i\left|M,p(1),\cdots ,p(n)\right.\right)}=1\end{array}$$

(51)

where $P_k^i$, ${\widehat{x}}_k^i$ are the variance and expectation predicted under the i th model and $M$ is the model set. According to the optimality criterion, when $P_k$ is smaller, its estimation is more effective. Therefore, the covariance matrices of different methods can be output, and the method with a smaller covariance matrix norm can be chosen as the optimal method. For example, $‖P‖=\sup (‖Px‖,‖x‖=1)$, as shown in the Figure 5.

3. Maneuvering Target Trajectory Datasets Construction

Three maneuvering target trajectory datasets were created for a more comprehensive algorithm performance validation. The first is the kinematic random maneuver trajectory dataset, which consists of trajectories simulated by kinematic equations for random maneuver. The second is the UAV maneuver trajectory dataset, which is based on the Gazebo simulation platform, and the trajectories are generated by manually manipulating a rotary-wing UAV. The third is the surface boat trajectory dataset, which is generated by measuring the trajectories of surface boats.

In the kinematic random maneuver trajectory dataset, the target maneuver parameters follow a known probability distribution. This dataset has a large number of trajectories, which allows for more effective neural network training. The trajectory data are used to evaluate the performance of various types of tracking algorithms in the presence of a substantial amount of trajectory data and adequate prior knowledge about the target maneuver.

In contrast to the kinematic random maneuver trajectory dataset, the UAV maneuver trajectory dataset possesses the following attributes: Firstly, target trajectories are dynamically constrained, making them more realistic. Second, there is limited previous knowledge about the maneuver features of the target, and the target’s maneuver is intricate and no longer follows a known probability distribution. This dataset is utilized to evaluate the tracking algorithm’s performance in scenarios where there is an abundance of data, but little a priori information regarding the target’s movement.

Compared to the first two datasets, the surface speedboat dataset exhibits a rather uncomplicated target maneuver pattern, and its maneuver features are also constrained by dynamics. The dataset contains a relatively modest amount of measured trajectory data due to limitations in the conditions of the experiment. The trajectory data are utilized to validate the tracking algorithm’s performance in scenarios with limited previous knowledge about the target’s maneuver, uncomplicated maneuver patterns, and a short dataset.

3.1. Kinematic Random Maneuver Trajectory Dataset

The differential equation of motion of the mass is established based on Equation (1) with a sampling period of 0.01s. In the direction of motion, the acceleration obeys a uniform distribution $U(a_{x\min },a_{x\max })$, in which the probability of both acceleration and deceleration is 0.5; In the horizontal plane perpendicular to the direction of motion, the acceleration obeys the uniform distribution $U(a_{y\min },a_{y\max })$. The probability that the direction of acceleration is to the left and the right is 0.5; The acceleration in the vertical direction is uniformly distributed as $U[a_{z\min },a_{z\max }]$, with positive values representing upward acceleration and negative values representing downward acceleration. The duration period of the acceleration follows an exponential distribution with parameters ${\lambda }_{\mathrm{m}}$. In addition, the drag factor $a_{x-}$ is added to prevent the target speed from being too large:

$$a_{x-}=-k_a{v}^2$$

(52)

where $v$ is the velocity rate of the target and $a_{x-}$ is the opposite direction of the target’s motion. The maximum velocity rate of the target is controlled by adjusting $k_a$.

The targets are categorized into three groups based on their acceleration distribution: strong maneuverability targets, medium maneuverability targets, and weak maneuverability targets. These categories are determined by the following parameters ($a$ has the measure of ${\mathrm{m}/\mathrm{s}}^2$):

(1) Strong maneuverability targets: maximum velocity rate: 30 m/s, $[a_{x\min },a_{x\max }]=[6,12]$, $[a_{y\min },a_{y\max }]=[6,12]$, $[a_{z\min },a_{z\max }]=[-4,4]$, ${\lambda }_{\mathrm{m}}=5 \mathrm{s}$.

(2) Medium maneuverability targets: maximum velocity rate: 20 m/s, $[a_{x\min },a_{x\max }]=[3,7]$, $[a_{y\min },a_{y\max }]=[3,7]$, $[a_{z\min },a_{z\max }]=[-2,2]$, ${\lambda }_{\mathrm{m}}=15 \mathrm{s}$.

(3) Weak maneuverability targets: maximum velocity rate: 10 m/s, $[a_{x\min },a_{x\max }]=[1,3]$, $[a_{y\min },a_{y\max }]=[1,3]$, $[a_{z\min },a_{z\max }]=[-1,1]$, ${\lambda }_{\mathrm{m}}=30 \mathrm{s}$.

Part of the generated maneuver target trajectory is shown in Figure 6.

3.2. UAV Maneuver Trajectory Dataset

The dataset is constructed using the Gazebo simulation tool, which offers accurate and realistic physics simulations. Utilizing high-fidelity sensor streams, this system enables many forms of robotic sensor simulations and allows for rapid iterative simulations in relatively authentic physical contexts. Thus, the data generated by the UAV in this platform are more accurate and realistic.

The rotor UAV is initially modelled with the Iris UAV data package supplied by PX4-AutoPilot. The UAV is a quadcopter model that shares the shape features of a quadcopter UAV. It is capable of transmitting data to and from the ground station via the Mavlink V2.0 software, as depicted in Figure 7. The accurate location of the UAV target is determined by disseminating gazebo/model states messages from the UAV.

Subsequently, the LiDAR sensor is simulated (as shown in Figure 8). The simulation model of the Velodyne multi-line LIDAR is created by the Gazebo simulation platform. The laser sensor’s functionality is implemented by specifying the .sdf file, as depicted in Figure 7. Upon detecting the UAV, the lidar will transmit the point cloud data of the UAV via the cloud points message topic, as seen in Figure 9b. The design parameters of the LiDAR include the following:

The detecting range is 400 m, the blind zone at the boundary is 2 m, the angular resolution in both horizontal and vertical directions is 0.2°, there are 64 laser lines, the scanning speed is 10 r/s, and the error parameter can be adjusted.

Ultimately, the point cloud data were standardized and refined to acquire the recorded trajectory of the unmanned aerial vehicle (UAV). Various maneuver modes were created by directing the UAV to perform turning, folding motion, and hovering, respectively. The trajectory data were then gathered, as depicted in Figure 10.

3.3. Surface Boat Trajectory Dataset

The steps to construct the dataset are as follows:

Firstly, Real-time kinematic (RTK) was placed on the boat, which has a centimetre-level positioning error and a sampling period of 0.05 s, and its measured position data were taken as the true value, as shown in Figure 11. The speedboat was then manoeuvred to accelerate, decelerate, straighten and “s-steer” on the lake. Due to the limitation of the water area, the recorded track data were in a single direction, and finally, the track data in each direction were obtained by rotating and turning operations, as shown in Figure 12.

4. Maneuvering Target Tracking Experiment

To verify the algorithm of this paper more comprehensively, tracking experiments are carried out based on the above three trajectory datasets respectively.

4.1. Target Tracking Experiment Based on Kinematic Random Maneuver Trajectory Dataset

Based on the dataset, the method of this paper is compared with IMM and HGMM (the traditional multi-mode algorithms), and DeepMTT (the deep learning algorithm) [17] respectively. Among them, DeepMTT was trained based on pytorch under the noise condition: azimuth:0.01 rad/ pitch: 0.01 rad/ distance: 15 m. The trajectory is sliced into trajectory segments of length 20, and then input into the neural network for training, the training epoch is 5000, the batchsize is set to 120, the batchsize is set to 60 after 2000 epochs of training, the loss function is the MSE function, and the optimizer is the SSD, and the learning rate is 10–3.

The model set of the traditional multi-model algorithm is set according to the motion law of the maneuver target in 3.1.1 respectively, so that it is aligned with the motion model set of the target as much as possible. The model set parameters are:

Strong maneuver target: [0, 0, 0], [−9, 0, 0], [9, 0, 0], [0, 9, 0], [0, −9, 0], [0, 0, 2], [0, 0, −2].

Medium maneuver targets: [0, 0, 0], [−5, 0, 0], [5, 0, 0], [0, 5, 0], [0, −5, 0], [0, 0, 1], [0, 0, −1].

Weak maneuver targets: [0, 0, 0], [−2, 0, 0], [2, 0, 0], [0, 2, 0], [0, −2, 0], [0, 0, 0.5], [0, 0, −0.5].

Further, the state noise matrix and observation noise matrix are adjusted iteratively to obtain better tracking accuracy. The tracking trajectories of the algorithms are shown in Figure 13. The tracking root mean square error (RMSE) results of the algorithms under different motion models and error conditions are shown in Table 2,

Table 3. Characteristics of the datasets.

Dataset	Reality	Maneuvering Information	Data Volume	Purpose
Kinematic	low	adequate	large	to validate the performance of various tracking methods when more comprehensive target motion information is available.
UAV	medium	medium	medium	to validate the ability of different methods to track the typical maneuvering trajectories of airborne targets in the absence of maneuvering information.
Surface Boat	high	deficient	small	to validate the performance of various tracking methods under the conditions of lack of target maneuvering information and lack of training data.

and

Table 4. RMSE of strong maneuverability targets (observation distance: 4 km).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.001 rad, 0.001 rad, 15 m	7.1034	6.5548	5.7835	6.5855
0.003 rad, 0.003 rad, 25 m	10.8474	10.0600	10.1511	10.4058
0.006 rad, 0.006 rad, 35 m	14.6329	13.7920	15.1262	15.6328

.

From Table 4,

Table 5. RMSE of medium maneuverability targets (observation distance: 4 km).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.001 rad, 0.001 rad, 15 m	6.4875	5.9776	5.4925	6.3487
0.003 rad, 0.003 rad, 25 m	9.4398	8.8448	9.6851	9.8823
0.006 rad, 0.006 rad, 35 m	14.1523	12.9611	13.8346	14.4856

and

Table 6. RMSE of weak maneuverability trajectories (observation distance: 4 km).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.001 rad, 0.001 rad, 15 m	5.9382	5.9053	5.1974	5.0989
0.003 rad, 0.003 rad, 25 m	9.8312	9.7763	9.4889	9.3737
0.006 rad, 0.006 rad, 35 m	12.7814	12.7350	14.7790	14.1203

, the following conclusions can be drawn.

(1)

Under the condition that the model set and noise parameters are set reasonably, the accuracy of the traditional multi-model algorithm is better than the proposed algorithm as a whole, and as the target maneuverability is stronger and the observation error is larger, the advantage is more obvious. As shown in Table 4, its tracking accuracy can be reduced by 13% compared with this paper’s algorithm under large error conditions. However, as the maneuvering strength and error decrease, the advantage of the traditional method is also weakening, and even in Table 6, its tracking accuracy is lower than that of this paper’s algorithm.

(2)

Under the condition of sufficient data volume, the deep learning algorithm can also show better tracking accuracy, however, with the increase of observation error, the accuracy of the algorithm decreases more obviously, as in Table 4, under the condition of small error, the tracking accuracy is in the optimal among the algorithms, and with the increase of error, its tracking accuracy is comparable with this paper’s algorithm.

(3)

In summary, when the target motion characteristics and system noise parameters are known, and the amount of data is sufficient, the traditional multi-model algorithm and the deep learning algorithm are better than the proposed algorithm as a whole.

4.2. Target Tracking Experiment Based on UAV Maneuver Trajectory Dataset

Target tracking experiments are conducted based on the dataset constructed in Section 3.1. The proposed algorithm is compared with IMM, HGMM and DeepMTT on turn maneuver, zigzagging maneuver and motion-hover maneuver trajectories respectively.

Under this experimental condition, the traditional multi-model algorithms are unable to establish a more accurate model set based on the motion characteristics of the target. At the same time, to simulate the actual situation, the system noise parameters of the traditional multi-mode algorithm are no longer adjusted as follows:

The model set is [0, 0, 0], [−5, 0, 0], [5, 0, 0], [0, 5, 0], [0, −5, 0], [0, 0, 1], [0, 0, −1], and the state noise matrix is 5 × I_6×6. The observation noise matrices are set to be 4 × I_3×3, 16 × I_3×3, 40 × I_3×3, and 100 × I_3×3, respectively.

The deep learning method is trained in the same way as in Section 4.1 under the noise condition: azimuth: 0.01 rad/pitch: 0.01 rad/distance: 0.5 m. The tracking curves of the different algorithms under different motion models and error conditions are shown in Figure 14, and the RMSE of the tracking algorithms are shown in

Table 7. Tracking RMSE for turning maneuver (observation distance 200 m).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.01 rad, 0.01 rad, 0.5 m	0.8167	0.8158	0.7094	0.8072
0.03 rad, 0.03 rad, 1.5 m	2.2619	2.2608	1.7638	1.6830
0.05 rad, 0.05 rad, 3.5 m	4.0122	4.0110	3.1262	2.9188
0.07 rad, 0.07 rad, 6.0 m	5.9267	5.9217	4.4565	3.9794

,

Table 8. Tracking RMSE for zigzagging maneuver (observation distance 200 m).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.01 rad, 0.01 rad, 0.5 m	0.9783	0.9802	0.7573	0.8146
0.03 rad, 0.03 rad, 1.5 m	2.8797	2.8766	2.2122	2.1440
0.05 rad, 0.05 rad, 3.5 m	5.1574	5.1518	3.9141	3.6573
0.07 rad, 0.07 rad, 6.0 m	7.2456	7.2540	5.4097	4.8844

and

Table 9. Tracking RMSE for motion-hover maneuver (observation distance 200 m).

Error Condition (Azimuth/Pitch/Distance)	IMM	HGMM	DeepMTT	Proposed
0.01 rad, 0.01 rad, 0.5 m	0.8656	0.8660	0.6854	0.8147
0.03 rad, 0.03 rad, 1.5 m	2.6387	2.6374	2.1030	2.0879
0.05 rad, 0.05 rad, 3.5 m	4.9555	4.9522	3.8335	3.5979
0.07 rad, 0.07 rad, 6.0 m	6.4603	6.4530	4.8063	4.4234

.

From Figure 14 and Table 7, Table 8 and Table 9, the following conclusions can be drawn.

(1) The tracking accuracy of the proposed algorithm is better than that of the traditional multi-model algorithm when there is insufficient a priori information about the target motion pattern. As the observation error increases, the tracking accuracy of the proposed algorithm improves more obviously, and the RMSE of the proposed algorithm can be reduced by 33% at most compared with the traditional multi-model algorithm. It indicates that the proposed algorithm has a higher performance when tracking targets that have less a priori motion information.

(2) Based on the dataset, the accuracy of the deep learning algorithm is comparable to that of the proposed algorithm. However, as the observation error increases, the tracking accuracy of the proposed algorithm is better than that of the deep learning algorithm.

The deep learning method was trained with the observation noise of azimuth: 0.01 rad/pitch: 0.01 rad/distance: 0.5 m. Thus, when the noisy condition of the test data becomes large, its data features are mismatched with the training set, which, in turn, deteriorates the performance of the neural network. The proposed method fits the target state with a polynomial of different observation lengths and orders. Firstly, the fitting polynomial is able to give a minimum mean square estimate of the target state under noisy conditions. Therefore, the fitted polynomials are better able to overcome the target state estimation problem in the absence of an a priori motion model of the target. Secondly, the proposed method employs multi-model polynomial parallel filtering and then fuses the filtering results based on likelihood probabilities. In comparison to deep learning-based methods, multiple polynomial models in a predefined model set are better adapted to noise than networks trained for specific noise conditions. When the observation noise increases, the accuracy of the proposed method is less affected than that of the deep learning-based method.

4.3. Target Tracking Experiment Based on Surface Boat Trajectory Dataset

Based on the surface boat trajectory dataset, the proposed algorithm and the comparison algorithm are verified for the tracking capability of maneuvering targets under real conditions. Three sections of speedboat trajectories, namely, slow steering, straight-line motion and steering acceleration, are intercepted and compared under different error conditions. The parameters of the comparison algorithms are set as follows:

The model sets of IMM and HGMM algorithms are designed as [0, 0, 0], [−9, 0, 0], [9, 0, 0], [0, 9, 0], [0, −9, 0]. The state noise matrix is 10 × I_6×6 The observation noise matrix is set to 120 × I_3×3, 400 × I_3×3 and 900 × I_3×3 based on the observation error, respectively. The training method of the deep learning method is the same as Section 3.2.

The experimental results are shown in Figure 15 and

Table 10. RMSE of slow steering.

Error Condition (Azimuth/Distance)	IMM	HGMM	DeepMTT	Proposed
0.003 rad, 6 m	4.8891	4.8521	6.0199	5.3103
0.007 rad, 14 m	9.2031	8.9553	13.6923	8.3379
0.010 rad, 20 m	12.9969	12.3169	20.4941	10.95414

,

Table 11. RMSE of straight-line motion.

Error Condition (Azimuth/Distance)	IMM	HGMM	DeepMTT	Proposed
0.003 rad, 6 m	4.5609	4.5585	6.1798	6.3821
0.007 rad, 14 m	9.2031	8.9553	15.80843	8.7176
0.010 rad, 20 m	13.2518	12.0387	23.7324	12.3036

and

Table 12. RMSE of steering acceleration.

Error Condition (Azimuth/Distance)	IMM	HGMM	DeepMTT	Proposed
0.003 rad, 6 m	5.6798	5.6723	7.8045	4.6865
0.007 rad, 14 m	9.6845	9.6685	17.6552	9.4866
0.010 rad, 20 m	11.5615	11.5356	25.5932	13.6073

.

The following conclusions can be drawn from Figure 14 and Table 10, Table 11 and Table 12.

(1) The deep learning algorithm is not adequately trained due to fewer data in the dataset, which, in turn, makes its tracking accuracy significantly lower than other algorithms. In contrast, the traditional multi-mode method and the proposed method maintain a high tracking accuracy as a whole, which reflects that the proposed method and the traditional multi-mode method do not rely on data training and have better practicality.

(2) Under the conditions of this experiment, the tracking accuracies of the traditional multi-model algorithm and the proposed algorithm are comparable. During linear motion (shown in Table 12), the tracking accuracy of the traditional multi-model algorithm is slightly higher than that of the proposed algorithm, and during curvilinear motion (shown in Table 10 and Table 12), the tracking accuracy of the proposed algorithm is slightly higher than that of the multi-model algorithm as a whole.

(3) In summary, the proposed algorithm and the multi-model algorithm are more suitable for tracking tasks with insufficient data, when the tracking accuracy of the traditional multi-model algorithm depends on the a priori maneuver characteristics of the target. The proposed algorithm is more suitable when there is a lack of a priori in-formation about the maneuver characteristics.

(4) It should be further noted that, in the above three experiments, we knew the observation noise characteristics. At the same time, we made the performance of the traditional method as good as possible by continuously adjusting the state noise matrix. However, in practice, it is often difficult to accurately estimate the observation noise and the state noise, so the performance of the traditional method may further deteriorate.

4.4. Comparison of Algorithm’s Computational Time Consumption

The experimental hardware is as follows: CPU: Intel(R) Xeon(R) Silver 4210; GPU: NVIDIA Quadro 4000; Memory: 8 GB, where IMM, HGMM, and the proposed algorithm are computed with CPU; the DeepMTT algorithm is computed with GPU. The software environments are MATLAB 2022b and Pytorch 1.11. The former was used to run traditional multi-model methods and the latter was used to run the deep learning method.

The IMM and HGMM algorithms presuppose seven models, and the proposed algorithm presupposes five polynomial models. The average computation time for the four algorithms is shown in

Table 13. Computational time consumption of different algorithms.

Algorithms	IMM	HGMM	DeepMTT	Proposed
Time (ms)	0.39064	0.7031	1.7587	0.7813

.

As shown in Table 13, the IMM algorithm takes the shortest time, and the proposed algorithm and the HGMM algorithm take a comparable amount of time, and both of them can satisfy the real-time target tracking task with a sampling frequency of less than 1200 HZ. The deep learning algorithm has the heaviest computational burden and can only support the target tracking task below 600 HZ sampling frequency.

5. Conclusions

A multi-model polynomials Kalman filtering algorithm is proposed for the problem of high maneuver frequency and large maneuver scale of targets, which makes it difficult to establish the target motion model in time and effectively. Firstly, time polynomials of multiple orders are used to fit observation sequences of different lengths for estimating target motion states. Secondly, under the multi-model framework, the parallel interactive filtering of multi-mode polynomials is achieved by presetting a fixed model set as well as generating an adaptive model set, which in turn improves the state estimation capability of complex maneuvering targets. The fixed model set is generated by simulated annealing with L2-centered-discrepancy as the optimization objective based on the distribution space of the model parameters, while the adaptive model set is generated by moment matching, which describes more accurately the distribution characteristics of the model set at each moment compared to the fixed model set.

To validate the algorithm, three types of maneuver target tracking datasets are constructed. Based on the above datasets, the proposed algorithm is compared and analyzed with the traditional multi-model algorithms IMM and HGMM, as well as the deep learning-based tracking algorithm DeepMTT. The results show that the algorithm in this paper and the two types of algorithms compared are suitable for different task scenarios, respectively. When the amount of data is large and the computing resources are sufficient, the deep learning-based tracking method can show better tracking accuracy. When there is more a priori information about the target motion pattern, the traditional multi-model algorithm can better perform the target tracking task. This method is suitable for scenarios with less data volume, and unknown and complex target motion patterns, and can somewhat solve the problem of insufficient data volume faced by deep learning and the problem of model set misalignment faced by traditional multi-model algorithms. Scenarios of application of the three types of algorithms are shown in Figure 16.

However, limited by the research conditions and other factors, this paper has some shortcomings in terms of experiments. Firstly, the collected data mainly come from simulation and semi-physical simulation. The construction and tracking experiments of manoeuvre target datasets in more realistic and complex scenarios have not yet been carried out. Secondly, the performance of the algorithm under other complex conditions (e.g., sensor failure, environmental factors) has not been comprehensively verified.

In our subsequent research, we will construct manoeuvre target datasets in more realistic and complex scenarios. For the algorithms, we will verify and improve the performance of the algorithms in complex environments. At the same time, the data-driven aspect will also be investigated. We will study how to use a large amount of observation data to improve the approximation accuracy of the target states and then enhance the tracking capability of the algorithm for complex manoeuvre targets.