Explainable Gated Bayesian Recurrent Neural Network for Non-Markov State
  Estimation
The optimality of Bayesian filtering relies on the completeness of prior
models, while deep learning holds a distinct advantage in learning models from
offline data. Nevertheless, the current fusion of these two methodologies
remains largely ad hoc, lacking a theoretical foundation. This paper presents a
novel solution, namely an explainable gated Bayesian recurrent neural network
specifically designed to state estimation under model mismatches. Firstly, we
transform the non-Markov state-space model into an equivalent first-order
Markov model with memory. It is a generalized transformation that overcomes the
limitations of the first-order Markov property and enables recursive filtering.
Secondly, by deriving a data-assisted joint state-memory-mismatch Bayesian
filtering, we design a Bayesian gated framework that includes a memory update
gate for capturing the temporal regularities in state evolution, a state
prediction gate with the evolution mismatch compensation, and a state update
gate with the observation mismatch compensation. The Gaussian approximation
implementation of the filtering process within the gated framework is derived,
taking into account the computational efficiency. Finally, the corresponding
internal neural network structures and end-to-end training methods are
designed. The Bayesian filtering theory enhances the interpretability of the
proposed gated network, enabling the effective integration of offline data and
prior models within functionally explicit gated units. In comprehensive
experiments, including simulations and real-world datasets, the proposed gated
network demonstrates superior estimation performance compared to benchmark
filters and state-of-the-art deep learning filtering methods.
Shi Yan, Yan Liang, Le Zheng, Mingyang Fan, Xiaoxu Wang, Binglu Wang
# Multi-level Gated Bayesian Recurrent Neural Network for State Estimation

###### Abstract

The optimality of Bayesian filtering relies on the completeness of prior models, while deep learning holds a distinct advantage in learning models from offline data. Nevertheless, the current fusion of these two methodologies remains largely ad hoc, lacking a theoretical foundation. This paper presents a novel solution, namely a multi-level gated Bayesian recurrent neural network specifically designed to state estimation under model mismatches. Firstly, we transform the non-Markov state-space model into an equivalent first-order Markov model with memory. It is a generalized transformation that overcomes the limitations of the first-order Markov property and enables recursive filtering. Secondly, by deriving a data-assisted joint state-memory-mismatch Bayesian filtering, we design a Bayesian multi-level gated framework that includes a memory update gate for capturing the temporal regularities in state evolution, a state prediction gate with the evolution mismatch compensation, and a state update gate with the observation mismatch compensation. The Gaussian approximation implementation of the filtering process within the gated framework is derived, taking into account the computational efficiency. Finally, the corresponding internal neural network structures and end-to-end training methods are designed. The Bayesian filtering theory enhances the interpretability of the proposed gated network, enabling the effective integration of offline data and prior models within functionally explicit gated units. In comprehensive experiments, including simulations and real-world datasets, the proposed gated network demonstrates superior estimation performance compared to benchmark filters and state-of-the-art deep learning filtering methods.

 state estimation, gated recurrent neural network, Bayesian filtering

## I Introduction

Bayesian filtering (BF) constructs a posterior probability density function (PDF) for the system state by leveraging available information [1], which is widely used in target tracking [2, 3, 4], navigation [5, 6], and localization [7] due to its solid theoretical foundation and effective computational design. During the filtering process, state-space models (SSMs) play a pivotal role in depicting state evolution and sensor observations. The optimality of traditional Bayesian filters depends on the correctness of prior SSMs [8], and hence their performance significantly deteriorates in the presence of model mismatches, while the complete and correct prior SSM is difficult to obtain in a complex and non-cooperative environment. Therefore, much attention has been paid to online Bayesian optimization, with representative approaches being the multi-model (MM) estimation and joint estimation identification (JEI). Its idea is the real-time inference of model parameters along with the recursive computation of density. The limited information available in online measurements necessitates certain prior conditions or assumptions for these parameters. Specifically, the MM estimation requires a prior model set and transition probability matrix [9, 10, 11]; JEI implemented by the expectation-maximization (EM) for parameter identification assumes that the estimated parameter is piecewise constant [12, 13, 14, 15]; and JEI implemented by the variational inference (VI) for covariance identification [16, 17, 18, 19] assume that the estimated covariance has a conjugate distribution. More importantly, the online Bayesian optimization approaches exclusively deal with real-time measurement data, while disregarding the abundant offline data.

With the widespread deployment of sensors and advancements in simulation technologies, there has been a notable increase in the available offline data, leading to a surge in the use of deep learning (DL). By capturing the temporal regularities in state evolution through their powerful nonlinear fitting capabilities and memory iteration mechanisms, various gated recurrent neural networks (RNNs), like long short-term memory (LSTM) [20] and gated recurrent unit (GRU) [21], have demonstrated remarkable performance in handling time series data [22, 23]. Gated RNNs are utilized to establish mapping functions from measurements to states, revealing the feasibility of extracting implicit modeling information from offline data [24, 25]. Nonetheless, pure deep learning approaches do not systematically integrate prior model knowledge, which necessitates a large number of learnable parameters and a substantial amount of labeled data, often resulting in a lack of interpretability.

To introduce prior model knowledge and improve interpretability, methods combining DL and BF have been widely studied in recent years. The RNN is utilized for dynamic model learning, i.e., it is integrated into the Bayesian filter as a mapping for the state transition [26, 27, 28], while this leads to the neglect of the prior dynamic model knowledge. This idea of utilizing RNN to directly learn the dynamic and measurement models of the Bayesian filter has more applications in the computer vision field [29, 30]. An improved approach employs the LSTM to correct the estimated state of a
model-based Bayesian filter [31]. When the model is severely mismatched, such an ad hoc combination faces failure as the filter diverges, and the covariance is not corrected so that it is inconsistent with the corrected state. Given the significance of filter gain in balancing prediction and measurement, a recent study utilized the GRU to estimate the gain [32]. By omitting second-order moment computations, this method accelerates the computation but sacrifices the utilization of prior noise statistical information and the ability to output covariance. In general, it is essential to derive an algorithmic framework based on a well-defined model and Bayesian filtering theory when introducing DL in the filtering. This is a crucial element currently lacking in existing methods, and it has the potential to enhance interpretability and estimation performance.

To this end, we propose a general state estimation problem for the non-Markov system, considering the non-Markov state evolution model, uncertain measurement model, and labeled offline data. For this problem, an underlying multi-level gated Bayesian RNN (MGBRNN) is designed. At first, it is discovered that the non-Markov SSM can be transformed into an equivalent first-order Markov model with memory through function nesting, which enables efficient recursive computation while overcoming the limitations of the first-order Markov property. Secondly, a data-assisted joint state-memory-mismatch BF framework is derived based on the transformed model, featuring three-level gated units: the memory update gate captures the temporal regularities in state evolution, and the state prediction and state update gates perform filtering while compensating for model mismatches. With its functionally explicit gated units, the framework efficiently integrates offline data and prior model knowledge with good interpretability. Considering computing efficiency, the framework is derived as a Gaussian approximation implementation. Finally, the internal neural network (NN) structure and end-to-end training method of each gated unit are designed. Notably, rather than artificially combing the DL and the Bayesian filter, the proposed MGBRNN is naturally derived from BF theory. In comprehensive experiments with benchmark simulations and real-world datasets, the MGBRNN outperforms various model-based (MB) filters and state-of-the-art filtering methods.

Notation: throughout this paper, the superscript \((\cdot)^{\top}\) denotes the transpose operation; \((\cdot)\) represents identical content to that stated in the preceding parenthesis; \(\mathrm{E}\left[\cdot|\cdot\right]\) denotes the conditional expectation.

## II Problem Formulation

The state estimation of dynamic systems typically involves the utilization of the following first-order Markov SSM with additive noises.

### System model:

\[\mathbf{x}_{k}=f_{k}\left(\mathbf{x}_{k-1}\right)+\mathbf{w}_{k} \tag{1}\]

### Measurement model:

\[\mathbf{z}_{k}=h_{k}\left(\mathbf{x}_{k}\right)+\mathbf{v}_{k} \tag{2}\]

where the subscript \(k\) represents the time index; \(\mathbf{x}_{k}\) is the state vector, such as the aircraft's position and velocity; \(\mathbf{z}_{k}\) is the measurement vector, such as the radial distance and azimuth of radar observation; \(f_{k}\) and \(h_{k}\) are the nominal linear/nonlinear state-evolution and state-measurement functions, respectively; the process noise \(\mathbf{w}_{k}\) and the measurement noise \(\mathbf{v}_{k}\) are independent of each other, and their PDFs are known and always being Gaussian distrubution.

