Filtering method based on probability data association in clutter environment with unknown observation noise

Abstract

The invention discloses a filtering method based on probability data association in a clutter environment with unknown observation noise, belongs to the field of adaptive filtering, and provides a filtering method based on a probability data association algorithm for solving the problem that the observation noise in the clutter environment is difficult to estimate. The method aims to solve the problems that observation noise is difficult to estimate and the accurate position of a target is difficult to determine in a clutter environment. The core technology of the method comprises the following steps: tracking a target in a clutter environment by using a probability data association algorithm to obtain a target position, and using obtained information for observing noise calculation and filtering updating. Simulation results show that the method has good estimation performance in a clutter environment in which observation noise is far greater than process noise.

Description

Technical Field

[0001] The invention belongs to the field of adaptive filtering, and proposes a filtering method based on a probability data association algorithm to solve the problem that observation noise is difficult to estimate in a clutter environment. Background Art

[0002] The Kalman filter provides a recursive optimal solution for the linear motion filtering problem. Due to its low computational cost and small storage space requirements, it is widely used in navigation, state estimation, target tracking and other fields. The excellent performance of the Kalman filter mainly depends on the accurate knowledge of the observed linear motion system, such as the initial state, state transition matrix, measurement matrix, process noise covariance matrix and observation noise covariance matrix. However, in many practical applications, it is difficult or even impossible to obtain prior knowledge of the motion system, which limits the actual use of the Kalman filter.

[0003] Research shows that the use of inaccurate process noise covariance and observation noise covariance by the Kalman filter will reduce the effectiveness of state estimation and even cause the Kalman filter to diverge. In different environments and for different targets, process noise is time-varying and difficult to obtain. In order to obtain better filtering performance, we hope to obtain observation noise values ​​that are closer to the true value for filtering prediction.

[0004] Since the estimation method of observation noise was proposed in the 1870s, the algorithm for calculating radar observation noise based on the continuous target position observed by radar in a stable environment has gradually matured. However, due to the interference of false alarm information in a clutter environment, it is difficult to determine the exact position of the radar detection target, and the observation noise estimation method in a clutter environment is lacking. When obtaining more accurate observation noise, the Kalman filter can achieve better filtering effect. Summary of the invention

[0005] The present invention proposes a filtering method based on probability data association in a clutter environment with unknown observation noise, aiming to solve the problem that the observation noise in a discrete-time linear time-invariant system in a clutter environment is difficult to estimate during radar target detection, resulting in a poor Kalman filter effect.

[0006] In order to achieve the above object, the technical solution of the present invention is:

[0007] A filtering method based on probabilistic data association in a clutter environment with unknown observation noise is proposed. The motion equation of the discrete linear time-invariant system targeted by this method is:

[0008] x k =F*x k-1 +q k-1

[0009] zk =H*x k +r k

[0010] Among them, x k ∈R n is the state of the object at time k, F∈R n is the object state transfer matrix, q k-1 is the process noise of the object at time k-1; z k ∈R n is the observable quantity of the object at time k, H∈R nxn is the observation matrix of the object, r k is the observation noise of the object at time k, where R n represents an n-dimensional real column vector, R nxn represents an n×n real matrix;

[0011] In a clutter environment, the state transfer matrix F, the process noise covariance Q are known, the observation matrix H is known, and it satisfies the conditions of complete controllability and observability and is a constant value. It is assumed that the observation matrix C is reversible and a quantity matrix; the observation quantity z of the system is k is measurable and bounded; assume that the process noise and observation noise are uncorrelated, their values ​​are completely unknown and are positive definite matrices, and that the noise values ​​are bounded; for the above discrete-time linear time-invariant system model, the steps at time k are as follows:

[0012] Step 1: Use the motion equation to estimate the state of the target at time k, and predict the corresponding observation value as the center of the relevant wave gate:

[0013] x(k|k-1)=F*x(k-1|k-1) Formula (1)

