Dense multi-unmanned aerial vehicle target real-time tracking method based on GM-PHD filtering

Abstract

The invention relates to the technical field of multi-target tracking, and discloses a dense multi-unmanned aerial vehicle target real-time tracking method based on GM-PHD filtering, and the method comprises the steps: constructing a system state space model for describing a target motion and observation mechanism; executing prediction based on the posterior strength of the previous moment to generate predicted target strength; receiving real-time measurement data, combining with the predicted target intensity and utilizing a weight normalization formula introducing clutter intensity to execute updating, and generating updated posterior intensity; performing Mahalanobis distance-based trimming and merging on the Gaussian mixture items to obtain a refined Gaussian component set; component information is extracted to serve as a state estimation value, state circulation updating is conducted on the free trace point, the experimental track and the formal track through association matching with historical data, and a target track set is output. According to the method, the weight is corrected through clutter intensity, so that the false alarm rate is reduced, the problem of target identity label missing is solved in combination with hierarchical track management, and the tracking precision of dense overlapping targets is improved by using judgment based on the weight.

Description

Technical Field

[0001] This invention relates to the field of multi-target tracking technology, specifically to a method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering. Background Technology

[0002] With the development of UAV swarm technology, real-time detection and tracking of dense multi-target targets has become a key aspect of radar and sensor data processing. Traditional multi-target tracking algorithms, such as joint probabilistic data interconnection algorithms or multi-hypothesis tracking algorithms, are effective in scenarios with a small number of targets and simple environments. However, when faced with complex scenarios with a large number of targets and densely interwoven motion trajectories, the computational load becomes excessive due to the need for complex measurement and trajectory data correlation. The computational dimension increases exponentially with the increase in the number of targets and measurements, making it difficult to meet the real-time requirements of the system.

[0003] To address the computational complexity issue, a Probability Hypothesis Density (PHD) filtering method based on random finite set theory has been proposed. This method replaces the complete multi-objective posterior probability density by passing the first moment of the multi-objective posterior probability density, avoiding the complex data association process in traditional methods. Among them, the Gaussian Mixture Probability Hypothesis Density (GM-PHD) filtering algorithm uses a Gaussian mixture model to recursively approximate the PHD in a closed-form form, and is widely used due to its ease of engineering implementation.

[0004] However, the standard GM-PHD filtering algorithm still has limitations in practical engineering applications for tracking dense UAV swarms. First, the algorithm is highly sensitive to environmental noise during intensity updates. In high-density Poisson clutter, relying solely on detection probability and traditional weight update mechanisms is insufficient to completely eliminate false measurements, leading to clutter accumulation and the formation of false target intensity peaks, thus increasing the system's false alarm rate. Second, the theoretical basis of PHD filtering focuses on estimating the spatial target density at each moment. Its output is essentially a disordered set of points, lacking an inherent target identification maintenance mechanism. This makes it impossible for the algorithm to directly distinguish different targets in the time dimension, making it difficult to form continuous and stable historical tracks. Once a target is obscured or missed, the track is prone to breakage. Furthermore, when UAVs perform dense formation flight or track intersections, due to the close physical distance between targets, the standard pruning and merging steps fuse multiple spatially adjacent Gaussian components into a single component. This causes the algorithm to underestimate the number of targets in overlapping areas, failing to accurately resolve the independent state of each target. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a real-time tracking method for dense multi-UAV targets based on GM-PHD filtering, which solves the problems of traditional probabilistic hypothesis density filtering algorithms in dense target environments, such as difficulty in maintaining target identity continuity, weak clutter suppression capability, and inaccurate multi-target state extraction.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a real-time tracking method for dense multi-UAV targets based on GM-PHD filtering, comprising the following steps: Construct a constant-velocity motion model describing the changes in the target's position and velocity, and a linear Gaussian measurement model describing the sensor's observation mechanism, and generate a system state-space model; Based on the posterior intensity information of the previous moment and the system state space model, intensity prediction is performed to generate the predicted target intensity at the current moment. The system receives real-time measurement data from sensors, which includes the real target and Poisson clutter. It combines the predicted target strength with the strength update using a weighted normalization formula that incorporates clutter strength, and generates an updated posterior strength. The Gaussian mixture terms contained in the updated posterior strength are pruned and merged based on Mahalanobis distance to obtain a refined Gaussian component set; Extract the component information from the refined Gaussian component set as the state estimate, and update the state transition of the free point track, experimental track and formal track by association and matching with historical data, and output the target track set.

[0007] Preferably, the construction of the system state space model specifically includes: defining a four-dimensional target motion state vector containing lateral and longitudinal position and corresponding velocity components for the constant-speed motion model; and establishing a state transition matrix and a process noise vector describing the evolution of the four-dimensional target motion state vector based on the linear assumption of the constant-speed motion model.

[0008] Preferably, the construction of the system state-space model further includes: performing differential processing on the position observation data of the sensor at continuous time intervals for the linear Gaussian measurement model to construct a four-dimensional measurement vector containing position and pseudo-velocity; constructing a measurement noise covariance matrix to describe the observation accuracy of the four-dimensional measurement vector, and setting the main diagonal elements of the noise component submatrix contained in the measurement noise covariance matrix to correspond to the variance of the sensor in the position observation. This setting enables the convergence speed of the target motion trend to be obtained using the pseudo-velocity measurement enhancement filtering algorithm.

