Uncertainty-Aware Variational Inference for Target Tracking

In the low Earth orbit, target tracking with ground based assets in the context of situational awareness is particularly difficult. Because of the nonlinear state propagation between the moments of measurement arrivals, the inevitably accumulated errors will make the target state prediction and the measurement likelihood inaccurate and uncertain. In this article, optimizable models with learned parameters are constructed to model the state and measurement prediction uncertainties. A closed-loop variational iterative framework is proposed to jointly achieve parameter inference and state estimation, which comprises an uncertainty-aware variational filter (UnAVF). The theoretical expression of the evidence lower bound and the maximization of the variational lower bound are derived without the need for the true states, which reflect the awareness and reduction of uncertainties. The evidence lower bound can also evaluate the estimation performance of other Gaussian density filters, not only the UnAVF. Moreover, two rules, estimation consistency and lower bound consistency, are proposed to conduct the initialization of hyperparameters. Finally, the superior performance of UnAVF is demonstrated over an orbit state estimation problem.


I. INTRODUCTION
Aerospace technologies have witnessed a significant development in recent years [1], especially for autonomous vehicles in the low Earth orbit. However, this creates a number of challenges, such as in collision avoidance, and in security due to the growing number of resident space objects [2], [3]. Detection and tracking are essential components of such systems. In response to such challenges, this article focuses on low Earth orbit tracking problems with ground-based assets.

A. Related Work
To cope with the nonlinear estimation problem, a popular classical method is the Gaussian density filter (GDF) [4], which comprises the state prediction stage and the measurement update stage within the Bayesian update framework [5], [6]. These two stages have a number of unknown parameters that need to be updated, but cannot be calculated analytically due to intractable nonlinear integrals. By using different numerical methods to approximate the posterior update parameters of the GDF framework, many nonlinear filters are proposed, such as extended Kalman filter (EKF), unscented Kalman filter [7], cubature Kalman Filter [4], Gaussian-Hermite quadrature filter [8], sparse-grid quadrature filter (SGQF) [9], and their improved variants [10]- [13]. However, due to numerical approximations, there is an inevitable error in the state prediction, which propagates further and impacts the estimation results. Moreover, in order to improve the approximation accuracy, the iterated posterior linearization filter (IPLF) [14] is proposed, which performs iterated statistical linear regression with respect to the posterior instead of the prior so that the approximation accuracy can be improved.
Within the Kalman filtering framework [15], between the moments of measurement arrivals, repeated state predictions are performed until the receipt of the measurements. This could lead to error accumulation in the likelihood, which affects the accuracy of the nonlinear state propagation.
Particle filters [16] have also been popular in solving nonlinear non-Gaussian state estimation problems. However, they face challenges with the state initialization and with high-dimensional spaces [17]. Next, Markov chain Monte Carlo (MCMC) methods [18] are proposed to enhance particle filters where the estimates of the initial state conditions are uncertain [19]. However, as the computational cost of MCMC methods [19] increases with time, conventional MCMC can perform well if the needed convergence conditions are met [18].
Compared with MCMC methods, the variational inference (VI) is a faster and powerful tool [20]. There are the following four kinds of VI, i.e., scalable VI, generic VI, accurate VI, amortized VI [21]. The idea of VI is to first posit a family of densities and then to find the member of that family close to the target, which is measured by Kullback-Leibler (KL) divergence. In VI-based methods [22], [23], the introduction of hyperparameters is inevitable, which need to be artificially set in advance. However, there are no well-established adaptation approaches to determine these hyperparameters. If the hyperparameters deviate from their true values, the estimation performance will be influenced greatly.