[0014] z(k|k-1)=H*x(k|k-1) Formula (2)

[0015] Among them, x(k|k-1) represents the prior estimate of the target state at time k based on the information at time k-1, F represents the state transfer matrix, x(k-1|k-1) represents the posterior estimate of the target state at time k-1; z(k|k-1) represents the observation value corresponding to the prior predicted state of the target at time k according to the observation equation;

[0016] Step 2: Update the covariance:

[0017] P(k|k-1)=F*P(k-1|k-1)F T +Q formula (3)

[0018] S(k)=H*P(k|k-1)H T +R(k-1) Formula (4)

[0019] P(k|k-1) is the covariance between the prior target state x(k|k-1) predicted by equation (1) at time k and the actual target state x(k), P(k-1|k-1) is the covariance between the posterior target state x(k-1|k-1) at time k-1 and the actual target state x(k-1), Q is the process error covariance of the target motion; S(k) is the innovation covariance matrix at time k, which is used to filter the observations in step 3, R(k-1) is the R value calculated at time k-1, [.] T represents the transpose of the matrix [.];

[0020] Step 3: Obtain multiple actual observations at time k, establish a correlation gate with the prior state observation value obtained by formula (2) as the center, and filter out the set of observations that fall into the correlation gate from the actual observations:

[0021]

[0022] Z(k) is the set of observations that fall into the relevant gate at time k, m means there are m observations that fall into the relevant gate, z i (k) represents the i-th observation in the set;

[0023] Step 4: Calculate the probability β that each observation value comes from the target according to the PDA algorithm i ;

[0024] Step 5: Using the calculated probability as the weight, perform weighted calculation on the observations to obtain the conditional mean as the observation value of the target at time k, and use this to calculate the new information V(k+1) at time k+1 of the combination;

[0025] Step 6: Let the state value corresponding to z(k) be x(k), according to the observation equation z(k) = H*x(k|k-1)+r k get:

[0026] α=x(k)-x(k|k-1) Formula (9)

[0027] V(k)=z(k)-z(k|k-1)=H*x(k)+r k -H*x(k|k-1)=H*α+r Formula (10)

[0028] Among them, α is the prior prediction state error, that is, P(k|k-1)=αα T , r is the noise;

[0029] Step 7: Update the calculation to obtain the R(k) value;

[0030] Step 8: Update the Kalman gain K(k), make a posterior estimate of the object state x(k|k), and update the posterior state estimate covariance P(k|k):

[0031] K(k)=P(k|k-1)H T S -1 (k) Formula (15)

[0032] x(k|k)=x(k|k-1)+K(k)V(k) Formula (16)

[0033]

[0034] Among them, x i (k|k) represents θ i (k) is the corresponding posterior estimated state value when β0(k) is true, β0(k) represents the probability that all measured values ​​at time k do not come from the target, S -1 (k) represents the inverse matrix of S(k) of the innovation covariance at time k, K(k) represents the Kalman gain at time k, and x i (k|k) represents the posterior estimate of the target state under the condition that the i-th measurement value at time k comes from the target, and x(k|k) represents the posterior estimate of the target state at time k.

[0035] Furthermore, the specific calculation method of step 4 is:

[0036]

[0037] Where Z k represents the cumulative set of confirmation measurements up to time k+1, It means that it is identical to, β i is identical to the event θ at time k i (k) probability of occurrence, θ i (k) represents the event "the set of measurements confirmed at this time is Z k Under the premise of measuring z i (k) is the measurement from the target", where i takes values ​​between 0 and m. When i takes 0, it indicates the event "none of the m measurements at time k comes from the target, that is, the target measurement is the prior prediction value z(k|k-1)".

[0038] Furthermore, the specific calculation method of step 5 is:

[0039]

[0040] Furthermore, the specific method of step 7 is:

[0041] V(k)V(k) T =(Hα+r)(Hα+r)T =Hαα T H T +Hαr T +rα T H T +rr T Formula (11)

