Target tracking method of phd filter with adaptive measurement noise and detection probability

Abstract

The application discloses a target tracking method of PHD filtering with adaptive measurement noise and detection probability, and mainly solves the problem that the traditional multi-target tracking filtering algorithm cannot accurately track targets, resulting in the decline of the performance of the filter under the condition that the measurement noise covariance and the detection probability are both unknown in the multi-target tracking scene. The method models the measurement noise covariance and the detection probability with inverse gamma distribution and beta distribution respectively on the basis of the traditional Gaussian mixture probability hypothesis density filter (GM-PHD), updates the weight values of the PHD of the detected and missed targets in the posterior update step by using the properties of the beta function, and then decouples the state and the measurement noise covariance by using the variational Bayesian method, so that the variational iteration approximates the real posterior distribution, and the joint estimation of the target state, the measurement noise covariance and the detection probability is realized. Compared with the traditional tracking algorithm, the application has better tracking performance under the condition that the measurement noise covariance and the detection probability are unknown.

Description

Technical Field

[0001] This invention relates to the fields of random finite set multi-target tracking in radar data processing and variational Bayesian filtering, specifically to a target tracking method using PHD filtering that adapts to measurement noise and detection probability. Background Technology

[0002] Multi-target tracking has always been an important research topic in the field of information fusion. With the development of various sensors such as radar and sonar, this technology has been widely applied in many fields such as aerospace, transportation, and guidance. Existing multi-target tracking methods are mainly divided into two categories. One category is based on traditional data association algorithms to track multiple targets. The principle of this method is to first match multiple measurement values ​​with multiple targets and then perform single-target Bayesian filtering. However, due to the clutter effect in actual tracking scenarios, a large number of real and false measurements are generated, and the computational load of data association increases exponentially, which cannot meet the requirements of real-time tracking.

[0003] In 2003, Mahler proposed the first multi-target tracking algorithm based on random finite set theory, called the Probabilistic Hypothesis Density (PHD) filter. This algorithm uses first-order statistical moments to approximate the posterior probability density of multiple targets, providing a closed-form solution for multi-target Bayesian filtering. Because it avoids complex data correlations, it has attracted widespread attention. Currently, multi-target tracking theory based on the RFS framework is becoming increasingly sophisticated, but many challenges remain. In practical tracking scenarios, the lack of prior information is a major challenge limiting target tracking performance. Prior information includes target emergence intensity, measurement noise, process noise, clutter rate, detection probability, etc. In standard multi-target tracking algorithms, all of this prior information is known, but in practical tracking scenarios, this information is unknown and time-varying. If manually set parameters are used, tracking performance will rapidly decline. Summary of the Invention

[0004] Technical Problem Solved: This invention addresses the performance degradation of traditional multi-target tracking filtering algorithms in scenarios where both measurement noise covariance and detection probability are unknown, leading to inaccurate target tracking. The present invention provides a target tracking method based on adaptive measurement noise and detection probability PHD filtering. This method uses inverse gamma distribution and beta distribution to model the measurement noise covariance matrix and detection probability, respectively. Variational Bayesian iteration is employed to continuously adjust the inverse gamma parameter to approximate the true value of the measurement noise covariance, thereby obtaining the optimal state estimate. Simultaneously, during posterior update, beta distribution parameters are used to model the detection probability and missed detection probability and participate in the weight update, ultimately achieving a joint estimate of the target state, measurement noise covariance, and detection probability.

[0005] Technical solution: The target tracking method with adaptive measurement noise and detection probability PHD filtering proposed in this invention includes the following steps:

[0006] S1, Initialize the target state and tracking scenario, with k = 0 at the initial time:

[0007] S11, initialize the target extended state X = (x, p, R), which can be further written as X = (w k m k P k u k v k α k ,β k ), where x represents the target state, p represents the detection probability, and R represents the measurement noise covariance matrix; w k These are the Gaussian component weights, m k and P k These are the mean and error covariance matrix of the target state, respectively. These three components constitute the traditional PHD components, which are approximations of the first moment u of the target state. k and v k Here are the parameters of the beta distribution, used to model the detection probability: α k and β k These are the degrees of freedom and scale parameters of the inverse gamma distribution, respectively, used to model the measurement noise covariance matrix;