The first-order Markov property of the system model in Eq. (1) has certain reasonable aspects. At first, the most recent state has the most effect on the evolution since the state evolution could be roughly described as an accumulative process. Secondly, the state transition process between adjacent frames is easy to analyze, making it possible to characterize the process with a first-order model. Finally, the first-order Markov property supports the recursive computation. The above advantages have led to the widespread application of this modeling form.

In fact, the real-world state evolution may be more complex than the model in Eq. (1). In practical state estimation tasks such as target tracking and navigation, the state evolution could be related to any of the historical states. Fig. 1 shows two representative examples of target tracking. In Fig. 1 (a), during target turning motion, determining the target's state at time \(k\) necessitates not only the state at time \(k-1\) but also parameters like turning radius and turning center, which are determined from state information across multiple frames. The traditional first-order Markov model, with no consideration of historical state, is difficult to get the turning parameters, resulting in modeling errors. Moreover, the target's motion involves a combination of recurring typical motion modes, allowing for the current state's determination based on historical states due to the invariance caused by the same motion mode. As shown in Fig. 1 (b), the invariant distance allows to obtain \(\mathbf{x}_{k}\) from \(\mathbf{x}_{k-N}\), but \(N\) is unknown. Such a non-Markov property is also prevalent in natural language processing. For

Fig. 1: **Examples of non-Markov state evolution model.** (a) The turning motion mode requires multiple frames of state to be recognized. (b) The current state can be obtained by analyzing historical states.
instance, in a sentence, the meaning of "bank" depends on the preceding context, with interpretations varying between a financial institution and the edge of a river or lake. Obviously, the state evolution could not be accurately depicted by the first-order Markov model, and it should be modeled in a non-Markov form [33], which is described as

\[\mathbf{x}_{k}=f_{k}^{\text{general}}\left(\mathbf{x}_{k-1},\mathbf{x}_{k-2}, \cdots,\mathbf{x}_{1}\right)+\mathbf{w}_{k} \tag{3}\]

where the general state-evolution function \(f_{k}^{\text{general}}\) reflects the effects of previous states and cannot be given a priori due to the increasing dimensionality. Note that Eq. (1) is a particular case of Eq. (3) when the state sequence in \(f_{k}^{\text{general}}\) contains only the previous single frame.

The measurement model in Eq. (2) is also uncertain in practical state estimation tasks, and the uncertainties mainly come from two aspects: 1) the model mismatch caused by observation anomalies, such as offsets or unknown rotations in the state observation function due to poor sensor calibration; 2) it is often necessary to approximate intractable complex mappings with computationally convenient functions, which causes errors, e.g., the linearization of the non-linear system. We define the actual measurement model as

\[\mathbf{z}_{k}=h_{k}^{\text{general}}\left(\mathbf{x}_{k}\right)+\mathbf{v}_ {k} \tag{4}\]

where the general state-measurement function \(h_{k}^{\text{general}}\) is unknown, which depicts the actual observation process.

For the filtering of the uncertain non-Markov SSMs presented in (3) and (4), we also consider access to labeled offline data. The abundant regularity information about state evolution and sensing observations is implicit in the offline data, which needs to be mined using time-series DL methods. Here, the available offline data consists of some sequence of observations and their corresponding ground truth states, i.e.,

\[\mathcal{D}=\left\{\mathbf{x}_{1:K}^{i},\mathbf{z}_{1:K}^{i}\right\}_{i=1 \cdots I} \tag{5}\]

where \(\mathbf{x}_{1:K}^{i}\) is the state ground-truth sequence obtained offline and \(\mathbf{z}_{1:K}^{i}\) is the corresponding measurement sequence; \(K\) is the duration of the system evolution; and \(I\) is the number of samples in the data set. Such an offline data set can be accumulated by collecting sensor readings or executing simulations. For instance, collecting location ground-truth and odometer readings for wheeled robots.

It is quite challenging for the filtering to SSMs in Eqs. (3) and (4) with considering the offline data set \(\mathcal{D}\). On the one hand, computational cost considerations require the derivation of efficient recursive estimation frameworks for non-Markov SSMs. On the other hand, it is required to extract regularity information about state evolution and sensing observations implicit in the offline data, which is used to learn uncertain SSMs. In addition, one needs to consider the prior information that can be available while utilizing the offline data, i.e., the prior nominal model in Eqs. (1) and (2). To this end, our goal is to construct a state estimation algorithm for uncertain non-Markov SSMs by synthesizing prior model information and offline data. This algorithm should have an efficient recursive estimation framework derived from BF theory.

## III Multi-level Gated Bayesian Recurrent Neural Network

The algorithm's design comprises the following three stages. At first, the non-Markov SSM is transformed into an equivalent first-order Markov model with memory. Secondly, a data-assisted joint state-memory-mismatch BF is derived based on the transformed model, leading to a multi-level gated framework. The filtering process within the gated framework is implemented using the Gaussian approximation. Finally, we design the structure of the NNs in each gated unit and the corresponding training method.

### _Problem transformation_

By introducing prior nominal models \(f_{k}\) and \(h_{k}\) in Eqs. (1) and (2), the problem of state estimation under uncertain and non-Markov models is transformed into the processing of the corresponding errors, thus we have

\[\mathbf{x}_{k}= f_{k}\left(\mathbf{x}_{k-1}\right)+\underbrace{f_{k}^{\text{ general}}\left(\mathbf{x}_{k-1},\cdots,\mathbf{x}_{1}\right)-f_{k}\left( \mathbf{x}_{k-1}\right)}_{\Delta_{k}^{f}}+\mathbf{w}_{k} \tag{6}\] \[\mathbf{z}_{k}= h_{k}\left(\mathbf{x}_{k}\right)+\underbrace{h_{k}^{\text{ general}}\left(\mathbf{x}_{k}\right)-h_{k}\left(\mathbf{x}_{k}\right)}_{ \Delta_{k}^{h}}+\mathbf{v}_{k} \tag{7}\]

where \(\Delta_{k}^{f}\) and \(\Delta_{k}^{h}\) are the corresponding state-dependent evolution mismatch error and observation mismatch error, respectively.

In Eq. (6), the non-Markov error \(\Delta_{k}^{f}\) is difficult to recognize directly since it continues to change as the states accumulate. Here, a two-layer nested function is constructed to approximate \(\Delta_{k}^{f}\), which is described as

\[\Delta_{k}^{f}=f_{k}^{\Delta}\left(\underbrace{g_{k}^{\Delta}\left(\mathbf{x} _{1},\cdots,\mathbf{x}_{k-1}\right)}_{\mathbf{c}_{k}}\right) \tag{8}\]

where \(g_{k}^{\Delta}\) encapsulates the effects of historical states, representing a regular memory pattern in state evolution, and hence denoted as memory \(\mathbf{c}_{k}\), which exhibits non-Markov characteristics; and \(f_{k}^{\Delta}\) is an output function. In complex state evolution, \(\mathbf{c}_{k}\) is essentially a hidden variable that is difficult to explicitly construct into a physical meaning, while the extraction regularity of the memory is implicit in a large amount of offline data.

For the memory \(\mathbf{c}_{k}\), considering the evolutionary regularities implicit in the states at different instants are similar, we also use the idea of nested functions to deal with it. Specifically, \(\mathbf{c}_{k}\) is approximated through the following nested function:

\[\mathbf{c}_{k}=\underbrace{g^{e}(g^{e}(g^{e}(\cdots),\mathbf{x}_{k-3}), \mathbf{x}_{k-2}),\mathbf{x}_{k-1})}_{k\text{ times}} \tag{9}\]

where, at each instant, \(g^{e}\) selectively retains information distilled from historical states and introduces information about new states. In fact, the function nesting in Eq. (9) reveals the idea of memory (or hidden variable) iteration in the LSTM
[20]. Such a nested form constructs an iterative update of \(\mathbf{c}_{k}\), which is equivalently expressed as

\[\mathbf{c}_{k}=g^{c}\left(\mathbf{c}_{k-1},\mathbf{x}_{k-1}\right) \tag{10}\]

For the state-dependent observation mismatch error \(\Delta_{k}^{h}\), we construct the mapping relation as

