Robust relative navigation method for aircraft based on hybrid distribution under non-gaussian noise

Abstract

This invention discloses a robust relative navigation method for aircraft based on a hybrid distribution under non-Gaussian noise, belonging to the technical field of computation, estimation, or counting. To address the problem of filter divergence caused by non-stationary heavy-tail noise in time-varying environments during relative navigation, this invention introduces a Dirichlet random mixture vector fusion of Gaussian, Student's t, and multivariate K-distributions, proposing a Gaussian-Student's t-multivariate K-distribution modeling of measurement likelihood. Then, by minimizing the KLD of the true posterior probability density function and the approximate posterior probability density function using variational Bayesian techniques, the approximate posterior estimates of the aircraft's relative motion state and filter parameters are obtained, yielding the target's state information relative to the aircraft and solving for relative position and velocity. Finally, a nonlinear filter based on the Gaussian-Student's t-multivariate K-distribution is derived to improve relative navigation accuracy for angle-only relative navigation of aircraft in time-varying environments.

Description

Technical Field

[0001] This invention pertains to aircraft navigation technology, specifically to measurement and relative navigation technology, and particularly to robust aircraft filtering technology. It discloses a robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, belonging to the technical field of calculation, estimation, or counting. Background Technology

[0002] Angle measurement navigation has become an important method for relative navigation between aircraft due to its advantages such as low cost, small size, and long detection range. Angle measurement navigation relies on optical cameras, which are susceptible to optical vignetting and stray light within the field of view. During actual flight, optical cameras experience environmental changes such as front lighting and backlighting. Measurement noise is low and close to Gaussian noise under front lighting, while outliers and heavy-tail noise characteristics exist under backlighting. Therefore, the corresponding noise characteristics during actual flight may also be non-stationary heavy-tail noise that slowly changes within Gaussian and heavy-tail patterns, leading to filter divergence. How to solve this problem has become one of the current research hotspots in the field of relative navigation.

[0003] Among the existing methods for suppressing heavy-tail noise, VBAKF (Variational Bayesian Adaptive Kalman Filter) and RSKF (Robust Student't Kalman Filter) are representative. The VBAKF algorithm is an adaptive variational Bayesian Kalman filter based on Gaussian modeling of measurement noise. This method relies on the Gaussian assumption and cannot handle non-Gaussian noise well. RSKF is a robust variational Bayesian filter algorithm based on Student't modeling. However, during target tracking with only angle measurement in an aircraft, time-varying environmental changes such as front lighting and backlighting cause variations in the degree to which the sensor is affected by outliers. That is, the measurement noise may be non-stationary heavy-tail noise that slowly changes within Gaussian and heavy-tail patterns. A single Student't distribution is insufficient to accurately characterize the non-stationary nature of the noise. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of the aforementioned background technology by providing a robust relative navigation method for aircraft based on a hybrid distribution under non-Gaussian noise. This method solves the technical problem of filter divergence caused by non-stationary heavy tail noise in time-varying environments during relative navigation of aircraft, thereby achieving the invention's objective of high-precision relative navigation for aircraft.

[0005] To achieve the above-mentioned objectives, the present invention employs the following technical solution:

[0006] This invention discloses a robust relative navigation method for aircraft based on a hybrid distribution in non-stationary heavy-tail noise environments. First, an angle measurement model based on optical sensors is established. The target's position vector is then calculated from the target's pitch and azimuth angles in the measurement model. The state equations and observation equations of the relative navigation system are input into a nonlinear navigation filter based on a hybrid distribution to obtain the target's position and velocity information relative to the aircraft, thereby achieving robust relative navigation of the aircraft based on a hybrid distribution in non-stationary heavy-tail noise environments.

[0007] A robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise includes the following steps:

[0008] Step 1: Introduce a Dirichlet random mixture vector representing the Gaussian distribution weights, the Student's t-distribution weights, and the multivariate K-distribution weights. Construct a joint posterior probability density function of the state variables, the mixture vector, and the observation vector under the constraint of measurement noise variance. The joint posterior probability density function follows a Gaussian-Student's t-multivariate K-mixture distribution. Model the statistical prior probability density function of the measurement noise variance as an inverse Wissaud distribution. Rewrite the joint posterior probability density function as a Gaussian hierarchical form.

[0009] Step 2: Obtain the state variables and state variance of the previous time step, the measurement and noise variance of the current time step, and the mixed distribution parameters;

[0010] Step 3: The capacitive Kalman filter algorithm is used for time update and measurement update. During the measurement update process, the relative operating state quantity and the posterior estimate of the mixed distribution parameter are approximated by the prior probability density function of the measurement noise variance statistics of the inverse Wissaud distribution and the joint posterior probability density function in Gaussian hierarchical form through variational Bayesian inference.

[0011] As a further optimization scheme for the robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, the joint posterior probability density function of the state variables, mixed vector, and observation vector under measurement noise variance constraints constructed in step one is: Among them, z k Let x be the observation vector at time k. k Let τ be the state variable at time k. k Let R be the mixture vector at time k. k Let p(z) be the measurement noise variance at time k. k |x k ,τ k ,R k ) is x k ,τ k ,R k z under constraints k The joint posterior probability density function, τ k,1 τ k,2 τk,3 Let N(x) be the weights of the Gaussian, Student t, and multivariate K distributions, respectively, and let N(x), St(x), and MK(x) represent the probability density functions of the Gaussian, Student t, and multivariate K distributions, respectively. k ) represents the mean vector composed of the average values ​​of the observation vectors at time k, and v and a represent the degrees of freedom parameters of the student t-distribution and the multivariate K-distribution, respectively.

