**1. Introduction**

State estimation is an important topic in maneuvering target tracking of traditional ships and surface autonomous ships. As a minimum mean square error estimator, the Kalman filter (KF) provides the optimal state estimation method for a linear Gaussian system (The abbreviations are summarized in Table 1). However, the system nonlinearity and system model uncertainty in jump Markov systems (JMSs) may lead to the KF's performances degrading dramatically. Additionally, the KF assumes that the measurement noises obey Gaussian distributions and that all measurements need to arrive in time. Both aspects may affect the filtering accuracy significantly when these assumptions are not satisfied. Thus, the extension for the KF under various assumptions has attracted considerable attention on account of its widely applied in engineering, for example, targets tracking, signal processing, integrated navigation, fault diagnosis, ect. [1–6].

The surface maneuvering target tracking is a problem of state estimation in JMSs. Since the state estimation in JMSs is well-known computational intractability and nondeterministic polynomial difficult, it is hard to find an optimal solution [7]. In the past few decades, some sub-optimal solutions were presented, for example, the interacting multiple model (IMM) approach [8], pseudo-Bayesian technique [9], particle filter [10], etc. Among these methods, the IMM approach is one of the most efficient since it reasonably balances the estimation performances and the computation complexity. In the process of IMM technique

**Citation:** Chen, C.; Zhou, W.; Gao, L. A Novel Robust IMM Filtering Method for Surface-Maneuvering Target Tracking with Random Measurement Delay. *J. Mar. Sci. Eng.* **2023**, *11*, 1047. https://doi.org/ 10.3390/jmse11051047

Academic Editor: Mohamed Benbouzid

Received: 17 April 2023 Revised: 9 May 2023 Accepted: 13 May 2023 Published: 14 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

development, many researchers contributed remarkable achievements. The stability and performances of the IMM filter are analyzed by [11,12]. A series of executable pseudocode is provided in [13]. IMM algorithms are widely used in various practical engineering [14–16].

The typical IMM filters parallelly perform a bank of KFs to estimate mode-conditional system state and probability at each cycle. The weighted sum of sub-filter outputs obtains the final state estimates [17]. An effective method to develop the performances of sub-filters (KFs) is the variational Bayesian (VB) inference, which implements approximations for the conjugate exponential model with acceptable computational costs [18]. Based on the combination of VB theory and the KF, the development of the KF has made significant progress and remarkable achievements [19–22]. Extended to IMM filters for the state estimation in JMSs, some VB-based IMM methods were designed in past few years. Ref. [23] presented an IMM filter to adaptively estimate the unknown process and measurement noise covariance matrices in JMS. The authors selected the conjugate prior distribution of noise covariance matrix as the inverse-Wishart distribution, the conditional system state vectors and the noise parameters can be estimated simultaneously by VB inference. This method achieved remarkable estimation accuracy when the noise covariances were unknown. However, this technique requires linear Gaussian systems, and the measurements need to be arrived in time, which is usually not satisfied in actual tracking scenarios affected by measurement outliers and communication channel latency.


**Table 1.** Definitions of notations in this paper.

Aiming to solve state estimation problems in JMSs with measurement outliers, Ref. [24] modeled the measurement noises as students' t-distributions (STDs), and estimated the system states and mode probabilities jointly by VB inference. Compared with the typical IMM algorithm, the estimation performance was significantly developed. Ref. [25] developed the IMM filter in [24] by reasonably selecting the conjugate prior distributions of covariances and degree of freedom parameters as inverse Wishart and Gamma distributions, respectively. The state vectors and the degree of freedom parameters were inferred simultaneously by VB method. This algorithm improved the estimation accuracy and can be utilized for nonlinear JMSs. The methods mentioned in [24,25] show superior estimation accuracy in the scenarios of state estimation in JMSs containing heavy-tailed measurement noises (HTMNs). However, the performances of these approaches may lose efficacy when the random measurement delay happens since they assume the measurements can be obtained in real-time. To deal with the one-step random measurement delay (OSRMD), Ref. [26] designed a particle filter by defining two variables to extend the original system state vector. This filter can effectively handle the systems with delayed measurements. Although this algorithm provided a solution to measurement delay and achieved satisfactory estimation accuracy, the particle filters have problems of high computational burden and curse of dimensionality. Ref. [27] proposed a Gaussian approximate filter with OSRMD, and measurement noise vectors are inferred by the Bayesian rules. Ref. [28] developed

the KF to deal with the OSRMD in nonlinear systems. In this method, the authors utilized a random Bernoulli variable (RBV) to describe the random measurement delay, and the system states were inferred by the VB technique. However, these methods failed to deal with the HTMN and system model jumping problems in JMSs. In complex JMSs engineering applications, the HTMNs and random measurement delay often exist simultaneously. It is necessary to propose a more general IMM filter for HTMNs coexist with random measurement delay.

This article designs a novel robust IMM filter to deal with the filtering problems in JMSs with HTMN and OSRMD. The HTMN and the OSRMD in JMSs are reasonably modeled in the designed filter. By transforming the measurement likelihood function to a new form, this algorithm constructed a new hierarchical Gaussian state space model (HGSSM). According to the VB inference, the state vectors, model probability, RBVs, and distribution parameters are inferred simultaneously. Target tracking experiment validates that the presented approach has better performance than existing IMM approaches.

The main contribution of our work is summarized as follows:


The rest parts of this paper can be arranged as the following sections. Section 2 provides the problem statement. We construct a new HGSSM in Section 3. Furtherly, we summarize the derivations of our designed IMM method in Section 4. The performances of exsiting IMM algorithms and the proposed are verified in Section 5. Finally, the conclusions of this paper is in Section 6.

#### **2. Problem Statement**

Considering a state-space model of nonlinear jump Markov system:

$$\mathbf{x}\_s = f\_{s-1}(\mathbf{x}\_{s-1}, M\_{s-1}) + \mathbf{g}\_{s-1}, \mathbf{s} \ge 1 \tag{1}$$

$$z\_s = h\_s(\mathbf{x}\_s, M\_s) + \varepsilon\_{s\_s} \mathbf{s} \ge 1 \tag{2}$$

$$y\_s = \begin{cases} \ (1 - \sigma\_s)z\_s + \sigma\_s z\_{s-1}, & s > 1 \\ z\_{s\prime} & s = 1 \end{cases} \tag{3}$$

where *<sup>s</sup>* is time index, *xs* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is the system state vector, *zs* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* and *ys* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* denote the ideal and real measurement vectors, *m* and *n* are their dimension numbers. *fs*−1(·) and *hs*(·) refer to the process function and measurement function. We denote *fs*−1(*xs*−1, *Ms*−1) as *f <sup>i</sup> <sup>s</sup>*−1(*xs*−1) and *hs*(*xs*, *Ms*) as *<sup>h</sup><sup>i</sup> <sup>s</sup>*(*xs*) conditioned on *Ms* = *i*, respectively. The system mode *Ms* indicates the state of the Markov chains. It selects a value from a limited set

{1, 2, 3, . . . , *<sup>k</sup>*} according to the transition probability matrix *<sup>T</sup>* = ' *λji*( *k*×*k* , which satisfies the equations as follows:

$$
\lambda\_{ji} \stackrel{\Delta}{=} P\{M\_{\mathfrak{s}} = i | M\_{\mathfrak{s}} = j\} \tag{4}
$$

$$\sum\_{i=1}^{k} \lambda\_{ji} = 1\tag{5}$$

*gs*−<sup>1</sup> represents the Gaussian process noise vector with nominal covariance matrix *Qs*, while *ε<sup>s</sup>* refers to the HTMN vector caused by measurement outliers with nominal covariance matrix *Rs*. *ys* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* denotes the real random delayed measurement vector. The RBV *σ<sup>s</sup>* ∈ {0, 1} with corresponding probability can be defined as follows:

$$p(\sigma\_s = 1) = \varphi\_s \tag{6}$$

$$p(\sigma\_s = 0) = 1 - \varphi\_s \tag{7}$$

where *ϕ<sup>s</sup>* ∈ (0, 1) denotes the probability of OSRMD. Note that random variables *xs*, *σs*, *gs*−1, and *rs* are independent of each others.

Utilizing Equation (3) and Equations (6) and (7), the measurement likelihood PDF is as follows:

$$\begin{aligned} p(y\_s|\mathbf{x}\_s, \mathbf{x}\_{s-1}, M\_s = i) &= \sum\_{\sigma\_s=0}^{1} p(y\_s, \sigma\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, M\_s = i) \\ p(y\_s = 1) p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 1, M\_s = i) &+ p(\sigma\_s = 0) p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 0, M\_s = i) \\ p = \mathfrak{g}\_s p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 1, M\_s = i) &+ (1 - \mathfrak{g}\_s) p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 0, M\_s = i) \end{aligned} \tag{8}$$

Utilizing Equations (2) and (3) yields

