Target Tracking Method under Unknown Probability Skew and Heavy-Tailed Noise

Abstract

The present invention discloses a target tracking method under unknown probability skew and heavy-tailed noise, comprising the following steps: first, in a distributed target tracking framework, prior modeling is performed based on the characteristics obeyed by the unknown parameters; second, the posterior distribution parameters such as the target state, system deviation, and noise covariance are jointly inferred through fixed-point iteration of the standard variational Bayesian method; and finally, the fusion and correction of the local state estimates are achieved based on the covariance cross-fusion strategy. The proposed method comprehensively considers the multiple effects of process noise that occurs randomly with unknown probability, abnormal heavy-tailed measurement noise, and unknown and time-varying system deviations during the target tracking process. It can effectively estimate unknown parameters such as the target state, noise covariance, and system deviation, 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 a target tracking method in which process noise comprises skew noise or heavy-tail noise occurring randomly with unknown probability, unknown and time-varying system deviation and heavy-tail measurement noise. Background Art

[0002] Today's battlefields are complex and ever-changing, with increasingly complex electromagnetic environments. Using single sensor information to guide missiles inevitably faces challenges such as inaccurate detection, susceptibility to enemy deception and jamming, and target information interruption. Therefore, in sensor networks, information collected from various sensors is collaboratively fused to improve overall system performance. However, this fusion process relies on precise sensor registration. Furthermore, in practical engineering, targets are often affected by the external environment, and severe process noise can cause targets to become outliers. Therefore, the impact of process noise cannot be ignored during target tracking. Furthermore, sensors are also subject to measurement noise when observing targets. Using erroneous / inaccurate process noise and measurement noise can lead to significant estimation errors and even filter divergence. Therefore, it is necessary to estimate noise information while performing sensor offset registration to improve target tracking accuracy.

[0003] Some traditional registration algorithms, such as least squares (LS) and maximum likelihood estimation (MLE), can effectively estimate system biases when noise is negligible or minimal. However, these algorithms are offline and lack real-time performance. In recent years, distributed sensor networks have emerged to reduce the communication burden of sensors and improve their efficiency. Okello et al. proposed an equivalent measurement method for distributed track-level registration, which involves expanding the sensor biases to a state vector. This method designs two parallel extended Kalman filters (EKFs) to estimate the state and sensor biases of the same target after the measurement tracks are associated. Although this algorithm is easy to implement and can be processed online, achieving improved real-time performance and accounting for noise, it suffers from poor convergence when the noise is simply Gaussian white noise and is time-invariant, and the expanded state dimension increases the communication burden. Y. Huang et al. applied the variational Bayesian (VB) method to target tracking. VB is an online processing algorithm that comprehensively considers the effects of process noise and measurement noise. First, the noise information covariance matrix is modeled a priori, and then a variational iteration is used to solve it. This consideration of noise improves target tracking accuracy, but does not account for the influence of system bias. Y. Huang et al. also applied the idea of variational iteration. In addition to considering the influence of measurement noise, they also modeled the system bias with a Gaussian distribution. They then used VB fixed-point iteration to estimate a posteriori parameters such as the system bias and the measurement noise covariance matrix. However, they did not consider the potential influence of external process noise on the target itself. Furthermore, no relevant scholars have yet addressed the issue of state process noise being heavy-tailed noise or skew noise that occurs randomly with unknown probability.

[0004] In summary, existing algorithms are rarely able to handle various problems that arise during target tracking, as follows:

[0005] 1. During target tracking, the target itself may be randomly affected by heavy-tailed process noise or skew process noise with unknown probability, which may cause the target itself to generate outliers and make target tracking difficult;

[0006] 2. Sensor measurements may be subject to the dual influence of abnormal heavy-tailed measurement noise and unknown and time-varying system bias, which may lead to inaccurate sensor measurement information or even loss of target;

[0007] 3. The working environment in actual engineering is complex and changeable, and the electromagnetic environment is becoming increasingly complex. The use of single sensor information will inevitably face problems such as inaccurate detection, susceptibility to enemy deception and interference, and easy interruption of target information. Summary of the Invention

[0008] The purpose of the present invention is to provide a target tracking method under unknown probability Skew and heavy-tailed noise, which can handle the process noise involved in multi-sensor target tracking, such as Skew noise or heavy-tailed noise that appears randomly with unknown probability, unknown and time-varying system deviation, and heavy-tailed measurement noise, thereby improving the target tracking accuracy and robustness.

[0009] The technical solution adopted in the present invention is:

[0010] Target tracking methods under unknown probability Skew and heavy-tailed noise, including

[0011] The following steps are involved:

[0012] S1: A physical scenario where multiple sensors jointly track the same target. During the target tracking process, the effects of target state process noise, sensor measurement noise, and system bias are comprehensively considered. Therefore, the target state model, measurement model, and system bias model are established.

[0013] S2: Set the initial values and initialization expectations of the variational iterative solution parameters at time k-1, k = 1, 2, 3..., ts, ts is the total simulation time;

[0014] S3: 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;

[0015] S4: Variational Bayes Approximate Posterior:

[0016] The model parameters constructed in step S3 are mutually coupled, and it is impossible to solve the posterior of the unknown parameters analytically.

