CN112528479A - Robust self-adaptive smoothing method based on Gibbs sampler - Google Patents

Abstract

The invention belongs to the technical field of state estimation, and particularly relates to a robust self-adaptive smoothing method. A robust self-adaptive smoothing method based on a Gibbs sampler firstly models process noise and observation noise in a linear state space model into Student's t distribution, and noise parameters of the distribution are completely unknown. The process noise and observation noise satisfying the Student's t distribution are decomposed into a combination of Gauss and Gamma distributions by introducing auxiliary parameters. The unknown noise parameters are also regarded as random variables, and the prior distribution of the random variables is modeled as conjugate prior corresponding to the parameters. And under the framework of a Gibbs sampler, carrying out iterative sampling on unknown noise parameters, auxiliary parameters and the system state at the same time. After multiple iterations are performed, the average value of the iteration period state samples after the steady state is reached is selected as the final state estimation value. The robust adaptive smoothing method provided by the invention can still obtain a better state estimation result when the model noise is thick tail and the initial noise parameter setting error is larger.

Description

Robust self-adaptive smoothing method based on Gibbs sampler

Technical Field

The invention belongs to the technical field of state estimation, and particularly relates to a robust self-adaptive smoothing method.

Background

The Rauch-Tung-striebel (rts) smoother, also known as Kalman smoother, is widely used in the field of linear state estimation. But the Kalman smoother can only guarantee its optimality when the process noise and the observation noise of the state space model are gaussian distributed and the noise parameters are completely known. However, in practical estimator applications, model noise parameters are often difficult to obtain accurately, and the sensor often has observation outliers. For example, in an underwater positioning system based on ranging, uncertainty parameters of underwater acoustic ranging are difficult to obtain accurately, and are affected by an underwater harsh communication environment, a ranging outlier often occurs, and the performance of the Kalman smoother is seriously deteriorated, and even the estimator may be diverged.

Some existing adaptive estimators can solve the problem of unknown noise parameters to some extent. The main adaptive estimation methods at present include a maximum likelihood method, a correlation method, a variance matching method and a Bayesian method. But these methods are all based on the gaussian distribution of model noise. The appearance of the outlier can make the original gaussian noise exhibit thick tail characteristics, that is, the probability density of the noise value in the region far from the mean value becomes large. Student's t distribution can better model thick tail noise, and therefore, is widely applied to the problem of state estimation with outliers. The current robust estimator based on Student's t distribution design mainly adopts the Variational Bayes (VB) approximation technology. A correlated VB approximation based robust adaptive estimator can simultaneously handle the thick tailness of noise and the unknown of noise parameters. However, the VB-approximation-based robust adaptive estimation method usually needs to perform a free factorization approximation of the joint probability density, and only an approximate solution of the target probability density distribution is obtained, so that the estimation accuracy is usually difficult to guarantee. The Markov Chain Monte Carlo (MCMC) method is a stochastic approximation method and is widely used in bayesian inference. MCMC can sample complex distributions by generating a Markov chain that is quasi-static with respect to the target distribution. When the number of samples is sufficiently large, the MCMC method can ensure that sufficiently accurate estimates are produced. The MCMC method has the potential to obtain better estimation accuracy than a method based on VB approximation when being used in the field of robust adaptive estimation, and no relevant research is found at present. Because of its simplicity and ease of use, the Gibbs sampler is the most widely used method in the MCMC framework.

Disclosure of Invention

The purpose of the invention is: aiming at the problems that the observation outlier of a sensor and the noise parameter of a model are unknown frequently in the actual linear state estimation application, the model noise is modeled into the Student's t distribution with unknown parameters, and a robust self-adaptive smoothing method based on a Gibbs sampler is designed, so that the thick tailability of the model noise and the unknown noise parameter can be processed simultaneously.

The technical scheme of the invention is as follows: the process noise and observation noise in the linear state space model are first modeled as Student's t distributions, and their noise parameters are completely unknown. The process noise and observation noise satisfying the Student's t distribution are decomposed into a combination of Gauss and Gamma distributions by introducing auxiliary parameters. The unknown noise parameters are also regarded as random variables, and the prior distribution of the random variables is modeled as conjugate prior corresponding to the parameters. And under the framework of a Gibbs sampler, carrying out iterative sampling on unknown noise parameters, auxiliary parameters and the system state at the same time. The flow of each iteration cycle is as follows: 1) calculating the state posterior distribution of the current iteration cycle based on a Kalman smoother under the condition of the noise parameter and the auxiliary parameter sampled in the previous iteration cycle, and sampling the system state of the current iteration cycle from the posterior distribution; 2) calculating posterior distribution of an observation noise mean value and a scale matrix under the condition of a system state, an observation variable and an observation noise auxiliary parameter sampled in the current iteration period, and sampling the observation noise mean value and the scale matrix in the current iteration period from the posterior distribution; 3) calculating posterior distribution of a process noise mean value and a scale matrix by taking a system state and a process noise auxiliary parameter sampled in a current iteration period as conditions, and sampling the process noise mean value and the scale matrix of the current iteration period from the posterior distribution; 4) calculating posterior distribution of the observation noise auxiliary parameters under the condition of the system state, the observation variable and the observation noise parameters sampled in the current iteration period, and sampling the observation noise auxiliary parameters in the current iteration period from the posterior distribution; 5) calculating posterior distribution of the process noise auxiliary parameters by taking the system state and the process noise parameters sampled in the current iteration period as conditions, and sampling the process noise auxiliary parameters in the current iteration period from the posterior distribution; 6) calculating posterior distribution of the observation noise freedom degree under the condition of observation noise auxiliary parameters sampled in the current iteration period, and sampling the observation noise freedom degree of the current iteration period from the posterior distribution; 7) and calculating posterior distribution of the process noise freedom degree by taking the process noise auxiliary parameters sampled in the current iteration period as conditions, and sampling the process noise freedom degree of the current iteration period from the posterior distribution. After multiple iterations are performed, the average value of the iteration period state samples after the steady state is reached is selected as the final state estimation value.

The invention comprises the following steps:

A. modeling the process noise and the observation noise of the linear state space model into Student's t distribution, wherein the parameters of the process noise and the observation noise are random variables, and the prior distribution of the parameters of the process noise and the observation noise is the corresponding conjugate prior distribution.

B. And introducing auxiliary parameters, and decomposing the process noise and the observation noise which meet the Student's t distribution into a combination of Gauss distribution and Gamma distribution by using a heuristic Gauss model.

C. And initially sampling process noise, observation noise parameters and auxiliary parameters of the linear state space model under the frame of a Gibbs sampler.

D. Setting a total number of iteration cycles N and a number of stationary cycles N_bIs required to satisfy N_b< N; and performing N iterations, wherein the flow of each iteration is as follows:

D1. and calculating the state posterior distribution of the current iteration period based on a Kalman smoother under the condition of the noise parameter and the auxiliary parameter sampled in the previous iteration period, and sampling the system state of the current iteration period from the posterior distribution.

D2. And calculating posterior distribution of the observation noise mean value and the scale matrix under the condition of the system state, the observation variable and the observation noise auxiliary parameter sampled in the current iteration period, and sampling the observation noise mean value and the scale matrix in the current iteration period from the posterior distribution.

D3. And calculating posterior distribution of the process noise mean value and the scale matrix by taking the system state and the process noise auxiliary parameters sampled in the current iteration period as conditions, and sampling the process noise mean value and the scale matrix in the current iteration period from the posterior distribution.

D4. And calculating posterior distribution of the observation noise auxiliary parameters under the condition of the system state, the observation variable and the observation noise parameters sampled in the current iteration period, and sampling the observation noise auxiliary parameters in the current iteration period from the posterior distribution.

D5. And calculating posterior distribution of the process noise auxiliary parameters by taking the system state and the process noise parameters sampled in the current iteration period as conditions, and sampling the process noise auxiliary parameters in the current iteration period from the posterior distribution.

D6. And calculating posterior distribution of the observation noise freedom degree under the condition of the observation noise auxiliary parameters sampled in the current iteration period, and sampling the observation noise freedom degree of the current iteration period from the posterior distribution.

D7. And calculating posterior distribution of the process noise freedom degree by taking the process noise auxiliary parameters sampled in the current iteration period as conditions, and sampling the process noise freedom degree of the current iteration period from the posterior distribution.

E. After N iterations, selecting N-N after reaching steady state_bAnd taking the average value of the state samples in the secondary iteration period as a final state estimation value.

On the basis of the above scheme, specifically, the method adopted in step a is as follows:

for a linear state space model:

x_k＝F_k-1x_k-1+w_k-1

z_k＝H_kx_k+v_k

wherein:

is a vector of the states of the system,

in order to observe the vector for the system,

in order to be a state transition matrix,

in order for the system to observe the matrix,

in order to be a noise of the process,

to observe noise; w is a_kAnd v_kBoth modeled as Student's t distributions, as follows:

p(w_k)＝St(w_k；ξ，Q，ω)

p(v_k)＝St(v_k；μ，R，v)

wherein: st (x; mu, sigma, omega) represents a random vector x which satisfies Student's t distribution and takes mu as a mean vector, sigma as a scale matrix and omega as a degree of freedom; xi, Q, omega, mu, R and v are random variables, and the prior distribution of the random variables is modeled as corresponding conjugate prior; wherein, the prior distribution of mu and R is modeled as Gaussian-Inverse Wishart (GIW) distribution:

wherein: GIW (X, X; a, B, c, D) represents a Gaussian-inverse Wishart distribution with a, B, c, D as parameters; IW (X; a, B) represents that the random matrix X meets inverse Wisharp distribution with a as the degree of freedom and B as the inverse scale matrix; n (x; a, A) represents a Gaussian distribution in which the random vector x satisfies a as a mean and A as a variance; likewise, the prior distribution of ξ and Q is also modeled as a GIW distribution:

the degree of freedom parameters v and ω are both modeled as a Gamma distribution, as follows:

p(v)＝Gamma(v；e₀，f₀)

p(ω)＝Gamma(ω；g₀，h₀)

wherein: gamma (x; a, b) represents a random variable x satisfying the Gamma distribution with a as a shape parameter and b as a rate parameter; setting nominal values of mean vectors of process noise and observation noise to be xi respectively₀And mu₀(ii) a Order:

E(ξ)＝ξ₀

E(μ)＝μ₀

wherein: e (-) is the expectation of a random variable; to obtain:

α₀and beta₀Is a modulation parameter; setting nominal values of process noise and observation noise scale matrix to be Q respectively₀And R₀(ii) a Order:

E(Q)＝Q₀

E(R)＝R₀

to obtain:

t₀＝ρ_t+2n+2

T₀＝ρ_tQ₀

u₀＝ρ_u+2m+2

U₀＝ρ_uR₀

wherein: rho_uAnd rho_tIs a modulation parameter;

setting nominal values of process noise and observation noise freedom as omega₀And v₀(ii) a Order:

E(ω)＝ω₀

E(v)＝v₀

to obtain:

on the basis of the above scheme, specifically, the method adopted in step B is:

introducing an auxiliary parameter gamma_kAnd lambda_kThe process noise and the observation noise are decomposed into:

on the basis of the above scheme, specifically, the method adopted in step C is:

according to the process noise scale matrix and the observation noise scale matrix prior distribution parameter t obtained in the step A₀，T₀，u₀，U₀Directly from IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Obtaining initial sampling of a process noise scale matrix and an observation noise scale matrix by intermediate sampling; namely:

Q⁽¹⁾～IW(Q；t₀，T₀)

R⁽¹⁾～IW(R；u₀，U₀)

wherein: a. the^(j)Representing the sampling value of the random variable A in the jth iteration;

the prior distribution parameter alpha of the process noise mean vector and the observation noise mean vector obtained according to the step A₀，β₀，

And slave IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Q of middle sampling⁽¹⁾And R⁽¹⁾Directly from

And

obtaining a process noise mean vector and an observation noise mean vector initial sample by intermediate sampling; namely:

according to the process noise degree of freedom and observation noise degree of freedom prior distribution parameter e obtained in the step A₀，f₀，g₀，h₀Directly from Gamma (omega; g)₀，h₀) And Gamma (v; e.g. of the type₀，f₀) Obtaining initial sampling of process noise freedom and observation noise freedom through intermediate sampling; namely:

ω⁽¹⁾～Gamma(ω；g₀，h₀)

v⁽¹⁾～Gamma(v；e₀，f₀)

according to the process noise auxiliary parameter and observation noise auxiliary parameter prior distribution parameter obtained in the step B, directly obtaining the noise auxiliary parameter

And

obtaining initial sampling of the process noise auxiliary parameter and the observation noise auxiliary parameter by intermediate sampling; namely:

on the basis of the above scheme, specifically, the method adopted in the step D1 is as follows:

noise mean vector sample value xi during previous iteration cycle⁽ⁱ⁾Scale matrix sample value Q⁽ⁱ⁾And auxiliary parameter sampling value gamma_k ⁽ⁱ⁾Observation of noise mean vector sample value mu⁽ⁱ⁾Scale matrix sample value R⁽ⁱ⁾And an auxiliary parameter lambda_k ⁽ⁱ⁾Under the known condition, calculating a state posterior distribution parameter by adopting a Kalman smoother;

1) initializing system state and variance:

wherein: the subscript i | j represents the estimate of the variable at time i conditioned on the system observations prior to time j;

2) forward recursion:

let T be the total number of times, for k 1, 2, 3.

a) Prediction

Calculating the prior state and the prior variance at the k moment according to the posterior state and the posterior variance at the k-1 moment as follows:

P_k|k-1＝F_k-1P_k-1|k-1F_k-1 ^T+Q/γ_k

b) updating

Calculating the posterior state and posterior variance at the moment k as follows:

K_k＝P_k|k-1H_k ^T(H_kP_k|k-1H_k ^T+R/λ_k)

P_k|k＝P_k|k-1-K_kH_kP_k|k-1

wherein: k_kIs Kalman gain;

3) recursion backwards

For k ═ T-1, T-2.. 0, the following steps are performed:

P_k|T＝P_k|k+G_k(P_k+1|T-P_k+1|k)G_k ^T

wherein: g_kIs the smoothing gain;

according to the above Kalman smoother, xi is obtained⁽ⁱ⁾，Q⁽ⁱ⁾，γ_1：T ⁽ⁱ⁾，μ⁽ⁱ⁾，R⁽ⁱ⁾And lambda_1：T ⁽ⁱ⁾A state posterior parameter of a condition

And

further, for k equal to 0, 1, 2, 3.. T, the following are performed in order:

sampling is carried out;

the method adopted in the step D2 is as follows:

according to the Bayes rule:

p(μ，R|z_1：T，x_0：T，λ_1：T)＝cp(μ，R)p(z_1：T|μ，R，x_0：T，λ_1：T)

its logarithmic form is:

log p(μ，R|z_1：T，x_0：T，λ_1：T)＝log p(μ，R)+log p(z_1：T|μ，R，x_0：T，λ_1：T)+log c

as shown in the step a, the prior distribution of the observation noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the gaussian-inverse Wishart distribution is as follows:

the log form of the likelihood distribution can be written according to the observation model:

and calculating to obtain:

wherein:

the posterior distribution of the mean vector of the observed noise is sufficient

Is an average value of

The posterior probability density of the scale matrix is based on the Gaussian distribution of the variance

For the degree of freedom, to

The inverse Wishart distribution of the inverse scale matrix is obtained; in the actual sampling process, the following steps are required:

sample R⁽ⁱ⁺¹⁾And further according to:

sampling mu⁽ⁱ⁺¹⁾Obtaining an observation noise mean vector and a scale matrix sampling value of the current iteration cycle;

the method adopted in the step D3 is as follows:

according to the Bayes rule:

p(ξ，Q|x_0：T)＝cp(ξ，Q)p(x_0：T|ξ，Q)

its logarithmic form is:

log p(ξ，Q|x_0：T)＝log p(ξ，Q)+log p(x_0：T|ξ，Q)+log c

as shown in the step a, the prior distribution of the process noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the prior distribution is as follows:

the log form of the likelihood distribution is written according to a kinematic model:

and calculating to obtain:

wherein:

process noise mean satisfied with

Is an average value of

The posterior probability density of the scale matrix is based on the Gaussian distribution of the variance

For the degree of freedom, to

The inverse Wishart distribution of the inverse scale matrix is obtained; in the actual sampling process, the following steps are required:

sample Q⁽ⁱ⁺¹⁾And further according to:

sample xi⁽ⁱ⁺¹⁾Obtaining a process noise mean vector and a scale matrix sampling value of the current iteration cycle;

the method adopted in the step D4 is as follows:

according to the Bayes rule:

p(λ_k|x_k，z_k，μ，R，v)＝cp(λ_k|v)p(z_k|x_k，λ_k，μ，R)

its logarithmic form is:

log p(λ_k|x_k，z_k，μ，R，v)＝log p(λ_k|v)+log p(z_k|x_k，λ_k，μ，R)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

the posterior distribution of the noise auxiliary parameter is observed sufficiently

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling lambda_k ⁽ⁱ⁺¹⁾Obtaining an observation noise auxiliary parameter sampling value of the current iteration period, wherein k is 1, 2,. T;

the method adopted in the step D5 is as follows:

according to the Bayes rule:

p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝cp(γ_k-1|ω)p(x_k|x_k-1，γ_k-1，ξ，Q)

its logarithmic form is:

log p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝log p(γ_k-1|ω)+log p(x_k|x_k-1，γ_k-1，ξ，Q)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

process noise auxiliary parameter posterior distribution satisfied by

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling gamma_k-1 ⁽ⁱ⁺¹⁾，k＝1，Obtaining a process noise auxiliary parameter sampling value of the current iteration period;

the method adopted in the step D6 is as follows:

according to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(v)＝(e₀-1)log v-f₀v+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

observing the posterior distribution of the degree of freedom of the noise to satisfy

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sample v⁽ⁱ⁺¹⁾Obtaining an observation noise freedom degree sampling value of the current iteration period;

the method adopted in the step D7 is as follows:

according to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(ω)＝(g₀-1)log ω-h₀ω+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

process noise degree of freedom posterior distribution satisfied by

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling omega⁽ⁱ⁺¹⁾And obtaining a process noise freedom degree sampling value of the current iteration period.

On the basis of the above scheme, specifically, the method adopted in step E is:

after N times of iterative sampling, a state variable sampling set is obtained as

Selecting N-N after reaching steady state_bAnd taking the average value of the state samples of the sub-iteration period as a final state estimation value, namely:

has the advantages that: the invention processes the observation outlier of the sensor by modeling the model noise of the linear state space model as Student's t distribution, and processes the problem that the model noise parameter is unknown by taking the noise parameter as a random variable and randomly sampling the noise parameter and the system state simultaneously under the frame of a Gibbs sampler. The method can still obtain a better state estimation result when the sensor observation field value exists and the initial noise parameter setting error is larger.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 is a comparison of the position estimation root mean square error for the four methods of the present invention, i.e., GIB-RSTS versus KS, GIB-KS, and VB-KS;

FIG. 3 is a comparison of the root mean square error of velocity estimates for the four methods of the present invention, i.e., GIB-RSTS versus KS, GIB-KS, and VB-KS.

Detailed Description

Embodiment 1, referring to fig. 1, a robust adaptive smoothing method based on a Gibbs sampler includes the following steps:

A. modeling the process noise and the observation noise of the linear state space model into Student's t distribution, wherein the parameters of the process noise and the observation noise are random variables, and the prior distribution of the parameters of the process noise and the observation noise is the corresponding conjugate prior distribution.

For a linear state space model:

x_k＝F_k-1x_k-1+W_k-1

z_k＝H_kx_k+v_k

wherein:

is a vector of the states of the system,

in order to observe the vector for the system,

in order to be a state transition matrix,

in order for the system to observe the matrix,

in order to be a noise of the process,

to observe noise; w is a_kAnd v_kAre all modeled as Student's t, the distribution is as follows:

p(w_k)＝St(w_k；ξ，Q，ω)

p(v_k)＝St(v_k；μ，R，v)

wherein: st (x; mu, sigma, omega) represents a random vector x which satisfies Student's t distribution and takes mu as a mean vector, sigma as a scale matrix and omega as a degree of freedom; xi, Q, omega, mu, R and v are random variables, and the prior distribution of the random variables is modeled as corresponding conjugate prior; wherein, the prior distribution of mu and R is modeled as Gaussian-Inverse Wishart (GIW) distribution:

wherein: GIW (X, X; a, B, c, D) represents a Gaussian-inverse Wishart distribution with a, B, c, D as parameters; IW (X; a, B) represents that the random matrix X meets inverse Wisharp distribution with a as the degree of freedom and B as the inverse scale matrix; n (x; a, A) represents a Gaussian distribution in which the random vector x satisfies a as a mean and A as a variance; likewise, the prior distribution of ξ and Q is also modeled as a GIW distribution:

the degree of freedom parameters v and ω are both modeled as a Gamma distribution, as follows:

p(v)＝Gamma(v；e₀，f₀)

p(ω)＝Gamma(ω；g₀，h₀)

wherein: gamma (x; a, b) represents a random variable x satisfying the Gamma distribution with a as a shape parameter and b as a rate parameter; setting nominal values of mean vectors of process noise and observation noise to be xi respectively₀And mu₀(ii) a Order:

E(ξ)＝ξ₀

E(μ)＝μ₀

wherein: e (-) is the expectation of a random variable; to obtain:

α₀and beta₀Is a modulation parameter; setting nominal values of process noise and observation noise scale matrix to be Q respectively₀And R₀(ii) a Order:

E(Q)＝Q₀

E(R)＝R₀

to obtain:

t₀＝ρ_t+2n+2

T₀＝ρ_tQ₀

u₀＝ρ_u+2m+2

U₀＝ρ_uR₀

wherein: rho_uAnd rho_tIs a modulation parameter;

setting nominal values of process noise and observation noise freedom as omega₀And v₀(ii) a Order:

E(ω)＝ω₀

E(v)＝v₀

to obtain:

B. and introducing auxiliary parameters, and decomposing the process noise and the observation noise which meet the Student's t distribution into a combination of Gauss distribution and Gamma distribution by using a heuristic Gauss model.

Introducing an auxiliary parameter gamma_kAnd lambda_kThe process noise and the observation noise are decomposed into:

C. and initially sampling process noise, observation noise parameters and auxiliary parameters of the linear state space model under the frame of a Gibbs sampler.

According to the process noise scale matrix and the observation noise scale matrix prior distribution parameter t obtained in the step A₀，T₀，u₀，U₀Directly from IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Obtaining initial sampling of a process noise scale matrix and an observation noise scale matrix by intermediate sampling; namely:

Q⁽¹⁾～IW(Q；t₀，T₀)

R⁽¹⁾～IW(R；u₀，U₀)

wherein: a. the^(j)Representing the sampling value of the random variable A in the jth iteration;

the prior distribution parameter alpha of the process noise mean vector and the observation noise mean vector obtained according to the step A₀，β₀，

And slave IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Q of middle sampling⁽¹⁾And R⁽¹⁾Directly from

And

obtaining a process noise mean vector and an observation noise mean vector initial sample by intermediate sampling; namely:

according to the process noise degree of freedom and observation noise degree of freedom prior distribution parameter e obtained in the step A₀，f₀，g₀，h₀Directly from Gamma (omega; g)₀，h₀) And Gamma (v; e.g. of the type₀，f₀) Obtaining initial sampling of process noise freedom and observation noise freedom through intermediate sampling; namely:

ω⁽¹⁾～Gamma(ω；g₀，h₀)

v⁽¹⁾～Gamma(v；e₀，f₀)

according to the process noise auxiliary parameter and observation noise auxiliary parameter prior distribution parameter obtained in the step B, directly obtaining the noise auxiliary parameter

And

obtaining initial sampling of the process noise auxiliary parameter and the observation noise auxiliary parameter by intermediate sampling; namely:

D. setting a total number of iteration cycles N and a number of stationary cycles N_bIs required to satisfy N_b< N; and performing N iterations, wherein the flow of each iteration is as follows:

D1. and calculating the state posterior distribution of the current iteration period based on a Kalman smoother under the condition of the noise parameter and the auxiliary parameter sampled in the previous iteration period, and sampling the system state of the current iteration period from the posterior distribution.

Noise mean vector sample value xi during previous iteration cycle⁽ⁱ⁾Scale matrix sample value Q⁽ⁱ⁾And auxiliary parameter sampling value gamma_k ⁽ⁱ⁾Observation of noise mean vector sample value mu⁽ⁱ⁾Scale matrix sample value R⁽ⁱ⁾And an auxiliary parameter lambda_k ⁽ⁱ⁾Under the known condition, calculating a state posterior distribution parameter by adopting a Kalman smoother;

1) initializing system state and variance:

wherein: the subscript i | j represents the estimate of the variable at time i conditioned on the system observations prior to time j;

2) forward recursion:

let T be the total number of times, for k 1, 2, 3.

a) Prediction

Calculating the prior state and the prior variance at the k moment according to the posterior state and the posterior variance at the k-1 moment as follows:

P_k|k-1＝F_k-1P_k-1|k-1F_k-1 ^T+Q/γ_k

b) updating

Calculating the posterior state and posterior variance at the moment k as follows:

K_k＝P_k|k-1H_k ^T(H_kP_k|k-1H_k ^T+R/λ_k)

P_k|k＝P_k|k-1-K_kH_kP_k|k-1

wherein: k_kIs Kalman gain;

3) recursion backwards

For k ═ T-1, T-2.. 0, the following steps are performed:

P_k|T＝P_k|k+G_k(P_k+1|T-P_k+1|k)G_k ^T

wherein: g_kIs the smoothing gain;

according to the above Kalman smoother, xi is obtained⁽ⁱ⁾，Q⁽ⁱ⁾，γ_1：T ⁽ⁱ⁾，μ⁽ⁱ⁾，R⁽ⁱ⁾And lambda_1：T ⁽ⁱ⁾A state posterior parameter of a condition

And

further, for k equal to 0, 1, 2, 3.. T, the following are performed in order:

sampling is performed.

D2. And calculating posterior distribution of the observation noise mean value and the scale matrix under the condition of the system state, the observation variable and the observation noise auxiliary parameter sampled in the current iteration period, and sampling the observation noise mean value and the scale matrix in the current iteration period from the posterior distribution.

According to the Bayes rule:

p(μ，R|z_1：T，x_0：T，λ_1：T)＝cp(μ，R)p(z_1：T|μ，R，x_0：T，λ_1：T)

its logarithmic form is:

log p(μ，R|z_1：T，x_0：T，λ_1：T)＝log p(μ，R)+log p(z_1：T|μ，R，x_0：T，λ_1：T)+log c

as shown in the step a, the prior distribution of the observation noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the gaussian-inverse Wishart distribution is as follows:

the log form of the likelihood distribution can be written according to the observation model:

and calculating to obtain:

wherein:

the posterior distribution of the mean vector of the observed noise is sufficient

Is an average value of

The posterior probability density of the scale matrix is based on the Gaussian distribution of the variance

For the degree of freedom, to

The inverse Wishart distribution of the inverse scale matrix is obtained; in the actual sampling process, the following steps are required:

sample R⁽ⁱ⁺¹⁾And further according to:

sampling mu⁽ⁱ⁺¹⁾And obtaining an observation noise mean vector and a scale matrix sampling value of the current iteration period.

D3. And calculating posterior distribution of the process noise mean value and the scale matrix by taking the system state and the process noise auxiliary parameters sampled in the current iteration period as conditions, and sampling the process noise mean value and the scale matrix in the current iteration period from the posterior distribution.

According to the Bayes rule:

p(ξ，Q|x_0：T)＝cp(ξ，Q)p(x_0：T|ξ，Q)

its logarithmic form is:

log p(ξ，Q|x_0：T)＝log p(ξ，Q)+log p(x_0：T|ξ，Q)+log c

as shown in the step a, the prior distribution of the process noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the prior distribution is as follows:

the log form of the likelihood distribution is written according to a kinematic model:

and calculating to obtain:

wherein:

process noise mean satisfied with

Is an average value of

Is a Gaussian distribution of variance, a scale matrixHas a posterior probability density of

For the degree of freedom, to

The inverse Wishart distribution of the inverse scale matrix is obtained; in the actual sampling process, the following steps are required:

sample Q⁽ⁱ⁺¹⁾And further according to:

sample xi⁽ⁱ⁺¹⁾And obtaining a process noise mean vector and a scale matrix sampling value of the current iteration cycle.

D4. And calculating posterior distribution of the observation noise auxiliary parameters under the condition of the system state, the observation variable and the observation noise parameters sampled in the current iteration period, and sampling the observation noise auxiliary parameters in the current iteration period from the posterior distribution.

According to the Bayes rule:

p(λ_k|x_k，z_k，μ，R，v)＝cp(λ_k|v)p(z_k|x_k，λ_k，μ，R)

its logarithmic form is:

log p(λ_k|x_k，z_k，μ，R，v)＝log p(λ_k|v)+log p(z_k|x_k，λ_k，μ，R)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

the posterior distribution of the noise auxiliary parameter is observed sufficiently

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling lambda_k ⁽ⁱ⁺¹⁾And k is 1, 2,. T, and the sampling value of the observation noise auxiliary parameter of the current iteration period is obtained.

D5. And calculating posterior distribution of the process noise auxiliary parameters by taking the system state and the process noise parameters sampled in the current iteration period as conditions, and sampling the process noise auxiliary parameters in the current iteration period from the posterior distribution.

According to the Bayes rule:

p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝cp(γ_k-1|ω)p(x_k|x_k-1，γ_k-1，ξ，Q)

its logarithmic form is:

log p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝log p(γ_k-1|ω)+log p(x_k|x_k-1，γ_k-1，ξ，Q)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

process noise auxiliary parameter posterior distribution satisfied by

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling gamma_k-1 ⁽ⁱ⁺¹⁾And k is 1, 2.. T, and a process noise auxiliary parameter sampling value of the current iteration cycle is obtained.

