A method for state estimation of mobile robot based on Kalman filtering under thick tail noise

Abstract

The application relates to the field of mobile robots, in particular to a mobile robot state estimation method based on Kalman filtering under thick tail noise, which comprises the following steps: establishing a discrete time linear state space model of a mobile robot, including a state equation and an observation equation; and realizing state updating of the mobile robot through robust Kalman filtering based on a gamma-multivariate Laplace mixed distribution thick tail noise; the application exhibits obvious performance advantages under a non-stationary heavy tail noise environment.

Description

Technical Field

[0001] This invention relates to the field of mobile robots, and more specifically to a method for state estimation of mobile robots based on Kalman filtering under thick-tailed noise. Background Technology

[0002] Kalman filters (KFs) have been widely used in mobile robot state estimation due to their high computational efficiency and ease of implementation. However, outliers commonly found in real-world measurements often cause noise to exhibit heavy-tailed characteristics. Traditional KFs, based on the Gaussian noise assumption, struggle to accurately model such non-Gaussian noise, significantly degrading state estimation performance and potentially causing filter divergence under heavy-tailed interference. This limitation restricts the practical application value of KFs in complex dynamic environments. Therefore, research on robust state estimation methods for systems with heavy-tailed noise is of significant engineering importance.

[0003] Existing patent CN113971752A proposes using the Student's t-distribution to model non-Gaussian noise with outliers. This is achieved by approximating the posterior probability density function to a Student's t-distribution with fixed degrees of freedom, thereby improving estimation accuracy. However, in practice, simply using the Student's t-distribution to fit heavy-tailed noise suffers from limited accuracy and insufficient robustness, and its performance is highly dependent on the choice of the t-distribution's degrees of freedom. Besides the Student's t-distribution, the multivariate Laplace (ML) distribution provides another model for characterizing heavy-tailed non-Gaussian noise without introducing degrees of freedom parameters. For example, Wang et al. proposed a robust nonlinear Kalman filter (MLRCK) based on multivariate Laplace distribution, which combines multivariate Laplace distribution with variational Bayesian inference to achieve joint estimation of system state and noise characteristics (G. Wang, C. Yang, and X. Ma, “A novelrobust nonlinear Kalman filter based on multivariate Laplace distribution,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol. 68, no. 7, pp. 2705–2709, Jul. 2021.). However, when there are few or no outliers in the measured noise, the accuracy of noise modeling using multivariate Laplace distribution decreases, which also leads to a reduction in overall estimation performance.

[0004] Existing methods generally do not consider the case where the observed noise has a non-stationary thick-tailed distribution. When the thick-tailed noise exhibits time-varying or non-stationary characteristics, ignoring its non-stationarity will significantly reduce the estimation accuracy of the Kalman filter. Therefore, it is necessary to carry out research on robust state estimation for non-stationary thick-tailed noise. Summary of the Invention

[0005] To address the above problems, this invention provides a method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise, comprising the following steps:

[0006] S1. Establish a discrete-time linear state-space model of the mobile robot, including state equations and observation equations;

[0007] S2. The state update of the mobile robot is achieved by using a heavy-tailed noise robust Kalman filter based on a gamma-multivariate Laplace hybrid distribution.

[0008] The beneficial effects of this invention are:

[0009] This invention addresses the problem of non-stationary, heavy-tailed distribution of measurement noise in mobile robot systems by proposing a robust Kalman filter method for heavy-tailed noise based on a gamma-multivariate Laplace mixture (GaML) distribution. First, based on the multivariate Laplace distribution, auxiliary variables are introduced to achieve dynamic parameter adjustment: a gamma prior distribution is used as a scaling factor to adaptively calibrate the covariance matrix, and the mean vector is compensated and corrected using the product of the scaling factor and the shape parameter of the multivariate Laplace distribution. Second, a hierarchical Gaussian state-space model is constructed based on this theory, where both measurement noise and likelihood functions are modeled using GaML distributions, and the joint estimation of system state and distribution parameters is achieved through variational Bayesian methods. Simulation experiments show that the algorithm exhibits significant performance advantages in non-stationary heavy-tailed noise environments. Attached Figure Description

