Robust Kalman filtering method based on student t mixed distribution

Abstract

The invention discloses a filtering method based on a student t mixed distribution model, and aims to effectively process non-stationary heavy tail noise. A traditional Kalman filter has great limitation when dealing with non-Gaussian noise, especially noise with a heavy tail characteristic. According to the RSTMKF provided by the invention, a plurality of student t distributions are introduced, so that the limitation of a traditional method is overcome, and modeling can be effectively carried out, and multi-modal and asymmetric noise distribution can be effectively processed. The invention further provides a variational Bayesian inference scheme for approaching the student t mixed distribution model, so that the balance between the calculation efficiency and the estimation precision is realized. A simulation experiment verifies that when the RSTMKF algorithm provided by the invention is used for processing non-stationary heavy-tail noise, compared with the prior art, the estimation precision can be remarkably improved, and the robustness can be enhanced.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of multi-source information fusion and filtering, and particularly relates to a robust Kalman filtering method based on student t mixed distribution. The method can be widely applied in unmanned navigation systems, aerospace fields, intelligent control systems of gas turbines and other fields. BACKGROUND

[0002] In complex environments such as unmanned navigation, aerospace and intelligent control of gas turbines, noise often exhibits non-Gaussian characteristics, which severely limits the performance of existing filtering methods. In particular, when the noise has heavy-tailed characteristics, the effectiveness of traditional Kalman filtering methods is significantly reduced, making it difficult to meet the needs of practical applications. Therefore, how to effectively handle complex noise models has become a technical problem to be solved.

[0003] The existence of non-Gaussian noise leads to a non-standard and complex form of the likelihood function. For example, when the noise has significant heavy-tailed characteristics (such as student t distribution), the posterior distribution may exhibit multi-modality, asymmetry or heavy-tailed characteristics. This complexity makes integral calculations (such as expectation or marginalization calculations) difficult, often unable to be directly solved by analytical methods. Therefore, approximation methods are particularly important in such problems. Existing solutions include Markov Chain Monte Carlo (MCMC), Variational Bayesian Inference, Laplace Approximation, Particle Filtering (PF), M-estimation and Gaussian Sum Approximation methods. These methods can address the above challenges to some extent, but each has its limitations.

[0004] Among them, MCMC and PF algorithms can approximate complex posterior distributions, but their high computational complexity is usually difficult to implement in practical applications, and they can only converge to the true probability density function when the number of particles tends to infinity. Laplace approximation involves a quadratic approximation of the posterior distribution at its maximum value, but since it relies on Gaussian assumptions, it may produce large errors when the posterior distribution has significant non-Gaussian characteristics (such as heavy-tailed, asymmetric or multi-modal). Gaussian Sum Approximation methods can theoretically handle heavy-tailed distributions, but they often require a large number of Gaussian components to accurately approximate the distribution, significantly increasing the computational complexity. M-estimation methods, such as robust Kalman filtering based on Huber loss function, maximum correlation entropy criterion Kalman filtering and covariance matching method, enhance the robustness of the system to outliers and noise. However, these methods usually assume that the noise is light-tailed, and may not fully capture extreme values in heavy-tailed distributions.

[0005] To address these challenges, the student t process is introduced as an alternative to the Gaussian process for noise modeling. For example, non-Gaussian measurement noise can be modeled by a t-distributed process, while t-distributed models are used for both measurement noise and process noise, and Gaussian-t mixture models are also used to describe the probability density function of noise. However, the use of a single t-distribution may require auxiliary strategies to improve matching accuracy in some cases. Although the Gaussian-t mixture model contains a t-distribution component, the Gaussian component in it may limit the flexibility of the model when dealing with extreme data. In addition, heavy-tailed distributions may not only occur in measurement noise, but also in process noise. Therefore, an innovative method is needed to more effectively address these issues.

[0006] To this end, the present application proposes a new robust Kalman filter method based on student t mixture distribution to effectively address the problem of non-stationary heavy-tailed noise. By combining multiple student t distributions, the model proposed by the present application can capture heavy-tailed characteristics and more robustly handle multimodal and asymmetric distributions in practical applications, while balancing computational complexity and estimation accuracy. Compared with existing Kalman filters based on Gaussian-t mixture distribution and robust student t Kalman filters with fixed distribution parameters, the present application has higher flexibility and adaptability, and can more effectively handle complex non-stationary heavy-tailed noise. SUMMARY

[0007] The present application proposes a new robust Kalman filter method based on student t mixture distribution to effectively address the problem of non-stationary heavy-tailed noise. By combining multiple student t distributions, the model proposed by the present application can capture heavy-tailed characteristics and more robustly handle multimodal and asymmetric distributions in practical applications, while balancing computational complexity and estimation accuracy.

