A Target Tracking Method Based on Variational Bayesian under Non-Gaussian and Non-Stationary System Noise

Abstract

The present invention proposes a target tracking method based on variational Bayes under non-Gaussian and non-stationary system noise. The method comprises the following steps: first, in a distributed target tracking framework, prior modeling is performed according to the probability density obeyed by the unknown parameters; second, the standard variational Bayes method is used to jointly estimate the posterior distribution parameters such as the target state, mixing coefficient and scale matrix through fixed-point iteration; finally, the local platform state estimation value is fused according to the covariance cross-fusion strategy in a multi-sensor distributed framework. The present invention comprehensively considers the influence of non-Gaussian and non-stationary heavy-tailed process noise and Skew measurement noise in the target tracking process, and can effectively estimate the unknown parameters such as the target state and noise covariance, thereby improving the target tracking accuracy, while having good adaptability and robustness.

Description

Technical Field

[0001] The present invention relates to the technical field of multi-sensor target tracking, and in particular to an application environment in which state process noise is non-Gaussian and non-stationary heavy-tailed process noise and non-Gaussian and non-stationary Skew measurement noise. Background Art

[0002] In complex environments, target tracking systems within a multi-sensor distributed framework can suffer performance degradation due to complex system noise caused by pulse interference in clutter, inaccurate modeling, and measurement outliers from unreliable sensors. This can lead to target loss. This system noise exhibits non-Gaussian and non-stationary characteristics, leading to a continuous increase in tracking noise interference. Therefore, addressing this complex system noise is crucial.

[0003] To address the state estimation problem for state-space models with non-Gaussian and non-stationary heavy-tailed state noise, the process noise is modeled using the Student t distribution. Examples include the Student t hybrid filter proposed by J. Loxam and the various Student t-based outlier robust Kalman filters proposed by J. Ting et al. However, the performance of these robust Kalman filters degrades dramatically for state-space models with heavy-tailed state noise because they are all based on Gaussian state noise models. To address the state estimation problem for state-space models with heavy-tailed state noise, CD Karlgaard proposed the Hong Kong Kalman Filter (HKF) and R. Izanloo proposed the Maximum Entropy Kalman Filter (MCKF). Both the HKF and the MCKF are able to suppress the increase in estimation error caused by heavy-tailed noise, thereby mitigating the negative impact. Unfortunately, the design of the HKF and the MCKF does not exploit the inherent heavy-tailed characteristics of state and measurement noise, resulting in limited estimation accuracy. To achieve better estimation performance, a reasonable approach is to improve the modeling of heavy-tailed non-Gaussian probability density functions (PDFs). M. Roth, Y.L. Huang, and others proposed the Student-t-based filter (STF) and the robust Student-t-based Kalman filter (RSTKF). Most of these filters share a common problem: they cannot be truly applied to non-stationary scenarios and are designed under the assumption that measurement noise follows a heavy-tailed distribution. Therefore, to comprehensively consider the non-stationary, non-Gaussian, heavy-tailed characteristics of state noise, the Heavy-Tailed Mixture (HTM) distribution in the HTM-RKF proposed by Y.L. Huang is used to address the state noise issue.

[0004] Johnson et al. proposed a two-component Gaussian mixture (GM2) model to model skew-t noise, compensating for the degradation of state estimation accuracy caused by non-Gaussian skew measurement noise. However, this model has an exponentially decaying tail and is overly sensitive to measurement outliers. More advanced, Loxam et al. proposed an improved Student's t mixture filter, which simulates skew measurement noise using multiple Student's t mixtures. This improves filtering accuracy, but it is difficult to strike a balance between computational complexity and filtering accuracy. Unfortunately, all of these filters for skew measurement noise are designed based on the assumption that the process noise is Gaussian distributed. When the process noise follows a heavy-tailed distribution, severe performance degradation will occur. To address this issue, Huang et al. proposed a KF based on the Gaussian scale mixture (GSM) distribution (GSMKF) and an outlier-robust fixed-interval smoothing KF, providing a general framework for filtering problems in the presence of heavy-tailed process noise and skew measurement noise. Although all of these filters can compensate for the erosion of estimation accuracy caused by skew measurement noise to varying degrees, Unfortunately, since the system uses fixed parameter information, there are usually large state estimation errors in scenarios where the noise has a Skew-t distribution and non-stationary changes. Therefore, the Gaussian Scaled Mixture Gamma Inverse Wishart (GSMGIW) distribution is proposed here to address the impact of non-Gaussian and non-stationary skew measurement noise.

[0005] Although system noise has been comprehensively considered, in actual target tracking systems, the sensor measurement accuracy of a single platform often fails to meet the system filtering requirements. Designing a multi-platform sensor tracking system based on data fusion mechanisms can effectively address this problem of insufficient filtering accuracy. Distributed fusion is a typical implementation of this concept. It feeds the sensor measurements of each platform into local filters for state estimation, and then sends the local state estimation results to the fusion center for processing, thereby achieving better target state estimation than using single-platform sensor measurement data. Here, a covariance intersection (CI) fusion strategy is adopted to achieve better fusion results. Therefore, distributed fusion is very necessary to solve the problem of target tracking in complex environments with non-Gaussian and non-stationary system noise. Summary of the Invention

[0006] The present invention proposes a target tracking method based on variational Bayes under non-Gaussian and non-stationary system noise in a distributed fusion framework. Considering the physical scenario in which multiple sensors jointly track the same target, the method not only solves the problem that the process noise is non-Gaussian and non-stationary heavy-tailed noise; but also reduces the communication burden of the sensor and improves the target tracking accuracy.