[0042] E(VV T )=HE(αα T )H T +HE(αr T )+E(rα T )H T +E(rr T ) Formula (12)

[0043]

[0044] Among them, E([.]) represents the expected value of the matrix [.], E k ([.]) represents the expected value of the matrix [.] at time k obtained based on the iterative method;

[0045] Calculate the estimated R value R(k) at time k:

[0046] R(k)=E k (VV T )-H*P(k|k-1)H T Formula (14).

[0047] The present invention aims to solve the problem of large error in radar tracking targets in a clutter environment with unknown observation noise, and optimizes the tracking effect through an adaptive Kalman filter method based on a probabilistic data association algorithm. The present invention iteratively calculates the combined new information obtained in the probabilistic data association algorithm, estimates the observation noise matrix R close to the true value, and uses it in the tracking algorithm. Simulation results show that the invention has good estimation performance in a clutter environment where the observation noise is much greater than the process noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 A general flow chart of the operation of the present invention.

[0049] Figure 2 Estimate the first moments of the observation noise matrix with respect to the true values ​​for the simulation.

[0050] Figure 3 Estimate the RMSE of the observed noise matrix and the true values ​​for the simulation.

[0051] Figure 4 This is a comparison chart of the filtering effects of the probabilistic data association algorithm when R is known and the algorithm when R is unknown under the same environment. DETAILED DESCRIPTION

[0052] For the discrete-time linear time-invariant system model, when the process noise covariance matrix is ​​known and the observation noise covariance matrix is ​​unknown, the motion model is:

[0053] x k =Fx k-1 +q k-1

[0054] z k =Hx k +r k

[0055] where x k ∈R n is the state of the object at time k, F∈R n is the object state transfer matrix, q k-1 is the process noise of the object at time k-1; z k ∈R n is the observable quantity of the object at time k, H∈R nxn is the observation equation of the object, r k is the observation noise of the object at time k, where R n represents an n-dimensional real column vector, R nxn Represents an nxn real matrix.

[0056] In a clutter-free environment, the observed value z at each time k is uniquely determined, and its corresponding state x k It is also the only certainty that α=x k -x k|k-1 , then the new information is

[0057] e k =z k -Hx k|k-1 =H(x k -x k|k-1 )+r=Hα+r

[0058] E(e k e k T )=HE(αα T )H T +HE(αr T )+E(rα T )H T +E(rr T )=HPH T +R

[0059] According to the law of large numbers

[0060]

[0061] Right now

[0062]

[0063] However, in a clutter environment, due to the interference of false alarm information, the observation value z at time k k It is not necessarily unique. When there are multiple possible observations or the target is not observed in the prediction area, it is difficult to select the correct observation for calculation. In order to solve this problem, the present invention uses the screening idea of ​​probabilistic data association (PDA), calculates the probability that the observation comes from the target based on the previous observations, and uses the probability as the weight to weight the observations to obtain the conditional mean, and replaces the target observation z with the conditional mean. k Perform calculations.

[0064] Specific motion model establishment

[0065] For discrete linear time-invariant systems with process noise and unknown observation noise, a motion model is built using Kalman filtering as the framework:

[0066] x k =Fx k-1 +q k-1

[0067] z k =Hx k +r k

[0068] Status Update

[0069] At time k, according to the state transfer equation, the target state x(k-1|k-1) at time k-1 is used to make a priori prediction x(k|k-1) of the target state x(k) at time k, and the observation value z(k|k-1) corresponding to x(k|k-1) is obtained according to the observation equation.

[0070] x(k|k-1)=Fx(k-1|k-1)#(1)

[0071] z(k|k-1)=Hx(k|k-1)#(2)

[0072] Among them, x(k|k-1) represents the prior estimate of the target state at time k based on the information at time k-1, F represents the state transfer matrix, x(k-1|k-1) represents the posterior estimate of the target state at time k-1; z(k|k-1) represents the observation value corresponding to the prior predicted state of the target at time k according to the observation equation, and H represents the observation matrix.