[0008] S12, Initialize basic parameters of the tracking scenario: Set the target state transition matrix F, measurement matrix H, process noise covariance matrix Q, and target survival probability p. s The clutter number λ, attenuation factor ρ, and Gaussian component weight w of newly formed targets within the monitoring area. γ,k ,Target Mean m γ，k The error covariance matrix P of the new target γ，k ;

[0009] S2, using the following formula for the PHD component v in the Gaussian beta inverse gamma mixture form at time k-1. k-1 Predict using (x, p, R):

[0010]

[0011] Among them, v k|k-1 (x,p,R) represents the predicted mixed PHD components, the subscript ξ represents the parameters of the target still alive at time k, and v γ，k (x, p, R) represents the newly generated mixed PHD component, the subscript γ represents the parameter of the newly generated target at time k, and d represents the dimension of the measurement noise covariance matrix.

[0012] The prediction results and the newly generated target PHD component v in step S12 γ，k (x, p, R) are combined to form the final predicted value. The formula can be viewed as three parts: predicting the target state components, beta distribution parameters, and inverse gamma parameters separately, specifically including the following three sub-steps:

[0013] S21, use the standard Kalman filter formula to predict the mean and covariance of the target state at time k-1:

[0014]

[0015]

[0016] S22, using the beta distribution parameters at time k-1, we obtain its mean and variance. Then, using these mean and variance, and the conjugate prior relationship between the beta distribution and the binomial distribution, we predict the beta distribution parameters using the following formula:

[0017]

[0018]

[0019]

[0020]

[0021] in, This represents the predicted mean of the beta distribution. Let k represent the prediction variance of the beta distribution. β The parameter is generally taken to a value greater than 1 to increase the prediction variance;

[0022] S23. The inverse gamma parameter at time k-1 is predicted using a variational Bayesian framework. This parameter is then multiplied by a decay factor ρ to obtain the predicted inverse gamma parameter, specifically using the following formula:

[0023]

[0024]

[0025] S3, the predicted Gaussian-Beta-Gamma mixture PHD component is iteratively updated within a variational Bayesian framework. The update process is divided into two parts: the detected target and the missed target, and specifically includes the following steps:

[0026] S31, For the detected target, based on the measurement set Z obtained at time k... k The predicted extended state set is updated, and variational Bayes iteration is used to obtain the optimal state parameters. The iteration parameters are initialized based on the predicted parameters obtained in step S2. The inverse gamma distribution parameters are iteratively adjusted, and the adjusted inverse gamma parameters are used to model the measurement noise covariance matrix. Then, intermediate parameters such as the innovation covariance and Kalman gain are calculated, followed by further calculation of the adjusted mean, covariance, and weights. During the iteration process, the detection probability is modeled using the beta distribution parameters and participates in updating the weights. After the iteration is complete, the beta distribution parameters are updated. The specific iterative update process can be described using the following formula:

[0027]

[0028]

[0029]

[0030]

[0031]

[0032]

[0033]

[0034]

[0035] After iteration, the beta distribution parameters are updated using the following formula:

[0036]

[0037]

[0038] In the above formula, the superscript n represents the number of iterations. This represents the measurement noise covariance matrix obtained by modeling with the adjusted inverse gamma parameters. Represents the new information covariance matrix. Indicates Kalman gain, Let B(·) represent the j-th measurement value, B(·) represent the beta function, N(·) represent the Gaussian distribution function, and κ represent the kinematic coefficient. k (z) represents the clutter intensity. This represents the mean of the iterative updates. This represents the covariance of the iterative update. This represents the weights updated iteratively. This represents the scale parameter adjusted in the l-th iteration. The scale parameter is continuously adjusted during the iteration process, while the degrees of freedom are initialized only once before the iteration.

[0039] when When the value is basically constant, stop the loop and output. Otherwise, continue the loop;

[0040] S32, for the missed detection target, the predicted beta distribution parameters are used to model the missed detection probability and participate in the weight update. Finally, the beta distribution parameters are updated, while other parameters remain unchanged and consistent with the predicted values. Specifically, the update can be performed using the following formula:

[0041]

[0042]

[0043]

[0044]

[0045] Among them, v ND,k|k (x, p, R) represents the mixed PHD component of the missed detection portion.