$$p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 1, M\_s = i) = p\_{\varepsilon\_{s-1}}[y\_s - h\_{s-1}(\mathbf{x}\_{s-1})] \tag{9}$$

$$p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s = 0, M\_s = i) = p\_{\varepsilon\_s}[y\_s - h\_s(\mathbf{x}\_s)] \tag{10}$$

In Equation (8), replaced by Equations (9) and (10), we can obtain:

$$p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, M\_\mathbf{s} = i) = \varphi\_\mathbf{s} p\_{\mathbf{f}\_{s-1}}[y\_s - h\_{s-1}(\mathbf{x}\_{s-1})] + (1 - \varphi\_\mathbf{s}) p\_{\mathbf{f}\_\mathbf{s}}[y\_\mathbf{s} - h\_\mathbf{s}(\mathbf{x}\_\mathbf{s})] \tag{11}$$

However, the measurement likelihood PDF in Equation (11), which is formulated by a weighted sum form, has neither conjugate property nor closed property. Therefore, it is not easy to introduce the VB technique directly. Aiming at this difficulty, in the next section, the weighted sum form of measurement likelihood PDF is converted to exponential multiplication, which can be further utilized to design a novel adaptive estimation method for JMSs.

#### **3. Construction of the HGSSM**

*3.1. Measurement Likelihood PDF Convertion*

The relationship between RBV *σ<sup>s</sup>* with fixed delay probability can be seen in Equations (6) and (7). The probability mass function *p*(*σs*|*ϕ<sup>s</sup>* ) is defined as the following formula:

$$p(\sigma\_{\mathfrak{s}}) = (1 - \varphi\_{\mathfrak{s}})^{(1 - \sigma\_{\mathfrak{s}})} \varphi\_{\mathfrak{s}}^{\sigma\_{\mathfrak{s}}} \tag{12}$$

according to Equations (11) and (12), reformulating Equation (8) by:

$$p\left(y\_s|\mathbf{x}\_s, \mathbf{x}\_{s-1}, M\_s = i\right) = \sum\_{\sigma\_i=0}^{1} p\left(y\_s, \sigma\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, M\_s = i\right)$$

$$\begin{split} p &= \sum\_{\sigma\_i=0}^{1} \left(\varphi\_s\right)^{\sigma\_i} p\_{\varepsilon\_{s-1}} (y\_s - h\_{s-1}(\mathbf{x}\_{s-1}))^{\sigma\_i} (1 - q\_s)^{(1 - \sigma\_i)} p\_{\varepsilon\_s} (y\_s - h\_s(\mathbf{x}\_s))^{(1 - \sigma\_i)} \\ &= \sum\_{\sigma\_i=0}^{1} p(\sigma\_s) p\_{\varepsilon\_{s-1}} (y\_s - h\_s(\mathbf{x}\_s))^{\sigma\_i} p\_{\varepsilon\_s} (y\_s)^{(1 - \sigma\_i)} \end{split} \tag{13}$$

then, transforming the conditional measurement likelihood PDF to exponential product form:

$$p(y\_s | \mathbf{x}\_s, \mathbf{x}\_{s-1}, \sigma\_s) = \left[ p\_{\mathbf{f}\_{s-1}}(y\_s - h\_{s-1}(\mathbf{x}\_{s-1})) \right]^{\sigma\_s} \left[ p\_{\mathbf{f}\_s}(y\_s - h\_s(\mathbf{x}\_s)) \right]^{(1-\sigma\_s)} \tag{14}$$

Equation (14) can be also denoted by the formula as follows:

$$p(y\_s | \eta\_{s\cdot}, \sigma\_s) = \left[p\_{\mathbb{E}\_{s-1}}(y\_s - h\_{s-1}(\mathbf{x}\_{s-1}))\right]^{\sigma\_s} \left[p\_{\mathbb{E}\_s}(y\_s - h\_s(\mathbf{x}\_s))\right]^{(1-\sigma\_s)}\tag{15}$$

where the extended state vector *η<sup>s</sup>* = ' *x<sup>T</sup> <sup>s</sup> x<sup>T</sup> s*−1 (*T* .

#### *3.2. Prior PDFs Selection*

First of all, the one-step predictive PDF of the extended system state *η<sup>s</sup>* is formulated by

$$p(\eta\_s | y\_{1:s-1}) = N(\eta\_s; \mathfrak{f}\_{s|s-1}, P^{\eta \eta}\_{s|s-1}) \tag{16}$$

where *<sup>η</sup>*ˆ*s*|*s*−<sup>1</sup> represents the predicted mean vector, *<sup>P</sup>ηη <sup>s</sup>*|*s*−<sup>1</sup> refers to the predicted covariance matrix, and they are expressed by:

$$\boldsymbol{\mathfrak{H}}\_{s|s-1} = \left[ \underset{\boldsymbol{r}}{\mathfrak{A}}\_{s|s-1}^{T} \quad \mathfrak{A}\_{s-1|s-1}^{T} \right]^{T} \tag{17}$$

$$P\_{s|s-1}^{\eta\eta} = \begin{bmatrix} P\_{s|s-1} & P\_{s-1,s|s-1} \\ \left(P\_{s-1,s|s-1}\right)^T & P\_{s-1|s-1} \end{bmatrix} \tag{18}$$

where *<sup>η</sup>*ˆ*s*|*s*−<sup>1</sup> , *Ps*|*s*−<sup>1</sup> , and *Ps*−1,*s*|*s*−<sup>1</sup> can be calculated by the standard Gaussian approximate filter [29], i.e.,

$$\hat{\mathbf{x}}\_{s|s-1} = \int f\_{s-1}^{i}(\mathbf{x}\_{s-1}) \mathcal{N}(\mathbf{x}\_{s-1}; \hat{\mathbf{x}}\_{s-1|s-1}, P\_{s-1|s-1}) d\mathbf{x}\_{s-1} \tag{19}$$

$$P\_{s|s-1} = \int f\_{s-1}^{i}(\mathbf{x}\_{s-1}) f\_{s-1}^{i(T)}(\mathbf{x}\_{s-1}) N\left(\mathbf{x}\_{s-1}; \mathbf{f}\_{s-1|s-1}, P\_{s-1|s-1}\right) d\mathbf{x}\_{s-1} \tag{20}$$

$$P\_{\mathbf{s}-1,\mathbf{s}|\mathbf{s}-1} = \int \mathbf{x}\_{\mathbf{s}-1} f\_{\mathbf{s}-1}^{i(T)}(\mathbf{x}\_{\mathbf{s}-1}) \mathbf{N} \Big(\mathbf{x}\_{\mathbf{s}-1}; \hat{\mathbf{x}}\_{\mathbf{s}-1|\mathbf{s}-1}, P\_{\mathbf{s}-1|\mathbf{s}-1}\Big) d\mathbf{x}\_{\mathbf{s}-1} - \hat{\mathbf{x}}\_{\mathbf{s}-1|\mathbf{s}-1} \hat{\mathbf{x}}\_{\mathbf{s}|\mathbf{s}-1}^{T} \tag{21}$$

Secondly, in the aspects of processing the measurement noises, since the STD has heavier tails compared with Gaussian distribution and is more robust to outlier [30–32], STD is selected to model HTMN, i.e.,

$$p(\varepsilon\_{s-1}) = St(\varepsilon\_{s-1}; 0, R\_{s-1}, \tau\_{s-1}) \tag{22}$$

$$p(\varepsilon\_s) = St(\varepsilon\_s; 0, R\_s, \tau\_s) \tag{23}$$

where the *p*(*εs*−1) and *p*(*εs*) represent the PDFs of previous and current measurement noise vector. Based on the properties of the STDs, PDFs *p*(*εs*−1) and *p*(*εs*) are formulated in hierarchical form as follows:

$$p(\varepsilon\_{s-1}) \int \mathcal{N}(\varepsilon\_{s-1}; 0, R\_{s-1}/v\_{s-1}) p(v\_{s-1}) dv\_{s-1} \tag{24}$$

$$p(v\_{s-1}) = G\left(v\_{s-1}; \frac{a\_{s-1}}{2}, \frac{a\_{s-1}}{2}\right) \tag{25}$$

$$p(\varepsilon\_s) \int N(\varepsilon\_{\rm s}; 0, R\_{\rm s}/v\_{\rm s}) p(v\_{\rm s}) dv\_{\rm s} \tag{26}$$

$$p(v\_s) = G\left(v\_s; \frac{a\_s}{2}, \frac{a\_s}{2}\right) \tag{27}$$

Utilizing Equations (15), (24)–(27), the conditional likelihood PDF *p*(*ys*|*xs*, *xs*−1, *vs*, *vs*−1, *σs*) can be obtained by:

$$p(y\_s | \mathbf{x}\_{\delta}, \mathbf{x}\_{s-1}, v\_{\delta}, v\_{s-1}, \sigma\_{\delta}) = \left[\mathcal{N}(y\_{\delta}; h\_{\delta}(\mathbf{x}\_{\delta}), R\_{\delta}/v\_{\delta})\right]^{(1-\sigma\_{\delta})} \left[\mathcal{N}(y\_{s-1}; h\_{s-1}(\mathbf{x}\_{s-1}), R\_{\delta-1}/v\_{s-1})\right]^{\sigma\_{\delta}} \tag{28}$$

Finally, according to Equation (28), the measurement vector *ys* depends on the extended state vector *ηs*, auxiliary variables *vs* and *vs*−1, and the RBV *σs*, the joint prior PDFs *p*(*ηs*, *vs*, *vs*−1, *σs*|*y*1:*s*−<sup>1</sup> ) can be obtained as follows:

$$\begin{cases} \begin{aligned} p(\boldsymbol{\Xi}|\boldsymbol{y}\_{1:s-1}) &= p(\boldsymbol{\eta}\_{s}|\boldsymbol{y}\_{1:s-1})p(\boldsymbol{v}\_{s})p(\boldsymbol{v}\_{s-1})p(\boldsymbol{v}\_{s})\\ &= N\Big(\boldsymbol{\eta}\_{s};\boldsymbol{\eta}\_{s\mid s-1},\boldsymbol{P}\_{s\mid s-1}^{\eta\eta}\big)G\big(\boldsymbol{v}\_{s};\frac{\boldsymbol{a}\_{1}}{2},\frac{\boldsymbol{a}\_{s}}{2}\big)G\big(\boldsymbol{v}\_{s-1};\frac{\boldsymbol{a}\_{s-1}}{2},\frac{\boldsymbol{a}\_{s-1}}{2}\big)(1-\boldsymbol{q}\_{s}\big)^{(1-\sigma\_{s})}\boldsymbol{q}\_{s}^{\sigma\_{s}}\\ \boldsymbol{\Xi} \overset{\Lambda}{=} \{\boldsymbol{\eta}\_{s},\boldsymbol{v}\_{s},\boldsymbol{v}\_{s-1},\sigma\_{s}\} \end{aligned} \end{cases} \tag{29}$$

Then, the HGSSM consists Equations (15), (17)–(21) and (24)–(29) is constructed, based on which a novel adaptive estimation algorithm will be designed to infer the state vector *η<sup>s</sup>* and BRV *σ<sup>s</sup>* in JMSs with OSRMD and HTMN.

#### **4. Design of the Proposed Filter**

Based on the IMM filtering framework, the designed filter as an iterative algorithm mainly consists of four recursive steps, i.e., interacting/mixing, mode-conditioned filtering, mode probability updating, and combination.

**Step 1:** Interacting/Mixing

To interact the state vector and noise parameters, the mixing probability needs to be computed first utilizing the mode transition probability *λji* as follows:

$$
\mu\_{s-1}^{ji} = \frac{1}{\mathcal{E}\_i} p\_{ji} \mu\_{s-1}^j \tag{30}
$$

where the mode transition probability has been defined in Equations (4) and (5), and the normalized constant *e*¯*<sup>i</sup>* is defined by Equation (31):

$$
\bar{e}\_i = \sum\_{j=1}^k \lambda\_{ji} \mu\_{s-1}^j \tag{31}
$$

where *μ<sup>j</sup> <sup>s</sup>*−<sup>1</sup> represents the mode probability, and the mode index *<sup>j</sup>*, *<sup>i</sup>* <sup>∈</sup> {1, 2, 3, . . . , *<sup>k</sup>*}. Suppose that the posterior PDF at *j*-th mode can be obtained by

$$p(\eta\_{s-1}|y\_{1:s-1}, M\_{s-1} = j) = N(\eta\_{s-1}; \mathfrak{h}\_{s-1|s-1}^j, P\_{s-1|s-1}^j) \tag{32}$$

Utilizing the total probability thorem and the Bayes rules to calculate the mixing PDFs *p*(*ηs*|*y*1:*s*−1, *Ms* = *i*), i.e.,

$$p(\eta\_s | y\_{1:s-1}, M\_s = i) = \sum\_{j=1}^k p(\eta\_s | y\_{1:s-1}, M\_{s-1} = j) p(M\_{s-1} = j | M\_s = i, y\_{1:s-1})$$

$$= \sum\_{j=1}^k \mu\_{s-1}^{ji} N\left(\eta\_{s-1}; \eta\_{s-1}^j, P\_{s-1}^j\right) \tag{33}$$

Using the Kullback-Leibler average fusion method [33], the summation term in Equation (33) is approximated as follows:

$$p(\eta\_{s-1}|y\_{1:s-1}, M\_s = i) \approx N\left(\eta\_{s-1}; \mathfrak{h}\_{s-1\mid s-1}^{0i}, P\_{s-1\mid s-1}^{0i}\right) \tag{34}$$

where the mixing mean vector and covariance matrix of the distribution of *ηs*−<sup>1</sup> can be obtained as follows:

$$\begin{cases} \begin{aligned} \eta\_{s-1|s-1}^{\text{iif}} &= \sum\_{j=1}^{k} \mu\_{s-1|s-1}^{\text{iif}} \eta\_{s-1|s-1}^{j} \\\ P\_{s-1|s-1}^{\text{p\eta}, \text{iif}} &= \sum\_{j=1}^{k} \mu\_{s-1|s-1}^{\text{iif}} \left( P\_{s-1|s-1}^{\text{p\eta}, j} + \left( \mathfrak{x}\_{s-1|s-1}^{j} - \mathfrak{x}\_{s-1|s-1}^{\text{0i}} \right) \left( \mathfrak{x}\_{s-1|s-1}^{j} - \mathfrak{x}\_{s-1|s-1}^{\text{0i}} \right)^{T} \right) \end{aligned} \tag{35} \end{cases} \tag{36}$$

**Step 2:** Mode-Conditioned Filtering

The posterior PDFs at the *i*-th mode is formulated by:

$$p(\eta\_{s}, v\_{s}, v\_{s-1}, \sigma\_{s} | y\_{1:s}, M\_{\mathfrak{s}} = i) = \frac{p(\eta\_{s}, v\_{s}, v\_{s-1}, \sigma\_{\mathfrak{s}} | y\_{\mathfrak{s}}, M\_{\mathfrak{s}} = i)}{p(y\_{\mathfrak{s}} | y\_{1:s-1}, M\_{\mathfrak{s}} = i)}\tag{36}$$

Since the presence of couplings between the state vectors and the noise covariances, the analytic solution can not be obtained for the mode conditional posterior PDFs *p*(*ηs*, *vs*, *vs*−1, *σs*|*ys*, *Ms* = *i*). However, by using the VB inference, the approximate solution can be available. Based on the VB technique, the abovementioned posterior PDFs is approximated as a multiplied factor form:

$$p(\eta\_{s\cdot\nu}v\_{s\cdot\nu}v\_{s-1}, \sigma\_s | y\_{1:s\cdot\nu}M\_{\mathfrak{s}} = i) \approx q\left(\eta\_{s}^{i}\right)q\left(\upsilon\_{s}^{i}\right)q\left(\upsilon\_{s-1}^{i}\right)q\left(\sigma\_{s}^{i}\right) \tag{37}$$

where *q*(·) is approximated posterior PDF of *p*(·). Using the VB technique to minimize the kullback-Leibler divergence (KLD) between the actual and approximated posterior PDFs, the optimal approximation PDFs can be obtained, i.e.,

$$\begin{aligned} &q\left(\eta\_s^i\right)q\left(\upsilon\_s^i\right)q\left(\upsilon\_{s-1}^i\right)q\left(\sigma\_s^i\right) \\ &=\arg\min KLD\left[q\left(\eta\_s^i\right)q\left(\upsilon\_s^i\right)q\left(\upsilon\_{s-1}^i\right)q\left(\sigma\_s^i\right)\right]\left|p\left(\eta\_s^i,\upsilon\_s^i,\upsilon\_{s-1}^i,\sigma\_s^i\middle| y\_{1:s\prime}M\_s=i\right)\right] \end{aligned} \tag{38}$$

where *KLD*(*q*(·)*p*(·) ) <sup>Δ</sup> <sup>=</sup> ) *<sup>q</sup>*(·)ln *<sup>q</sup>*(*y*) *<sup>p</sup>*(*y*) *dy*, Equation (38) is the optimal method for approximation:

$$\ln q(\phi\_s) = E\_{\Xi\_s^{-\Phi\*}}[\ln p(\Xi|y\_{1:s}) | M\_s = i] + \alpha t\_{\Phi\*} \tag{39}$$

where *φ<sup>s</sup>* denotes one element from Ξ, Ξ−*φ<sup>s</sup> <sup>s</sup>* represents the rest elements after removing element *φ<sup>s</sup>* in Ξ*s*. *cstφ<sup>s</sup>* refers to the constant about *φs*.

The joint PDF *p*(Ξ, *y*1:*s*|*Ms* = *i*) in Equation (39) can be formulated by