[0017] The Bayesian method approximates the posterior of the unknown parameters in S3, and the solution formula is as follows:

[0018]

[0019] 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, Expressed as a natural constant

[0020] The logarithm with base is the approximate posterior probability density;

[0021] S5: Set the total number of variational iterations to i represents the i-th iteration of the variational method,

[0022] S6: Update the state prediction covariance matrix P k|k-1,s and the measurement noise covariance matrix R k,s ;

[0023] S7: Update the solution shape parameter β k , system deviation η k,s and the target state x k,s ;

[0024] S8: Update the auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ;

[0025] S9: Update the solution for the random variable ξ k and Skew noise occurrence probability π k ;

[0026] S10: Repeat steps S6-S9 until the preset total number of variational iterations is reached or the iteration termination condition is met:

[0027] ε is generally set to 1e-15, and the iteration ends;

[0028] S11: Output variational iteration results:

[0029] Output iteration results: and Σ k|k,s , and assign it as the k+1 moment variable

[0030] Iteration initial value;

[0031] S12: Multi-sensor distributed fusion feedback:

[0032] In a distributed network, it is often difficult to achieve the expected accuracy by relying on a single sensor to track the target after the local filter completes the deviation registration. Therefore, it is necessary to solve the state values of each local filter obtained by the variational iterative solution in step S11. And its covariance matrix P k|k,s The information is fused and the optimal target state after fusion is obtained and its covariance matrix Feedback to each local filter.

[0033] The Bayesian statistical formulas of the discrete target state model, measurement model, and system deviation model in step S1 are specifically as follows:

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

[0035] z k,s =H k x k,s +η k,s +v k,s

[0036] η k,s =η k-1,s +n k-1,s s=1,2,3,...,S

[0037] Among them, x k,s is the target state of sensor s at time k, k = 1, 2, 3..., ts, (x k ,y k ) is the position coordinate, Indicates the corresponding speed, F k-1 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 ;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 ;η k,s is the system deviation of sensor s at time k, n k-1,s is the system deviation noise of sensor s at time k-1, and its noise covariance matrix is Q η,k-1 .

[0038] The step S2 is to set the variational parameters and initialize the expectations, as shown below:

[0039] S2.1: Set the initial value of the variational parameter at time k-1 as follows:

[0040] Initial estimated state and its error covariance matrix P k-1|k-1,s ; Measurement noise covariance matrix R k-1,s The degrees of freedom parameter And the scale matrix is Degree of freedom parameter ω k-1 ,ν k-1Shape parameters Sum rate parameters and Skew noise occurrence probability π k-1 Shape parameters System deviation η k-1,s The mean and its covariance matrix Σ k-1|k-1,s ; Nominal, i.e., non-real process noise covariance matrix and the nominal measurement noise covariance matrix Shape parameter β k-1 The mean and covariance matrix of P g,k-1|k-1 ; The measurement value z of sensor s at time k k,s ;

[0041] S2.2: k-time parameter time update:

[0042] The parameter time is updated as follows: 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 T k|k-1,s =τ p P k|k-1,s ; Measurement noise covariance matrix R k,s The degrees of freedom parameter and the scale matrix U k|k-1,s =ρU k-1|k-1,s ; degrees of freedom ω k 、ν k Shape parameters Sum rate parameters Skew noise occurrence probability π k Shape parameters Mean of systematic deviations and its covariance matrix Σ k|k-1,s =(Σ k-1|k-1,s +Q η,k-1 ) / ρ;

[0043] in,(·) T represents the transposition operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter and takes 1-exp(-4), τ p represents the tuning factor of the forecast error covariance;

[0044] S2.3: Initialize the expected value at time k:

[0045] Initial expected state and its covariance matrix: Initialize the following expected values: auxiliary variable λ k, γ k Expected E 0 [λ k ]=E 0 [γ k ]=1; degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ]; Bernoulli variable ξ k Expectations 0 [ξ k ]=1; the expectation of the system deviation and its covariance: E 0 [η k ]、E 0 [P η,k ]; superscript E 0 [·] means finding the expectation of the 0th variational iteration.

[0046] In step S3, P(x k,s ∣z 1:k-1 ) and P(z k,s ∣x k,s ) Prior model, specifically including the following steps:

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

[0048] In practical engineering applications, various non-Gaussian noises caused by pulse interference, outliers and modeling artifacts usually have Student-t distribution or skew distribution. k-1 Following the Skew distribution, its probability density function P(w k-1 ) is represented as a Gaussian-Gamma mixture distribution as follows:

[0049]

[0050] where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, G(·; a, b) represents a gamma distribution with shape parameter a and rate parameter b, s(·) and δ(·) are positive skew and scale functions, respectively, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter, shape parameter β k Controls the symmetry and skewness of the Skew distribution; β k ≠0 is asymmetric, if s(λ k )=δ(λ k )=λ kis a typical Skew distribution; β k = 0 is symmetric, if δ(λ k )=λ k , the Skew distribution degenerates into the typical Student-t distribution;

[0051] In order to determine the target state x at time k k Obey the Student-t distribution or the Skew distribution, introduce the binary Bernoulli random variable ξ k According to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, the state prediction covariance matrix P k|k-1,s , random variable ξ k , Skew noise occurrence probability π k and the degree of freedom parameter ω k Establish a priori model and combine the process noise probability density function P(w k-1 ), then P(x k,s ∣z 1:k-1 ) is represented as:

[0052]

[0053] Where IW(·; t, T) represents the inverse Wishart distribution, with t degrees of freedom and T as the scaling matrix; Bn(·,π k ) represents the Bernoulli distribution, and the probability of its value being 1 is π k Be(·;a,b) represents the Beta distribution, whose shape parameters are a and b, and is often used as the conjugate prior for the bivariate Bernoulli distribution; ξ k The values are 0 and 1, ξ k =1 indicates state x k Obey the Skew distribution, otherwise ξ k =0 indicates state x k Obey the Student-t distribution;

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

[0055] Considering that the measurement noise is abnormal heavy-tailed noise, the Student-t distribution is used to model it and the auxiliary variable γ is introduced. k Written in the form of Gaussian-Gamma mixed distribution, according to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, the measurement noise covariance matrix R k,s , system deviation η k,s and the degree of freedom parameter ν k Establish a priori model, then P(z k,s ∣x k,s ) is represented as:

[0056]

[0057] S3.3: Constructing a logarithmic 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 S3.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 is represented as:

[0059]

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

[0061] In step S6, the solution P is updated. k|k-1,s and R k,s , as shown below:

[0062] Separate orders Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, the approximate posterior q(P k|k-1,s ) and q(R k,s ) are updated to IW distribution respectively:

[0063]

[0064] Solve the approximate posterior q separately i+1 (P k|k-1,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

[0065]

[0066] Among them, E i [·] means finding the expectation of the i-th variational iteration, and the superscript i means the i-th variational iteration.

[0067] In step S7, the updated solution β k ,η k,s and x k,s , as shown below:

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

[0069]

[0070] Similar to the Kalman filter measurement update step, the approximate posterior q is solved separately i+1 (β k ) and q i+1 (η k,s ) mean: and its corresponding covariance matrix According to the properties of Gaussian distribution, we can calculate the expectation: E i+1 [β k ]、E i +1 [P β,k ]; Calculate expected E i+1 [η k,s ] and the expected covariance matrix of the corrected system deviation

[0071] Calculated by Kalman filter measurement update step and

[0072]

[0073] In step S8, update and solve λ k , γ k 、ω k 、ν k , as shown below:

[0074] S8.1: Update auxiliary variable λ k and γ k :

[0075] make Applying the variational iteration formula in step S4, the approximate posterior is modeled as a generalized inverse Gaussian distribution GIG, that is:

[0076]

[0077] Solve for the approximate posterior q i+1 (λ k ) shape parameters and Calculate the expected E by solving the existing standard GIG distribution formula i+1 [λ k ]、E i+1 [logλ k ];

[0078] make Apply the variational iteration formula in step S4, and based on the conjugation of the exponential distribution, approximate the posterior q(γ k ) is updated to be gamma distributed:

[0079]

[0080] 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 ];

[0081] S8.2: Update the degree of freedom parameter ω k and ν k :

[0082] Separate orders and Apply the variational iteration formula in step S4, and based on the conjugation of the exponential distribution, approximate the posterior q(ω k ) and q(ν k ) is updated to be gamma distributed:

[0083]

[0084]

[0085] 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 ].

[0086] In step S9, the updated solution ξ is obtained. k and π k , as shown below:

[0087] S9.1: Update the random variable ξ k :

[0088] make Apply the variational iteration formula in step S4, then the approximate posterior q(ξ k ) is updated to Bernoulli distribution, and ξ is calculated respectively. k The probability of taking 0 Pr i+1 (ξ k =0) and ξ k Take Pr of 1 i+1 (ξ k =1), calculate the expected E according to the distribution properties of Bn i+1 [ξ k ];

[0089] S9.2: Update Skew noise occurrence probability π k :

[0090] make Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, q(π k ) is updated to a Beta distribution:

[0091]

[0092] Solve for the approximate posterior q i+1 (π k )’s shape parameters: According to the Be distribution properties, calculate the natural logarithm expectation E i+1 [logπ k ]、E i+1 [log(1-π k )].

[0093] The step S12 of sensor distributed fusion feedback specifically includes:

[0094] 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:

[0095]

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

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

[0098]

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

[0100]

[0101] Among them, ||·|| F represents the Frobenius norm; that is, for any matrix D,

[0102] The beneficial effects of the present invention are as follows: through the above technical solution, the present invention proposes a new target tracking method under unknown probability skew and heavy-tailed noise in a distributed fusion framework to address the existing problems:

[0103] First, the variational Bayesian method is introduced to consider the influence of the target itself on the process noise, which is heavy-tailed noise or skew noise that appears randomly with unknown probability. This method is more suitable for the uncertainty of the external environment in actual engineering applications.

[0104] Secondly, the dual effects of abnormal heavy-tailed measurement noise and unknown and time-varying system bias on the target tracking process are considered, that is, the measurement information value of the sensor is corrected, which indirectly improves the target tracking accuracy;