[0046] S4, Pruning and Absorption: The Gaussian-Beta-Inverse-Gamma Hybrid PHD components will continuously increase over time. A pruning threshold T is set, and components with weights less than T are directly removed. A merging threshold U is set; if the mean component satisfies... Then, the components corresponding to this mean component are merged and absorbed, specifically using the following formula:

[0047]

[0048]

[0049]

[0050]

[0051]

[0052]

[0053]

[0054] in, This represents the weights after the merger. This represents the combined mean. This represents the combined covariance. and This represents the beta distribution parameters after merging. and These represent the degrees of freedom and scale parameters after merging, respectively.

[0055] S5, State Extraction: The weights of all mixed components are merged and extracted to obtain the final estimated number of targets N. The N components with the largest weights are the final estimated target states.

[0056] Beneficial effects:

[0057] (1) The target tracking method with adaptive measurement noise and detection probability PHD filtering proposed in this invention can effectively track the target and accurately estimate the target base and detection probability.

[0058] (2) The target tracking method with adaptive measurement noise and detection probability PHD filtering proposed in this invention uses inverse gamma distribution and beta distribution to model the measurement noise covariance matrix and detection probability, respectively. Variational Bayesian iteration is used to continuously optimize the inverse gamma parameters to obtain the optimal state estimate. In the posterior update, the beta distribution parameters are used to model the detection and missed detection probabilities and participate in the weight update. This realizes the joint estimation of target state, measurement noise covariance and detection probability. Compared with traditional multi-target tracking algorithms, this invention has better tracking performance in scenarios where both measurement noise covariance and detection probability are unknown. Attached Figure Description

[0059] Figure 1 This is a flowchart of an embodiment of the present invention;

[0060] Figure 2 This is a diagram of the actual motion trajectory of the target in an embodiment of the present invention;

[0061] Figure 3 This is a target state estimation diagram according to an embodiment of the present invention;

[0062] Figure 4 This is a comparison chart of the average GOSPA error of embodiments of the present invention;

[0063] Figure 5 This is a comparison chart of target base estimation in embodiments of the present invention;

[0064] Figure 6 This is a detection probability estimation diagram according to an embodiment of the present invention; Detailed Implementation

[0065] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings:

[0066] like Figure 1 As shown in the figure, the target tracking method with adaptive measurement noise and detection probability PHD filtering proposed in this embodiment of the invention mainly includes the following steps:

[0067] Step 1, initialize the target state and basic parameters of the tracking scene, with k = 0 at the initial time:

[0068] (1a) Initialize the target extended state X = (x, p, R), also written as X = (w k m k P k u k v k α k ,β k ), where x represents the target state, p represents the detection probability, and R represents the measurement noise covariance matrix; w k These are the Gaussian component weights, m k and P k These are the mean and error covariance matrices of the target state, respectively; u k and v k α represents the parameters of the beta distribution, used to model the detection probability. k and β k These are the degrees of freedom and scale parameters of the inverse gamma distribution, respectively, used to model the measurement noise covariance matrix;

[0069] (1b) Initialize the basic parameters of the tracking scenario: Set the target state transition matrix F, measurement matrix H, process noise covariance matrix Q, and target survival probability p. s The clutter number λ, attenuation factor ρ, and Gaussian component weight w of newly formed targets within the monitoring area. γ，k Target mean m for new students γ，k The error covariance matrix P of the new target γ，k ;

[0070] Step 2, use the following formula to calculate the PHD component v in the Gaussian-Beta inverse gamma mixture form at time k-1. k-1 Predict using (x, p, R):

[0071]

[0072] Among them, v k|k-1 (x,p,R) represents the predicted mixed PHD components, the subscript ξ represents the parameters of the target still alive at time k, and v γ，k (x, p, R) represents the newly generated mixed PHD component, the subscript γ represents the parameter of the newly generated target at time k, and d represents the dimension of the measurement noise covariance matrix.

[0073] The prediction results and the newborn target PHD component v in step (1b) γ，k (x, p, R) are combined as the final predicted value, where the parameter prediction can be divided into three parts:

[0074] (2a) Predict the mean and covariance of the target state at time k-1:

[0075]

[0076]

[0077] (2b) Prediction of the beta distribution parameters at time k-1:

[0078]