\[\Delta_{k}^{h}=h^{\Delta}\left(\mathbf{x}_{k}\right) \tag{11}\]

Through the above function approximation, we discover that the non-Markov system model can be equivalently transformed into a first-order Markov model with memory iteration. In this context, the recursion of memory is affected by the state, while the evolution of the state is also affected by the memory. If \(\mathbf{c}_{k}\) is considered an extended state, the new state composed of \(\mathbf{x}_{k}\) and \(\mathbf{c}_{k}\) remains Markov, demonstrating that non-Markov systems can be equivalently represented using the first-order Markov system with an appropriately extended state dimension. Such a model preserves historical state information while supporting recursive computation. Considering the offline data, the estimation of the model mismatch errors is also transformed into the learning of the functions \(f^{\Delta}\), \(g^{c}\) and \(h^{\Delta}\) from \(\mathcal{D}\), which is equivalent to learning the corresponding conditional PDFs. Unlike the traditional Bayesian filtering framework, a triple PDF recursive framework with state, memory, and mismatch is required for this model.

### _Bayesian multi-level gated framework design_

By introducing the offline data set \(\mathcal{D}\), we derive the joint state-memory-mismatch Bayesian filtering as shown in Theorem 1.

**Theorem 1**.: _(Joint state-memory-mismatch Bayesian filtering) For the SSM in Eqs. (6) - (8), (10), and (11), given the previous measurements \(\mathbf{z}_{1:k-1}\), the offline data set \(\mathcal{D}\), and the previous joint state-memory-mismatch posterior density \(p\left(\mathbf{x}_{k-1},\mathbf{c}_{k-1}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right)\), the joint state-memory density prediction is_

\[p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k-1}, \mathcal{D}\right.\right)\] \[=\iint P_{k}^{1}p\left(\mathbf{x}_{k-1},\mathbf{c}_{k-1}\left| \mathbf{z}_{1:k-1},\mathcal{D}\right.\right)d\mathbf{x}_{k-1}d\mathbf{c}_{k-1} \tag{12}\]

_with_

\[P_{k}^{1}= \int p\left(\mathbf{x}_{k}\left|\mathbf{x}_{k-1},\Delta_{k}^{f} \right.\right)P_{k}^{2}d\Delta_{k}^{f} \tag{13}\] \[P_{k}^{2}= p\left(\Delta_{k}^{f}\left|\mathbf{c}_{k},\mathcal{D}\right. \right)p\left(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1},\mathcal{ D}\right.\right) \tag{14}\]

_The joint state-memory density update is_

\[p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k}, \mathcal{D}\right.\right)= \frac{\int p\left(\mathbf{z}_{k}\left|\Delta_{k}^{h},\mathbf{x}_ {k}\right.\right)p\left(\Delta_{k}^{h}\left|\mathbf{x}_{k},\mathcal{D}\right. \right)d\Delta_{k}^{h}}{p\left(\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1}, \mathcal{D}\right.\right)}\] \[\times p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k-1}, \mathcal{D}\right.\right) \tag{15}\]

Proof.: See Appendix A. 

The probabilistic transfer of the joint state-memory-mismatch BF is shown in Fig. 2, and its iteration process contains three gated units: 1) the memory update gate (MUG) determines the current memory based on the previous state, previous memory, and offline data; 2) the state prediction gate (SPG) determines the joint state-memory prediction PDF; 3) the state update gate (SUG) determines the joint update PDF of the state and memory. To get the state prediction and posterior PDFs, the outputs of the state prediction and update gates are computed as

\[p\left(\mathbf{x}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right. \right)= \int p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k-1}, \mathcal{D}\right.\right)d\mathbf{c}_{k} \tag{16}\] \[p\left(\mathbf{x}_{k}\left|\mathbf{z}_{1:k},\mathcal{D}\right. \right)= \int p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k}, \mathcal{D}\right.\right)d\mathbf{c}_{k} \tag{17}\]

Such a multi-level gated framework utilizes prior models in the recursive filtering process while integrating offline data for learning the unknown PDFs, including the memory update density \(p(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1},\mathcal{D}\right.)\), the evolution mismatch density \(p(\Delta_{k}^{f}\left|\mathbf{c}_{k},\mathcal{D}\right.)\), and the observation mismatch density \(p(\Delta_{k}^{h}\left|\mathbf{x}_{k},\mathcal{D})\). It is required to design NN modules for the learning of these unknown distributions from offline data, while the analytical implementation of the filtering process can be derived separately in the case that these distributions are known. To fully introduce the proposed gated framework, the filtering implementation is first derived, with the corresponding network module to be designed in the next subsection.

The filtering process of the gated framework has various implementation methods, such as Monte Carlo sampling [34] and Gaussian approximation [35]. The gated framework of our MGBRNN is implemented by Gaussian approximation, considering that it is a fast and stable implementation method. In fact, constructing the internal weights or the output of the network in the form of a Gaussian distribution is also widespread in the field of Bayesian NNs [36, 37]. The assumptions required to perform the Gaussian approximation are given as follows.

Fig. 2: Diagram of the joint state-memory-mismatch Bayesian filtering
Assumption 1. The process noise \(\mathbf{w}_{k}\) and the measurement noise \(\mathbf{v}_{k}\) obey Gaussian distributions as \(\mathbf{w}_{k}\sim N\left(0,\mathbf{Q}_{k}\right)\) and \(\mathbf{v}_{k}\sim N\left(0,\mathbf{R}_{k}\right)\), respectively.

Assumption 2. The measurement prediction PDF is a Gaussian distribution with first- and second-order moments of \(\mathbf{\hat{z}}_{k|k-1}\) and \(\mathbf{P}_{k|k-1}^{x}\), respectively, namely,

\[p\left(\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right)=N\left( \mathbf{z}_{k};\mathbf{\hat{z}}_{k|k-1},\mathbf{P}_{k|k-1}^{x}\right) \tag{18}\]

Assumption 3. The state posterior PDF is a Gaussian distribution with first- and second-order moments of \(\mathbf{\hat{x}}_{k|k}\) and \(\mathbf{P}_{k|k}\) respectively, namely,

\[p\left(\mathbf{x}_{k}\left|\mathbf{z}_{1:k},\mathcal{D}\right)=N\left( \mathbf{x}_{k};\mathbf{\hat{x}}_{k|k},\mathbf{P}_{k|k}\right) \tag{19}\]

Assumption 4. \(\mathbf{c}_{k}\) obeys a Gaussian distribution with first- and second-order moments of \(\mathbf{\hat{c}}_{k}\) and \(\mathbf{P}_{k}^{c}\), respectively, namely,

\[p\left(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1}, \mathcal{D}\right)=N\left(\mathbf{c}_{k};\mathbf{\hat{c}}_{k},\mathbf{P}_{k}^ {c}\right) \tag{20}\]

Assumption 5. \(\Delta_{k}^{f}\) obeys a Gaussian distribution with first- and second-order moments of \(\hat{\Delta}_{k}^{f}\) and \(\mathbf{P}_{k}^{f}\), respectively; and \(\Delta_{k}^{h}\) obeys a Gaussian distribution with first- and second-order moments of \(\hat{\Delta}_{k}^{h}\) and \(\mathbf{P}_{k}^{h}\), respectively, namely,

\[p\left(\Delta_{k}^{f}\left|\mathbf{c}_{k},\mathcal{D}\right)=N \left(\Delta_{k}^{f};\hat{\Delta}_{k}^{f},\mathbf{P}_{k}^{f}\right) \tag{21}\] \[p\left(\Delta_{k}^{h}\left|\mathbf{x}_{k},\mathcal{D}\right. \right)=N\left(\Delta_{k}^{h};\hat{\Delta}_{k}^{h},\mathbf{P}_{k}^{h}\right) \tag{22}\]

_Remark 1:_General Gaussian filtering without mismatch usually assumes that \(\mathbf{w}_{k}\), \(\mathbf{v}_{k}\), \(p\left(\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right)\) and \(p\left(\mathbf{x}_{k}\left|\mathbf{z}_{1:k},\mathcal{D}\right)\) are Gaussian [35], which corresponds to the Assumptions 3, 4, and 5.

The Gaussian approximation implementation of the filtering process in the proposed multi-level gated framework is shown in Theorem 2.

**Theorem 2**.: _(Gaussian approximation implementation) Under Assumptions 1-5, the state and covariance prediction are_

