A filtering method based on probabilistic data association in clutter environment with unknown observation noise

Abstract

The present invention discloses a filtering method based on probabilistic data association in a clutter environment where the observation noise is unknown. The method belongs to the field of adaptive filtering and addresses the problem that the observation noise in a clutter environment is difficult to estimate. A filtering method based on a probabilistic data association algorithm is proposed. The purpose is to solve the problem that the observation noise in a clutter environment is difficult to estimate and the accurate position of the target is difficult to determine. The core technology includes using a probabilistic data association algorithm to track the target in the clutter environment to obtain the target position, and using the obtained new information for observation noise calculation and filter update. Simulation results show that the invented method has good estimation performance in a clutter environment where the observation noise is much greater than the process noise.

Description

Technical Field

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

[0002] The Kalman filter provides a recursive optimal solution to linear motion filtering problems. Due to its low computational cost and minimal storage requirements, it is widely used in fields such as navigation, state estimation, and target tracking. The excellent performance of the Kalman filter relies primarily on 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, obtaining prior knowledge of the motion system is difficult or even impossible, limiting the practical use of the Kalman filter.

[0003] Research has shown that using inaccurate process and observation noise covariances in Kalman filters can reduce state estimation performance and even cause the Kalman filter to diverge. Process noise is time-varying and difficult to obtain for different targets in different environments. To achieve better filtering performance, we need to obtain observation noise values ​​that are closer to the true value for filtering predictions.

[0004] Since the observation noise estimation method 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 clutter environments, it is difficult to determine the exact position of the radar detection target. The observation noise estimation method in clutter environments is lacking. If more accurate observation noise is obtained, the Kalman filter can achieve better filtering effect. Summary of the Invention

[0005] The present invention proposes a filtering method based on probabilistic 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 poor Kalman filter performance.

[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 object's observation matrix, 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, and the observation matrix H are known and meet the conditions of complete controllability and observability and are constant values. Assume 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] Where 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 transition matrix, x(k-1|k-1) represents the posterior estimate of the target state at time k-1, and z(k|k-1) represents the observation value corresponding to the prior predicted state of the target at time k based on 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 noise;

[0029] Step 7: Update the calculated 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 true, and β0(k) represents the probability that all measured values ​​at time k do not come from the target. -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 confirmed measurements up to time k+1, It means 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 measurement 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 the value of 0-m. When i is 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)".

[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] This invention addresses the problem of large radar target tracking errors in clutter environments with unknown observation noise. It optimizes tracking performance through an adaptive Kalman filter method based on a probabilistic data association algorithm. The invention iteratively calculates the combined innovations obtained from the probabilistic data association algorithm to estimate an observation noise matrix R close to the true value, which is then used in the tracking algorithm. Simulation results demonstrate that this invention achieves good estimation performance in clutter environments where the observation noise is significantly 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 value for the simulation.

[0050] Figure 3 Estimate the RMSE of the observed noise matrix and the true value 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 noise-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, let α=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 is 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. To solve this problem, the present invention uses the screening idea of ​​probabilistic data association (PDA). The probability of each observation coming from the target is calculated based on the previous observations. The probability is used as the weight to weight the observations to obtain the conditional mean, and the conditional mean is used to replace the target observation value z. 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 the Kalman filter as a 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) for 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 innovation 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] Obtaining and filtering observations

[0079] Obtain 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 considered to be measurements that may come 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, it takes the prior state observation value z(k|k-1).

[0082] PDA algorithm calculates event probability

[0083] The probability that each observation value 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 measurement 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 the value of 0-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)".

[0086] Weighted conditional mean

[0087] Using the calculated probability as the weight, the observation value is weighted to obtain the conditional mean as the observation value of the target at time k, and the combined new information is obtained based on 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] The simulation results and the effectiveness of the method are compared. 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). 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 ]), calculate the initial value of R value 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 over 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 The following table compares the filtering performance of the two methods in this simulation environment. The vertical axis of the figure plots the distance between the two algorithms' posterior predictions and the actual target location. As can be seen from the figure, the filtering performance of the two algorithms is relatively stable. After the 10th second, the deviation between the posterior predictions and the actual target location of both algorithms remains stable within a range of 5.5-7 meters. The results after filtering stabilizes (after 50 seconds of observation) are averaged and presented in the table below. As can be seen from the table, when the basic PDA algorithm is used for filtering with a known R value, the mean distance between the posterior estimate and the actual target location is 6.1056 meters. When the proposed algorithm is used with an unknown R value, the mean distance between the posterior estimate and the actual target location is 6.2404 meters, a deviation of 0.14 meters. This is because the proposed algorithm fails to obtain a realistic noise matrix. When the process noise is unknown, the filtering performance is comparable to that of the PDA algorithm when the process noise is known, demonstrating the superior performance of the proposed algorithm 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 object's observation matrix, 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, and the observation matrix H are known and meet the conditions of complete controllability and observability and are constant values. Assume 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) Where 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 transition matrix, x(k-1|k-1) represents the posterior estimate of the target state at time k-1, and z(k|k-1) represents the observation value corresponding to the prior predicted state of the target at time k based on 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 noise; Step 7: Update the calculated 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 true, and β0(k) represents the probability that all measured values ​​at time k do not come from the target. -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 according to claim 1, characterized in that: The specific calculation method of step 4 is: where Z k represents the cumulative set of confirmed measurements up to time k+1, It means 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 measurement 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 the value of 0-m. When i is 0, it means that the event "the m measurements at time k do not come 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 according to 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 according to 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