Robust Tracking of Multi-Station Underwater Maneuvering Targets Under Skewed Heavy-Tailed Noise

1. Introduction

Advances in unmanned underwater platforms and acoustic sensing technologies have increased the need for accurate underwater maneuvering target tracking. By recursively estimating target states from noisy acoustic measurements, tracking algorithms support underwater vehicle localization, navigation, and marine monitoring. Maintaining estimation accuracy in the presence of target maneuvers and complex measurement noise remains a key challenge.

Maneuvering target tracking requires accurate state estimation from noisy measurements while accommodating changes in target motion. Within the Bayesian estimation framework, various filtering methods have been developed, including the Kalman filter (KF) [1,2], extended Kalman filter (EKF) [3,4], unscented Kalman filter (UKF) [5,6,7], and particle filter (PF) [8]. The UKF uses the unscented transform to approximate the mean and covariance of random variables undergoing nonlinear transformations and has been applied to underwater maneuvering target tracking. The PF can address more general nonlinear, non-Gaussian estimation problems, although its estimation performance and computational cost depend on the number of particles and the resampling strategy[9]. These methods provide a foundation for target state estimation; however, reliance on a single fixed motion model may lead to model mismatch when the target switches between motion modes, thereby degrading tracking performance.

To address model mismatch caused by changes in target motion, the interacting multiple model (IMM) algorithm improves adaptability through model interaction, parallel filtering, and model-probability-weighted fusion. Li et al. [10] developed a Bayesian filtering framework incorporating delayed measurements to address random measurement delays and model switching in underwater target tracking, with adaptive adjustment of model transition probabilities. Yang et al. [11] combined adaptive and robust filtering within the IMM framework to address simultaneous process and measurement modeling errors. Guo et al. [12] introduced a maneuver-adaptive model set to improve adaptability to target maneuvers compared with a fixed model set. More recently, Zheng et al. [13] combined variational Bayesian estimation with a multiple-model approach for underwater extended target tracking with inaccurate noise covariances. Although these studies enhance model adaptability in different ways, the choice of process noise covariance still affects the trade-off between responsiveness to maneuvers and steady-state estimation error. Adapting the process noise covariance to changes in target motion therefore offers a means of improving maneuvering target tracking performance.

In addition to model mismatch caused by target maneuvers, non-Gaussian measurement noise also limits underwater tracking performance. To address estimation errors caused by measurement outliers, Stojanovic and Nedic [14] investigated robust nonlinear Kalman filtering under non-Gaussian noise. Whereas robust update schemes reduce the influence of measurement outliers, heavy-tailed distribution models explicitly characterize the probability of large errors. Roth et al. [15] modeled both process and measurement noise using Student’s t distributions and developed a robust filtering method for linear systems. Huang et al. [16] extended Student’s t filtering and smoothing to nonlinear systems, using the unscented transform to approximate the associated weighted integrals. Xue et al. [17] developed a robust filtering framework based on elliptically contoured distributions and employed a multiple-model mechanism to adaptively select the noise distribution, improving adaptability to different noise conditions. In underwater acoustic ranging, however, range estimates are affected by errors in both the effective sound speed and propagation-delay estimates. Sound speed variations, multipath propagation, and erroneous delay estimates may introduce systematic biases and large measurement errors[18]. When measurement errors exhibit both skewness and heavy tails, symmetric heavy-tailed distributions or symmetric noise mixture models may not adequately characterize their statistical properties, and the resulting noise model mismatch may further degrade tracking accuracy. Huang et al. [19] proposed a robust Kalman filtering framework based on Gaussian scale mixture distributions, employing variational Bayesian inference to estimate the state and noise-related variables under heavy-tailed and skewed noise. For skew-t measurement noise, Nurminen et al. [20] improved the treatment of dependencies among variables in the posterior approximation and incorporated expectation propagation to improve covariance estimation. To address nonstationary noise in underwater target tracking, Zhou et al. [21] modeled process noise and measurement noise using heavy-tailed and skewed distributions, respectively, and combined variational Bayesian estimation with distributed fusion for state estimation across multiple platforms. These studies provide a methodological foundation for state estimation under skewed heavy-tailed noise, while the update strategies for noise-related parameters remain important determinants of model adaptability and estimation performance.

Sliding-window estimation has recently been incorporated into adaptive robust filtering to make greater use of recent observations in noise parameter estimation. Wang et al. [22] modeled heavy-tailed measurement noise using a Pearson type VII distribution and jointly updated the state and noise-related parameters through state smoothing and variational iterations within the window. Based on a Bernoulli-Gaussian measurement model, Jiang et al. [23] developed centralized and distributed variational Bayesian moving horizon estimation methods to jointly estimate the state trajectory and unknown noise parameters within the window. For vessel tracking under heavy-tailed noise and random transmission delays, Liu et al. [24] combined a Student’s t noise model, measurement time alignment, and variational Bayesian moving horizon estimation. These studies improve the use of information from multiple time steps through state and parameter estimation within a moving window. For recursive tracking, further investigation is warranted into how recent local estimates can improve the stability of parameter updates without repeatedly solving for the joint posterior of the state trajectory and parameters within the window.

To address the above issues, this article proposes an Adaptive Interacting Multiple Model Sliding-Window Variational Bayesian Unscented Filter (AIMM-SWVBUF). The main contributions are summarized as follows:

1.

A recursive estimation method combining per-time-step variational updates with sliding-window parameter fusion is developed under a generalized hyperbolic skewed Student’s t (GHSST) measurement noise model. The target state and noise-related variables are estimated using the unscented transform, while estimates of the location parameter and mixing variable are fused within the window to reduce the influence of fluctuations in local estimates on subsequent filtering updates.

2.

A sliding-window maximum a posteriori (MAP) estimation method is developed to address the experimentally observed tendency of shape parameter estimates to remain near their initial values with limited updates. By combining observations from multiple time steps within the window with the parameter prior, the method updates the shape parameter separately to improve its estimation.

3.

An adaptive process noise covariance adjustment strategy based on changes in model probabilities is designed. The strategy dynamically adjusts the process noise covariance in response to changes in target motion, balancing responsiveness during maneuvers with tracking accuracy during steady motion.

2. Materials and Methods

2.1. Problem Formulation

Establishment of a Nonlinear State-Space Model:

$$\left\{\begin{array}{c}{\mathit{x}}_k={\mathit{F}}_{k-1}{\mathit{x}}_{k-1}+{\mathit{G}}_{k-1}{\mathit{w}}_{k-1}\hfill \\ {\mathit{z}}_k=\mathit{h}\left({\mathit{x}}_k\right)+{\mathit{v}}_k\hfill \end{array}\right.$$

(1)

Where, k denotes the discrete time index, and ${\mathit{x}}_k=[x_{p,k},{\dot{x}}_{p,k},$ $y_{p,k},{\dot{y}}_{p,k},{\omega }_k{]}^{\mathrm{T}}\in {\mathbb{R}}^n$ is the state vector, $(x_{p,k},y_{p,k})$ and $({\dot{x}}_{p,k},{\dot{y}}_{p,k})$ represent the target position and velocity components in a two-dimensional Cartesian coordinate system, respectively, and ${\mathit{\omega }}_k$ denotes the target turn rate. The measurement vector is denoted by ${\mathit{z}}_k\in {\mathbb{R}}^m$. ${\mathit{F}}_{k-1}$ is the state transition matrix corresponding to the current motion model, and ${\mathit{G}}_{k-1}$ is the process noise weighting matrix. The function $\mathit{h}(·)$ denotes the nonlinear measurement function. The vectors ${\mathit{w}}_{k-1}\in {\mathbb{R}}^3$ and ${\mathit{v}}_k\in {\mathbb{R}}^m$ represent the process noise and the measurement noise, respectively, and they are assumed to be mutually independent.

In complex underwater environments, the combined effects of acoustic propagation mechanisms and highly variable environmental conditions cause the measurement noise to exhibit both skewness and heavy-tailed characteristics. First, the skewness of the noise mainly originates from systematic biases and asymmetries in acoustic propagation paths. On the one hand, the sound speed profile, which is influenced by temperature, salinity, and depth, shows significant spatiotemporal inhomogeneity [25], leading to systematic errors in the estimation of the effective sound speed. As a result, range measurements statistically deviate from the zero-mean symmetric assumption. On the other hand, multipath propagation introduces asymmetries in both energy distribution and arrival times between the direct path and reflected paths [18], which further yields asymmetric components in the measurement errors. Second, the heavy-tailed nature of the noise is mainly attributed to complex environmental interference and outliers. During active sonar operation, reverberation severely interferes with time-delay estimation, inducing measurement outliers. Meanwhile, impulsive noise, ship-radiated noise, and biological noise that are ubiquitous in underwater environments [26] can easily introduce outliers into the measurement data, causing the noise distribution to exhibit heavy-tailed behavior. In summary, multiple physical factors jointly break the symmetry and stationarity of the noise distribution, giving rise to a dual non-Gaussian nature characterized by both nonzero-mean skewness and heavy tails. Under such conditions, the conventional zero-mean Gaussian assumption or purely symmetric heavy-tailed models become seriously mismatched to the actual noise statistics. Therefore, in this article, the process noise is modeled as Gaussian to represent random perturbations in the target motion, while the contaminated measurement noise is modeled by a Generalized Hyperbolic Skewed Student’s t (GHSST) distribution to accommodate the complex underwater acoustic observation environment. This distribution can be expressed in the form of a Gaussian-inverse-gamma mixture as follows [19].

$$p\left({\mathit{w}}_k\right)=\mathcal{N}({\mathit{w}}_k;0,{\mathit{Q}}_k)$$

(2)

$$\begin{array}{cc}\hfill p\left({\mathit{v}}_k\right)& =\mathrm{GHSST}({\mathit{v}}_k;{\mathit{u}}_k^{\mathrm{noise}},{\mathit{R}}_k,\mathit{\beta },{\xi }_k)\hfill \\ \hfill & =\int \mathcal{N}({\mathit{v}}_k;{\mathit{u}}_k^{\mathrm{noise}}+{\lambda }_k\mathit{\beta },{\lambda }_k{\mathit{R}}_k)\mathrm{IG}({\lambda }_k;{\xi }_k,{\xi }_k)d{\lambda }_k\hfill \end{array}$$

(3)

Where ${\mathit{Q}}_k$ denotes the process noise covariance matrix, ${\mathit{u}}_k^{\mathrm{noise}}$ is the location parameter of the measurement noise, accounting for an additive bias, ${\mathit{R}}_k$ is the measurement noise scale matrix, ${\lambda }_k>0$ is the latent mixing variable, $\mathit{\beta }$ is the shape parameter governing the distribution’s skewness, ${\xi }_k$ is the degrees-of-freedom parameter, and $\mathrm{IG}(·)$ denotes the inverse-gamma distribution.

The measurement equation can be expressed as:

$$\begin{array}{cc}\hfill {\mathit{z}}_k& =\left[\begin{array}{c}\sqrt{{(x_{p,k}-x_1)}^2+{(y_{p,k}-y_1)}^2}\\ \begin{array}{c}\sqrt{{(x_{p,k}-x_2)}^2+{(y_{p,k}-y_2)}^2}\hfill \\ \sqrt{{(x_{p,k}-x_3)}^2+{(y_{p,k}-y_3)}^2}\hfill \end{array}\end{array}\right]+{\mathit{v}}_k\hfill \end{array}$$

(4)

Where $(x_i,y_i),(i=1,2,3)$ denote the locations of the three observation points, respectively. The observation geometry is shown in Figure 1.

2.2. Overall Architecture

The AIMM-SWVBUF consists of three key components. First, the proposed Sliding-Window Variational Bayesian Unscented Filter (SWVBUF) is employed as the core estimation algorithm to perform joint inference of the state and noise parameters under skewed heavy-tailed measurement noise. Second, an IMM mechanism is introduced to achieve maneuvering target tracking through parallel model-conditioned filtering and probabilistic model updating. Finally, an adaptive process noise covariance adjustment is designed, which dynamically adjusts the noise parameters based on the changes in model probabilities, thereby achieving a balance between rapid maneuver response and high steady-state accuracy. Its schematic structure is shown in Figure 2.

2.3. Sliding-Window Variational Bayesian Unscented Filter

The following derivation focuses on the single-model formulation. The extension to the IMM-based multi-model framework will be detailed in Section 2.4. For notational simplicity, the model index is omitted in the subsequent derivations.

2.3.1. Single-Step Variational Inference

According to (1)-(3), the one-step predictive probability density function $p\left({\mathit{x}}_k\right|{\mathit{z}}_{1:k-1})$ and the likelihood function $p\left({\mathit{z}}_k\right|{\mathit{x}}_k)$ can be expressed as follows:

$$p\left({\mathit{x}}_k\right|{\mathit{z}}_{1:k-1})=\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})$$

(5)

$$\begin{array}{c}\hfill p({\mathit{z}}_k\mid {\mathit{x}}_k)=\int \mathcal{N}({\mathit{z}}_k;\mathit{h}\left({\mathit{x}}_k\right)+{\mathit{u}}_k^{\mathrm{noise}}+{\lambda }_k\mathit{\beta },{\lambda }_k{\mathit{R}}_k)\times \mathrm{IG}({\lambda }_k;{\xi }_k,{\xi }_k)\mathrm{d}{\lambda }_k\end{array}$$

(6)

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{k|k-1}={\mathit{F}}_{k-1}{\widehat{\mathit{x}}}_{k-1|k-1}\end{array}$$

(7)

$$\begin{array}{c}\hfill {\mathit{P}}_{k|k-1}={\mathit{F}}_{k-1}{\mathit{P}}_{k-1|k-1}{\mathit{F}}_{k-1}^{\mathrm{T}}+{\mathit{G}}_{\mathrm{k}-1}{\mathit{Q}}_{k-1}{\mathit{G}}_{k-1}^{\mathrm{T}}\end{array}$$

(8)

Where ${\mathit{x}}_k$ denotes the system state at time k, ${\mathit{z}}_{1:k-1}$ represents the measurement set from time 1 to $k-1$, ${\widehat{\mathit{x}}}_{k|k-1}$ is the one-step predicted state mean, and ${\mathit{P}}_{k|k-1}$ is the corresponding one-step predicted state covariance matrix.

To compensate for the systematic measurement errors caused by biases in the effective sound speed in underwater environments, ${\mathit{u}}_k^{\mathrm{noise}}$ is introduced as an unknown parameter to characterize the persistent nonzero-mean bias. Meanwhile, to further describe the heavy-tailed and skewed characteristics of the dual noise sources in complex underwater acoustic environments, ${\lambda }_k$ is introduced, and the joint estimation of ${\mathit{u}}_k^{\mathrm{noise}}$ and ${\lambda }_k$ is carried out within the variational inference framework. Specifically, ${\mathit{u}}_k^{\mathrm{noise}}$ is assumed to follow a Gaussian prior, while ${\lambda }_k$ follows an inverse-gamma prior distribution. Combined with (6), the modified conditional likelihood function and parameter prior distributions are defined as follows:

$$p\left({\mathit{z}}_k\right|{\mathit{x}}_k,{\mathit{u}}_k^{\mathrm{noise}},{\lambda }_k)=\mathcal{N}({\mathit{z}}_k;\mathit{h}\left({\mathit{x}}_k\right)+{\mathit{u}}_k^{\mathrm{noise}}+{\lambda }_k\mathit{\beta },{\lambda }_k{\mathit{R}}_k)$$

(9)

$$p\left({\mathit{u}}_k^{\mathrm{noise}}\right)=\mathcal{N}({\mathit{u}}_k^{\mathrm{noise}};{\mathit{\phi }}_k,{\mathsf{\Phi }}_k)p\left({\lambda }_k\right)=\mathrm{IG}({\lambda }_k;{\xi }_k,{\xi }_k)$$

(10)

Let the set of latent variables at time k be denoted by ${\mathsf{\Theta }}_k\stackrel{\Delta }{=}\left\{{\mathit{x}}_k,{\mathit{u}}_k^{\mathrm{noise}},{\lambda }_k\right\}$. Since direct analytical evaluation of the posterior distribution $p\left({\mathsf{\Theta }}_k\right|{\mathit{z}}_{1:k})$ is intractable, a variational Bayesian approach is adopted to seek a tractable approximation $q\left({\mathsf{\Theta }}_k\right)$ to the true posterior. This problem is equivalent to minimizing the Kullback-Leibler divergence (KLD) between the two distributions [27].

$$\begin{array}{c}\hfill q^{*}\left({\mathsf{\Theta }}_k\right)=\underset{q\left({\mathsf{\Theta }}_k\right)}{\arg \min }\mathrm{KLD}(q\left({\mathsf{\Theta }}_k\right)\left|\right|p\left({\mathsf{\Theta }}_k\right|{\mathit{z}}_{1:k}))\end{array}$$

(11)

$$\begin{array}{c}\hfill \mathrm{KLD}(q\left({\mathsf{\Theta }}_k\right)\left|\right|p\left({\mathsf{\Theta }}_k\right|{\mathit{z}}_{1:k}))=\int q\left({\mathsf{\Theta }}_k\right)\ln \left({\displaystyle \frac{q\left({\mathsf{\Theta }}_k\right)}{p({\mathsf{\Theta }}_k,{\mathit{z}}_{1:k})}}\right)\mathrm{d}{\mathsf{\Theta }}_k+\ln p\left({\mathit{z}}_{1:k}\right)\end{array}$$

(12)

To obtain a tractable analytical solution, based on the mean-field theory[28], the variational approximation $q\left({\mathsf{\Theta }}_k\right)$ is assumed to factorize into a product of independent components:

$$q\left({\mathsf{\Theta }}_k\right)=q\left({\mathit{x}}_k\right)q\left({\mathit{u}}_k^{\mathrm{noise}}\right)q\left({\lambda }_k\right)$$

(13)

Since $\ln p\left({\mathit{z}}_{1:k}\right)$ is a constant, minimizing the KLD is equivalent to maximizing the Evidence Lower Bound (ELBO)[28]:

$$\mathcal{L}\left(q\left({\mathsf{\Theta }}_k\right)\right)={\mathrm{E}}_{{\mathsf{\Theta }}_k}\left[\ln \left(p(\mathsf{\Theta },{\mathit{z}}_{1:k})\right)-\ln \left(q\left({\mathsf{\Theta }}_k\right)\right)\right]$$