$$p(\Xi, y\_{1:s} | M\_{\sf s} = i) = p(y\_{1:s-1} | M\_{\sf s} = i) p(y\_s | \eta\_{\sf s}, v\_{\sf s}, v\_{\sf s-1}, \sigma\_{\sf s}, M\_{\sf s} = i) p(\eta\_{\sf s} | y\_{1:s-1}, M\_{\sf s} = i)$$

$$p(v\_{\sf s} | M\_{\sf s} = i) p(v\_{\sf s-1} | M\_{\sf s} = i) p(\sigma\_{\sf s} | M\_{\sf s} = i) \tag{40}$$

Substituting Equations (12) and (28), (29) into (40):

$$p(\Xi, y\_{1:s} | M\_s = i) = \left[ N\left( y\_s^i, h\_s^i \left( \mathbf{x}\_s^i \right), R\_s^i / v\_s^i \right) \right]^{(1 - \sigma\_s^i)} \left[ N\left( y\_{s-1}^i, h\_{s-1}^i \left( \mathbf{x}\_{s-1}^i \right), R\_{s-1}^i / v\_{s-1}^i \right) \right]^{(1 - \sigma\_s^i)}$$

$$\cdot N\left( \eta\_{1:s}^i \eta\_{1:s-1}^i, P\_{s \mid s-1}^{\eta \eta, i} \right) G\left( v\_{s:}^i, \frac{d\_s^i}{2}, \frac{d\_s^i}{2} \right) G\left( v\_{s-1}^i, \frac{d\_{s-1}^i}{2}, \frac{d\_{s-1}^i}{2} \right) \left( 1 - q\_s^i \right)^{(1 - \sigma\_s^i)} $$

$$\cdot \left( q\_s^i \right)^{\sigma\_s^i} p(y\_{1:s-1}) \tag{41}$$

Utilizing Equation (41), ln *p*(Ξ, *y*1:*s*|*Ms* = *i*) is further derived:

$$\begin{split} \ln p(\Xi, y\_{1:s}|M\_{5} = i) &= 0.5 \Big( m \Big( 1 - \sigma\_{s}^{i} \Big) + a\_{s}^{i} - 2 \Big) \ln \sigma\_{s}^{i} + 0.5 \Big( m \upsilon\_{s}^{i} + a\_{s-1}^{i} - 2 \Big) \ln \upsilon\_{s-1}^{i} \\ &- 0.5 \Big( 1 - \sigma\_{s}^{i} \Big) \upsilon\_{s}^{i} \Big( y\_{s}^{i} - h\_{s}^{i} \Big( \mathbf{x}\_{s}^{i} \Big) \Big)^{T} R\_{s}^{i(-1)} \Big( y\_{s}^{i} - h\_{s}^{i} \Big( \mathbf{x}\_{s}^{i} \Big) \Big) - 0.5 \sigma\_{s}^{i} \upsilon\_{s-1}^{i} \Big( y\_{s}^{i} - h\_{s-1}^{i} \Big( \mathbf{x}\_{s-1}^{i} \Big) \Big)^{T} \\ &- \sigma\_{s-1}^{i(-1)} \Big( y\_{s}^{i} - h\_{s-1}^{i} \Big( \mathbf{x}\_{s-1}^{i} \Big) \Big) - 0.5 \Big( \eta\_{s}^{i} - \eta\_{s|s-1}^{i} \Big)^{T} R\_{s|s-1}^{\eta \eta \prime - 1} \Big( \eta\_{s}^{i} - \eta\_{s|s-1}^{i} \Big) + \sigma\_{s}^{i} \ln \eta\_{s}^{i} \\ &+ \left( 1 - \sigma\_{s}^{i} \right) \ln \Big( 1 - \sigma\_{s}^{i} \Big) - 0.5 \upsilon\_{s-1}^{i} - 0.5 \upsilon\_{s}^{i} + \text{cst}\_{\eta,\nu,\sigma} \end{split} \tag{42}$$

The approximated PDFs togather with relevant expectations are obtained after *N* iterations:

$$q^{(N)}\left(\eta\_s^i\right) \approx N\left(\eta\_s^i; \left.\eta\_{s\mid s}^{(N)i}{}\_{\prime}P\_{s\mid s}^{(N)\eta\eta, i}\right)\right.\tag{43}$$

$$q^{(N)}\left(v\_s^i\right) \approx G\left(v\_s^i; a\_{s\mid s}^{(N)i}, \beta\_{s\mid s}^{(N)i}\right) \tag{44}$$

$$q^{(N)}\left(v\_{s-1}^{i}\right) \approx G\left(v\_{s-1}^{i}; \mathfrak{a}\_{s-1\mid s-1}^{(N)i}, \mathfrak{b}\_{s-1\mid s-1}^{(N)i}\right) \tag{45}$$

$$E^{(N)}\left[v\_s^i\right] = \frac{a\_s^{(N)i}}{\beta\_s^{(N)i}}\tag{46}$$

$$E^{(N)}\left[v\_{s-1}^i\right] = \frac{\alpha\_{s-1}^{(N)i}}{\beta\_{s-1}^{(N)i}}\tag{47}$$

$$E^{(N)}\left[\ln \upsilon\_s^i\right] = \psi\left(a\_s^{(N)i}\right) - \ln\left(\beta\_s^{(N)i}\right)\_{\cdot \cdot} \tag{48}$$

$$E^{(N)}\left[\ln v\_{s-1}^{i}\right] = \psi\left(\alpha\_{s-1}^{(N)i}\right) - \ln\left(\beta\_{s-1}^{(N)i}\right) \tag{49}$$

$$E^{(N)}\left[\sigma\_s^i\right] = \frac{\Pr^{(N)}\left(\sigma\_s^i = 1\right)}{\Pr^{(N)}\left(\sigma\_s^i = 1\right) + \Pr^{(N)}\left(\sigma\_s^i = 0\right)}\tag{50}$$

$$E^{(N)}\left[1-\sigma\_s^i\right] = 1-E^{(N)}\left[\sigma\_s^i\right] \tag{51}$$

where *ψ*(·) refers to digamma function, and the derivation details of Equations (43)–(51) can be found in Appendix A.

**Step 3:** Mode Probability Updating

In this step, the mode probability Pr{*Ms* = *i*|*y*1:*<sup>s</sup>* } can be updated by combining the total probability equation with the Bayes rule, i.e.,

$$\mu\_s^i = \Pr\{M\_s = i | y\_{1:s} \} = \frac{p(y\_s | M\_s = i, y\_{1:s-1}) \Pr(M\_s = i | y\_{1:s-1})}{\sum\_{j=1}^k p(y\_s | M\_s = j, y\_{1:s-1}) \Pr(M\_s = j | y\_{1:s-1})} = \frac{\bar{\mu}\_{s-1}^i \Gamma\_s^i}{\sum\_{j=1}^k \bar{\mu}\_{s-1}^i \Gamma\_s^i} \tag{52}$$

where

$$
\mu\_{s-1}^i = \sum\_{j=1}^k \lambda\_{ji} \mu\_{s-1}^j \tag{53}
$$

Γ*i <sup>s</sup>* refers to the likelihood function, and the ln *p*(*ys*|*y*1:*s*−1, *Ms* = *i*) is rewritten by

$$\ln p(y\_s | y\_{1:s-1}, M\_s = i) = L(\Lambda) + KL(\Lambda \| \, p(\Xi | y\_{1:s}, M\_s = i) \,) \tag{54}$$

In Equation (54), *L*(Λ) represents the evidence low bound of the function Λ(·). Utilizing the VB technique, the term of *KLD*(·) can be minimized and approximated to zero, therefore, we can obtain:

$$
\Gamma^i\_s \cong \exp\{L(\Lambda)\}\tag{55}
$$

where

$$L(\Lambda) = \ln \frac{p(y\_{1:s}, \eta\_{s\cdot s}, v\_{s-1}, \sigma\_s | M\_s = i)}{p(\eta\_{s\cdot}, v\_{s\cdot}, v\_{s-1}, \sigma\_s | y\_{1:s}, M\_s = i)} - \ln p(y\_{1:s-1} | M\_s = i) \tag{56}$$

and

$$p(\eta\_{s}, \upsilon\_{s}, \upsilon\_{s-1}, \sigma\_{s} | y\_{1:s}, M\_{s} = i) \simeq q\left(\eta\_{s}^{i}\right)q\left(\upsilon\_{s}^{i}\right)q\left(\iota\_{s-1}^{i}\right)q\left(\sigma\_{s}^{i}\right) \tag{57}$$

After approximate operation, the *L*(Λ) can be newly derived:

$$L(\Lambda) = E\_{\Lambda}[\ln p(y\_{1:s}, \eta\_s, v\_{s\cdot}, v\_{s-1}, \sigma\_s | M\_s = i)] - E\_{\Lambda} \left[ \ln q \left( \eta\_s^i \right) q \left( v\_s^i \right) q \left( v\_{s-1}^i \right) q \left( \sigma\_s^i \right) \right] \tag{58}$$
  $\left. - \ln p \left( y\_{1:s-1} | M\_{\tilde{s}} = i \right) \right| \tag{59}$ 