[0073] Prior covariance update

[0074] Update the state covariance and new information covariance at time k:

[0075] P(k|k-1)=FP(k-1|k-1)F T +Q#(3)

[0076] S(k)=HP(k|k-1)H T +R(k-1)#(4)

[0077] P(k|k-1) is the covariance between the prior target state x(k|k-1) predicted by (1) at time k and the actual target state x(k), P(k-1|k-1) is the covariance between the posterior target state x(k-1|k-1) at time k-1 and the actual target state x(k-1), Q is the process error covariance of the target motion; S(k) is the innovation covariance matrix at time k, which is used to filter the observations in step 3, R(k-1) is the R value calculated at time k-1, [■] T Represents the transpose of the matrix [■].

[0078] Get and filter observations

[0079] Get multiple actual observations at time k, and establish a correlation gate with the prior state observation value obtained in (2) as the center: Calculate the Mahalanobis distance between each observation value and the prior state observation value z(k|k-1):

[0080] D 2 =v T S -1 v#(5)

[0081] v is the coordinate difference between the measured value and the prior state observation value z(k|k-1), S -1 is the inverse matrix of the innovation covariance matrix. Compare the Mahalanobis distance with the value γ selected from the chi-square distribution table, D 2 The observations <γ are regarded as possible measurements from the target point, and the set of observations that may come from the target at time k is obtained. m is the number of measured values, z i (k) represents the i-th possible measurement from the target at time k. When i is 0, its value is the prior state observation z(k|k-1).

[0082] PDA algorithm calculates event probability

[0083] The probability that each observation comes from the target is calculated according to the PDA algorithm:

[0084]

[0085] Where Z k represents the cumulative set of confirmation measurements up to time k+1, βi is identical to the event θ at time k i (k) probability of occurrence, θ i (k) represents the event "the set of measurements confirmed at this time is Z k Under the premise of measuring z i (k) is the measurement from the target" where i takes values ​​between 0 and m. When i takes 0, it indicates the event "none of the m measurements at time k comes from the target, that is, the target measurement is the prior prediction value z(k|k-1)".

[0086] Weighted conditional mean

[0087] Using the calculated probability as the weight, the observations are weighted to obtain the conditional mean as the observation value of the target at time k, and the combined new information is obtained from this:

[0088]

[0089] Substitute the calculated conditional mean as the only observation value in the noise-free environment into z(k), and its corresponding state value is x(k). According to the observation equation z k =Hx k +r k get

[0090] α=x(k)-x(k|k-1)#(9)

[0091] V(k)=z(k)-z(k|k-1)=H*x(k)+rH*x(k|k-1)=H*α+r#(10)

[0092] α is the prior prediction state error, that is, P(k|k-1)=αα T

[0093] Update the calculated R(k) value

[0094] Find the variance of V(k):

[0095] V(k)V(k) T =(Hα+r)(Hα+r) T =Hαα T H T +Hαr T +rα T H T +rr T #(11)

[0096] E(VV T )=HE(αα T )H T +HE(αr T )+E(rα T )HT +E(rr T )=HPH T +R#(12)

[0097] According to the law of large numbers, we get the recursive formula:

[0098]

[0099] Calculate the estimated R value at time k:

[0100] R(k)=E k (VV T )-HP(k|k-1)H T #(14)

[0101] Using formula (14), the covariance matrix estimation sequence {R(k)} of the real-time measurement noise is obtained. The noise covariance matrix gradually approaches the theoretical value over time, providing noise information for applying the Kalman filter algorithm to generate system state estimation.

[0102] Update the posterior state

[0103] Update the Kalman gain K:

[0104] K(k)=P(k|k-1)H T S -1 (k)#(15)

[0105] Make a posterior estimate of the object state x(k|k):

[0106] x(k|k)=x(k|k-1)+K(k)V(k)#(16)

[0107] Update the posterior state estimate covariance P(k|k):

[0108]