\[\mathbf{\hat{x}}_{k|k-1}= \int f_{k}\left(\mathbf{x}_{k-1}\right)P_{k-1}^{x+}d\mathbf{x}_{k -1}+\hat{\Delta}_{k}^{f} \tag{23}\] \[\mathbf{P}_{k|k-1}= \int f_{k}\left(\mathbf{x}_{k-1}\right)[f_{k}\left(\mathbf{x}_{k -1}\right)]^{\top}P_{k-1}^{x+}d\mathbf{x}_{k-1}+\mathbf{P}_{k}^{f}\] \[- \mathbf{\hat{x}}_{k|k-1}\mathbf{\hat{x}}_{k|k-1}^{\top}+\mathbf{Q }_{k} \tag{24}\]

_and the state and covariance update are_

\[\mathbf{\hat{x}}_{k|k}= \mathbf{\hat{x}}_{k|k-1}+\mathbf{P}_{k|k-1}^{xz}\left(\mathbf{P} _{k|k-1}^{z}\right)^{-1}\left(\mathbf{z}_{k}-\mathbf{\hat{z}}_{k|k-1}\right) \tag{25}\] \[\mathbf{P}_{k|k}= \mathbf{P}_{k|k-1}-\mathbf{P}_{k|k-1}^{xz}\left(\mathbf{P}_{k|k -1}^{z}\right)^{-1}\left(\mathbf{P}_{k|k-1}^{xz}\right)^{\top} \tag{26}\]

_where_

\[\mathbf{\hat{z}}_{k|k-1}= \int h_{k}\left(\mathbf{x}_{k}\right)P_{k}^{x-}d\mathbf{x}_{k}+ \hat{\Delta}_{k}^{h} \tag{27}\] \[\mathbf{P}_{k|k-1}^{z}= \int h_{k}\left(\mathbf{x}_{k}\right)\left(h_{k}\left(\mathbf{x}_{k }\right)\right)^{\top}P_{k}^{x-}d\mathbf{x}_{k}+\mathbf{P}_{k}^{h}\] \[- \mathbf{\hat{z}}_{k|k-1}\mathbf{\hat{z}}_{k|k-1}^{\top}+\mathbf{R }_{k}\] (28) \[\mathbf{P}_{k|k-1}^{xz}= \int\mathbf{x}_{k}(h_{k}\left(\mathbf{x}_{k}\right))^{\top}P_{k}^{ x-}d\mathbf{x}_{k}-\mathbf{\hat{x}}_{k|k-1}\mathbf{\hat{z}}_{k|k-1}^{\top}\] (29) \[P_{k-1}^{x+}= N(\mathbf{x}_{k-1};\mathbf{\hat{x}}_{k-1|k-1},\mathbf{P}_{k-1|k-1})\] \[P_{k}^{x-}= N(\mathbf{x}_{k};\mathbf{\hat{x}}_{k|k-1},\mathbf{P}_{k|k-1})\]

Proof.: See Appendix B. 

In the above implementation process, Gaussian-weighted integrals are crucial. There are various methods for fast calculation of integrals in Gaussian approximate filtering [35]. Here, we give the computational process under the first-order Taylor expansion. Let \(\mathbf{F}_{k}\) and \(\mathbf{H}_{k}\) be the Jacobian matrices of \(f_{k}\) and \(h_{k}\), respectively. Based on the idea of function approximation used in extended KF (EKF) [38], Eqs. (23) and (24) are replaced by:

\[\mathbf{\hat{x}}_{k|k-1}= f_{k}\left(\mathbf{\hat{x}}_{k-1|k-1}\right)+\hat{\Delta}_{k}^{f} \tag{30}\] \[\mathbf{P}_{k|k-1}= \mathbf{F}_{k}\mathbf{P}_{k|k-1}\mathbf{F}_{k}^{\top}+\mathbf{Q}_{ k}+\mathbf{P}_{k}^{f} \tag{31}\]

and the Eqs. (27)-(29) are replaced by:

\[\mathbf{\hat{z}}_{k|k-1}= h_{k}\left(\mathbf{\hat{x}}_{k|k-1}\right)+\hat{\Delta}_{k}^{h} \tag{32}\] \[\mathbf{P}_{k|k-1}^{z}= \mathbf{H}_{k}\mathbf{P}_{k|k-1}\mathbf{H}_{k}^{\top}\text{ + }\mathbf{R}_{k}+\mathbf{P}_{k}^{h}\] (33) \[\mathbf{P}_{k|k-1}^{xz}= \mathbf{P}_{k|k-1}\mathbf{H}_{k}^{\top} \tag{34}\]

Under the Gaussian approximation, the learning of the unknown distributions is transformed into the estimation of their corresponding first- and second-order moments. It is necessary to construct the corresponding mappings by learning from offline data to realize the estimation. In the following, we design the corresponding NN structure for each gated unit in the three-level gated framework and train them as the corresponding estimation mappings via offline data.

### _Gated network structure design_

The overall structure of the proposed network is shown in Fig. 3, which contains six NN modules in its three gated units for estimating the first- and second-order moments of memory and model mismatch errors, respectively. The prior model information is introduced through the filtering formulas, and the offline data is used for the learning of the NN, while such a structure is rigorously derived based on the BF theory. Benefiting from this, the network modules in each gated unit do not need complex structures. Therefore, we construct the same network with a conventional form for these network modules, i.e., a two-layer NN with the hyperbolic tangent activation function, which is described as a function \(F\):

\[F\left(\mathbf{i},j\right)=\omega^{j2}\left(\tanh\left(\omega^{j1}\mathbf{i}+ \mathbf{b}^{j1}\right)\right)+\mathbf{b}^{j2} \tag{35}\]

where \(\mathbf{i}\) is the input vector and \(j\) is the index; \(\omega\) and \(\mathbf{b}\) are the weight and bias of the NN, respectively, and \(\omega,\mathbf{b}\in\mathbf{\Phi}\); and tanh denotes the hyperbolic tangent activation function.

The inputs to the NN modules of each gating unit are introduced as follows, with descriptions of some of the involved notations provided here: \(\sigma\) denotes the sigmoid function, \(con[\mathbf{a},\mathbf{b}]\) denotes the concatenation of vectors \(\mathbf{a}\) and \(\mathbf{b}\), flat denotes the flattening of the matrix into vectors, and \(\Psi\) denotes the maximal normalization operation. The normalization, flattening, and other operations we conducted are standard procedures in DL.

#### Iii-C1 Input to the NN in MUG

In the MUG, NN modules are utilized for the estimation of \(\mathbf{\hat{c}}_{k}\) and \(\mathbf{P}_{k}^{c}\). Therefore, based on the conditions of PDF \(p\left(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1},\mathcal{D}\right)\), the corresponding input vector is constructed as \(\mathbf{\hat{l}}_{k}^{c}=\text{con}\left[\sigma\left(\text{con}\left[\mathbf{ \hat{c}}_{k-1},\text{flat}\left(\mathbf{P}_{k-1}^{c}\right)\right]\right), \Psi\left(\mathbf{\hat{x}}_{k-1|k-1}\right)\right]\).
#### Iii-C2 Input to the NN in SPG

In the SPG, NN modules are utilized for the estimation of \(\hat{\Delta}_{k}^{f}\) and \(\mathbf{P}_{k}^{f}\). Therefore, based on the conditions of the PDF \(p(\Delta_{k}^{f}\,|\,\mathbf{c}_{k},\mathcal{D})\), the corresponding input vector is constructed as \(\mathbf{i}_{k}^{\mathbf{f}}=\sigma\left(\text{con}\left[\mathbf{\hat{c}}_{k}, \text{flat}\left(\mathbf{P}_{k}^{\mathbf{c}}\right)\right]\right)\).

#### Iii-C3 Input to the NN in SUG

In the SUG, NN modules are utilized for the estimation of \(\hat{\Delta}_{k}^{h}\) and \(\mathbf{P}_{k}^{h}\). Therefore, based on the conditions of the PDF \(p\left(\Delta_{k}^{h}\,|\,\mathbf{x}_{k},\mathcal{D}\right)\), the corresponding input vector is constructed as \(\mathbf{i}_{k}^{\mathbf{i}}=\Psi\left(\mathbf{\hat{x}}_{k|k-1}\right)\).

Based on the structures and inputs of the gated units described above, their internal NN modules are constructed as in Tab. I. Obviously, each gated unit has its own specified function, and the errors caused by model mismatches are targeted to be compensated.