B. Contributions
This article has the following main contributions.
1) A state prior model (SPM) and a measurement likelihood model (MLM) are proposed to cope with inaccuracies during the state and measurement predictions. Prior parameters and fitting parameters are introduced in the SPM and MLM, respectively, to reflect the uncertainty in predictions. 2) A closed-loop coordinate ascent variational iterative framework is designed to infer states and introduced parameters by maximizing the variational lower bound (VLBO). It includes model optimization and state estimation, which comprises a novel uncertainty-aware variational filter (UnAVF). 3) Theoretical expressions for the evidence lower bound (ELBO) and the maximization of VLBO are derived, by which the higher estimation accuracy of UnAVF compared with GDF can be theoretically explained. 4) Two rules, estimation and lower bound consistencies, are proposed, which can determine the values of hyperparameters reasonably and adaptively in the VI approach. There is no need to artificially and randomly set hyperparameters in advance.
In the UnAVF, the objective function, (VLBO for parameter inference and state estimation) and the ELBO are known and can be calculated without the need for the true states. However, in conventional Kalman filters, the calculation of the root mean square error (RMSE) needs the true states. Hence, the UnAVF can be adaptively aware and reduce the impact of uncertainties in the state and measurement predictions by monitoring the ELBO and maximizing the VLBO. The proposed UnAVF is validated and tested over an orbit state estimation problem.

C. Organization and Notation
This rest of this article is organized as follows. The considered problem is formulated and the motivation of UnAVF is discussed in Section II. In Section III, the theoretical derivation of the modeling and optimization process in the UnAVF are shown. Section IV proposes a quantitative evaluation indicator for different filters and two rules, estimation consistency and lower bound consistency, to initialize hyperparameters. Section V reports simulation results to demonstrate the superior performance of UnAVF. Finally, Section VI concludes this article.
We will use the following notations. The superscripts "−1" and " " represent the matrix inverse and transpose operations, respectively; | · | and Tr(·) denote the matrix determinant and trace, respectively; N (x|μ, P) denotes that variable x obeys Gaussian distribution with mean μ and covariance P; E [ · ] denotes mathematical expectation; p(x|z) represents the conditional density of x on z; I m denotes the unit matrix; the superscripts "∧" and "∼," used as the hat of random variables, represent the estimate and the estimation error, respectively; for example,x denotes the estimate of variablex and its estimation error isx = x −x; C(·) and D(·) denotes two functions with expressions C(x) = xx and D(x,P) = x Px, respectively.

II. PROBLEM FORMULATION
We consider the nonlinear state-space system given by where x τ ∈ R n and z k ∈ R m denote the system state and the measurement, respectively; τ and k represent the state propagation time and the measurement sampling time, respectively. The sampling intervals of state and measurement are t x = t τ − t τ −1 and t z = t k − t k−1 , respectively. The system noise w τ and the measurement noise v k are zeromean Gaussian white noises with covariances Q and R, respectively. Note that the sampling interval t z of measurements is always larger than the sampling interval t x of states. Hence, states have a long nonlinear propagation during t z without any measurement to correct the accumulated error in mean and uncertainty in covariance of state prediction. When the next measurement arrives, the state and measurement predictions have been influenced by the difference in the sampling rates. This will affect the posterior state estimation processes.

A. Gaussian Density Filter
GDF is a powerful framework for solving the nonlinear state estimation problems. In GDF, the joint state and measurement predictive probability density function (PDF) is assumed to be Gaussian where the state prior PDF and the measurement likelihood PDF are given by where Z 1:k−1 =[z 1 , z 2 , . . . , z k−1 ] . The parameters ξ x , xx , ξ z , zz are calculated via nonlinear integrals in (8)- (11). Based on the assumption in (2), the tractable update formula of the state posterior PDF q GDF (x k ) can be derived naturally [24] as However, the basic GDF framework has some disadvantages. During the long measurement sampling interval t z , GDF can only approximately calculate the state prior PDF p(x k |Z 1:k−1 ) under the minimum mean square error (MMSE) criterion in (8) and (9). The error in mean and uncertainty in covariance of the state prior PDF will accumulate and increase with the nonlinear state propagation. Then, the measurement likelihood PDF will also become inaccurate, because its calculation in (10) and (11) also relies on the corrupted state prior PDF. Finally, the joint Gaussian PDF assumption in (2) cannot reflect the true prior information when the next measurement arrives In spite of this, these inaccurate assumption, state prior, and measurement likelihood PDFs in (8)- (11) are still directly used to update the state posterior PDF in the measurement update stage (5)-(7) without any optimization process. As a result, the estimation performance of GDF will be influenced greatly.