[0007] The present invention solves the technical problem by providing a target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise, comprising the following steps:

[0008] S1: Establish the one-step prediction probability density P(x k,s ∣z 1:k-1 ) model, measure likelihood probability density P(z k,s ∣x k,s ) model and log joint probability density Model;

[0009] S2: Variational Bayes approximate posterior:

[0010] The model parameters constructed in step S1 are mutually coupled, and the posterior of the unknown parameters cannot be solved analytically. The variational Bayesian method is introduced to approximately solve the posterior of the unknown parameters in S1. The solution formula is as follows:

[0011]

[0012] in, Express about Ask for the expected operation. Indicates taking An element of a set, yes Centralized removal The remaining elements, Relatively unknown parameters is an irrelevant constant.

[0013] S3: Use the solution formula of S2 to estimate the shape parameter β k,s and the target state x k,s ;

[0014] S4: Using the solution formula of S2, estimate the state prediction covariance matrix Σ k,s and the measurement noise covariance matrix R k,s ;

[0015] S5: Use the solution formula of S2 to estimate the mixed random vector τ k and species distribution vector ξ k ;

[0016] S6: Use the solution formula of S2 to estimate the auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ;

[0017] S7: Repeat steps S3-S6 until the total number of variational iterations is reached or the iteration termination condition is met: ε is usually set to 10 -10 , the iteration ends; the final iteration result is: and Σ k|k,s , and assign it as the initial value of the variational iteration at time k+1.

[0018] S8: Multi-sensor distributed fusion:

[0019] In a distributed network framework, it is often difficult to achieve the expected accuracy by relying on a single local platform for target tracking after the local filter completes the deviation registration. Therefore, it is necessary to use the state values of each local filter obtained by the variational iteration of each platform in step S7. And its covariance matrix P k|k,s The information is fused using the covariance intersection (CI) fusion method to obtain the optimal target state after fusion. and its covariance matrix

[0020] Compared with the physical scenario of traditional multi-sensor target tracking, the linear system represented by the discrete-time linear state space model is as follows, taking into account the influence of target state process noise and sensor measurement noise.

[0021] x k,s =F k x k-1,s +w k-1

[0022] z k,s =H k x k,s +v k,s

[0023] Among them, x k,s is the target state of sensor s at time k, k = 1, 2, 3..., ts, F k is the state transfer matrix from time k-1 to time k, w k-1 is the k-1 moment process noise, and its noise covariance matrix is Q k-1 , the process noise presents non-Gaussian and non-stationary heavy-tail characteristics; z k,s is the measurement value of sensor s at time k, H k is the measurement transfer matrix at time k; v k,s is the measurement noise of sensor s at time k, and its noise covariance matrix is R k,s , the measurement noise presents non-Gaussian and non-stationary Skew characteristics.

[0024] The variational parameter settings and initialization expectations are as follows:

[0025] (1) Set the initial value of the variational parameter at time k-1 as follows:

[0026] Initial estimated state and its error covariance matrix P k-1|k-1,s ;

[0027] Measurement noise covariance matrix R k-1,s The degree of freedom parameter and the scale matrix is U k-1|k-1,s ;

[0028] Degree of freedom parameter ω k-1 ,ν k-1 Shape parameters Sum rate parameters and

[0029] Mixed random vector τ k The prior concentration parameter α k|k-1 and species distribution vector ξ k The mean vector of

[0030] Nominal (not true) process noise covariance matrix and the nominal measurement noise covariance matrix R k-1,s ;

[0031] Shape parameter β k-1 The mean and covariance matrix of P β,k|k,s ;

[0032] and the measurement value z of sensor s at time k k,s and other related parameters.

[0033] (2) k-time parameter time update:

[0034] State one-step prediction: And its predicted state error covariance matrix:

[0035] P k|k-1,s The degree of freedom parameter and the scale matrix

[0036] Measurement noise covariance moment R k,s The degree of freedom parameter and the scale matrix U k|k-1,s =U k|k-1,s +E[f(γ k )]B k ;

[0037] Degrees of freedom ωk 、ν k Shape parameters Sum rate parameters

[0038] Shape parameter β k,s The covariance matrix P β,k|k,s =σ k (I m -E j [γ k ]K β,k ) and mean

[0039] Mixed random vector τ k Concentration parameters and species distribution vector ξ k Mixed vector in Where E[·] represents the expectation operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter and is usually 1-exp(-4), τ p represents the tuning factor for the forecast error covariance.

[0040] (3) Initialize the expected value at time k:

[0041] Initial expected state and its covariance matrix:

[0042] Initialize the following expected values: shape parameter β k,s Expectations

[0043] Auxiliary variable λ k , γ k Expected E 0 [λ k ]=E 0 [γ k ]=1;

[0044] Degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ];

[0045] Measurement noise covariance R k,s Expectations

[0046] Mixed random vector τ k Expectations and species distribution vector ξ k Expectations

[0047] Superscript E 0 [·] means finding the expectation of the 0th variational iteration.

[0048] In step S1, P(x k,s ∣z 1:k-1 ) and P(z k,s ∣x k,s ) Prior model, logarithmic joint probability density model The model specifically includes the following steps:

[0049] S1.1: Establish the state one-step prediction probability density P(x k,s ∣z 1:k-1 ) Prior model:

[0050] An HTM distribution is proposed to describe the non-Gaussian and non-stationary heavy-tailed characteristics of state process noise. The heavy-tailed and non-stationary noise are modeled by a hybrid method, in which the mixing coefficient is adaptively learned according to the sequential measurement. k,s ∣z 1:k-1 ) can be expressed in Gaussian hierarchical form, introducing the classification distribution vector ξ k , specifically expressed as:

[0051]

[0052] Where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, and G(·; a, b) represents a gamma distribution with shape parameter a, rate parameter b, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter; IW(·; t, T) represents the inverse Wishart distribution, whose degrees of freedom are t and the scaling matrix is T; Cat(·; m, M) represents the species distribution, whose mixing random vector is m and the number of mixing elements is M; Dir(·; θ, M) represents the Dirichlet distribution, whose concentration parameter is θ and the number of mixing elements is M.

[0053] S1.2: Establish the measurement likelihood probability density P(z k,s ∣x k,s ) Prior model:

[0054] Propose GSMGIW distribution and introduce auxiliary variable γ k Written in the form of Gaussian-Gamma mixed distribution, the measurement noise covariance matrix R k,s Establish a priori model, then P(z k,s ∣x k,s ) can be expressed as:

[0055]

[0056] Among them, Hk x k,s and R k,s denote the mean vector and covariance matrix respectively, U k|k-1,s Denote the degrees of freedom and the scale matrix respectively. k >0 is the mixing parameter, ν k represents the degree of freedom, f(γ k ) and IG(γ k ν k / 2,ν k / 2) are expressed as positive scaling function and inverse gamma function respectively. k,s is a shape parameter that controls the degree of skewness. k,s ≠0 is asymmetric, otherwise it is symmetric. In particular, when β k,s =0, and f(γ k )=γ k , GSMGIW degenerates into the Normal-Gamma-inverse Wishart (NGIW) distribution model.

[0057] S1.3: Constructing a joint probability density model

[0058] According to the state one-step probability density P(x k,s ∣z 1:k-1 ) prior model and step S1.2 to establish the measurement likelihood probability density P(z k,s ∣x k,s ) prior model, rewritten as a layered Gaussian, then the logarithmic joint probability density can be expressed as:

[0059]

[0060] in, represents the unknown parameter set, Z represents the measurement data set, represents the logarithm of the joint probability density with respect to a natural constant.

[0061] In step S3, β is estimated k,s and x k,s The details are as follows:

[0062] Separate orders and Apply the variational iteration formula in step S2, and based on the conjugation of the exponential distribution, approximate the posterior q(β k,s ) and q(x k,s ) are updated to conform to the Gaussian distribution:

[0063]

[0064] Similar to the Kalman filter measurement update step to solve the approximate posterior q i+1 (β k,s ) mean: and its corresponding covariance matrix By calculating the expectation separately: E i+1 (β k,s ), E i+1 (P β,k|k,s ).

[0065] Calculated by Kalman filter measurement update step and

[0066]

[0067] In step S4, the estimated Σ k,s and R k,s The details are as follows:

[0068] Separate orders Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, the approximate posterior q(Σ k,s ) and q(R k,s ) are updated to conform to the IW distribution:

[0069]

[0070] Solve the approximate posterior q separately i+1 (Σ k,s ) and q i+1 (R k,s ) degrees of freedom parameter: And the scale matrix parameters: Calculate the expected E from the IW distribution properties i+1 [P k|k-1,s ], E i+1 [R k,s ], calculate the modified state prediction covariance matrix and the measurement noise covariance matrix The state one-step prediction probability density P(x k,s ∣z 1:k-1 ) The prior model shows that the state noise is non-stationary, and here it is designed as a mixture of multiple heavy tails, so P k|k-1,s and Σ k,s For a one-to-many relationship:

[0071]

[0072] Among them, E i [·] means finding the expectation of the i-th variational iteration, i=0:Nm -1.

[0073] In step S5, the estimated value ξ is updated. k and τ k , as shown below:

[0074] S5.1: Update the mixed random vector τ k :

[0075] make Apply the variational iteration formula in step S2, and the approximate posterior q(τ k ) is updated to conform to the Dirichlet distribution:

[0076]

[0077] Solve for the approximate posterior q j+1 (τ k ) concentration parameters: According to the properties of Dirichlet distribution, calculate the natural logarithm expectation E j+1 [τ k ] and E j+1 [logτ k,i ].

[0078] S5.2: Update the species distribution vector ξ k :

[0079] make Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, q(ξ k ) is updated to match the category distribution:

[0080]

[0081] Solve for the approximate posterior q j+1 (ξ k ) of the mixed vector: According to the species distribution properties, calculate the natural logarithm expectation E i+1 [ξ k ].

[0082] In step S6, update and solve λ k , γ k 、ω k 、ν k , as follows:

[0083] S6.1: Update auxiliary variable λ k and γ k :

[0084] make Applying the variational iteration formula in step S2, the approximate posterior conforms to the gamma distribution, that is:

[0085]

[0086] Solve for the approximate posterior q i+1 (λ k )’s shape and rate parameters: Calculate the expectation E from the properties of the gamma distribution i+1 [λ k ]、E i+1 [logλ k ].

[0087] make Apply the variational iteration formula of step S2, and according to the conjugation of the exponential distribution, approximate the posterior q(γ k ) conforms to the generalized inverse Gaussian distribution (GIG):

[0088]

[0089] Solve for the approximate posterior q i+1 (γ k )’s shape parameters: and Calculate the expectation E by the generalized inverse Gaussian distribution property i+1 [γ k ]、E i+1 [logγ k ].