Compared to the existing algorithmic frameworks that incorporate gated RNNs like LSTM and GRU into Bayesian filters, MGBRNN is an underlying gated RNN tailored for state estimation and designed based on BF theory. The gated units in MGBRNN possess finely partitioned functionality, rendering them physically interpretable and potentially reducing the demand for learnable parameters, thus diminishing the need for a large amount of training data. Despite the extensive research on task-specific underlying gated RNNs in recent years [39, 40, 41, 42], gated RNNs tailored specifically for state estimation remain largely unexplored, with MGBRNN introducing a novel perspective in this research direction.

### _Algorithm training_

As depicted in Fig. 3, in a supervised manner, the gated network within MGBRNN undergoes end-to-end training by extracting the underlying regularities governing state evolution and observation from \(\mathcal{D}\). Despite the fact that the NN computes intermediate parameters instead of directly generating estimated states, we employ the loss function computed from the estimated states and their corresponding ground-truth in \(\mathcal{D}\) to simultaneously train all the learnable parameters within MGBRNN. This training approach has been demonstrated to be viable in [32].

In line with the commonly used minimum mean square error (MMSE) criterion in BF, we employ the mean square error (MSE) loss function with the addition of L2 regularization for training MGBRNN. Consider the offline data set described in Eq. (5), the MSE loss is described as

\[l_{i}\left(\Phi\right)=\frac{1}{K}\sum_{k=1}^{K}\left\|\mathbf{\hat{x}}_{k|k} ^{i}\left(\Phi\right)-\mathbf{x}_{k}^{i}\right\|^{2}+\tau\left\|\Phi\right\|^ {2} \tag{36}\]

Fig. 3: **The overall structure of the MGBRNN and its training process. The multi-level gated structure with the memory update gate, the state prediction gate, and the state update gate is rigorously derived from Bayesian filtering theory and integrates prior model knowledge and offline data efficiently.**
where \(\tau\) is a coefficient used to adjust the proportion of the regular term; \(i\in\{1,2,\cdots,I\}\) represents the index of the training sample; and \(\Phi\) is the learnable parameter set. A variation of mini-batch stochastic gradient descent is utilized to optimize \(\Phi\). Specifically, for each batch indexed by \(n\), we select \(J<I\) trajectories denoted as \(i_{1}^{n},i_{2}^{n},\ldots,i_{J}^{n}\). The mini-batch loss is then computed as follows:

\[L_{n}\left(\Phi\right)=\frac{1}{J}\sum_{j=1}^{J}l_{j^{-}_{\tau}}\left(\Phi\right) \tag{37}\]

Given that MGBRNN features a recursive architecture, the backpropagation through time (BPTT) algorithm [43] is utilized for its training, which is a technique commonly used in training RNNs. The BPTT algorithm efficiently captures temporal dependencies by unfolding the network over time and computing gradient estimates in both forward and backward directions, making it effective for handling sequential data. Concretely, we utilize shared network parameters to unfold MGBRNN over time and calculate the forward and backward gradient estimations across the network.

## IV Experiments

In this section, we comprehensively evaluate the MGBRNN in four state estimation tasks, which encompass two benchmark simulation scenarios and two real-world data sets. At first, we conduct an ablation experiment on a real-world data set to illustrate the functions and effectiveness of each gated unit in the MGBRNN. Subsequently, we compared the MGBRNN with MB- and DL-based representative methods across all four tasks. This comparison showcases the outstanding performance and the lightweight parameter scale of the MGBRNN.

Throughout the experiment, we employ two units to quantify estimation errors: MSE in [dB], commonly used in signal processing, is utilized for simulation scenarios to narrow down the range of values; and root MSE (RMSE), a standard metric in target tracking and navigation, is utilized for the real-world data sets. Here, we describe the calculation for these two units. The MSE of the estimated state is calculated as

\[M_{\mathbf{x}}=\frac{1}{KN}\sum_{i=1}^{N}\sum_{k=1}^{K}\left(\mathbf{x}_{i,k}- \mathbf{\hat{x}}_{i,k|k}\right)\left(\mathbf{x}_{i,k}-\mathbf{\hat{x}}_{i,k|k }\right)^{\top} \tag{38}\]

where \(N\) is the number of test samples and \(K\) is the time duration. The MSE in [dB] is calculated as \(10\times\log_{10}\left(M_{\mathbf{x}}\right)\) and the RMSE is calculated as \(\sqrt{M_{\mathbf{x}}}\).

### _Ablation Experiment_

The Michigan NCLT data set is used to evaluate each gated unit of MGBRNN in the ablation experiment. The labeled trajectory data of the wheeled robots is gathered in the NCLT data set, which contains the robots' ground-truth locations and noisy sensor readings (e.g., GPS and odometer). We follow the setup in [32] to perform robot localization using measurements from these noisy readings. Our specific goal is to track the two-dimensional position using only odometry data, which means that the observations come from noisy velocity readings while modeling information is not available in the dataset. Odometer-based localization is available in locations with weak GPS signals or dense obstructions, such as inside or in canyons, but sensor reading drift and the inability to directly monitor positions make it difficult.

**Experimental setups.** In the NCLT data set, the trajectory at 2012-01-22 with outliers removed is segmented into 48 trajectories, each containing a sequence of ground-truth positions and corresponding noisy velocity readings with \(K=100\) and a sampling frequency of \(1\mathrm{\,Hz}\). In these trajectories, we set 34 for training (\(71\%\)), 9 for validation (\(19\%\)), and 5 for testing (\(10\%\)). The state vector is \(\mathbf{x}_{k}=[x_{k},y_{k},\dot{x}_{k},\dot{y}_{k}]^{\top}\) with \([x_{k},y_{k}]^{\top}\) and \([\dot{x}_{k},\dot{y}_{k}]^{\top}\) being the position and velocity, respectively, and the measurement is the noisy velocity. This data set contains only state ground-truth and corresponding measurements, so one need to construct the nominal SSM manually. The prior nominal state-evolution matrix \(\mathbf{F}^{\text{prior}}\) is a constant velocity (CV) model with the sampling interval \(\Delta t=1\), and the nominal state-measurement matrix \(\mathbf{H}^{\text{prior}}\) is a linear model, i.e.,

\[\mathbf{F}^{\text{prior}}=\begin{bmatrix}1&0&\Delta t&0\\ 0&1&0&\Delta t\\ 0&0&1&0\\ 0&0&0&1\end{bmatrix},\mathbf{H}^{\text{prior}}=\begin{bmatrix}0&0&1&0\\ 0&0&0&1\end{bmatrix}^{\top}\]

The process noise covariance and measurement noise covariance are computed as

\[\mathbf{Q}=q^{2}\mathbf{I}_{d_{x}},\quad\mathbf{R}=r^{2}\mathbf{I}_{d_{x}} \tag{39}\]

where we set nominal values of \(q=0.1\) and \(r=1\).

To evaluate the effectiveness of each gated unit, we set four schemes by masking the different gated units in MGBRNN: 1) full MGBRNN with three gated units; 2) MGBRNN without the MUG; 3) MGBRNN without the SPG; and 4) MGBRNN without the SUG. In Scheme 2, the input of the SPG does not contain the memory but only the state estimation of the previous instant.

**Experimental results and analysis** Fig. 4 depicts the convergence of localization errors on the validation set during the training process of the four schemes. MGBRNN with a full-gated structure exhibits the quickest convergence to achieve the lowest localization error. Although the scheme without the MUG converges quickly, it is constrained by the first-order Markov property, resulting in greater localization error. The scheme without the SUG cannot compensate for the observation mismatch, leading to slow convergence and a high loss. The scheme without the SPG trains the worst since learning the dynamic model is the most critical of state estimation.

The localization RMSE of the four schemes on the test set is shown in Tab. II, and the localization results of a single test sample are displayed in Fig. 5. Masking out each gated unit will increase localization errors, reflecting the effectiveness of each gated unit in MGBRNN.

### _Linear System_

In this scenario, the MGBRNN is compared with the Kalman filter (KF) and LSTM under a linear SSM to evaluate
its performance in response to model mismatch under linear conditions.

**Experimental setups.** We follow the setup in [32] to estimate the state of a linear system as

\[\mathbf{x}_{k}= \mathbf{F}\mathbf{x}_{k-1}+\mathbf{w}_{k} \tag{40}\] \[\mathbf{z}_{k}= \mathbf{H}\mathbf{x}_{k}+\mathbf{v}_{k} \tag{41}\]