[0079]

[0080]

[0081]

[0082] in, This represents the predicted mean of the beta distribution. Let k represent the prediction variance of the beta distribution. β The parameter is generally taken to a value greater than 1 to increase the prediction variance;

[0083] (2c) Prediction of the inverse gamma parameter at time k-1:

[0084]

[0085]

[0086] Step 3 involves iteratively updating the predicted Gaussian-Beta-Gamma mixture PHD components within a variational Bayesian framework. The update process is divided into two parts: the detected targets and the missed targets. Specifically, it includes the following steps:

[0087] (3a) For the detected target, based on the measurement set Z obtained at time k k The predicted extended state set is updated by using variational Bayes iteration to obtain the optimal state parameters, and the iteration parameters are initialized based on the predicted parameters obtained in step 2. The optimal estimate is obtained by iteratively adjusting the inverse gamma distribution parameters to approximate the true measurement noise covariance matrix. During the iteration process, the detection probability is modeled using the beta distribution parameters and used to update the weights. After each iteration, the beta distribution parameters are updated. The specific iterative update process can be described by the following formula:

[0088]

[0089]

[0090]

[0091]

[0092]

[0093]

[0094]

[0095]

[0096] After iteration, the beta distribution parameters are updated using the following formula:

[0097]

[0098]

[0099] In the above formula, the superscript n represents the number of iterations. This represents the measurement noise covariance matrix obtained by modeling with the adjusted inverse gamma parameters. Represents the new information covariance matrix. Indicates Kalman gain, Let B(·) represent the j-th measurement value, B(·) represent the beta function, N(·) represent the Gaussian distribution function, and κ represent the kinematic coefficient. k (z) represents the clutter intensity. This represents the mean of the iterative updates. This represents the covariance of the iterative update. This represents the weights updated iteratively. This represents the scale parameter adjusted in the l-th iteration. The scale parameter is continuously adjusted during the iteration process, while the degrees of freedom are initialized only once before the iteration.

[0100] when When the value is basically constant, stop the loop and output. Otherwise, continue the loop;

[0101] (3b) For the missed target, the predicted beta distribution parameters are used to model the missed probability and participate in the weight update. Finally, the beta distribution parameters are updated using the following formula:

[0102]

[0103]

[0104]

[0105]

[0106] Among them, v ND,k|k (x, p, R) represents the mixed PHD component of the missed detection portion.

[0107] Step 4, Pruning and Absorption: The Gaussian-Beta-Inverse Gamma mixture will continuously increase over time. A pruning threshold T is set, and components with weights less than T are directly removed. A merging threshold U is set; if the mean component satisfies... Then, the components corresponding to this mean component are merged and absorbed, specifically using the following formula:

[0108]

[0109]

[0110]

[0111]

[0112]

[0113]

[0114]

[0115] in, This represents the weights after the merger. This represents the combined mean. This represents the combined covariance. and This represents the beta distribution parameters after merging. and These represent the degrees of freedom and scale parameters after merging, respectively.

[0116] Step 5, State Extraction: Merge the weights of all mixed components to extract the final estimated number of targets N. The N components with the largest weights are the final estimated target states.

[0117] The algorithm of this invention has been verified through simulation and has achieved relatively good results:

[0118] 1. Simulation parameter settings:

[0119] The multi-target tracking scenario in the simulation experiment is set in a two-dimensional XY observation area of ​​[0, 300]. There are four targets within the area, moving at a constant velocity in a straight line. The initial positions of all targets are {150, 150}, and their birth and death times are {1, 40}, {1, 70}, {20, 70}, and {20, 81}, respectively. The sampling time interval T = 1s, and the observation time is 81s. The clutter number λ in the monitored area is 10, and the survival probability p of the targets is... s =0.99, detection probability p d =0.95. The initial position of the new target is {100, 100}, and the initial covariance P is 0.95. γ=diag{150 2 1,150 2 ,1}, initial inverse gamma parameter α γ =β γ =[5;5], attenuation factor ρ=1-e -4 The initial beta distribution parameters are u = 8 and v = 1. The state transition matrix F, measurement matrix H, measurement noise covariance matrix R, and process noise covariance matrix Q are set as follows:

[0120]

[0121]

[0122]

[0123]