[0010] Figure 1 This is a flowchart of a mobile robot state estimation method based on Kalman filtering under thick-tailed noise according to the present invention.

[0011] Figure 2 This is a schematic diagram of the root mean square error of various algorithms in the embodiments of the present invention. Detailed Implementation

[0012] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0013] To address the issue that existing Kalman filters neglect noise non-stationarity, thus affecting estimation accuracy, this invention proposes a Kalman filter-based state estimation method for mobile robots under heavy-tailed noise. This method employs a robust Kalman filter for heavy-tailed noise based on a gamma-multivariate Laplace mixture (GaML) distribution. By introducing a gamma distribution with adjustable shape and scale parameters into the multivariate Laplace model, this robust Kalman filter method adaptively adjusts the mean vector and covariance matrix of the multivariate Laplace distribution. To support online parameter updates, a hierarchical Gaussian representation of GaML is constructed, and variational Bayesian inference is used to jointly estimate the system state and model parameters. This method significantly improves the state estimation accuracy of mobile robots under non-stationary conditions and heavy-tailed noise interference.

[0014] Please see Figure 1 In this embodiment of the invention, a method for estimating the state of a mobile robot based on Kalman filtering under thick-tailed noise includes the following steps:

[0015] S1. Establish a discrete-time linear state-space model of the mobile robot, including state equations and observation equations.

[0016] In some embodiments, the state equations and observation equations of the discrete-time linear state-space model are expressed as follows:

[0017] (1)

[0018] In the formula, Indicates the mobile robot at time The state vector, Indicates the mobile robot at time Given the state transition matrix, This indicates that it conforms to a Gaussian distribution. Process noise; This represents a multivariate normal distribution, also known as a Gaussian distribution. Indicates the mobile robot at time The observation vector, Indicates the mobile robot at time The known observation matrix; This indicates that the measurement noise follows a heavy-tailed characteristic with uncertainty. The dimension of the observation vector. Represents the set of real numbers. The dimension of the state vector. This represents the process noise covariance.

[0019] S2. The state update of the mobile robot is achieved by using a heavy-tailed noise robust Kalman filter based on a gamma-multivariate Laplace hybrid distribution.

[0020] Preferably, in engineering practice, measurement noise often exhibits non-stationary, heavy-tailed statistical characteristics due to the influence of time-varying interference factors. To address this issue, this invention innovatively proposes using the GaML distribution to model this type of noise. The probability density function of the GaML distribution is expressed as:

[0021] (2)

[0022] In the formula, This represents the position vector (i.e., the mean vector). Represents the scale matrix (i.e., the covariance matrix). Indicates shape parameters, Indicates the first auxiliary parameter. Indicates the rate parameter. Indicates the DOF parameter. Represents a random vector. Represents the gamma distribution, where , This is a gamma function. The multivariate Laplace distribution has the following form:

[0023] (3)

[0024] Indicates having a rate parameter The exponential distribution with a value of 1, , This indicates the second auxiliary parameter.

[0025] Combining equation (3), the probability density function of the GaML distribution can be further expressed as:

[0026] (4)

[0027] This represents a multivariate normal distribution.

[0028] This invention employs a gamma distribution. First auxiliary parameter As a scaling factor, the scaling matrix is ​​realized. Adaptive calibration, and utilizing the first auxiliary parameter Shape parameters of multivariate Laplace distribution The product of position vectors Compensation and corrections will be made.

[0029] The embodiments of the present invention will use shape parameters The model is a random variable, and it is assigned a Gaussian prior distribution, i.e.:

[0030] (5)

[0031] , These are the mean hyperparameter and variance hyperparameter of the Gaussian prior distribution, respectively.

[0032] Will , , Treating them as latent variables, based on equations (4) and (5), GaML can be expressed as the following Gaussian hierarchical model:

[0033] (6)

[0034] In the formula, Represents observed variables Given latent variables Multivariate normal distribution under certain conditions Representing latent variables The prior distribution, Representing latent variables The multivariate normal prior distribution, Representing latent variables The gamma prior distribution.

[0035] Preferably, step S2 specifically includes:

[0036] S21. Prediction Step: Obtain the mobile robot's position in the previous time step. posterior state estimation vector and the posterior error covariance matrix The one-step prediction probability density function It approximates a Gaussian distribution, and from this, we can deduce the mobile robot's current state at the current moment. Prior state estimation vector and prior error covariance matrix .

[0037] Specifically, at the current moment The one-step prediction probability density function can be approximated by a Gaussian distribution, that is:

[0038] (7)

[0039] In the formula, This represents the time from moment 1 to the previous moment. All observation data.

[0040] Wherein, the prior state estimation vector and prior error covariance matrix Based on the state equations, it is expressed as:

[0041] (8)

[0042] In the formula, This represents the process noise covariance matrix.

[0043] S22. Bayesian Update: Get the Current Moment observation vector By combining the likelihood function constructed based on the gamma-multivariate Laplace mixture distribution, a joint posterior probability density function containing state and model noise parameters is formed. ; Represents the set of parameters to be estimated. Indicates the mobile robot at time The state vector, , Indicates auxiliary parameters, Indicates shape parameters.

[0044] Preferably, the GaML distribution proposed in this invention can be used to characterize measurement noise with uncertain heavy-tailed statistical properties, expressed as:

[0045] (9)

[0046] In the formula, Indicates measurement noise The probability density function, Indicates time The measurement noise vector, Indicates time Measurement noise covariance, This represents the DOF parameter.

[0047] According to equations (1) and (9), at time... The likelihood function constructed based on the gamma-multivariate Laplace mixture distribution is expressed as:

[0048] (10)

[0049] In the formula, Indicates the current time The likelihood probability density function.

[0050] Preferably, the joint posterior probability density function Represented as:

[0051] (11)

[0052] Furthermore,

[0053] (12)

[0054] In the formula, Represents the time from time 1 to the current time. All observation data; Indicates having a rate parameter The exponential distribution; Indicates having a rate parameter The exponential distribution with a value of 1, The vector representing the mean of a normal distribution is a representation of the mean of the distribution. The initial guess value, Represents constant real numbers. express An identity matrix of order 1.

[0055] S23. Variational Inference: Using variational Bayesian methods, the joint posterior probability density function is... It can be decomposed into the product of independent factors of the state and model noise parameters, i.e. And through alternating iterative optimization, it approximates the true posterior. This represents the approximate posterior probability density function.

[0056] Preferably, the joint posterior probability density function is typically difficult to solve directly analytically, while the Variational Bayesian (VB) method can efficiently approximate it. Its core idea lies in decomposing the joint posterior into the product of several independent factors and obtaining an approximate solution through iterative optimization. The specific expression is as follows:

[0057] (13)

[0058] Wherein, the approximate posterior probability density function The following constraints must be met:

[0059] (14)

[0060] In the formula, Indicates except All parameters except Expectations Denotes the set of parameters to be estimated. One of the parameters, It represents a constant. This represents the joint posterior probability density function.

[0061] According to formulas (7), (12) and (14), The formal expression is:

[0062] (15)

[0063] In the formula, Indicates and Irrelevant constant terms This represents the dimension of the observation vector.

[0064] Preferably, in step S23, after initializing the initial values ​​for the iteration, Z approximate posterior probability density functions are performed respectively. The update process includes:

[0065] S2301. Initialization: , , , This represents the initial expected value.

[0066] S2302. Calculation and :

[0067] (16)

[0068] In the formula, This represents the corrected measurement noise covariance. This represents the adjusted mean. Let represent the mathematical expectation calculated based on the approximate posterior distribution of the z-th iteration.

[0069] S2303 Update .

[0070] In one specific embodiment, let Based on formulas (14) and (15), we can obtain formula (17):

[0071] (17)

[0072] In the formula, This represents the mathematical expectation calculated based on the approximate posterior distribution of the z-th iteration. Indicates and Irrelevant constant terms This represents the mean of the equivalent Gaussian observations. This represents the covariance of equivalent Gaussian observations.

[0073] From the above formula, we can see that It is approximated by a Gaussian distribution, as shown below:

[0074] (18)