where

\[\mathbf{F}\in \mathbb{R}_{d_{x}\times d_{x}},\quad F_{ij}=\begin{cases}1,& \text{if }i=j\text{ or }i=1\\ 0,&\text{else}\end{cases} \tag{42}\] \[\mathbf{H}\in \mathbb{R}_{d_{x}\times d_{z}},\quad H_{ij}=\begin{cases}1,& \text{if }i=d_{x}-j+1\text{ or }i=1\\ 0,&\text{else}\end{cases} \tag{43}\]

with \(d_{x}\times d_{z}=2\times 2\) and \(K=20\).

In the filtering process, our MGBRNN only acquires the mismatched model, i.e., the rotated state evolution matrix \(\mathbf{F}_{\theta}\) and the rotated state observation matrix \(\mathbf{H}_{\theta}\) with \(\theta=10^{\circ}\), which are described as

\[\mathbf{H}_{\theta}=\mathbf{T}_{\theta}\mathbf{H},\quad\mathbf{F}_{\theta}= \mathbf{T}_{\theta}\mathbf{F},\quad\mathbf{T}_{\theta}=\left[\begin{array}{ cc}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{array}\right] \tag{44}\]

The process noise covariance and the measurement noise covariance are calculated as Eq. (39), and the same \(q\) and \(r\) are used for data generation and filtering. Under several inverse observation noise levels as \(\frac{1}{r^{2}}\left[dB\right]\in\{-10,0,10,20\}\), we compare the MGBRNN with three methods: a single-layer LSTM directly estimates \(\mathbf{x}_{k}\) by input \(\mathbf{z}_{k}\); a KF with the accurate model; and a KF with the mismatched model.

**Experimental results and analysis.** The test result is shown in Fig. 6. When the model is mismatched, the KF without learning capability causes an extreme estimation error. Although LSTM learns state estimation by mining offline data, it lacks a specific learning architecture and cannot introduce prior models, leading to its poor estimated performance. The effect of the model on the estimation performance is particularly significant at lower observation noise, at which time the estimation performance of the LSTM is degraded. The estimation error of MGBRNN is close to MMSE (KF with the accurate model), demonstrating that it effectively compensates for model mismatches by learning from offline data.

### _Lorenz Attractor_

The Lorenz attractor is a chaotic system generated by a set of differential equations used to simulate atmospheric convection, characterized by ever-changing three-dimensional trajectories. In this scenario, the tracking ability of the MGBRNN in complex and highly nonlinear settings under model mismatch is validated.

**Experimental setups.** The Lorentz attractor scenario here has the same parameter settings as in [32] and [44]. The Lorenz attractor is mathematically described by a system of three coupled ordinary differential equations, and the noiseless state-evolution of the continuous time process

Fig. 4: **Training loss on the validation set.** The absence of each gate affects the convergence speed and final convergence result of MGBRNN training.

Fig. 5: **Localization results on test data with different ablation schemes.** The absence of each gate affects the localization performance of MGBRNN.

Fig. 6: Test results for linear systems
\(\left[x_{t}^{1},x_{t}^{2},x_{t}^{3}\right]^{\top}\) is

\[\frac{\partial\mathbf{x}_{t}}{\partial t}=\mathbf{A}\left(\mathbf{x}_{t}\right) \mathbf{x}_{t},\mathbf{A}\left(\mathbf{x}_{t}\right)=\left[\begin{array}{ccc} -\sigma&\sigma&0\\ \beta&-1&-x_{t}^{1}\\ 0&x_{t}^{1}&-\rho\end{array}\right] \tag{45}\]

where \(\sigma=10\), \(\rho=8/3\), and \(\beta=28\) are parameters governing the system's dynamics. We follow the processing in [32] to get a discrete-time state-evolution model, and the resulting discrete-time evolution process is \(\mathbf{x}_{k+1}=\mathbf{F}\left(\mathbf{x}_{k}\right)\mathbf{x}_{k}\), where \(\mathbf{x}_{k}=\left[x_{k}^{1},x_{k}^{2},x_{k}^{3}\right]^{\top}\) is a discrete-time state vector, and the discrete-time state-evolution matrix \(\mathbf{F}\left(\mathbf{x}_{k}\right)\) is obtained by the Taylor series expansion and a finite series approximation (with \(J\) coefficients). The covariance of the process noise and the measurement noise is computed as in Eq. (39).

Similar to the linear system scenario, we consider both evolution and observation model mismatches. Here, the observation model is an identity matrix rotated by merely \(\theta=10^{\circ}\), and the evolution model is \(\mathbf{F}\left(\mathbf{x}_{k}\right)\) with \(J=1\) while the actual data is generated by \(\mathbf{F}\left(\mathbf{x}_{k}\right)\) with \(J=5\). Under different noise levels of \(1/r^{2}\left[\text{dB}\right]\in\left\{-10,0,10\right\}\), our MGBRNN is compared with several benchmark nonlinear filters, i.e., the EKF, unscented KF (UKF), and particle filter (PF).

**Experimental results and analysis.** The tracking MSE in [dB] of each method is shown in Tab. III. MGBRNN exhibits the lowest estimation error, highlighting its capacity to significantly improve performance in scenarios with severe nonlinearities through offline learning. Other MB nonlinear filters all degrade in performance under model mismatch, especially under lower measurement noise. The tracking trajectories of each method for a noise level of \(1/r^{2}=0\text{[dB]}\) are shown in Fig. 7, with the trajectory estimated by MGBRNN being more closely fitting the ground-truth.

### _Arceraft trajectory tracking_

In this experiment, our goal is to track the real trajectories of civil aircraft landing at the airport. The tracking of landing aircraft is an important aspect of airspace surveillance, which is critical to flight safety and aircraft scheduling. Due to their complex maneuvers, tracking airborne targets presents a challenging state estimation task. Here, we focus on compensation for uncertain motion modes.

**Experimental setups.** The civil aviation trajectory data comes from [45], which contains the flight records covering Jan-Mar 2006 in Northern California TRACON. We extracted the trajectories of aircraft landing in the east direction on a designated runway and performed cubic spline interpolation on the trajectories to achieve a fixed sampling interval of \(\Delta t=4\). The trajectory length is set as \(K=200\)s, and shorter trajectories are discarded or longer ones are truncated accordingly. This resulted in 260 trajectories, with 182 (\(70\%\)) for training, 52 (\(20\%\)) for validation, and 26 (\(10\%\)) for testing. We consider tracking targets in the x-y two-dimensional coordinate system, and the form of the target state is the same as that in IV-A. The measurement model used corresponds to a general radar observation, i.e., the measurement \(\mathbf{z}_{k}\!=\!\!\left[d_{k},\mu_{k}\right]^{\top}\) with the radial distance \(d_{k}\) and the azimuth angle \(\mu_{k}\), the state-to-measurement function is

\[h(\mathbf{x}_{k})\!=\!\!\left[\sqrt{x_{k}^{2}+{y_{k}}^{2}},\arctan\left(\frac{y _{k}}{x_{k}}\right)\right]^{\top} \tag{46}\]

and the observation noise \(\mathbf{v}=\left[v_{d},v_{\mu}\right]^{\top}\) with \(v_{d}\sim\mathcal{N}\left(v_{d};0,\sigma_{d}{}^{2}\right)\) and \(v_{\mu}\sim\mathcal{N}\left(v_{\mu};0,\sigma_{\mu}{}^{2}\right)\).

Under various radar accuracy settings of \(\sigma_{\mu}\) and \(\sigma_{d}\), the MGBRNN is compared with various representative filtering methods under various radar accuracies. For setups of these methods, the MGBRNN, EKF, UKF, and PF adopt the CV model, while the model set of the IMM estimator includes one CV model and two coordinate turn (CT) models with \(\omega=\pm 3^{\circ}/s\). Another method is a combination method that employs LSTM to real-time compensate for the estimated state of EKF. The process noise covariance \(\mathbf{Q}\) is computed as Eq. (39) with \(q^{2}=1\) for IMM, and \(q^{2}=10\) for other methods. Since airport air traffic control radars are typically cooperative sensors with known accuracy, the measurement model is known, and the measurement noise of each filter is given as \(\mathbf{R}=\left[\begin{array}{cc}\sigma_{d}^{2}&0\\ 0&\sigma_{\mu}^{2}\end{array}\right]\).