[0009] Preferably, generating the predicted target intensity at the current moment specifically includes: extracting the surviving target component from the posterior intensity information of the previous moment; calculating the prediction parameters using the state transition matrix and process noise vector in the system state-space model and performing weight attenuation processing to generate a surviving Gaussian component set; initializing the newly generated Gaussian component according to the boundary of the monitoring area or a preset beam direction to generate a newly generated Gaussian component set; and superimposing and merging the surviving Gaussian component set and the newly generated Gaussian component set to generate the predicted target intensity.

[0010] Preferably, the intensity update specifically includes: combining the prediction covariance in the predicted target intensity with the measurement noise characteristics corresponding to the real-time measurement data to calculate the covariance of the prediction error and obtain an innovative covariance matrix; calculating a Kalman gain matrix for correcting the state prediction value based on the innovative covariance matrix; and using the Kalman gain matrix to calculate and update the covariance matrix.

[0011] Preferably, the step of performing intensity update using the weight normalization formula that incorporates clutter intensity specifically includes: calculating the correlation likelihood and combining it with the detection probability parameter for each measurement vector in the real-time measurement data and each Gaussian component in the predicted target intensity to obtain unnormalized update weights; setting clutter intensity parameters, and based on the mutual exclusion principle of measurement sources, substituting the unnormalized update weights and the clutter intensity into the weight normalization formula to calculate normalized update weights.

[0012] Preferably, obtaining the refined Gaussian component set specifically includes: traversing all Gaussian components contained in the updated posterior strength, setting a truncation threshold, removing components with weights lower than the truncation threshold, and obtaining a remaining component set; selecting the principal component with the largest weight in the remaining component set, calculating the Mahalanobis distance between the other components in the remaining component set and the principal component; performing a weighted average calculation on the components with Mahalanobis distances less than a merging threshold to fuse them into new components, and using the new components to constitute the refined Gaussian component set. Wherein, the truncation threshold is a preset lower weight limit value used to filter out low-weight invalid components; the merging threshold is a preset upper distance limit value used to determine whether two Gaussian components belong to the same target.

[0013] Preferably, the extraction of the state estimate specifically includes: traversing the refined Gaussian component set, filtering components whose weights are greater than or equal to an extraction threshold, and adding the corresponding state mean to the estimated target state set; if the weight of the target component in the refined Gaussian component set indicates that there are multiple targets at the location corresponding to the target component, extracting a corresponding number of overlapping points according to the rounding principle and adding them to the estimated target state set, and using the points in the estimated target state set as the state estimate. Wherein, the extraction threshold is a preset weight probability value used to determine whether a target exists.

[0014] Preferably, the state transition update for free tracks, experimental tracks, and formal tracks specifically includes: acquiring the state estimate and the free tracks, experimental tracks, and formal tracks established at historical times; calculating the Euclidean distance between the state estimate and the predicted position of the existing track; if the Euclidean distance is less than the association gate, the association is considered successful; determining that the state estimate that fails to be successfully associated with any existing track is a free track; upgrading the state estimate to which the free track from the previous time step is successfully associated at the current time step is an experimental track; when the accumulated value of the hit counter of the experimental track reaches the initial confirmation frame number, upgrading the experimental track to a formal track. The association gate is a preset distance threshold used to determine whether the observed track and the track belong to the same target.

[0015] Preferably, updating the state transition of the formal track further includes: for the established formal track, if the association is successful at the current time, updating the state of the formal track using Kalman filtering; if the formal track is not associated with any point, performing a pure prediction update and accumulating the loss counter; when the loss counter value of the formal track exceeds the maximum loss deletion threshold, deleting the formal track. The maximum loss deletion threshold is a preset maximum number of frames that a track can remain undetected for consecutive periods.

[0016] This invention provides a method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering. It has the following advantages: 1. This invention adopts a weighted normalization formula that incorporates clutter intensity during the intensity update stage. When calculating the Gaussian component weight, clutter intensity is used as a denominator term in the calculation. By utilizing the mutual exclusivity principle of measurement sources, the posterior weight of false measurements generated in high-density Poisson clutter regions is reduced. This effectively distinguishes real targets from environmental noise under dense clutter interference, thereby reducing the false alarm rate during the tracking process.

[0017] 2. This invention combines the intensity estimation of random finite set theory with hierarchical track management logic, mapping the disordered points output by filtering to a three-level circulation system of free points, experimental tracks and formal tracks. By setting the initial confirmation frame number and the loss deletion threshold, it eliminates instantaneous interference and continuously associates stable targets, solving the problem that traditional probability hypothesis density filtering algorithms cannot provide target identification, and ensuring the continuity of multi-UAV target tracks.

[0018] 3. In the state extraction stage, this invention adopts a weight-based rounding decision principle. By utilizing the characteristic that the sum of the Gaussian component weights corresponds to the expected number of targets, when the weight of a single Gaussian component indicates the existence of multiple overlapping targets, the corresponding number of overlapping points can be directly parsed out. This enables accurate extraction of the state of multiple targets in dense formations or trajectory intersection scenarios, thereby improving the ability to analyze and distinguish overlapping targets. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the system functional modules according to an embodiment of the present invention; Figure 2 This is the main flowchart of the tracking method according to an embodiment of the present invention; Figure 3 This is a flowchart of the Gaussian component trimming and merging sub-process of an embodiment of the present invention; Figure 4 This is a comparison chart of the simulated and predicted positions in an embodiment of the present invention; Figure 5 This is a diagram illustrating the trajectory tracking effect of an embodiment of the present invention. Figure 6 This is a distribution diagram of the original measurement data on the X-axis in an embodiment of the present invention; Figure 7 This is a distribution diagram of the original measurement data on the Y-axis in an embodiment of the present invention; Figure 8 This is a graph showing the OSPA distance evaluation results according to an embodiment of the present invention; Figure 9 This is a comparison chart of multi-objective potential estimation in an embodiment of the present invention.

