Aircraft robust relative navigation method based on mixed distribution under non-Gaussian noise

Abstract

The invention discloses an aircraft robust relative navigation method based on mixed distribution under non-Gaussian noise, and belongs to the technical field of calculation, reckoning or counting. In order to solve the problem that non-stationary heavy tail noise caused by a time-varying environment in a relative navigation process of an aircraft generally causes filtering divergence, Dirichlet random mixed vectors are introduced to fuse Gaussian, student t and multivariate K distribution, so that a Gaussian-student t-multivariate K mixed distribution modeling measurement likelihood is provided; and then the KLD of the real posterior probability density function and the approximate posterior probability density function is minimized through the variational Bayesian technology to solve the relative motion state of the aircraft and approximate posterior estimation of filter parameters, the state quantity information of the target relative to the aircraft is obtained, and the relative position and speed are solved. Finally, a nonlinear filter based on Gaussian-student t-multivariate K mixed distribution is deduced, the relative navigation precision is improved, and the nonlinear filter is used for angle-measurement-only relative navigation of an aircraft in a time-varying environment.

Description

Technical Field

[0001] The present invention belongs to aircraft navigation technology, specifically to measurement and relative navigation technology, and especially to aircraft robust filtering technology. It discloses an aircraft robust relative navigation method based on mixed distribution under non-Gaussian noise, and belongs to the technical field of calculation, estimation or counting. Background Art

[0002] Angular navigation has become an important method of relative navigation between aircraft due to its advantages such as low cost, small size and long detection distance. Angular navigation relies on optical cameras, which are easily affected by optical vignetting and stray light in the field of view when imaging. In actual flight, optical cameras are subject to environmental changes such as front light and back light. When the light is front light, the measurement noise is small and the noise is close to Gaussian. When the light is back light, there are wild values ​​in the measurement and the noise presents heavy-tail characteristics. Therefore, the corresponding noise characteristics in the actual flight process may also be non-stationary heavy-tail noise that changes slowly in Gaussian and heavy-tail forms, resulting in 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-tailed 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's modeling. However, during the tracking of an aircraft's only angle-measuring target, changes in the time-varying environment, such as front light and back light, cause the degree to which the sensor is affected by abnormal values ​​to change, that is, the measurement noise may be a non-stationary heavy-tailed noise that changes slowly in Gaussian and heavy-tailed forms. A single Student's t distribution is difficult to accurately characterize the non-stationary characteristics of the noise. Summary of the invention

[0004] The purpose of the present invention is to address the deficiencies of the above-mentioned background technology and to provide a method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise, so as to solve the technical problem of filtering divergence caused by non-stationary heavy-tail noise caused by the time-varying environment during the relative navigation of aircraft, and to achieve the purpose of the invention of high-precision relative navigation of aircraft.

[0005] The present invention adopts the following technical solutions to achieve the above-mentioned invention object:

[0006] The present invention discloses a method for robust relative navigation of an aircraft based on mixed distribution in a non-stationary heavy-tail noise environment. First, an angle measurement model based on an optical sensor is established. The position vector of the target is solved according to the pitch angle and azimuth angle of the target in the measurement model. The state equation and observation equation in the relative navigation system are input into a nonlinear navigation filter based on mixed distribution to obtain the position and velocity information of the target relative to the aircraft, thereby realizing robust relative navigation of the aircraft based on mixed distribution in a non-stationary heavy-tail noise environment.

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

[0008] Step 1: introduce the Dirichlet random mixture vector that represents the Gaussian distribution weight, the Student t distribution weight, and the multivariate K distribution weight, construct the joint posterior probability density function of the state quantity, the mixture vector, and the observation vector under the constraint of the measurement noise variance, the joint posterior probability density function obeys the Gaussian-Student t-multivariate K mixture distribution, the statistical prior probability density function of the measurement noise variance is modeled as an inverse Wishart distribution, and the joint posterior probability density function is rewritten into a Gaussian stratified form;

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

[0010] Step three, the cubature Kalman filter algorithm is used for time update and measurement update. During the measurement update process, the prior probability density function of the measurement noise variance statistics of the inverse Wishart distribution and the joint posterior probability density function in the Gaussian stratified form are used to approximate the posterior estimates of the relative operating state quantity and the mixed distribution parameters through variational Bayesian inference.