**Experimental results and analysis.** The mean position and velocity RMSE of each method on the test set are shown in Tabs. IV and V, respectively. It can be seen that the MGBRNN has the lowest tracking error in both position and velocity. In comparison to the MB method, MGBRNN is capable of extracting implicit state evolution modes from offline data, i.e., multiple aircraft motion modes during landing. In comparison to the combination method, the gated structure of MGBRNN is tailored to state estimation under model mismatches, which adequately compensates for tracking errors caused by model mismatches.

Fig. 8 shows the tracking trajectories on some samples of the test set. MGBRNN provides accurate tracking during target maneuvering (see zoomed-in section in figure), demonstrating

Fig. 7: Tracking result of Lorenz attractor under model mismatch
that it effectively mines the implicit target motion regularities in the offline data.

### _Michigan NCLT Data Set_

On the Michigan NCLT data set mentioned in the ablation experiment, we compare MGBRNN with several representative filtering algorithms to demonstrate its performance and the small learnable parameter scale.

**Experimental setups.** With the same data setups as the ablation experiment, we compare MGBRNN with four types of representative methods, namely, the MB KF, the pure DL method LSTM, the state-of-the-art online Bayesian optimization method variational-Bayesian-based adaptive KF (VBKF), and the state-of-the-art DL-based filtering method KalmanNet. In these methods, KF and MGBRNN use the same model and noise settings as in the ablation experiments; LSTM has a single-layer structure and directly maps velocity measurements to positions; and the settings of VBKF and KalmanNet are referenced to [17] and [32], respectively.

**Experimental results and analysis.** Tab. VI shows the localization RMSE of each method on the test set. KF without learning capability has the highest localization error, and VBKF slightly outperforms KF by online inference to the covariance. KalmanNet performs offline learning while introducing prior model knowledge, but it only modifies the filter gain, making it incapable of compensating for model mismatches effectively. The localization trajectories of some test samples are shown in Fig. VI, where we can see that MGBRNN consistently maintains high localization accuracy with low termination error. In MGBRNN, each gated unit has its own explicit function, which is targeted to compensate for the model mismatch, so it achieves the lowest localization error in such a challenging real-world localization task.

In Fig. 10, we present the localization RMSE results for MGBRNN and KalmanNet with varying numbers of learnable parameters. It can be seen that MGBRNN outperforms KalmanNet, achieving lower localization errors despite having

Fig. 8: **Tracking results on the test data.** We show the tracking performance of different filtering methods under various radar accuracies, with these black boxes indicating the target maneuver instants. When the target undergoes maneuvers, the proposed MGBRNN responds promptly and attains high tracking accuracy.
fewer learnable parameters. The MGBRNN is an underlying gated RNN with a computationally efficient gated structure, which allows it to excel in state estimation tasks while maintaining simplicity. In fact, in such a task with very limited training samples, the small parameter scale of MGBRNN confers a significant advantage.

## V Conclusion

In this paper, a multi-level gated Bayesian RNN is specifically designed for state estimation under model mismatch, i.e., MGBRNN. The gated structure of MGBRNN is designed by introducing the memory mechanism and offline data into the BF theory, which provides good interpretability, low complexity, and data efficiency. Each gated unit in the multi-level gated framework has its own explicit functionality, which effectively compensates for the errors resulting from model mismatches by integrating prior models and offline data. For computational efficiency, the Gaussian approximation is used to implement the filtering process of the gated framework, along with the development of an internal neural network structure and the corresponding end-to-end training method. The MGBRNN has been verified by both simulations and real-world data sets. It is demonstrated that the proposed algorithm outperforms the state-of-the-art DL-based filtering methods and could achieve high estimation accuracy with a small number of learnable parameters.

## Appendix A Proof of the joint state-memory-mismatch BF

Denote that \(p\left(\mathbf{x}_{1:K},\mathbf{c}_{1:K}\left|\mathbf{z}_{1:K},\mathcal{D}\right.\right)\) is the joint posterior density of 1 to \(K\) time states and memories. Expanding this density under the Bayesian theorem as