[0090] S6.2: Update the degree of freedom parameter ω k and ν k :

[0091] Separate orders and Apply the variational iteration formula of step S2, and according to the conjugation of the exponential distribution, approximate the posterior q(ω k ) and q(ν k ) is updated to conform to the gamma distribution:

[0092]

[0093] Solve for the approximate posterior q i+1 (ω k ) and q i+1 (ν k )’s shape parameters: and rate parameters: Calculate the expectation E from the properties of the gamma distribution i+1 [ω k ] and E i+1 [v k ].

[0094] The step S8 of sensor distributed fusion feedback specifically includes:

[0095] First, the fusion center receives the state estimates of each local filter. According to the covariance cross fusion criterion, the fusion calculation method is as follows:

[0096]

[0097] in, Represents the optimal state estimate after fusion and its corresponding covariance matrix.

[0098] Secondly, the state and its covariance estimate are fed back to each local filter, and the feedback is calculated as follows:

[0099]

[0100] Among them, κ k,s is the feedback weight coefficient, which is adaptively adjusted with the local state estimation covariance matrix to satisfy:

[0101]

[0102] in,· F represents the Frobenius norm. That is, for any matrix D,

[0103] The beneficial effects of the present invention are as follows: first, by introducing a variational Bayesian method, the present invention considers that the target itself is affected by non-Gaussian and non-stationary heavy-tailed noise in the process noise, and is more adaptable to the uncertainty of the external environment in actual engineering applications; secondly, the influence of non-Gaussian and non-stationary skew measurement noise on the target tracking process is considered, that is, the measurement information value of the sensor is corrected, and the accuracy of target tracking is indirectly improved; finally, a multi-sensor distributed fusion feedback algorithm is introduced, the target state of each local filter and its corresponding covariance matrix are online fused and fed back to each local filter, the target state and its corresponding covariance matrix are corrected in real time, which makes up for the shortcomings of a single platform sensor, makes the system have a certain robustness, and improves the accuracy of target tracking. BRIEF DESCRIPTION OF THE DRAWINGS

[0104] Figure 1 It is an experimental comparison diagram of the present invention.

[0105] Figure 2 It is a flow chart of the method of the present invention. DETAILED DESCRIPTION

[0106] The following will be combined with the accompanying drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0107] In actual use, the present invention further includes step S0: initializing the parameters of the linear system. Compared with the physical scenario of traditional multi-sensor target tracking, the present application comprehensively considers the influence of target state process noise and sensor measurement noise. The linear system represented by the discrete-time linear state space model is as follows:

[0108] x k,s =F k x k-1,s +w k-1

[0109] z k,s =H k x k,s +v k,s

[0110] Among them, x k,s is the target state of sensor s at time k, k = 1, 2, 3..., ts, F k is the state transfer matrix from time k-1 to time k, w k-1 is the k-1 moment process noise, and its noise covariance matrix is Q k-1 , the process noise presents non-Gaussian and non-stationary heavy-tail characteristics; z k,s is the measurement value of sensor s at time k, H k is the measurement transfer matrix at time k; v k,s is the measurement noise of sensor s at time k, and its noise covariance matrix is R k,s , the measurement noise presents non-Gaussian and non-stationary Skew characteristics.

[0111] The variational parameter setting and initialization expectations are as follows:

[0112] (1) Set the initial value of the variational parameter at time k-1 as follows:

[0113] Initial estimated state and its error covariance matrix P k-1|k-1,s ;

[0114] Measurement noise covariance matrix R k-1,s The degrees of freedom parameter And the scale matrix is U k-1|k-1,s ;

[0115] Degree of freedom parameter ω k-1 ,ν k-1 Shape parameters Sum rate parameters and

[0116] Mixed random vector τ k The prior concentration parameter α k|k-1 and species distribution vector ξ k The mean vector of

[0117] Nominal (not true) process noise covariance matrix and the nominal measurement noise covariance matrix R k-1,s ;

[0118] Shape parameter β k-1 The mean and covariance matrix of P β,k|k,s ;

[0119] and the measurement value z of sensor s at time k k,s and other related parameters.

[0120] (2) k-time parameter time update:

[0121] State one-step prediction: And its predicted state error covariance matrix:

[0122] P k|k-1,s The degree of freedom parameter and the scale matrix

[0123] Measurement noise covariance moment R k,s The degree of freedom parameter and the scale matrix U k|k-1,s =U k|k-1,s +E[f(γ k )]B k ;

[0124] Degrees of freedom ω k 、ν k Shape parameters Sum rate parameters

[0125] Shape parameter β k,s The covariance matrix P β,k|k,s =σ k (I m -E j [γ k ]K β,k ) and mean

[0126] Mixed random vector τ k Concentration parameters and species distribution vector ξ k Mixed vector in Where E[·] represents the expectation operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter and is usually 1-exp(-4), τ p represents the tuning factor for the forecast error covariance.

[0127] (3) Initialize the expected value at time k:

[0128] Initial expected state and its covariance matrix:

[0129] Initialize the expected value as follows: shape parameter β k,s Expectations

[0130] Auxiliary variable λ k , γ k Expected E 0 [λ k ]=E 0 [γ k ]=1;

[0131] Degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ];

[0132] Measurement noise covariance R k,s Expectations

[0133] Mixed random vector τ k Expectations and species distribution vector ξ k Expectations

[0134] Superscript E 0 [·] means finding the expectation of the 0th variational iteration.