[0011] As a further optimization scheme for the robust relative navigation method of aircraft based on mixed distribution under non-Gaussian noise, the joint posterior probability density function of the state quantity, mixed vector, and observation vector under the constraint of measurement noise variance constructed in step 1 is: Among them, z k is the observation vector at time k, x k is the state quantity at time k, τ k is the mixing vector at time k, R k is the measurement noise variance at time k, p(z k |x k ,τ k ,R k ) is x k ,τ k ,R k Under the constraint z k The joint posterior probability density function, τ k,1 , τ k,2 , τk,3 is the weight of Gaussian distribution, Student t distribution and multivariate K distribution, N(;), St(;) and MK(;) represent the probability density functions of Gaussian, Student t and multivariate K respectively, h(x k ) represents the mean vector composed of the average values ​​of the observation vectors at time k, and v and a represent the degree of freedom parameters of Student’s t distribution and multivariate K distribution, respectively.

[0012] As a further optimization scheme for the robust relative navigation method of aircraft based on mixed distribution under non-Gaussian noise, in step 1, the statistical prior probability density function of the measurement noise variance is modeled as an inverse Wishart distribution, specifically: p(R k )=IW(R k ;u k|k-1 ,U k|k-1 ), where p(R k ) is the statistical prior probability density function of the measurement noise variance at time k, IW(;) is the probability density function of the inverse Wishart distribution, u k|k-1 represents the one-step forecast value of the degree of freedom parameter of the inverse Wishart distribution at time k-1, U k|k-1 represents the one-step prediction of the inverse scale matrix of the inverse Wishart distribution at time k-1, For the forgetting factor, u k-1|k-1 represents the degree of freedom parameter of the inverse Wishart distribution at time k-1, U k-1|k-1 Represents the inverse scaling matrix of the k-1 inverse Wishart distribution.

[0013] As a further optimization scheme for the robust relative navigation method of aircraft based on mixed distribution under non-Gaussian noise, step 1 rewrites the joint posterior probability density function into a Gaussian hierarchical form, specifically: Among them, λ k is the vector that follows the classification distribution at time k, λ k,1 , k,2 , k,3 is k The first, second, and third elements of k,1 , k,2 , k,3 It is composed of a unit vector with only one element being 1 and all other elements being 0, ω k is the density generator for generating the scaling function of the Student's t distribution at time k, δ k is the density generator for generating the multivariate K distribution scaling function at time k, G(;) and IG(;) represent the probability density function of the gamma distribution and the inverse gamma distribution, v1, v2 and a1, a2 represent the shape parameters of the gamma distribution and the inverse gamma distribution, respectively, p(λ k |τ k ) is τ k Under the constraintk The probability density function of the distribution of the categorical distribution, Cat(;) is the probability density function of the categorical distribution, p(τ k ) is the mixing vector τ at time k k The probability density function of Dir(;) is the probability density function of Dirichlet distribution, α k|k-1 is the one-step prediction of the concentration parameter matrix at time k-1, and N is the number of mixing elements.

[0014] As a further optimization scheme for the robust relative navigation method of aircraft based on mixed distribution under non-Gaussian noise, the mixed distribution parameters obtained in step 2 include: the one-step prediction value u of the degree of freedom parameter of the inverse Wishart distribution at time k-1 k|k-1 , one-step prediction of the inverse scale matrix of the inverse Wishart distribution at time k-1 U k|k-1 , Forgetting Factor Gamma distribution with shape parameters v1, v2 and inverse gamma distribution with shape parameters a1, a2.

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

[0016] The cubature Kalman filter algorithm is used for time update, specifically:

[0017] Calculate the state volume point:

[0018] One-step prediction of state quantity: Variance one-step forecast: The volumetric Kalman filter algorithm is used in combination with variational Bayesian reasoning for measurement update, specifically: Calculate the volumetric point of the state quantity one-step prediction:

[0019]

[0020] Auxiliary matrix A at time k k renew: The variance of the measurement noise at time k renew: R k =U k|k / u k|k ,

[0021]

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

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

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

[0025] in, is the lth volume point of the state quantity at time k-1, chol() represents Cholesky decomposition, n represents the dimension of the state vector, e l represents the lth column vector of the n-dimensional identity matrix, is the one-step prediction of the state quantity at time k-1, f(;) represents the state transfer function, P k|k-1 is the one-step prediction of the state variance at time k-1, is the lth volume point of the one-step predicted state quantity at time k-1, is the estimated value of the state vector at time k, h(;) represents the nonlinear measurement function, E[λ k,1 ], E[λ k,2 ], E[λ k,3 ]、E[ω k ], E[δ k ]

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