where

$$\begin{aligned} p(y\_{1:s}, \eta\_s, \upsilon\_s, \upsilon\_{s-1}, \upsilon\_s | M\_s = i) &= p(y\_s | \eta\_s, \upsilon\_s, \upsilon\_{s-1}, \upsilon\_s, M\_s = i) p(\upsilon\_s | M\_s = i) p(y\_{1:s-1} | M\_s = i) \\ &\cdot p(\eta\_s | y\_{1:s-1}, M\_s = i) p(\upsilon\_s | M\_s = i) p(\upsilon\_{s-1} | M\_s = i) \end{aligned} \tag{59}$$

and Equation (60) is obtained:

$$\begin{aligned} L(\Lambda) &= E\_{\Omega} \left[ \ln p(y\_s | \eta\_{s\cdot} v\_{s\cdot} v\_{s-1}, \sigma\_{s\cdot} M\_{\sf s} = j) \right] + E\_{\Omega} \left[ \ln p(\eta\_{s\cdot} | y\_{1:s-1}, M\_{\sf s} = j) \right] \\ &+ E\_{\Omega} \left[ \ln p(v\_s | M\_{\sf s} = j) \right] + E\_{\Omega} \left[ \ln p(v\_{s-1} | \, M\_{\sf s} = j) \right] + E\_{\Omega} \left[ \ln p(\sigma\_{s} | M\_{\sf s} = j) \right] \\ &- E\_{\Omega} \left[ \ln q \left( \eta\_{s}^{j} \right) \right] - E\_{\Omega} \left[ \ln q \left( \boldsymbol{v}\_{s}^{j} \right) \right] - E\_{\Omega} \left[ \ln q \left( \boldsymbol{v}\_{s-1}^{j} \right) \right] - E\_{\Omega} \left[ \ln q \left( \boldsymbol{v}\_{s}^{j} \right) \right] \end{aligned} \tag{60}$$

based on Equation (60), the likelihood function in Equation (55) can be calculated. **Step 4:** Combination

The estimated system states and covariances of sub-filters are combined in this step and can be defined as follows:

$$\pounds\_{\mathfrak{s}|\mathfrak{s}} = \sum\_{i=1}^{k} \mu\_{\mathfrak{s}}^{i} \pounds\_{\mathfrak{s}|\mathfrak{s}}^{i} \tag{61}$$

$$P\_{\mathbf{s}|\mathbf{s}} = \sum\_{i=1}^{k} \mu\_{\mathbf{s}}^{i} \left( P\_{\mathbf{s}|\mathbf{s}}^{i} + \left( \mathfrak{A}\_{\mathbf{s}|\mathbf{s}}^{i} - \mathfrak{A}\_{\mathbf{s}|\mathbf{s}} \right) \left( \mathfrak{A}\_{\mathbf{s}|\mathbf{s}}^{i} - \mathfrak{A}\_{\mathbf{s}|\mathbf{s}} \right)^{T} \right) \tag{62}$$

Based on the abovementioned derivations, the design of the proposed method consists of the time update Equations (16)–(21) and the recursive measurement update Equations (43)–(51). Executable pseudocode of the presented algorithm can be found in Algorithm 1.


#### **5. Simulation Results of Maneuvering Target Tracking**

To verify the estimation performances of the designed method, this paper utilizes moving target tracking simulation experiments. The target moves by following the coordinate

turning (CT) and the constant velocity (CV) model alternately, where the dynamics of the CV model can be described as:

$$\mathbf{x}\_{\mathfrak{s}} = \begin{bmatrix} 1 & D & 0 & 0 \\ & 0 & 1 & 0 & 0 \\ & 0 & 0 & 1 & D \\ & 0 & 0 & 0 & 1 \end{bmatrix} \mathbf{x}\_{\mathfrak{s}-1} + \mathcal{g}\_{\mathfrak{s}-1} \tag{63}$$

where *D* represents sampling interval, the system state in CV model **x***<sup>s</sup>* = *px*, *vx*, *py*, *vy* comprises the positions *px*, *py* and the velocities *vx*, *vy* . The nominal process noise covariance matrix is provided as follows:

$$Q = \begin{bmatrix} & b\_1 \mathbf{U}\_{2 \times 2} & \mathbf{0}\_{2 \times 2} \\ & \mathbf{0}\_{2 \times 2} & b\_1 \mathbf{U}\_{2 \times 2} \end{bmatrix} \tag{64}$$

where the matrix

$$\mathbf{U}\_{2\times 2} = \begin{bmatrix} & \frac{D^3}{3} & \frac{D^2}{2} \\ & \frac{D^2}{2} & D \end{bmatrix} \tag{65}$$

and the CT model can be described by:

$$\mathbf{x}\_{s} = \begin{bmatrix} 1 & \frac{\sin(\varpi D)}{\mathcal{Q}} & 0 & \frac{\cos(\varpi D) - 1}{\mathcal{Q}} & 0\\ 0 & \cos(\varpi D) & 0 & -\sin(\varpi D) & 0\\ 0 & \frac{1 - \cos(\varpi D)}{\mathcal{Q}} & 1 & \frac{\sin(\varpi D)}{\mathcal{Q}} & 0\\ 0 & \sin(\varpi D) & 0 & \cos(\varpi D) & 0\\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \mathbf{x}\_{s-1} + \mathcal{g}\_{s-1}, \ \mathcal{Q} \neq 0 \tag{66}$$

where the system state in CT model **x***<sup>s</sup>* = *px*, *vx*, *py*, *vy*, comprises positions *px*, *py* , the velocities *vx*, *vy* , and the turning rate . The nominal process noise covariance matrix is in Equation (67)

$$Q = \begin{bmatrix} b\_1 \mathbf{U}\_{2 \times 2} & \mathbf{0}\_{2 \times 2} & \mathbf{0} \\ \mathbf{0}\_{2 \times 2} & b\_1 \mathbf{U}\_{2 \times 2} & \mathbf{0} \\ \mathbf{0}\_{1 \times 2} & \mathbf{0}\_{1 \times 2} & b\_2 D \end{bmatrix} \tag{67}$$

where *b*<sup>1</sup> and *b*<sup>2</sup> represent the power spectral densities.

One sensor locates in *px*0, *py*<sup>0</sup> receives the noise ranges and bearing measurements by the following formula [25]:

$$h(\mathbf{x}\_s) = \begin{bmatrix} \sqrt{\left(y\_s - p\_{y0}\right)^2 + \left(\mathbf{x}\_s - p\_{x0}\right)^2} \\ \arctan\left(\left(y\_s - p\_{y0}\right)/\left(\mathbf{x}\_s - p\_{x0}\right)\right) \end{bmatrix} \tag{68}$$

According to [20,34], the HTMNs are generated by

$$\varepsilon\_{\varepsilon} \sim \begin{cases} \begin{array}{l} N(0, R) \\ N(0, 100R) \end{array} & w.p. 0.9 \\ N(0, 100R) & w.p. 0.1 \end{cases} \tag{69}$$

where *w*.*p*. denotes "with probabilities of", the nominal covariance matrices of the measurement noises are assumed to be *R* = *diag* (10 m) 2 ,(0.1◦) 2 . The meaning of Equation (69) is that the measurement noise mostly obeys Gaussian distributions in a probability of 0.9, and follows Gaussian distributions with severely large covariances in a probability of 0.1.

The simulation time is from 0*s* to 1000 s. In detail, the moving target switches dynamic models alternately between the CT model and the CV model, i.e., the target executes the CV model with coordinate turn −5◦/s in 0–200 s, 401–600 s and 801–1000 s, while moves by the CT model in 201–400 s and 601–800 s. The sensor is positioned at (0 m, 0 m). Initialize the state vector *x*<sup>0</sup> = (10,000 m, 20 m/s, 10,000 m, 20 m/s, −5◦)/s and the covariance matrix *P*<sup>0</sup> = diag(100 m2, 50 m2/s2, 100 m2, 50 m2/s2, 100 m rad2/s2).

Four IMM approaches, containing the IMM cubature KF (IMM-CKF) [8], the robust STD-based IMM unscented KF (IMM-RSTUKF) [25], the Gaussian approximate filter with OSRMD [27] in the IMM framework (IMM-DGAF), and the designed IMM filter are compared simultaneously. All these algorithms are implemented with MATLAB 2020b in Windows 10, and the simulation experiments are run on a computer with Intel Core i5- 10400F CPU at 2.9 GHz and 16GB 2666MHz memory. The simulation parameters are summarized in Table 2. 1000 Monte Carlos are carried out for each filter. Two evaluation indicators are utilized to validate the performances of all algorithms, one is root mean square errors (RMSEs), the other one is averaged RMSEs (ARMSEs). According to [35,36], the RMSE and the ARMSE can be defined as follows:

$$RMSE\_{\text{pos}}(s) = \sqrt{\frac{1}{K\_{\text{mc}}} \sum\_{j=1}^{K\_{\text{mc}}} \left( \left( \hat{\mathbf{x}}\_s^{(j)} - \mathbf{x}\_s^{(j)} \right)^2 + \left( \hat{y}\_s^{(j)} - y\_s^{(j)} \right)^2 \right)} \tag{70}$$

$$ARMSE\_{\text{pos}} = \frac{1}{t\_s} \sum\_{s=1}^{t\_s} RMSE\_{\text{pos}}(s) \tag{71}$$

where *Kmc* denotes the total Monte Carlo run times, *xj <sup>s</sup>*, *y j s* and *x*ˆ *j <sup>s</sup>*, *y*ˆ *j s* represent the real and estimated positions in *j*-th Monte Carlo cycle, and *ts* refers to the total simulation time. Additionally, the RMSEs and the ARMSEs of the velocities and turning rates can be defined similarly to the abovementioned positions.

**Table 2.** Filtering parameters.


The surface maneuvering target tracking simulation can be divided into four parts to evaluate the performances of different algorithms. The first part considers the designed algorithm and existing ones regarding their estimation accuracy. Figures 1–3 displays the RMSEs of position, velocity, and turning rate of all algorithms. The IMM-CKF performs poorly due to the coexistence of the HTMNs with random measurement delay, significantly affecting filtering accuracy. It can be seen in Figures 1–3 that the proposed method has smaller RMSEs than existing filters. When the measurement outliers occur, the RMSE fluctuation of the proposed method is relatively small. This is because the proposed method is robust to the measurement outliers and adaptive to random measurement delay. As illustrated in Figures 1–3, the designed algorithm outperforms existing algorithms in estimation accuracy. Figures 4 and 5 provide the estimated and actual model probabilities for the two models. The estimated model probabilities from the designed filter are closer to the actual than the other filters. Figures 4 and 5 illustrate that the proposed method can better match the system model. Table 3 shows that the designed filter has smaller ARMSEs in each period in the scenarios of HTMNs coexisting with OSRMD. Compared with IMM-RSTUKF, which utilizes VB approach, the accuracy of ARMSE in position, velocity and

turning rate of the designed filter is improved by 41.46%, 47.96% and 6.92%, respectively. Moreover, the introduction of the VB inference increases computational costs.

**Figure 1.** RMSEs of position for different filters.

**Figure 2.** RMSEs of velocity for different filters.

**Figure 3.** RMSEs of turning rate for different filters.

**Figure 4.** The CV model probability for different filters.

**Figure 5.** The CT model probability for different filters.

**Table 3.** ARMSEs and SSRTs from all algorithms in different time periods. Filter 1, 2, 3, and 4 respectively represent the IMM cubature KF [8], the robust STD-based IMM unscented KF [25], the GAF with OSRMD in IMM filtering framework [27], and our designed robust IMM filter with HTMNs and OSRMD.


The second experiment part validates the impact of different OSRMD probabilities on estimation accuracy. The range of the measurement delay probability *ϕ<sup>s</sup>* is set to be from 0.1 to 0.9. In Figures 6–8, the ARMSEs of the positions, velocities, and turning rate from the designed filtering method and existing algorithms are provided. The horizontal axises of Figures 6–8 denote different measurement delay probabilities, and the longitudinal axises of Figures 6–8 represent the value of ARMSEs. From Figures 6–8, it can be seen that the IMM-CKF shows the worst estimation performance among all filters. The reason is that IMM-CKF deals with neither measurement outliers nor random measurement delay, leading to unsatisfactory estimation performance. As the probability of delay increases, the ARMSEs of IMM-RSTUKF has a relatively large change, the reason is that the IMM-RSTUKF lacks of adaptive estimation ability to random measurement delay. For the IMM-DGAF, although it has good adaptability to random measurement delay, it still has large estimation error due to inaccurate modeling of HTMNs. Since the proposed method can address the HTMNs and OSRMD simultaneously, this algorithm shows better adaptability to different measurement delays and obtain higher estimation accuracy in the presence of measurement outliers than existing algorithms.

**Figure 6.** ARMSEs of position with different measurement delay probabilities.

**Figure 7.** ARMSEs of velocity with different measurement delay probabilities.

**Figure 8.** ARMSEs of turning rate with different measurement delay probabilities.

The third part of the experiment aimed to verify the performances of all filters under different measurement outlier probabilities. The fixed measurement outlier probability value in Equation (69) is transformed into a range from 0.05 to 0.15. Figures 9–11 exhibits the ARMSEs from all algorithms under different measurement outlier probabilities. It can be seen from Figures 9–11 that the IMM-CKF still performs poorly, the conclusions about the IMM-CKF in the second experiment are further verified. In addition, with the increasing of the measurement outlier probability, the ARMSE value of the IMM-DGAF shows great fluctuation, this is because the IMM-DGAF assumes the measurement noises to obey Gaussian distributions, however, the Gaussian distribution is sensitive to outliers. As a result, the IMM-DGAF shows serious performance degradation to the increase of the measurement outlier probability. The reason for the IMM-RSTUKF performs worse than the proposed method is that the IMM-RSTUKF assumes all the measurement vectors can arrive in time, which is not satisfied when the OSRMD exists. The results in Figures 9–11 indicate that the presented algorithm can achieve better adaptive estimation performances than existing algorithms under different measurement outlier probabilities.

**Figure 9.** ARMSEs of position with different measurement outlier probabilities.

**Figure 10.** ARMSEs of velocity with different measurement outlier probabilities.

**Figure 11.** ARMSEs of turning rate with different measurement outlier probabilities.

In the fourth part of the experiment, the effects of the variational iterations on filtering accuracy are validated individually. In Figures 12–14, the ARMSEs from all filters under different variational iterations *N* = 1, 2, ..., 15 are shown. The simulation results in Figures 12–14 indicate that the presented algorithm converges when iteration number *N* = 2, and performs better than existing filters in estimation accuracy when *N* ≥ 2. The proposed algorithm converges faster than the IMM-RSTUKF which is also based on the VB framework. We can see from Figures 12–14 that the estimation accuracy of the IMM-CKF and IMM-DGAF is relatively stable since these two algorithms do not depend on the variational iterations. As the variational iteration *N* increases, the estimation performance of the designed filter improves. However, computational efficiency must be taken into account. Considering the balance between computing efficiency and estimation performances, the recommended variational iteration number is 4 to 8.

**Figure 12.** ARMSEs of position under different iteration numbers.

**Figure 13.** ARMSEs of velocity under different iteration numbers.

**Figure 14.** ARMSEs of turning rate under different iteration numbers.

### **6. Conclusions**

In the applications of surface maneuvering target tracking, due to the adverse effects of unreliable sensors and communication channel latency, the measurement outliers and random measurement delay often occur simultaneously. The state estimation issue of surface maneuver target faces severe challenges. Aiming at solving the problem of state estimation in JMSs with HTMNs and OSRMD, this paper designed a novel robust IMM filtering method and applied to surface maneuvering target tracking. Firstly, the HTMNs are assumed to follow STDs, and we introduced an RBV to characterize the random measurement delay. Secondly, this paper established a new HGSSM to utilize VB method. Finally, by using VB inference, the system state, RBV, and model probability are estimated simultaneously. The surface maneuvering target tracking simulation results indicated that the designed filtering method has better estimation accuracy than existing algorithms. In addition, compared with existing filters, the proposed filtering method showed better robustness to different outlier probabilities, and achieved better adaptive estimation performances under various of random measurement delay probabilities. Although the computational costs of the designed method increased slightly, the accuracy of ARMSE in position, velocity and turning rate is improved by 41.46%, 47.96% and 6.92% than the state-of-the-art, respectively. In our future work, based on the theoretical content and results of the proposed filtering method, more complex measurement environments will be considered, including the cases of multi-step random measurement delay and random measurement loss with unknown probability. The effect of heavy-tailed process noises on filtering accuracy and computational efficiency are also focuses of our future research.

**Author Contributions:** Conceptualization, C.C.; methodology, C.C.; software, C.C. and L.G.; data curation, W.Z.; writing—original draft preparation, C.C.; writing—review and editing, C.C. and L.G.; visualization, C.C.; supervision, W.Z.; project administration, W.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grants 61573113.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data was created or analyzed in this research work. Data sharing is not applicable to this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Proposition A1.** *Let φ* = *η<sup>i</sup> <sup>s</sup> in Equation (39), and* ln *q*(*L*+1) *ηi s is derived by:*

