A variational bayes-based carrier integrated navigation method

Abstract

The application relates to a carrier combined navigation method based on variational Bayes, and relates to the technical field of combined navigation.The application aims to solve the problem that the freedom degree parameter of the existing student t distribution is usually set as a fixed value, and the robustness and positioning accuracy of a carrier navigation system cannot be guaranteed when the carrier navigation system processes abnormal values.The process is as follows: S1: constructing a MEMS and GPS combined navigation system; letting k=1; S2: updating the state quantity at the k moment; S3: constructing a measurement noise student t distribution model at the k moment; updating the state quantity, the auxiliary variable and the freedom degree parameter in the model based on the updated state quantity at the k moment, and obtaining optimal estimation values of the state quantity, the auxiliary variable and the freedom degree parameter; S4: obtaining an optimal combined navigation system based on the optimal estimation values, and navigating the carrier based on the optimal combined navigation system; and S5: letting k=k+1, and repeatedly executing S2 to S4 until the combined navigation system of the carrier is closed.

Description

Technical Field

[0001] This invention relates to a carrier-based integrated navigation method based on variational Bayes, and relates to the field of integrated navigation technology. Background Technology

[0002] In integrated navigation systems utilizing both the Global Positioning System (GPS) and the Inertial Navigation System (INS), inertial auxiliary information can effectively improve GPS tracking performance and enhance the stability of the navigation system. Micro-Electro-Mechanical Systems (MEMS) inertial measurement units, due to their small size and low cost, can provide the inertial information required by aircraft in low-cost integrated navigation systems. However, due to manufacturing errors, in random vibration environments, low-cost sensors are simultaneously affected by internal noise and external vibration noise, resulting in unknown non-Gaussian noise in the output measurement noise, which significantly impacts the performance of the integrated navigation system. How to address the unknown non-Gaussian noise in low-cost integrated navigation systems to improve positioning accuracy is a problem urgently needing to be solved by those skilled in the art.

[0003] Currently, to estimate unknown noise in real time, the Variational Bayesian (VB) method is commonly used for adaptive estimation. However, this method typically assumes that the noise follows a Gaussian distribution, and the system performance will be significantly affected when outliers appear in the measurements. To improve the robustness of low-cost integrated navigation systems to non-Gaussian noise, Kalman filtering based on Huber, the maximum correlation entropy criterion, and the Student t-distribution is mainly used to handle situations where outliers appear in the system's measurements. However, the degrees of freedom parameters of the Student t-distribution are usually set to fixed values, which makes it impossible to guarantee the performance of robust filters in practical applications. Summary of the Invention

[0004] The purpose of this invention is to address the problem that the degree-of-freedom parameters of the existing student t-distribution are usually set to fixed values, which leads to the inability of the aircraft navigation system to guarantee robustness and positioning accuracy when handling outliers. Therefore, this invention proposes an aircraft integrated navigation method based on variational Bayes.

[0005] The specific process of a carrier-based integrated navigation method based on variational Bayes is as follows:

[0006] Step S1: Construct a MEMS and GPS integrated navigation system, which includes: state equations and measurement equations;

[0007] MEMS stands for Micro-Electro-Mechanical Systems;

[0008] GPS is the Global Positioning System.

[0009] The MEMS and GPS integrated navigation system is installed on the carrier aircraft;

[0010] Let k = 1;

[0011] Step S2: Update the state variables at time k;

[0012] Step S3: Construct the student t-distribution model of measurement noise at time k;

[0013] Based on the updated state variables at time k, the state variables X in the model are... k Auxiliary variable u k The degree-of-freedom parameter λk is updated to obtain the state variable X. k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimate;

[0014] Step S4: Based on state variable X k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimated value is used to obtain the optimal integrated navigation system, and the carrier aircraft uses the optimal integrated navigation system for navigation;

[0015] Step S5: Let k = k + 1, and repeat steps S2 to S4 until the aircraft's integrated navigation system is turned off.

[0016] Preferably, in step S1, a MEMS and GPS integrated navigation system is constructed, which includes: a state equation and a measurement equation.

[0017] MEMS stands for Micro-Electro-Mechanical Systems;

[0018] GPS is the Global Positioning System.

[0019] The MEMS and GPS integrated navigation system is installed on the carrier aircraft;

[0020] The specific process is as follows:

[0021] Step S11: In the launch inertial frame, select the aircraft's attitude error angle δα, position error Δp, velocity error Δv, and gyro constant drift error Δb. g and accelerometer constant drift error Δb a As a state variable X k-1 ,Right now

[0022] X k-1 =[δα,Δp,Δv,Δb] g ,Δba ] T (1)

[0023] in,

[0024] The superscript T denotes the transpose of the matrix, and k is the time k.

[0025] Based on state variable X k-1 The state equation for the integrated navigation system is as follows:

[0026]

[0027] In the formula,

[0028] Represents the state variable X k-1 The first derivative;

[0029] F represents the integrated navigation system matrix, and G represents the noise input matrix;

[0030] w is the noise of the combined navigation system with zero mean, and the variance of w is Q;

[0031] Step S12: Take the difference between the position and velocity information output by GPS and the position and velocity information output by MEMS as the measurement, i.e.

[0032]

[0033] In the formula,

[0034] Z k-1 Indicative measurement;

[0035] x GPS ,y GPS ,z GPS These represent the position coordinate components of the carrier aircraft in the x, y, and z directions in the launch inertial frame, as output by GPS.

[0036] These represent the velocity information components of the carrier aircraft in the x, y, and z directions in the launch inertial frame, respectively, output by GPS;

[0037] x INS ,y INS ,z INS These represent the position coordinate components of the carrier aircraft in the x, y, and z directions in the launch inertial frame output by the MEMS;