[0008] The technical solution for achieving the object of the present application is: a robust Kalman filter method based on student t mixture distribution, comprising the following steps:

[0009] Step 1: input and initialize the state vector of the previous time , the measurement vector , the state vector dimension n, the measurement vector dimension m, the error covariance matrix , the state transition matrix , the process noise covariance matrix , the Gamma distribution parameters , , , and , , , Beta distribution parameters and Maximum number of iterations ;

[0010] Step 2: Solve for the one-step predicted state vector and the one-step predicted error covariance matrix based on the Kalman filter time update equation;

[0011] Step 3: Iterate and recursively update the variational Bayesian measurement;

[0012] Step 4: Output the variational Bayesian measurement update results to obtain the final posterior state vector. And error covariance matrix ;

[0013] Step 5: Write and debug the program, verify it through simulation, and finally apply the algorithm to the established state-space model to achieve robust filtering for non-stationary heavy-tailed noise problems.

[0014] In step two, the one-step predicted state vector and the one-step predicted error covariance matrix are solved according to the Kalman filter time update equation. The specific method is as follows:

[0015]

[0016] in, and Let these represent the one-step predicted state vector and the one-step predicted error covariance matrix, respectively. and Given the state estimate and error covariance matrix from the previous step, The process noise covariance matrix is... Let be the state transition matrix.

[0017] In step three, variational Bayesian measurement updates are performed iteratively, specifically as follows:

[0018] (1) Initialize the expected values ​​of the following parameters, including the Gamma random variable. , , , Each other's expectations , , , Bernoulli random variable , Each other's expectations , and mixed weights , Each other's expectations , .

[0019] (2) The error covariance matrix after introducing the mixed t-distribution and observation noise covariance matrix It is given by the following formula:

[0020]

[0021]

[0022] in, Indicates the first iteration Indicates the first The next iteration.

[0023] (3) The corrected error covariance matrix and observation noise covariance matrix Substituting the measurement update equation into the standard Kalman filter framework yields the corrected filter gain. Corrected state estimates and the corrected state estimation error covariance matrix It is given by the following formula:

[0024]

[0025]

[0026] in, It is a measurement matrix. Represents the measurement matrix Transpose.

[0027] (4) Calculate the state prediction error covariance and new information covariance It is given by the following formula:

[0028]

[0029]

[0030] (5) Calculate the probability , , and It is given by the following formula:

[0031]

[0032]

[0033]

[0034]

[0035] in, and As the normalization factor, This indicates that the trace of the matrix expression within the parentheses is being calculated.

[0036] (6) Calculate the Bernoulli random variable , Expectations and It is given by the following formula:

[0037]

[0038]

[0039] (7) Update the Gamma distribution parameters , , , , , , , It is given by the following formula:

[0040]

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047]

[0048] in, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameters, which affect the noise scale and The uncertainty in the model was modeled; similarly, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameters, which affect the measurement noise scale and The uncertainty in the model was modeled.

[0049] (8) Update the expected value , , , , , , , It is given by the following formula:

[0050]

[0051]

[0052]

[0053]

[0054]

[0055]

[0056]

[0057]

[0058] Among them, for random variables with Gamma distribution , , , Their expected values ​​are directly given by the parameters of the Gamma distribution, and the expected value of the logarithm can be expressed using the digamma function. It is calculated to be the derivative of the logarithmic Gamma function.

[0059] (9) Update the Beta-distributed random variable and Logarithmic expected value , , , It is given by the following formula:

[0060]

[0061]

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068] in, and It is a Beta distribution Shape parameters, and It is a Beta distribution The shape parameters, the parameters input and initialized in step one. and These are Beta distributions. and The shape parameters, and and It is both a Beta-distributed random variable and a mixed weight in the specific method (1) of step three.