[0105] Finally, a multi-sensor distributed fusion feedback algorithm is introduced to online fuse the target state of each local filter and its corresponding covariance matrix and feed it back to each local filter, thereby correcting the target state and its corresponding covariance matrix in real time, making up for the shortcomings of a single sensor, making the system more robust, and improving the accuracy of target tracking. BRIEF DESCRIPTION OF THE DRAWINGS

[0106] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0107] Figure 1 Flowchart of the present invention. DETAILED DESCRIPTION

[0108] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without creative work are within the scope of protection of the present invention.

[0109] like Figure 1 As shown, the present invention includes the following steps:

[0110] S1: Considering the physical scenario of multiple sensors tracking the same target, the influence of target state process noise, sensor measurement noise and system deviation is comprehensively considered during the target tracking process. Therefore, the target state model, measurement model and system deviation model are established.

[0111] S2: Set the initial values and initialization expectations of the variational iterative solution parameters at time k-1, k = 1, 2, 3..., ts; ts is the total simulation time;

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

[0113] S4: Variational Bayes Approximate Posterior:

[0114] The model parameters constructed in step S3 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 S3. The solution formula is as follows:

[0115]

[0116] 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, represents the logarithm approximated posterior probability density based on a natural constant.

[0117] S5: Set the total number of variational iterations to i represents the i-th iteration of the variational method,

[0118] S6: Update the state prediction covariance matrix P k|k-1,s and the measurement noise covariance matrix R k,s ;

[0119] S7: Update the solution shape parameter β k , system deviation η k,s and the target state xk,s ;

[0120] S8: Update the auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ;

[0121] S9: Update the solution for the random variable ξ k and Skew noise occurrence probability π k ;

[0122] S10: Repeat steps S6-S9 until the preset total number of variational iterations is reached or the iteration termination condition is met: ε is generally set to 1e-15, and the iteration ends;

[0123] S11: Output variational iteration results:

[0124] Output iteration results: and Σ k|k,s , and assign it as the initial value of the variational iteration at time k+1;

[0125] S12: Multi-sensor distributed fusion feedback:

[0126] In a distributed network, it is often difficult to achieve the expected accuracy by relying on a single sensor to track the target after the local filter completes the deviation registration. Therefore, it is necessary to solve the state values of each local filter obtained by the variational iterative solution in step S11. And its covariance matrix P k|k,s The information is fused and the optimal target state after fusion is obtained and its covariance matrix Feedback to each local filter.

[0127] The Bayesian statistical formulas of the discrete target state model, measurement model, and system deviation model in step S1 are specifically as follows:

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

[0129] z k,s =H k x k,s +η k,s +v k,s

[0130] η k,s =η k-1,s +n k-1,s s=1,2,3,...,S

[0131] Among them, x k,s is the target state of sensor s at time k, k = 1, 2, 3..., ts, (x k ,y k ) is the position coordinate, Indicates the corresponding speed, F k-1 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 ;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 ;η k,s is the system deviation of sensor s at time k, n k-1,s is the system deviation noise of sensor s at time k-1, and its noise covariance matrix is Q η,k-1 .

[0132] The step S2 is to set the variational parameters and initialize the expectations, as shown below:

[0133] S2.1: Set the initial value of the variational parameter at time k-1 as follows:

[0134] Initial estimated state and its error covariance matrix P k-1|k-1,s ; Measurement noise covariance matrix R k-1,s The degrees of freedom parameter And the scale matrix is Degree of freedom parameter ω k-1 ,ν k-1 Shape parameters Sum rate parameters and Skew noise occurrence probability π k-1 Shape parameters System deviation η k-1,s The mean and its covariance matrix Σ k-1|k-1,s ; Nominal (not real) process noise covariance matrix and the nominal measurement noise covariance matrix Shape parameter β k-1 The mean and covariance matrix of P g,k-1|k-1 The measurement value z of sensor s at time k k,s and other related parameters;

[0135] S2.2: k-time parameter time update:

[0136] The parameter time is updated as follows: 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 T k|k-1,s =τ p P k|k-1,s ; Measurement noise covariance matrix R k,s The degrees of freedom parameter and the scale matrix U k|k-1,s =ρU k-1|k-1,s ; degrees of freedom ω k 、ν k Shape parameters Sum rate parameters Skew noise occurrence probability π k Shape parameters Mean of systematic deviations and its covariance matrix Σ k|k-1,s =(Σ k-1|k-1,s +Q η,k-1 ) / ρ.

[0137] in,(·) T represents the transpose operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter, usually 1-exp(-4), τ p represents the tuning factor for the forecast error covariance.

[0138] S2.3: Initialize the expected value at time k:

[0139] Initial expected state and its covariance matrix: Initialize the following expected values: auxiliary variable λ k , γ k Expected E 0 [λ k ]=E 0 [γ k ]=1; degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ]; Bernoulli variable ξ k Expectations 0 [ξ k ]=1; the expectation of the system deviation and its covariance: E 0 [η k ]、E 0 [P η,k ]; superscript E 0[·] means finding the expectation of the 0th variational iteration.

[0140] In step S3, P(x k,s ∣z 1:k-1 ) and P(z k,s ∣x k,s ) Prior model, specifically including the following steps:

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

[0142] In practical engineering applications, various non-Gaussian noises caused by pulse interference, outliers and modeling artifacts usually have Student-t distribution or skew distribution. k-1 Following the Skew distribution, its probability density function P(w k-1 ) can be expressed as a Gaussian-Gamma mixture distribution as follows:

[0143]

[0144] where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, G(·; a, b) represents a gamma distribution with shape parameter a and rate parameter b, s(·) and δ(·) are positive skew and scale functions, respectively, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter, shape parameter β k Controls the symmetry and skewness of the Skew distribution. k ≠0 is asymmetric, if s(λ k )=δ(λ k )=λ k is a typical Skew distribution; β k = 0 is symmetric, if δ(λ k )=λ k , the Skew distribution degenerates into the typical Student-t distribution.

[0145] In order to determine the target state x at time k k Obey the Student-t distribution or the Skew distribution, introduce the binary Bernoulli random variable ξ k According to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, select the appropriate distribution to predict the state covariance matrix P k|k-1,s , random variable ξ k , Skew noise occurrence probability π k and the degree of freedom parameter ω kEstablish a priori model and combine the process noise probability density function P(w k-1 ), then P(x k,s ∣z 1:k-1 ) can be expressed as:

[0146]

[0147] Where IW(·; t, T) represents the inverse Wishart distribution, and its degrees of freedom are t , the scale matrix is T; Bn(·,π k ) represents the Bernoulli distribution, and the probability of its value being 1 is π k Be(·;a,b) represents the Beta distribution, whose shape parameters are a and b, and is often used as the conjugate prior for the bivariate Bernoulli distribution; ξ k The values are 0 and 1, ξ k =1 indicates state x k Obey the Skew distribution, otherwise ξ k =0 indicates state x k Obey the Student-t distribution.

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

[0149] Considering that the measurement noise is abnormal heavy-tailed noise, the Student-t distribution is used to model it and the auxiliary variable γ is introduced. k Written in the form of Gaussian-Gamma mixed distribution, according to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, select the appropriate distribution for the measurement noise covariance matrix R k,s , system deviation η k,s and the degree of freedom parameter ν k Establish a priori model, then P(z k,s ∣x k,s ) can be expressed as:

[0150]

[0151] S3.3: Constructing a logarithmic joint probability density model

[0152] According to the state one-step probability density P(x k,s ∣z 1:k-1 ) prior model and step S3.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:

[0153]

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

[0155] In step S6, the solution P is updated. k|k-1,s and R k,s , specifically as follows: In step S6, update and solve P k|k-1,s and R k,s , as shown below:

[0156] Separate orders Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, the approximate posterior q(P k|k-1,s ) and q(R k,s ) are updated to IW distribution respectively:

[0157]

[0158] S6.1: Update the state prediction covariance matrix P k|k-1,s :

[0159] Since the coefficients of the same terms are equal, the approximate posterior q(P k|k-1 )'s degrees of freedom and its scaling matrix can be calculated as:

[0160]

[0161]

[0162] Among them, i represents the i-th variational iteration, and the auxiliary parameters Auxiliary parameters Then the expectation E i+1 [P k|k-1 ]、E i [A k ]、E i [B k ] is calculated as:

[0163]

[0164] Among them, E i [·] represents the expectation of the i-th variational iteration, Definition E i+1 [q k ]=E i [β k ] / E i [λ k], the modified state prediction covariance matrix can be calculated as:

[0165]

[0166] S6.2: Update the measurement noise covariance matrix R k,s :

[0167] Since the coefficients of the same terms are equal, the approximate posterior q(R k,s )'s degrees of freedom and its scaling matrix can be calculated as:

[0168]

[0169]

[0170] Among them, the auxiliary parameter C k =(z k,s -H k x k,s -η k,s )(z k,s -H k x k,s -η k,s ) T , then the expectation E i+1 [R k,s ] and E i [C k ] is calculated as:

[0171]

[0172] Among them, the definition Corrected measurement noise covariance matrix

[0173]

[0174] In step S7, the updated solution β k ,η k,s and x k,s , as shown below:

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

[0176]

[0177] S7.1: Update shape parameter β k :

[0178] Similar to the measurement update step of Kalman filtering, the approximate posterior q i+1 (β k ) and its covariance matrix can be calculated as:

[0179]

[0180] Expected E i+1 [β k ]、E i+1 [P β,k ]and Can be calculated as:

[0181]

[0182] S7.2: Update system deviation η k,s :

[0183] Similar to the measurement update step of Kalman filtering, the approximate posterior q i+1 (η k,s ) and its covariance matrix can be calculated as:

[0184]

[0185] Expected E i+1 [η k,s ]、E i+1 [P η,k ] can be calculated as:

[0186]

[0187] S7.3: Update target state x k,s :

[0188] According to the measurement update step of the Kalman filter, the approximate posterior q i+1 (x k,s ) and its covariance matrix can be calculated as:

[0189]

[0190] In step S8, update and solve λ k , γ k 、ω k 、ν k , as shown below:

[0191] S8.1: Update auxiliary variable λ k and γ k :

[0192] make Applying the variational iteration formula in step S4, the approximate posterior is modeled as a generalized inverse Gaussian distribution (GIG), that is:

[0193]

[0194] Among them, h k ,g k ,p k is the shape parameter of the generalized inverse Gaussian. Then the formula is updated as:

[0195]