D6. And calculating posterior distribution of the observation noise freedom degree under the condition of the observation noise auxiliary parameters sampled in the current iteration period, and sampling the observation noise freedom degree of the current iteration period from the posterior distribution.

According to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(v)＝(e₀-1)log v-f₀v+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

observing the posterior distribution of the degree of freedom of the noise to satisfy

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sample v⁽ⁱ⁺¹⁾And obtaining an observation noise freedom degree sampling value of the current iteration period.

D7. And calculating posterior distribution of the process noise freedom degree by taking the process noise auxiliary parameters sampled in the current iteration period as conditions, and sampling the process noise freedom degree of the current iteration period from the posterior distribution.

According to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(ω)＝(g₀-1)log ω-h₀ω+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

process noise degree of freedom posterior distribution satisfied by

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling omega⁽ⁱ⁺¹⁾And obtaining a process noise freedom degree sampling value of the current iteration period.

E. After N iterations, selecting N-N after reaching steady state_bAnd taking the average value of the state samples in the secondary iteration period as a final state estimation value.

After N times of iterative sampling, a state variable sampling set is obtained as

Selecting N-N after reaching steady state_bAnd taking the average value of the state samples of the sub-iteration period as a final state estimation value, namely:

embodiment 2, the pseudo code implemented by the present invention is:

embodiment 3, verification is performed by simulation data using the robust adaptive smoothing method described in embodiment 1.

The simulation case is two-dimensional target tracking and defines a state vector x_k＝[x_k y_k v_x，k v_y，k]^TWherein x is_k，y_k，v_x，kAnd v_y，kRespectively, the x and y direction positions and the x and y direction velocities. The observed variables are the positions in the x and y directions. The system kinematics model and the observation model are respectively as follows:

wherein: Δ t is a discrete time interval, chosen to be 1 second in this simulation. The total simulation time T was chosen to be 200 seconds. To simulate the observed outliers of the sensor, we use:

to generate true process noise and observation noise, wherein:

the nominal process noise scale matrix and the observation noise scale matrix are respectively set to be Q_n＝5Q_t，R_n＝R_t/5. The nominal process noise and the observed noise mean vector are both set to zero vectors. The simulated modulation parameters are set as follows: rho_t＝1，ρ_u＝1，α₀＝0.1，β₀＝0.1，e₀＝5，f₀＝1，g₀＝5，h₀1, total iteration cycle N5000, stationary cycle N_b1000. By way of comparison, the present embodiment simultaneously shows the estimation performance of the Kalman smoother estimation result (KS) solved by the nominal noise parameters, the adaptive Kalman smoother (gibb-KS) based on the Gibbs sampler and the noise gaussian distribution, and the robust adaptive Kalman smoother estimation result (VB-RSTS) based on the variational bayesian approximation. The method of the present invention is abbreviated as "GIB-RSTS".

The results of 500 independent sub-Monte Carlo simulations were used to verify the proposed method. The Root Mean Square Error (RMSE) of position and velocity and the Average Root Mean Square Error (ARMSE) are used to evaluate the performance of the different estimators. Several evaluation indices were calculated as follows:

wherein:

and

the real coordinate of the aircraft at the moment k in the ith Monte Carlo simulation;

and

position coordinates of the aircraft estimated at the moment k in the ith Monte Carlo simulation;

and

the real speed of the aircraft at the moment k in the ith Monte Carlo simulation;

and

the estimated speed of the aircraft at the moment k in the ith Monte Carlo simulation. M-500 indicates the total number of simulations.

FIG. 2 is a comparison of position estimation RMS errors for four methods, while FIG. 3 is a comparison of velocity RMS errors for four methods. The mean root mean square error for the four methods is shown in table 1.

TABLE 1

From fig. 2, fig. 3 and table 1, it can be seen that the GIB-RSTS proposed by the present invention can obtain better state estimation accuracy, including position estimation and velocity estimation, than KS, GIB-KS and VB-KS. KS does not consider noise parameter setting error and observation field value, so that the state estimation error is larger. The method provided by the invention is superior to GIB-KS, mainly because GIB-KS only considers the unknown of noise parameters and ignores the influence of the thick tail of noise. The main reason why the method of the invention is superior to VB-RSTS is that VB-RSTS is based on the free factorization approximation of posterior distribution and has certain approximation error. Generally speaking, the proposed GIB-RSTS can still obtain a more ideal state estimation result under the condition of large field value and noise parameter setting errors, and has better practical application potential.

Although the invention has been described in detail above with reference to a general description and specific examples, it will be apparent to one skilled in the art that modifications or improvements may be made thereto based on the invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.

Claims (6)

1. A robust adaptive smoothing method based on a Gibbs sampler is characterized by comprising the following steps:

A. modeling process noise and observation noise of a linear state space model into Student's t distribution, wherein parameters of the process noise and the observation noise are random variables, and the prior distribution of the parameters of the process noise and the observation noise is the corresponding conjugate prior distribution;

B. introducing auxiliary parameters, and decomposing process noise and observation noise meeting the Student's t distribution into a combination of Gauss distribution and Gamma distribution by using a heuristic Gauss model;

C. initially sampling process noise, observation noise parameters and auxiliary parameters of a linear state space model under the frame of a Gibbs sampler;

D. setting a total number of iteration cycles N and a number of stationary cycles N_bIs required to satisfy N_b< N; and performing N iterations, wherein the flow of each iteration is as follows:

D1. calculating the state posterior distribution of the current iteration cycle based on a Kalman smoother under the condition of the noise parameter and the auxiliary parameter sampled in the previous iteration cycle, and sampling the system state of the current iteration cycle from the posterior distribution;

D2. calculating posterior distribution of an observation noise mean value and a scale matrix under the condition of a system state, an observation variable and an observation noise auxiliary parameter sampled in the current iteration period, and sampling the observation noise mean value and the scale matrix in the current iteration period from the posterior distribution;

D3. calculating posterior distribution of a process noise mean value and a scale matrix by taking a system state and a process noise auxiliary parameter sampled in a current iteration period as conditions, and sampling the process noise mean value and the scale matrix of the current iteration period from the posterior distribution;