[0135] like Figure 1 As shown, the present invention includes the following steps: S1: Establishing the one-step prediction probability density P(x k,s ∣z 1:k-1 ) model, measure likelihood probability density P(z k,s ∣x k,s ) model and log joint probability density Model;

[0136] In step S1, P(x k,s ∣z 1:k-1 ) and P(z k,s∣x k,s ) Prior model, logarithmic joint probability density model The model specifically includes the following steps:

[0137] S1.1: Establish the state one-step prediction probability density P(x k,s ∣z 1:k-1 ) Prior model:

[0138] An HTM distribution is proposed to describe the non-Gaussian and non-stationary heavy-tailed characteristics of state process noise. The heavy-tailed and non-stationary noise are modeled by a hybrid method, in which the mixing coefficient is adaptively learned according to the sequential measurement. k,s ∣z 1:k-1 ) can be expressed in Gaussian hierarchical form, introducing the classification distribution vector ξ k , specifically expressed as:

[0139]

[0140] Where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, and G(·; a, b) represents a gamma distribution with shape parameter a, rate parameter b, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter; IW(·; t, T) represents the inverse Wishart distribution, whose degrees of freedom are t and the scaling matrix is T; Cat(·; m, M) represents the species distribution, whose mixing random vector is m and the number of mixing elements is M; Dir(·; θ, M) represents the Dirichlet distribution, whose concentration parameter is θ and the number of mixing elements is M.

[0141] S1.2: Establish the measurement likelihood probability density P(z k,s ∣x k,s ) Prior model:

[0142] Propose GSMGIW distribution and introduce auxiliary variable γ k Written in the form of Gaussian-Gamma mixed distribution, the measurement noise covariance matrix R k,s Establish a priori model, then P(z k,s ∣x k,s ) can be expressed as:

[0143]

[0144] Among them, H k x k,s and R k,s denote the mean vector and covariance matrix respectively, U k|k-1,s Denote the degrees of freedom and the scale matrix respectively. k >0 is the mixing parameter, ν krepresents the degree of freedom, f(γ k ) and IG(γ k ν k / 2,ν k / 2) are expressed as positive scaling function and inverse gamma function respectively. k,s is a shape parameter that controls the degree of skewness. k,s ≠0 is asymmetric, otherwise it is symmetric. In particular, when β k,s = 0, and f(γ k )=γ k , GSMGIW degenerates into the Normal-Gamma-inverse Wishart (NGIW) distribution model.

[0145] S1.3: Constructing a joint probability density model

[0146] According to the state one-step probability density P(x k,s ∣z 1:k-1 ) prior model and step S1.2 to establish the measurement likelihood probability density P(z k,s ∣x k,s ) prior model, rewritten as a layered Gaussian, then the logarithmic joint probability density can be expressed as:

[0147]

[0148] in, represents the unknown parameter set, Z represents the measurement data set, represents the logarithm of the joint probability density with respect to a natural constant.

[0149] In step S3, β is estimated k,s and x k,s The details are as follows:

[0150] Separate orders and Apply the variational iteration formula in step S2, and based on the conjugation of the exponential distribution, approximate the posterior q(β k,s ) and q(x k,s ) are updated to conform to the Gaussian distribution:

[0151]

[0152] Similar to the Kalman filter measurement update step to solve the approximate posterior q i+1 (β k,s ) mean: and its corresponding covariance matrix By calculating the expectation separately: E i+1 (βk,s ), E i+1 (P β,k|k,s ).

[0153] Calculated by Kalman filter measurement update step and

[0154]

[0155] S2: Variational Bayes approximate posterior:

[0156] The model parameters constructed in step S1 are mutually coupled, and the posterior of the unknown parameters cannot be solved analytically. The variational Bayesian method is introduced to approximately solve the posterior of the unknown parameters in S1. The solution formula is as follows:

[0157]

[0158] in, Indicates about Ask for the expected operation. Indicates taking An element of a set, yes Centralized removal The remaining elements, Relatively unknown parameters is an irrelevant constant.

[0159] S3: Use the solution formula of S2 to estimate the shape parameter β k,s and the target state x k,s ;

[0160] S4: Using the solution formula of S2, estimate the state prediction covariance matrix Σ k,s and the measurement noise covariance matrix R k,s ; In step S4, the update estimate Σ k,s and R k,s The details are as follows:

[0161] Separate orders Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, the approximate posterior q(Σ k,s ) and q(R k,s ) are updated to conform to the IW distribution:

[0162]

[0163] Solve the approximate posterior q separately i+1 (Σ k,s ) and q i+1 (R k,s ) degrees of freedom parameter: And the scale matrix parameters: Calculate the expected E from the IW distribution properties i+1 [P k|k-1,s ], E i+1 [R k,s ], calculate the modified state prediction covariance matrix and the measurement noise covariance matrix The state one-step prediction probability density P(x k,s ∣z 1:k-1 ) The prior model shows that the state noise is non-stationary, and here it is designed as a mixture of multiple heavy tails, so P k|k-1,s and Σ k,s For a one-to-many relationship:

[0164]

[0165] Among them, E i [·] means finding the expectation of the i-th variational iteration, i=0:N m -1.

[0166] S5: Use the solution formula of S2 to estimate the mixed random vector τ k and species distribution vector ξ k ;

[0167] In step S5, the estimated value ξ is updated. k and τ k , as shown below:

[0168] S5.1: Update the mixed random vector τ k :