[0038] These represent the velocity information components of the carrier aircraft in the x, y, and z directions in the launch inertial frame output by the MEMS, respectively.

[0039] V GPS This refers to GPS measurement noise.

[0040] The measurement equations for the integrated navigation system are:

[0041] Z k-1 =H k-1 X k-1 +V GPS (4)

[0042] Among them, H k-1 This represents the measurement matrix.

[0043] Preferably, the integrated navigation system matrix F in step S11 is represented as:

[0044]

[0045] In the formula, Let be the direction cosine matrix between the carrier aircraft's own frame and the launch inertial frame. The specific force is obtained after correction by the accelerometer. Indicates calculation The antisymmetric matrix, F 32 Indicates intermediate variables;

[0046] The origin O of the carrier's OXYZ system is located at the center of mass of the carrier aircraft. The X-axis is the axis of symmetry of the carrier's outer shell, pointing towards the head as positive; the Y-axis is perpendicular to the X-axis, pointing upwards as positive; and the Z-axis is perpendicular to the XY plane.

[0047] I 3×3 Represents a 3×3 identity matrix, 0 3×3 Represents a 3×3 dimensional zero matrix;

[0048]

[0049] In the formula,

[0050] x, y, z are the position coordinate components of the carrier aircraft in the launch inertial frame, respectively;

[0051] R0 is the average Earth radius. Where GM is the distance from the aircraft to the Earth's center, and GM is the Earth's gravitational constant.

[0052] The noise input matrix G in step S11 is represented as follows:

[0053]

[0054] In step S11, the variance of w, Q, is expressed as:

[0055]

[0056] in, are the random drift noise variances of the gyroscope and accelerometer, respectively, and I represents the identity matrix;

[0057] In step S12, the measurement matrix H k-1 Represented as:

[0058]

[0059] Preferably, in step S2, the state quantity at time k is updated; the specific process is as follows:

[0060] Step S21: Based on the inertial information and state equation obtained from MEMS measurements, determine the state variable X at time k-1. k-1 Update the time;

[0061] The inertial information obtained from MEMS measurements includes and F 32 ;

[0062] The state equation is

[0063] The specific process is as follows:

[0064] Based on the state variable X at time k-1 k-1 Obtain the one-step prediction X of the state quantity at time k k|k-1 ; indicates as:

[0065] X k|k-1 =Φ k|k-1 X k-1 (10)

[0066] Where, Φ k|k-1 Here is the state transition matrix.

[0067] T′ represents the sampling time interval;

[0068] I 15×15 Represents a 15×15 dimensional identity matrix;

[0069] Step S22: Based on the state error covariance P at time k-1 k-1 One-step prediction P for calculating the state error covariance at time k k|k-1 ; indicates as:

[0070]

[0071] Among them, Γ k|k-1 For noise driving matrix,

[0072] T′ represents the sampling time interval; the superscript T indicates transpose;

[0073] P k-1 Let be the state error covariance at time k-1.

[0074] Preferably, in step S3, a student t-distribution model of measurement noise at time k is constructed;

[0075] Based on the updated state variables at time k, the state variables X in the model are... k Auxiliary variable u k Degrees of freedom parameter λ k Update and obtain the state variable X. k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimate;

[0076] The specific process is as follows:

[0077] Step S31: Construct the student t-distribution model of measurement noise at time k; initial iteration number j = 1;

[0078] Step S32: Obtain the set of parameters to be estimated Θ at time k. k A parameter in Bayesian updated Gaussian distribution

[0079] The set of parameters to be estimated Θ k Includes state variable X k Auxiliary variable u k and degree of freedom parameter λ k ;

[0080] Step S33: Based on the set of parameters to be estimated Θ k A parameter in Bayesian updated Gaussian distribution Update the state variables at time k and in the j-th iteration of the measurement noise student t-distribution model.

[0081] Step S34: Based on the set of parameters to be estimated Θ k A parameter in Bayesian updated Gaussian distribution Update the auxiliary variable in the k-th iteration of the student t-distribution model for measurement noise.

[0082] Step S35: Determine the degree of freedom parameters at time k and in the j-th iteration of the constructed measurement noise Student t-distribution model. Adaptive adjustment;

[0083] Step S36: Let the iteration number j = j + 1, and repeat steps S33 to S35 until the parameters to be estimated in the measurement noise student t-distribution model remain unchanged, thus obtaining the state variable X at time k. k Auxiliary variable uk and degree of freedom parameter λ k The optimal estimate;

[0084] The parameter to be estimated is state X. k Auxiliary variable u k and degree of freedom parameter λ k .

[0085] Preferably, in step S31, a student t-distribution model of measurement noise at time k is constructed; the specific process is as follows:

[0086] The probability density function of a random variable x that follows a Student's t-distribution is:

[0087]

[0088] in,