[0196] Then E i+1 [λ k ]、E i+1 [logλ k The expectation of ] is approximately calculated by the following formula:

[0197]

[0198] Among them, K p (·) represents a second-kind modified Bessel function, and the subscript p represents an index.

[0199] make Apply the variational iteration formula in step S4, and based on the conjugation of the exponential distribution, approximate the posterior q(γ k ) is updated to be gamma distributed:

[0200]

[0201] Approximate posterior q i+1 (γ k )’s shape and rate parameters: can be calculated as:

[0202]

[0203]

[0204] Expected E i+1 [γ k ]、E i+1 [logγ k ] can be calculated as:

[0205]

[0206] Here, ψ(·) represents the digamma function.

[0207] S8.2: Update the degree of freedom parameter ω k and ν k :

[0208] Separate orders and Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, the approximate posterior q(ω k ) and q(ν k ) is updated to be gamma distributed:

[0209]

[0210]

[0211] Since the coefficients corresponding to the same terms are equal, the approximate posterior q i+1 (ω k ) and q i+1 (ν k )’s shape parameters: and rate parameters: can be calculated as:

[0212]

[0213] Expected E i+1 [ω k ]、E i+1 [ν k ] can be calculated as:

[0214]

[0215] In step S9, the updated solution ξ is obtained. k and π k , as shown below:

[0216] S9.1: Update the random variable ξ k :

[0217] make Apply the variational iteration formula in step S4, and according to the conjugation property of the exponential family distribution, the approximate posterior q(ξ k ) is updated to the Bernoulli distribution, ξ k The probability of taking 0 and 1 can be calculated as:

[0218] Pr i+1 (ξ k =0) = Λ i+1 exp{E i [log(1-π k )]

[0219] +0.5nE i+1 [logλ k ]-0.5E i+1 [λ k ]E i [A k ](E i [P k|k-1,s ]) -1}

[0220] Pr i+1 (ξ k =1) =Λ i+1 exp{E i [logπ k ]

[0221] +0.5nE i+1 [logλ k ]-0.5E i+1 [λ k ]E i [B k ](E i [P k|k-1,s ]) -1}

[0222] Among them, exp(·) represents the exponential function with the natural constant e as the base, and the auxiliary parameter expectation E i [A k ]、E i [B k ] is given in step S6.1, is the normalization constant, the expected E i+1 [ξ k ] can be calculated as:

[0223]

[0224] S9.2: Update Skew noise occurrence probability π k :

[0225] make Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, q(π k ) is updated to a Beta distribution:

[0226]

[0227] Since the coefficients of the same terms are equal, the shape parameter for:

[0228]

[0229]

[0230] Expected E i+1 [logπ k ]、E i+1 [log(1-π k )] can be calculated as:

[0231]

[0232] The step S12 of sensor distributed fusion feedback specifically includes:

[0233] 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:

[0234]

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

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

[0237]

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

[0239]

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

[0241] Through the above technical solution, the present invention first introduces a variational Bayesian method, considers that the target itself is affected by heavy-tailed noise or Skew noise that appears randomly with unknown probability during process noise, and is more adaptable to the uncertainty of the external environment in actual engineering applications; secondly, considers the dual influence of abnormal heavy-tailed measurement noise and unknown and time-varying system deviation on the target tracking process, that is, corrects the measurement information value of the sensor, and indirectly improves the accuracy of target tracking; finally, introduces a multi-sensor distributed fusion feedback algorithm, online fuses the target state of each local filter and its corresponding covariance matrix and feeds it back to each local filter, corrects the target state and its corresponding covariance matrix in real time, makes up for the shortcomings of a single sensor, makes the system have a certain robustness, and improves the accuracy of target tracking.

[0242] The present invention proposes a new target tracking method under unknown probability Skew and heavy-tailed noise under the distributed fusion framework. The forgetting factor ρ is introduced in step S2.2 to play the role of adaptive adjustment parameter. Considering the physical scenario in which multiple sensors jointly track the same target, not only the process noise is identified and solved as heavy-tailed noise or Skew noise that appears randomly with unknown probability in steps S6, S7 and S9; the problem of measurement noise being abnormal heavy-tailed noise and unknown and time-varying system deviation is also solved in steps S6 and S7; and a distributed sensor fusion feedback algorithm is proposed in step S12, which not only reduces the communication burden of the sensor but also improves the target tracking accuracy.

[0243] In the description of the present invention, it should be noted that, for directional words, such as the terms "center", "horizontal", "longitudinal", "length", "width", "thickness", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", "clockwise", "counterclockwise" and the like, indicating directions and positional relationships, are based on the directions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and cannot be understood as limiting the specific scope of protection of the present invention.

[0244] It should be noted that the terms "including" and "having" and any variations thereof in the specification and claims of this application are intended to cover non-exclusive inclusions. For example, a process, method, system, product or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units that are not explicitly listed or are inherent to these processes, methods, products or apparatuses.

[0245] Note that the above are only preferred embodiments of the present invention and the principles of the technology used. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and that various obvious changes, readjustments, and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention is described in detail through the above embodiments, the present invention is not limited to the specific embodiments described herein. Without departing from the concept of the present invention, it may also include many other effective embodiments, and the scope of the present invention is determined by the scope of the appended claims.