[0109] where x i (k|k) is θ i (k) corresponds to the posterior estimated state value when it is true.

[0110] Simulation result image comparison and method effectiveness demonstration The initial state of the target is [0m, 30m / s, 0m, -100m / s] (initial x coordinate, initial x direction speed, initial y coordinate, initial y direction speed, respectively), and the noise matrix Q value of the simulation process is diag ([0.01m 2 ,0.01m 2 / s 2 ,0.01m 2 ,0.01m 2 / s 2]), the observation noise R value is diag([100m 2 ,100m 2 ]), the initial value of R is calculated as diag([500m 2 ,500m 2 ]), target detection probability Pd = 0.99, clutter is randomly distributed, and the clutter density is 30 / km 2 ,Perform Monte Carlo simulation, and calculate the first-order moment and RMSE between the obtained R value and the true value.

[0111] Figure 2 is the change of the first-order moment of the elements at each position with time. It can be seen that the deviations of the diagonal elements converge to -8.95% and -9.76% respectively, which are about 9.5% different from the actual R value.

[0112] Figure 3 They are the changes of elements at each position of the RMSE matrix. The diagonal elements converge to 12.34 and 11.6, which are close to 12% deviation from the actual R value.

[0113] Figure 4 This is a comparison of the filtering effects of the two methods in this simulation environment. The vertical axis of the figure is the distance between the two algorithms' posterior prediction results and the actual position of the target. As can be seen from the figure, the algorithm filtering effect is relatively stable. Since the 10th second, the distance deviation between the two algorithms' posterior prediction results and the actual position of the target is stable in the range of 5.5m-7m. The average value of the filtering effect after stabilization (after 50 seconds of observation time) is taken, and the results are shown in the table below. As can be seen from the table, when the basic PDA algorithm is used for filtering when the R value is known, the average distance between the posterior estimate and the actual target position is 6.1056 meters; when this algorithm is used when the R value is unknown, the average distance between the posterior estimate and the actual target position is 6.2404 meters, which is 0.14 meters more than the former. This is because this algorithm fails to obtain the real noise matrix. When the process noise is unknown, the filtering effect is not much different from the PDA when the process noise is known, indicating that this algorithm has excellent effect in this simulation environment.

[0114]

Claims