(14)

$$q^{*}\left({\mathsf{\Theta }}_k\right)=\underset{q\left({\mathsf{\Theta }}_k\right)}{\arg \max }\mathcal{L}\left(q\left({\mathsf{\Theta }}_k\right)\right)$$

(15)

The factors are iteratively updated using the Coordinate Ascent Variational Inference algorithm. For any element ${\mathit{\theta }}_k$ in ${\mathsf{\Theta }}_k$ , its optimal variational posterior update is given by:

$$\log q^{*}\left({\mathit{\theta }}_k\right)={\mathrm{E}}_{{\mathsf{\Theta }}_k^{(-\mathit{\theta })}}\left[\log p({\mathsf{\Theta }}_k,{\mathit{z}}_{1:k})\right]+{\mathrm{C}}_{\mathsf{\Theta }}$$

(16)

Where ${\mathrm{E}}_{{\mathsf{\Theta }}_k^{(-\mathit{\theta })}}[·]$ denotes the expectation with respect to all variables in ${\mathsf{\Theta }}_k$ except ${\mathit{\theta }}_k$, and ${\mathrm{C}}_{\mathsf{\Theta }}$ is a normalization constant.

Next, taking the $(i+1)$-th iteration as an example, the update procedures for the posterior estimate of the target state, the posterior estimate of the measurement-noise mean, and the posterior estimate of the mixing parameter are provided.

Let ${\mathit{\theta }}_k={\mathit{x}}_k$. As shown in Appendix B, its posterior distribution follows a Gaussian distribution. Therefore, the VB update can be written as:

$$q^{(i+1)}\left({\mathit{x}}_k\right)=\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k}^{(i+1)},{\mathit{P}}_{k|k}^{(i+1)})$$

(17)

Where ${\widehat{\mathit{x}}}_{k|k}^{(i+1)}$ and ${\mathit{P}}_{k|k}^{(i+1)}$ can be computed as follow:

$$\begin{array}{c}\hfill {\mathit{b}}_k=\int \mathit{h}\left({\mathit{x}}_k\right)\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_{\mathrm{k}}\end{array}$$

(18)

$$\begin{array}{c}\hfill {\tilde{\mathit{q}}}_k^{\left(i\right)}={\mathrm{E}}^{\left(i\right)}\left[{\mathit{u}}_k^{\mathrm{noise}}\right]+{\mathrm{E}}^{\left(i\right)}\left[{\lambda }_k\right]\mathit{\beta }\end{array}$$

(19)

$$\begin{array}{c}\hfill {\widehat{\mathit{z}}}_{k|k-1}^{(i+1)}={\mathit{b}}_k+{\tilde{\mathit{q}}}_k^{\left(i\right)}\end{array}$$

(20)

$$\begin{array}{c}\hfill {\tilde{\mathit{R}}}_k^{\left(i\right)}={\mathit{R}}_k/{\mathrm{E}}^{\left(i\right)}\left[{\lambda }_k^{-1}\right]\end{array}$$

(21)

$$\begin{array}{c}\hfill {\mathit{P}}_{k|k-1}^{zz(i+1)}=\int \mathit{h}\left({\mathit{x}}_k\right){\mathit{h}}^{\mathrm{T}}\left({\mathit{x}}_k\right)\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_{\mathrm{k}}-{\mathit{b}}_{\mathrm{k}}{\left({\mathit{b}}_{\mathrm{k}}\right)}^{\mathrm{T}}+{\tilde{\mathit{R}}}_{\mathrm{k}}^{\left(\mathrm{i}\right)}\end{array}$$

(22)

$$\begin{array}{c}\hfill {\mathit{P}}_{k|k-1}^{xz}=\int {\mathit{x}}_k{\mathit{h}}^{\mathrm{T}}\left({\mathit{x}}_k\right)\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1})\mathrm{d}{\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-1}{\left({\mathit{b}}_k\right)}^{\mathrm{T}}\end{array}$$

(23)

$$\begin{array}{c}\hfill {\mathit{K}}^{(i+1)}={\mathit{P}}_{k|k-1}^{xz}{\left({\mathit{P}}_{k|k-1}^{zz(i+1)}\right)}^{-1}\end{array}$$

(24)

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{k|k}^{(i+1)}={\widehat{\mathit{x}}}_{k|k-1}+{\mathit{K}}^{(i+1)}({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1}^{(i+1)})\end{array}$$

(25)

$$\begin{array}{c}\hfill {\mathit{P}}_{k|k}^{(i+1)}={\mathit{P}}_{k|k-1}-{\mathit{K}}^{(i+1)}{\mathit{P}}_{k|k-1}^{zz(i+1)}{\left({\mathit{K}}^{(i+1)}\right)}^{\mathrm{T}}\end{array}$$

(26)

Let ${\mathit{\theta }}_k={\mathit{u}}_k^{\mathrm{noise}}$. As shown in Appendix C, its posterior distribution follows a Gaussian distribution. Therefore, the VB update can be written as:

$$\begin{array}{c}\hfill q^{(i+1)}\left({\mathit{u}}_k^{\mathrm{noise}}\right)=\mathcal{N}({\mathit{u}}_k^{\mathrm{noise}};{\mathit{\phi }}_k^{(i+1)},{\mathsf{\Phi }}_k^{(i+1)})\end{array}$$

(27)

Where ${\mathit{\phi }}_k^{(i+1)}$ and ${\mathsf{\Phi }}_k^{(i+1)}$ can be computed as follows:

$$\begin{array}{c}\hfill {\mathit{K}}_u^{(i+1)}={\mathsf{\Phi }}_k^{\left(i\right)}{({\mathsf{\Phi }}_k^{\left(i\right)}+{\tilde{\mathit{R}}}_k^{\left(i\right)})}^{-1}\end{array}$$

(28)

$$\begin{array}{c}\hfill {\mathit{c}}_k={\mathit{z}}_k-\int \mathit{h}\left({\mathit{x}}_k\right)\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k}^{(i+1)},{\mathit{P}}_{k|k}^{(i+1)})\mathrm{d}{\mathit{x}}_k-{\mathrm{E}}^{\left(i\right)}\left[{\lambda }_k\right]\mathit{\beta }\end{array}$$

(29)

$$\begin{array}{c}\hfill {\mathit{\phi }}_k^{(i+1)}={\mathit{\phi }}_k^{\left(i\right)}+{\mathit{K}}_u^{(i+1)}({\mathit{c}}_k-{\mathit{\phi }}_k^{\left(i\right)})\end{array}$$

(30)

$${\mathsf{\Phi }}_k^{(i+1)}={\mathsf{\Phi }}_k^{\left(i\right)}-{\mathit{K}}_u^{(i+1)}{\mathsf{\Phi }}_k^{\left(i\right)}$$

(31)

Let ${\mathit{\theta }}_k={\lambda }_k$. As shown in Appendix D, its posterior distribution follows the generalized inverse gaussian (GIG) distribution. Therefore, the VB update can be written as:

$$\begin{array}{c}\hfill q^{(i+1)}\left({\lambda }_k\right)=\mathrm{GIG}\left({\lambda }_k;{\xi }_k^1,{\xi }_k^2,{\xi }_k^3\right)\end{array}$$

(32)

$$\begin{array}{c}\hfill {\xi }_k^1=-{\xi }_k-{\displaystyle \frac{m}{2}}{\xi }_k^2=\mathrm{tr}\left(\mathit{\beta }{\mathit{\beta }}^{\mathrm{T}}{\mathit{R}}_k^{-1}\right){\xi }_k^3=2{\xi }_k+\mathrm{tr}\left({\mathit{d}}_k{\mathit{R}}_k^{-1}\right)\end{array}$$

(33)

Where ${\mathit{d}}_k$ can be computed as:

$$\begin{array}{c}\hfill {\mathit{d}}_k=\int {\mathit{d}}_k^1{\left({\mathit{d}}_k^1\right)}^{\mathrm{T}}\mathcal{N}({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k}^{(i+1)},{\mathit{P}}_{k|k}^{(i+1)})\mathrm{d}{\mathit{x}}_k+{\mathsf{\Phi }}_k^{(i+1)}+{\mathit{d}}_k^2{\left({\mathit{d}}_k^2\right)}^{\mathrm{T}}\end{array}$$

(34)

$$\begin{array}{c}\hfill {\mathit{d}}_k^1={\mathrm{E}}^{\left(\mathrm{i}+1\right)}\left[\mathit{h}\left({\mathit{x}}_k\right)\right]-\mathit{h}\left({\mathit{x}}_k\right){\mathit{d}}_k^2={\mathit{z}}_k-{\mathrm{E}}^{\left(\mathrm{i}+1\right)}\left[\mathit{h}\left({\mathit{x}}_k\right)\right]-{\mathit{\phi }}_k^{(i+1)}\end{array}$$

(35)

Since computing the moments of the GIG distribution requires evaluating modified Bessel functions of the second kind, its numerical implementation is relatively complicated and introduces additional computational burden. To ensure the real-time implementation of the algorithm, a MAP approach is adopted in this paper, where the complete posterior distribution is approximated by its mode, thereby yielding the update of ${\lambda }_k$ :