[0075] In the formula, Indicates the first The posterior state estimation vector under the next iteration. Indicates the first The posterior error covariance matrix under the next iteration is given by the following equation:

[0076] (19)

[0077] In the formula, express An identity matrix of order 1. Indicates the first Kalman gain in the next iteration.

[0078] S2304 Update .

[0079] In one specific embodiment, let Based on formulas (14) and (15), we can obtain formula (20):

[0080] (20)

[0081] In the formula, Indicates and Irrelevant constant terms , , These represent the third, fourth, and fifth auxiliary parameters, respectively.

[0082] From the above formula, we can see that It is approximated as a point distribution, as shown below:

[0083] (twenty one)

[0084] In the formula, Indicates the first The first auxiliary parameter in the next iteration is given by the following formula:

[0085] (twenty two)

[0086] (twenty three)

[0087] S2305. Calculation .

[0088] S2306 Update .

[0089] In one specific embodiment, let Based on formulas (14) and (15), we can obtain formula (24):

[0090] (twenty four)

[0091] In the formula, Indicates and Irrelevant constant terms.

[0092] From (24), we can see that It is approximated by a Gaussian distribution, as shown below:

[0093] (25)

[0094] (26)

[0095] In the formula, Indicates the sixth auxiliary parameter No. The value of the next iteration. This represents the seventh auxiliary parameter. Indicates the first Shape parameters in the next iteration.

[0096] S2307. Calculation , .

[0097] S2308 Update :

[0098] (27)

[0099] S2309 Update .

[0100] In one specific embodiment, let Based on formulas (14) and (15), we can obtain formula (28):

[0101] (28)

[0102] In the formula, Indicates and Irrelevant constant terms.

[0103] From (28), we can see that It is approximated by a generalized inverse Gaussian distribution, as shown below:

[0104] (29)

[0105] (30)

[0106] In the formula, This represents the generalized inverse Gaussian distribution. This represents the shape parameter of the GIG distribution. The scale parameter representing the GIG distribution. This represents the location parameters of the GIG distribution. This indicates finding the trace of a matrix;

[0107] S2310. Calculation :

[0108] (31)

[0109] In the formula, express The second-order modified Bessel function of the second kind express The second-order modified Bessel function.

[0110] S2311. If the maximum number of iterations has not been reached, return to step S2302; otherwise, proceed to step S24.

[0111] S24. State Extraction: After the iteration converges, extract the state from the approximate distribution. Extract the mobile robot at the current moment The final posterior state estimation vector and the posterior error covariance matrix .

[0112] Specifically, , .

[0113] In some embodiments, the simulation experiment is based on a linear mobile robot system. The state transition matrix of the system model is set to 0.5, and the observation matrix is ​​set to 1. Initial state estimation. The initial covariance matrix is ​​100. The value is set to 10, the simulation duration is t = 2000 seconds, and the number of Monte Carlo simulations is 1000. In the experiment, the state noise is set to a normal distribution, and the observation noise is set to non-stationary, heavy-tailed noise. The corresponding noise covariance matrix values ​​are as follows: =100, =100. Parameter , , , These correspond to the probability of anomalies occurring within different time intervals. "wp" means "by probability".

[0114]

[0115]

[0116] This invention compares the root mean squared error (RMSE) and average root mean squared error (ARMSE) of its algorithm in state estimation with those of KFTCM, KFTPAC, RSTKF, GSTMKF, MLKF, and GaSTKF. Figure 2 As shown in Table 1.

[0117] Table 1 ARMSE of various algorithms

[0118]

[0119] Experimental results show that the proposed algorithm significantly outperforms the benchmark methods in terms of RMSE. Especially under time-varying noise covariance conditions, the proposed algorithm achieves the smallest RMSE fluctuation, demonstrating outstanding robustness. Quantitative results show that the proposed algorithm achieves a minimum RMSE of 8.076, significantly lower than methods such as RSTKF, GSTMKF, GaSTKF, and MLKF.