\[p \left(\mathbf{x}_{1:K},\mathbf{c}_{1:K}\left|\mathbf{z}_{1:K}, \mathcal{D}\right.\right)\] \[=p\left(\mathbf{x}_{1:K},\mathbf{c}_{1:K}\left|\mathbf{z}_{K}, \mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[=\frac{p\left(\mathbf{z}_{K}\left|\mathbf{x}_{1:K},\mathbf{c}_{1 :K},\mathbf{z}_{1:K-1}\right.\right)p\left(\mathbf{x}_{1:K},\mathbf{c}_{1:K} \left|\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)}{p\left(\mathbf{z}_{K} \left|\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)} \tag{47}\]

According to Eq. (11), the measurement is independent of the previous state given the current state. In addition, the state information contained in previous measurements has been accurately characterized by the current state. Then we have

\[p\left(\mathbf{z}_{K}\left|\mathbf{x}_{1:K},\mathbf{c}_{1:K}, \mathbf{z}_{1:K-1},\mathcal{D}\right.\right)=p\left(\mathbf{z}_{K}\left| \mathbf{x}_{K},\mathcal{D}\right.\right) \tag{48}\]

According to the total probability formula, we have

\[p\left(\mathbf{x}_{1:K},\mathbf{c}_{1:K}\left|\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[=p\left(\mathbf{x}_{K},\mathbf{c}_{K}\left|\mathbf{x}_{K-1}, \mathbf{c}_{K-1},\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[\times p\left(\mathbf{x}_{1:K-1},\mathbf{c}_{1:K-1}\left|\mathbf{ z}_{1:K-1},\mathcal{D}\right.\right) \tag{49}\]

Thus, (A) is further written as

\[p\left(\mathbf{x}_{1:K},\mathbf{c}_{1:K}\left|\mathbf{z}_{1:K},\mathcal{D}\right.\right)\] \[=\frac{p\left(\mathbf{z}_{K}\left|\mathbf{x}_{K},\mathcal{D} \right.\right)p\left(\mathbf{x}_{K},\mathbf{c}_{K}\left|\mathbf{x}_{K-1}, \mathbf{c}_{K-1},\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)}{p\left(\mathbf{z} _{K}\left|\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)}\] \[\times p\left(\mathbf{x}_{1:K-1},\mathbf{c}_{1:K-1}\left|\mathbf{ z}_{1:K-1},\mathcal{D}\right.\right) \tag{50}\]

For Eq. (50), the conditional independence expressed by Eqs. (6) and (10) leads to

\[p\left(\mathbf{x}_{K},\mathbf{c}_{K}\left|\mathbf{x}_{K-1}, \mathbf{c}_{K-1},\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[=p\left(\mathbf{x}_{K}\left|\mathbf{c}_{K},\mathbf{x}_{K-1}, \mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[\times p\left(\mathbf{c}_{K}\left|\mathbf{x}_{K-1},\mathbf{c}_{K-1},\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[=\int p\left(\mathbf{x}_{K}\left|\mathbf{x}_{K-1},\Delta_{k}^{f}, \mathcal{D}\right.\right)P_{K-1}^{3}d\Delta_{k}^{f} \tag{51}\]

Fig. 10: The localization RMSE (y-axis) comparison of MGBRNN and KalmanNet on the test set under different numbers of learnable parameters (x-axis).

Fig. 9: **Localization results on the NCLT test data.** We present the localization results of the MGBRNN alongside several representative methods, which include KF, VBKF, and KalmanNet. MGBRNN demonstrates the highest accuracy with stable localization.
For \(P_{K-1}^{3}\), considering the transfer relations in Eq. (8) and Eq. (6) yields

\[P_{K-1}^{3}= p\left(\Delta_{k}^{f},\mathbf{c}_{K}\left|\mathbf{x}_{K-1}, \mathbf{c}_{K-1},\mathbf{z}_{1:K-1},\mathcal{D}\right.\right)\] \[= p\left(\Delta_{k}^{f}\left|\mathbf{c}_{K},\mathcal{D}\right. \right)p\left(\mathbf{c}_{K}\left|\mathbf{c}_{K-1},\mathbf{x}_{K-1},\mathbf{z }_{1:K-1},\mathcal{D}\right.\right) \tag{52}\]

Based on Chapman-Kolmogorov equations for Eq. (51), the joint state-memory prediction density \(p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right)\) can be obtained as Eqs. (12) - (14).

Based on the Bayesian theorem, the joint state-memory update density can be obtained as

\[p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k}, \mathcal{D}\right.\right)=\frac{p\left(\mathbf{z}_{k}\left|\mathbf{x}_{k}, \mathcal{D}\right.\right)p\left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z }_{1:k-1},\mathcal{D}\right.\right)}{p\left(\mathbf{z}_{k}\left|\mathbf{z}_{1 :k-1},\mathcal{D}\right.\right)} \tag{53}\]

From Eq.(7) it follows that \(\mathbf{z}_{k}\) is only related to the current state \(\mathbf{x}_{k}\) and not to the memory \(\mathbf{c}_{k}\), thus

\[p\left(\mathbf{z}_{k}\left|\mathbf{x}_{k},\mathcal{D}\right.\right)= p\left(\mathbf{z}_{k}\left|\mathbf{x}_{k},\mathcal{D}\right.\right)\] \[= \!\!\int p\left(\mathbf{z}_{k}\left|\Delta_{k}^{h},\mathbf{x}_{k} \right.\right)p\left(\Delta_{k}^{h}\left|\mathbf{x}_{k},\mathcal{D}\right. \right)d\Delta_{k}^{h} \tag{54}\]

## Appendix B Proof of the Gaussian approximation implementation

The mean of the state prediction is calculated as:

\[\mathbf{\hat{x}}_{k|k-1}= \mathrm{E}\left[\mathbf{x}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right]\] \[= E\left[\bar{f}_{k}\left(\mathbf{x}_{k-1}\right)+\Delta_{k}^{f}+ \mathbf{w}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right]\] \[= \!\!\int\!\!\!\int\!\!\!\int\left[\bar{f}_{k}\left(\mathbf{x}_{k -1}\right)+\Delta_{k}^{f}\right]P_{k}^{1}d\Delta_{k}^{f}d\mathbf{c}_{k}d \mathbf{c}_{k-1}d\mathbf{x}_{k-1} \tag{55}\]

with

\[P_{k}^{1}= p\left(\Delta_{k}^{f}\left|\mathbf{c}_{k},\mathcal{D}\right. \right)p\left(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1},\mathcal{D }\right.\right)\] \[\times p\left(\mathbf{x}_{k-1},\mathbf{c}_{k-1}\left|\mathbf{z}_{1:k-1 },\mathcal{D}\right.\right) \tag{56}\]

and

\[\int p\left(\mathbf{x}_{k-1},\mathbf{c}_{k-1}\left|\mathbf{z}_{1:k -1},\mathcal{D}\right.\right)d\mathbf{c}_{k-1}\] \[= p\left(\mathbf{x}_{k-1}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right)\] \[= N\left(\mathbf{x}_{k-1};\mathbf{\hat{x}}_{k-1\left|k-1}, \mathbf{P}_{k-1\left|k-1\right.}\right) \tag{57}\]

Substituting \(p(\mathbf{c}_{k}\left|\mathbf{x}_{k-1},\mathbf{c}_{k-1},\mathcal{D}\right.)= N\left(\mathbf{c}_{k};\mathbf{\hat{c}}_{k},\mathbf{P}_{k}^{c}\right)\) and \(p(\Delta_{k}^{f}\left|\mathbf{c}_{k},\mathcal{D}\right.)=N(\Delta_{k}^{f}; \hat{\Delta}_{k}^{f},\mathbf{P}_{k}^{f})\) into Eq. (56), Eq. (55) can be calculated as Eq. (23).

Considering the prediction residual is \(\varepsilon_{k}^{\tau}=\mathbf{x}_{k}-\mathbf{\hat{x}}_{k\left|k-1}\), the state prediction covariance is calculated as

\[\mathbf{P}_{k|k-1}= \mathrm{E}\left[\varepsilon_{k}^{\tau}\varepsilon_{k}^{\tau\top} \left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right]\] \[= \!\!\int\!\!\!\int\!\!\!\int M_{k}^{f}P_{k}^{1}d\Delta_{k}^{f}d \mathbf{c}_{k}d\mathbf{c}_{k-1}d\mathbf{x}_{k-1} \tag{58}\]

with

\[M_{k}^{f}=\left(f_{k}\left(\mathbf{x}_{k-1}\right)+\Delta_{k}^{f}+\mathbf{w}_{ k}-\mathbf{\hat{x}}_{k\left|k-1}\right)\left(\cdot\right)^{\top} \tag{59}\]

Substituting \(N(\mathbf{c}_{k};\mathbf{\hat{c}}_{k},\mathbf{P}_{k}^{c})\) and \(N(\Delta_{k}^{f};\hat{\Delta}_{k}^{f},\mathbf{P}_{k}^{f})\) into Eq. (58) while considering Eq. (57), then we have the covariance prediction as in Eq. (24).

In the state update stage, we compute the state posterior according to the Bayesian rule:

\[p\left(\mathbf{x}_{k}\left|\mathbf{z}_{1:k},\mathcal{D}\right.\right)=\frac{p \left(\mathbf{x}_{k},\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right. \right)}{p\left(\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right)} \tag{60}\]

The joint distribution of state prediction and measurement prediction is Gaussian [35], i.e.,

\[p\left(\mathbf{x}_{k},\mathbf{z}_{k}\left|\mathbf{z}_{1:k-1}, \mathcal{D}\right.\right)\] \[=N\left[\left(\begin{array}{c}\mathbf{\hat{x}}_{k\left|k-1}\\ \mathbf{\hat{z}}_{k\left|k-1}\end{array}\right.\right),\left(\begin{array}{c} \mathbf{P}_{k\left|k-1}&\mathbf{P}_{k\left|k-1}^{xz}\\ \left(\mathbf{P}_{k\left|k-1\right.}^{xz}\right)^{\top}&\mathbf{P}_{k\left|k-1 \right.}^{z}\end{array}\right)\right] \tag{61}\]

Following the conclusions in [35], we substitute Eqs. (18) and (61) into Eq. (60) to obtain updates of the state and covariance as shown in Eqs. (25) and (26). The required \(\mathbf{\hat{z}}_{k\left|k-1\right.}\), \(\mathbf{P}_{k\left|k-1\right.}^{z}\), and \(\mathbf{P}_{k\left|k-1\right.}^{xz}\) are obtained by the following procedure.

The mean value of the Gaussian measurement prediction is calculated as

\[\mathbf{\hat{z}}_{k\left|k-1\right.}= \mathrm{E}\left[h_{k}\left(\mathbf{x}_{k}\right)+\Delta_{k}^{h} \left|\mathbf{z}_{1:k-1}\right.\right]\] \[= \!\!\int\!\!\int\!\!\!\int\left(h_{k}\left(\mathbf{x}_{k}\right)+ \Delta_{k}^{h}\right)\!P_{k}^{2}d\Delta_{k}^{h}d\mathbf{c}_{k}d\mathbf{x}_{k} \tag{62}\]

with

\[P_{k}^{2}=p\left(\Delta_{k}^{h}\left|\mathbf{x}_{k},\mathcal{D}\right.\right)p \left(\mathbf{x}_{k},\mathbf{c}_{k}\left|\mathbf{z}_{1:k-1},\mathcal{D}\right.\right) \tag{63}\]

Substituting \(p\left(\Delta_{k}^{h}\left|\mathbf{x}_{k-1},\mathcal{D}\right.\right)=N( \Delta_{k}^{h};\hat{\Delta}_{k}^{h},\mathbf{P}_{k}^{h})\) into Eq. (63), Eq. (62) is computed as Eq. (27).

Considering the measurement residuals is \(\sigma_{k}^{\tau}=\mathbf{z}_{k}-\mathbf{\hat{z}}_{k\left|k-1}\), the measurement prediction covariance is calculated as

\[\mathbf{P}_{k\left|k-1\right.}^{z}= \mathrm{E}\left[\sigma_{k}^{z}(\sigma_{k}^{z})^{\top}\left| \mathbf{z}_{1:k-1}\right.\right]\] \[= \int\int\int M_{k}^{z}P_{k}^{2}d\Delta_{k}^{h}d\mathbf{c}_{k}d \mathbf{x}_{k} \tag{64}\]

with

\[M_{k}^{z}=\left(h_{k}\left(\mathbf{x}_{k}\right)+\Delta_{k}^{h}+\mathbf{v}_{k}- \mathbf{\hat{z}}_{k\left|k-1}\right.\right)\left(\cdot\right)^{\top}\]