[0124] This invention employs the Generalized Optimal Subpattern Assignment Metric (GOSPA) as the algorithm's evaluation metric, and further verifies the algorithm's performance through 100 Monte Carlo experiments. The simulation parameters in GOSPA are set as p = 2, c = 10, and α = 2. A control experiment is also set up, where the detection probability p is artificially set under the condition that the measurement noise covariance and detection probability are unknown. d =0.6, measurement noise covariance R = diag{5, 5}, observe the tracking performance differences among the three: the present invention, the PHD algorithm under artificially set mismatched conditions, and the PHD algorithm under matched conditions.

[0125] 2. Simulation Result Analysis

[0126] Figure 2 The actual trajectory of the target. Figure 3 The diagram shows a comparison between the target's actual trajectory and the target state estimation in this embodiment of the invention. The target state estimation represented by the circles basically overlaps with the actual trajectory, indicating that the invention can effectively track the target.

[0127] like Figure 4 The figure shown is a comparison of the average GOSPA error of the embodiments of the present invention. It can be seen that the root mean square error is the largest of the PHD filtering algorithm, which has a mismatch between the measurement noise covariance and the detection probability and the actual value. However, the root mean square error of the embodiment of the present invention (VB-BGM-PHD) is almost the same as that of the conditionally matched PHD filtering algorithm. This shows that the present invention can still track the target relatively accurately when the measurement noise covariance and the detection probability are unknown, and solves the problem of the decline in tracking performance of traditional algorithms when both are unknown.

[0128] like Figure 5The figure shown is a comparison of target number estimation in an embodiment of the present invention. It can be seen that the PHD filtering algorithm with mismatched conditions generally overestimates the number of targets, while the target number estimation in the embodiment of the present invention is almost identical to that of the PHD filtering algorithm with matching conditions, and fits the true number of targets well. Therefore, the present invention can effectively estimate the number of targets.

[0129] like Figure 6 The figure shown is a detection probability estimation diagram of an embodiment of the present invention. It can be seen that the present invention can effectively estimate the detection probability and fits well with the actual detection probability.

[0130] The above description is only a partial embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A target tracking method using adaptive measurement noise and detection probability PHD filtering, characterized in that, The target tracking method includes the following steps: S1, Initialize the target state and tracking scenario, with k = 0 at the initial time: S11, initialize the target extended state X = (x, p, R), which can be further written as X = (w k m k P k u k v k α k ,β k ), where x represents the target state, p represents the detection probability, and R represents the measurement noise covariance matrix; w k These are the Gaussian component weights, m k and P k These are the mean and error covariance matrix of the target state, respectively. These three components constitute the traditional PHD components, i.e., the first-order moment approximation of the target state; u k and v k α represents the parameters of the beta distribution, used to model the detection probability. k and β k These are the degrees of freedom and scale parameters of the inverse gamma distribution, respectively, used to model the measurement noise covariance matrix; S12, Initialize basic parameters of the tracking scenario: Set the target state transition matrix F, measurement matrix H, process noise covariance matrix Q, and target survival probability p. s The clutter number λ, attenuation factor ρ, and Gaussian component weight w of newly formed targets within the monitoring area. γ,k ,Target Mean m γ，k The error covariance matrix P of the new target γ，k ; S2, predicting the target extended state, specifically includes the following steps: S21, use the standard Kalman filter algorithm to predict the mean and covariance of the target state at time k-1, and merge the prediction result with the newly generated target component in step S12 as the predicted Gaussian mixture component. S22, use the beta distribution parameters at time k-1 to obtain its mean and variance, and use the mean and variance as well as the relationship that the beta distribution and the binomial distribution are conjugate priors to make a prediction and obtain the predicted beta distribution parameters. S23, use the variational Bayesian framework to predict the inverse gamma parameter at time k-1, and multiply the parameter by the decay factor ρ to obtain the predicted inverse gamma parameter; S3, the predicted Gaussian-Beta-Gamma mixture components are iteratively updated within a variational Bayesian framework, specifically including the following steps: S31, For the detected target, based on the measurement set Z obtained at time k... k The predicted extended state set is updated, and the variational Bayesian method is used to decouple the state and measurement noise covariance. This process iterates continuously to approximate the true posterior distribution until the state set stabilizes. In each iteration, the predicted beta distribution parameters are used to model the detection probability, which is then used to update the weight components. S32, For the target that is missed, the probability of missed detection is modeled using the predicted beta distribution parameters, which is used to update the weights of the missed part; S4, Pruning and Absorption: The Gaussian-Beta-Gamma mixture will continue to grow over time. Set a pruning threshold T to directly remove components with weights less than the pruning threshold T; set a merging threshold U to merge components with weights less than the merging threshold U. S5, State Extraction: The weights of all mixed components are merged and extracted to obtain the final estimated number of targets N. The N components with the largest weights are the final estimated target states.