[0169] make Apply the variational iteration formula in step S2, and the approximate posterior q(τ k ) is updated to conform to the Dirichlet distribution:

[0170]

[0171] Solve for the approximate posterior q j+1 (τ k ) concentration parameters: According to the properties of Dirichlet distribution, calculate the natural logarithm expectation E j+1 [τ k ] and E j+1 [logτ k,i ].

[0172] S5.2: Update the species distribution vector ξ k :

[0173] make Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, q(ξ k ) is updated to match the category distribution:

[0174]

[0175] Solve for the approximate posterior q j+1 (ξ k ) of the mixed vector: According to the species distribution properties, calculate the natural logarithm expectation E i+1 [ξ k ].

[0176] S6: Use the solution formula of S2 to estimate the auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ;

[0177] In step S6, update and solve λ k , γ k 、ω k 、ν k , as follows:

[0178] S6.1: Update auxiliary variable λ k and γ k :

[0179] make Applying the variational iteration formula in step S2, the approximate posterior conforms to the gamma distribution, that is:

[0180]

[0181] Solve for the approximate posterior q i+1 (λ k )’s shape and rate parameters: Calculate the expectation E from the properties of the gamma distribution i+1 [λ k ]、E i+1 [logλ k ].

[0182] make Apply the variational iteration formula of step S2, and according to the conjugation of the exponential distribution, approximate the posterior q(γ k ) conforms to the generalized inverse Gaussian distribution (GIG):

[0183]

[0184] Solve for the approximate posterior q i+1 (γ k )’s shape parameters: and Calculate the expectation E by the generalized inverse Gaussian distribution property i+1 [γ k ]、E i+1 [logγ k ].

[0185] S6.2: Update the degree of freedom parameter ω k and ν k :

[0186] Separate orders and Apply the variational iteration formula of step S2, and according to the conjugation of the exponential distribution, approximate the posterior q(ω k ) and q(ν k ) is updated to conform to the gamma distribution:

[0187]

[0188] Solve for the approximate posterior q i+1 (ω k ) and q i+1 (ν k )’s shape parameters: and rate parameters: Calculate the expectation E based on the properties of the gamma distribution i+1 [ω k ] and E i+1 [v k ].

[0189] S7: Repeat steps S3-S6 until the total number of variational iterations is reached or the iteration termination condition is met: ε is usually set to 10 -10 , the iteration ends; the final iteration result is: and Σ k|k,s , and assign it as the initial value of the variational iteration at time k+1.

[0190] S8: Multi-sensor distributed fusion:

[0191] In a distributed network framework, it is often difficult to achieve the expected accuracy by relying on a single local platform for target tracking after the local filter completes the deviation registration. Therefore, it is necessary to use the state values of each local filter obtained by the variational iteration of each platform in step S7. And its covariance matrix P k|k,s The information is fused using the covariance intersection (CI) fusion method to obtain the optimal target state after fusion. and its covariance matrix

[0192] The step S8 of sensor distributed fusion feedback specifically includes:

[0193] First, the fusion center receives the state estimates of each local filter. According to the covariance cross fusion criterion, the fusion calculation method is as follows:

[0194]

[0195] in, Represents the optimal state estimate after fusion and its corresponding covariance matrix.

[0196] Secondly, the state and its covariance estimate are fed back to each local filter, and the feedback is calculated as follows:

[0197]

[0198] Among them, κ k,s is the feedback weight coefficient, which is adaptively adjusted with the local state estimation covariance matrix to satisfy:

[0199]

[0200] Among them, ||·|| F represents the Frobenius norm. That is, for any matrix D,

[0201] The present invention considers the physical scenario in which multiple sensors jointly track the same target. It not only solves the problem of non-Gaussian and non-stationary heavy-tailed process noise in steps S3, S4, S5 and S6; it also solves the problem of non-Gaussian and non-stationary Skew heavy-tailed measurement noise in steps S3, S4 and S6; and proposes a distributed sensor fusion feedback algorithm in step S8, which not only reduces the communication burden of the sensors but also improves the target tracking accuracy.

[0202] according to Figure 1 The RMSE of the state position estimation shows that, excluding the Reference Value, the following are ranked from highest to lowest: VB-ADFKF and HTM-RKF. Theoretically, the Reference Value method offers the best estimation accuracy because it uses real system noise information and performs CI fusion under the same conditions. However, in practical engineering, real system noise is often difficult or impossible to obtain, making it only an optimal estimate under ideal conditions. Although the HTM-RKF considers the effects of non-Gaussian and non-stationary heavy-tailed system noise, its inappropriate modeling of skew-t measurement noise adds additional state estimation error. Therefore, VB-ADFKF is considered the optimal method.

[0203] The present invention aims to solve the existing problems and proposes a new target tracking method based on variational Bayes under non-Gaussian and non-stationary system noise in a distributed fusion framework. Firstly, the variational Bayes method is introduced to consider the influence of the process noise of the target itself, which is non-Gaussian and non-stationary heavy-tailed noise, and is more adaptable to the uncertainty of the external environment in actual engineering applications. Secondly, the influence of non-Gaussian and non-stationary skew measurement noise on the target tracking process is considered, that is, the measurement information value of the sensor is corrected, which indirectly improves the accuracy of target tracking. Finally, a multi-sensor distributed fusion feedback algorithm is introduced to online fuse the target state and the corresponding covariance matrix of each local filter and feed it back to each local filter, so that the target state and the corresponding covariance matrix are corrected in real time, which makes up for the shortcomings of the single platform sensor, makes the system have a certain robustness, and improves the accuracy of target tracking.