B. Main Idea
In order to accommodate the accumulated error in mean and uncertainty in covariance of the state prediction discussed earlier, we propose an optimizable SPM as where η and are prior parameters.
REMARK 1: In (12), the state prediction PDF p(x k |Z 1:k−1 ) is approximated by the parameterized SPM N (x|η, −1 ). The prior parameters are introduced to reflect the error in mean and uncertainty in the covariance in the state prediction caused by the long nonlinear propagation during t z . Its advantage is that the inaccurate state prediction is only the initial value of the prior parameters at the beginning of the variational iteration. In detail, based on the VI theory, the convergency of the coordinate ascent variational iteration can be guaranteed theoretically. When the variational iteration converges, the optimized SPM will only slightly depend on the inaccurate initial state prediction. Hence, through inferring the prior parameters η and , SPM can accommodate the propagation error in mean and uncertainty in covariance of the state prediction.
Moreover, in order to characterize the uncertainty in the measurement likelihood, an optimizable MLM with fitting parameters is constructed as follows: whereR = R −1 ; H k ∈ R m×n , and u k ∈ R m denote the Jacobian matrix and the first-order constant term in the Taylor expansion of the measurement function, respectively, The fitting parameters μ ∈ R m and the scalar λ need to be optimized.
REMARK 2: Indeed, the MLM (13) essentially fits the latent variable x and the measured data z as a linear parametric Gaussian regression process. Due to the Gaussianity in the MLM, the simple and explicit results of parameter inference and state posterior estimation can be obtained. Moreover, through minimizing the KL divergence to infer the fitting parameters μ and λ, the error in mean and uncertainty in covariance caused by the inaccurate state prediction will be decreased. Specially, the fitting parameter μ is introduced to adjust the mean of p(z k |x k , μ, λ). The fitting parameter λ is designed as a ratio to adjust the covariance. If the approximate error of mean is large, λ will also amplify the covariance to reflect the uncertainty of measurements.
The Gaussian-Wishart and Gaussian-Gamma, distributions [25] are commonly adopted to represent parameters in the VI framework.
where η 0 , β 0 , W 0 , ν 0 , μ 0 , M 0 , c 0 , and d 0 are hyperparameters [25], which denote the prior initialization values for depicting the prior distributions of the introduced parameters. Here, the prior distributions of μ and η are typically Gaussian, which is an usual expression in using VI for regression analysis [26]. The prior parameters and λ are used to manage the covariance accuracy in (12) and (13). It is a standard processing to consider the prior PDFs of and λ as the Wishart and Gamma distributions [27], [28]. Abovementioned constructed forms in (14) and (15) guarantee that the posterior distributions of the fitting and prior parameters conjugate with their prior ones for easing the following maximization of the VLBO. Note that the correlated characteristics can be expressed in other forms, as discussed in [21]. REMARK 3: In the iterative EKF (IEKF) and other Kalman filter types, only the mean and covariance of the state posterior PDF are estimated. There is no optimization of the state model and of the measurement model parameters to reduce the approximation error (such as linearization) of measurement likelihood and state PDFs. Hence, even if the nonlinear functions are approximated via the updated posterior state, the approximation error still exists and is not compensated. In the IPLF, even if there are unknown parameters to be optimized in approximation, only the first moment of their unknown parameters is considered, which cannot capture and reflect the uncertainty and inaccuracy of the state and measurement predictions well, caused by the difference in the sampling rates in this article. Then, in both IPLF and IEKF, based on the approximation which does not contain optimizable parameters or only considers the first moment of parameters, their posterior states will become less informative and more conservative, i.e., their covariances of posterior states will become large. Accordingly, their posterior estimation accuracy will also be influenced. (12) and (13), we distribute the unknown parameters η, , μ, λ, which contain both first and second moments information. Hence, we need to iteratively infer the posterior distributions of the introduced parameters η, , μ, λ and update state posterior estimation, instead of directly calculating the values of unknown parameters (first moment) and posterior states. By inferring the posterior distributions of the introduced parameters, the approximation (12) and (13) in UnAVF can capture and aware the uncertainty and inaccuracy of the state and measurement predictions. Then, based on the approximation considering both the first and second moments of the introduced parameters, the posterior state estimation will become more informative and less conservative, i.e., the covariance of posterior states will become small and estimation accuracy will also be guaranteed. This is also the main reason why our proposed method is called UnAVF. The abovementioned analysis is also demonstrated in the simulations in Figs. 9-10