D4. calculating posterior distribution of the observation noise auxiliary parameters under the condition of the system state, the observation variable and the observation noise parameters sampled in the current iteration period, and sampling the observation noise auxiliary parameters in the current iteration period from the posterior distribution;

D5. calculating posterior distribution of the process noise auxiliary parameters by taking the system state and the process noise parameters sampled in the current iteration period as conditions, and sampling the process noise auxiliary parameters in the current iteration period from the posterior distribution;

D6. calculating posterior distribution of the observation noise freedom degree under the condition of observation noise auxiliary parameters sampled in the current iteration period, and sampling the observation noise freedom degree of the current iteration period from the posterior distribution;

D7. calculating posterior distribution of the process noise freedom degree under the condition of the process noise auxiliary parameters sampled in the current iteration period, and sampling the process noise freedom degree of the current iteration period from the posterior distribution;

E. after N iterations, selecting N-N after reaching steady state_bAnd taking the average value of the state samples in the secondary iteration period as a final state estimation value.

2. The method of claim 1, wherein the step a comprises the following steps:

for a linear state space model:

x_k＝F_k-1x_k-1+w_k-1

z_k＝H_kx_k+v_k

wherein:

is a vector of the states of the system,

in order to observe the vector for the system,

in order to be a state transition matrix,

in order for the system to observe the matrix,

in order to be a noise of the process,

to observe noise; w is a_kAnd v_kBoth modeled as Student's t distributions, as follows:

p(w_k)＝St(w_k；ξ，Q，ω)

p(v_k)＝St(v_k；μ，R，v)

wherein: st (x; mu, sigma, omega) represents a random vector x which satisfies Student's t distribution and takes mu as a mean vector, sigma as a scale matrix and omega as a degree of freedom; xi, Q, omega, mu, R and v are random variables, and the prior distribution of the random variables is modeled as corresponding conjugate prior; wherein, the prior distribution of mu and R is modeled as Gaussian-Inverse Wishart (GIW) distribution:

wherein: GIW (X, X; a, B, c, D) represents a Gaussian-inverse Wishart distribution with a, B, c, D as parameters; IW (X; a, B) represents that the random matrix X meets inverse Wisharp distribution with a as the degree of freedom and B as the inverse scale matrix; n (x; a, A) represents a Gaussian distribution in which the random vector x satisfies a as a mean and A as a variance; likewise, the prior distribution of ξ and Q is also modeled as a GIW distribution:

the degree of freedom parameters v and ω are both modeled as a Gamma distribution, as follows:

p(v)＝Gamma(v；e₀，f₀)

p(ω)＝Gamma(ω；g₀，h₀)

wherein: gamma (x; a, b) represents a random variable x satisfying the Gamma distribution with a as a shape parameter and b as a rate parameter; setting nominal values of mean vectors of process noise and observation noise to be xi respectively₀And mu₀(ii) a Order:

E(ξ)＝ξ₀

E(μ)＝μ₀

wherein: e (-) is the expectation of a random variable; to obtain:

α₀and beta₀Is a modulation parameter; setting nominal values of process noise and observation noise scale matrix to be Q respectively₀And R₀(ii) a Order:

E(Q)＝Q₀

E(R)＝R₀

to obtain:

t₀＝ρ_t+2n+2

T₀＝ρ_tQ₀

u₀＝ρ_u+2m+2

U₀＝ρ_uR₀

wherein: rho_uAnd rho_tIs a modulation parameter;

setting nominal values of process noise and observation noise freedom as omega₀And v₀(ii) a Order:

E(ω)＝ω₀

E(v)＝v₀

to obtain:

3. the Gibbs sampler-based robust adaptive smoothing method as claimed in claim 2, wherein the step B comprises:

introducing an auxiliary parameter gamma_kAnd lambda_kThe process noise and the observation noise are decomposed into:

4. the Gibbs sampler-based robust adaptive smoothing method as claimed in claim 3, wherein the step C comprises:

according to the process noise scale matrix and the observation noise scale matrix prior distribution parameter t obtained in the step A₀，T₀，u₀，U₀Directly from IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Obtaining initial sampling of a process noise scale matrix and an observation noise scale matrix by intermediate sampling; namely:

Q⁽¹⁾～IW(Q；t₀，T₀)

R⁽¹⁾～IW(R；u₀，U₀)

wherein: a. the^(j)Representing the sampling value of the random variable A in the jth iteration;

the prior distribution parameter alpha of the process noise mean vector and the observation noise mean vector obtained according to the step A₀，β₀，

And slave IW (Q; t)₀，T₀)，IW(R；u₀，U₀) Q of middle sampling⁽¹⁾And R⁽¹⁾Directly from

And

obtaining a process noise mean vector and an observation noise mean vector initial sample by intermediate sampling; namely:

obtaining the prior distribution parameter e of the process noise freedom degree and the observation noise freedom degree of the mountain according to the step A₀，f₀，g₀，h₀Directly from Gamma (omega; g)₀，h₀) And Gamma (v; e.g. of the type₀，f₀) Obtaining initial sampling of process noise freedom and observation noise freedom through intermediate sampling; namely:

ω⁽¹⁾～Gamma(ω；g₀，h₀)

v⁽¹⁾～Gamma(v；e₀，f₀)

the process noise auxiliary parameter obtained according to the step BAnd observing the noise-aided parameter prior distribution parameter directly from

And

obtaining initial sampling of the process noise auxiliary parameter and the observation noise auxiliary parameter by intermediate sampling; namely:

5. the Gibbs sampler-based robust adaptive smoothing method as claimed in claim 4, wherein the step D1 is implemented by:

noise mean vector sample value xi during previous iteration cycle⁽ⁱ⁾Scale matrix sample value Q⁽ⁱ⁾And auxiliary parameter sampling value gamma_k ⁽ⁱ⁾Observation of noise mean vector sample value mu⁽ⁱ⁾Scale matrix sample value R⁽ⁱ⁾And an auxiliary parameter lambda_k ⁽ⁱ⁾Under the known condition, calculating a state posterior distribution parameter by adopting a Kalman smoother;