[0204] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise, characterized by: The following steps are involved: S1: Establish the one-step prediction probability density P(x k,s ∣z 1:k-1 ) model, measure likelihood probability density P(z k,s ∣x k,s ) model and log joint probability density Model; The step S1 specifically includes the following steps: S1.1: Establish the state one-step prediction probability density P(x k,s ∣z 1:k-1 ) Prior model: Adaptively learn the mixing coefficients based on sequential measurements to transform P(x k,s ∣z 1:k-1 ) can be expressed in Gaussian hierarchical form, introducing the classification distribution vector ξ k , specifically expressed as: Where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, and G(·; a, b) represents a gamma distribution with shape parameter a, rate parameter b, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter; IW(·; t, T) represents the inverse Wishart distribution, with t degrees of freedom and T as the scale matrix; Cat(·; m, M) represents the species distribution, with m as the mixing random vector and M as the number of mixing elements; Dir(·; θ, M) represents the Dirichlet distribution, with θ as the concentration parameter and M as the number of mixing elements; S1.2: Establish the measurement likelihood probability density P(z k,s ∣x k,s ) Prior model: Propose GSMGIW distribution and introduce auxiliary variable γ k Written in the form of Gaussian-Gamma mixed distribution, the measurement noise covariance matrix R k,s Establish a priori model, then P(z k,s ∣x k,s ) can be expressed as: Among them, H k x k,s and R k,s denote the mean vector and covariance matrix respectively, denote the degrees of freedom and scale matrix respectively; γ k >0 is the mixing parameter, ν k represents the degree of freedom, f(γ k ) and IG(γ k ν k / 2,ν k / 2) are expressed as positive scaling function and inverse gamma function respectively; β k,s is a shape parameter that controls the degree of skewness. k,s ≠0 is asymmetric, otherwise it is symmetric; when β k,s = 0, and f(γ k )=γ k , GSMGIW degenerates into the normal gamma inverse Wishart distribution model; S1.3: Constructing a joint probability density model According to the state one-step probability density P(x k,s ∣z 1:k-1 ) prior model and step S1.2 to establish the measurement likelihood probability density P(z k,s ∣x k,s ) prior model, rewritten as a layered Gaussian, then the logarithmic joint probability density can be expressed as: in, represents the unknown parameter set, Z represents the measurement data set, represents the logarithmic joint probability density with a base of natural constants S2: Variational Bayes approximate posterior: The variational Bayesian method is introduced to approximately solve the posterior of the unknown parameters in S1. The solution formula is as follows: in, Indicates about Ask for the expected operation. Indicates taking An element of a set, yes Centralized removal The remaining elements, Relatively unknown parameters is an irrelevant constant; S3: Use the solution formula of S2 to estimate the shape parameter β k,s , target state x k,s , state prediction covariance matrix Σ k,s , measurement noise covariance matrix R k,s , mixed random vector τ k , species distribution vector ξ k , auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ; S4: Repeat step S3 until the iteration termination condition is reached: The iteration ends; the final iteration result is obtained: and Σ k|k,s , and assign it as the initial value of the variational iteration at time k+1; S5: Multi-sensor distributed fusion: The state values of each local filter obtained by iterative variation of each platform in step S4 are And its covariance matrix P k|k,s The information is fused using the covariance cross CI fusion method to obtain the optimal target state after fusion. and its covariance matrix 2. The target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise according to claim 1, characterized in that: The system also includes step S0, initializing the linear system parameters. Compared with the physical scenario of traditional multi-sensor target tracking, the effects of target state process noise and sensor measurement noise are comprehensively considered. The linear system represented by the discrete-time linear state space model is as follows: x k,s =F k x k-1,s +w k-1 z k,s =H k x k,s +v k,s Among them, x k,s is the target state of sensor s at time k, F k is the state transfer matrix from time k-1 to time k, w k-1 is the k-1 moment process noise, and its noise covariance matrix is Q k-1 , the process noise presents non-Gaussian and non-stationary heavy-tail characteristics; z k,s is the measurement value of sensor s at time k, H k is the measurement transfer matrix at time k; v k,s is the measurement noise of sensor s at time k, and its noise covariance matrix is R k,s , the measurement noise presents non-Gaussian and non-stationary Skew characteristics.

3. The target tracking method based on variational Bayesian in the non-Gaussian and non-stationary system noise according to claim 2, characterized in that: The variational parameter settings and initialization expectations in the linear system initialization are as follows: (1) Set the initial value of the variational parameter at time k-1 as follows: Initial estimated state and its error covariance matrix Measurement noise covariance matrix R k-1,s The degrees of freedom parameter And the scale matrix is U k-1|k-1,s ; Degree of freedom parameter ω k-1 ,ν k-1 Shape parameters Sum rate parameters and Mixed random vector τ k The prior concentration parameter α k|k-1 and species distribution vector ξ k The mean vector of Nominal (not true) process noise covariance matrix and the nominal measurement noise covariance matrix R k-1,s ; Shape parameter β k-1 The mean and covariance matrix of and the measurement value z of sensor s at time k k,s and other related parameters; (2) k-time parameter time update: State one-step prediction: And its predicted state error covariance matrix: P k|k-1,s The degrees of freedom parameter and the scale matrix Measurement noise covariance moment R k,s The degrees of freedom parameter and the scale matrix U k|k-1,s =U k|k-1,s +E[f(γ k )]B k ; Degrees of freedom ω k 、ν k Shape parameters Sum rate parameters Shape parameter β k,s The covariance matrix P β,k|k,s =σ k (I m -E j [γ k ]K β,k ) and mean Mixed random vector τ k Concentration parameters and species distribution vector ξ k Mixed vector in Where E[·] represents the expectation operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter and is usually 1-exp(-4), τ p represents the tuning factor of the forecast error covariance; (3) Initialize the expected value at time k: Initial expected state and its covariance matrix: Initialize the expected value as follows: shape parameter β k,s Expectations Auxiliary variable λ k , γ k Expected E 0 [λ k ]=E 0 [γ k ]=1; Degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ]; Measurement noise covariance R k,s Expectations Mixed random vector τ k Expectations and species distribution vector ξ k Expectations Superscript E 0 [·] means finding the expectation of the 0th variational iteration.