REMARK 4: In UnAVF in the approximation
The core idea of the proposed UnAVF is to construct an interaction between the model optimization and the state estimation via the maximization of the VLBO, as shown in Fig. 1. Based on the optimized SPM and MLM, the accurate state estimation results are obtained in the VI measurement update. Given the state estimation results, prior and fitting parameters are inferred to optimize the SPM and MLM in the VI state and measurement predictions. The prediction and update comprise the coordinate ascent variational iteration in the proposed UnAVF. In the variational iteration, the inaccurate and uncertain state and measurement information only provide initial values for the optimization of the SPM and MLM. As the convergency of the variational iteration, the finial iterative estimation results will not be influenced largely by the uncertain initial information. In  (14) and (15) need to be initialized at the beginning of the variational iteration. Hence, is there a rule to rationally conduct the adaptive initialization, instead of setting their values artificially and casually?

III. UNCERTAINTY-AWARE VARIATIONAL FILTER
In Section II, the motivation and core idea of the Un-AVF have been discussed. Through answering the first question, the complete UnAVF algorithm will be proposed in Section III. Based on the VI framework, the VI state and measurement predictions and the VI measurement update will be derived to iteratively maximize the VLBO for jointly optimizing the SPM and MLM and estimating states.
To estimate the state x k and inferring the prior parameters η, and the fitting parameters μ, λ, the joint posterior PDF p(x k , η, , μ, λ|Z 1:k ) needs to be calculated. Because there is no analytical solution to the joint posterior PDF in the nonlinear system (1), the VI approach is, therefore, employed to obtain a suboptimal approximation for the joint posterior PDF. Based on the VI approach, we are going to look for an approximate solution by making the following variational approximation: where = {x k , η, , μ, λ}; q(.) denotes the variational posterior PDF of p(.). The variational posterior PDFs q(x k ), q(η, ), q(μ, λ) can be calculated by minimizing the following KL divergence between the factorized variational posterior PDFs q(x k )q(η, )q(μ, λ) and the true joint posterior PDF p( |Z 1:k ) [21], [25], [26] as p(x) dx is the KL divergence between q(x) and p(x). Based on the VI theory [25], the minimization of the KL divergence (18) is equal to the maximization of the VLBO as follows: The optimal solution to (19) satisfies the equation [25] as where θ is an arbitrary element of and ( =θ) is a subset of with θ ∪ ( =θ) = . The operator E ( =θ) [· · · ] denotes an expectation with respect to the variational posterior PDF q( ( =θ) ). Because the calculations of q(x k )q(η, )q(μ, λ) are coupled with each other, the coordinate ascent variational iteration is needed to solve (21). In detail, t denotes the number of variational iteration. Based on the results of the tth variational iteration q t ( ( =θ) ), the variational posterior PDF q(θ) of an arbitrary element θ is calculated as q t+1 (θ) at the t + 1th variational iteration by solving the expectation in (21).
In the following, we will calculate the variational posterior PDFs q(η, )q(μ, λ)q(x k ) at the t + 1th variational iteration, denoted as q t+1 (η, )q t+1 (μ, λ)q t+1 (x k ) given by Theorems 1-3, respectively. To this end, using the conditional independence properties of the Gaussian Gamma state-space model in (1) and (14) and (15), the joint complete-data likelihood PDF can be factored as The direct probability graph is given in Fig. 2 to illustrate the factorized complete-data likelihood PDF. The state x k is controlled by the prior parameters η, , and the measurement z k is controlled by the fitting parameters μ, λ and the state x k . There is no direct connection between the fitting parameters and the prior parameters so they are independent with each other. Now, given the factorized complete-data likelihood PDF in (22), based on the solution (21), the VLBO can be maximized gradually by our derived coordinate ascent variational iteration process. At the same time, the model optimization and state estimation can also be iteratively achieved in Theorems 1-3.
Theorem 1 (VI state prediction): Let θ = {η, } and accordingly the variational posterior PDF of ( =θ) at the previous tth variational iteration is then, based on (21) and (22), the variational posterior PDF of the prior parameters at the t + 1th iteration are calculated as PROOF: See Appendix A.
Theorem 3 (VI Measurement Update): Let θ = {x k } and accordingly the variational posterior PDF of ( =θ) at the t + 1th variational iteration is then, based on (21) and (22), the variational posterior PDF of states at the t + 1th variational iteration is calculated as PROOF: See Appendix C.
REMARK 5: In Theorems 1 and 2, the optimization of the SPM and MLM is achieved. Then, based on the optimized models, the state is estimated in Theorem 3. Accordingly, the state posterior estimation will also contribute to the optimization in Theorems 1 and 2. The VI state and measurement prediction and the VI measurement update in Theorems 1-3 comprise the variational iteration by maximizing the VLBO in the UnAVF. The uncertain and inaccurate state prior and measurement likelihood PDFs are only the initial value of the variational iteration. As the increase of the VLBO in the variational iteration, the prior and fitting parameters can be inferred so that the impact of uncertainties on the state estimation will be reduced gradually.
REMARK 6: In our proposed UnAVF, the optimization of the MLM, i.e., Theorem 2 (VI measurement prediction), can be considered as relinearization. Specifically, in Theorem 2, the distributions of the learned parameters μ, λ in the MLM (13) is updated. Besides, in each iteration, statesx k is also updated so H k will also be recalculated using the  The terminal condition of the variational iteration at each sampling time is to measure the difference between the tth and t + 1th variational iteration results. If the difference is too small, we can assess that the coordinate ascent variational iteration has converged and should be stopped. Hence, in the proposed UnAVF, the KL divergence with respect to the state posterior PDFs in the tth and t + 1th variational iteration is considered as the terminal index. The setting of the threshold value δ will be given in the simulation part.