[0027] As a further optimization scheme for the robust relative navigation method of aircraft based on mixed distribution under non-Gaussian noise, E[λ k ]、E[ω k ], E[δ k ], E[logτ k,j ]According to the measurement noise variance statistical prior probability density function of the inverse Wishart distribution and the joint posterior probability density function of the Gaussian stratified form and approximated by variational Bayesian inference, specifically:

[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 is the observation vector from time 1 to time k, z 1:k-1 is the observation vector from time 1 to time k-1, p(x k |z 1:k-1 ) is z 1:k-1 Under the constraint x k The probability density function of

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

[0030] The VB method is used to obtain the general solution of the mixed distribution parameters: logq(ζ k )=Ε[logp(Δ k ,z 1:k )]+c ζk , k For the set Δ k The elements in, log(;) is the natural logarithm function, c ζk is a constant, ψ(;) is the degamma function, K a (;), K a+1 (;) represent the modified Bessel function of the a-order moment and the modified Bessel function of the a+1-order moment, τ k,j is the jth element of the mixing vector at time k, α k|k,j is the jth element of the concentration parameter matrix at time k, α k|k =α k|k-1 +Ε[λ k ],

[0031] An electronic device comprises a memory and a processor. The memory stores a computer program running on the processor. When the processor runs the computer program, the steps of the relative navigation method are executed.

[0032] A computer-readable storage medium stores a computer program, which executes the steps of the relative navigation method when the computer program is running.

[0033] The present invention adopts the above technical solution and has the following beneficial effects:

[0034] (1) The present invention aims to solve the problem of filter divergence caused by non-stationary heavy-tailed noise in a time-varying environment, and improves the existing heavy-tailed noise suppression method. Specifically, in the measurement update of the volumetric Kalman filter, the non-stationary measurement noise variance 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 the relative navigation accuracy.

[0035] (2) The robust relative navigation method for aircraft based on mixed distribution under non-Gaussian noise disclosed in the present invention has excellent generalization ability, because the present invention proposes a universal filtering framework for nonlinear state space models with non-stationary heavy-tailed measurement noise, which embeds different Gaussian weighted integration techniques to obtain different computational complexities and linearization accuracies. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 Schematic diagram of relative navigation of aircraft angle measurement.

[0037] Figure 2 A graphical model diagram of the hierarchical Gaussian state-space model based on mixture distribution.

[0038] Figure 3(a)-Figure 3(f) This is a comparison chart of simulation results of six state quantities of the present invention and other navigation filtering algorithms under non-stationary heavy-tail noise conditions.

[0039] Figure 4 These are simulation results of weight changes of the three distributions, Gaussian, Student's t, and multivariate K, in the mixed distribution proposed by the present invention under non-stationary heavy-tailed noise conditions.

[0040] Figure 5 The present invention is a single-step operation flow chart of the aircraft robust relative navigation filter based on mixed distribution. DETAILED DESCRIPTION

[0041] The technical solution of the invention is described in detail below in conjunction with the accompanying drawings.

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

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

[0044] Step 1: Create an angle measurement model:

[0045]

[0046] In angle-only relative navigation, the observation vector z at each moment consists of two angle measurements, here represented as the azimuth angle θ and the pitch angle φ, which correspond to the linear relative position vector δr from the aircraft to the target in the aircraft's own system RTN: δr = [δr R δr T δr N ] T ,δr R ,δr T ,δr N are the components of the linear relative position vector of the aircraft pointing to the target in the aircraft's own system, namely, the R-axis, T-axis, and N-axis. Figure 1 As shown, in the aircraft's original system RTN, the R axis points from the center of the earth to the aircraft, the T axis is along the tangent direction of the flight path, and the N axis is in a right-hand system with the R axis and the T axis along the normal direction.

[0047] Solve the relationship between measurement and state quantity change:

[0048]

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

[0050] Step 2: Construct a nonlinear filter for the GSK mixture distribution

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

[0052]

[0053] In formula (3), z k is the observation vector at time k, x k is the state quantity at time k, τ k is the mixing vector at time k, R k is the measurement noise variance at time k, p(z k |x k ,τ k ,R k ) is x k ,τ k ,R k Under the constraint z k The joint posterior probability density function, τ k,1 , τ k,2 , τ k,3are the first, second, and third elements of the mixed vector at time k, τ k,1 , τ k,2 , τ k,3 represents the weights of Gaussian distribution, Student t distribution, and multivariate K distribution, N(;), St(;) and MK(;) represent the probability density functions of Gaussian, Student t, and multivariate K, respectively, and h(x k ) 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 degree of freedom (dof) parameters of the Student t distribution and the multivariate K distribution, respectively. The mixing vector τ at time k k The probability density function p(τ k ) follows a Dirichlet distribution:

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

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

[0056] The elliptical distribution is a natural extension of the Gaussian distribution, which includes a variety of non-Gaussian distributions with heavy-tail characteristics. It can be expressed in Gaussian form by introducing density functions and scale functions. The following table gives the corresponding density functions and scale functions of Student's t and multivariate K.

[0057] Table 1 Density functions and scale functions of Student's t and multivariate K distributions

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

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

[0060]

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

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

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

[0064] (3) Rewrite the joint posterior probability density function into a Gaussian hierarchical form, such as Figure 2 As shown, we introduce a vector λ=[λ1,λ2…λ N ], λ is a unit vector with only one element being 1 and all other elements being 0, λ1,λ2…λ N is the 1st, 2nd, ..., Nth element.

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

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

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

[0068]

[0069] In formula (8), λ k is the vector that follows the classification distribution at time k, λ k,1 , k,2 , k,3 is k The first, second, and third elements of p(λ k |τ k ) is τ k Under the constraint k The probability density function of . It can be seen that through λ k,1 , k,2 , k,3 The weighted sum form of GSK mixed distribution in formula (3) is transformed into the product form of GSK mixed distribution in formula (8).

[0070] (4) Approximating the posterior estimation of relative motion state quantities and parameters through variational Bayesian techniques:

[0071] Using the free factor distribution to approximate the joint distribution, rewrite the joint posterior PDF:

[0072]

[0073] In formula (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 is the observation vector from time 1 to time k-1, p(x k |z 1:k-1 ) is z 1:k-1 Under the constraint x k The probability density function of

[0074] Minimizing the Courbet–Léibler divergence (KLD) cost function yields a suboptimal approximation to the joint posterior probability density function:

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

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

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

[0078]

[0079] In formula (11), is a constant, ζ k is the set Δ k The elements in , log(;) represents the natural logarithm function, and Ε[;] represents the expectation.

[0080] Step 3: Input the state equation and measurement equation in the relative navigation system into the navigation filter in 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 relative navigation of the aircraft.

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

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

[0083]

[0084] in, is the estimate of the state quantity at time k-1, P k-1|k-1 represents the state variance at time k-1, z k represents the measurement at time k, Q k-1 represents the process noise variance at time k-1, R k represents the measurement noise variance at time k-1, u k|k-1 represents the prior value of the degree of freedom parameter of the inverse Wishart distribution at time k, U k|k-1 represents the inverse scale prior matrix of the inverse Wishart distribution at time k, is the forgetting factor, v1, v2 and a1, a2 represent the shape parameters of the gamma distribution and inverse gamma distribution respectively.

[0085] (2) Time update:

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

[0087] Calculate the state volume point:

[0088]

[0089] In formula (12) and formula (13), is the lth volume point of the state estimation at time k-1, chol() represents Cholesky decomposition, n represents the dimension of the state vector, and e l Represents the lth column vector of the n-dimensional identity matrix.

[0090] State one-step prediction:

[0091]

[0092] In formula (14), is the one-step prediction of the state quantity at time k-1, and f(;) represents the state transfer function.

[0093] Variance one-step forecast:

[0094]

[0095] In formula (15), P k|k-1 represents the 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 one-step prediction

[0099]

[0100] In formula (16), It is the lth volume point of the state quantity predicted in one step at time k-1.

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

[0102]

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

[0104] According to the variational Bayesian inference process from equation (9) to equation (11), we get k , τ k ,ω k , δ k The update formula is as follows:

[0105]

[0106] In formula (18), ψ(;) represents the degamma function, K a (;), K a+1 (;) represent the modified Bessel function of the a-order moment and the modified Bessel function of the a+1-order moment, τ k,j is the jth element of the mixing vector at time k, α k|k,j is the jth element of the concentration parameter matrix at time k, and N is the number of elements in 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 formula (20), m represents the number of observations, and tr(;) represents the trace of the matrix.

[0110]

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

[0112]

[0113] Get the measurement noise variance update:

[0114]

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

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

[0117]

[0118] In formula (24), is the measurement noise variance update at time k.

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

[0120]

[0121] In formula (25), is the one-step prediction of the observation vector at time k-1, It is the lth volume point of the one-step prediction of the state quantity at time k-1.

[0122] Calculate the one-step forecast measurement variance matrix and the covariance matrix

[0123]

[0124]

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

[0126]

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

[0128]

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

[0130] parameter Numeric Aircraft initial state (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 Measuring noise 18arcsec Process noise [10m,0.1m / s]

[0131] Non-stationary heavy-tailed measurement noise v k Produces the following:

[0132]

[0133] In formula (31), wp means “with probability”, N represents 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] The simulation results are as follows Figure 3(a) to Figure 3(f) As shown, the root mean square error of the three filters in estimating the six state quantities of the aircraft in the case of non-Gaussian non-stationary noise is plotted. The six state quantities represent the flight trajectory parameters and uniquely determine the relative position and speed of the aircraft and the target. It can be seen from the figure that the robust relative navigation method based on mixed distribution proposed in the present 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 three distributions in the GSK mixed distribution under non-stationary heavy-tailed noise are plotted. It can be seen from the figure that the relative navigation method proposed in the present invention can adaptively adjust the weights of the mixed distribution to cope with the time-varying noise environment.

[0135] In summary, the present invention proposes a method for robust relative navigation of aircraft based on mixed distribution for non-stationary heavy-tailed noise caused by aircraft operation in complex time-varying environments. The effective application and implementation of this technology has important theoretical and practical significance in reducing the demand for platform power consumption, adapting to more complex and changeable 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: The steps include: Step 1, introducing a Dirichlet random mixture vector representing a Gaussian distribution weight, a Student t distribution weight, and a multivariate K distribution weight, constructing a joint posterior probability density function of the state quantity, the mixture vector, and the observation vector under the constraint of the measurement noise variance, wherein the joint posterior probability density function obeys a Gaussian-Student t-multivariate K mixture distribution, and the statistical prior probability density function of the measurement noise variance is modeled as an inverse Wishart distribution, and the joint posterior probability density function is rewritten into a Gaussian hierarchical form; Step 2: Obtain the state quantity and state variance of the previous moment, the measurement and noise variance of the current moment, and the mixed distribution parameters; Step three, using the cubature Kalman filter algorithm to perform time update and measurement update. In the measurement update process, the measurement noise variance statistics prior probability density function of the inverse Wishart distribution and the joint posterior probability density function of the Gaussian stratified form are used, and the posterior estimates of the relative operating state quantity and the mixed distribution parameters are approximated through variational Bayesian inference.

2. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 1, characterized in that: The joint posterior probability density function of the state quantity, mixed vector, and observation vector under the constraint of measurement noise variance constructed in step 1 is: Among them, z k is the observation vector at time k, x k is the state quantity at time k, τ k is the mixing vector at time k, R k is the measurement noise variance at time k, p(z k |x k ,τ k ,R k ) is x k ,τ k ,R k Under the constraint z k The joint posterior probability density function, τ k,1 , τ k,2 , τ k,3 is the weight of Gaussian distribution, Student t distribution and multivariate K distribution, N(;), St(;) and MK(;) represent the probability density functions of Gaussian, Student t and multivariate K respectively, h(x k ) represents the mean vector composed of the average values ​​of the observation vectors at time k, and v and a represent the degree of freedom parameters of Student’s t distribution and multivariate K distribution, respectively.

3. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 2, characterized in that: The first step is to model the statistical prior probability density function of the measurement noise variance as an inverse Wishart distribution, specifically: p(R k )=IW(R k ;u k|k-1 ,U k|k-1 ), where p(R k ) is the statistical prior probability density function of the measurement noise variance at time k, IW(;) is the probability density function of the inverse Wishart distribution, u k|k-1 represents the one-step forecast value of the degree of freedom parameter of the inverse Wishart distribution at time k-1, U k|k-1 represents the one-step prediction of the inverse scale matrix of the inverse Wishart distribution at time k-1, For the forgetting factor, u k-1|k-1 represents the degree of freedom parameter of the inverse Wishart distribution at time k-1, U k-1|k-1 Represents the inverse scaling matrix of the k-1 inverse Wishart distribution.

4. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 3, characterized in that: The step 1 rewrites the joint posterior probability density function into a Gaussian hierarchical form, specifically: Among them, λ k is the vector that follows the classification distribution at time k, λ k,1 , k,2 , k,3 is k The first, second, and third elements of k,1 , k,2 , k,3 It is composed of a unit vector with only one element being 1 and all other elements being 0, ω k is the density generator for generating the scaling function of the Student's t distribution at time k, δ k is the density generator for generating the multivariate K distribution scaling function at time k, G(;) and IG(;) represent the probability density function of the gamma distribution and the inverse gamma distribution, v1, v2 and a1, a2 represent the shape parameters of the gamma distribution and the inverse gamma distribution, respectively, p(λ k |τ k ) is τ k Under the constraint k The probability density function of the distribution of the categorical distribution, Cat(;) is the probability density function of the categorical distribution, p(τ k ) is the mixing vector τ at time k k The probability density function of Dir(;) is the probability density function of Dirichlet distribution, α k|k-1 is the one-step prediction of the concentration parameter matrix at time k-1, and N is the number of mixing elements.

5. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 4, characterized in that: The mixed distribution parameters obtained in step 2 include: the one-step predicted value u of the degree of freedom parameter of the inverse Wishart distribution at time k-1 k|k-1 , one-step prediction of the inverse scale matrix of the inverse Wishart distribution at time k-1 U k|k-1 , Forgetting Factor Gamma distribution with shape parameters v1, v2 and inverse gamma distribution with shape parameters a1, a2.

6. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 5, characterized in that: In the step three: The cubature Kalman filter algorithm is used for time update, specifically: Calculate the state volume point: One-step prediction of state quantity: Variance one-step forecast: The volumetric Kalman filter algorithm is used in combination with variational Bayesian reasoning for measurement update, specifically: Calculate the volumetric point of the state quantity one-step prediction: Auxiliary matrix A at time k k renew: The variance of the measurement noise at time k renew: R k =U k|k / u k|k , One-step prediction of observation vector: One-step prediction of the measurement variance matrix and measurement covariance matrix: Calculate the filter gain K at time k k : Update the state quantity at time k and the state variance P at time k k|k : in, is the lth volume point of the state quantity at time k-1, chol() represents Cholesky decomposition, n represents the dimension of the state vector, e l represents the lth column vector of the n-dimensional identity matrix, is the one-step prediction of the state quantity at time k-1, f(;) represents the state transfer function, P k|k-1 is the one-step prediction of the state variance at time k-1, is the lth volume point of the one-step predicted state quantity at time k-1, is the estimated value of the state vector at time k, h(;) represents the nonlinear measurement function, E[λ k,1 ], E[λ k,2 ]、E[λ k,3 ]、E[ω k ], E[δ k ] is k,1 , k,2 , k,3 ,ω k , δ k The expectation of , tr(;) represents the trace of the matrix, is the one-step prediction of the observation vector at time k-1, is the one-step forecast of the k-1 moment measurement variance matrix, is the one-step forecast of the k-1-time measurement covariance matrix.

7. The method for robust relative navigation of aircraft based on mixed distribution under non-Gaussian noise according to claim 6, characterized in that: E[λ k ]、E[ω k ], E[δ k ], E[logτ k,j ]According to the measurement noise variance statistical prior probability density function of the inverse Wishart distribution and the joint posterior probability density function of the Gaussian stratified form and approximated by variational Bayesian inference, specifically: 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 is the observation vector from time 1 to time k, z 1:k-1 is the observation vector from time 1 to time k-1, p(x k |Z 1:k-1 ) is z 1:k-1 Under the constraint x k The probability density function of By minimizing the Courbet-Leibler divergence cost function, we obtain a suboptimal approximation of the joint posterior probability density function: q(Δ k )=argminKLD{p(Δ k ,z 1:k )||q(Δ k ,z 1:k )}, q(;) represents the approximate posterior probability distribution of p(;); Use VB method to obtain the general solution of mixed distribution parameters: ζ k For the set Δ k The elements in , log(;) is the natural logarithm function, is a constant, ψ(;) is the degamma function, K a (;), K a+1 (;) represent the modified Bessel function of the a-order moment and the modified Bessel function of the a+1-order moment, τ k,j is the jth element of the mixing vector at time k, α k|k,j is the jth element of the concentration parameter matrix at time k, α k|k =α k|k-1 +Ε[λ k ], 8. 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 according to claim 1 when running the computer program.

9. A computer-readable storage medium having a computer program stored thereon, wherein the computer program executes the steps of the relative navigation method according to claim 1 when executed.

Images