4. The target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise according to claim 1, characterized in that: In step S3, the shape parameter β is estimated. k,s and the target state x k,s The specific steps include: Separate orders and Apply the variational iteration formula in step S2, and based on the conjugation of the exponential distribution, approximate the posterior q(β k,s ) and q(x k,s ) are updated to conform to the Gaussian distribution: Similar to the Kalman filter measurement update step: Solve for the approximate posterior q i+1 (β k,s ) mean: and its corresponding covariance matrix E i+1 (β k,s ), E i+1 (P β,k|k,s ) are calculated as 5. The target tracking method based on variational Bayesian in the non-Gaussian and non-stationary system noise according to claim 1, characterized in that: In step S3, the state prediction covariance matrix Σ is updated and estimated. k,s and the measurement noise covariance matrix R k,s The specific steps include: Separate orders Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, the approximate posterior q(Σ k,s ) and q(R k,s ) are updated to conform to the IW distribution: Solve the approximate posterior q separately i+1 (Σ k,s ) and q i+1 (R k,s ) degrees of freedom parameter: And the scale matrix parameters: Calculate the expected E from the IW distribution properties i+1 [P k|k-1,s ], E i+1 [R k,s ], calculate the modified state prediction covariance matrix and the measurement noise covariance matrix The state one-step prediction probability density P(x k,s ∣z 1:k-1 ) The prior model shows that the state noise is non-stationary, and here it is designed as a mixture of multiple heavy tails, so P k|k-1,s and Σ k,s For a one-to-many relationship: Among them, E i [·] means finding the expectation of the i-th variational iteration, i=0:N m -1.

6. The target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise according to claim 1, characterized in that: In step S3, the mixed random vector ξ is updated and estimated. k and the species distribution vector τ k , specifically including the following steps: S5.1: Update the mixed random vector τ k : make Apply the variational iteration formula in step S2, and the approximate posterior q(τ k ) is updated to conform to the Dirichlet distribution: Solve for the approximate posterior q j+1 (τ k ) concentration parameters: According to the properties of Dirichlet distribution, calculate the natural logarithm expectation E j+1 [τ k ] and E j+1 [logτ k,i ]; S5.2: Update the species distribution vector ξ k : make Apply the variational iteration formula in step S2. According to the conjugation property of the exponential distribution, q(ξ k ) is updated to match the category distribution: Solve for the approximate posterior q j+1 (ξ k ) of the mixed vector: According to the species distribution properties, calculate the natural logarithm expectation E i+1 [ξ k ].

7. The target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise according to claim 1, characterized in that: In step S3, the auxiliary variable λ is updated and solved. k , γ k and its corresponding degree of freedom parameter ω k 、ν k , the specific steps are as follows: S6.1: Update auxiliary variable λ k and γ k : make Applying the variational iteration formula in step S3, the approximate posterior conforms to the gamma distribution, that is: Solve for the approximate posterior q i+1 (λ k )’s shape and rate parameters: Calculate the expectation E from the properties of the gamma distribution i +1 [λ k ]、E i+1 [logλ k ]; make Apply the variational iteration formula in step S3, and based on the conjugation of the exponential distribution, approximate the posterior q(γ k ) conforms to the generalized inverse Gaussian distribution GIG: Solve for the approximate posterior q i+1 (γ k )’s shape parameters: and Calculate the expectation E by the generalized inverse Gaussian distribution property i+1 [γ k ]、 E i+1 [logγ k ]; S6.2: Update the degree of freedom parameter ω k and ν k : Separate orders and Apply the variational iteration formula in step S3, and based on the conjugation of the exponential distribution, approximate the posterior q(ω k ) and q(ν k ) is updated to conform to the gamma distribution: Solve for the approximate posterior q i+1 (ω k ) and q i+1 (ν k )’s shape parameters: and rate parameters: Calculate the expectation E from the properties of the gamma distribution i+1 [ω k ] and E i+1 [v k ].

8. The target tracking method based on variational Bayesian in the presence of non-Gaussian and non-stationary system noise according to claim 1, characterized in that: The step S5 of sensor distributed fusion feedback specifically includes the following steps: First, the fusion center receives the state estimates of each local filter, and according to the covariance cross fusion criterion, the fusion calculation method is as follows: in, Represents the optimal state estimate after fusion and its corresponding covariance matrix; Secondly, the state and its covariance estimate are fed back to each local filter, and the feedback is calculated as follows: Among them, κ k,s is the feedback weight coefficient, which is adaptively adjusted with the local state estimation covariance matrix to satisfy: Among them, ||·|| F represents the Frobenius norm; that is, for any matrix D,

Images