[0020] Among them, 10 is the parameter configuration module; 20 is the filtering processing module; 30 is the track management module; and 40 is the visualization evaluation module. Detailed Implementation

[0021] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] Reference Figure 1 and Figure 2This invention provides a method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering. This method is executed by a computer system or embedded processing device, which includes a parameter configuration module 10, a filtering module 20, a trajectory management module 30, and a visualization evaluation module 40. The parameter configuration module 10 is used to set environmental and sensor parameters and generate simulation scene data; the filtering module 20 is configured to perform the core operation of the Gaussian mixture probability hypothesis density filtering algorithm; the trajectory management module 30 is used to handle the association between point traces and trajectories and their lifecycle management; and the visualization evaluation module 40 is used to output trajectory results and performance indicators.

[0023] This invention provides a method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering, comprising the following steps: S1. Construct the target motion model and measurement model. For dense UAV swarms, establish a constant-speed motion model describing the changes in target position and velocity, and a linear Gaussian measurement model describing the sensor observation mechanism. Define the target's state vector, state transition matrix, measurement matrix, and related noise covariance matrix to construct the system state space model. S2, Execute intensity prediction; The filtering module 20 performs Kalman prediction on surviving targets and adaptive intensity initialization on new targets based on the posterior intensity information of the previous time step and the system state space model, calculates the state mean and covariance, and generates the predicted target intensity at the current time step. S3, perform intensity update; receive real-time measurement data including real target and Poisson clutter collected by the sensor, combine the predicted target intensity to calculate the predicted observation value, innovative covariance and Kalman gain, use the weight normalization formula introduced by the clutter intensity to correct the weight of each Gaussian component, and generate the updated posterior intensity. S4, perform Gaussian component pruning and merging; set truncation threshold and merging threshold, perform low-weight removal and Mahalanobis distance-based merging operations on the Gaussian mixture terms contained in the updated posterior strength, count the number of remaining components, and obtain the refined Gaussian component set. S5, execute state extraction and track management; extract the mean of the components whose weights meet the extraction threshold in the refined Gaussian component set as the state estimate, and input the state estimate into the track management module 30. Through association and matching with historical data, the state transition of free track, experimental track and formal track is updated, and the smoothed target track set is output. S6, to perform performance evaluation; the visualization evaluation module 40 calculates the optimal sub-mode allocation distance and potential estimation error based on the output target track set and preset real trajectory, and quantifies the convergence and robustness of the tracking algorithm in dense clutter environment.

[0024] The specific implementation details of steps S1 to S6 above will be elaborated below.

[0025] Reference Figure 1 In step S1, the parameter configuration module 10 initializes the system state-space model based on preset UAV swarm characteristics and sensor performance indicators. This initialization of the system state-space model forms the mathematical basis for subsequent state estimation by the filtering module 20. For the motion characteristics of dense UAV swarms maintaining formation or stable flight during the cruise or penetration phase, this embodiment uses a constant-speed motion model to describe the target dynamics and constructs a sensor measurement model based on linear assumptions. The specific model construction process includes the following sub-steps: S11, Define the target motion state vector and state transition equation. To fully characterize the target's motion trend in the two-dimensional plane, the lateral and longitudinal positions and their corresponding velocity components are selected to constitute the system state. A four-dimensional state vector is constructed using the state vector definition formula: ; in, Represents the target motion state vector. Indicates the horizontal position coordinates. Indicates the vertical position coordinates. Represents the lateral velocity component. This represents the longitudinal velocity component.

[0026] Assuming the target moves at a constant velocity in a straight line for a short period, and the state evolution conforms to the characteristics of a linear time-invariant system, the current state is obtained by linearly predicting the state at the previous time step and adding process noise. The process is described using a linear state transition formula: ; in, express The target motion state vector at time t. express The target motion state vector at time t. Represents the process noise vector. This represents the state transition matrix.

[0027] S12, set the state transition matrix and process noise parameters.

[0028] The state transition matrix describes the coupling relationship between position and velocity over time. Based on the constant velocity model assumption, the state transition matrix is ​​defined by the following formula: ; in, Represents a second-order identity matrix. Represents a second-order zero matrix. This represents the time sampling interval. In this embodiment, the time sampling interval is set to 1 second. Considering the impact of airflow disturbances or slight maneuvers on the target during actual movement, a process noise vector following a zero-mean Gaussian distribution is introduced, i.e. The noise distribution is determined using the formula for calculating the process noise covariance matrix: ; in, Represents the process noise covariance matrix. This represents the standard deviation of process noise. This reflects the degree of uncertainty in the target's acceleration. In the simulation scenario of this embodiment, the standard deviation of the process noise is set to 0.8.

[0029] S13, Constructing a linear measurement equation. To improve the convergence of the filtering algorithm for target velocity estimation, this embodiment constructs a four-dimensional measurement vector containing position and pseudo-velocity by performing differential processing on the position observation data of the sensor at continuous time points. The measurement process is described using a linear measurement equation: ; in, express The measurement vector at time includes the observed values ​​of the horizontal and vertical coordinate positions and the observed values ​​of the velocity obtained by differential calculation; for The target motion state vector at any given time; Represents the measurement noise vector; The measurement matrix is ​​used to describe the mapping from the state space to the observation space. Since the dimensions of the constructed measurement vectors are the same as those of the state vectors, the measurement matrix is ​​set to a fourth-order identity matrix.