IV. PERFORMANCE ANALYSIS
As discussed in Section III, the hyperparameters in the UnAVF need to be initialized at beginning of variational iteration. In Section IV, two rules will be proposed to conduct the adaption initialization of hyperparameters. Moreover, we need to derive an index to evaluate the KL divergences of different filters for quantitatively comparing their accuracy performance. To this end, it is necessary to transform the KL divergence evaluation into the corresponding lower bound evaluation.

A. Quantitative Assessment by the ELBO
The VLBO in (20) is the objective function of the UnAVF with respect to states and parameters η, , μ, λ.
Through iteratively maximizing the VLBO, states and parameters can be joint inferred. This means that the increase of the VLBO is due to the mutual effect of states and parameters. Consequently, the VLBO cannot measure the accuracy of other filters, which do not introduce these parameters. Hence, we need to derive a new lower bound only with respect to states to fairly and theoretically evaluate the performance of different filters. To this end, the ELBO is derived from the following equation: where The subscript "E" aims to distinguish with the VLBO in (20) and q(x k ) denotes the variational posterior PDF of states at the tth variational iteration.
REMARK 7: From (44), q(x k ) will approximate the true posterior PDF p(x k |Z 1:k ) if the ELBO in (46) increases. Thus, the ELBO has the ability to quantitatively evaluate the accuracy performance of different filters, i.e., that the higher the filter's ELBO is, the better the estimation performance is.
Calculating the ELBOs of the GDF and UnAVF starts from According to the modeling process in the GDF and Un-AVF, the expression forms of the measurement likelihood PDF, the state prior PDF and the state posterior PDF can be summarized, respectively, as The difference of (48)-(50) in the GDF and UnAVF depends on the expressions of these parameters A k , B k , C k ,x k/k−1 , P k/k−1 ,x k/k ,P k/k . Given (48)-(50), the completely uniform expression of the ELBO in both GDF and UnAVF is where m and n represent the dimensions of measurement and state vectors, respectively.