$$\begin{array}{c}\hfill {\lambda }_k^{(i+1)}=-(2{\xi }_k+2+m)/2{\rho }_2+\sqrt{{(2{\xi }_k+2+m)}^2+4{\rho }_2(2{\xi }_k+{\rho }_1)}/2{\rho }_2\end{array}$$

(36)

$$\begin{array}{c}\hfill {\mathrm{E}}^{(i+1)}\left[{\mathit{u}}_k^{\mathrm{noise}}\right]={\mathit{\phi }}_k^{(i+1)}{\mathrm{E}}^{(i+1)}\left[{\lambda }_k\right]={\lambda }_k^{(i+1)}\end{array}$$

(37)

To reduce computational complexity and avoid unnecessary repeated calculations, a convergence criterion is established for the variational Bayesian fixed-point iteration process. The iterative update is terminated when the relative change between the posterior state estimates obtained in two successive iterations satisfies the following condition:

$$\begin{array}{c}\hfill {\displaystyle \frac{∥{\widehat{\mathit{x}}}_{k|k}^{\left(i\right)}-{\widehat{\mathit{x}}}_{k|k}^{(i-1)}∥}{∥{\widehat{\mathit{x}}}_{k|k}^{(i-1)}∥}}\le \tau \end{array}$$

(38)

Where $\tau >0$ denotes the convergence threshold, and $∥·∥$ represents the Euclidean norm.

This criterion uses the relative variation of the state estimate as the basis for convergence. When the update magnitude falls below the predefined threshold, the algorithm is considered to have reached a stable state, and the iteration is terminated. Furthermore, to ensure numerical stability under extreme conditions, a maximum number of iterations $N_m$ is specified, and the iteration process is forcibly stopped once this upper bound is reached.

2.3.2. Sliding-Window Parameter Fusion

In recursive Bayesian estimation, information from past observations can be propagated to the current time step through the prior distribution. Therefore, a recursive update at a single time step does not imply that only the current observation is used. In practical variational inference, however, parameter estimation performance depends on the approximate posterior structure, the prior propagation scheme, and the quality of the current measurement. Under skewed heavy-tailed measurement noise, measurement outliers may cause substantial fluctuations in local parameter estimates, thereby compromising the stability of target state estimation.

To mitigate these fluctuations, a sliding-window parameter fusion mechanism is introduced following the per-time-step variational update. For noise statistics that vary slowly over short time intervals, local estimates from recent time steps are fused to obtain parameter estimates at the current time step.

This subsection considers the fusion of estimates of the location parameter ${\mathit{u}}_k^{\mathrm{noise}}$ and the mixing variable ${\lambda }_k$. To distinguish these estimates from those obtained during the iterations, let ${\mathit{u}}_k^{\mathrm{noise}\left(\mathrm{loc}\right)}$ and ${\lambda }_k^{\mathrm{noise}\left(\mathrm{loc}\right)}$ denote , respectively, the posterior mean of the location parameter and the MAP estimate of the mixing variable after termination of the variational iterations, but before sliding-window fusion. Define

$${\mathit{\zeta }}_k^{\left(\mathrm{loc}\right)}=\left[\begin{array}{c}{\mathit{u}}_k^{\mathrm{noise}\left(\mathrm{loc}\right)}\hfill \\ {\lambda }_k^{\left(\mathrm{loc}\right)}\hfill \end{array}\right]\in {\mathbb{R}}^m\times (0,\infty )$$

(39)

Let L denote the prescribed window length. For $k\ge L$, the sliding window comprises the L most recent time steps:

$${\mathit{W}}_k=\left\{k-L+1,\dots ,k\right\}$$

(40)

The effective window length is defined as the number of local estimates included in the fusion:

$${\mathit{L}}_k^{\mathrm{eff}}=\left|{\mathit{W}}_k\right|$$

(41)

During initialization, the effective window length ${\mathit{L}}_k^{\mathrm{eff}}$ is used as the normalization factor in the averaging operation. Once the window is full, ${\mathit{L}}_k^{\mathrm{eff}}=L$.

For the location parameter, slow variation within the window is assumed:

$${\mathit{u}}_t^{\mathrm{noise}}\approx {\mathit{u}}_k^{\mathrm{noise}}\forall t\in {\mathit{W}}_k$$

(42)

Under this assumption, local estimates of the location parameter at different time steps can be approximately regarded as estimates of a common offset. Fusing these estimates can mitigate fluctuations in per-time-step estimates. The mixing variable ${\lambda }_t$ characterizes the random scale of the measurement at each time step and is allowed to vary across observations. Its estimates are also smoothed over the window as an additional algorithmic approximation. The resulting estimate is used in subsequent filtering computations and is not interpreted as a single underlying mixing variable shared across all time steps within the window.

To fuse the local estimates of the location parameter and the mixing variable within a unified formulation, the following objective function is constructed using a weighted least-squares criterion:

$${\mathcal{J}}_k^{\mathrm{fus}}\left(\mathit{\zeta }\right)={\displaystyle \frac{1}{2}}\sum _{t\in {\mathit{W}}_k}{(\mathit{\zeta }-{\mathit{\zeta }}_t^{\left(\mathrm{loc}\right)})}^{\mathrm{T}}{\mathit{M}}_t(\mathit{\zeta }-{\mathit{\zeta }}_t^{\left(\mathrm{loc}\right)})$$

(43)

Where, ${\mathit{M}}_t\in {\mathbb{R}}^{(m+1)\times (m+1)}$ is a positive definite diagonal weighting matrix.

Taking the gradient of (43) with respect to $\mathit{\zeta }$ yields

$${\nabla }_{\mathit{\zeta }}{\mathcal{J}}_k^{\mathrm{fus}}=\sum _{t\in {\mathit{W}}_k}{\mathit{M}}_t(\mathit{\zeta }-{\mathit{\zeta }}_t^{\left(\mathrm{loc}\right)})$$

(44)

Since the Hessian $\sum _{t\in {\mathit{W}}_k}{\mathit{M}}_t$ is positive definite, the objective function is strictly convex. Setting the gradient in (44) to zero gives the unique stationary point:

$${\mathit{\zeta }}_k^{\left(\mathrm{SW}\right)}={\left(\sum _{t\in {\mathit{W}}_k}{\mathit{M}}_t\right)}^{-1}\sum _{t\in {\mathit{W}}_k}{\mathit{M}}_t{\mathit{\zeta }}_t^{\left(\mathrm{loc}\right)}$$

(45)

The positive definite diagonal structure of the weighting matrices ensures that each component of the fused estimate is a convex combination of the corresponding local estimates. Because ${\lambda }_t^{\left(\mathrm{loc}\right)}>0$, the fused mixing variable remains strictly positive. Therefore, (45) is feasible and constitutes the unique global minimizer of (43) over the specified domain.

For computational simplicity, a common weighting matrix is assigned to all time steps within the window:

$${\mathit{M}}_t=\mathit{M}\forall t\in {\mathit{W}}_k$$

(46)

Where, $\mathit{M}$ is a fixed positive definite diagonal matrix. Substituting (46) into (45) reduces the fusion rule to an arithmetic average:

$${\mathit{\zeta }}_k^{\left(\mathrm{SW}\right)}={\displaystyle \frac{1}{{\mathit{L}}_k^{\mathrm{eff}}}}\sum _{t\in {\mathit{W}}_k}{\mathit{\zeta }}_t^{\left(\mathrm{loc}\right)}$$

(47)

The resulting parameter updates are

$${\mathit{u}}_k^{\mathrm{noise}\left(\mathrm{SW}\right)}={\displaystyle \frac{1}{{\mathit{L}}_k^{\mathrm{eff}}}}\sum _{t\in {\mathit{W}}_k}{\mathit{u}}_t^{\mathrm{noise}\left(\mathrm{loc}\right)}{\lambda }_k^{\left(\mathrm{SW}\right)}={\displaystyle \frac{1}{{\mathit{L}}_k^{\mathrm{eff}}}}\sum _{t\in {\mathit{W}}_k}{\lambda }_t^{\left(\mathrm{loc}\right)}$$

(48)

Equation (48) performs sliding-window averaging of the local posterior means of the location parameter and the local MAP estimates of the mixing variable. For $L=1$, the updates reduce to the corresponding local estimates at the current time step.

2.3.3. Sliding-Window MAP Estimation

As shown in Figure 3, the $\mathit{\beta }$ estimates obtained by the joint variational inference method in [19] remain essentially around the initial value over the tracking sequence, exhibiting only slight fluctuations and showing no clear tendency to converge toward the true parameter value. This indicates that, under the joint variational inference framework, the posterior update of $\mathit{\beta }$ is subject to strong coupling constraints and suffers from insufficient identifiability, thereby limiting its correction capability. To address this issue, $\mathit{\beta }$ is no longer treated as a complete latent variable in the joint variational inference at the current time step; instead, a sliding-window MAP scheme is adopted for its separate estimation. By incorporating observations from multiple adjacent time steps within the window, the proposed method strengthens the statistical support for $\mathit{\beta }$ estimation and thereby improves its identifiability. The results show that the proposed method enables the estimate of $\mathit{\beta }$ to move away from the initial value and gradually approach the true value, indicating that the proposed estimation strategy effectively overcomes the stagnation problem in variational updating and achieves an effective update of $\mathit{\beta }$. Furthermore, a more accurate estimate of $\mathit{\beta }$ helps improve the modeling accuracy of skewed heavy-tailed noise and ultimately enhances state tracking performance.

