Unmanned aerial vehicle composite non-Gaussian attitude estimation method for dynamic instruction packet loss
By establishing a combined drone malicious/polarization pose system model with embedded dynamics and designing a composite non-Gaussian attitude estimation method, the problem of insufficient accuracy and robustness under composite non-Gaussian interference is solved, real-time estimation of the posterior distribution of the drone attitude and high-precision autonomous navigation in complex environments is realized.
Patent Information
- Application Number
- CN202411696434.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-11-25
- Publication Date
- 2025-03-04
AI Technical Summary
The existing particle filtering algorithms lack accuracy and robustness in composite non-Gaussian interference situations, making it difficult to deal with the uncertainty of state transfer distribution caused by gust interference and unknown dynamics command packet loss probability.
Establish a combined attitude system model of the inertia/polarization of UAV with embedded dynamics, and carefully characterize the composite non-Gaussian interference such as dynamic command packet loss, gust interference, inertial device error and polarization measurement noise, establish a hierarchical probability model of packet loss indicator factors, and design a composite non-Gaussian attitude estimation method to realize real-time estimation of the posterior distribution of UAV.
Improve the identification and compensation ability of the UAV inertia/polarization combined attitude system for composite non-Gaussian interference, and improve the accuracy and reliability of UAV attitude resolution under complex unknown environments and intermittent communication conditions.
Smart Images