[0012] As a further optimization of the robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, step one models the statistical prior probability density function of the measurement noise variance as an inverse Wissaud distribution, specifically: p(R k ) = IW(R k ;u k|k-1 U k|k-1 ), where p(R) k Let IW(k) be the statistical prior probability density function of the measurement noise variance at time k, and let IW(k) be the probability density function of the inverse Wissaud distribution. k|k-1 U represents the one-step predicted value of the degree-of-freedom parameter of the inverse Wissaud distribution at time k-1. k|k-1 This represents a one-step prediction of the inverse scaling matrix of the inverse Wissaud distribution at time k-1. Forgetting factor, u k-1|k-1 U represents the degree of freedom parameter of the inverse Wieshard distribution at time k-1. k-1|k-1 Let represent the inverse scaling matrix of the inverse Wissaud distribution at time k-1.

[0013] As a further optimization of the robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, step one rewrites the joint posterior probability density function into a Gaussian hierarchical form, specifically as follows: Where, λ k Let λ be a vector that follows a categorical distribution at time k. k,1 , λ k,2 , λ k,3 For λ k The first element, the second element, the third element, λ k,1 , λ k,2 , λ k,3 The vector ω is composed of a single element that is 1 and all others are 0. k To generate a density generator for the scaling function of the student's t-distribution at time k, δ k To generate a density generator for the scaling function of a multivariate K-distribution at time k, G(;) and IG(;) represent the probability density functions of the gamma and inverse gamma distributions, respectively, v1, v2 and a1, a2 represent the shape parameters of the gamma and inverse gamma distributions, respectively, and p(λ k |τ k ) is τ k λ under constraintsk The probability density function of τ is Cat(;), which is the probability density function of the classification distribution. k ) is the mixing vector τ at time k. k The probability density function of the Dirichlet distribution is α, where Dir(;) is the probability density function of the Dirichlet distribution. k|k-1 It is a one-step prediction of the concentration parameter matrix at time k-1, where N is the number of mixture elements.

[0014] As a further optimization scheme for the robust relative navigation method for aircraft based on hybrid distribution under non-Gaussian noise, the hybrid distribution parameters obtained in step two include: the one-step predicted value u of the degree of freedom parameter of the inverse Wissaud distribution at time k-1. k|k-1 One-step prediction of the inverse scaling matrix of the inverse Wissaud distribution at time k-1 k|k-1 Forgetting factors The shape parameters v1, v2 of the gamma distribution and the shape parameters a1, a2 of the inverse gamma distribution.