On this basis, a sliding-window MAP estimation model for $\mathit{\beta }$ is further established. Specifically, within a sliding window of length L, denoted by ${\mathcal{W}}_k$, $\mathit{\beta }$ is assumed to remain constant and to follow a Gaussian prior distribution:

$$\begin{array}{c}\hfill p\left({\mathit{\beta }}_k^{\left(\mathrm{SW}\right)}\right)=\mathcal{N}({\mathit{\beta }}_k^{\left(\mathrm{SW}\right)};{\mathit{\beta }}_0,{\Sigma }_0)\end{array}$$

(49)

Where ${\mathit{\beta }}_0$ and ${\Sigma }_0$ denote the prior mean and prior covariance, respectively.

Then, the sliding-window MAP estimate of $\mathit{\beta }$ is defined as:

$$\begin{array}{c}\hfill {\mathit{\beta }}_k^{\left(\mathrm{SW}\right)}=\underset{\mathit{\beta }}{\arg \max }\mathcal{L}\left(\mathit{\beta }\right)\end{array}$$

(50)

$$\begin{array}{c}\hfill \mathcal{L}\left(\mathit{\beta }\right)={\mathrm{E}}_{q\left({\mathsf{\Theta }}_k\right)}\left[\sum _{t\in {\mathit{W}}_k}\log p\left({\mathit{z}}_t\right|{\mathsf{\Theta }}_t,\mathit{\beta })+\log p\left(\mathit{\beta }\right)\right]\end{array}$$

(51)

Where ${\mathsf{\Theta }}_k$ denotes the set of latent variables other than $\mathit{\beta }$, which have already been estimated by variational Bayesian inference.

According to (9), $\log p\left({\mathit{z}}_t\right|{\mathsf{\Theta }}_t,\mathit{\beta })$ can be written as:

$$\begin{array}{c}\hfill \log p\left({\mathit{z}}_t\right|{\mathsf{\Theta }}_t,\mathit{\beta })={\mathit{\beta }}^{\mathrm{T}}{\mathit{R}}_t^{-1}{\mathit{\sigma }}_t^1-{\displaystyle \frac{1}{2}}{\lambda }_t{\mathit{\beta }}^{\mathrm{T}}{\mathit{R}}_t^{-1}\mathit{\beta }+C_{\mathit{\beta }}\end{array}$$

(52)

Where ${\mathit{\sigma }}_t^1={\mathit{z}}_t-\mathit{h}\left({\mathit{x}}_t\right)-{\mathit{u}}_t^{\mathrm{noise}}$. (51) can be further expressed as:

$$\begin{array}{c}\hfill \mathcal{L}\left(\mathit{\beta }\right)=-{\displaystyle \frac{1}{2}}{\mathit{\beta }}^{\mathrm{T}}\left({\Sigma }_0^{-1}+\sum _{t\in {\mathcal{W}}_k}\mathrm{E}\left[{\lambda }_t\right]{\mathit{R}}_t^{-1}\right)\mathit{\beta }+{\mathit{\beta }}^{\mathrm{T}}\left({\Sigma }_0^{-1}{\mathit{\beta }}_0+\sum _{t\in {\mathcal{W}}_k}{\mathit{R}}_t^{-1}\mathrm{E}\left[{\mathit{\sigma }}_t^1\right]\right)+C_{\mathit{\beta }}\end{array}$$

(53)

Taking the derivative of (53) with respect to $\mathit{\beta }$ and setting it to zero yields:

$$\begin{array}{c}\hfill \left({\Sigma }_0^{-1}+\sum _{t\in {\mathcal{W}}_k}\mathrm{E}\left[{\lambda }_t\right]{\mathit{R}}_t^{-1}\right)\mathit{\beta }-\left({\Sigma }_0^{-1}{\mathit{\beta }}_0+\sum _{t\in {\mathcal{W}}_k}{\mathit{R}}_t^{-1}\mathrm{E}\left[{\mathit{\sigma }}_t^1\right]\right)=0\end{array}$$

(54)

Therefore, the closed-form solution is obtained as:

$$\begin{array}{c}\hfill {\mathit{\beta }}_k^{\left(\mathrm{SW}\right)}={\left({\Sigma }_0^{-1}+\sum _{t\in {\mathcal{W}}_k}\mathrm{E}\left[{\lambda }_t\right]{\mathit{R}}_t^{-1}\right)}^{-1}\times \left({\Sigma }_0^{-1}{\mathit{\beta }}_0+\sum _{t\in {\mathcal{W}}_k}{\mathit{R}}_t^{-1}\mathrm{E}\left[{\mathit{\sigma }}_t^1\right]\right)\end{array}$$

(55)

During initialization, single-time-step MAP estimation is used to update the shape parameter until the sliding window is full. The corresponding update is given by

$$\begin{array}{c}\hfill {\mathit{\beta }}_k={\left({\Sigma }_0^{-1}+\mathrm{E}\left[{\lambda }_k\right]{\mathit{R}}_k^{-1}\right)}^{-1}\left({\Sigma }_0^{-1}{\mathit{\beta }}_0+{\mathit{R}}_k^{-1}\mathrm{E}\left[{\mathit{\sigma }}_k^1\right]\right)\end{array}$$

(56)

The proposed sliding-window variational Bayesian unscented filter (SWVBUF) is composed of (7)-(8), (18) -(26), (28)-(38), (48), (55). To evaluate the required weighted integrals, the UT is employed for numerical approximation.

2.4. IMM-Based Multi-Model Extension

2.4.1. Model Interaction

The model interaction step mixes the posterior state estimates of all models at time $k-1$ using the Markov transition probability matrix, thereby providing the initial conditions for each model filter at time k. Assume that the system consists of N models, and the transition probability from model i to model j is denoted by ${\mathsf{\Pi }}_{ij}$.

First, the mixing probability ${\mathit{\mu }}_{k-1|k-1}^{ij}$ from model i to model j is computed as:

$$\begin{array}{c}\hfill {\mathit{\mu }}_{k|k-1}^j={\displaystyle \sum _{i=1}^N}{\mathsf{\Pi }}_{ij}{\mathit{\mu }}_{k-1}^i{\mathit{\mu }}_{k-1|k-1}^{ij}={\displaystyle \frac{{\mathsf{\Pi }}_{ij}{\mathit{\mu }}_{k-1}^i}{{\mathit{\mu }}_{k|k-1}^j}}\end{array}$$

(57)

Where ${\mathit{\mu }}_{k-1}^i$ denotes the probability of model i at time $k-1$, and ${\mathit{\mu }}_{k|k-1}^j$ denotes the predicted probability of model j at time k, which also serves as the normalization constant.

Based on the above mixing probabilities, the mixed initial state ${\widehat{\mathit{x}}}_{k-1|k-1}^{0j}$ and the mixed initial covariance ${\mathit{P}}_{k-1|k-1}^{0j}$ for model j at time k can be computed as:

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{k-1|k-1}^{0j}={\displaystyle \sum _{i=1}^N}{\widehat{\mathit{x}}}_{k-1|k-1}^i{\mathit{\mu }}_{k-1|k-1}^{ij}\end{array}$$

(58)

$$\begin{array}{c}\hfill {\mathit{P}}_{k-1|k-1}^{0j}={\displaystyle \sum _{i=1}^N}{\mathit{\mu }}_{k-1|k-1}^{ij}\left({\mathit{P}}_{k-1|k-1}^i+{\mathit{d}}_{k-1}^{ij}{\left({\mathit{d}}_{k-1}^{ij}\right)}^{\mathrm{T}}\right)\end{array}$$

(59)

Where ${\mathit{d}}_{k-1}^{ij}={\widehat{\mathit{x}}}_{k-1|k-1}^i-{\widehat{\mathit{x}}}_{k-1|k-1}^{0j}$. ${\widehat{\mathit{x}}}_{k-1|k-1}^i$ and ${\mathit{P}}_{k-1|k-1}^i$ denote the posterior state estimate and covariance estimate of model i at time $k-1$, respectively. The mixed initial state ${\widehat{\mathit{x}}}_{k-1|k-1}^{0j}$ and the mixed initial covariance ${\mathit{P}}_{k-1|k-1}^{0j}$ are then used as inputs to the filter matched to model j at time k.

2.4.2. Model Filtering

For each motion model, the SWVBUF described in Section 2.3 is executed independently using the corresponding mixed initial state and covariance obtained from the model interaction step. Specifically, conditioned on the current motion model, the state and noise parameters are jointly estimated through the sliding-window variational Bayesian inference procedure.

It is worth emphasizing that the internal structure of the SWVBUF remains unchanged across different models; the only difference lies in the model-specific state transition matrix and process noise configuration.

2.4.3. Model Probability Update

After completing the variational iterations, the likelihood function of each model is updated using the corrected innovation residual $e_k^j$.

$$\begin{array}{c}\hfill {\mathit{e}}_k^j={\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1}^j\end{array}$$

(60)

