Installation error calibration method based on Student's T distribution and variational Bayes

Abstract

The invention discloses an installation error calibration method based on Student's T distribution and variational Bayes, which comprises the following steps: firstly, constructing an installation error calibration geometric model of an SINS / USBL system, defining coordinate systems, establishing an attitude transfer matrix between the coordinate systems, and designing a state equation and a measurement equation by taking installation error angles in three directions as state variables; time updating and measurement updating are carried out based on a Kalman filtering framework; the method comprises the following steps of: embedding Student's T distribution into a variational Bayesian filtering framework, alternately updating distribution parameters of a state variable, a noise covariance and an auxiliary variable through a variational iterative optimization process, maximizing a variational lower bound until convergence, and outputting an optimized installation error angle estimated value. According to the method, acoustic measurement noise is modeled by using the heavy tail characteristic of the SINS / USBL combined system, the interference of outliers on installation error angle estimation is remarkably inhibited, the positioning accuracy of the SINS / USBL combined system in a complex underwater environment can be effectively improved, and the calibration robustness is enhanced.

Description

Technical Field

[0001] The present invention belongs to the technical field of underwater navigation, and particularly relates to a calibration method and system for installation errors based on Student's T distribution and variational Bayesian, which is applicable to suppressing the interference of acoustic measurement outliers and improving the calibration accuracy of integrated navigation systems in complex underwater environments. Background Technique

[0002] Traditional SINS / USBL calibration methods are mostly based on Kalman filtering under Gaussian assumptions. However, in actual underwater environments, acoustic measurements are easily contaminated by outliers, resulting in a significant decline in calibration accuracy. Existing technologies such as the least squares method and the extended Kalman filter have problems such as being sensitive to outliers and relying on the prior distribution of noise. For example, traditional calibration methods use the least squares method to solve the installation error angle, but multiple navigation iterations are required; some research has proposed iterative extended Kalman filtering, which is still limited by the Gaussian noise assumption. The present invention breaks through the limitations of traditional methods by integrating Student's T distribution and variational Bayesian framework. Summary of the Invention

[0003] To solve the above problems, the present invention discloses a calibration method for installation errors based on Student's T distribution and variational Bayesian, which can effectively improve the positioning accuracy of the SINS / USBL integrated system in complex underwater environments and enhance the calibration robustness.

[0004] To achieve the above object, the technical solution of the present invention is as follows:

[0005] A calibration method for installation errors based on Student's T distribution and variational Bayesian includes the following steps:

[0006] Step 1: Construct a geometric model for calibrating the installation errors of the SINS / USBL system, define the Earth coordinate system (e-system), the navigation coordinate system (n-system), the vehicle body coordinate system (b-system), and the acoustic array coordinate system (u-system), and establish the attitude transformation matrix between the coordinate systems;

[0007] Step 2: Take the installation error angles in three directions as state variables, and design the state equation and the measurement equation;

[0008] Step 3: Implement the filtering algorithm, perform time update and measurement update based on the Kalman filtering framework, and embed the Student's T distribution into the variational Bayesian filtering framework. During the measurement update process, use its heavy-tailed characteristic to model the measurement noise, and combine variational Bayesian inference to dynamically jointly estimate the noise covariance matrix and the auxiliary variables;

[0009] Step 4: Through the variational iteration optimization process, alternately update the distribution parameters of the state variables, noise covariance, and auxiliary variables, maximize the variational lower bound (ELBO) until convergence, and output the optimized installation error angle estimation value.

[0010] The specific process of Step 1 is as follows:

[0011] Step 1.1: Construct a geometric model for installation error calibration and complete the definition of coordinate systems, including the Earth coordinate system (e-system): the origin is the center of the Earth, the axis points to the prime meridian, the equatorial plane, and the Earth's axis of rotation; the navigation coordinate system (n-system): the origin is the GNSS center, and the axes point to the east, north, and sky directions; the vehicle coordinate system (b-system): the origin is the center of gravity of the vehicle, and the axes point to the right, front, and upward directions; the acoustic array coordinate system (a-system): the origin is the center of gravity of the acoustic array, and the axes point to the right, front, and upward directions of the array.