1) initializing system state and variance:

wherein: the subscript i | j represents the estimate of the variable at time i conditioned on the system observations prior to time j;

2) forward recursion:

let T be the total number of times, for k 1, 2, 3.

a) Prediction

Calculating the prior state and the prior variance at the k moment according to the posterior state and the posterior variance at the k-1 moment as follows:

P_k|k-1＝F_k-1P_k-1|k-1F_k-1 ^T+Q/γ_k

b) updating

Calculating the posterior state and posterior variance at the moment k as follows:

K_k＝P_k|k-1H_k ^T(H_kP_k|k-1H_k ^T+R/λ_k)

P_k|k＝P_k|k-1-K_kH_kP_k|k-1

wherein: k_kIs Kalman gain;

3) recursion backwards

For k ═ T-1, T-2.. 0, the following steps are performed:

P_k|T＝P_k|k+G_k(P_k+1|T-P_k+1|k)G_k ^T

wherein: g_kIs the smoothing gain;

according to the above Kalman smoother, xi is obtained⁽ⁱ⁾，Q⁽ⁱ⁾，γ_1：T ⁽ⁱ⁾，μ⁽ⁱ⁾，R⁽ⁱ⁾And lambda_1：T ⁽ⁱ⁾A state posterior parameter of a condition

And

further, for k equal to 0, 1, 2, 3.. T, the following are performed in order:

sampling is carried out;

the method adopted in the step D2 is as follows:

according to the Bayes rule:

p(μ，R|z_1：T，x_0：T，λ_1：T)＝cp(μ，R)p(z_1：T|μ，R，x_0：T，λ_1：T)

its logarithmic form is:

log p(μ，R|z_1：T，x_0：T，λ_1：T)＝log p(μ，R)+log p(z_1：T|μ，R，x_0：T，λ_1：T)+log c

as shown in the step a, the prior distribution of the observation noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the gaussian-inverse Wishart distribution is as follows:

the log form of the likelihood distribution can be written according to the observation model:

and calculating to obtain:

wherein:

the posterior distribution of the mean vector of the observed noise is sufficient

Is an average value of

The posterior probability density of the scale matrix is based on the Gaussian distribution of the variance

For the degree of freedom, to

Is of inverse scaleInverse Wishart distribution of the matrix; in the actual sampling process, the following steps are required:

sample R⁽ⁱ⁺¹⁾And further according to:

sampling mu⁽ⁱ⁺¹⁾Obtaining an observation noise mean vector and a scale matrix sampling value of the current iteration cycle;

the method adopted in the step D3 is as follows:

according to the Bayes rule:

p(ξ，Q|x_0：T)＝cp(ξ，Q)p(x_0：T|ξ，Q)

its logarithmic form is:

log p(ξ，Q|x_0：T)＝log p(ξ，Q)+log p(x_0：T|ξ，Q)+log c

as shown in the step a, the prior distribution of the process noise mean and the scale matrix is modeled as gaussian-inverse Wishart distribution, and the logarithmic form of the prior distribution is as follows:

the log form of the likelihood distribution is written according to a kinematic model:

and calculating to obtain:

wherein:

process noise mean satisfied with

Is an average value of

The posterior probability density of the scale matrix is based on the Gaussian distribution of the variance

For the degree of freedom, to

The inverse Wishart distribution of the inverse scale matrix is obtained; in the actual sampling process, the following steps are required:

sample Q⁽ⁱ⁺¹⁾And further according to:

sample xi⁽ⁱ⁺¹⁾Obtaining a process noise mean vector and a scale matrix sampling value of the current iteration cycle;

the method adopted in the step D4 is as follows:

according to the Bayes rule:

p(λ_k|x_k，z_k，μ，R，v)＝cp(λ_k|v)p(z_k|x_k，λ_k，μ，R)

its logarithmic form is:

log p(λ_k|x_k，z_k，μ，R，v)＝log p(λ_k|v)+log p(z_k|x_k，λ_k，μ，R)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

the posterior distribution of the noise auxiliary parameter is observed sufficiently

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling lambda_k ⁽ⁱ⁺¹⁾Obtaining an observation noise auxiliary parameter sampling value of the current iteration period, wherein k is 1, 2,. T;

the method adopted in the step D5 is as follows:

according to the Bayes rule:

p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝cp(γ_k-1|ω)p(x_k|x_k-1，γ_k-1，ξ，Q)

its logarithmic form is:

log p(γ_k-1|x_k，x_k-1，ξ，Q，ω)＝log p(γ_k-1|ω)+log p(x_k|x_k-1，γ_k-1，ξ，Q)+log c

defined by the auxiliary parameters:

according to the observation model, the following results are obtained:

and further:

wherein:

process noise auxiliary parameter posterior distribution satisfied by

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling gamma_k-1 ⁽ⁱ⁺¹⁾Obtaining a process noise auxiliary parameter sampling value of the current iteration period, wherein k is 1, 2,. T;

the method adopted in the step D6 is as follows:

according to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(v)＝(e₀-1)log v-f₀v+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

observing the posterior distribution of the degree of freedom of the noise to satisfy

As a shape parameter, of

Is the Gamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sample v⁽ⁱ⁺¹⁾Obtaining an observation noise freedom degree sampling value of the current iteration period;

the method adopted in the step D7 is as follows:

according to the Bayes rule:

its logarithmic form is:

obtaining the following according to the prior distribution of the degree of freedom parameters:

log p(ω)＝(g₀-1)logω-h₀ω+c₁

according to the auxiliary parameter definition and the Stirling approximation:

and further:

wherein:

process noise degree of freedom posterior distribution satisfied by

As a shape parameter, of

To be fastGamma distribution of the rate parameter; in the actual sampling process, the following steps are required:

sampling omega⁽ⁱ⁺¹⁾And obtaining a process noise freedom degree sampling value of the current iteration period.

6. The method of claim 5, wherein the step E comprises the following steps:

after N times of iterative sampling, a state variable sampling set is obtained as

Selecting N-N after reaching steady state_bAnd taking the average value of the state samples of the sub-iteration period as a final state estimation value, namely:

Images