[0030] S14, construct the measurement noise covariance matrix.

[0031] The measurement noise vector follows a zero-mean Gaussian distribution. To describe the sensor's observation accuracy in each dimension, a measurement noise covariance matrix is ​​constructed. The formula for constructing the measurement noise covariance matrix is ​​as follows: ; in, Represents the measurement noise covariance matrix; The noise component submatrix is ​​a product of a second-order identity matrix and the observation variance. In this embodiment, the main diagonal elements of the noise component submatrix are set to 100, which corresponds to the variance of the sensor's position observations and reflects the sensor's measurement accuracy level. The zero matrix represents a second-order zero matrix, indicating that the position measurement noise and velocity measurement noise are statistically independent.

[0032] In step S2, the filtering module 20 performs intensity prediction. This step, based on the Bayesian filtering principle, extrapolates the posterior multi-target intensity from the previous moment along the time axis. The basic principle is that the spatial target intensity distribution at the current moment is composed of two superimposed parts: one part is the intensity of surviving targets from the previous moment after migration, and the other part is the intensity of newly formed targets entering the monitoring range at the current moment. The specific prediction process includes the following sub-steps: S21, calculate the Gaussian component prediction parameters of the surviving targets.

[0033] The filtering module iterates through the previous time step 20 times. The output is a posterior Gaussian mixture intensity set, which consists of multiple Gaussian components, each containing weights, mean, and covariance. For each Gaussian component in the set, a one-step state prediction is performed using the time update equation of the Kalman filter.

[0034] The expected position of each component at the current time is calculated using the state prediction mean formula: ; in, express The mean of the predicted state at time t, Represents the state transition matrix. express The mean of the posterior state at time t.

[0035] The uncertainty diffusion range of each component is calculated using the formula for the predictive covariance matrix: ; in, express The prediction covariance matrix at time 1, Represents the process noise covariance matrix. express The posterior covariance matrix at time t, with superscript This represents the matrix transpose operation.

[0036] Simultaneously, the weights of the survival components are attenuated. A survival probability parameter is set, which represents the probability that the target has not disappeared from the previous moment to the current moment and is still within the field of view. In this embodiment, the survival probability is set to 0.99. The updated prediction weight is calculated as the product of the previous weight and the survival probability.

[0037] S22, construct the Gaussian components of the new target and generate the overall predicted intensity.

[0038] To capture newly emerging targets within the field of view, the filtering module 20 employs a Gaussian mixture model to construct the newborn intensity. Based on the boundary of the monitoring area or a preset beam orientation, a set of newborn Gaussian components is initialized.

[0039] For each newly generated Gaussian component, the initial parameters are set as follows: Newborn weight: set to 0.03, representing the prior probability density of a new target appearing at this position; Average value of newborn status: set as a preset coordinate point of the monitoring area boundary or a coordinate point determined based on prior information; The newborn covariance matrix: is set as a diagonal matrix, for example... , representing the initial uncertainty range of the position and velocity of the new target.

[0040] Finally, the filtering module 20 superimposes and merges the surviving Gaussian component set calculated in step S21 with the newly formed Gaussian component set constructed in this step to form the predicted target intensity at the current moment. This predicted target intensity contains the spatial distribution information of all potential targets (surviving and newly formed), and serves as the prior input for intensity updating in step S3.

[0041] In step S3, the filtering module 20 performs an intensity update. This step uses real-time measurement data acquired by the sensor to correct the predicted target intensity generated in step S2. Its core logic is to calculate the posterior probability weight of each measurement vector originating from the real target or background clutter. The specific update process includes the following sub-steps: S31, calculate the updated Kalman gain and state covariance. Filtering module 20 obtains the current time step. Given a set of sensor measurements, for each measurement vector in the set, calculate its statistical distance to each Gaussian component in the predicted target intensity.

[0042] First, calculate the covariance matrix of the prediction error. The innovative covariance calculation formula is as follows: ; in, Represent the innovation covariance matrix; Represents the measurement noise covariance matrix; Represents the measurement matrix; express The predictive covariance matrix at time t; superscript This represents the matrix transpose operation.

[0043] Based on this, the gain matrix used to correct the predicted state values ​​is calculated. The Kalman gain calculation formula is as follows: ; in, Represents the Kalman gain matrix; This represents the inverse matrix of the innovation covariance matrix.

[0044] By correcting the prediction covariance using the Kalman gain matrix, a more convergent state distribution shape is obtained. The state covariance update formula is as follows: ; in, express The covariance matrix is ​​updated at each time step; This represents the fourth-order identity matrix corresponding to the dimension of the state vector.

[0045] Simultaneously, the filtering module 20 calculates the updated mean of each Gaussian component, using the following calculation method: The mean of the predicted state at time step is added to the product of the Kalman gain matrix and the innovation vector (i.e., the measurement vector minus the predicted measurement vector).

[0046] S32, calculate the unnormalized update weights for each Gaussian component.

[0047] For the first in the measurement set Measurement vectors The first in the predicted target intensity The correlation likelihood between the two is calculated using Gaussian components. A detection probability parameter is set to characterize the probability that the real target is detected by the sensor; in this embodiment, the detection probability is set to 0.98. The Gaussian component weight update formula is used to calculate: ; in, This indicates that the update weights have not been normalized. Indicates the detection probability; Indicates the first Predicted weights for each component; Let represent the Gaussian likelihood function. The Gaussian likelihood function is used to calculate the _th ... The measurement vector falls within the first... Each component The predicted state at time t is the probability density value within a Gaussian distribution centered at the mean and with the innovation covariance matrix as the covariance.