[0069] (10) Each time the variational Bayesian measurement update as shown above is performed, a convergence judgment of the state estimation is made. If the following conditions are met... ( If the convergence threshold of the relative change in the monitored state estimate is reached, the iterative iteration ends; otherwise, the iterative iteration of variational Bayesian measurement updates continues until the convergence condition is met or the maximum number of iterations is reached. .

[0070] Compared with the prior art, the beneficial effects of the present invention are:

[0071] This invention proposes a novel robust Kalman filter method based on a mixture of Student's t distributions to effectively address non-stationary heavy-tailed noise problems. The proposed model captures heavy-tailed features and handles multimodal and asymmetric distributions more robustly in practical applications. Furthermore, it employs a variational Bayesian approximation method, transforming the problem of solving complex integrals into an optimization problem, thus avoiding the direct computation of these difficult integrals and achieving a balance between computational complexity and estimation accuracy. Compared to existing Kalman filters based on Gaussian-t mixture distributions and robust Student's t Kalman filters with fixed distribution parameters, this invention offers greater flexibility and adaptability, and can more effectively handle complex non-stationary heavy-tailed noise. Attached Figure Description

[0072] The advantages of the present invention described above will become apparent and readily understood from the following description of the simulation experiments in conjunction with the accompanying drawings, wherein:

[0073] Figure 1 The results are for simulation experiment scenario 1;

[0074] Figure 2 The results are for simulation experiment scenario 2;

[0075] Figure 3 The results are for simulation experiment scenario 3;

[0076] Figure 4 The results are for simulation experiment scenario 4; Detailed Implementation

[0077] The invention will now be further described with reference to the accompanying drawings.

[0078] This paper, based on the Student's t-mixture distribution, overcomes the problem of significant finiteness reduction in traditional Kalman filters when facing heavy-tailed noise by combining multiple Student's t-distributions. This enables the Kalman filter to effectively model and manage multimodal and asymmetric noise distributions. Furthermore, a comprehensive derivation of the variational Bayesian inference formula for approximating the Student's t-mixture distribution is provided, ensuring an optimal trade-off between computational efficiency and estimation accuracy. Finally, a robust Kalman filtering method based on the Student's t-mixture distribution is derived, comprising the following steps:

[0079] Step 1: Input and initialize the state vector from the previous time step. Measurement vector The state vector dimension is n, the measurement vector dimension is m, and the error covariance matrix is... State transition matrix Process noise covariance matrix Gamma distribution parameters , , , and , , , Beta distribution parameters and Maximum number of iterations ;

[0080] Step 2: Solve for the one-step predicted state vector and the one-step predicted error covariance matrix based on the Kalman filter time update equation;

[0081] Step 3: Iterate and recursively update the variational Bayesian measurement;

[0082] Step 4: Output the variational Bayesian measurement update results to obtain the final posterior state vector. And error covariance matrix ;

[0083] Step 5: Write and debug the program, verify it through simulation, and finally apply the algorithm to the established state-space model to achieve robust filtering for non-stationary heavy-tailed noise problems.

[0084] In step two, the one-step predicted state vector and the one-step predicted error covariance matrix are solved according to the Kalman filter time update equation. The specific method is as follows:

[0085]

[0086] in, and Let these represent the one-step predicted state vector and the one-step predicted error covariance matrix, respectively. and Given the state estimate and error covariance matrix from the previous step, The process noise covariance matrix is... Let be the state transition matrix.

[0087] In step three, variational Bayesian measurement updates are performed iteratively, specifically as follows:

[0088] (1) Initialize the expected values ​​of the following parameters, including the Gamma random variable. , , , Each other's expectations , , , Bernoulli random variable , Each other's expectations , and mixed weights , Each other's expectations , .

[0089] (2) The error covariance matrix after introducing the mixed t-distribution and observation noise covariance matrix It is given by the following formula:

[0090]

[0091]

[0092] in, Indicates the first iteration Indicates the first The next iteration.

[0093] (3) The corrected error covariance matrix and observation noise covariance matrix Substituting the measurement update equation into the standard Kalman filter framework yields the corrected filter gain. Corrected state estimates and the corrected state estimation error covariance matrix It is given by the following formula:

[0094]

[0095]

[0096] in, It is a measurement matrix. Represents the measurement matrix Transpose.

[0097] (4) Calculate the state prediction error covariance and new information covariance It is given by the following formula:

[0098]

[0099]

[0100] (5) Calculate the probability , , and It is given by the following formula:

[0101]

[0102]

[0103]

[0104]

[0105] in, and As the normalization factor, This indicates that the trace of the matrix expression within the parentheses is being calculated.

[0106] (6) Calculate the Bernoulli random variable , Expectations and It is given by the following formula:

[0107]

[0108]

[0109] (7) Update the Gamma distribution parameters , , , , , , , It is given by the following formula:

[0110]

[0111]

[0112]

[0113]

[0114]

[0115]

[0116]

[0117]

[0118] in, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameter, and It is both a random variable of Gamma and a measure of process noise; similarly, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameter, and It is both a random variable of Gamma and a measure of noise.

[0119] (8) Update the expected value , , , , , , , It is given by the following formula:

[0120]

[0121]

[0122]

[0123]

[0124]

[0125]

[0126]

[0127]

[0128] in, Let denote the digamma function, which is the derivative of the logarithmic Gamma function.

[0129] (9) Update the Beta-distributed random variable and Logarithmic expected value , , , It is given by the following formula:

[0130]

[0131]

[0132]

[0133]

[0134]

[0135]

[0136]

[0137]

[0138] in, and It is a Beta distribution Shape parameters, and It is a Beta distribution The shape parameters, the parameters input and initialized in step one. and These are Beta distributions. and The shape parameters, and and It is both a Beta-distributed random variable and a mixed weight in the specific method (1) of step three.

[0139] (10) Each time the variational Bayesian measurement update as shown above is performed, a convergence judgment of the state estimation is made. If the following conditions are met... ( If the convergence threshold of the relative change in the monitored state estimate is reached, the iterative iteration ends; otherwise, the iterative iteration of variational Bayesian measurement updates continues until the convergence condition is met or the maximum number of iterations is reached. .

[0140] Finally, to verify the effectiveness and robustness of the robust Kalman filter method based on the Student's t mixture distribution proposed in this invention in handling non-stationary heavy-tailed noise, simulation experiments were set up and compared with the following filters: the Maximum Correntropy Criterion Kalman Filter (MCCKF), the Huber-based Robust Kalman Filter (HRKF), the Robust Student's t-based Kalman Filter (RSTKF), the Gaussian–Student's t Mixture Distribution Based Kalman Filter (GSTMKF), and the Kalman Filter with dynamically adjusted noise covariance matrices (KFANM). For KFANM, the noise was dynamically adjusted based on the noise variation probability and scaling factors U1 and U2 in the simulation. For other filtering methods, only the initial noise covariance is set. Furthermore, the optimal Kalman filter (OKF) with the true noise covariance matrix is ​​used as a benchmark.

[0141] The linear discrete state-space model used in the simulation can be represented as:

[0142]

[0143]

[0144] in, Let k be the state vector at time k. For the corresponding measurement vector, and These are the state transition matrix and the measurement matrix, respectively. and These are process noise and measurement noise, respectively.

[0145] The relevant parameters involved in the experiment are as follows: Assume the coefficient matrix of the state-space model is as follows: , The nominal process noise and measurement noise covariances are respectively... and The initial state and the initial estimation error covariance matrix are set as follows: , The number of Monte Carlo simulations was [number missing]. .

[0146] Scenario 1: The non-stationary heavy-tailed process and measurement noise change over time as follows:

[0147]

[0148]

[0149] in , In scenario 1, the process noise and measurement noise follow different Gaussian mixture distributions with varying variances at different time intervals. This setting is designed to simulate heavy-tailed phenomena, thereby testing the estimation performance of various algorithms.

[0150] Scenario 2: In Scenario 2, the heavy-tailed phenomenon is further aggravated. The noise follows a Gaussian distribution mixture, with the same probability in each time interval as in Scenario 1, but... and .

[0151] Scenario 3: In Scenario 3, compared to Scenario 1, the covariance of the process noise is significantly amplified, while the measurement noise is amplified relatively less. At each time interval, the noise follows a Gaussian distribution mixture with the same probability as in Scenario 1, but the covariance scaling factor is different. , .

[0152] Case 4: In Case 4, compared to Case 3, except for the covariance scaling factor , Apart from that, all other conditions remain unchanged.

[0153] To evaluate the estimation performance of different methods, the root mean square error (RMSE) was chosen as the evaluation criterion.

[0154]

[0155] in, This is the state estimate. This represents the actual state value. The maximum number of iterations is set to... Simulated measurement data was collected over 400 seconds. To improve the differentiation in the comparison graphs, a moving smoothing method with a window size of 10 was applied to the data. The simulation results are as follows: Figures 1 to 4 As shown.

[0156] In scenarios 1 and 2, the RMSE curves of KFANM and MCCKF are quite similar within the time intervals of 1–100 seconds and 301–400 seconds, respectively, indicating a very low probability of anomalies causing mild heavy-tailed noise. The method proposed in this invention performs slightly better than other methods and is closer to the true value, while the performance of the HRKF method is moderate. RSTKF and GSTMKF have the highest RMSE values.

[0157] In the time intervals of 101–200 s and 201–300 s, the heavy-tailed distribution probabilities of process noise and measurement noise increased from 0.05 to 0.15 and from 0.1 to 0.2, respectively. Under these conditions, the performance of all filtering methods deteriorated, with the degree of difference proportional to the severity of the heavy-tailed phenomenon. The RMSE curve of the proposed method deviated more significantly from that of other methods, further demonstrating its superior performance in handling heavy-tailed noise.

[0158] For scenarios 3 and 4, the parameters are reconfigured to obtain better estimation results. The purpose of this reconfiguration is to demonstrate that the parameters of the filtering method should be adjusted according to the actual situation. Scenarios 1 and 2 mainly involve severe measurement noise tailing problems, while scenarios 3 and 4 involve severe process noise tailing problems. Obviously, it is highly unlikely that both scenarios will occur simultaneously, but it is possible to predict which scenario is closer to reality and adjust the parameters accordingly in advance. Figure 3 and Figure 4 Simulation results show that MCCKF performs the worst in terms of estimation accuracy. RSTKF and GSTMKF exhibit performance close to KFANM. However, the method proposed in this invention has higher estimation accuracy and remains robust and stable even under more severe heavy-tailed noise conditions.

Claims

1. A robust Kalman filter method based on a mixed distribution of student t, characterized in that, Includes the following steps: Step 1: Input and initialize the state vector from the previous time step. Measurement vector The state vector dimension is n, the measurement vector dimension is m, and the error covariance matrix is... State transition matrix Process noise covariance matrix Gamma distribution parameters , , , and , , , Beta distribution parameters and Maximum number of iterations ; Step 2: Solve for the one-step predicted state vector and the one-step predicted error covariance matrix based on the Kalman filter time update equation; Step 3: Iterate and recursively update the variational Bayesian measurement; Step 4: Output the variational Bayesian measurement update results to obtain the final posterior state vector. And error covariance matrix ; Step 5: Write and debug the program, verify it through simulation, and finally apply the algorithm to the established state-space model to achieve robust filtering for non-stationary heavy-tailed noise problems.

2. According to claim 1, the second step of a robust Kalman filtering method based on a mixed distribution of student t, which involves solving for the one-step predicted state vector and the one-step predicted error covariance matrix based on the Kalman filter time update equation, is characterized in that... The specific method is as follows: in, and Let these represent the one-step predicted state vector and the one-step predicted error covariance matrix, respectively. and Given the state estimate and error covariance matrix from the previous step, The process noise covariance matrix is... Let be the state transition matrix.

3. According to claim 1, the variational Bayesian measurement update performed by the cyclic recursion in step three of a robust Kalman filter method based on a student t-mixed distribution is characterized in that... The specific method is as follows: (1) Initialize the expected values ​​of the following parameters, including the Gamma random variable. , , , Each other's expectations , , , Bernoulli random variable , Each other's expectations , and mixed weights , Each other's expectations , . (2) The error covariance matrix after introducing the mixed t-distribution and observation noise covariance matrix It is given by the following formula: in, Indicates the first iteration Indicates the first The next iteration. (3) The corrected error covariance matrix and observation noise covariance matrix Substituting the measurement update equation into the standard Kalman filter framework yields the corrected filter gain. Corrected state estimates and the corrected state estimation error covariance matrix It is given by the following formula: in, It is a measurement matrix. Represents the measurement matrix Transpose. (4) Calculate the state prediction error covariance and new information covariance It is given by the following formula: (5) Calculate the probability , , and It is given by the following formula: in, and As the normalization factor, This indicates that the trace of the matrix expression within the parentheses is being calculated. (6) Calculate the Bernoulli random variable , Expectations and It is given by the following formula: (7) Update the Gamma distribution parameters , , , , , , , It is given by the following formula: in, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameter, and It is both a random variable of Gamma and a measure of process noise; similarly, and It is a Gamma distribution and Shape parameters, and It is a Gamma distribution and The scale parameter, and It is both a random variable of Gamma and a measure of noise. (8) Update the expected value , , , , , , , It is given by the following formula: in, Let denote the digamma function, which is the derivative of the logarithmic Gamma function. (9) Update the Beta-distributed random variable and Logarithmic expected value , , , It is given by the following formula: in, and It is a Beta distribution Shape parameters, and It is a Beta distribution The shape parameters, the parameters input and initialized in step one. and These are Beta distributions. and The shape parameters, and and It is both a Beta-distributed random variable and a mixed weight in the specific method (1) of step three. (10) Each time the variational Bayesian measurement update as shown above is performed, a convergence judgment of the state estimation is made. If the following conditions are met... ( If the convergence threshold of the relative change in the monitored state estimate is reached, the iterative iteration ends; otherwise, the iterative iteration of variational Bayesian measurement updates continues until the convergence condition is met or the maximum number of iterations is reached. .

Images