Claims

1. A target tracking method under unknown probability skew and heavy-tail noise, characterized by: The steps include: S1: A physical scenario where multiple sensors jointly track the same target. During the target tracking process, the effects of target state process noise, sensor measurement noise, and system bias are comprehensively considered. Therefore, the target state model, measurement model, and system bias model are established. S2: Set the initial values and initialization expectations of the variational iterative solution parameters at time k-1, k = 1, 2, 3..., ts, ts is the total simulation time; S3: 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; S4: Variational Bayes Approximate Posterior: The model parameters constructed in step S3 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 S3. The solution formula is as follows: 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, represents the logarithmic approximation of the posterior probability density with a natural constant as the base; S5: Set the total number of variational iterations to i represents the i-th iteration of the variational method, S6: Update the state prediction covariance matrix P k|k-1,s and the measurement noise covariance matrix R k,s ; S7: Update the solution shape parameter β k , system deviation η k,s and the target state x k,s ; S8: Update the auxiliary variable λ k , γ k and its corresponding degree of freedom parameter ω k 、ν k ; S9: Update the solution for the random variable ξ k and Skew noise occurrence probability π k ; S10: Repeat steps S6-S9 until the preset total number of variational iterations is reached or the iteration termination condition is met: If 1e-15 is taken, the iteration ends; S11: Output variational iteration results: Output iteration results: and Σ k|k,s , and assign it as the initial value of the variational iteration at time k+1; S12: Multi-sensor distributed fusion feedback: In a distributed network, it is often difficult to achieve the expected accuracy by relying on a single sensor to track the target after the local filter completes the deviation registration. Therefore, it is necessary to solve the state values of each local filter obtained by the variational iterative solution in step S11. And its covariance matrix P k|k,s The information is fused and the optimal target state after fusion is obtained and its covariance matrix Feedback to each local filter.

2. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: The Bayesian statistical formulas of the discrete target state model, measurement model, and system deviation model in step S1 are specifically as follows: x k,s =F k-1 x k-1,s +w k-1 z k,s =H k x k,s +η k,s +v k,s or k,s =the k-1,s +n k-1,s s=1,2,3,...,S Among them, x k,s is the target state of sensor s at time k, k = 1, 2, 3..., ts, (x k ,y k ) is the position coordinate, Indicates the corresponding speed, F k-1 is the state transfer matrix from time k-1 to time k, w k-1 k-1 The noise of the moment process is Q k-1 ;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 ;η k,s is the system deviation of sensor s at time k, n k-1,s is the system deviation noise of sensor s at time k-1, and its noise covariance matrix is Q η,k-1 .

3. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: The step S2 is to set the variational parameters and initialize the expectations, as shown below: S2.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 degree of freedom parameter And the scale matrix is Degree of freedom parameter ω k-1 ,ν k-1 Shape parameters Sum rate parameters and Skew noise occurrence probability π k-1 Shape parameters System deviation η k-1,s The mean and its covariance matrix Σ k-1|k-1,s ; Nominal process noise covariance matrix and the nominal measurement noise covariance matrix Shape parameter β k-1 The mean and covariance matrix of P g,k-1|k-1 ; The measurement value z of sensor s at time k k,s ; S2.2: k-time parameter time update: The parameter time is updated as follows: 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 T k|k-1,s =τ p P k|k-1,s ; Measurement noise covariance matrix R k,s The degrees of freedom parameter and the scale matrix U k|k-1,s =ρU k-1|k-1,s ; degrees of freedom ω k 、ν k Shape parameters Sum rate parameters Skew noise occurrence probability π k Shape parameters Mean of systematic deviations and its covariance matrix Σ k|k-1,s =(Σ k-1|k-1,s +Q η,k-1 ) / ρ; in,(·) T represents the transposition operation, ρ represents the variational forgetting factor, which plays the role of adaptive adjustment parameter and takes 1-exp(-4), τ p represents the tuning factor of the forecast error covariance; S2.3: Initialize the expected value at time k: Initial expected state and its covariance matrix: Initialize the following expected values: auxiliary variable λ k , γ k Expected E 0 [λ k ]=E 0 [γ k ]=1; degree of freedom parameter ω k 、ν k Expectations 0 [ω k ]、E 0 [ν k ]; Bernoulli variable ξ k Expectations 0 [ξ k ]=1; the expectation of the system deviation and its covariance: E 0 [η k ]、E 0 [P η,k ]; superscript E 0 [·] means finding the expectation of the 0th variational iteration.

4. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: In step S3, P(x k,s ∣z 1:k-1 ) and P(z k,s ∣x k,s ) Prior model, specifically including the following steps: S3.1: Establish the state one-step prediction probability density P(x k,s ∣z 1:k-1 ) Prior model: In practical engineering applications, various non-Gaussian noises caused by pulse interference, outliers and modeling artifacts usually have Student-t distribution or skew distribution. k-1 Following the Skew distribution, Its probability density function P(w k-1 ) is represented as a Gaussian-Gamma mixture distribution as follows: where N(·; μ, Σ) represents a Gaussian distribution with mean μ and covariance matrix Σ, G(·; a, b) represents a gamma distribution with shape parameter a and rate parameter b, s(·) and δ(·) are positive skew and scale functions, respectively, and λ k >0 is the mixing parameter, ω k is the degree of freedom parameter, shape parameter β k Controls the symmetry and skewness of the Skew distribution; β k ≠0 is asymmetric, if s(λ k )=δ(λ k )=λ k is a typical Skew distribution; β k = 0 is symmetric, if δ(λ k )=λ k , the Skew distribution degenerates into the typical Student-t distribution; In order to determine the target state x at time k k Obey the Student-t distribution or the Skew distribution, introduce the binary Bernoulli random variable ξ k According to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, the state prediction covariance matrix P k|k-1,s , random variable ξ k , Skew noise occurrence probability π k and the degree of freedom parameter ω k Establish a priori model and combine the process noise probability density function P(w k-1 ), then P(x k,s ∣z 1:k-1 ) is represented as: Where IW(·; t, T) represents the inverse Wishart distribution, with t degrees of freedom and T as the scaling matrix; Bn(·,π k ) represents the Bernoulli distribution, and the probability of its value being 1 is π k Be(·;a,b) represents the Beta distribution, whose shape parameters are a and b, and is often used as the conjugate prior for the bivariate Bernoulli distribution; ξ k The values are 0 and 1, ξ k =1 indicates state x k Obey the Skew distribution, otherwise ξ k =0 indicates state x k Obey the Student-t distribution; S3.2: Establish the measurement likelihood probability density P(z k,s ∣x k,s ) Prior model: Considering that the measurement noise is abnormal heavy-tailed noise, the Student-t distribution is used to model it and the auxiliary variable γ is introduced. k Written in the form of Gaussian-Gamma mixed distribution, according to the characteristics of the unknown parameters and the Bayesian statistical theory formula in step S1, the measurement noise covariance matrix R k,s , system deviation η k,s and the degree of freedom parameter ν k Establish a priori model, then P(z k,s ∣x k,s ) is represented as: S3.3: Constructing a logarithmic joint probability density model According to the state one-step probability density P(x k,s ∣z 1:k-1 ) prior model and step S3.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 is represented as: in, represents the unknown parameter set, Z represents the measurement data set of sensor s, represents the logarithm of the joint probability density with respect to a natural constant.

5. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: In step S6, the solution P is updated. k|k-1,s and R k,s , as shown below: Separate orders Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, the approximate posterior q(P k|k-1,s ) and q(R k,s ) are updated to IW distribution respectively: Solve the approximate posterior q separately i+1 (P k|k-1,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 Among them, E i [·] means finding the expectation of the i-th variational iteration, the superscript i means the i-th variational iteration, i=0:

6. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: In step S7, the updated solution β k ,η k,s and x k,s , as shown below: Separate orders and Apply the variational iteration formula in step S4, based on the conjugation property of the exponential distribution, Approximate posterior q(β k ),q(η k,s ) and q(x k,s ) are updated to Gaussian distribution respectively: Similar to the Kalman filter measurement update step, the approximate posterior q is solved separately i+1 (β k ) and q i+1 (η k,s ) mean: and its corresponding covariance matrix According to the properties of Gaussian distribution, we can calculate the expectation: E i+1 [β k ]、E i+1 [P β,k ]; Calculate expected E i+1 [η k,s ] and the expected covariance matrix of the corrected system deviation Calculated by Kalman filter measurement update step and 7. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: In step S8, the updated solution λ k , γ k 、ω k 、ν k , as shown below: S8.1: Update auxiliary variable λ k and γ k : make Applying the variational iteration formula in step S4, the approximate posterior is modeled as a generalized inverse Gaussian distribution GIG, that is: Solve for the approximate posterior q i+1 (λ k ) shape parameters and Calculate the expected E by solving the existing standard GIG distribution formula i+1 [λ k ]、E i+1 [logλ k ]; make Apply the variational iteration formula in step S4, and based on the conjugation of the exponential distribution, approximate the posterior q(γ k ) is updated to be gamma distributed: 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 ]; S8.2: Update the degree of freedom parameter ω k and ν k : Separate orders and Apply the variational iteration formula in step S4, and based on the conjugation of the exponential distribution, approximate the posterior q(ω k ) and q(ν k ) is updated to be gamma distributed: 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 under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: In step S9, the updated solution ξ is obtained. k and π k , as shown below: S9.1: Update the random variable ξ k : make Apply the variational iteration formula in step S4, then the approximate posterior q(ξ k ) is updated to Bernoulli distribution, and ξ is calculated respectively. k The probability of taking 0 Pr i+1 (ξ k =0) and ξ k Take Pr of 1 i+1 (ξ k =1), calculate the expected E according to the distribution properties of Bn i+1 [ξ k ]; S9.2: Update Skew noise occurrence probability π k : make Apply the variational iteration formula in step S4. According to the conjugation property of the exponential distribution, q(π k ) is updated to a Beta distribution: Solve for the approximate posterior q i+1 (π k )’s shape parameters: According to the Be distribution properties, calculate the natural logarithm expectation E i+1 [logπ k ]、E i+1 [log(1-π k )].

9. The target tracking method under unknown probability skew and heavy-tail noise according to claim 1, characterized in that: The step S12 of sensor distributed fusion feedback specifically includes: 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: 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