[0015] As a further optimization scheme for the robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, in step three:

[0016] The capacitive Kalman filter algorithm is used for time updates, specifically:

[0017] Calculate the volume point of the state variables:

[0018] One-step prediction of state variables: One-step variance prediction: The measurement update is performed using a capacitive Kalman filter algorithm combined with variational Bayesian inference, specifically: calculating the capacitive points for one-step prediction of the state variables.

[0019]

[0020] Auxiliary matrix A at time k k renew: Noise variance measured at time k renew: R k =U k|k / u k|k ,

[0021]

[0022] One-step prediction of observation vectors: One-step prediction of measurement variance matrix and measurement covariance matrix:

[0023] Calculate the filter gain K at time k. k :

[0024] Update the state at time k The state variance P at time k k|k :

[0025] in, Let be the l-th volume point of the state variables at time k-1, where 'chol()' denotes the Cholesky decomposition, 'n' represents the dimension of the state vector, and 'e' represents the state vector dimension. l This represents the l-th column vector of an n-dimensional identity matrix. For one-step prediction of the state variables at time k-1, f(;) represents the state transfer function, P k|k-1 This is a one-step prediction of the state variance at time k-1. For the l-th volume point of the state variable predicted in one step at time k-1, Let h(;) be the estimated state vector at time k, and let E[λ] be the nonlinear measurement function. k,1 ]、E[λ k,2 ]、E[λ k,3 ]、E[ω k ]、E[δ k ]

[0026] For λ k,1 , λ k,2 , λ k,3 ω k δ k The expectation, tr(); represents the trace of the matrix. This is a one-step prediction of the observation vector at time k-1. This is a one-step prediction of the variance matrix measured at time k-1. This is a one-step prediction of the covariance matrix measured at time k-1.

[0027] As a further optimization scheme for the robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, E[λ k ]、E[ω k ]、E[δ k ]、E[logτ k,j Based on the prior probability density function of the measurement noise variance statistics of the inverse Wittsch distribution and the joint posterior probability density function in Gaussian hierarchical form, and approximated by variational Bayesian inference, the specific details are as follows:

[0028] The joint posterior probability density function is approximated using the free factor distribution. Δ k ={x k ,R k ,λ k ,τ k ,ω k,δ k}, z 1:k Let z be the observation vector from time 1 to time k. 1:k-1 Let p(x) be the observation vector from time 1 to time k-1. k |z 1:k-1 ) for z 1:k-1 x under constraints k The probability density function;

[0029] By minimizing the Coulbert-Leibler divergence cost function, a suboptimal approximation of the joint posterior probability density function is obtained: q(Δ k )=argminKLD{p(Δ k ,z 1:k )||q(Δ k ,z 1:k )}, q(;) represents the approximate posterior probability distribution of p(;);

[0030] The general solution for the mixed distribution parameters is obtained using the VB method: logq(ζ) k )=Ε[logp(Δ k ,z 1:k )]+c ζk , ζ k For set Δ k The elements in the array, log(), are the natural logarithm function, and c ζk As a constant, ψ(;) is the degamma function, K a (;), K a+1 (;) represent the corrected Bessel functions for the a-th moment and the (a+1)-th moment, respectively, τ k,j Let α be the j-th element of the mixture vector at time k. k|k,j Let α be the j-th element of the concentration parameter matrix at time k. k|k =α k|k-1 +Ε[λ k ],

[0031] An electronic device includes a memory and a processor, wherein the memory stores a computer program that runs on the processor, and the processor executes the steps of the relative navigation method described above when running the computer program.

[0032] A computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed, performs the steps of the relative navigation method described above.

[0033] The present invention, by adopting the above technical solution, has the following beneficial effects:

[0034] (1) This invention addresses the problem of filter divergence caused by non-stationary heavy tail noise in time-varying environments. It improves existing heavy tail noise suppression methods. Specifically, in the volume Kalman filter measurement update, the variance of non-stationary measurement noise is updated by updating the mixed distribution parameters of the measurement posterior probability density function, so that the aircraft can adapt to the time-varying environment and improve relative navigation accuracy.