The likelihood value ${\mathsf{\Lambda }}_k^j$ and the model probability ${\mathit{\mu }}_k^j$ of model j are then computed as:

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_k^j=\mathcal{N}\left({\mathit{e}}_k^j;0,{\mathit{P}}_{k\mid k-1}^{\left(zz\right)j}\right){\mathit{\mu }}_k^j={\displaystyle \frac{\left({\mathsf{\Lambda }}_k^j{\mathit{\mu }}_{k|k-1}^j\right)}{\left({\displaystyle \sum _{i=1}^N}{\mathsf{\Lambda }}_k^i{\mathit{\mu }}_{k|k-1}^i\right)}}\end{array}$$

(61)

2.4.4. State Estimate Fusion

Finally, the global state estimate ${\widehat{\mathit{x}}}_{k|k}$ and the corresponding covariance matrix ${\mathit{P}}_{k|k}$ are obtained by a weighted sum of the filtering results of all models as follows:

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{k|k}={\displaystyle \sum _{j=1}^N}{\mathit{\mu }}_k^j{\widehat{\mathit{x}}}_{k|k}^j\end{array}$$

(62)

$$\begin{array}{c}\hfill {\mathit{P}}_{k|k}={\displaystyle \sum _{j=1}^N}{\mathit{\mu }}_k^j\left[{\mathit{P}}_{k|k}^j+({\widehat{\mathit{x}}}_{k|k}-{\widehat{\mathit{x}}}_{k|k}^j){({\widehat{\mathit{x}}}_{k|k}-{\widehat{\mathit{x}}}_{k|k}^j)}^{\mathrm{T}}\right]\end{array}$$

(63)

$$\begin{array}{c}\hfill {\mathit{u}}_k^{\mathrm{noise}}={\displaystyle \sum _{j=1}^N}{\mathit{\mu }}_k^j{\mathit{u}}_k^{\mathrm{noise},\mathrm{j}}{\lambda }_k={\displaystyle \sum _{j=1}^N}{\mathit{\mu }}_k^j{\lambda }_k^j{\mathit{\beta }}_k={\displaystyle \sum _{j=1}^N}{\mathit{\mu }}_k^j{\mathit{\beta }}_k^j\end{array}$$

(64)

2.5. Adaptive Process Noise Covariance Adjustment

In underwater maneuvering target tracking, the setting of the process noise covariance matrix ${\mathit{Q}}_{k-1}$ involves a trade-off between response speed and tracking accuracy. When ${\mathit{Q}}_{k-1}$ is chosen too large, excessive random perturbations are introduced during steady-state phases, degrading tracking accuracy. Conversely, when ${\mathit{Q}}_{k-1}$ is too small, the filter fails to promptly adapt to rapid state variations during target maneuvers, resulting in delayed response. To resolve this dilemma, an adaptive regulation strategy based on changes in model probabilities is proposed to achieve a dynamic balance between the two.

2.5.1. Construction of the Modification Function

Within the IMM framework, the model probability ${\mathit{\mu }}_k^j$ can reflect in real time the degree of matching between the target motion mode and model j. When the target is in steady motion, the probabilities of all models remain relatively stable; when a maneuver occurs, the probability of the corresponding maneuvering model exhibits a pronounced increasing trend. Therefore, the changes in model probabilities can serve as a sensitive indicator for maneuver detection.

The increment of the maneuvering model probability $\Delta {\mathit{\mu }}_{k-1}^j$ at time k is defined as:

$$\begin{array}{c}\hfill \Delta {\mathit{\mu }}_{k-1}^j={\mathit{\mu }}_{k-1}^j-{\mathit{\mu }}_{k-2}^j\end{array}$$

(65)

Based on this increment, an adaptive modification function for the process noise is constructed as $f(k-1)$:

$$\begin{array}{c}\hfill f(k-1)={\displaystyle \frac{1}{1-\Delta {\mathit{\mu }}_{k-1}^j}}\end{array}$$

(66)

Figure 4 illustrates an example of the variation of the modification function.

2.5.2. Process Noise Covariance Regulation

On this basis, the modification term of the process noise covariance matrix is defined as:

$$\begin{array}{c}\hfill {\tilde{\mathit{Q}}}_{k-1}=\psi (k-1){\mathit{Q}}_{k-1}\end{array}$$

(67)

where $\psi \left(k\right)$ is the adaptive scaling matrix, defined as:

$$\begin{array}{c}\hfill \psi (k-1)=\left\{\begin{array}{c}\mathrm{diag}(1,1,1)\left|f\left(k-1\right)-1\right|\le \delta \hfill \\ \mathrm{diag}(\mathrm{c},\mathrm{c},\mathrm{c})\left|f\left(k-1\right)-1\right|>\delta \hfill \end{array}\right.\end{array}$$

(68)

Where $\delta$ denotes the tolerance threshold and c is the process noise covariance inflation factor.

When the variation of model probabilities is small, i.e., the value of the modification function lies within $\left[-\delta ,\delta \right]$, the target is regarded as being in a steady-state phase. In this case, $\psi (k-1)=\mathrm{diag}(1,1,1)$, and the nominal process noise covariance ${\mathit{Q}}_{k-1}$ is retained to avoid introducing unnecessary disturbances and to ensure high estimation accuracy. When the model probabilities exhibit abrupt variations such that $\left|f\left(k-1\right)-1\right|>\delta$, the target is considered to be in a maneuvering phase. Then, the process noise covariance is enlarged by a factor c, so as to enhance the filter responsiveness to rapid state changes.

By integrating the derivations presented in Section 2.3, Section 2.4 and Section 2.5. the complete procedure of the AIMM-SWVBUF algorithm can be obtained. The corresponding implementation pseudocode is summarized in Algorithm 1.

Algorithm 1 Pseudocode of the Proposed AIMM-SWVBUF for One Time Step

Require: ${\widehat{\mathit{x}}}_{k-1|k-1}$, ${\mathit{P}}_{k-1|k-1}$, ${\tilde{\mathit{Q}}}_{k-1}$, ${\mathit{R}}_k$, ${\mathsf{\Pi }}_{ij}$, ${\mathit{\mu }}_{k-1}$, ${\xi }_k$, ${\lambda }_k$, ${\phi }_k$, N, $N_m$, ${\mathit{\beta }}_0$, L, c, $\delta$, $\tau$ Ensure:  ${\widehat{\mathit{x}}}_{k|k}$, ${\mathit{P}}_{k|k}$ 1: Initialization: ${\mathrm{E}}^{\left(0\right)}\left[{\lambda }_k\right]={\lambda }_k$, ${\mathrm{E}}^{\left(0\right)}\left[{\mathit{u}}_k^{\mathrm{noise}}\right]={\mathit{\phi }}_k$ 2: Compute ${\mathit{\mu }}_{k-1|k-1}^{ij}$, ${\widehat{\mathit{x}}}_{k-1|k-1}^{0j}$, ${\mathit{P}}_{k-1|k-1}^{0j}$ using (57)-(59). 3: Compute $\mathit{\psi }(k-1)$ using (66) and (68). 4: for  $j=1,\dots ,N$do 5:     Compute ${\widehat{\mathit{x}}}_{k|k-1}^{\left(0\right)}$ and ${\mathit{P}}_{k|k-1}^{\left(0\right)}$ using (7) and (8). 6:     for $i=1,\dots ,N_m$ do 7:         Compute $q^{\left(i\right)}\left({\mathit{x}}_k\right)$ using (17)-(26) 8:         Compute $q^{\left(i\right)}\left({\mathit{u}}_k^{\mathrm{noise}}\right)$ using (27)-(31) 9:         Compute $q^{\left(i\right)}\left({\lambda }_k\right)$ using (32)-(36) 10:         Compute ${\mathrm{E}}^{\left(i\right)}\left[{\mathit{u}}_k^{\mathrm{noise}}\right]$ and ${\mathrm{E}}^{\left(i\right)}\left[{\lambda }_k\right]$ using (37) 11:         if $∥{\widehat{\mathit{x}}}_{k|k}^{\left(i\right)}-{\widehat{\mathit{x}}}_{k|k}^{(i-1)}∥/∥{\widehat{\mathit{x}}}_{k|k}^{(i-1)}∥\le \tau$ then 12:            break; 13:         end if 14:     end for 15:     Compute ${\mathit{u}}_k^{\mathrm{noise}\left(\mathrm{SW}\right)}$, ${\lambda }_k^{\left(\mathrm{SW}\right)}$ and ${\mathit{\beta }}_k^{\left(\mathrm{SW}\right)}$ using (48) and (55) 16:     Update ${\mathit{q}}^{(i+1)}\left({\mathit{x}}_k\right)$ using (17)-(26) 17: end for 18: Compute ${\mathit{\mu }}_k^j$ using (61) 19: Compute ${\widehat{\mathit{x}}}_{k|k}$, ${\mathit{P}}_{k|k}$, ${\mathit{u}}_k^{\mathrm{noise}}$, ${\mathit{\beta }}_k$ and ${\lambda }_k$ using (62)-(64)

2.6. Computational Complexity Analysis

This subsection analyzes the per-time-step computational complexity of the proposed AIMM-SWVBUF algorithm. The main computational costs arise from filtering updates within each model, sliding-window parameter fusion, sliding-window MAP estimation of the shape parameter, and model interaction and estimate fusion within the IMM framework. Let n and m denote the state and measurement dimensions. Let N denote the number of models, L the sliding-window length, and $S=2n+1$ the number of sigma points used in the unscented transform. The computational complexity of each main step is summarized in Table 1.

For the single-model SWVBUF, let $N_d$ denote the number of variational iterations performed at the current time step. The corresponding computational complexity is

$$\begin{array}{cc}\hfill C_{\mathrm{SWVBUF}}=& O\left(N_d(n^3+n{m}^2+n^{2}m+m^3)+N_{d}S(n^2+nm+m^2)+L{m}^2\right)\hfill \end{array}$$

(69)

Where $N_d$ is the number of iterations required to satisfy the convergence criterion. For the tracking problem considered in this study, $n>m$. With $S=2n+1$, the above expression simplifies to

$$\begin{array}{c}\hfill C_{\mathrm{SWVBUF}}=O\left(N_d{n}^3+L{m}^2\right)\end{array}$$

(70)

Within the IMM framework, filtering updates are performed separately for N models. Model interaction involves computing the mixed state mean and covariance for each model, with a total computational complexity of $O\left(N^2{n}^2\right)$. Model probability updating and final state estimate fusion incur computational costs of $O\left(N^2\right)$ and $O\left(N{n}^2\right)$, respectively. Adaptive process noise adjustment primarily involves evaluating changes in model probabilities and scaling the process noise covariance; it therefore does not alter the asymptotic order of the dominant computations. Consequently, the per-time-step computational complexity of AIMM-SWVBUF is bounded by

$$C_{\mathrm{AIMM}-\mathrm{SWVBUF}}=O\left(N(N_d{n}^3+L{m}^2)+N^2{n}^2\right)$$

(71)

Equation (71) shows that the computational cost of the proposed algorithm depends primarily on the state dimension, the number of models, the number of variational iterations, and the window length. Figure 5 presents the number of variational iterations performed by AIMM-SWVBUF at each time step in the considered experiment, with a mean of approximately 6.4. A comparison of the measured execution times is provided in Section 3.4.

3. Results and Discussion

To evaluate the tracking performance of the proposed algorithm in complex underwater acoustic environments, a three-stage validation scheme is designed.

First, a parameter sensitivity analysis is conducted. Specifically, the sliding-window length and the process noise amplification factor are systematically swept to investigate their impacts on tracking accuracy. This enables the identification of an optimal parameter configuration that balances accuracy and computational efficiency, thereby establishing a baseline for the subsequent experiments.

Second, an ablation study is performed using the selected parameter configuration. A baseline algorithm and variants obtained by removing individual modules are compared with the complete algorithm to assess the contributions of sliding-window parameter fusion, adaptive process noise covariance adjustment, and sliding-window MAP estimation of the shape parameter to tracking performance.

Finally, comprehensive comparative experiments are conducted. Building upon the optimal parameter configuration and the validated effectiveness of each module, the proposed algorithm is compared with the Adaptive Interacting Multiple Model Unscented Kalman Filter (AIMM-UKF)[12], the Interacting Multiple Model Unscented Kalman Filter (IMM-UKF)[6], the Robust Student’s t Unscented Kalman Filter (STUKF)[16], the Student’s t -Generalized Hyperbolic Skewed Student’s t Kalman Filter (ST-GHSST-KF)[19], and Robust Distributed Variational Bayesian Moving Horizon Estimation with a Bernoulli-Gaussian Measurement Noise Model (R-DVB-MHE-BG)[23]. It is worth noting that STUKF, ST-GHSST-KF and R-DVB-MHE-BG were originally developed as single-model algorithms. To ensure fairness and validity in maneuvering target tracking scenarios, they are extended to the IMM framework in this study, thereby enabling a comprehensive evaluation of the proposed algorithm’s overall advantages in handling skewed heavy-tailed measurement noise and target maneuvers.

3.1. Simulation Setup

This simulation considers a two-dimensional underwater target tracking scenario. The target motion is modeled using the state-space formulation introduced in Section 2.1, and tracking is performed using a Constant Velocity (CV) model and a Constant Turn Rate (CT) model. The corresponding state transition matrices ${\mathit{F}}_{k-1}$ and noise weighting matrices ${\mathit{G}}_{k-1}$ follow those in [6]. The total sampling duration is 800 s, with a sampling interval of $T=4\mathrm{s}$. The target motion models over different time periods are listed in Table 2.

The process noise and measurement noise are specified as follows:

$$\begin{array}{c}\hfill {\mathit{Q}}_k=\mathrm{diag}\left[{\left(0.01\mathrm{m}\right)}^2,{\left(0.01\mathrm{m}/\mathrm{s}\right)}^2,{\left(0.01^{\mathrm{o}}/\mathrm{s}\right)}^2\right]\end{array}$$

(72)

$$\left\{\begin{array}{l}{\mathit{v}}_k\sim \mathcal{N}\left({\gamma }_k+\mathsf{\Omega }{\mathit{\eta }}_k,{\mathsf{\Gamma }}_k^{-1}\overline{\mathit{R}}\right)\\ {\mathit{\eta }}_k\sim {\mathcal{N}}_{+}\left(0,{\mathsf{\Gamma }}_k^{-1}\right)\\ {\left[{\mathsf{\Gamma }}_k\right]}_{ii}\sim \mathrm{G}\left({\displaystyle \frac{\alpha }{2}},{\displaystyle \frac{\alpha }{2}}\right)\end{array}\right.$$

(73)

Where $\overline{\mathit{R}}=100{\mathrm{I}}_3$ is the noise scale matrix, ${\gamma }_k={[5,5,5]}^{\mathrm{T}}$ denotes the location parameter, $\mathsf{\Omega }=5{\mathrm{I}}_3$ denotes the shape parameter , ${\mathsf{\Gamma }}_k$ is a $3\times 3$ diagonal matrix whose diagonal elements ${\left[{\mathsf{\Gamma }}_k\right]}_{ii}$ are independent and identically distributed, ${\mathit{\eta }}_k$ is an auxiliary random vector, and the degrees of freedom are set to $\alpha =5$.

In this simulation, the parameters of all compared methods are configured according to the settings reported in the corresponding references. For the proposed algorithm, the prior parameters are set as follows: ${\mathit{\beta }}_0={[2,2,2]}^{\mathrm{T}}$, ${\xi }_k=5$, ${\lambda }_k=1$, ${\mathit{\phi }}_k={[0,0,0]}^{\mathrm{T}}$, $\tau ={10}^{-3}$ and $\delta =0.05$. The initial target state is ${\mathit{x}}_0={[2000,-4,2000,-3,0]}^{\mathrm{T}}$, the state covariance matrix is ${\mathit{P}}_0=\mathrm{diag}\left(\left[{20}^2,5^2,{20}^2,5^2,{(0.3\mathrm{pi}/180)}^2\right]\right)$, The transition probability matrix is $\mathsf{\Pi }=\left[\begin{array}{c}0.980.02\hfill \\ 0.020.98\hfill \end{array}\right]$, the initial model probability is ${\mathit{\mu }}_0=\left[0.80.2\right]$, The number of iteration is set as $N_m=10$. All algorithms were implemented in MATLAB, and all simulation were performed on a computer equipped with an Intel Core i9-13900HX CPU @ 2.2 GHz and 64 GB of RAM.

Figure 6 shows the ground-truth trajectory of the maneuvering target and the locations of the three buoys. Figure 7 presents the range measurements obtained from the three buoys, and Figure 8 illustrates the probability density distributions of the measurement noise for the three buoys.

To evaluate the filtering performance, the root mean square error (RMSE) and the average RMSE (ARMSE) are employed as performance metrics. All results are obtained from 1000 independent Monte Carlo runs.

3.2. Parameter Sensitivity Analysis

3.2.1. Sliding-Window Length

Figure 9 illustrates the effect of the sliding-window length L on the position and velocity ARMSE. As the window length increases, the position ARMSE initially decreases markedly, while the velocity ARMSE decreases more gradually. The position ARMSE is lowest around $L=4$, with relatively low values over $L\in [3,5]$, whereas the velocity ARMSE reaches its minimum at $L=5$. Beyond $L=5$, both metrics show an overall upward trend, although the velocity ARMSE exhibits minor local fluctuations. These results indicate that an appropriate window length helps reduce tracking errors by fusing local estimates from several recent time steps, thereby mitigating the influence of fluctuations in individual estimates on subsequent filtering updates. However, an excessively long window may incorporate historical information that no longer adequately represents the current state or noise characteristics, reducing the responsiveness of the estimator and increasing estimation errors. Therefore, considering both position and velocity estimation accuracy, $L=5$ is selected as the optimal window length among the configurations examined and is used in the subsequent experiments.

3.2.2. Process Noise Amplification Factor

Table 3 presents the impact of the process noise amplification factor c on the tracking performance of the algorithm. The results indicate that, as c increases within the range [1, 100], the ARMSE of both position and velocity decreases and reaches its minimum at $c=100$. This phenomenon confirms the rationality of the proposed adaptive process noise covariance adjustment: when a target maneuver is detected, moderately increasing the process noise covariance can effectively increase the filter gain, thereby compensating for the lag induced by model mismatch and improving the responsiveness to maneuvering dynamics. However, when c is further increased to the range [100, 200] the ARMSE gradually rises, suggesting that an excessively large amplification factor makes the algorithm overly sensitive to measurement noise in steady-state phases, leading to unnecessary errors. In summary, to balance maneuver responsiveness and steady-state tracking accuracy, $c=100$ is selected as the optimal setting in the subsequent experiments.

3.3. Algorithm Ablation Experiments

To assess the contribution of each proposed module to tracking performance, an ablation study is conducted using five algorithm configurations. The complete AIMM-SWVBUF algorithm is denoted as Full, while IMM-VBUF, which incorporates none of the three proposed enhancements, is denoted as Baseline. Three variants—w/o Fusion, w/o MAP, and w/o Adapt—are obtained by individually removing sliding-window parameter fusion, sliding-window MAP estimation of the shape parameter, and adaptive process noise covariance adjustment from the complete algorithm, respectively. Table 4 summarizes the modules included in each configuration. Sliding-window parameter fusion, sliding-window MAP estimation of the shape parameter, and adaptive process noise covariance adjustment are described in Section 2.3.2, Section 2.3.3, and Section 2.5, respectively.

All algorithms are evaluated under the same simulation conditions and parameter settings. Figure 10 and Figure 11 present the position and velocity tracking RMSE curves, respectively, for the algorithms in the ablation study, while Table 5 reports the corresponding ARMSE values.

Compared with w/o Fusion, Full reduces the position and velocity ARMSE by approximately 26.3% and 4.8%, respectively. Sliding-window parameter fusion integrates local estimates from recent time steps, reducing the direct influence of fluctuations in individual estimates on subsequent filtering updates. These results indicate that the module contributes to lower tracking errors, with a more pronounced improvement in position estimation accuracy.

Compared with w/o MAP, Full reduces the position ARMSE by approximately 11.7%, while the velocity ARMSE remains essentially unchanged. Sliding-window MAP estimation of the shape parameter combines observations from multiple time steps within the window with the parameter prior to refine the prescribed shape parameter, allowing the noise model to adapt to the observations. These results suggest that the module primarily improves position tracking accuracy.

Compared with w/o Adapt, Full reduces the position and velocity ARMSE by approximately 9.6% and 7.2%, respectively, and also reduces the dominant error peaks following the onset of both turns. Adaptive process noise covariance adjustment responds to changes in model probabilities, improving the filter’s responsiveness to state changes during maneuvers and thereby mitigating tracking errors caused by motion model mismatch.

Overall, Full reduces the position and velocity ARMSE by approximately 32.3% and 19.6%, respectively, relative to Baseline. These results indicate that the combined use of the three enhancements reduces the average tracking error in the scenario considered.

3.4. Comprehensive Comparative Experiments

Figure 12 and Figure 13 show the RMSE curves of different algorithms for position and velocity, respectively, while Table 6 reports the corresponding average RMSE values.

Both IMM-UKF and AIMM-UKF are formulated under the Gaussian noise assumption, which does not fully capture the skewed heavy-tailed measurement noise considered in this study. AIMM-UKF employs a maneuver-adaptive model set to improve adaptability to changes in target motion. As shown in Table 6 AIMM-UKF achieves slightly lower position and velocity ARMSE values than IMM-UKF. However, both methods yield higher errors than the proposed algorithm under the considered noise conditions, indicating that improved motion-model adaptability alone may be insufficient to achieve robust tracking in the presence of skewed heavy-tailed measurement noise.

Although STUKF, ST-GHSST-KF and R-DVB-MHE-BG are designed for non-Gaussian noise and can suppress outliers to some extent, their performance is still limited by inherent bottlenecks. STUKF employs a prespecified degrees-of-freedom parameter, which limits its adaptability to changes in noise characteristics. Although ST-GHSST-KF models noise using Gaussian scale mixture distributions, its shape parameter estimates remain close to their initial values under the configuration considered, exhibiting only minor adjustments. This behavior may restrict the noise model’s ability to adapt to the observed measurements. R-DVB-MHE-BG combines a Bernoulli-Gaussian measurement model with moving horizon estimation to handle measurement outliers; however, its position and velocity ARMSE values remain higher than those of ST-GHSST-KF and the proposed algorithm.

In contrast, AIMM-SWVBUF achieves the best overall performance. Benefiting from the sliding-window mechanism that provides smoothed estimation of noise parameters and the adaptive process noise covariance adjustment that enables rapid response to maneuvers, the proposed method attains a superior balance between suppressing heavy-tailed noise and capturing strong maneuvers. As shown in Table 6, relative to the best-performing benchmark for each metric, the proposed algorithm reduces the position and velocity ARMSE by approximately 24.8% and 20.2%, respectively. These results demonstrate the high accuracy and strong robustness of AIMM-SWVBUF for maneuvering target tracking in complex underwater acoustic environments.

Table 7 compares the average execution times of three methods involving variational Bayesian iterations. To facilitate a fair comparison, convergence-based stopping criteria are also applied to the variational iterations of ST-GHSST-KF and R-DVB-MHE-BG. Under the same experimental conditions, AIMM-SWVBUF has a lower average execution time than the other two methods.

3.5. Validation Using Lake Trial Data

To validate the tracking performance of the proposed algorithm in a real underwater acoustic environment, an underwater target tracking experiment was conducted in June 2025 at an artificial lake. In the experiment, an acoustic transceiver was mounted beneath a small boat to serve as the underwater moving target, and a high-precision GPS positioning system was installed on the boat to provide the ground-truth target trajectory. Within an area of approximately $1{\mathrm{km}}^2$, three positioning buoys (labeled S1-S3) were deployed. Each buoy was equipped with an acoustic transceiver and a GPS receiver, forming a distributed observation network. The GPS positioning information was updated at 1 Hz. The co-located transmit-receive transducers on the buoys transmitted probing signals to the target boat and received echoes every 2s, with a sampling interval of $T=2\mathrm{s}$. All collected data were transmitted in real time to the mothership platform for logging and processing. The total sampling duration was 200s.

The experiment established a two-dimensional Cartesian coordinate system with the mothership as the origin. After coordinate transformation, the coordinates of the three positioning buoys are S1(-442.7085m, -192.0706m), S2(-783.0345, -579.1814m), S3(-38m, -318.8523m). For filter initialization, the process noise covariance is set to ${\mathit{Q}}_k=\mathrm{diag}\left(\left[0.02,0.02,0.01\right]\right)$ and the measurement noise scale matrix is $\mathit{R}=10{\mathit{I}}_3$. The GPS-measured initial position of the boat is (-230m, -253.08m). Accordingly, the initial state is set to ${\mathit{x}}_0=[-230,0,-253.08,0,0]$ and the initial state covariance matrix is ${\mathit{P}}_0=\mathrm{diag}\left(\left[5,2,5,2,0.3\right]\right)$. The remaining parameters are kept consistent with those used in the numerical simulations.

Figure 14 shows the ground-truth target trajectory and the buoy deployment layout, while Figure 15 presents the range measurements from the three buoys. By subtracting the GPS-derived geometric ranges from the acoustic range measurements at each time instant, the true measurement errors are obtained. Kernel density estimation (KDE) is then applied to the error sequence, producing the probability density curves shown in Figure 16. It can be observed that the measurement noise in the real underwater acoustic channel exhibits pronounced heavy-tailed and skewed characteristics, providing strong experimental evidence for adopting the GHSST distribution for noise modeling in this article.

Figure 17 and Figure 18 show the position and velocity estimation errors of different algorithms in the lake-trial experiment, respectively. The position average error values of the IMM-UKF, STUKF, ST-GHSST-KF, R-DVB-MHE-BG and AIMM-SWVBUF algorithms are 9.696 m, 9.398 m, 7.630 m, 10.681 m and 3.468 m, respectively, while the corresponding velocity average error values are 0.373 m/s, 0.412 m/s, 0.367 m/s, 0.308 m/s and 0.275 m/s. It can be seen that the proposed algorithm achieves the best performance in both metrics. Compared with the other algorithms, it reduces the average position error by more than 54.5% and the average velocity error by more than 10.7%. These results indicate that, in real underwater acoustic environments, the proposed algorithm can more effectively suppress the adverse effects of non-Gaussian skewed heavy-tailed noise and target maneuvering, thereby achieving higher tracking accuracy and robustness.

4. Conclusions

This paper proposes an adaptive interacting multiple-model sliding-window variational Bayesian unscented filtering algorithm to address fluctuations in local parameter estimates and insufficient updates of certain parameters under skewed heavy-tailed measurement noise, as well as motion model mismatch caused by target maneuvers. The method models measurement noise using a generalized hyperbolic skewed Student’s t distribution. Building on per-time-step updates, it applies sliding-window fusion to estimates of the location parameter and temporal smoothing to local MAP estimates of the mixing variable. It also updates the shape parameter through sliding-window MAP estimation and adaptively adjusts the process noise covariance based on changes in model probabilities. Simulation results demonstrate reductions in position and velocity tracking errors in the scenarios considered, while lake-trial results further validate the algorithm’s effectiveness and robustness in a complex underwater acoustic environment. Future work will investigate adaptive window-length selection, allowing the window length to vary with noise statistics and target motion to further improve parameter estimation stability and tracking performance.