$$\begin{split} \ln q^{(L+1)} \left( \eta\_{s}^{i} \right) &= -0.5 \left( \eta\_{s}^{i} - \eta\_{s|s-1}^{i} \right)^{T} \left( P\_{s|s-1}^{\eta \eta, i} \right)^{-1} \left( \eta\_{s}^{i} - \eta\_{s|s-1}^{i} \right) \\ &- 0.5 E^{(L)} \left[ \sigma\_{s}^{i} \right] E^{(L)} \left[ \boldsymbol{v}\_{s-1}^{i} \right] \left( \boldsymbol{y}\_{s}^{i} - h\_{s-1}^{i} \left( \boldsymbol{x}\_{s-1}^{i} \right) \right)^{T} R\_{s-1}^{-1} \left( \boldsymbol{y}\_{s}^{i} - h\_{s-1}^{i} \left( \boldsymbol{x}\_{s-1}^{i} \right) \right) \\ &- 0.5 E^{(L)} \left[ 1 - \sigma\_{s}^{i} \right] E^{(L)} \left[ \boldsymbol{v}\_{s}^{i} \right] \left( \boldsymbol{y}\_{s}^{i} - h\_{s}^{i} \left( \boldsymbol{x}\_{s}^{i} \right) \right)^{T} R\_{s}^{-1} \left( \boldsymbol{y}\_{s}^{i} - h\_{s}^{i} \left( \boldsymbol{x}\_{s}^{i} \right) \right) + c\_{\eta} \tag{A1} \end{split} \tag{A1}$$

*where <sup>q</sup>*(*L*+1)(·) *is the approximate posterior PDF of <sup>p</sup>*(·) *when variational iteration is* (*<sup>L</sup>* <sup>+</sup> <sup>1</sup>) *. The likelihood PDF is modified as:*

$$p^{(L+1)}\left(\vec{y}\_s^i \middle| \eta\_s^i\right) = N\left(\vec{y}\_s^i; \bar{h}\_s^i \middle| \eta\_s^i\right), \mathcal{R}\_s^{(L+1)i}\right) \tag{A2}$$

*where y*¯*<sup>i</sup> s,* ¯ *hi s ηi s and <sup>R</sup>*¯(*L*+1)*<sup>i</sup> <sup>s</sup> represent the extended measurement vector, mean vector, and updated measurement noise covariances, they are defined as follows:*

$$
\Psi\_s^i = \begin{bmatrix} y\_s^i \\ y\_s^i \end{bmatrix} \tag{A3}
$$

$$
\bar{h}\_s^i \left( \eta\_s^i \right) = \left[ \begin{array}{c} \bar{h}\_s^i \left( \mathbf{x}\_s^i \right) \\ \bar{h}\_{s-1}^i \left( \mathbf{x}\_{s-1}^i \right) \end{array} \right] \tag{A4}
$$

$$
\bar{R}\_s^{(L+1)i} = \begin{bmatrix}
\frac{R\_s^i}{E^{(L)}\left[1 - \sigma\_s^i\right]E^{(L)}\left[\sigma\_s^i\right]} & \mathbf{0}\_{m \times m} \\
& \mathbf{0}\_{m \times m} & \frac{R\_{s-1}^i}{E^{(L)}\left[\sigma\_s^i\right]E^{(L)}\left[\sigma\_{s-1}^i\right]}
\end{bmatrix} \tag{A5}
$$

*Then, the q*(*L*+1) *ηi s can be updated by*

$$q^{(L+1)}\left(\eta\_s^i\right) = N\left(\eta\_s^i; \mathfrak{h}\_{s|s}^{(L+1)i}, P\_{s|s}^{(L+1)\eta\eta, i}\right) \tag{A6}$$

*where η*ˆ (*L*+1)*i <sup>s</sup>*|*<sup>s</sup> and P*(*L*+1)*ηη*,*<sup>i</sup> <sup>s</sup>*|*<sup>s</sup> are calculated through the following equations:*

$$K\_s^{(L+1)i} = P\_{\eta \bar{y}, s|s-1}^i \left( P\_{\mathcal{gg}, s|s-1}^{(L+1)i} \right)^{-1} \tag{A7}$$

$$\boldsymbol{\eta}\_{s|s}^{(L+1)i} = \boldsymbol{\eta}\_{s|s-1}^{i} + \boldsymbol{K}\_{s}^{(L+1)i} \left( \boldsymbol{\mathcal{Y}}\_{s}^{i} - \boldsymbol{\mathcal{Y}}\_{s|s-1}^{i} \right) \tag{A8}$$

$$P\_{s|s}^{(L+1)\eta\eta,i} = P\_{s|s-1}^{\eta\eta,i} - K\_s^{(L+1)i} P\_{gg|s|s-1}^{(L+1)i} \left(K\_{s-1}^{(L+1)i}\right)^T \tag{A9}$$

*where*

$$\mathbf{\hat{y}}\_{s|s-1}^{i} = \int \bar{h}\_{s}^{i} \left( \eta\_{s}^{i} \right) \mathbf{N} \left( \eta\_{s}^{i}; \boldsymbol{\eta}\_{s|s-1}^{i} \boldsymbol{P}\_{s|s-1}^{\eta \eta, i} \right) d\eta\_{s}^{i} \tag{A10}$$

$$\boldsymbol{P}\_{gg,s|s-1}^{\dot{\boldsymbol{r}}} = \int \ddot{\boldsymbol{h}}\_{s}^{i} \left(\boldsymbol{\eta}\_{s}^{i}\right) \left(\ddot{\boldsymbol{h}}\_{s}^{i} \left(\boldsymbol{\eta}\_{s}^{i}\right)\right)^{T} \boldsymbol{N} \left(\boldsymbol{\eta}\_{s}^{i}; \boldsymbol{\eta}\_{s|s-1}^{i}, \boldsymbol{P}\_{s|s-1}^{\eta \eta, i}\right) d\boldsymbol{\eta}\_{s}^{i} - \boldsymbol{\mathfrak{Y}}\_{s|s-1}^{i} \left(\boldsymbol{\mathfrak{Y}}\_{s|s-1}^{i}\right)^{T} + \boldsymbol{\mathcal{R}}\_{s}^{(L+1)i} \tag{A11}$$

$$P\_{\eta|\mathcal{S},s|s-1}^{i} = \int \eta\_s^i \left(\bar{h}\_s^i \left(\eta\_s^i\right)\right)^T \mathcal{N} \left(\eta\_s^i; \mathfrak{h}\_{s|s-1}^i, P\_{s|s-1}^{\eta\eta,i}\right) d\eta\_s^i - \eta\_{s|s-1}^i \left(\mathfrak{h}\_{s|s-1}^i\right)^T \tag{A12}$$

**Proposition A2.** *Let φ* = *v<sup>i</sup> <sup>s</sup> in Equation (39), and* ln *q*(*L*+1) *vi s is updated by Equation (A13):*

$$\begin{split} \ln q^{(L+1)} \left( v\_s^i \right) &= 0.5 \left( m E^{(L)} \left[ 1 - \sigma\_s^i \right] + a\_s^i - 2 \right) \ln v\_s^i \\ &- 0.5 \left\{ E^{(L)} \left[ 1 - \sigma\_s^i \right] \left( y\_s^i - h\_s^i \left( \mathbf{x}\_s^i \right) \right)^T R\_s^{-1} \left( y\_s^i - h\_s^i \left( \mathbf{x}\_s^i \right) \right) + a\_s^i \right\} v\_s^i + c\_v \quad \text{(A13)} \end{split} \tag{A13}$$

*and q*(*L*+1) *vi s can be updated by*

$$q^{(L+1)}\left(v\_s^i\right) = G\left(v\_s^i; a\_s^{(L+1)i}, \beta\_s^{(L+1)i}\right) \tag{A14}$$

*where <sup>α</sup>*(*L*+1)*<sup>i</sup> <sup>s</sup> and <sup>β</sup>*(*L*+1)*<sup>i</sup> <sup>s</sup> are calculated as:*