2. The target tracking method with adaptive measurement noise and detection probability PHD filtering according to claim 1, characterized in that, In step S2, the following formula is used to evaluate the Gaussian-Beta inverse gamma mixing PHD component v at time k-1. k-1 Predict using (x, p, R): Among them, v k|k-1 (x, p, R) represents the predicted mixed PHD components, the subscript ξ represents the parameters of the target still alive at time k, and v γ，k (x, p, R) represents the newly generated mixed PHD component, the subscript γ represents the parameter of the newly generated target at time k, and d represents the dimension of the measurement noise covariance matrix. Both the PHD components before and after the prediction are in the form of a mixture of Gaussian, beta, and inverse gamma distributions. Essentially, they are predictions of the target state parameters, beta distribution parameters, and inverse gamma distribution parameters. In step S21, the mean and covariance of the target at time k-1 are predicted using the following standard Kalman filter formula: In step S22, the beta distribution parameters at time k-1 are predicted using the following formula: in, This represents the predicted mean of the beta distribution. Let k represent the prediction variance of the beta distribution. β The parameter is generally taken to a value greater than 1 to increase the prediction variance; In step S23, the inverse gamma distribution parameters are predicted using the following formula:

3. The target tracking method with adaptive measurement noise and detection probability PHD filtering according to claim 1, characterized in that, In step S3, the detected targets and missed targets need to be updated respectively, corresponding to steps S31 and S32 respectively: In step S31, for the detected target, variational Bayes iteration is used to obtain the optimal state parameters, and the iteration parameters are initialized based on the prediction parameters obtained in step S2. The inverse gamma distribution parameters are iteratively adjusted, and the adjusted inverse gamma parameters are used to model the measurement noise covariance matrix. Then, intermediate parameters such as the innovation covariance and Kalman gain are calculated, followed by further calculation of the adjusted mean, covariance, and weights. During the iteration process, the detection probability is modeled using the beta distribution parameters and participates in updating the weights. After the iteration is complete, the beta distribution parameters are updated. The specific iterative update process can be described using the following formula: After iteration, the beta distribution parameters are updated using the following formula: In the above formula, the superscript n represents the number of iterations. This represents the measurement noise covariance matrix obtained by modeling with the adjusted inverse gamma parameters. Represents the new information covariance matrix. Indicates Kalman gain, Let B(·) represent the j-th measurement value, B(·) represent the beta function, N(·) represent the Gaussian distribution function, and κ represent the kinematic coefficient. k (z) represents the clutter intensity. This represents the mean of the iterative updates. This represents the covariance of the iterative update. This represents the weights updated iteratively. Let represent the scale parameter adjusted in the l-th iteration. During the iteration process, the scale parameter is continuously adjusted, while the degrees of freedom are initialized only once before each iteration: when When the value is basically constant, stop the loop and output. Otherwise, continue the loop; In step S32, for the missed target, the predicted beta distribution parameters are used to model the missed probability and participate in the weight update. Finally, the beta distribution parameters are updated, while other parameters remain unchanged and consistent with the predicted values. Specifically, the update can be performed using the following formula: Among them, V ND，k|k (x, p, R) represents the mixed PHD component of the missed detection portion.

4. The target tracking method with adaptive measurement noise and detection probability PHD filtering according to claim 1, characterized in that, In step S4, if a merging threshold U is given, and the mean value corresponding to the largest weight component is recorded as... If the mean component satisfy The following formula can then be used for merging and absorption: in, This represents the weights after the merger. This represents the combined mean. This represents the combined covariance. and This represents the beta distribution parameters after merging. and These represent the degrees of freedom and scale parameters after merging, respectively.

Images