1) The ELBO of GDF
In GDF, expressions of parameters A k , B k , C k ,x k/k−1 , P k/k−1 ,x k/k ,P k/k are given by Then, we can obtain the following complete expression of the ELBO in GDF: 2

) The ELBO of the UnAVF Based on the modeling of the SPM and MLM in Section II, expressions of parameters
Then, the complete expression of the ELBO in the UnAVF can be calculated as Summarizing the above, computing ELBOs of different filters is achievable and tractable, without nonlinear integrals. Moreover, it provides a quantitative assessment to the accuracy performance of different filters. In (44), the ELBO increase reflects the decrease of the KL divergence so that the variational approximate distribution q(x k ) approximates well the true posterior p(x k |Z k ). If the ELBO of the UnAVF is higher than that of GDF, then the UnAVF can outperform GDF.

B. Initialization of Hyperparameters at the Beginning of the Variational Iteration
The principle of initialization is to ensure that the Un-AVF can at least outperform GDF. To this end, two rules, estimation consistency and lower bound consistency, are proposed to conduct the initialization of hyperparameters. 1) Estimation Consistency: the state estimation of the UnAVF at the beginning of variational iteration t = 0 at least is identical with that of GDF. According to Theorem 3, only if the hyperparameters μ 0 , c 0 , d 0 , η 0 , W 0 , ν 0 are initialized as the state posterior estimation of the UnAVF q t=1 (x k ) will be identical with that of GDF q GDF (x k ) in (5). The similar proof can be seen in [29] in Appendix B. In other words, GDF only provides a prior estimation to the proposed UnAVF. Then, due to the variational iteration, the UnAVF's performance will become better than that of GDF.
2) Lower Bound Consistency: Based on the estimation consistency, the ELBO of the UnAVF at the beginning of the variational iteration should also be identical with that of GDF. The aim of the lower bound consistency is to initialize the hyperparameters and to further theoretically explain why the UnAVF can outperform GDF At the beginning of the variational iteration t = 0, the ELBOs of the UnAVF and GDF are calculated in (57) and (58), respectively. Based on the estimation consistency, the third terms in (57) and (58) are the same. To achieve the lower bound consistency, the key point is to guarantee that the first and second terms in (57) equal to that in (58), which are calculated in (61), (62) and (63), (64), respectively. As a result, two conditions are yielded as Summarizing abovementioned derivation and analysis, we can obtain the following proposition for initializing hyperparameters at the beginning of the variational iteration.
Proposition 1: Based on both estimation consistency and lower bound consistency, the hyperparameters in the UnAVF should be initialized as follows: where ξ x and xx are calculated in the same way as GDF in (8) and (9). Although there are many hyperparameters, based on Proposition 1, only c 0 and ν 0 need to be initialized. Frankly speaking, there may exist other initialization methods, but at least Proposition 1 provides an available initialization scheme with both operability and practicability. REMARK 8: Note that with the condition of the two rules in Proposition 1 to determine the hyperparameters, at the beginning of the variational iteration t = 0, the ELBO of UnAVF is the same as that of GDF. As shown in Fig. 3, then, with the variational iteration proceeding, the iterative optimization of the SPM and MLM can further promote the ELBO of the UnAVF. Finally, the ELBO promotion in the UnAVF implies its accuracy improvement and theoretically Fig. 3. Increase of the ELBO explains why the posterior estimation of the UnAVF is more accurate than that of GDF.