[0048] S33 introduces clutter intensity and performs weighted normalization.

[0049] To eliminate spurious signals, the filtering module 20 introduces a clutter model to correct the weights. A clutter intensity parameter is set, its value equal to the average number of clutter particles in the monitoring area divided by the total area of ​​the monitoring area. In the simulation environment of this embodiment, the clutter intensity is set to... (corresponding to) An average of 50 clutter points exist within a square meter area.

[0050] Based on the principle of mutual exclusion of measurement sources (i.e., the measurement either comes from clutter or from the target), the final weight is calculated using a weight normalization formula: ; in, This indicates that the weights are updated using normalization. Represents measurement vector Clutter intensity at the location; For summation index variables; This represents the total number of predicted Gaussian components; : Represents the clutter intensity parameter, the value of which is equal to the average number of clutter particles in the monitoring area divided by the total area of ​​the monitoring area; Represents all predicted Gaussian components relative to the first... The sum of the correlation contributions of each measurement vector; Indicates the first The predicted Gaussian component is for the first... Unnormalized update weights for each measurement vector; : Represents the summation operator, used to calculate sums from the first digit to the last digit. One to the first The sum of the weights of each component; This represents the normalized denominator term, that is, the term in the 1st division. The total intensity density at each measurement vector location encompasses the sum of intensity contributions from clutter sources and from all potential real targets.

[0051] This formula allows for the denominator to be increased when the clutter intensity at a certain location is high, thereby suppressing the normalized update weights generated at that location and effectively filtering out clutter interference. Finally, the filtering module 20 outputs a set containing all update components (components obtained for each measurement update) and missed detection components (i.e., predicted components multiplied by the missed detection probability), which serves as the updated posterior strength.

[0052] Reference Figure 3 In step S4, since the updated posterior strength output in step S3 is generated based on the combination of each measurement and each predicted component, the number of Gaussian components increases exponentially over time, leading to a surge in computational load. Therefore, the filtering module 20 performs adaptive pruning and merging of the Gaussian components, aiming to control the computational load and extract the single-target peak value while preserving effective target information. This process corresponds to... Figure 3 The logic shown, in its specific implementation process, includes the following sub-steps: S41, perform low-weight component truncation. Filtering module 20 executes. Figure 3The operation of "pruning targets with weights below the threshold" is described. A truncation threshold is set, defining the minimum weight lower bound that the system considers a Gaussian component to be valid. In this embodiment, the truncation threshold is set to 10. -5 The system iterates through the Gaussian component set output in step S3 and checks the weight value of each component. Components with weights below the truncation threshold are treated as noise or extremely low-probability assumptions and directly discarded, while only components whose weights meet the criteria are retained for subsequent processing.

[0053] S42, perform nearest neighbor search based on Mahalanobis distance. For the remaining component set after truncation, the filtering module 20 performs... Figure 3 The process involves initializing and merging related variables, and iterating through the target. First, the Gaussian component with the largest current weight is identified as the principal component. Then, the Mahalanobis distance between each other component and the principal component is calculated. Mahalanobis distance is a statistical distance that considers the covariance distribution characteristics of Gaussian components and effectively measures the degree of overlap between the distributions of two random variables. A merging threshold is set; in this embodiment, it is set to 4 (this value is determined based on the statistical characteristics of the chi-square distribution, representing the statistical similarity tolerance between components). If the calculated Mahalanobis distance is less than the merging threshold, the component and the principal component are considered to describe the same target and are marked as a group to be merged; if the distance is greater than or equal to the merging threshold, no marking is performed.

[0054] S43, perform the merging calculation of Gaussian components. For multiple Gaussian components marked as the group to be merged in step S42, the filtering module 20 performs... Figure 3 The operation of "merging similar targets" merges them into a new Gaussian component. The specific calculation logic is as follows: Merged weights: The sum of the weights of all components in the group is used as the probability of the existence of the new component. The mean state after merging is calculated by weighting the mean state of each component in the group using the weight of each component as a weighting factor, and the center position of the new component is obtained. The merged covariance matrix is ​​obtained by calculating the weighted average of the covariance matrices of each component within the group, and then adding the dispersion of each component mean relative to the merged state mean (i.e., the weighted sum of the outer product terms), thus obtaining a new covariance matrix that includes the uncertainty range of all sub-components. Through the above steps, multiple loosely defined assumptions describing the same objective are coalesced into a compact unimodal Gaussian distribution.

[0055] S44, perform component quantity control. After completing one round of merging, add the newly generated merged components to the result list and remove the components that participated in the merging from the original set. Repeat steps S42 to S43 until the original set is empty. Finally, the filtering module 20 executes... Figure 3The operation of "calculating the number of pruned components" is described. A maximum number of Gaussian components is set, which is 100 in this embodiment. If the final total number of merged components exceeds this limit, the system will sort the components in descending order according to their weights, retaining only the top 100 components with the highest weights and discarding the rest. The refined Gaussian component set obtained after the above processing is the optimal approximate representation of the multi-target intensity at the current moment.