1. A filtering method based on probabilistic data association in a clutter environment with unknown observation noise. The motion equation of the discrete linear time-invariant system targeted by this method is: x k =F*x k-1 +q k-1 z k =H*x k +r k in, x k ∈R n is the state of the object at time k, F∈R n is the object state transfer matrix, q k-1 is the process noise of the object at time k-1; z k ∈R n is the observable quantity of the object at time k, H∈R nxn is the observation matrix of the object, r k is the observation noise of the object at time k, where R n represents an n-dimensional real column vector, R nxn represents an n×n real matrix; In a clutter environment, the state transfer matrix F, the process noise covariance Q are known, the observation matrix H is known, and it satisfies the conditions of complete controllability and observability and is a constant value. It is assumed that the observation matrix C is reversible and a quantity matrix; the observation quantity z of the system is k is measurable and bounded; assume that the process noise and observation noise are uncorrelated, their values ​​are completely unknown and are positive definite matrices, and that the noise values ​​are bounded; for the above discrete-time linear time-invariant system model, the steps at time k are as follows: Step 1: Use the motion equation to estimate the state of the target at time k, and predict the corresponding observation value as the center of the relevant wave gate: x(k|k-1)=F*x(k-1|k-1) Formula (1) z(k|k-1)=H*x(k|k-1) Formula (2) Among them, x(k|k-1) represents the prior estimate of the target state at time k based on the information at time k-1, F represents the state transfer matrix, x(k-1|k-1) represents the posterior estimate of the target state at time k-1; z(k|k-1) represents the observation value corresponding to the prior predicted state of the target at time k according to the observation equation; Step 2: Update the covariance: P(k|k-1)=F*P(k-1|k-1)F T +Q formula (3) S(k)=H*P(k|k-1)H T +R(k-1) Formula (4) P(k|k-1) is the covariance between the prior target state x(k|k-1) predicted by equation (1) at time k and the actual target state x(k), P(k-1|k-1) is the covariance between the posterior target state x(k-1|k-1) at time k-1 and the actual target state x(k-1), Q is the process error covariance of the target motion; S(k) is the innovation covariance matrix at time k, which is used to filter the observations in step 3, R(k-1) is the R value calculated at time k-1, [.] T represents the transpose of the matrix [.]; Step 3: Obtain multiple actual observations at time k, establish a correlation gate with the prior state observation value obtained by formula (2) as the center, and filter out the set of observations that fall into the correlation gate from the actual observations: Z(k) is the set of observations that fall into the relevant gate at time k, m means there are m observations that fall into the relevant gate, z i (k) represents the i-th observation in the set; Step 4: Calculate the probability β that each observation value comes from the target according to the PDA algorithm i ; Step 5: Using the calculated probability as the weight, perform weighted calculation on the observations to obtain the conditional mean as the observation value of the target at time k, and use this to calculate the new information V(k+1) at time k+1 of the combination; Step 6: Let the state value corresponding to z(k) be x(k), according to the observation equation z(k) = H*x(k|k-1)+r k get: α=x(k)-x(k|k-1) Formula (9) V(k) = z(k) - z(k|k - 1) = H * x(k) + r k -H * x(k|k - 1) = H * α + r Equation (10) Among them, α is the prior prediction state error, that is, P(k|k-1)=αα T , r is the noise; Step 7: Update the calculation to obtain the R(k) value; Step 8: Update the Kalman gain K(k), make a posterior estimate of the object state x(k|k), and update the posterior state estimate covariance P(k|k): K(k)=P(k|k-1)H T S -1 (k) Formula (15) x(k|k)=x(k|k-1)+K(k)V(k) Formula (16) Among them, x i (k|k) represents θ i (k) is the corresponding posterior estimated state value when β0(k) is true, β0(k) represents the probability that all measured values ​​at time k do not come from the target, S -1 (k) represents the inverse matrix of S(k) of the innovation covariance at time k, K(k) represents the Kalman gain at time k, and x i (k|k) represents the posterior estimate of the target state under the condition that the i-th measurement value at time k comes from the target, and x(k|k) represents the posterior estimate of the target state at time k.

2. The filtering method based on probabilistic data association in a clutter environment with unknown observation noise as claimed in claim 1, characterized in that: The specific calculation method of step 4 is: Where Z k represents the cumulative set of confirmation measurements up to time k+1, It means that it is identical to, β i is identical to the event θ at time k i (k) probability of occurrence, θ i (k) represents the event "The set of measured values ​​confirmed at this k moment is Z k Under the premise of measuring z i (k) is the measurement from the target", where i takes values ​​between 0 and m. When i takes 0, it means that the event "none of the m measurements at time k comes from the target, that is, the target measurement is the prior prediction value z(k|k-1)".

3. The filtering method based on probabilistic data association in a clutter environment with unknown observation noise as claimed in claim 1, characterized in that: The specific calculation method of step 5 is: Where z0(k)=z(k|k-1) Formula (7) 4. The filtering method based on probabilistic data association in a clutter environment with unknown observation noise as claimed in claim 1, characterized in that: The specific method of step 7 is: V(k)V(k) T =(Hα+r)(Hα+r) T =Haa T H T +Har T +rα T H T +rr T formula(11) E(VV T ) = HE(αα T )H T + HE(αr T ) + E(rα T )H T + E(rr T ) Equation (12) Among them, E([.]) represents the expected value of the matrix [.], E k ([.]) represents the expected value of the matrix [.] at time k obtained based on the iterative method; Calculate the estimated R value R(k) at time k: R(k) = E k (VV T ) - H * P(k|k - 1)H T Equation (14).

Images