[0035] (2) The robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise disclosed in this invention has excellent generalization ability because this invention proposes a general filtering framework for nonlinear state-space models with non-stationary heavy-tailed measurement noise. By embedding different Gaussian weighted integral techniques, different computational complexities and linearization accuracies can be obtained. Attached Figure Description

[0036] Figure 1 This is a schematic diagram of relative navigation for aircraft using angle measurement.

[0037] Figure 2 This is a graphical representation of a hierarchical Gaussian state-space model based on a mixed distribution.

[0038] Figures 3(a)-3(f) This is a comparison of the simulation results of the present invention with those of other navigation filtering algorithms under non-stationary heavy-tail noise conditions, involving six state variables.

[0039] Figure 4 The simulation results show the weight changes of the Gaussian-Student t-multivariate K-type distributions in the mixed distribution proposed in this invention under non-stationary heavy-tailed noise conditions.

[0040] Figure 5 This is a flowchart of the single-step operation of the robust relative navigation filter for aircraft based on hybrid distribution according to the present invention. Detailed Implementation

[0041] The technical solution of the invention will now be described in detail with reference to the accompanying drawings.

[0042] like Figure 1 As shown, this example focuses on the relative navigation scenario of aircraft angle measurement in a time-varying environment. It uses the measurement information of an optical camera and employs a nonlinear filter based on a GSK (Gaussian-Student's-Multivariate K) hybrid distribution for filtering and solving. Fixed-point iteration is performed inside the filter to adaptively adjust the parameters of different distributions, thereby achieving high-precision estimation of the target position and velocity of the aircraft and completing high-precision relative navigation under non-stationary and non-Gaussian noise conditions.

[0043] The specific implementation method of this example is as follows:

[0044] Step 1: Establish the angle measurement model:

[0045]

[0046] In angular relative navigation only, the observation vector z at each moment consists of two angle measurements, denoted here as azimuth θ and pitch φ, which correspond to the straight-line relative position vector δr from the aircraft to the target under the aircraft's own system RTN: δr=[δr R δr T δr N ] T ,δr R δr T δr N These are the components of the relative position vector of the aircraft pointing towards the target along the straight line in the aircraft's own system, corresponding to the R-axis, T-axis, and N-axis. For example... Figure 1 As shown, in the RTN system of the aircraft, the R-axis points from the Earth's center to the aircraft, the T-axis is along the tangent direction of the flight path, and the N-axis is along the normal direction, forming a right-handed relationship with the R-axis and T-axis.

[0047] Solve for the relationship between measurements and changes in state variables:

[0048]

[0049] In equation (2), C is the measurement matrix and the state variable x is the flight trajectory parameter of the relative motion of the aircraft.

[0050] Step 2: Construct a nonlinear filter with GSK hybrid distribution

[0051] (1) Introducing the Dirichlet random mixture vector τ, a new GSK mixture distribution is obtained, and the joint probability density function (PDF) of GSK is:

[0052]

[0053] In equation (3), z k Let x be the observation vector at time k. k Let τ be the state variable at time k. k Let R be the mixture vector at time k. k Let p(z) be the measurement noise variance at time k. k |x k ,τ k ,R k ) is x k ,τ k ,R k z under constraints k The joint posterior probability density function, τ k,1 τ k,2 τ k,3Let τ be the first, second, and third elements of the mixing vector at time k. k,1 τ k,2 τ k,3 Let N(x), St(t), and MK(t) represent the weights of the Gaussian, Student t, and multivariate K distributions, respectively. k The vector τ represents the mean vector composed of the average values ​​of the observation vectors at time k. The mean vector provides the average value of each observation. v and a represent the degrees of freedom (dof) parameters of the Student's t-distribution and the multivariate K-distribution, respectively. The mixture vector τ at time k is... k The probability density function p(τ) k It follows a Dirichlet distribution:

[0054] Dir(τ k ;α k|k-1 ,N), N=3 (4)