[0012] Step 1.2: Establish the attitude transfer matrix. The installation error angles include the pitch angle θ x , roll angle θ y , and heading angle θ y . The conversion relationship from the u-system to the b-system is described by the rotation matrix :

[0013]

[0014] Preferably, the specific process of Step 2 is as follows:

[0015] Based on the geometric model, taking the installation error angles in three directions as state variables, the state equation is linear, and the measurement equation is a nonlinear discrete system. At this time, the state space model is:

[0016]

[0017] In the state equation, X k is the state of the system at time k, φ k|k-1 is the system state transition matrix, X k-1 is the state of the system at time k-1, Γ k|k-1 is the system noise matrix, and W k-1 is the system noise vector. In the measurement equation, Z k is the measurement at time k, H k is the measurement function with respect to the state variables, and V k is the measurement noise vector. Taking the installation error angle θ as the state variable X k , it is expressed as:

[0018] X k = [θ x θ y θ z

[0019] Since the installation error angle of the SINS / USBL integrated device no longer changes after installation, φ k|k-1 is a three-dimensional identity matrix.

[0020] Preferably, step 3 specifically includes the following process:

[0021] Step 3.1: Perform time update based on the Kalman filter framework:

[0022]

[0023] where is the state value at the previous moment k - 1; is the prior estimate value, that is, the optimal one-step prediction of X k-1 ; P k|k-1 is the prior covariance matrix of state prediction, P k-1 is the posterior covariance matrix.

[0024] Step 3.2: Initialize the state quantity, covariance, auxiliary random variable λ k the expectation of and the measurement noise

[0025]

[0026] where and are the state quantity and the covariance matrix estimate value at the initial moment respectively, E(0)[λ k is the mathematical expectation of λ k at the initial moment, is the measurement noise estimate value at the kth moment.

[0027] The probability density function of the Student's T distribution is:

[0028] p(ν k ) = St(ν k ; 0, R k , υ)

[0029] where p(·) represents the probability density function, v k is the measurement noise with zero mean and covariance R k and having the thick-tailed property, and υ is the degree-of-freedom parameter. When υ → ∞, the Student's t distribution degenerates into a Gaussian distribution; a smaller υ value gives the algorithm stronger resistance to outliers. The measurement noise whose probability density function conforms to the Student's T distribution is updated as follows:

[0030]

[0031] where, is the measurement noise estimation value of the i-th iteration, is the prior estimation value of the measurement noise. and are the state quantity and covariance matrix at the (i - 1)-th iteration respectively, is the observation matrix H k transpose.

[0032] Step 3.4: Update the unknown measurement noise covariance R based on the variational Bayesian algorithm k :

[0033]

[0034] where m is the dimension of the state quantity, and are the measurement noise covariance estimation and the observed value of the measurement noise covariance at the i-th iteration respectively, E (i) [λ k is the mathematical expectation of λ k at the i-th iteration.

[0035] Preferably, the specific steps of step 4 include the following process:

[0036] Step 4.1: Assume that the state variable follows a Gaussian distribution, and the update steps of the state variable, covariance, and Kalman gain are consistent with the traditional Kalman filter update.

[0037]

[0038] and are the state quantity and covariance matrix at the (i + 1)-th iteration respectively, is the Kalman gain at the (i + 1)-th iteration, is the observed value of the measurement noise covariance at the (i + 1)-th iteration.

[0039] Step 4.2: Assume that the noise covariance follows an inverse Wishart distribution in the variational Bayesian inference framework, and the probability density function of the auxiliary variable λ k follows a gamma distribution:

[0040]

[0041] where the shape parameter controls the skewness, kurtosis, and tail characteristics of the noise distribution, and the scale parameter controls the scale or dispersion degree of the distribution.

[0042] Let the variable be expressed as follows:

[0043]

[0044] Update according to m and the degree of freedom parameter the parameters, the process is as follows:

[0045]

[0046] According to update the proportionality parameter

[0047]

[0048] According to and perform the expectation update of the auxiliary variable:

[0049]

[0050] The above steps achieve the adaptive matching of state estimation and noise parameters.

[0051] Step 4.3: The objective function of variational iterative optimization is:

[0052] {q(X k ), q(λ k )} = argmin KLD(q(X k ), q(λ k ) || p(X k , λ k | Z 1:k ))

[0053] where the difference between two probability distributions is measured by KLD, and q(X k , λ k | Z 1:k ) and q(λ k ) q(λ k ) are iteratively optimized by minimizing the KLD distance between the probability density functions of p(X k ), q(λ k ). When KLD reaches the minimum, at this time:

[0054]

[0055] N is the set number of iterations, are the estimated values at the termination of the corresponding estimated value iterations respectively. Finally, the estimated value of the optimized installation error angle is output

[0056] The beneficial effects of the present invention are:

[0057] 1. Strong anti-outlier interference ability: By embedding the Student's T distribution into the variational Bayesian framework and using its heavy-tailed characteristics to model the acoustic measurement noise, the interference of outliers to the installation error angle estimation is significantly suppressed. Compared with the traditional Gaussian hypothesis method, in a contaminated environment where the noise variance is amplified by 100 times, the calibration error fluctuation is reduced by 52.1%, and the robustness is significantly enhanced.

[0058] 2. Adaptive estimation of dynamic noise parameters: Combining variational Bayesian inference, jointly estimate the noise covariance matrix and auxiliary variables in real time, avoiding the dependence on prior assumptions about the statistical characteristics of the noise, and adapting to the non-Gaussian and time-varying characteristics of the noise in complex underwater environments.

[0059] 3. Greatly improved calibration accuracy: Simulation experiments show that compared with the traditional Kalman filter (KF), the method of the present invention reduces the average error mean of the installation error angle from 0.3814° to 0.1206°, with the accuracy improved by 68.38%; compared with the recursive least squares method (RLS), the error mean is reduced by 91.87%, and the standard deviation is reduced by more than 75%.

[0060] 4. Reduced need for repeated calibration: Through the closed-loop parameter update mechanism and iterative optimization strategy, the adaptive matching of state estimation and noise model is realized. Experiments show that only a single voyage calibration is required to achieve the accuracy of the traditional method after multiple iterations, significantly improving the engineering efficiency and reducing the operation cost. Brief Description of the Drawings

[0061] Figure 1 is the schematic diagram of the installation error calibration method based on the Student's T distribution and variational Bayesian of the present invention;

[0062] Figure 2 is the flowchart of the installation error calibration method based on the Kalman filter of the present invention. Detailed Embodiment

[0063] The following further clarifies the present invention in conjunction with the drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and not to limit the scope of the present invention.

[0064] An installation error calibration scheme based on the Student's T distribution and variational Bayesian of the present invention, the implementation principle is as Figure 1 shown, and the experimental process is as Figure 2 shown. Its process mainly includes the following steps:

[0065] Step 1: Construct a geometric model for installation error calibration;

[0066] Step 1.1: Coordinate system definition and establishment of attitude transfer matrix:

[0067] Earth coordinate system (e - system): The origin is the center of the earth, O e x e 、O e y e and O e z e axes point to the prime meridian, the orthogonal axis in the equatorial plane, and the earth's axis of rotation (ECEF coordinate system) respectively; Navigation coordinate system (n - system): The origin is the GNSS center, O n x n 、O n y n and O n z n axes point to the east, north, and sky directions; Vehicle body coordinate system (b - system): The origin is the center of gravity of the vehicle body, O b x b 、O b y b and O b z b axes point to the right side, front, and upward direction of the vehicle body; Acoustic array coordinate system (u - system): The origin is the center of gravity of the acoustic array, O u x u 、O u y u and O u z u axes point to the right side, front, and upward direction of the array.

[0068] Step 1.2: Define the installation error angles. The installation error angles include the pitch angle θ x 、roll angle θ y 、heading angle θ z , which characterize the attitude deviation between the b - coordinate system and the u - coordinate system; The attitude transformation matrix describes the transformation relationship from the u - system to the b - system:

[0069]

[0070] Step 2: Selection of state variables and derivation of equations.

[0071] Step 2.1: Complete the definition of state variables. The installation error angles are used as state variables, i.e., X k =[θ x θ y θ z ; Since the installation error angles are constant after system installation, the state equation is in linear form:

[0072] X k =Φ k / k-1 X k-1 +Γ k / k-1 W k-1

[0073] Where, Φk / k-1 is a three-dimensional identity matrix, w k-1 is the state at time k-1 of the system, Γ k / k-1 is the system noise matrix, W k-1 is the system noise vector.

[0074] Step 2.1: Complete the construction of the measurement equation. Based on the geometric relationship between the acoustic array and the transponder, the measurement equation is a non-linear discrete system:

[0075] Z k = H k X k + V k

[0076] where H k is the observation matrix with respect to the state variables. The position of the transponder in the e-frame is subtracted from the position of the SINS / GNSS integrated system in the e-frame to form a position vector, and then the position vector is projected onto the b-frame. The position of the transponder in the u-frame is subtracted from the position of the USBL in the u-frame to form a position vector. Taking the subtraction of these two position vectors as the observation quantity, the non-linear functional expression of the observation quantity with respect to the state variables can be sorted out. In the formula

[0077] The specific expression of the observation quantity is:

[0078]

[0079] Step 3: Implement the variational Bayesian filtering algorithm;

[0080] The variational Bayesian algorithm is an approximate inference method based on probabilistic graphical models. Its core goal is to approximate the true posterior distribution by optimizing the variational distribution

[0081] Step 3.1: Modeling and parameter setting of the Student's T distribution:

[0082] Traditional KF models the system state and measurement noise based on the Gaussian distribution assumption. However, in practical applications, the measurement noise may contain outliers or spike noise, resulting in the failure of the Gaussian assumption. The Student's t distribution has heavier tails than the Gaussian distribution and can better capture extreme values or abnormal noise. Specifically, the statistical data from the noise sequence including pulse signals is used to reduce the impact of measurement noise. Model the measurement noise, and the measurement noise V k follows the Student's T distribution, and its probability density function is:

[0083] p(ν k ) = St(ν k ; 0, Rk , υ)

[0084] where v k is measurement noise with zero mean and covariance R k , having a heavy - tailed characteristic. The degrees - of - freedom parameter is set to υ = 5. Experimental verification shows that there is a slight improvement in accuracy compared to υ = 10, and it can be dynamically adjusted according to the situation in actual use. The measurement dimension m = 3, and the noise covariance matrix R k is dynamically updated.

[0085] Step 3.2: Variational Bayesian inference process:

[0086] Initialize, set the number of iterations N = 10, and the initial covariance is set to 0

[0087] Adopt the mean - field assumption and decompose the joint variational distribution into independent factors.

[0088] p(X k , R k , λ k ) ≈ q(X k )q(R k )q(λ k )

[0089] The state variable follows a Gaussian distribution; the noise covariance R k follows an inverse Wishart distribution; the auxiliary variable λ k follows a gamma distribution.

[0090]

[0091] where the shape parameter controls the skewness, kurtosis, and tail characteristics of the noise distribution, and the scale parameter controls the scale or dispersion degree of the distribution.

[0092] Update alternately:

[0093] {q(X k ), q(λ k )} = argmin KLD(q(X k ), q(λ k ) || p(X k , λ k |Z 1:k ))

[0094] Step 3.3: Iterative convergence and adaptive mechanism:

[0095] The adaptive mechanism is reflected in the dynamic estimation of noise statistics and the adjustment of robust weights. The former updates R based on the current residual k, to achieve online calibration of the noise covariance; the latter reduces the expectation E[λ k at outliers, automatically reducing the weight of the corresponding measurement to suppress outlier interference. k at outliers, automatically reducing the weight of the corresponding measurement to suppress outlier interference.

[0096] Step 4: Complete the simulation verification experiment,

[0097] Step 4.1: Experimental setup:

[0098] The starting point longitude and latitude of the underwater vehicle's trajectory are (28°, 112°), the altitude is 0m, and the simulation trajectory lasts for 1900s; noise model: the measurement noise variances of slant range and azimuth are amplified by 100 times, and the probability of outliers is 10%; comparison methods: Recursive Least Squares (RLS), Kalman Filter (KF), and the method of the present invention (ST-VBKF).

[0099] Step 4.2: Analysis of experimental results:

[0100] The calibration error of the RLS calibration algorithm θ x is 1.8656°, the calibration error of θ y is -0.5654°, the calibration error of θ z is -0.5654°, and the average error is 1.4837°; the calibration error of the KF calibration algorithm θ x is 0.3441°, the calibration error of θ y is -0.3623°, the calibration error of θ z is 0.4378°, and the average error is 0.3814°; the calibration error of the ST-VBKF calibration algorithm θ x is 0.0425°, the calibration error of θ y is -0.2160°, the calibration error of θ z is 0.1034°, and the average error is 0.1206°.

[0101] The simulation experiment shows that compared with the traditional Kalman filter (KF), the method of the present invention reduces the average error mean of the installation error angle from 0.3814° to 0.1206°, and the accuracy is improved by 68.38%; compared with the Recursive Least Squares method (RLS), the error mean is reduced by 91.87%, and the standard deviation is reduced by more than 75%.

[0102] Compared with the traditional Gaussian hypothesis method, in an outlier-contaminated environment where the noise variance is amplified by 100 times, the calibration error fluctuation is reduced by 52.1%, and the robustness is significantly enhanced.

[0103] It should be noted that the above content only illustrates the technical idea of the present invention and cannot be used to limit the protection scope of the present invention. For those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can still be made, and these improvements and refinements all fall within the protection scope of the claims of the present invention.

Claims

1. An installation error calibration method based on Student's T distribution and variational Bayesian, characterized in that: It includes the following steps: S1: Construct the installation error calibration geometric model of the SINS / USBL system, define the coordinate systems, and establish the attitude transfer matrix between the coordinate systems; S2: Based on the geometric model, with the installation error angles in three directions as state variables, design the state equation and the measurement equation; S3: Perform time update based on the Kalman filter framework: S4: Perform measurement update based on the Kalman filter framework, and embed the Student's T distribution into the variational Bayesian filter framework; during the measurement update process, use its heavy-tailed property to model the measurement noise, and combine variational Bayesian inference to dynamically jointly estimate the noise covariance matrix and the auxiliary variables; S5: Variational iterative optimization process, alternately update the state variables, covariance, and Kalman gain; S6: Reach the iteration termination condition and output the optimized estimated value of the installation error angle.

2. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, wherein: The specific process of step S1 is as follows: Step S1.1: Construct the installation error calibration geometric model and complete the definition of the coordinate systems, including the Earth coordinate system (e-system): the origin is the center of the Earth, and the axes point to the prime meridian, the equatorial plane, and the Earth's axis of rotation; the navigation coordinate system (n-system): the origin is the GNSS center, and the axes point to the east, north, and up directions; the vehicle coordinate system (b-system): the origin is the center of gravity of the vehicle, and the axes point to the right, front, and up directions; the acoustic array coordinate system (a-system): the origin is the center of gravity of the acoustic array, and the axes point to the right, front, and up directions of the array; Step S1.2: Establish the attitude transfer matrix The installation error angle is defined as the rotational offset between the b-frame and the u-frame, expressed as θ = [θ x θ y θ z T , which are respectively the pitch installation error angle θ x on the x-axis, the roll installation error angle θ y on the y-axis, and the heading installation error angle θ z on the z-axis; the conversion relationship from the u-frame to the b-frame is described by the rotation matrix :​ 3. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, characterized in that: The specific process of step S2 is as follows: Based on the geometric model, with the installation error angles in three directions as state variables, the state equation is linear, and the measurement equation is a nonlinear discrete system. At this time, the state space model is: X in the state equation k is the state of the system at time k, and φ k|k-1 is the system state transition matrix, X k-1 is the state of the system at time k-1, Γ k|k-1 is the system noise matrix, W k-1 is the system noise vector. In the measurement equation, Z k is the measurement at time k, H k is the measurement function with respect to the state variables, V k is the measurement noise vector. Taking the installation error angle θ as the state variable X k , it is expressed as: X k = [θ x θ y θ z ​ Since the installation error angle of the SINS / USBL integrated device does not change once it is installed, so φ k|k-1 is a three-dimensional identity matrix.

4. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, characterized in that: The time update described in step S3 based on the Kalman filter framework: where is the state value at the previous moment k - 1; is the prior estimate value, that is, the optimal one-step prediction of X k-1 ; P k|k-1 is the prior covariance matrix of state prediction, and P k-1 is the posterior covariance matrix.

5. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, characterized in that: The specific process of step S4 is as follows: S4.1: Initialize the state quantity, covariance, auxiliary random variable λ k expectation and measurement noise where and are the state quantity and the covariance matrix estimated value at the initial moment, respectively, and E (0) [λ k is the λ at the initial moment k mathematical expectation, is the measurement noise estimated value at the k-th moment; S4.2: The probability density function of the Student's T distribution is: p(ν k ) = St(ν k ; 0, R k , υ) where p(·) represents the probability density function, and v k is the measurement noise with zero mean and covariance R k , which has the characteristic of heavy tails, and υ is the degree-of-freedom parameter; when υ → ∞, the Student's t-distribution degenerates into a Gaussian distribution; smaller values of υ endow the algorithm with stronger ability to resist outliers; the measurement noise whose probability density function conforms to the Student's T-distribution is updated as follows: Among them, is the measurement noise estimation value of the i-th iteration, is the prior estimation value of the measurement noise; and are the state quantity and covariance matrix at the (i - 1)-th iteration respectively, is the observation matrix H k transpose; S4.3: Update the unknown measurement noise covariance R based on the variational Bayesian algorithm k : where m is the dimension of the state quantity, and are the estimated measurement noise covariance and the observed value of the measurement noise covariance at the i-th iteration, respectively, E (i) [λ k is the mathematical expectation of λ k at the i-th iteration.

6. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, characterized in that: The specific implementation of the variational iterative optimization in step S5 includes: S5.1: Assume that the state variables follow a Gaussian distribution, and the update steps of the state variables, covariance, and Kalman gain are the same as those of the traditional Kalman filter update; and are the state quantity and covariance matrix at the (i + 1)-th iteration respectively, is the Kalman gain at the (i + 1)-th iteration, is the observed value of the measurement noise covariance at the (i + 1)-th iteration; S5.2: Assume that the noise covariance follows an inverse Wishart distribution in the variational Bayesian inference framework, and the probability density function of the auxiliary variable λ k follows a gamma distribution: where the shape parameter controls the skewness, kurtosis, and tail characteristics of the noise distribution, and the scale parameter controls the scale or dispersion of the distribution; Let the variable The expression is as follows: Update according to the m and degrees of freedom parameters The parameters are updated as follows: According to Update ratio parameter According to and perform the expected update of the auxiliary variable: The above steps achieve the adaptive matching of state estimation and noise parameters.

7. The installation error calibration method based on Student's T distribution and variational Bayesian according to claim 1, characterized in that The objective function of the variational iterative optimization in step S6 is: {q(X k ),q(λ k )} = argmin KLD(q(X k ),q(λ k ) || p(X k ,λ k |Z 1:k )) where the difference between two probability distributions is measured by KLD, and by minimizing the KLD distance between the probability density function of p(X k ,λ k |Z 1:k ) and the probability density functions of the approximate posterior q(X k )q(λ k ), q(X k ) and q(λ k ) are iteratively optimized. When the KLD reaches minimization, then: N is the set number of iterations, which are the estimated values corresponding to the termination of the estimated value iteration respectively; finally, the optimized estimated value of the installation error angle is output

Images