[0089] St(x|μ,Δ,λ k Let f(x) represent the probability density function of a random variable x that follows a Student t-distribution;

[0090] St() represents the probability density function;

[0091] d is the dimension of the random variable x;

[0092] μ,Δ,λ k Let represent the mean, variance, and degrees of freedom parameters of the student t-distribution, respectively;

[0093] Γ(·) represents the Gamma function;

[0094] Introduce auxiliary variable u k Rewrite equation (12) as follows:

[0095]

[0096] in,

[0097] N(·) and C(·) represent the Gaussian distribution and the Gamma distribution, respectively;

[0098] N(x|μ,(u k Δ) -1 () represents the Gaussian distribution of the random variable x;

[0099] Represents the auxiliary variable u k The Gamma distribution;

[0100] Then the likelihood function of the integrated navigation system at time k is expressed as:

[0101]

[0102] In the formula,

[0103] R k Let be the measurement noise variance at time k.

[0104] Z represents a quantity that follows a student t-distribution. k The probability density function;

[0105] p(Z k |X k ) represents the likelihood function;

[0106] Z k Indicative measurement;

[0107] According to equation (13), the likelihood function of the integrated navigation system at time k is rewritten as follows:

[0108]

[0109] in,

[0110] Representing variable Z k Gaussian distribution;

[0111] p(u k |λ k ) represents the auxiliary variable u k The probability density function.

[0112] Preferably, in step S32, the set of parameters to be estimated Θ at time k is obtained. k A parameter in Bayesian updated Gaussian distribution The specific process is as follows:

[0113] The state X is achieved by minimizing the following index function. k Auxiliary variable u k and degree of freedom parameter λ k The estimate; the expression is:

[0114] J=KL(q(Θ k )||p(Θ k |Z 1:k (17)

[0115] In the formula,

[0116] J represents the index function;

[0117] KL(·) represents the calculation of KL divergence;

[0118] q(Θ k ) represents the set of parameters to be estimated Θ k Approximation of the distribution, the set of parameters to be estimated Θk State X k Auxiliary variable u k and degree of freedom parameter λ k ;

[0119] p(Θ k |Z 1:k ) represents the set of parameters to be estimated Θ k The posterior distribution of Z; 1:k Represents the set of all measurements from time 1 to time k;

[0120] The solution to equation (17) is

[0121]

[0122] In the formula,

[0123] express The Gaussian distribution updated by Bayes;

[0124] E[·] represents the expectation operation. The set of parameters to be estimated Θ k One of the parameters in, Represents the set of parameters to be estimated, Θ k Except Other parameters;

[0125] p(Θ k Z 1:k ) represents the parameter set Θ k and measurement set Z 1:k The joint probability distribution of ;

[0126] To and Relevant constants.

[0127] Preferably, in step S33, the set of parameters to be estimated Θ is used. k A parameter in Bayesian updated Gaussian distribution Update the state variables at time k and in the j-th iteration of the measurement noise student t-distribution model. The specific process is as follows:

[0128] make Substituting into equation (18), we obtain the state. Bayesian updated Gaussian distribution

[0129] Gaussian distribution The mean is (19), Gaussian distribution. The variance is (20);

[0130]

[0131]

[0132] in,

[0133] This represents the state quantity at time k and in the j-th iteration of the student t-distribution model for measurement noise;

[0134] This represents the error covariance estimate for the j-th iteration at time k;

[0135] X k|k-1 This represents the one-step prediction of the state at time k.

[0136] P k|k-1 This represents a one-step prediction of the state error covariance at time k.

[0137] H k This represents the measurement matrix at time k.

[0138] This represents the Kalman filter gain matrix at time k in the j-th iteration;

[0139] δZ k The information at time k is calculated as follows:

[0140]

[0141] In the formula,

[0142] This represents the actual measurement value at time k.

[0143] The superscript T indicates transpose;

[0144] Let the new information covariance matrix be the matrix of the j-th iteration at time k.

[0145]

[0146] In the formula,

[0147] This represents the auxiliary variable for the (j-1)th iteration at time k.

[0148] R k This represents the measurement noise variance at time k.

[0149] Preferably, in step S34, the set of parameters to be estimated Θ is used as the basis. k A parameter in Bayesian updated Gaussian distribution Update the auxiliary variable in the k-th iteration of the student t-distribution model for measurement noise. The specific process is as follows:

[0150] make Substituting into equation (18), we get Bayesian updated Gamma distribution

[0151] Gamma distribution The update formula is (24);

[0152]

[0153] in,

[0154] This represents the degree of freedom parameter at time k, during the (j-1)th iteration;

[0155] d represents Z k The dimension;

[0156] This represents the auxiliary variable for the j-th iteration at time k.

[0157] γ k Indicates intermediate variables;

[0158]

[0159] In the formula,

[0160] trace is the matrix trace operation;

[0161] This represents the error covariance estimate for the (j-1)th iteration at time k;

[0162] R k This represents the measurement noise variance at time k.

[0163] Preferably, in step S35, the degree of freedom parameters at time k and the j-th iteration in the constructed measurement noise student t-distribution model are... Adaptive adjustment; the specific process is as follows:

[0164] The information δZ of ​​the j-th iteration at time k in the integrated navigation system k Mahalanobis distance for:

[0165]

[0166] In the formula,

[0167] Let the covariance matrix be the new information of the (j-1)th iteration at time k.

[0168] Based on the new information δZ of ​​the j-th iteration at time k k Mahalanobis distance Adjusting the degree of freedom parameter at time k in the student t-distribution model Represented as:

[0169]

[0170] In the formula, λ0 is the prior degrees of freedom parameter of the student t distribution at the initial time.

[0171] The beneficial effects of this invention are as follows:

[0172] This invention proposes a carrier-based integrated navigation method based on variational Bayesian methods. It models non-Gaussian measurement noise as a Student t-distribution and uses the variational Bayesian method to estimate the unknown statistical characteristics of the noise. At the same time, it uses Mahalanobis distance to dynamically adjust the degree of freedom parameters of the Student t-distribution, thereby improving the robustness and positioning accuracy of the low-cost integrated navigation method when subjected to unknown non-Gaussian noise interference. Attached Figure Description

[0173] Figure 1 This is a flowchart of the present invention;

[0174] Figure 2 This is a flight trajectory diagram of the carrier aircraft;

[0175] Figure 3 A location settlement error diagram;

[0176] Figure 4 This is a speed calculation error diagram. Detailed Implementation

[0177] Specific Implementation Method 1: The specific process of this implementation method for aircraft integrated navigation based on variational Bayes is as follows:

[0178] Step S1: Construct a MEMS and GPS integrated navigation system, which includes: state equations and measurement equations;

[0179] MEMS stands for Micro-Electro-Mechanical Systems;

[0180] GPS is the Global Positioning System.

[0181] The MEMS and GPS integrated navigation system is installed on the carrier aircraft;

[0182] Let k = 1;

[0183] Step S2: Update the state variables at time k;

[0184] Step S3: Construct the student t-distribution model of measurement noise at time k;

[0185] Based on the updated state variables at time k, the state variables X in the model are... k Auxiliary variable u k Degrees of freedom parameter λ k Update and obtain the state variable X. k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimate;

[0186] Step S4: Based on state variable X k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimated value is used to obtain the optimal integrated navigation system, and the carrier aircraft uses the optimal integrated navigation system for navigation;

[0187] Step S5: Let k = k + 1, and repeat steps S2 to S4 until the aircraft's integrated navigation system is turned off.

[0188] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that: in step S1, a MEMS and GPS integrated navigation system is constructed. The MEMS and GPS integrated navigation system includes: state equation (Formula 2) and measurement equation (Formula 4).

[0189] MEMS stands for Micro-Electro-Mechanical Systems;

[0190] GPS is the Global Positioning System.

[0191] The MEMS and GPS integrated navigation system is installed on the carrier aircraft;

[0192] The specific process is as follows:

[0193] Step S11: In the launch inertial frame, select the aircraft's attitude error angle δα, position error Δp, velocity error Δv, and gyro constant drift error Δb. g and accelerometer constant drift error Δb a As a state variable X k-1 ,Right now

[0194] X k-1 =[δα,Δp,Δv,Δb] g ,Δb a ] T (1)

[0195] in,

[0196] The superscript T denotes the transpose of the matrix, and k is the time k.

[0197] The launch inertial frame is the same as the launch coordinate system at the moment of takeoff, and then the launch inertial frame is stationary relative to inertial space.

[0198] Attitude error angle δα includes heading angle error, pitch angle error, and roll angle error;

[0199] Based on state variable X k-1 The state equation for the integrated navigation system is as follows:

[0200]

[0201] In the formula,

[0202] Represents the state variable X k-1 The first derivative;

[0203] F represents the integrated navigation system matrix, and G represents the noise input matrix;

[0204] w is the noise of the combined navigation system with zero mean, and the variance of w is Q;

[0205] Step S12: Take the difference between the position and velocity information output by GPS and the position and velocity information output by MEMS as the measurement, i.e.

[0206]

[0207] In the formula,

[0208] Z k-1 Indicative measurement;

[0209] x GPS ,y GPS ,z GPS These represent the position coordinate components of the carrier aircraft in the x, y, and z directions in the launch inertial frame, as output by GPS.

[0210] These represent the velocity information components of the carrier aircraft in the x, y, and z directions in the launch inertial frame, respectively, output by GPS;

[0211] x INS ,y INS ,z INS These represent the position coordinate components of the carrier aircraft in the x, y, and z directions in the launch inertial frame output by the MEMS;

[0212] These represent the velocity information components of the carrier aircraft in the x, y, and z directions in the launch inertial frame output by the MEMS, respectively.

[0213] V GPS This refers to GPS measurement noise.

[0214] The measurement equations for the integrated navigation system are:

[0215] Z k-1 =H k-1 X k-1 +V GPS (4)

[0216] Among them, H k-1 This represents the measurement matrix.

[0217] The other steps and parameters are the same as in Specific Implementation Method 1.

[0218] Specific Implementation Method Three: This implementation method differs from Specific Implementation Method One or Two in that the integrated navigation system matrix F in step S11 is represented as follows:

[0219]

[0220] In the formula, Let be the direction cosine matrix between the carrier aircraft's own frame and the launch inertial frame. The specific force is obtained after correction by the accelerometer. Indicates calculation The antisymmetric matrix, F 32 Indicates intermediate variables;

[0221] The origin O of the carrier system OXYZ is located at the center of mass of the carrier aircraft. The X-axis is the axis of symmetry of the carrier aircraft shell (the direction of the line connecting the carrier aircraft nose and fuselage is the X-axis), pointing towards the nose is positive; the Y-axis is perpendicular to the X-axis, pointing upwards is positive; the Z-axis is perpendicular to the XY plane.

[0222] I 3×3 Represents a 3×3 identity matrix, 0 3×3 Represents a 3×3 dimensional zero matrix;

[0223]

[0224] In the formula,

[0225] x, y, z are the position coordinate components of the carrier aircraft in the launch inertial frame, respectively;

[0226] R0 is the average Earth radius. Where GM is the distance from the aircraft to the Earth's center, and GM is the Earth's gravitational constant.

[0227] The noise input matrix G in step S11 is represented as follows:

[0228]

[0229] In step S11, the variance of w, Q, is expressed as:

[0230]

[0231] in, are the random drift noise variances of the gyroscope and accelerometer, respectively, and I represents the identity matrix;

[0232] In step S12, the measurement matrix H k-1 Represented as:

[0233]

[0234] Other steps and parameters are the same as in specific implementation method one or two.

[0235] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that: the state quantity at time k is updated in step S2; the specific process is as follows:

[0236] Step S21: Based on the inertial information and state equation obtained from MEMS measurements, determine the state variable X at time k-1. k-1 Update the time;

[0237] The inertial information obtained from MEMS measurements includes and F 32 ;

[0238] The state equation is

[0239] The specific process is as follows:

[0240] Based on the state variable X at time k-1 k-1 Obtain the one-step prediction X of the state quantity at time k k|k-1 ; indicates as:

[0241] X k|k-1 =Φ k|k-1 X k-1 (10)

[0242] Where, Φ k|k-1 Here is the state transition matrix.

[0243] T′ represents the sampling time interval;

[0244] I 15×15 Represents a 15×15 dimensional identity matrix;

[0245] Step S22: Based on the state error covariance P at time k-1 k-1 One-step prediction P for calculating the state error covariance at time k k|k-1 ; indicates as:

[0246]

[0247] Among them, Γ k|k-1 For noise driving matrix,

[0248] T′ represents the sampling time interval; the superscript T indicates transpose;

[0249] P k-1 Let be the state error covariance at time k-1.

[0250] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0251] Specific Implementation Method Five: This implementation method differs from one of the specific implementation methods one to four in that: in step S3, a student t-distribution model of measurement noise at time k is constructed;

[0252] Based on the updated state variables at time k, the state variables X in the model are... k Auxiliary variable u k Degrees of freedom parameter λ k Update and obtain the state variable X. k Auxiliary variable u k Degrees of freedom parameter λ k The optimal estimate;

[0253] The specific process is as follows:

[0254] To address the issue of unknown measurement noise and non-Gaussian statistical characteristics in low-cost integrated navigation systems, this invention models the non-Gaussian measurement noise as a Student's t-distribution and uses the variational Bayesian method to jointly estimate the unknown noise and the system state.

[0255] Step S31: Construct the student t-distribution model of measurement noise at time k; initial iteration number j = 1;

[0256] Step S32: Obtain the set of parameters to be estimated Θ at time k. k A parameter in Bayesian updated Gaussian distribution

[0257] The set of parameters to be estimated Θ k Includes state variable X k Auxiliary variable u k and degree of freedom parameter λ k ;

[0258] Step S33: Based on the set of parameters to be estimated Θ k A parameter in Bayesian updated Gaussian distribution Update the state variables at time k and in the j-th iteration of the measurement noise student t-distribution model.

[0259] Step S34: Based on the set of parameters to be estimated Θ kA parameter in Bayesian updated Gaussian distribution Update the auxiliary variable in the k-th iteration of the student t-distribution model for measurement noise.

[0260] Step S35: Determine the degree of freedom parameters at time k and in the j-th iteration of the constructed measurement noise Student t-distribution model. Adaptive adjustment;

[0261] Step S36: The variational Bayesian method iterates to make the estimated parameters gradually approach the true values.

[0262] Let the iteration number j = j + 1, and repeat iteration steps S33 to S35 until the parameters to be estimated in the measurement noise student t-distribution model remain unchanged, thus obtaining the state variable X at time k. k Auxiliary variable u k and degree of freedom parameter λ k The optimal estimate;

[0263] The parameter to be estimated is state X. k Auxiliary variable u k and degree of freedom parameter λ k ;

[0264] The noise mean of the t-distribution is usually assumed to be 0 and does not need to be estimated. The variance R is prior information, and the auxiliary variable u is used. k For (1 / u) k To make adaptive adjustments, the parameters that need to be estimated are the state, auxiliary variables, and degrees of freedom parameters.

[0265] The other steps and parameters are the same as those in specific implementation methods one through four.

[0266] Specific Implementation Method Six: This implementation method differs from Specific Implementation Methods One to Five in that: in step S31, the student t-distribution model of measurement noise at time k is constructed (Formula 13); the specific process is as follows:

[0267] This invention introduces the Student t-distribution to model non-Gaussian measurement noise;

[0268] The probability density function of a random variable x that follows a Student's t-distribution is:

[0269]

[0270] in,

[0271] St(x|μ,Δ,λ k Let f(x) represent the probability density function of a random variable x that follows a Student t-distribution;

[0272] St() represents the probability density function;

[0273] d is the dimension of the random variable x;

[0274] μ,Δ,λ k Let represent the mean, variance, and degrees of freedom parameters of the student t-distribution, respectively;

[0275] Γ(·) represents the Gamma function;

[0276] To obtain a closed-form solution to the student t-distribution, this invention introduces an auxiliary variable u. k Rewrite equation (12) as follows:

[0277]

[0278] in,

[0279] N(·) and C(·) represent the Gaussian distribution and the Gamma distribution, respectively;

[0280] N(x|μ,(u k Δ) -1 () represents the Gaussian distribution of the random variable x;

[0281] Represents the auxiliary variable u k The Gamma distribution;

[0282] Then the likelihood function of the integrated navigation system at time k is expressed as:

[0283]

[0284] In the formula,

[0285] R k Let be the measurement noise variance at time k.

[0286] Z represents a quantity that follows a student t-distribution. k The probability density function;

[0287] p(Z k |X k ) represents the likelihood function;

[0288] Z k Indicative measurement;

[0289] According to equation (13), the likelihood function of the integrated navigation system at time k is rewritten as follows:

[0290]

[0291] in,

[0292] Representing variable Z k Gaussian distribution;

[0293] p(u k |λ k ) represents the auxiliary variable u k The probability density function.

[0294] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0295] Specific Implementation Method Seven: This implementation method differs from Specific Implementation Methods One through Six in that: in step S32, the set of parameters to be estimated at time k is obtained. k A parameter in Bayesian updated Gaussian distribution The specific process is as follows:

[0296] To address the problem of estimating noise parameters and state, this invention employs a variational Bayesian method for joint estimation. The variational Bayesian method achieves the estimation of state X by minimizing the following index function. k Auxiliary variable u k and degree of freedom parameter λ k The estimate; the expression is:

[0297] J=KL(q(Θ k )||p(Θ k |Z 1:k (17)

[0298] In the formula,

[0299] J represents the index function;

[0300] KL(·) represents the calculation of KL divergence;

[0301] q(Θ k ) represents the set of parameters to be estimated Θ k Approximation of the distribution, the set of parameters to be estimated Θ k State X k Auxiliary variable u k and degree of freedom parameter λ k ;

[0302] p(Θ k |Z 1:k ) represents the set of parameters to be estimated Θ k The posterior distribution of Z; 1:k Represents the set of all measurements from time 1 to time k;

[0303] The solution to equation (17) is

[0304]

[0305] In the formula,

[0306] express The Gaussian distribution updated by Bayes;

[0307] E[·] represents the expectation operation. The set of parameters to be estimated Θ k One of the parameters in, Represents the set of parameters to be estimated, Θ k Except Other parameters;

[0308] p(Θ k Z 1:k ) represents the parameter set Θ k and measurement set Z 1:k The joint probability distribution of ;

[0309] To and Relevant constants.

[0310] The other steps and parameters are the same as those in specific implementation methods one through six.

[0311] Specific Implementation Method Eight: This implementation method differs from Specific Implementation Methods One to Seven in that: in step S33, the set of parameters to be estimated Θ is used... k A parameter in Bayesian updated Gaussian distribution Update the state variables at time k and in the j-th iteration of the measurement noise student t-distribution model. The specific process is as follows:

[0312] make Substituting into equation (18), we obtain the state. Bayesian updated Gaussian distribution

[0313] Gaussian distribution The mean is (19), Gaussian distribution. The variance is (20);

[0314]

[0315] in,

[0316] This represents the state quantity at time k and in the j-th iteration of the student t-distribution model for measurement noise;

[0317] This represents the error covariance estimate for the j-th iteration at time k;

[0318] Xk|k-1 This represents the one-step prediction of the state at time k.

[0319] P k|k-1 This represents a one-step prediction of the state error covariance at time k.

[0320] H k This represents the measurement matrix at time k.

[0321] This represents the Kalman filter gain matrix at time k in the j-th iteration;

[0322] δZ k The information at time k is calculated as follows:

[0323]

[0324] In the formula,

[0325] This represents the actual measurement value at time k.

[0326] The superscript T indicates transpose;

[0327] Let the new information covariance matrix be the matrix of the j-th iteration at time k.

[0328]

[0329] In the formula,

[0330] This represents the auxiliary variable for the (j-1)th iteration at time k.

[0331] R k This represents the measurement noise variance at time k.

[0332] The other steps and parameters are the same as those in specific implementation methods one through seven.

[0333] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One to Eight in that: in step S34, the set of parameters to be estimated Θ is used... k A parameter in Bayesian updated Gaussian distribution Update the auxiliary variable in the k-th iteration of the student t-distribution model for measurement noise. The specific process is as follows:

[0334] make Substituting into equation (18), we get Bayesian updated Gamma distribution

[0335] Gamma distribution The update formula is (24);

[0336]

[0337] in,

[0338] This represents the degree of freedom parameter at time k, during the (j-1)th iteration;

[0339] d represents Z k The dimension;

[0340] This represents the auxiliary variable for the j-th iteration at time k.

[0341] γk represents an intermediate variable;

[0342]

[0343] In the formula,

[0344] trace is the matrix trace operation;

[0345] This represents the error covariance estimate for the (j-1)th iteration at time k;

[0346] R k This represents the measurement noise variance at time k.

[0347] The other steps and parameters are the same as those in specific implementation methods one through eight.

[0348] Specific Implementation Method Ten: This implementation method differs from Specific Implementation Methods One through Nine in that: in step S35, the degree of freedom parameters of the j-th iteration at time k in the constructed measurement noise student t-distribution model are... Adaptive adjustment; the specific process is as follows:

[0349] Traditional methods typically model the Student t-distribution using fixed degrees of freedom parameters. However, when the prior degrees of freedom parameters do not match the actual situation, the estimation accuracy of the algorithm decreases. To address this issue, this invention adaptively adjusts the degrees of freedom parameters in the Student t-distribution using Mahalanobis distance. Mahalanobis distance is commonly used to describe the distance between a point and a distribution. Considering the innovation δZ of ​​the j-th iteration at time k in a low-cost integrated navigation system... k Mahalanobis distance for:

[0350]

[0351] In the formula,

[0352] Let the covariance matrix be the new information of the (j-1)th iteration at time k.

[0353] As can be seen, the larger the Mahalanobis distance at the current time, the greater the difference between the measured value and the Gaussian distribution, and therefore the greater the probability of measurement outliers. The corresponding degree-of-freedom parameter of the Student's t-distribution should be adjusted to be smaller. Then, based on the information δZ from the j-th iteration at time k... k Mahalanobis distance Adjusting the degree of freedom parameter at time k in the student t-distribution model Represented as:

[0354]

[0355] In the formula, λ0 is the prior degrees of freedom parameter of the student t distribution at the initial time.

[0356] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0357] The beneficial effects of the present invention are verified using the following embodiments:

[0358] Example 1:

[0359] To verify the effectiveness of the proposed method, a navigation scenario was first constructed for simulation. It was assumed that the vehicle was initially positioned at longitude 120°, latitude 70°, and altitude 2km. The accelerometer constant drift was 1mg, and the random drift was 25mg. The gyroscope constant drift was 20° / h, and the random drift was... The simulation duration is 100 seconds, with a step size of 0.1 seconds. The simulation scenario is as follows: Figure 2 .

[0360] The GPS navigation position error is 10m and the velocity error is 2.5m / s. However, its measurement is affected by the environment, which causes outliers in the GPS navigation error. Tracking simulations were performed using the extended Kalman filter (EKF), the variational Bayes-based Kalman filter (VBEKF), and the algorithm proposed in this invention.

[0361] Tracking results for example Figure 3 , Figure 4It can be seen that the traditional extended Kalman filter cannot handle situations where the statistical characteristics of the measurement are unknown, and therefore the navigation estimation accuracy is the worst when outliers appear in the measurement. While the traditional variational Bayesian Kalman filter adaptively estimates the unknown measurement noise, it assumes that it follows a Gaussian distribution, so its accuracy is still affected when the measurement noise exhibits non-Gaussian characteristics. The algorithm proposed in this invention, by modeling the unknown non-Gaussian noise as a Student's t-distribution and using the variational Bayesian method to adaptively estimate the unknown statistical characteristic parameters, and by employing the Mahalanobis distance enhancement system to adapt to the unknown parameters, achieves the highest navigation estimation accuracy. The method proposed in this invention effectively improves the robustness and positioning accuracy of low-cost integrated navigation methods when subjected to unknown non-Gaussian noise interference.

[0362] This invention may have other embodiments. Without departing from the spirit and essence of this invention, those skilled in the art can make various corresponding changes and modifications according to this invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A carrier-based integrated navigation method based on variational Bayes, characterized in that: The specific process of the method is as follows: Step S1: Construct a MEMS and GPS integrated navigation system, which includes: state equations and measurement equations; MEMS stands for Micro-Electro-Mechanical Systems; GPS is the Global Positioning System. The MEMS and GPS integrated navigation system is installed on the carrier aircraft; make ; Step S2: Update State quantity at any given moment; Step S3: Build Time measurement noise for students Distribution model; Based on the update The state variables at time step versus the state variables in the model Auxiliary variables Degrees of freedom parameters Update and obtain the state. Auxiliary variables Degrees of freedom parameters The optimal estimate; Step S4: Based on state variables Auxiliary variables Degrees of freedom parameters The optimal estimated value is used to obtain the optimal integrated navigation system, and the carrier aircraft uses the optimal integrated navigation system for navigation; Step S5: Let Repeat steps S2 to S4 until the aircraft's integrated navigation system is turned off; The adaptive adjustment process of the degree-of-freedom parameters in step S3 is as follows: Integrated navigation system Time of the first The next iteration's new information Mahalanobis distance for: (26) In the formula, For the present Time of the first The new information covariance matrix of the next iteration; the superscript T denotes the transpose of the matrix; according to Time of the first The next iteration's new information Mahalanobis distance Adjusting students In the distribution model Degrees of freedom parameters at time ; indicates as: (27) In the formula, For students at the initial moment The prior degrees of freedom parameters of the distribution.

2. The aircraft integrated navigation method based on variational Bayes as described in claim 1, characterized in that: In step S1, a MEMS and GPS integrated navigation system is constructed. The MEMS and GPS integrated navigation system includes: state equation and measurement equation. MEMS stands for Micro-Electro-Mechanical Systems; GPS is the Global Positioning System. The MEMS and GPS integrated navigation system is installed on the carrier aircraft; The specific process is as follows: Step S11: Select the attitude error angle of the carrier aircraft in the launch inertial frame. Position error Speed ​​error Gyroscope constant drift error and accelerometer constant drift error As a state quantity ,Right now (1) in, The superscript T denotes the transpose of a matrix. for time; Based on state variables The state equation for the integrated navigation system is as follows: (2) In the formula, State variables The first derivative; Represents the matrix of the integrated navigation system. Represents the noise input matrix; The noise of the integrated navigation system is zero-mean. The variance is ; Step S12: Take the difference between the position and velocity information output by GPS and the position and velocity information output by MEMS as the measurement, i.e. (3) In the formula, Indicative measurement; These represent the GPS output of the carrier aircraft in the launching inertial frame, respectively. Position coordinate components in three directions; These represent the GPS output of the carrier aircraft in the launching inertial frame, respectively. Velocity information components in three directions; These represent the carrier aircraft output by the MEMS in the launch inertial frame, respectively. Position coordinate components in three directions; These represent the carrier aircraft output by the MEMS in the launch inertial frame, respectively. Velocity information components in three directions; This refers to GPS measurement noise. The measurement equations for the integrated navigation system are: (4) in, This represents the measurement matrix.

3. The aircraft integrated navigation method based on variational Bayes as described in claim 2, characterized in that: The integrated navigation system matrix in step S11 Represented as: (5) In the formula, Let be the direction cosine matrix between the carrier aircraft's own frame and the launch inertial frame. The specific force is obtained after correction by the accelerometer. Indicates calculation antisymmetric matrix, Indicates intermediate variables; Carrier System The origin Located at the center of gravity of the carrier aircraft, The axis is the axis of symmetry of the aircraft's outer shell, with the direction pointing towards the head being positive; Axis perpendicular The axis, pointing upwards, is positive; Axis perpendicular flat; express An identity matrix of dimension 1 express A zero-dimensional matrix; (6) In the formula, These are the position coordinate components of the carrier aircraft in the launch inertial frame; The average Earth radius, The distance from the aircraft to the Earth's center. The gravitational constant of Earth; The noise input matrix in step S11 Represented as: (7) In step S11 The variance is Represented as: (8) in, These are the variances of random drift noise for the gyroscope and accelerometer, respectively. Represents the identity matrix; The measurement matrix in step S12 Represented as: (9)。 4. The aircraft integrated navigation method based on variational Bayes as described in claim 3, characterized in that: Update in step S2 The state quantity at time t; the specific process is as follows: Step S21: Based on the inertial information and state equations obtained from MEMS measurements, perform... state quantity at time 1 Update the time; The inertial information obtained from MEMS measurements includes , , and ; The state equation is ; The specific process is as follows: based on state quantity at time 1 get One-step prediction of state quantities at time step ; indicates as: (10) in, Here is the state transition matrix. ; Indicates the sampling time interval; express An identity matrix of dimensionality; Step S22: Based on State error covariance at time 1 calculate One-step prediction of state error covariance at time step ; indicates as: (11) in, For noise driving matrix, ; Indicates the sampling time interval; the superscript T indicates transpose. for The state error covariance at time t.

5. The aircraft integrated navigation method based on variational Bayes as described in claim 4, characterized in that: The construction in step S3 Time measurement noise for students Distribution model; Based on the update The state variables at time step versus the state variables in the model Auxiliary variables Degrees of freedom parameters Update and obtain the state. Auxiliary variables Degrees of freedom parameters The optimal estimate; The specific process is as follows: Step S31: Build Time measurement noise for students Distribution model; initial number of iterations ; Step S32: Obtain The set of parameters to be estimated at time 1 A parameter in Bayesian updated Gaussian distribution ; Set of parameters to be estimated Includes state variables Auxiliary variables and degrees of freedom parameters ; Step S33: Based on the set of parameters to be estimated A parameter in Bayesian updated Gaussian distribution Update student noise measurement In the distribution model Time of the first The state variables of the next iteration ; Step S34: Based on the set of parameters to be estimated A parameter in Bayesian updated Gaussian distribution Update student noise measurement In the distribution model Time of the first Auxiliary variables for the next iteration ; Step S35: For the constructed measurement noise student In the distribution model Time of the first The degree of freedom parameters of the next iteration Adaptive adjustment; Step S36: Let the number of iterations be... Repeat steps S33 to S35 until the measurement noise is reached. In the distribution model, the parameters to be estimated remain unchanged, and the results are obtained. Time-state quantity Auxiliary variables and degrees of freedom parameters The optimal estimate; The parameter to be estimated is the state. Auxiliary variables and degrees of freedom parameters .

6. The aircraft integrated navigation method based on variational Bayes as described in claim 5, characterized in that: The construction in step S31 Time measurement noise for students Distribution model; the specific process is as follows: Obeying students Distributed random variables The probability density function is: (12) in, Showing obedience to students Distributed random variables The probability density function; Represents the probability density function; For random variables dimensionality; They represent students respectively. The mean, variance, and degrees of freedom parameters of the distribution; Represents the Gamma function; Introducing auxiliary variables Rewrite equation (12) as follows: (13) in, and Represent the Gaussian distribution and the Gamma distribution, respectively; Represents random variables Gaussian distribution; Representing auxiliary variables The Gamma distribution; but The likelihood function of the combined navigation system at time t is expressed as: (14) In the formula, For the present Measurement noise variance at any given time; Showing obedience to students Distribution measurement The probability density function; Represents the likelihood function; Indicative measurement; Based on equation (13), rewrite The likelihood function of the combined navigation system at time t is (15) (16) in, Representing variables Gaussian distribution; Representing auxiliary variables The probability density function.

7. The aircraft integrated navigation method based on variational Bayes as described in claim 6, characterized in that: In step S32, the following is obtained The set of parameters to be estimated at time 1 A parameter in Bayesian updated Gaussian distribution The specific process is as follows: The state is achieved by minimizing the following index function. Auxiliary variables and degrees of freedom parameters The estimate; the expression is: (17) In the formula, Indicates the index function; This indicates the calculation of KL divergence; Represents the set of parameters to be estimated Approximation of the distribution, set of parameters to be estimated For state Auxiliary variables and degrees of freedom parameters ; The set of parameters to be estimated The posterior distribution of; Indicates time 1 to The set of all measurements at any given time; The solution to equation (17) is (18) In the formula, express The Gaussian distribution updated by Bayes; This indicates the operation of finding the expected value. The set of parameters to be estimated One of the parameters in, Represents the set of parameters to be estimated Except Other parameters; Represents the parameter set and measurement set The joint probability distribution of ; To and Relevant constants.

8. The aircraft integrated navigation method based on variational Bayes as described in claim 7, characterized in that: In step S33, the set of parameters to be estimated is used. A parameter in Bayesian updated Gaussian distribution Update student noise measurement In the distribution model Time of the first The state variables of the next iteration The specific process is as follows: make Substituting into equation (18), we obtain the state. Bayesian updated Gaussian distribution ; Gaussian distribution The mean is (19), Gaussian distribution. The variance is (20); (19) (20) in, Indicates the student's measurement noise In the distribution model Time of the first The state variables of the next iteration; express Time of the first The error covariance estimate for the next iteration; Indicates the current A one-step prediction of the state at a given time; Indicates the current One-step prediction of the state error covariance at time step; Indicates the current The measurement matrix at time; Indicates the current Time of the first The Kalman filter gain matrix for the next iteration; For the present The information at any given moment is calculated as follows: (21) (22) In the formula, Indicates the current The actual measured value at that moment; The superscript T indicates transpose; For the present Time of the first The new covariance matrix of the next iteration; (23) In the formula, Indicates the current Time of the first Auxiliary variables for the next iteration; Indicates the current Measurement noise variance at any given time.

9. The aircraft integrated navigation method based on variational Bayes as described in claim 8, characterized in that: In step S34, the set of parameters to be estimated is used. A parameter in Bayesian updated Gaussian distribution Update student noise measurement In the distribution model Time of the first Auxiliary variables for the next iteration The specific process is as follows: make Substituting into equation (18), we get Bayesian updated Gamma distribution ; Gamma distribution The update formula is (24); (24) in, Indicates the current Time of the first The degree of freedom parameters for the next iteration; express The dimension; Indicates the current Time of the first Auxiliary variables for the next iteration; Indicates intermediate variables; (25) In the formula, This is a matrix trace operation; express Time of the first The error covariance estimate for the next iteration; express Measurement noise variance at any given time.

Images