V. SIMULATION
In Section V, the proposed UnAVF is compared with the EKF, SGQF, IEKF, and IPLF. The performance of the UnAVF is demonstrated in the orbit estimation problem.
In order to evaluate the performance of different filters, we run different filters with 1000 Monte Carlo simulations and use the RMSE where N mc = 1000 is the number of Monte Carlo simulations; K is the simulation length in time steps; x n k [i] is the true value of state andx n k [i] is the estimated value of state in nth simulation at time k. Furthermore, we used the average RMSE (ARMSE) with respect to times k 1 s to k 2 s We also considered the mean absolute error (MAE) over time k of the nth simulation run The dynamic model of the low Earth orbit satellite [9] is given byr where r = [x, y, z] is the position of satellite in the inertial coordinate frame (I-J-K); scalar r = x 2 + y 2 + z 2 . Vector ν is the zero-mean Gaussian state noise [9]. Vector a G is the acceleration caused by the J 2 perturbation [30]. Vector a D is the atmospheric drag [31]. The measurement model is described by where the azimuth (az), the elevation (el), and the range ρ = [ρ u , ρ e , ρ n ] can be obtained by the radar site on the ground with respect to the local observer coordinate system. n az , n el , and n ρ are the white Gaussian noise. The transformation relationship between range ρ of measurement and position r of state in the inertial coordinate frame is described by where L = 6378.1363 km is the Earth radius; ε and ϑ are the latitude and local sidereal time of the observer, respectively. Generally speaking, a single radar cannot track a satellite with the entire orbit. In this simulation, the rational tracking time is 300 s. The measurement update interval t z is 5 s. For the accurately describing the state propagation, the dynamic model of the low Earth orbit satellite is discretized by fourth-order Runge-Kutta algorithm with the step size t x = 0.1 s. Hence, the satellite states have long nonlinear propagation without any measurement. Other information about the radar is identical with that in [9]. The trajectory of the low Earth orbit satellite is shown in Fig. 4.
In the simulation of the low Earth orbit satellite, the reasonable measurement noise covariance is assumed to be [9] R k = diag (0.015 • ) 2 (0.015 • ) 2 0.025 2 (km) 2 . (72) The six dimensional state is  .841182]km/s In each Monte Carlo simulation, the initial states of different filters are both generated randomly from the Gaussian distribution N (x 0|0 ,P 0|0 ), wherê P p 0|0 = 10 4 , 10 4 , 10 4 km 2 P v 0|0 = 10 −2 , 10 −2 , 10 −2 (km/s) 2 . (77) According to the Proposition 1 in Section IV, only hyperparameters c 0 and ν 0 in the UnAVF need to be initialized as In the following, we will evaluate and compare performance of different filters from estimation error and estimation credibility aspects. Moreover, we will demonstrate the validity of the variational iteration in the UnAVF.