Figure CN119555078A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to a composite non-Gaussian attitude estimation method for unmanned aerial vehicles for dynamics command packet loss. For an unmanned aerial vehicle inertial / polarization combined attitude determination system with embedded dynamics, the influence of dynamics command packet loss, gust interference, inertial device error and polarization measurement noise is considered. First, a model of an unmanned aerial vehicle inertial / polarization combined attitude determination system with embedded dynamics is established, and composite non-Gaussian interferences such as dynamics command packet loss, gust interference, inertial device error and polarization measurement noise are finely characterized, and a hierarchical probability model of packet loss indicator factors is established; secondly, a composite non-Gaussian attitude estimation method containing gust interference compensation and packet loss indicator factor posterior distribution tracking is designed to achieve real-time estimation of the posterior distribution of the unmanned aerial vehicle attitude under composite non-Gaussian interference. The present invention can improve the recognition and compensation capabilities of the unmanned aerial vehicle inertial / polarization combined attitude determination system for composite non-Gaussian interference, improve the accuracy and reliability of the unmanned aerial vehicle attitude solution under complex unknown environments and intermittent communication conditions, and can be applied to outdoor autonomous navigation tasks of various types of rotorcraft unmanned aerial vehicles. Background Art
[0002] For the UAV inertial / polarization combined attitude determination system, particle filtering is currently the most widely used type of algorithm. It does not rely on the Gaussian assumption and can theoretically achieve arbitrarily high accuracy. However, most existing particle filtering algorithms rely on precisely known model information, making it difficult to handle the uncertainty of state transition distribution caused by gust interference and the probability of unknown dynamic command packet loss. Therefore, it is necessary to improve the adaptability of the UAV inertial / polarization combined attitude determination system to gust interference and the probability of unknown dynamic command packet loss by improving the existing particle filtering algorithm.
[0003] The above-mentioned UAV inertial / polarization combined attitude determination system filtering problem can be attributed to the design of an improved particle filter algorithm under dynamic interference and unknown parameters. At present, there are few related research results, and the existing results mainly include robust particle filtering and parameter adaptive particle filtering. Among them, the robust particle filter can only process interference signals with dynamic characteristics, and cannot adaptively identify model parameters (such as the probability of packet loss of unknown dynamic instructions); the parameter adaptive particle filter can update the unknown model parameters in real time, but the estimation performance is difficult to guarantee when facing dynamic interference signals (such as gust interference). Summary of the invention
[0004] The problem to be solved by the present invention is: considering the unknown dynamic instruction packet loss probability, gust interference, inertial device error and polarization measurement noise and other composite non-Gaussian interferences in the UAV inertial / polarization combined attitude determination system, overcoming the problem of insufficient accuracy and robustness of the existing particle filter algorithm under the composite non-Gaussian interference situation, and proposing a composite non-Gaussian filter design method for UAVs with dynamic instruction packet loss. The method is based on the UAV inertial / polarization combined attitude determination system model with embedded dynamics, and finely characterizes composite non-Gaussian interferences such as dynamic instruction packet loss, gust interference, inertial device error and polarization measurement noise, establishes a hierarchical probability model of packet loss indicator factors, and designs a composite non-Gaussian attitude estimation method containing gust interference compensation and packet loss indicator factor posterior distribution tracking, so as to realize real-time estimation of the UAV attitude posterior distribution under the composite non-Gaussian interference situation.
[0005] The technical solution of the present invention is: a composite non-Gaussian attitude estimation method for unmanned aerial vehicles for dynamic instruction packet loss, and its implementation steps are as follows:
[0006] Step 1: Establish a UAV inertial / polarization combined attitude determination system model with embedded dynamics, and finely characterize the complex non-Gaussian interference such as dynamic command packet loss, gust interference, inertial device error and polarization measurement noise, and establish a hierarchical probability model of packet loss indication factors:
[0007] First, the state equation is established as follows:
[0008]
[0009] Among them, φ(t), θ(t), ψ(t) are the roll angle, pitch angle and yaw angle of the UAV in the navigation coordinate system "n-system" respectively; p(t), q(t), r(t) are the components of the angular velocity of the UAV in the body coordinate system "b-system" on the x-axis, y-axis and z-axis respectively; M x (t), M y (t), M z (t) are the components of the UAV torque on the x-axis, y-axis and z-axis respectively; I xx ,I yy ,I zz is the moment of inertia of the drone, usually given in the drone manual; w p (t), w q (t), w r (t) are all Gaussian noise items with zero mean and known variance; γ(t) is the packet loss indicator factor, which takes the value of 1 or 0, where 1 represents that the dynamics command is normal and 0 represents that the dynamics command has packet loss; d p (t), d q (t), d r(t) are the components of gust disturbance in the x-axis, y-axis and z-axis of the “b-system” respectively, and d(t) = [d p (t)d q (t)d r (t)] T , and its dynamic characteristics are described as the following first-order Markov process:
[0010]
[0011] in, is the state variable of the gust interference subsystem, V=diag{[a p 0],[a q 0],[a r 0]}, and σ p , σ q , σ r is the angular frequency of the harmonic interference, a p 、a q 、a r is the amplitude of harmonic interference; n(t) is a Gaussian noise term with zero mean and known variance;
[0012] Secondly, the gyroscope and accelerometer measurement equations are established as follows:
[0013]
[0014] in, are the x-axis, y-axis and z-axis angular velocity components measured by the gyroscope respectively; are the x-axis and y-axis acceleration components measured by the accelerometer; g is the local gravity acceleration; v p (t), v q (t), v r (t), v x (t), v y (t) is the inertial navigation sensor measurement noise, which obeys a Gaussian distribution with zero mean and known variance;
[0015] Next, the polarization sensor measurement equation is established as follows:
[0016]
[0017] in, is the real-time measurement value of the polarization sensor; v pol (t) is the polarization sensor measurement noise, which obeys Gaussian distribution; K I is the known polarization sensor light intensity gain coefficient; K d is the known optical path coupling coefficient of the polarization sensor; A is the known solar direction angle; is the known polarization detection angle; I in is the known incident light intensity coefficient;
[0018] Furthermore, by combining equations (1), (2), (3), and (4), we can obtain:
[0019]
[0020] in, x(t)=[φ(t)θ(t)ψ(t)p(t)q(t)r(t)] T , u(t)=[000M x (t)M y (t)M z (t)] T , w(t)=[000w p (t)w q (t)w r (t)] T , v(t)=[v p (t)v q (t)v r (t)a x (t)a y (t)I out (t)] T ; and They are respectively the nonlinear system model and measurement equation consisting of the terms in formula (1) and formula (3) that are only related to x(t);
[0021] Finally, the model described by equation (5) is discretized to obtain the discrete state space model of the UAV inertial / polarization combined attitude determination system with embedded dynamics as follows:
[0022]
[0023] Among them, W k =I+WT s , V k =V; the nonlinear system model and measurement equations f(·) and h(·) can be respectively calculated according to and And the sampling time interval T s Calculated; w k 、n k and v k is a Gaussian random noise term, satisfying: w k ~N(0,Σ w,k ), n k ~N(0,Σ n,k ), and v k ~N(0,Σv,k ), and Σ n,k ,Σ w,k and Σ v,k are all known; γ k is the packet loss indicator variable at the kth moment, and its probability distribution obeys the following hierarchical probability model:
[0024] γ k ~Bern(γ k ; 1-τ k ),τ k ~Beta(τ k ; α k ,β k ) (7)
[0025] Where Bern(p) represents the Bernoulli distribution with success probability p; τ k is the unknown packet loss probability, subject to the parameter α k With β k Beta distribution of
[0026] Step 2: For the case where the system structure parameters are known, a composite non-Gaussian attitude estimation method with gust interference compensation and packet loss indicator factor posterior distribution tracking is designed to achieve real-time estimation of the posterior distribution of the UAV attitude under composite non-Gaussian interference:
[0027] At the kth moment, firstly, according to the particle set obtained at the k-1th moment and the parameter estimates corresponding to each particle The parameter prediction value corresponding to each particle at the kth moment is calculated as follows:
[0028]
[0029] Where ρ is a pre-specified forgetting factor, satisfying 0<ρ<1;
[0030] Next, using the above parameter prediction values and the distribution type of each parameter distribution, sample the parameters corresponding to each particle at the kth moment:
[0031]
[0032] Where N is the number of particles specified in advance; N(·; μ, Σ) represents a Gaussian distribution with mean μ and variance Σ; E[·] represents the expectation of a random variable; Then, the particle set calculated by the filter at the k-1th moment is combined And the system model (6), generating particle set
[0033]
[0034] Secondly, for each particle definition E <0> [γ k ]=E i- [γ k ], E <0) [log(1-τ k )]=E i- [log(1-τ k )],
[0035] E (0) [logτ k ]=E i- [logτ k ], where E i- [γ k ]、E i- [log(1-τ k )] and E i- [logτ k ] is calculated by the following formula:
[0036]
[0037] where ψ(·) represents the inverse Gamma function; then, for j = 1, 2, ..., N max , the loop performs the following calculation steps:
[0038]
[0039] Among them, N max is the pre-specified maximum number of iterations; is the virtual quantity measurement corresponding to the i-th particle; <j>< / j> To ensure Pr <j>< / j> (γ k =1)+Pr <j>< / j> (γ k =0)=1; after the cycle ends, the parameter estimate is updated to
[0040]
[0041] Next, for i=1,2,...,N, calculate the particle The corresponding weights are as follows:
[0042]
[0043] Among them, w k,i is the particle weight; according to the weight, a resampling sequence {r1,r2,...,r N}, where ri Represents the sequence number of the i-th particle after resampling, satisfying Where l(x) is the indicator function, satisfying According to the resampling sequence, the particles and their associated parameters are updated as follows:
[0044]
[0045] Finally, output the particle set Packet loss indicator variable estimate and parameter estimates and As the input variable of the filter at the k+1th moment.
[0046] The advantages of the present invention compared with the prior art are:
[0047] (1) The present invention takes the bionic autonomous integrated navigation system of unmanned aerial vehicles as the research object, introduces accelerometers, gyroscopes and polarization sensors to realize the observation of the Euler angles and angular velocities of the unmanned aerial vehicles, and considers the dynamics command packet loss, gust interference, inertial device error and polarization measurement noise that have a great impact on the attitude dynamics stability of the quad-rotor unmanned aerial vehicle. Based on physical reality, a model of the unmanned aerial vehicle inertial / polarization combined attitude determination system with embedded dynamics is established, and the composite non-Gaussian interferences such as dynamics command packet loss, gust interference, inertial device error and polarization measurement noise are finely characterized, and a hierarchical probability model of packet loss indication factor is established;
[0048] (2) The present invention adopts a composite non-Gaussian filter to simultaneously realize gust interference compensation, online tracking of the probability distribution of the dynamic command loss indicator factor, and real-time estimation of the posterior distribution of the quadrotor UAV attitude, thereby overcoming the problems of insufficient estimation accuracy and robustness of the existing robust inertial / polarization combined attitude determination filtering algorithm under the conditions of dynamic interference (gust interference) and unknown random parameters (dynamic command loss indicator factor). BRIEF DESCRIPTION OF THE DRAWINGS
[0049] Figure 1 It is a design flow chart of the present invention. DETAILED DESCRIPTION
[0050] like Figure 1 As shown, the specific implementation steps of the present invention are as follows:
[0051] Step 1: Establish a UAV inertial / polarization combined attitude determination system model with embedded dynamics, and finely characterize the complex non-Gaussian interference such as dynamic command packet loss, gust interference, inertial device error and polarization measurement noise, and establish a hierarchical probability model of packet loss indication factors:
[0052] First, the state equation is established as follows:
[0053]
[0054] Among them, φ(t), θ(t), ψ(t) are the roll angle, pitch angle and yaw angle of the UAV in the navigation coordinate system "n-system" respectively; p(t), q(t), r(t) are the components of the angular velocity of the UAV in the body coordinate system "b-system" on the x-axis, y-axis and z-axis respectively; M x (t), M y (t), M z (t) are the components of the UAV torque on the x-axis, y-axis and z-axis respectively; I xx ,I yy ,I zz is the moment of inertia of the drone, usually given in the drone manual; w p (t), w q (t), w r (t) are all Gaussian noise items with zero mean and known variance; γ(t) is the packet loss indicator factor, which takes the value of 1 or 0, where 1 represents that the dynamics command is normal and 0 represents that the dynamics command has packet loss; d p (t), d q (t), d r (t) are the components of gust disturbance in the x-axis, y-axis and z-axis of the “b-system” respectively, and d(t) = [d p (t)d q (t)d r (t)] T , and its dynamic characteristics are described as the following first-order Markov process:
[0055]
[0056] in, is the state variable of the gust interference subsystem, V=diag{[a p 0],[a q 0],[a r 0]}, and σ p , σ q , σ r is the angular frequency of the harmonic interference, a p 、a q 、a r is the amplitude of harmonic interference; n(t) is a Gaussian noise term with zero mean and known variance;
[0057] Secondly, the gyroscope and accelerometer measurement equations are established as follows:
[0058]
[0059] in, are the x-axis, y-axis and z-axis angular velocity components measured by the gyroscope respectively; are the x-axis and y-axis acceleration components measured by the accelerometer; g is the local gravity acceleration; v p (t), v q (t), v r (t), v x (t), v y (t) is the inertial navigation sensor measurement noise, which obeys a Gaussian distribution with zero mean and known variance;
[0060] Next, the polarization sensor measurement equation is established as follows:
[0061]
[0062] in, is the real-time measurement value of the polarization sensor; v pol (t) is the polarization sensor measurement noise, which obeys Gaussian distribution; K I is the known polarization sensor light intensity gain coefficient; K d is the known optical path coupling coefficient of the polarization sensor; A is the known solar direction angle; is the known polarization detection angle; I in is the known incident light intensity coefficient;
[0063] Furthermore, by combining equations (1), (2), (3), and (4), we can obtain:
[0064]
[0065] in, x(t)=[φ(t)θ(t)ψ(t)p(t)q(t)r(t)] T , u(t)=[000M x (t)M y (t)M z (t)] T , w(t)=[000w p (t)w q (t)w r (t)] T , v(t)=[v p (t)v q (t)v r (t)a x (t)a y (t)I out (t)] T ; and They are respectively the nonlinear system model and measurement equation consisting of the terms in formula (1) and formula (3) that are only related to x(t).
[0066] Finally, the model described by equation (5) is discretized to obtain the inertial / polarization combined attitude-fixed quadrotor UAV system model as follows:
[0067]
[0068] Among them, W k =I+WT s , V k =V; the nonlinear system model and measurement equations f(·) and h(·) can be respectively calculated according to and And the sampling time interval T s Calculated; w k 、n k and v k is a Gaussian random noise term, satisfying: w k ~N(0,Σ w,k ), n k ~N(0,Σ n,k ), and v k ~N(0,Σ v,k ), and Σ n,k ,Σ w,k and Σ v,k are all known; γ k is the packet loss indicator variable at the kth moment, and its probability distribution obeys the following hierarchical probability model:
[0069] γ k ~Bern(γ k ; 1-τ k ),τ k ~Beta(τ k ; α k ,β k ) (7) Where Bern(p) represents the Bernoulli distribution with success probability p; τ k is the unknown packet loss probability, subject to the parameter α k With β k Beta distribution of
[0070] Step 2: For the case where the system structure parameters are known, a composite non-Gaussian filter with gust interference compensation and online identification of control signal packet loss probability is designed:
[0071] At the kth moment, firstly, according to the particle set obtained at the k-1th moment and the parameter estimates corresponding to each particle The parameter prediction value corresponding to each particle at the kth moment is calculated as follows:
[0072]
[0073] Where ρ is a pre-specified forgetting factor, satisfying 0<ρ<1;
[0074] Next, using the above parameter prediction values and the distribution type of each parameter distribution, sample the parameters corresponding to each particle at the kth moment:
[0075]
[0076] Where N is the number of particles specified in advance; N(·; μ, Σ) represents a Gaussian distribution with mean μ and variance Σ; E[·] represents the expectation of a random variable; Then, the particle set calculated by the filter at the k-1th moment is combined And the system model (6), generating particle set
[0077]
[0078] Secondly, for each particle definition E <0> [γ k ]=E i- [γ k ], E <0> [log(1-τ k )]=E i- [log(1-τ k )],E <0> [logτ k ]=E i- [logτ k ], where E i- [γ k ]、E i- [log(1-τ k )] and E i- [logτ k ] is calculated by the following formula:
[0079]
[0080] where ψ(·) represents the inverse Gamma function; then, for j = 1, 2, ..., N max , the loop performs the following calculation steps:
[0081]
[0082] Among them, Nmax is the pre-specified maximum number of iterations; is the virtual quantity measurement corresponding to the i-th particle; <j>< / j> To ensure Pr <j>< / j> (γ k =1)+Pr <j>< / j> (γ k =0)=1; after the cycle ends, the parameter estimate is updated to
[0083]
[0084] Next, for i=1,2,...,N, calculate the particle The corresponding weights are as follows:
[0085]
[0086] Among them, w k,i is the particle weight; according to the weight, a resampling sequence {r1,r2,...,r N}, where r i Represents the sequence number of the i-th particle after resampling, satisfying Where l(x) is the indicator function, satisfying According to the resampling sequence, the particles and their associated parameters are updated as follows:
[0087]
[0088] Finally, output the particle set Packet loss indicator variable estimate and parameter estimates and As the input variable of the filter at the k+1th moment.
[0089] The contents not described in detail in the specification of the present invention belong to the prior art known to the professional and technical personnel in this field.
Claims
1. A composite non-Gaussian attitude estimation method for unmanned aerial vehicles for dynamic command packet loss, characterized in that: The implementation steps are as follows: Step 1: Establish a UAV inertial / polarization combined attitude determination system model with embedded dynamics, and finely characterize the complex non-Gaussian interference such as dynamic command packet loss, gust interference, inertial device error and polarization measurement noise, and establish a hierarchical probability model of packet loss indication factors; Step 2: Design a composite non-Gaussian attitude estimation method that includes gust interference compensation and packet loss indicator factor posterior distribution tracking to achieve real-time estimation of the UAV attitude posterior distribution under composite non-Gaussian interference.
2. The method for estimating composite non-Gaussian attitude of a UAV for dynamic instruction packet loss according to claim 1 is characterized in that: The specific implementation of the first step is as follows: First, the UAV dynamics equation is established as follows: Among them, φ(t), θ(t), ψ(t) are the roll angle, pitch angle and yaw angle of the UAV in the navigation coordinate system "n-system" respectively; p(t), q(t), r(t) are the components of the angular velocity of the UAV in the body coordinate system "b-system" on the x-axis, y-axis and z-axis respectively; M x (t), M y (t), M z (t) are the components of the UAV torque on the x-axis, y-axis and z-axis respectively; I xx ,I yy ,I zz is the moment of inertia of the drone, usually given in the drone manual; w p (t), w q (t), w r (t) are all Gaussian noise items with zero mean and known variance; γ(t) is the packet loss indicator factor, which takes the value of 1 or 0, where 1 represents that the dynamics command is normal and 0 represents that the dynamics command has packet loss; d p (t), d q (t), d r (t) are the components of gust disturbance in the x-axis, y-axis and z-axis of the "b-system", respectively. Define d(t) = [d p (t)d q (t)d r (t)] T , and its dynamic characteristics are described as the following first-order Markov process: in, is the state variable of the gust interference subsystem, V = diag{[a p 0],[a q 0],[a r 0]} and σ p , σ q , σ r is the angular frequency of harmonic interference, a p 、a q 、a r is the amplitude of harmonic interference; n(t) is a Gaussian noise term with zero mean and known variance; Secondly, the gyroscope and accelerometer measurement equations are established as follows: in, are the x-axis, y-axis and z-axis angular velocity components of the "b-frame" measured by the gyroscope respectively; are the x-axis and y-axis acceleration components measured by the accelerometer; g is the local gravity acceleration; v p (t), v q (t), v r (t), v x (t), v y (t) is the inertial navigation sensor measurement noise, which obeys a Gaussian distribution with zero mean and known variance; Next, the polarization sensor measurement equation is established as follows: in, is the real-time measurement value of the polarization sensor; v pol (t) is the polarization sensor measurement noise, which obeys Gaussian distribution; K I is the known polarization sensor light intensity gain coefficient; K d is the known optical path coupling coefficient of the polarization sensor; A is the known solar direction angle; is the known polarization detection angle; I in is the known incident light intensity coefficient; Furthermore, by combining equations (1), (2), (3), and (4), we can obtain: in, x(t)=[φ(t)θ(t)ψ(t)p(t)q(t)r(t)] T , u(t)=[000M x (t)M y (t)M z (t)] T , w(t)=[000w p (t)w q (t)w r (t)] T , v(t)=[v p (t)v q (t)v r (t)a x (t)a y (t)I out (t)] T ; and They are respectively the nonlinear system model and measurement equation consisting of the terms in formula (1) and formula (3) that are only related to x(t); Finally, the continuous system model described by equation (5) is discretized to obtain the discrete state space model of the UAV inertial / polarization combined attitude determination system with embedded dynamics as follows: Among them, W k =I+WT s , V k =V; the nonlinear system model and measurement equations f(·) and h(·) can be respectively calculated according to and And the sampling time interval T s Calculated; w k 、n k and v k is a Gaussian random noise term, satisfying: w k ~N(0,Σ w,k ), n k ~N(0,Σ n,k ), and v k ~N(0,Σ v,k ), and Σ n,k ,Σ w,k and Σ v,k are all known; γ k is the packet loss indicator variable at the kth moment, and its probability distribution obeys the following hierarchical probability model: γ k ~Bern(γ k ; 1-τ k ),τ k ~Beta(τ k ; α k ,β k )(7)Wherein, Bern(p) represents the Bernoulli distribution with success probability p; τ k is the unknown packet loss probability, subject to the parameter α k With β k Beta distribution.
3. The method for estimating composite non-Gaussian attitude of a UAV for dynamic instruction packet loss according to claim 1, characterized in that: The specific implementation of the second step is as follows: At the kth moment, firstly, according to the particle set obtained at the k-1th moment and the parameter estimates corresponding to each particle The parameter prediction value corresponding to each particle at the kth moment is calculated as follows: Where ρ is a pre-specified forgetting factor, satisfying 0<ρ<1; Next, using the above parameter prediction values and the distribution type of each parameter distribution, sample the parameters corresponding to each particle at the kth moment: Where N is the number of particles specified in advance; N(·; μ, Σ) represents a Gaussian distribution with mean μ and variance Σ; E[·] represents the expectation of a random variable; Then, the particle set calculated by the filter at the k-1th moment is combined And the system model (6), generating particle set Secondly, for each particle definition E <0> [γ k ]=E i- [γ k ], E (0) [log(1-τ k )]=E i- [log(1-τ k )],E <0> [logτ k ]=E i- [logτ k ], where E i- [γ k ]、E i- [log(1-τ k )] and E i- [logτ k ] is calculated by the following formula: where ψ(·) represents the inverse Gamma function; then, for j = 1, 2, ..., N max , the loop performs the following calculation steps: Among them, N max is the pre-specified maximum number of iterations; is the virtual quantity measurement corresponding to the i-th particle; <j>< / j> To ensure Pr <j>< / j> (γ k =1)+Pr <j>< / j> (γ k =0)=1; after the cycle ends, the parameter estimate is updated to Next, for i=1,2,...,N, calculate the particle The corresponding weights are as follows: Among them, w k,i is the particle weight; according to the weight, a resampling sequence {r1,r2,...,r N }, where r i Represents the sequence number of the i-th particle after resampling, satisfying Where l(x) is the indicator function, satisfying According to the resampling sequence, the particles and their associated parameters are updated as follows: Finally, output the particle set Packet loss indicator variable estimate and parameter estimates and As the input variable of the filter at the k+1th moment.