$$\begin{cases} \begin{aligned} a\_s^{(L+1)i} &= 0.5 \Big( m E^{(L)} \left[ 1 - \sigma\_s^i \right] + a\_s^i \Big) \\ \mathcal{B}\_s^{(L+1)i} &= 0.5 \Big\{ E^{(L)} \left[ 1 - \sigma\_s^i \right] tr \left( A\_s^{(L+1)} Y\_1 \left( B\_s^{(L+1)} \right)^{-1} \right) + a\_s^i \end{aligned} \end{cases} \tag{A15}$$

*where necessary matrices are defined as follows:*

$$\begin{cases} \begin{aligned} A\_s^{(L+1)i} &= \int \left( y\_s^i - h\_s^i(\mathbf{x}\_s^i) \right) \left( y\_s^i - h\_s^i(\mathbf{x}\_s^i) \right)^T N \left( \eta\_s^i; \mathfrak{H}\_{s|s}^{(L+1)i}, P\_{s|s}^{(L+1)\eta\eta; i} \right) d\eta\_s^i \\\ B\_s^{(L+1)i} &= \begin{bmatrix} \frac{R\_s^i}{E(L+1)\left[1-\sigma\_s^i\right]} & \mathbf{0}\_{m \times m} \\\ \mathbf{0}\_{m \times m} & \frac{R\_{s-1}^i}{E(L+1)\left[\sigma\_s^i\right]} \end{bmatrix} \\\ Y\_1 &= \begin{bmatrix} I\_{m \times m} & \mathbf{0}\_{m \times m} \\\ \mathbf{0}\_{m \times m} & \mathbf{0}\_{m \times m} \end{bmatrix}, Y\_2 = \begin{bmatrix} \mathbf{0}\_{m \times m} & \mathbf{0}\_{m \times m} \\\ \mathbf{0}\_{m \times m} & I\_{m \times m} \end{bmatrix} \end{aligned} \tag{A16}$$

*Im*×*<sup>m</sup> refers the identity matrix, and the necessary expectations are provided by:*

$$E^{(L+1)}\left[v\_s^i\right] = \frac{a\_s^{(L+1)i}}{\beta\_s^{(L+1)i}}\tag{A17}$$

$$E^{(L+1)}\left[\ln \upsilon\_s^i\right] = \psi\left(\alpha\_s^{(L+1)i}\right) - \ln\left(\beta\_s^{(L+1)i}\right) \tag{A18}$$

**Proposition A3.** *Let φ* = *v<sup>i</sup> <sup>s</sup>*−<sup>1</sup> *in Equation (39), and* ln *<sup>q</sup>*(*L*+1) *vi s*−1 *can be derived by the following equations:*

$$\begin{split} \ln q^{(L+1)} \left( \boldsymbol{v}\_{s-1}^{i} \right) &= 0.5 \left( m E^{(L)} \left[ \boldsymbol{o}\_{s}^{i} \right] + \boldsymbol{a}\_{s-1}^{i} - 2 \right) \ln \boldsymbol{v}\_{s-1}^{i} + \boldsymbol{c}\_{v} \\ &- 0.5 \left\{ E^{(L)} \left[ \boldsymbol{o}\_{s}^{i} \right] \left( \boldsymbol{y}\_{s}^{i} - \boldsymbol{h}\_{s-1}^{i} \left( \boldsymbol{x}\_{s-1}^{i} \right) \right)^{T} \boldsymbol{R}\_{s-1}^{-1} \left( \boldsymbol{y}\_{s}^{i} - \boldsymbol{h}\_{s-1}^{i} \left( \boldsymbol{x}\_{s-1}^{i} \right) \right) + \boldsymbol{a}\_{s-1}^{i} \right\} \boldsymbol{v}\_{s-1}^{i} \end{split} \tag{A19}$$

*and q*(*L*+1) *vi s*−1 *can be updated by*

$$q^{(L+1)}\left(v\_{s-1}^i\right) = G\left(v\_{s-1}^i; a\_{s-1}^{(L+1)i}, \beta\_{s-1}^{(L+1)i}\right) \tag{A20}$$

*where <sup>α</sup>*(*L*+1)*<sup>i</sup> <sup>s</sup>*−<sup>1</sup> *and <sup>β</sup>*(*L*+1)*<sup>i</sup> <sup>s</sup>*−<sup>1</sup> *are calculated as:*

$$\begin{cases} \boldsymbol{a}\_{s-1}^{(L+1)i} = 0.5 \Big( \boldsymbol{m} \boldsymbol{E}^{(L)} \left[ \sigma\_s^i \right] + \boldsymbol{a}\_{s-1}^i \Big) \\\ \boldsymbol{\rho}\_{s-1}^{(L+1)i} = 0.5 \Big\{ \boldsymbol{E}^{(L)} \left[ \sigma\_s^i \right] \boldsymbol{r} \boldsymbol{r} \Big( \boldsymbol{A}\_s^{(L+1)} \boldsymbol{Y}\_2 \left( \boldsymbol{B}\_s^{(L+1)} \right)^{-1} \Big) + \boldsymbol{a}\_{s-1}^i \Big\} \end{cases} \tag{A21}$$

*and the matrices uesd in Equation (A21) can be calculated by*

$$E^{(L+1)}\left[v\_{s-1}^i\right] = \frac{\alpha\_{s-1}^{(L+1)i}}{\beta\_{s-1}^{(L+1)i}}\tag{A22}$$

$$E^{(L+1)}\left[\ln v\_{s-1}^i\right] = \psi\left(\alpha\_{s-1}^{(L+1)i}\right) - \ln\left(\beta\_{s-1}^{(L+1)i}\right) \tag{A23}$$

**Proposition A4.** *Let φ* = *σ<sup>i</sup> <sup>s</sup> in Equation (39), and* ln *q*(*L*+1) *σi s is updated as follows:*

$$\begin{aligned} \ln q^{(L+1)}\left(\sigma\_{s}^{i}\right) &= -0.5\sigma\_{s}^{i}E^{(L+1)}\left[v\_{s-1}^{i}\right]\left(y\_{s}^{i}-h\_{s-1}^{j}\left(\mathbf{x}\_{s-1}^{i}\right)\right)^{T}\left(\mathbf{R}\_{s-1}^{i}\right)^{-1}\left(y\_{s}^{i}-h\_{s-1}^{j}\left(\mathbf{x}\_{s-1}^{i}\right)\right) \\ &+ 0.5mv\_{s}^{i}E^{(L+1)}\left[\ln v\_{s-1}^{i}\right]-0.5\left(1-\sigma\_{s}^{i}\right)E^{(L+1)}\left[v\_{s}^{i}\right]\left(y\_{s}^{i}-h\_{s}^{j}\left(\mathbf{x}\_{s}^{i}\right)\right)^{T}\left(\mathbf{R}\_{s}^{i}\right)^{-1}\left(y\_{s}^{i}-h\_{s}^{j}\left(\mathbf{x}\_{s}^{i}\right)\right) \\ &+ 0.5m\left(1-\sigma\_{s}^{i}\right)E^{(L+1)}\left[\ln v\_{s}^{i}\right]+\left(1-\sigma\_{s}^{i}\right)\ln\left(1-q\_{s}^{j}\right)+\sigma\_{s}^{i}\ln q\_{s}^{j}+\text{cst}\_{\sigma} \end{aligned} \tag{A24}$$

*the RBV σ<sup>s</sup> selecting values from* {0, 1} *is updated as follows:*

$$Pr^{(L+1)}\left(\sigma\_s^i = 1\right) = \Omega^{(L+1)i} \exp\left\{ S\_1^{(L+1)i} \right\} \tag{A25}$$

$$Pr^{(L+1)}\left(\sigma\_s^i = 0\right) = \Omega^{(L+1)i} \exp\left\{ S\_1^{(L+1)i} \right\} \tag{A26}$$

*where the* <sup>Ω</sup>(*L*+1)*<sup>i</sup> denotes the normalized constant, and <sup>S</sup>*(*L*+1)*<sup>i</sup>* <sup>1</sup> *and <sup>S</sup>*(*L*+1)*<sup>i</sup>* <sup>2</sup> *are calculated as follows:*

$$S\_1^{(L+1)i} = E^{(L+1)} \left[ \ln \sigma\_s^i \right] - 0.5tr \left( A\_s^{(L+1)i} Y\_1 \left( \mathbb{C}\_s^{(L+1)i} \right)^{-1} \right) + \ln \left( 1 - q\_s^i \right) \tag{A27}$$

$$S\_1^{(L+1)i} = E^{(L+1)} \left[ \ln v\_{s-1}^i \right] - 0.5tr \left( A\_s^{(L+1)i} Y\_2 \left( C\_s^{(L+1)i} \right)^{-1} \right) + \ln \left( q\_s^i \right) \tag{A28}$$

*and the matrix C*(*L*+1)*<sup>i</sup> <sup>s</sup> is formulated as:*

$$\mathbf{C}\_{s}^{(L+1)i} = \begin{bmatrix} \frac{R\_{s}}{E^{(L+1)}\left[v\_{s}^{i}\right]} & \mathbf{0}\_{m \times m} \\ \mathbf{0}\_{m \times m} & \frac{R\_{s-1}}{E^{(L+1)}\left[v\_{s-1}^{i}\right]} \end{bmatrix} \tag{A29}$$

*the necessary expectations E*(*L*+1) ' *σi s* ( *and E*(*L*+1) ' <sup>1</sup> − *<sup>σ</sup><sup>i</sup> s* ( *can be obtained by:*

$$E^{(L+1)}\left[\sigma\_s^i\right] = \frac{\Pr^{(L+1)}\left(\sigma\_s^i = 1\right)}{\Pr^{(L+1)}\left(\sigma\_s^i = 1\right) + \Pr^{(L+1)}\left(\sigma\_s^i = 0\right)}\tag{A30}$$

$$E^{(L+1)}\left[1-\sigma\_s^i\right] = 1 - E^{(L+1)}\left[\sigma\_s^i\right] \tag{A31}$$

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