A. Estimation Accuracy
First, the ARMSEs of different filters with respect to different times are shown in Tables I and II. The ARMSE of the UnAVF is always the smallest from 1 to 300 s. From times 1 to 100 s, the ARMSEs of different filters are similar. However, when the results of different filters converge (from 201 to 300 s), the accuracy of UnAVF is increased by 51.52 and 40% in the position and velocity, respectively, compared with the IEKF.
Then, the RMSEs of different filters and the posterior Cramer-Rao lower bound (PCRLB) [32] are shown in Fig. 5. At each sampling time of measurements, the UnAVF is better than other filters in the RMSEs of the position and  Tables I and II with different times. Note that we do not show the EKF's RMSE curves, because its RMSE is quite large compare with other filters. Hence, we only report its ARMSE in Tables I and II. As we discussed in Section IV, only two hyperparameters ν 0 , c 0 need to be initialized at the beginning of the variational iteration at each sampling time.  Tables I and II and Fig. 5. Hence, for clearly reporting the difference of RMSE curves, we show most of the curves. Moreover, there is a trend that the larger c 0 and ν 0 is, the lower the estimation error of the UnAVF is.
The trend of the results matches the theoretical analysis. According to Proposition 1, we always set d 0 = c 0 and η 0 = ξ x , which means that E (λ) = c 0 (d 0 ) −1 = 1 and E (η) = ξ x [ξ x is the GDF's state prediction mean, which is calculated by (8)]. However, we do not know any information about  λ and ξ x is also inaccurate. A direct and efficient method to make the UnAVF understand that E (η) and E (λ) are not reliable is to increase the variances of λ and η. The larger variance means the lower credibility and validity of mean. Hence, the parameters c 0 and ν 0 should be small so that the variances of λ and η will be large where V(·) denotes a variance. This can clearly explain why smaller c 0 and ν 0 are more rational and can yield more accurate estimation results.

B. Estimation Credibility
Besides the estimation error, the estimation credibility is significant as well. For evaluating different filters' estimation credibility, the contour of the state prior and posterior PDFs at times 1 s, 150 s, and 300 s are shown in Figs. 8-10. At time 1 s, the initial state prior PDFs of all filters are almost the same and inaccurate. However, at times 150 s and 300 s, the state prior PDF's mean in the UnAVF is obviously more accurate and its covariance is also smaller. This is because the flexible SPM in the UnAVF can be iteratively optimized to fit the true state situation by inferring prior parameters using Theorem 1, which outperforms the unadjustable state prior PDF in GDF.  Moreover, based on a more accurate state prior PDF and the approximation considering both the first and second moments of introduced parameters, the calculation of the state posterior PDF will be more credible and informative, which means a smaller covariance. In Figs. 9 and 10, the posterior PDF's contours of the UnAVF are more tight and the contours' centers of the UnAVF are more close to the true state value as well. These simulation results can further demonstrate that the UnAVF can handle the inaccurate and uncertain state and measurement predictions.

C. Validity of Variational Iteration in the UnAVF
In the proposed UnAVF, the variational iteration consists of the VI state and measurement predictions and the VI measurement update. The state estimation in the VI measurement update has been reported in above simulation results. Hence, to demonstrate the validity of the variational iteration, we will focus on the optimization of the SPM and MLM in the VI state and measurement predictions, respectively, and the maximization of the ELBO.  As discussed in Section IV, the ELBO is an evaluation about the accuracy performance of different algorithms. The higher the ELBO is, the better the estimation performance is. In Fig. 11, ELBOs of different filters are shown. The final ELBO of the UnAVF is obvious higher than that of other filters. The higher ELBO of the UnAVF can further explain why the state posterior PDF of the UnAVF is closer to the true one than that of other filters. Moreover, the increase from the initial to final ELBOs can also demonstrate the validity and contribution of the variational iteration at each sampling times.
The KL divergence of the estimated and true measurement likelihood PDFs at each sampling time are shown in Fig. 12. The final KL divergence of the UnAVF is smaller than that of other filters. Thus, the optimized MLM in the UnAVF can fit the true measurement model better. The decrease from the initial to final KL divergence further demonstrates the validity of the variational iteration in the UnAVF.
In addition, for more directly exhibiting the variation of the MLM and SPM in the variational iteration, the changes of the measurement likelihood PDF and the state prior error are shown in Figs. 13 and 14, respectively. As the variational iteration proceeding, the optimized measurement likelihood PDF can gradually fit the true one and the state prior error is decreased step by step. Hence, these changes of the ELBO, the KL divergence and the state prior error can clarify the impact and validity of the variational iteration in the UnAVF.

VI. CONCLUSION
Based on the VI framework, this article develops a novel nonlinear estimation method, which dynamically optimizes the parameterized state prior and MLMs by maximizing the VLBO. In the VI state and measurement predictions, the prior parameters and fitting parameters are inferred so that the SPM and MLM will be self-adaptively adjusted. Correspondingly, based on the optimized SPM and MLM, the state posterior PDF can be calculated more accurately in the VI measurement update. The uncertainties and inaccuracies in the state priori and measurement likelihood PDFs can only have effect at the beginning of the variational iteration. As the maximization of the VLBO by the variational iteration proceeding, the uncertainties' effect will be decreased gradually. Moreover, the estimation and lower bound consistency are proposed, which can rationally guide the initialization of hyperparameters at the beginning of the variational iteration at each sampling time. In the simulation, we have shown that the performance of the UnAVF is better and the validity of the variational iteration is demonstrated.