[0055] In equation (4), Dir(;) represents the PDF of the Dirichlet distribution, α k|k-1 It is a one-step prediction of the concentration parameter matrix at time k-1, where N is the number of mixture elements.

[0056] The elliptic distribution is a natural extension of the Gaussian distribution and includes various non-Gaussian distributions with heavy-tailed characteristics. It can be represented in Gaussian form by introducing a density function and a scaling function. The table below gives the corresponding density function and scaling function for Student t and multivariate K.

[0057] Table 1. Density function and scaling function of student t and multivariate K distribution.

[0058] distributed Scale function Density function ST <![CDATA[ω k ]]> <![CDATA[G(v1,v2)]]> MK <![CDATA[δ k ]]> <![CDATA[IG(a1,a2)]]>

[0059] Where, ω k To generate a density generator for the scaling function of the student's t-distribution at time k, δ k To generate the density generator for the scaling function of the multivariate K-distribution at time k, G(;) and IG(;) represent the gamma distribution PDF and the inverse gamma distribution PDF, respectively, and v1, v2 and a1, a2 represent the shape parameters of the gamma and inverse gamma distributions, respectively. Therefore, the Student's t-distribution and the multivariate K-distribution can be expressed as:

[0060]

[0061] (2) The statistical prior of the measurement noise variance is modeled as an inverse Wishart distribution;

[0062] p(R k ) = IW(R k ;u k|k-1 U k|k-1 (6)

[0063] In equation (6), p(R) k Let IW(;) be the statistical prior probability density function of the variance of the observation vector measurement noise at time k, and let IW(;) be the inverse Wishart PDF. k|k-1 U represents the one-step prediction of the degree of freedom parameter of the inverse Wishart distribution at time k-1. k|k-1 This represents a one-step prediction of the inverse scaling matrix of the inverse Wishart distribution at time k-1. Considering that the measurement noise covariance matrix is ​​relatively stable or changes slowly, a forgetting factor is used. in u k-1|k-1 U represents the degree of freedom parameter of the inverse Wishart distribution at time k-1. k-1|k-1 Let represent the inverse scaling matrix of the inverse Wishart distribution at time k-1.

[0064] (3) Rewrite the joint posterior probability density function in Gaussian hierarchical form, such as Figure 2 As shown, a vector λ = [λ1, λ2, ..., λ] that follows a categorical distribution is introduced. N ], where λ is a unit vector with only one element equal to 1 and all other elements equal to 0, λ1, λ2, ..., λ N For the 1st, 2nd, ..., Nth element.

[0065] p(λ|τ)=Cat(λ;τ,N) (7)

[0066] In equation (7), Cat(;) is the PDF of the classification distribution;

[0067] like Figure 2 As shown, the joint posterior probability density function is rewritten in Gaussian hierarchical form:

[0068]

[0069] In equation (8), λ k Let λ be a vector that follows a categorical distribution at time k. k,1 , λ k,2 , λ k,3 For λ k The first element, the second element, the third element, p(λ) k |τ k ) is τ k λ under constraints k The probability density function. It can be seen that through λ k,1 , λ k,2 , λ k,3 The weighted summation form of the GSK mixture distribution in equation (3) is transformed into the product form of the GSK mixture distribution in equation (8).

[0070] (4) Approximating the posterior estimates of relative motion state variables and parameters using variational Bayesian techniques:

[0071] By approximating the joint distribution using the free factor distribution, the joint posterior PDF is rewritten:

[0072]

[0073] In equation (9): Δ k ={x k ,R k ,λ k ,τ k ,ω k ,δ k}, z1:k is the observation vector from time 1 to time k, z 1:k-1 Let p(x) be the observation vector from time 1 to time k-1. k |z 1:k-1 ) for z 1:k-1 x under constraints k The probability density function;

[0074] Minimizing the Courbet-Leibler divergence (KLD) cost function yields a suboptimal approximation of the joint posterior probability density function:

[0075] q(Δ k )=argminKLD{p(Δ k ,z 1:k )||q(Δ k ,z 1:k )} (10)