[0120] Simulation results show that RSTKF has the lowest estimation accuracy. This is because although the Student's t distribution can adaptively learn scale and degree-of-freedom parameters to model non-stationary heavy-tailed noise, it is limited by the choice of degree-of-freedom parameters and cannot adjust the mean vector or covariance matrix, thus its ability to model non-stationary measurement noise is weak. MLKF does not require the selection of degree-of-freedom parameters and outperforms RSTKF, but it still lacks correction for the mean vector and covariance matrix, resulting in poor modeling performance. In contrast, GSTMKF achieves better estimation performance by modeling non-stationary noise as a GSTM distribution, but it still suffers from the inherent limitations of the Student's t distribution in selecting degree-of-freedom parameters; furthermore, this method only uses the current measurement values ​​for mode selection and lacks dynamic inference of the noise covariance, which limits its effectiveness in anomaly detection. GaSTKF, on the other hand, can adjust the mean vector and covariance matrix, which further improves the estimation accuracy. However, as a variant of the Student's t-distribution, this method still requires the selection of a degree-of-freedom parameter; its adaptive estimation of this parameter introduces two degrees of freedom, increasing filtering uncertainty and thus limiting performance. In contrast, the proposed GaMLKF extends the standard machine learning model by introducing shape and scale parameters. This method overcomes the limitation of traditional machine learning distributions, which only adjust the scale matrix, by eliminating the need to select multiple degrees of freedom. GaMLKF can adaptively correct the mean vector and covariance matrix, thus more effectively modeling non-stationary measurement noise.

[0121] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise, characterized in that, Includes the following steps: S1. Establish a discrete-time linear state-space model of the mobile robot, including state equations and observation equations; S2. The state update of the mobile robot is achieved through a robust Kalman filter based on a gamma-multivariate Laplace hybrid distribution for heavy-tailed noise, including: S21. Prediction Step: Obtain the mobile robot's position in the previous time step. posterior state estimation vector and the posterior error covariance matrix The one-step prediction probability density function It approximates a Gaussian distribution, and from this, we can deduce the mobile robot's current state at the current moment. Prior state estimation vector and prior error covariance matrix ; S22. Bayesian Update: Get the Current Moment observation vector By combining the likelihood function constructed based on the gamma-multivariate Laplace mixture distribution, a joint posterior probability density function containing state and model noise parameters is formed. ; Represents the set of parameters to be estimated. Indicates the mobile robot at time The state vector, , Indicates auxiliary parameters, Indicates shape parameters; S23. Variational Inference: Using variational Bayesian methods, the joint posterior probability density function is... It can be decomposed into the product of independent factors of the state and model noise parameters, i.e. And approximate the true posterior through fixed-point iterative optimization; This represents the approximate posterior probability density function; S24. State Extraction: After the iteration converges, extract the state from the approximate distribution. Extract the mobile robot at the current moment The final posterior state estimation vector and the posterior error covariance matrix .

2. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 1, characterized in that, In step S1, Equations of state: , In the formula, Indicates the mobile robot at time The state vector, Indicates the mobile robot at time Given the state transition matrix, This represents process noise that conforms to a Gaussian distribution; Represents the set of real numbers. The dimension of the state vector; Observation equation: , In the formula, Indicates the mobile robot at time The observation vector, Indicates the mobile robot at time The known observation matrix; This indicates that the measurement noise follows a heavy-tailed characteristic with uncertainty. This represents the dimension of the observation vector.

3. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 1, characterized in that, Current moment One-step prediction probability density function Represented as: ， In the formula, Represents a multivariate normal distribution; Wherein, the prior state estimation vector and prior error covariance matrix Based on the state equations, it is expressed as: ， ， In the formula, Indicates the mobile robot at time Given the state transition matrix, This represents the process noise covariance matrix.

4. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 1, characterized in that, Based on the multivariate Laplace distribution, a gamma-multivariate Laplace mixture distribution is constructed by introducing a gamma distribution with adjustable shape and scale parameters. This gamma-multivariate Laplace mixture distribution is then used to model noise, and its probability density function is expressed as: ， In the formula, Represents a position vector. Represents the scale matrix, Indicates shape parameters, Indicates the first auxiliary parameter. Indicates the rate parameter. Indicates the DOF parameter. Represents a multivariate normal distribution; Indicates having a rate parameter The exponential distribution with a value of 1, Represents a random vector. Indicates the second auxiliary parameter. Represents a gamma distribution; where, a gamma distribution is used. First auxiliary parameter As a scaling factor, the scaling matrix is ​​realized. Adaptive calibration, and utilizing the first auxiliary parameter Shape parameters of multivariate Laplace distribution The product of position vectors Make compensation and corrections; Furthermore, the shape parameters Model it as a random variable and assign it a Gaussian prior distribution, defining... , This represents the mean hyperparameter. This represents the variance hyperparameter; Therefore, the gamma-multivariate Laplace mixture distribution can be represented by the following Gaussian hierarchical model: In the formula, Represents observed variables Given latent variables Multivariate normal distribution under certain conditions Representing latent variables The prior distribution, Representing latent variables The multivariate normal prior distribution, Representing latent variables The gamma prior distribution.

5. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 4, characterized in that, At any moment The likelihood function constructed based on the gamma-multivariate Laplace mixture distribution is expressed as: ， In the formula, Indicates the current time The likelihood probability density function, Indicates the mobile robot at time The observation vector, Indicates the mobile robot at time The known observation matrix; Indicates the mobile robot at time The state vector, Indicates time The measurement noise covariance matrix, This represents the DOF parameter.

6. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 1, characterized in that, Joint posterior probability density function Represented as: ， In the formula, Represents the time from time 1 to the current time. All observation data, This represents the time from moment 1 to the previous moment. All observation data, Represents a multivariate normal distribution; Indicates the mobile robot at time The observation vector, Indicates the mobile robot at time The known observation matrix; Indicates the mobile robot at time The state vector, Indicates having a rate parameter The exponential distribution with a value of 1, Represents a position vector. Indicates time The measurement noise covariance matrix, Indicates the gamma distribution; The vector representing the mean of a normal distribution is a representation of the mean of the distribution. The initial guess value, Represents constant real numbers. express An identity matrix of order 1.

7. The method for state estimation of a mobile robot based on Kalman filtering under thick-tailed noise according to claim 1, characterized in that, Step S23 involves initializing the initial values ​​for the iteration, and then performing Z approximate posterior probability density functions. The update process includes: S2301. Initialization: , , , The vector representing the mean of a normal distribution is a representation of the mean of the distribution. The initial guess value, Represents constant real numbers. express An identity matrix of order 1. Indicates the initial expected value; S2302. Calculation and : ， ， In the formula, This represents the corrected measurement noise covariance. This represents the adjusted mean. Indicates time Measurement noise covariance, Indicates the mobile robot at time The known observation matrix; Indicates the mobile robot at time The state vector; Indicates based on the first The mathematical expectation calculated from the approximate posterior distribution in the next iteration; S2303 Update : ， ， ， In the formula, Represents a multivariate normal distribution; Indicates the first The posterior state estimation vector under the next iteration. Indicates the mobile robot at time The observation vector, Indicates the first The posterior error covariance matrix under the next iteration express An identity matrix of order 1. Indicates the first Kalman gain in the next iteration; S2304 Update : ， ， ， ， In the formula, Indicates the first The first auxiliary parameter in the next iteration Indicates the DOF parameter. Represents point distribution. , , These represent the third, fourth, and fifth auxiliary parameters, respectively. Indicates the dimension of the observation vector; S2305. Calculation ; S2306 Update : ， ， ， ， In the formula, Indicates the sixth auxiliary parameter No. The value of the next iteration. express An identity matrix of order 1. Indicates the seventh auxiliary parameter No. The value of the next iteration. Represents constant real numbers. Indicates the first Shape parameters in the next iteration; S2307. Calculation , ; S2308 Update : ， S2309 Update : ， ， ， In the formula, This represents the generalized inverse Gaussian distribution. This represents the shape parameter of the GIG distribution. The scale parameter representing the GIG distribution. This represents the location parameters of the GIG distribution. This indicates finding the trace of a matrix; S2310. Calculation : ， In the formula, express The second-order modified Bessel function of the second kind express The second-order modified Bessel function; S2311. If the maximum number of iterations has not been reached, return to step S2302.

Images