[0056] In step S5, the track management module 30 performs state extraction and hierarchical track management. Since the refined Gaussian component set output in step S4 is essentially still a spatial probability density distribution, lacking a continuous temporal identifier (i.e., no track ID), and may contain residual strong clutter interference, this embodiment constructs a hard-decision mechanism based on an intensity threshold, combined with three-level lifecycle management logic, aiming to extract definite target points from the probability field and establish stable tracks. The specific implementation process includes the following sub-steps: S51, Target state estimation based on extraction threshold. The trajectory management module 30 receives the refined Gaussian component set containing several Gaussian components output from step S4. According to the theory of random finite sets, the sum of the weights of the Gaussian components corresponds to the expected number of targets within the monitoring area. To extract the determined target position, a threshold decision needs to be made on the weights. The effective components are screened using the state target extraction formula: ; in, express The estimated target state set at time, i.e., the final extracted list of measurement points; express The first moment The state mean of a Gaussian component, that is, the center coordinate of the component in the state space; express The first moment The weights of each Gaussian component are determined. In this embodiment, the extraction threshold is set to 0.5. The system iterates through each Gaussian component in the set, and if its weight... Then the corresponding state mean join in The estimated set of target states at each time step. Furthermore, considering the potential target overlap in dense clusters, if the weight of a certain component... (This means that there may be multiple targets at this location), and the system uses rounding principles (i.e., Extract the corresponding number of overlapping points, and take the coordinate values ​​of the components as the state mean.

[0057] S52 performs three-level flow management of free tracks, experimental tracks, and official tracks. The track management module 30 maintains a dynamic track list and is based on... The points in the estimated target state set at each time point are correlated and matched with existing trajectories. The nearest neighbor algorithm is used to calculate the Euclidean distance between the predicted positions of the points and the trajectories. An association gate is set, which is set to 50 meters in this embodiment (determined based on the UAV's maximum maneuver speed and sampling interval). If the distance is less than the association gate, the association is considered successful. Based on the association results, the system executes the following hierarchical flow logic: Free Tracking (Level 1) Management: If If a point in the estimated target state set at a given time fails to be successfully associated with any existing track (including experimental and official tracks), it is determined to be a free point. The system assigns a temporary ID to this point and initializes the hit counter to 1. Free points represent potential new targets or transient clutter and are not output to the outside world at this time.

[0058] Experimental Track (Level 2) Management: If a free track from the previous time step successfully associates with a new track in the current time step, it is upgraded to an experimental track. The system continuously monitors the continuous association status of experimental tracks. An initial confirmation frame count is set; in this embodiment, it is set to 3. If the hit counter of the experimental track reaches the initial confirmation frame count, it indicates that the target has a stable motion pattern, eliminating the possibility of clutter interference, and it is upgraded to a formal track.

[0059] Formal Track (Level 3) Management: Formal tracks are published externally through the visualization evaluation module 40. For established formal tracks, if a successful association is achieved at the current time, the state is updated using Kalman filtering, and the loss counter is reset to 0; if no points are associated (i.e., a missed detection or target occlusion occurs), a pure prediction update is performed using the state transition matrix in step S1 to maintain track projection, and the loss counter is incremented by 1. A maximum loss deletion threshold is set, which is set to 5 in this embodiment. If the loss counter value of a formal track exceeds this threshold, the system determines that the target has disappeared or flown out of the monitoring area, permanently deletes the track from the list, and reclaims storage resources.

[0060] After completing the association and status transition update of all points and tracks at the current moment, the track management module 30 summarizes the target information of all targets marked "official track" in the current list and constructs a target track set. This target track set, as the final output of this step, contains the unique track ID, position coordinates, and velocity estimate of each confirmed target, and is passed to the visualization evaluation module 40 for subsequent processing.

[0061] To verify the effectiveness of the method proposed in this embodiment of the invention, a multi-target simulation scenario under a high clutter background was constructed, and the output data was analyzed in detail. The following, in conjunction with the appendix... Figure 4 - Appendix Figure 7The simulation environment settings, the distribution characteristics of the raw data, and the consistency results of trajectory tracking are described.

[0062] S61, Construct a high-clutter simulation environment and analyze the distribution of raw measurement data. This step is performed by parameter configuration module 10, which aims to generate test data that conforms to the characteristics of real radar or sonar detection. Parameter configuration module 10 first sets the basic parameters of the simulation scene: the monitoring area is set to a two-dimensional plane. (Unit: meters), the total simulation time is 100 seconds. Three real targets moving in uniform linear motion are preset in the scene: target 1 and target 2 appear simultaneously at time 0 and their trajectories intersect near the coordinates (450, -400); target 3 enters the monitoring area at time 20 to simulate the scene of a newly emerging target.

[0063] To simulate harsh detection environments, the parameter configuration module 10 generates clutter data based on a Poisson distribution. The clutter intensity parameters are then set. for This value means that, on average, there are approximately 50 spurious measurement points within the field of view in each scan. Simultaneously, the sensor's detection probability is set. The value is set to 0.98, and the standard deviation of the measurement noise is set to 10 meters. This parameter is used to construct the measurement noise covariance matrix in the formula. .

[0064] See attached document Figure 6 and attached Figure 7 By using the above markings, the relationship between the actual measurement and the filtered output under high-density clutter background can be intuitively compared.

[0065] From the data distribution characteristics, the trajectory signal of the real target is submerged in a large amount of environmental noise, resulting in an extremely low signal-to-noise ratio per frame. In particular, when target 3 appears in the 20th frame, its initial signal is spatially indistinguishable from clutter. This indicates that signal strength detection based on a single frame alone cannot effectively extract the target. It is necessary to adopt the intensity update algorithm that introduces clutter intensity in step S3 of this embodiment. This algorithm utilizes the spatiotemporal continuity and correlation of multi-frame data and suppresses discretely distributed clutter (i.e., clutter / all measurements (gray ×)) through a weight normalization formula, thereby enhancing the posterior probability weight of the target signal (i.e., the real target measurement (white circle)).

[0066] S62 verifies the accuracy of trajectory tracking and the consistency of position.

[0067] After completing the filter update and track management in step S5, the visualization evaluation module 40 receives the target track set output in step S5 and compares and analyzes it with the preset real trajectory data in the spatiotemporal dimension.

[0068] See attached document Figure 4In this embodiment, the predicted position output by the algorithm closely overlaps with the actual trajectory line, indicating that the filtered state estimation has extremely high positional accuracy. Particularly in the intersection region near coordinates (450, -400), the two trajectories intersect in an "X" shape. In this complex region, thanks to the Gaussian component weight update formula in step S3, which independently calculates the likelihood of each Gaussian component, the algorithm does not exhibit target identity swapping or trajectory divergence; the GM-PHD estimated trajectory still maintains good linear continuity.

[0069] Figure 5 To highlight the tracking effectiveness of valid tracks and improve readability, candidate points / short track segments that did not meet the track confirmation criteria have been removed during the track management process. Figure 5 The text only shows the comparison between the final valid trajectory and the corresponding true trajectory.

[0070] Experimental tracks or free points determined as invalid by track management module 30 (corresponding to...) Figure 6 The clutter points that occasionally appear consecutively in the path (though temporarily associated as candidate targets by the system) are ultimately identified as false paths and removed by the path management logic because they fail to meet the initial confirmation frame number (3 frames in this embodiment) or the continuous association condition set in step S5. Figure 5 The result is no longer displayed. This result visually demonstrates the engineering practicality of the three-level flow mechanism, which includes free track, experimental track, and official track, in suppressing false tracks and ensuring the purity of the final output track.

[0071] See attached document Figure 8 and appendix Figure 9 S63, perform a comprehensive error assessment based on the optimal sub-pattern allocation distance.

[0072] To quantitatively verify the overall performance of the tracking method proposed in this embodiment in a dense clutter environment, the visualization evaluation module 40 uses the target track set output in step S5 as the estimation set and the real target state in the simulated scenario as the ground truth set to calculate the Optimal Sub-pattern Assignment (OSPA) distance between the two. The OSPA distance, as a core evaluation metric, can simultaneously measure the position error of the target track set relative to the real trajectory and the cardinality (i.e., the number of targets) estimation error.

[0073] See attached document Figure 8 , attached Figure 8 The horizontal axis “simulation step” represents the number of frames advanced on the simulation timeline, and the vertical axis “OSPA distance” represents the overall error value. The lower the value, the closer the algorithm’s estimation result is to the real scene.

[0074] In the initial stage of the simulation (frames 0 to 10), the OSPA distance remained at a high level of approximately 200 because the filter had not yet fully established stable tracking of the surviving target and the clutter density was high.

[0075] Entering the interval from frame 10 to frame 20, as the weighted normalization formula for clutter intensity introduced in step S3 continues to take effect, the weight of clutter components is effectively suppressed, the posterior probability of the target component is significantly improved, and the OSPA distance shows a sharp downward trend, indicating that the algorithm has successfully captured the target features. The peaks that appear during this period correspond to the unstable state in the early stage of track establishment.

[0076] Around frame 20, the entry of a third newly formed target into the monitoring area disrupted the previous steady state, causing a brief jump in the OSPA distance. However, thanks to the Gaussian mixture model used in step S2, the algorithm quickly acquired the newly formed target in the following frames, and the OSPA distance subsequently dropped back and stabilized within the range of 75 to 80.

[0077] During the stable tracking phase from frame 20 to frame 100, despite the persistent background clutter, the OSPA distance remained at a low level and showed no divergence trend. Figure 8 The occasional spike pulses indicate that even if a false alarm or missed detection occurs at an individual moment, the track management mechanism in step S5 can quickly correct the state estimate through the correlation and matching of historical information, preventing the accumulation and spread of errors. This demonstrates that the method in this embodiment has good convergence and robustness in complex environments.

[0078] S64, perform accuracy verification of multi-objective potential estimation.

[0079] Potential estimation reflects the accuracy of the algorithm in judging the number of targets within the field of view and is a key performance indicator of multi-target tracking systems.

[0080] See attached document Figure 9 , attached Figure 9 The horizontal axis, "Number of Iterations," corresponds to the number of time frames, while the vertical axis, "Number of Targets," represents the statistical value of the number of targets within the field of view. Specifically, the red dashed squares represent the actual number of targets set in the simulation scene (2 targets for frames 0-20, 3 targets for frames 20-100), and the blue solid circles represent the estimated number of targets output by the algorithm in this embodiment.

[0081] from Figure 9 It can be seen that in the initial stage of frames 0 to 5, the estimated number of targets gradually increases from 0. This is because step S5 sets the initial confirmation frame number (3 frames in this embodiment). The system needs to accumulate enough confidence to upgrade the experimental track to a formal track. This lag is in line with the design expectation of the track management logic.

[0082] At frame 20, the number of real targets abruptly changes from 2 to 3. The blue solid line follows this change almost synchronously, demonstrating that the newly constructed Gaussian component set in step S2 can sensitively cover the boundary of the monitoring area and respond promptly to the emergence of new targets.

[0083] Throughout the simulation, the blue solid line and the red dashed line largely overlapped. Although the estimated curve experienced momentary downward fluctuations near frames 30, 65, and 90 (i.e., the estimated number was briefly less than the actual number), this was typically caused by a single-frame missed detection due to strong clutter obscuring the target signal. However, thanks to the maximum loss deletion threshold set in step S5 (5 frames in this embodiment), the system did not immediately delete tracks due to a single-frame detection failure, but instead performed a pure prediction update. Therefore, the estimated curve quickly recovered to the true value level, without permanent loss of target identity or track breakage. This demonstrates that this embodiment effectively solves the common problem of unstable target number estimation in dense multi-UAV tracking by combining the probabilistic update advantages of GM-PHD filtering with deterministic track lifecycle management.

Claims

1. A method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering, characterized in that, Includes the following steps: Construct a constant-velocity motion model describing the changes in the target's position and velocity, and a linear Gaussian measurement model describing the sensor's observation mechanism, and generate a system state-space model; Based on the posterior intensity information of the previous moment and the system state space model, intensity prediction is performed to generate the predicted target intensity at the current moment. The system receives real-time measurement data from sensors, which includes the real target and Poisson clutter. It combines the predicted target strength with the strength update using a weighted normalization formula that incorporates clutter strength, and generates an updated posterior strength. The Gaussian mixture terms contained in the updated posterior strength are pruned and merged based on Mahalanobis distance to obtain a refined Gaussian component set; Extract the component information from the refined Gaussian component set as the state estimate, and update the state transition of the free point track, experimental track and formal track by association and matching with historical data, and output the target track set.

2. The method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The generated system state-space model includes: For the constant-speed motion model, a four-dimensional target motion state vector is defined, which includes the horizontal and vertical positions and their corresponding velocity components; Based on the linear assumption of the constant-speed motion model, a state transition matrix and a process noise vector are established to describe the evolution of the four-dimensional target motion state vector.

3. The method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The generated system state-space model also includes: For the linear Gaussian measurement model, the position observation data of the sensor at continuous time points are differentially processed to construct a four-dimensional measurement vector containing position and pseudo-velocity. A measurement noise covariance matrix is ​​constructed to describe the observation accuracy of the four-dimensional measurement vector. The main diagonal elements of the noise component submatrix contained in the measurement noise covariance matrix correspond to the variance of the sensor in the position observation.

4. The method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The generation of the predicted target intensity at the current moment includes: The surviving target component is extracted from the posterior intensity information of the previous time step. The prediction parameters are calculated using the state transition matrix and process noise vector in the system state space model and weight decay processing is performed to generate a set of surviving Gaussian components. The newly generated Gaussian components are initialized according to the boundary of the monitoring area or the preset bundle direction, and a set of newly generated Gaussian components is generated. The surviving Gaussian component set and the newly formed Gaussian component set are superimposed and merged to generate the predicted target intensity.

5. A method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The execution strength update also includes: By combining the prediction covariance in the predicted target intensity with the measurement noise characteristics corresponding to the real-time measurement data, the covariance of the prediction error is calculated to obtain the innovative covariance matrix; Calculate the Kalman gain matrix for correcting the state prediction value based on the innovative covariance matrix; The updated covariance matrix is ​​calculated using the Kalman gain matrix.

6. The method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The intensity update using the weighted normalization formula that incorporates clutter intensity includes: For each measurement vector in the real-time measurement data and each Gaussian component in the predicted target intensity, the correlation likelihood is calculated and combined with the detection probability parameter to obtain the unnormalized update weight; Set clutter intensity parameters, and based on the mutual exclusion principle of measurement sources, substitute the unnormalized update weights and the clutter intensity into the weight normalization formula to calculate the normalized update weights.

7. The method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The refined Gaussian component set obtained includes: Traverse all Gaussian components contained in the updated posterior strength, set a truncation threshold, and remove components with weights lower than the truncation threshold to obtain the set of remaining components; Select the principal component with the largest weight from the set of remaining components, and calculate the Mahalanobis distance between the other components in the set of remaining components and the principal component; The components whose Mahalanobis distance is less than the merging threshold are weighted and averaged to generate new components, and the new components constitute the refined Gaussian component set. The truncation threshold is a preset lower weight limit used to filter out low-weight invalid components; the merging threshold is a preset upper distance limit used to determine whether two Gaussian components belong to the same target.

8. A method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The step of extracting component information from the refined Gaussian component set as a state estimate includes: Traverse the refined Gaussian component set, filter out the components whose weights are greater than or equal to the extraction threshold, and add the corresponding state mean to the estimated target state set. If the weight of the target component in the refined Gaussian component set indicates that there are multiple targets at the position corresponding to the target component, the corresponding number of overlapping points are extracted according to the rounding principle and added to the estimated target state set, and the points in the estimated target state set are used as the state estimate. The extraction threshold is a preset weighted probability value used to determine whether a target exists.

9. A method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 1, characterized in that, The state transition update for free track, experimental track, and official track includes: Obtain the state estimate and the free point track, experimental track and formal track established at historical time points; Calculate the Euclidean distance between the state estimate and the predicted position of the existing track. If the Euclidean distance is less than the associated gate, the association is considered successful. The state estimate that fails to be successfully associated with any existing track is determined to be a free point; The state estimate that is successfully associated with the free point trace in the previous time step in the current time step is upgraded to the experimental track. When the hit counter of the experimental track reaches the initial confirmation frame number, the experimental track will be upgraded to a formal track. The associated gate is a preset distance threshold used to determine whether the observation point track and the flight track belong to the same target.

10. A method for real-time tracking of dense multi-UAV targets based on GM-PHD filtering according to claim 9, characterized in that, The state transition update of the formal track also includes: If the established formal track is successfully associated at the current moment, the state of the formal track is updated using Kalman filtering. If the formal track is not associated with any point, perform a pure prediction update and increment the loss counter; When the loss counter value of the formal track exceeds the maximum loss deletion threshold, the formal track is deleted; The maximum loss deletion threshold is a preset maximum number of frames that a track can remain undetected for an extended period.

Images