[0076] In equation (10), q(;) represents the approximate posterior probability distribution of p(;).

[0077] Using the VB method, the general solution is as follows:

[0078]

[0079] In equation (11), It is a constant, ζ k It is a set Δ k In the array, log(;) represents the natural logarithm function, and Ε[;] represents the expectation.

[0080] Step 3: Input the state equation and measurement equation of the relative navigation system into the navigation filter of Step 2 to obtain a high-precision estimate of the relative motion state, and then obtain the position and velocity information of the aircraft relative to the target, thereby realizing the relative navigation of the aircraft.

[0081] like Figure 5 As shown, the specific steps are as follows:

[0082] (1) Initialization: Input the state variables and state variance of the system at the previous time step, the measurement and noise variance at the current time step, and various parameters required for the mixed distribution:

[0083]

[0084] in, For the estimation of the state variables at time k-1, P k-1|k-1 Let z represent the state variance at time k-1. k Q represents the measurement at time k. k-1 R represents the variance of the process noise at time k-1. k U represents the measurement noise variance at time k-1. k|k-1 U represents the prior value of the degrees of freedom parameter of the inverse Wishart distribution at time k. k|k-1 Let represent the inverse scaling prior matrix of the inverse Wishart distribution at time k. , where v1, v2 and a1, a2 are forgetting factors, and v1, v2 and a1, a2 are shape parameters of the gamma distribution and the inverse gamma distribution, respectively.

[0085] (2) Time update:

[0086] The weighted integral is achieved using the third-order radial spherical radial volume rule:

[0087] Calculate the state volume point:

[0088]

[0089] In equations (12) and (13), Let be the l-th volume point for the state quantity estimation at time k-1, where 'chol()' denotes the Cholesky decomposition, 'n' represents the dimension of the state vector, and 'e' represents the state vector dimension. l Let l represent the l-th column vector of the n-dimensional identity matrix.

[0090] State prediction in one step:

[0091]

[0092] In equation (14), For the one-step prediction of the state variables at time k-1, f(;) represents the state transfer function.

[0093] One-step variance prediction:

[0094]

[0095] In equation (15), P k|k-1 This represents a one-step prediction of the state variance at time k-1.

[0096] (3) Measurement update:

[0097] Iteration parameter initialization:

[0098] Calculate the volume point of the state in one step.

[0099]

[0100] In equation (16), This is the l-th volume point for predicting the state variables at time k-1.

[0101] Auxiliary matrix A at time k k renew:

[0102]

[0103] In equation (17), h(;) represents a nonlinear measurement function.

[0104] Based on the variational Bayesian inference process of equations (9) to (11), we obtain λ. k τ k ω k δ k The update formula is as follows:

[0105]

[0106] In equation (18), ψ(;) represents the degamma function, K a (;), K a+1 (;) represent the corrected Bessel functions for the a-th moment and the (a+1)-th moment, respectively, τ k,j Let α be the j-th element of the mixture vector at time k. k|k,j Let be the j-th element of the concentration parameter matrix at time k, and N be the number of elements in the concentration parameter matrix at time k. The concentration parameter matrix α at time k... k|k The calculation formula is as follows:

[0107] α k|k =α k|k-1 +Ε[λ k (19)

[0108]

[0109] In equation (20), m represents the number of observations, and tr(;) represents the trace of the matrix.

[0110]

[0111] In equation (21), exp(;) represents the natural base.

[0112]

[0113] The measurement noise variance update is obtained:

[0114]

[0115] In equation (23), u k|k U represents the degree of freedom parameter of the inverse Wishart distribution at time k. k|k Let represent the inverse scaling matrix of the inverse Wishart distribution at time k.

[0116] R k =U k|k / u k|k

[0117]

[0118] In equation (24), Update the measurement noise variance at time k.

[0119] Calculate the one-step predicted measurement value:

[0120]

[0121] In equation (25), This is a one-step prediction of the observation vector at time k-1. This is the l-th volume point predicted in one step for the state variables at time k-1.

[0122] Calculate the variance matrix of one-step prediction measurement Covariance Matrix

[0123]

[0124]

[0125] Calculate the filter gain K at time k. k :

[0126]

[0127] Update the state at time k The state variance P at time k k|k :

[0128]

[0129] The above steps complete one single-step operation of the nonlinear filter based on a hybrid distribution. The final state quantity at time k is obtained through fixed-point iteration, and then the relative motion state of the aircraft is calculated from the state quantity at time k. The simulation parameters are set as follows:

[0130] parameter numerical values Initial state of the aircraft (7200000m,0.5,30deg,60deg,120deg,180deg) Initial relative motion state (-150,-20000,300,0-300,0) / m Initial state variance (2500,20000,1000,1000,1000,1000) / m Total filtering time 3000s Measurement noise 18arcsec Process noise [10m, 0.1m / s]

[0131] Non-stationary heavy-tailed measurement noise v k The following results were generated:

[0132]

[0133] In equation (31), wp means "there is a probability", N represents the normal distribution, R represents the noise mean, and p1 and p2 represent the abnormal probabilities in the time intervals [1001,2000] and [2001,3000], respectively.

[0134] Simulation results are as follows Figures 3(a) to 3(f) As shown, the root mean square error of the estimation of the six state variables of the aircraft by three filters under the condition of non-Gaussian non-stationary noise is plotted. The six state variables represent flight trajectory parameters and uniquely determine the relative position and velocity of the aircraft and the target. It can be seen from the figure that the robust relative navigation method based on hybrid distribution proposed in this invention has better accuracy in estimating the relative motion state of the target than VBAKF and RSKF. Figure 4 The changes in the weights of the three distributions in the GSK hybrid distribution under non-stationary heavy-tailed noise were plotted. As can be seen from the figure, the relative navigation method proposed in this invention can adaptively adjust the weights of the hybrid distribution to cope with time-varying noise environments.

[0135] In summary, this invention addresses the non-stationary heavy tail noise caused by aircraft operating in complex time-varying environments by proposing a robust relative navigation method based on a hybrid distribution. The effective application and implementation of this technology has significant theoretical and practical implications for reducing platform power consumption, adapting to more complex and variable operating environments, and improving navigation accuracy.

Claims

1. A robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise, characterized in that, Includes the following steps: Step one involves introducing a Dirichlet random mixture vector representing the Gaussian, Student's t-distribution, and multivariate K-distribution weights. This constructs a joint posterior probability density function for the state variables, the mixture vector, and the observation vector under the constraint of measurement noise variance. This joint posterior probability density function follows a Gaussian-Student's t-multivariate K-mixture distribution. The statistical prior probability density function of the measurement noise variance is modeled as an inverse Wissaud distribution. Finally, the joint posterior probability density function is rewritten in a Gaussian hierarchical form. The joint posterior probability density function of the constructed state variables, mixture vector, and observation vector under the constraint of measurement noise variance is: , Let k be the observation vector at time k. Let k be the state variable at time k. Let be the mixing vector at time k. Let k be the measurement noise variance. for Under constraints The joint posterior probability density function, , , The weights are for the Gaussian distribution, Student's t-distribution, and multivariate K-distribution. , and Let represent the probability density functions of Gaussian, Student t, and multivariate K, respectively. This represents the mean vector composed of the average values ​​of the observation vectors at time k. , Let represent the degrees of freedom parameters of the student's t-distribution and multivariate K-distribution, respectively. The statistical prior probability density function of the measurement noise variance is modeled as an inverse Wissaud distribution, specifically: , Let be the statistical prior probability density function of the measurement noise variance at time k. Let be the probability density function of the inverse Wissaud distribution. This represents the one-step prediction of the degree-of-freedom parameter of the inverse Wissaud distribution at time k-1. This represents a one-step prediction of the inverse scaling matrix of the inverse Wissaud distribution at time k-1. , , Forgetting factor, , This represents the degree of freedom parameter of the inverse Wissaud distribution at time k-1. Let the inverse scaling matrix of the inverse Wissaud distribution at time k-1 be denoted as . The joint posterior probability density function is rewritten in Gaussian hierarchical form as follows: , Let be a vector that follows a categorical distribution at time k. , , for The first element, the second element, the third element, , , This forms a unit vector with only one element being 1 and all others being 0. To generate a density generator for the scaling function of the student's t-distribution at time k, To generate a density generator for the scaling function of the multivariate K-distribution at time k, and Represent the probability density function of the gamma distribution and the probability density function of the inverse gamma distribution. and Let these represent the shape parameters of the gamma distribution and the inverse gamma distribution, respectively. for Under constraints The probability density function, Let be the probability density function of the classification distribution. The mixing vector at time k The probability density function, Let be the probability density function of the Dirichlet distribution. It is a one-step prediction of the concentration parameter matrix at time k-1. It is the number of hybrid elements; Step two: Obtain the state variables and state variance of the previous time step, the measurement and noise variance of the current time step, and the mixed distribution parameters. The obtained mixed distribution parameters include: the one-step predicted value of the degree of freedom parameters of the inverse Wissaud distribution at time k-1. One-step prediction of the inverse scaling matrix of the inverse Wissaud distribution at time k-1 Forgetting factors Gamma distribution shape parameters and inverse gamma distribution shape parameters ; Step three involves using a capacitive Kalman filter algorithm for time and measurement updates. During the measurement update process, the prior probability density function of the measurement noise variance statistics based on the inverse Wissaud distribution and the joint posterior probability density function in Gaussian hierarchical form are used, along with variational Bayesian inference to approximate the posterior estimates of the relative operating state quantities and the mixed distribution parameters. The capacitive Kalman filter algorithm is used for time updates, specifically: Calculate the volume point of the state variables: , , One-step prediction of state variables: , One-step variance prediction: ; The measurement update is performed using a capacitive Kalman filter algorithm combined with variational Bayesian inference, specifically as follows: Calculate the volume point for one-step prediction of state variables: ， Timing Auxiliary Matrix renew: , Time-based noise variance renew: , , , , , One-step prediction of observation vectors: One-step prediction of measurement variance matrix and measurement covariance matrix: , , Calculate the filter gain at time k : , Update the state at time k and the variance of the state at time k : , ; for The l-th volume point of the state quantity at time step, This represents the Cholesky decomposition. This represents the dimension of the state vector. express The first dimension of the identity matrix Column vectors This is a one-step prediction of the state variables at time k-1. Represents the state transfer function. This is a one-step prediction of the state variance at time k-1. for Predict the l-th volume point of the state variables at time -1. for State vector estimate at time step Represents a nonlinear measurement function. , , , , for , , , , The expectation, tr(); represents the trace of the matrix. This is a one-step prediction of the observation vector at time k-1. This is a one-step prediction of the variance matrix measured at time k-1. This is a one-step prediction of the covariance matrix measured at time k-1. The , , , Based on the prior probability density function of the measurement noise variance statistics of the inverse Wieshart distribution and the joint posterior probability density function in Gaussian hierarchical form, and approximated by variational Bayesian inference, the specific approximation is as follows: The joint posterior probability density function is approximated using the free factor distribution. , , Let k be the observation vector from time 1 to time k. Let be the observation vector from time 1 to time k-1. for Under constraints The probability density function, By minimizing the Courbet-Leibler divergence cost function, a suboptimal approximation of the joint posterior probability density function is obtained: , express The approximate posterior probability distribution, Obtain the general solution for the mixed distribution parameters using the VB method: , For set The elements in It is the natural logarithm function. As a constant, , For the degamma function, , Let represent the corrected Bessel functions for the a-th moment and the (a+1)-th moment, respectively. for The j-th element of the mixing vector at time step 1. Let j be the j-th element of the concentration parameter matrix at time k. , , , Indicates the number of observations.

2. An electronic device comprising a memory and a processor, wherein the memory stores a computer program that runs on the processor, and the processor executes the steps of the relative navigation method of claim 1 when running the computer program.

3. A computer-readable storage medium having a computer program stored thereon, the computer program executing the steps of the relative navigation method of claim 1 when run.

Images