CN117075631A - Detection tracking identification method based on unmanned aerial vehicle cluster target - Google Patents

Abstract

The invention discloses a detection tracking identification method based on unmanned aerial vehicle cluster targets, and belongs to the radar cluster target tracking technology. The method comprises the following steps: firstly, predicting the motion state of a target of a unmanned aerial vehicle cluster by using different cluster target detection tracking modes according to different types of targets in the current unmanned aerial vehicle cluster; then, the predicted motion state of the unmanned aerial vehicle cluster target is updated by using tolerance multimode, and the updated result is used as a final tracking and identifying result; the tolerance multi-model is constructed by using a multi-model extended Kalman filter, the predicted motion state is calculated by using a plurality of linear models, and then weighted summation is carried out to obtain the updated motion state of the unmanned plane cluster target. The method can adapt to the change of the formation of the unmanned aerial vehicle group, obtain the estimation of the more real group expansion state, and improve the tracking precision of the group target.

Description

Detection tracking identification method based on unmanned aerial vehicle cluster target

Technical Field

The invention belongs to the technical field of radar swarm target tracking, and particularly relates to a detection tracking identification method based on an unmanned aerial vehicle swarm target.

Background

The distinction between extended target tracking and point target tracking arises from the nature of the radar, particularly the resolution of the radar. If the resolution is high enough relative to the object size, one object can occupy several resolution cells, thus representing an extended object. As shown in fig. 1, each resolution cell occupied by a target is defined as a measurement source of the target, in which case the target generates multiple measurements at each instant, i.e., the measurement source generates measurements with a probability of detection. The measurement is fully utilized, and the position and the appearance of the extended target are tracked at the same time, so that the method is beneficial to tracking different types of targets. In addition to a shaped expanded object, a group object is composed of a plurality of sub-objects having the same motion state, which are not individually tracked but are integrated. Thus, the target as a whole occupies a plurality of metrology resolution cells, with sub-targets occupying one or more of the resolution cells.

Along with the rapid development and application of the large-scale unmanned aerial vehicle formation flight technology, the traditional radar is limited by time resources, processing resources, resolving power and tracking strategies, and the problems of insufficient target tracking quantity, mixed batch, batch mixing, batch loss and the like can occur when dense group targets are tracked; meanwhile, the radar has the problems of missed detection, false alarm and the like, and random change of the number of targets is caused. Therefore, when the conventional radar is used for tracking the group targets, the problem of lower precision exists, and the requirements cannot be met.

Disclosure of Invention

Aiming at the problem that the tracking positioning accuracy is low when the traditional radar method tracks the group targets, the invention provides a novel detection tracking recognition method based on the unmanned aerial vehicle group targets, realizes accurate estimation of the group number, the group centroid motion state, the group shape and the intra-group target state, and improves the group target tracking positioning accuracy.

The invention provides a detection tracking identification method based on an unmanned aerial vehicle cluster target, which comprises the following steps:

firstly, according to different types of targets in the current unmanned aerial vehicle cluster, different cluster target detection tracking modes are used, as described in the steps one to four:

when the targets in the group are indistinguishable, a centroid group tracking algorithm is used for detecting and tracking the targets in the group;

secondly, when the targets in the group are distinguishable and are non-rigid targets, detecting and tracking the group targets by using a group target and multi-target information interaction tracking algorithm;

step three, when the target in the group can be resolved and is a rigid target, if the group configuration is known, using an extended group target tracking method A to detect and track the group target; the extended group target tracking method A obtains the motion centroid of each group of tracks based on an extended Kalman filter;

Fourthly, when the targets in the group can be resolved and are rigid targets, if the configuration of the group is unknown, using an extended group target tracking method B to detect and track the group targets; the extended group target tracking method B uses a Bayesian estimation framework based on a random matrix extended state to detect and track group targets; wherein the random matrix is described by inverse Weisal distribution, gaussian distribution is used for the target motion state, and gamma distribution is used for the measurement rate state;

secondly, constructing a tolerance multi-model, namely, updating the predicted motion state of the unmanned aerial vehicle cluster target by using a tolerance multi-model, and taking the updated result as a final tracking and identifying result;

fifthly, constructing a tolerance multi-model by using a multi-model extended Kalman filter, wherein the tolerance multi-model comprises the following steps:

the detection tracking recognition problem of the unmanned aerial vehicle cluster target is a nonlinear estimation problem, and a posterior probability distribution function of the motion state of the unmanned aerial vehicle cluster target is set to obey Gaussian distribution; setting M particles, wherein each particle uses a linear model to approximate the Gaussian distribution, and the probability accumulated value of the M linear models is approximate to a true value; m is a positive integer;

and calculating the corresponding expansion states of the predicted movement states of the unmanned aerial vehicle cluster targets by using M linear models, and carrying out weighted summation on the M expansion states to obtain updated movement states of the unmanned aerial vehicle cluster targets.

Compared with the prior art, the method has the advantages that:

(1) The method realizes the switching and scheduling of various tracking modes, constructs a novel tolerance multi-model of a multi-model extended Kalman filtering algorithm (MMEKF), establishes a unified group target tracking frame through multi-model estimation, approximates a continuous state space by using the concept of multi-probability weight and designing a limited probability weight point, can obtain an extended point which is closer to a true value, can adapt to the change of an unmanned aerial vehicle group formation, and obtains more real group extended state estimation. The method realizes joint estimation and decision, has the capability of obtaining the global optimal solution, and improves the group target tracking precision.

(2) The method of the invention uses an extended group target tracking method to identify and track targets of unmanned aerial vehicle clusters with known group formations, improves the processing of point cloud information, performs gating, scoring and distributing processing on the point cloud information, and completes the state extraction of the point cloud information based on an extended Kalman filtering algorithm. Aiming at the group target tracking with unknown group formation, the method provides a random matrix method for evaluating the expansion state by inverse Weisal distribution description, and more accurately describes the complex expansion state of the group target so as to improve the group target tracking precision.

Drawings

FIG. 1 is a schematic diagram of tracking an extended target (left) and a group target (right) using radar;

FIG. 2 is a diagram of an overall implementation framework of the detection tracking identification method of the present invention;

FIG. 3 is a flow chart of a centroid group target tracking algorithm in the method of the present invention;

FIG. 4 is a flow chart of a group information and multi-target interactive tracking algorithm in the method of the present invention;

FIG. 5 is a schematic diagram of classification of a drone swarm configuration;

FIG. 6 is a schematic diagram of an extended target tracking algorithm for the known case of group configuration in the method of the present invention;

FIG. 7 is a flow chart of an extended target tracking algorithm with known group configuration in the method of the present invention;

FIG. 8 is a schematic diagram of an ellipse represented by a random matrix;

FIG. 9 is a diagram of a ladder quantization method of the present invention with a given tolerance;

fig. 10 is a flow chart of the MMEKF filtering of the present invention.

Detailed Description

The invention will be described in further detail with reference to the drawings and examples.

Firstly, the invention analyzes the research target and determines the group tracking mode adopted for different situations. The group targets may be classified into group targets that are distinguishable from the group targets and group targets that are not distinguishable from the group targets from the perspective of radar detection. Group targets with fixed group configuration and group targets with unfixed group configuration can be classified from the characteristics of the targets themselves. The method of centroid group target tracking can be adopted for the indistinguishable targets, and the distinguishable targets can be divided into rigid targets (targets with rigid connection inside) and non-rigid targets such as bird groups, machine groups, shrapnel groups and the like. The non-rigid body targets can adopt a group information and multi-target interaction algorithm, and each target in the group can be tracked while the group information is reserved. For rigid objects, it essentially corresponds to a group configuration fixed object. In this case, two types of targets are classified, one is a rigid body target with a known configuration, and the other is a target with an unknown configuration. The target with known configuration can adopt a general extended target tracking algorithm, and the target with unknown configuration needs to be tracked after the extended state is estimated through a random matrix. And finally, dispatching the related methods by using a tolerance multi-model method to form a unified group target tracking framework. The overall implementation framework of the unmanned aerial vehicle cluster target-based detection tracking identification method is shown in fig. 2, and is described in five steps below.

Step one, detecting and tracking the group targets by using a centroid tracking mode under the condition that the targets in the group are indistinguishable.

In many realistic scenarios, the target exhibits characteristics of cluster motion. Such as migratory bird candidates, traveling fleet, sailing fleet, etc. When the target is relatively close to the sensor, the measurement data is more accurate, so that the target can be tracked as multiple targets. However, when the target is far from the sensor, the sensor cannot distinguish a plurality of targets due to insufficient resolution of the sensor radar. In this case we are more concerned with the motion state of the whole object group, such as the centroid position, velocity, acceleration, etc. of the object group, i.e. centroid tracking of the object group. As the targets move from far to near, more targets are observed by the radar, and there is a possibility that targets originally observed as one group actually become a plurality of groups, so that group separation is required. When the targets show cluster movement and the number of targets is many, the targets can be divided according to clusters and calculated in parallel because the calculated amount of data association and the number of targets are in an exponential relation, so that the calculated amount is reduced, and the cluster target tracking under multi-target tracking is realized.

And 1.1, establishing a group target model.

When tracking the group target, the main judgment basis is the group track center and the wave gate selection. And (3) carrying out data association on the measurement and the group track center by selecting the measurement falling into the group wave gate, updating the group centroid track by using the successfully-associated measurement, and estimating the related information such as the group speed.

Let the sensor track t groups in clutter environment, the state equation of the system is:

X ^t (k+1)＝F(k)X ^t (k)+G(k)u ^t (k)+V(k),k＝1,2,…,T

wherein X is ^t (k+1) is a global state vector of the group t at the time of k+1, including the position and the speed of the center of the group t, in a two-dimensional coordinate systemIn three-dimensional coordinate system->Acceleration dimensions, i.e. global state vector +.>(x, y) and (x, y, z) respectively represent two-dimensional position coordinates and three-dimensional position coordinates of the group t at the time k+1; f (k) represents a state transition matrixG (k) is an input control matrix acting on state variables, u ^t (k) V (k) represents zero-mean Gaussian white noise for the input signal at time k. Here T represents the time of observation T. The superscript' indicates the transpose.

The sensor measurement equation is:

z ^t (k)＝H(k)X ^t (k)+W(k)

wherein H (k) represents the measurement matrix, W (k) represents zero-mean Gaussian white noise, z ^t (k) The set of measurements representing group t at time k requires the use of a group gate to determine which of the total measurements received by the sensor are associated with group t, the criteria for determining are:

Wherein,representing the normalized distance between the ith measurement and the center of the wave gate and the center of the cluster, +.>Representing the size of the wave gate, and the maximum value of the normalized distance; />Representing information between one-step prediction and measurement, wherein the one-step prediction refers to a predicted value obtained according to a sensor measurement equation, and the information represents a difference value between the measured predicted value and a measured actual value;is a innovation covariance; the upper corner mark-1 indicates the inverse matrix. />The calculation is as follows:

wherein,related to the distribution of the measurements, representing an estimate of the distribution of the measurements; />Representing the measurement error covariance of group G; h' (k) represents a transpose of the metrology matrix; p (P) ^t (k|k-1) represents a covariance matrix of the state variables; group G represents a group after estimating the measurement distribution of group t.

Since the sensor usually uses a spherical coordinate system, i.e. azimuth, pitch and radial distances, to observe the target, and the system state space coordinate system is usually a rectangular coordinate system, it is necessary to convert the measurement coordinate system and the system state space coordinate system as follows:

wherein C represents a transformation matrix between the measurement coordinate system and the system state space coordinate system;the measurement covariance of group t in the measurement coordinate system is represented.

After the relevant measurements are selected, the center of the cluster can be obtained by the measurements as follows:

Wherein L is _t A state vector representing the center of group t; n (N) _max Representing the number of measurements falling into the wave gate; z is Z _t,i Representing the i-th measurement vector within group t.

The above-described group judgment mechanism is established by distinguishing the measurements using the wave gate, and then calculating the group center using the measurements in turn.

And 1.2, tracking the group target by using a centroid group tracking algorithm.

As shown in fig. 3, the centroid group tracking algorithm (Centroid Group Tracking) first determines whether the measurement belongs to the group, i.e., determines whether the measurement falls within the gates of the group, and all the measurements falling into the gates can be used to derive the measurement center. And (3) associating the predicted states of the measurement center and the mass center of the group, and updating the mass center state by using the measurement center to obtain the information such as the position, the speed and the like of the mass center of the group if the association is successful.

Set existing group target track X ^t (k) A one-step prediction of the state of the group object can be obtained by kalman filtering. It is calculated whether the measurements received by the sensor at time k+1 fall within the gates of the group. The present embodiment employs an elliptic wave gate. According to different shapes of group targets, different wave gate judgment can be selected. The shape of the group object is generally elliptical, rectangular, star-shaped, or any other shape.

The initial time is tracked without a determined group target track, so that measurement can be directly clustered, a measurement center is built through measurement of the same group, and the initial state of the group is built through kalman filtering. The speed is obtained from at least 2 moments of measurement information. Acceleration requires at least 3 moments of measurement information to obtain.

After establishing the gates of the group trajectories, the measurements meeting the requirements are selected for each group. However, the determination of the measurement is not simply a wave gate determination. Often there are cases where false measurements fall into the gate or where measurements from real targets do not fall into the gate due to cluster maneuvers or predictive inaccuracies. Thus, further discrimination of the measurement is required.

Measurement set of gates falling into group GZ _G,i Representing the i-th measurement vector within group G. Selecting the measurement of which the normalized distance is the smallest +.>For seed measurement, seedThe sub-measurement is taken as the center, and the radius is built againIs a wave gate of (c). And judging all the measurements in the measurement set again, reserving the measurements in the newly built wave gate, and deleting the measurements not in the newly built wave gate. After the determination, part of the false measurement is measured from the measurement set Z _G Kicking out the middle part. Again from the measurement set Z _G The seed measurement is selected and taken as the center, and the rest measurement is judged again. And re-judging whether the measurement which is not divided into groups is positioned in the wave gate of the existing target group and dividing the measurement, if the existing measurement does not belong to any existing group, establishing a new group until all the measurements are divided into the groups. Thus, measurements divided into clusters satisfy at least two conditions, the gate of one-step prediction of cluster centroid trajectories and the gate of seed measurements. The method can ensure the accuracy of the division measurement to the maximum extent.

Data association may begin after the measurements are divided into clusters. The subject of data correlation is one-step prediction of cluster centroid trajectories and measurement sets within clusters. The cluster centroid trace is updated with the associated metrology set by the kalman filter. When the measurement coordinate system is different from the system state coordinate system, the track can be updated by nonlinear Extended Kalman Filtering (EKF) or Conversion Measurement Kalman Filtering (CMKF).

Since multiple measurements are used to update a group track, the noise covariance consists of all the measurements

Wherein N is _i For the number of measurements in group G,for predicting the number of targets within group t, < > j- >To weightAnd the weighting parameter is 1 when only one measurement exists in the group. />For the distribution estimation parameter of the intra-group measurement, the subscript D represents the distribution parameter of the group G track after updating.

Prediction of the number of targets within a groupDistribution estimation parameter for intra-group measurement>Is that

Where, here, alpha denotes the filter gain constant,from the group distribution. DG (differential g) ₀ Representing initialized distribution parameters of group structure, N ₀ Representing the number of initialisation targets within a group, +.>Representing the distribution estimation parameters of the initial intra-cluster measurements.

And secondly, detecting and tracking the group target by using a group target and multi-target information interaction tracking algorithm under the condition that the target in the group can be resolved and is a non-rigid target.

The multi-target tracking technology is the basis of group target tracking, but the group target tracking has its own characteristics. The biggest difference between the two is that the target states in group target tracking cannot be regarded as independent of each other, so that the conventional multi-target system model cannot well describe the movement of the group targets. The existing mainstream group target tracking algorithm mainly comprises two types: the first type uses group structure update as a framework, and uses group structure information to assist multi-target tracking information, wherein the group structure information is updated preferentially and independently. And secondly, taking multi-target tracking as a framework, and extracting group target information from the multi-target information after each multi-target tracking update is finished. Thus, the multi-target information and the group structure information can be mutually assisted and mutually influenced. Therefore, the invention uses a tracking algorithm of closely interacting group information and multi-target information based on the logic of group target tracking and multi-target tracking to improve the group target tracking precision.

As shown in FIG. 4, the interactive tracking algorithm for group target and multi-target information adopted by the invention mainly comprises 9 steps, which are denoted by numbers and mainly divided into three parts of measurement grouping, group structure updating and group structure correction. Wherein the multi-target portion uses a track-oriented MHT (multi-hypothesis tracking) algorithm. Steps (1), (2), (3) and (9) are group structure updating steps, and steps (4) to (8) are multi-objective updating steps. Meanwhile, the multi-target information is used for assisting in updating the group structure, and the group information is used for assisting in updating the multi-target information.

Step (1) obtaining a group structure at the previous moment, namely the k-1 moment; performing the measurement grouping of the step (2), wherein the measurement grouping can use the mode of the step one; updating the group structure in the step (3); step (4) obtaining multi-target state information at the moment k-1; step (5) multi-objective state prediction is carried out; step (6) carrying out the association of the track and the measurement; step (7) generating a hypothesis; step (8) updating the multi-target track; and (9) performing group structure correction.

Step (7) generating hypotheses refers to generating hypotheses for different possible types of tracks, and then screening. Specifically, according to different multi-target tracking algorithms, a K-best algorithm is used for an assumed MHT algorithm, a sequential probability ratio test algorithm is used for screening local track scores for the MHT algorithm for tracks, tracks with extremely low probability are deleted, and then a track hypothesis is generated through global track probability.

The step (9) of group structure correction refers to updating the group structure by volume measurement. For example, when two groups of structures cross, the measurement packet uses the normalized distance as the criterion of the wave gate, so that the measurement of two groups that are closer to each other will be divided into one group, resulting in the two crossed groups being merged into one group or one group being terminated. The subsequent measurement continuously causes the state of the group structure to deviate seriously, the group structure is changed into a structure with 2 groups combined, the structure is not in accordance with the actual situation, and a large error is caused to the subsequent target tracking. In practical situations, the group target stably tracked only generates larger error at the initial moment, and the group target is stably tracked in the subsequent process, so that the group structure can be updated by using the stably tracked multi-target state in the group structure.

In each group, the multi-objective state is predicted, and the group structure state affects the multi-objective state prediction. Meanwhile, due to the fact that the distance between the group targets is small, the relative speed is small, shielding is easy to occur, and the track is enabled to terminate and start for many times. Therefore, for tracks without measurement update, the cluster structure state is used as the track state to continuously update N times, namely, the cluster structure is used for assisting the updating of the weak and small targets. If after N times the track is still not associated with any measurements, the track is deemed to have terminated. The method can effectively solve the problem of shielding among targets and interruption of measurement of weak and small targets.

And thirdly, if the group configuration is known, performing group target detection tracking by using an extended group target tracking method under the condition that the targets in the group can be distinguished and are rigid targets.

The unmanned aerial vehicle group configurations may be categorized from different angles. The configuration of the unmanned aerial vehicle group is classified based on the coupling between unmanned aerial vehicles, as shown in fig. 5, such as a physical coupling type, b formation type, c group type, d intentional cooperation type, etc. The physical coupling type means that unmanned aerial vehicles are connected through physical links, and the movement of the unmanned aerial vehicles is constrained by the force of the movement of other unmanned aerial vehicles. The formation type, the unmanned aerial vehicles are not physically coupled, but rather their relative movement is severely limited to maintain formation. The exercise planning problem is then considered and formulated as a whole. To avoid problems with collisions in the team, it is necessary to embed them into the formation control strategy. Group types, which are homogeneous teams composed of many unmanned aerial vehicles, create emerging collective behaviors. The final movement of the drone group does not necessarily lead to formation. Scalability is a major problem because of the large number of unmanned aerial vehicles involved, and then the imperative of a purely decentralized control architecture. The intentional collaboration type, the movement track of the unmanned aerial vehicle of the team is defined by each executing global task. Therefore, these unmanned trajectories are typically geometrically uncorrelated. In this case, problems such as multi-drone task allocation, advanced planning, resolution planning, and conflict resolution should be solved, and global tasks to be performed and the different drones involved should be considered. In this case, both centralized and decentralized decision architectures can be applied.

For unmanned aerial vehicle formation such as physical coupling, the configuration is known a priori, so that the unmanned aerial vehicle formation can be processed by adopting an extended target tracking method, as shown in fig. 6, and the specific flow is shown in fig. 7. Wherein the function in the predicting and updating steps is the operation of a classical Extended Kalman Filter (EKF), and the operation in the associating and assigning steps is used to support the processing of multi-point data.

The processing steps of the extended group target tracking method of the invention are as follows: 1) Predicting, namely calculating an estimated value of the centroid state of the target track of the current time group through the motion centroid of the target track of the previous time group by an EKF (extended Kalman Filter) predicting step; 2) And (3) gating and scoring, wherein all sensor measurement values at the current moment are input in the association process, and the measurement falling into the current track wave gate is screened out through a gating function and scored. Because one measurement can fall into a plurality of track wave gates at the same time, each measurement only keeps the optimal score finally, and the measurement is associated with the track with the optimal score; 3) The distribution, for the measurement which is not associated with any track, the distribution process executes clustering treatment on the measurement, and a new track is distributed for the measurement when the track generation condition is met; 4) Updating, namely updating each track, and calculating a group innovation amount and a group innovation covariance matrix through the mass center estimated amount of each track and the measurement related to the mass center estimated amount to obtain the motion mass center of the group at the current moment.

Three states exist for each sensor measurement point: free state, detected state and active state. There are two possibilities for tracking the track per frame: if there are measurement points associated with the track, then hit, and if there are no measurement points associated with the track, then track is lost. The state of the measurement point is converted by continuously hitting or losing frames as a threshold. The detection state and the activation state are unique states of the track centroid and can be regarded as the states of the track, and the target tracking algorithm only predicts and updates the measuring points in the non-free state. Only the track in the active state is considered to be the detected target.

At the beginning, all measuring points are in a free state, and at the moment, all measuring points generate tracks through an allocation process and convert the tracks into a detection state. The initial process is to use the measuring point with the number of 1 as a centroid, and use the centroid to perform speed detection first and then perform distance detection. And comparing the absolute value of the speed difference between the two points with a threshold value by setting the speed threshold value, and detecting if the absolute value is smaller than the threshold value. And similarly, setting a distance threshold value, comparing the square of the Euclidean distance of the two measuring points with the distance threshold value, and detecting if the square is smaller than the distance threshold value. And (3) re-calculating the mass center (detecting the average value of qualified points) of the point set detected by the method, adding the measuring points into the point cloud, detecting the next measuring point if the measured point does not pass the point cloud, and repeating the process until all the measuring points are detected. Once completed, detection of points and signal to noise ratios is performed for the point cloud. And a set point number threshold value, wherein if the number of the measuring points in the point cloud is larger than the threshold value, the set point number threshold value is regarded as passing detection. And setting a signal-to-noise ratio threshold, and summing all measured signal-to-noise ratios in the point cloud, wherein if the sum is larger than the threshold, the sum is regarded as passing detection. After the process is finished, the algorithm forms a new target track and assigns numbers to the tracks. If the initial centroid does not pass, the initial centroid is replaced and the process is repeated.

If the number of frames of the new track hit continuously meets the threshold requirement of 'detect activation', namely the duration of the new track meets the threshold requirement, the track (centroid) is converted into an activated state, otherwise, the waiting is continued.

For the motion centroid in the active state, the following three conditions may occur in the updating process, so that tracking is lost:

(1) If the target is very slowly approaching stationary due to stationary or moving objects, the assumption is treated as stationary clutter filtering.

(2) If the target is not stationary but is not detected, the target is considered to have left the observation area.

(3) If this is the case without detecting the target, it is presumed that the target is occluded by another object.

The processing method corresponding to the above situation comprises the following steps: if the x, y, z axis velocities are simultaneously less than the velocity resolution in the state estimator, then the object is considered stationary and a large "stationary to free" threshold (longer time threshold) is used to extend its observation time for objects that are within the boundary and stationary. The state quantity of the previous moment is reserved at the next moment of prediction. If the boundary is exceeded, a "leave to free" threshold (a shorter time threshold) is used to quickly remove the out-of-bounds target. If this is the other case, the "activate to free" threshold is used to continue the motion of the object according to the object's motion model. For the first and third cases, if the number of frames lost continuously exceeds the corresponding threshold, the track is deleted.

The method of the invention requires the prior information of the group target formation to be utilized to calculate the gating function. The gating function represents the amount of innovation present in the current track that is acceptable for the current uncertainty state, in order to be able to measure this uncertainty, a group innovation covariance matrix C is defined _G The method comprises the following steps:

wherein S is _apr (n) state vector, J of unmanned aerial vehicle cluster G at time n _H Representing the Jacobian matrix, P _apr (n) represents an a priori estimate matrix of the state vector covariance matrix at time n.Representing centroid uncertainty due to target maneuvers, R _G Is a measurement error covariance matrix, C _D Is the track variance estimation matrix. The upper corner T here represents the transpose of the matrix, the lowerThe corner mark D represents the estimated error of clustering to form the track centroid. G in FIG. 6 _apr (n) represents the amount of innovation acceptable to the n-time gating function.

For each group track, a distance function may be defined for all measurement vectors obtained at the time of track updateThe following are provided:

wherein j represents the obtained j observation point, i represents the stored i-th track, C _i Representing the variance of track i, y _ij An innovation vector representing an observation point j on a track i, inv represents inverting a matrix, u _j (n) represents the j-th measurement vector at time n, each measurement vector u (n) contains the distance r and azimuth angle And pitch angle θ information, C _Gi Representing covariance of track i, +.>Representing the observation matrix of the state vector.

The distance function represents the amount of innovation that the new measurement adds to the existing track. Then limiting the accepted information amount by chi-square detection, and meeting the requirementd _max Indicating an acceptable maximum distance, may be achieved by pre-processingThe set gating parameters, including the length, width, height and speed limits, are calculated and these parameters need to be obtained from a priori information on the formation of the unmanned aerial vehicle cluster.

And step four, an extended group target tracking method with unknown group configuration.

The extended target tracking is used as a technology widely applied to modern high-precision sensors, and provides more accurate and multidimensional target state estimation, such as extended state, measurement rate state, class state estimation and the like, than traditional point target tracking, and provides rich target information for post-processing of tracking system classification, identification, situation estimation and the like. In particular, the object expansion state estimation is a difficulty in tracking an expansion object due to the arbitrary shape and complexity of the object. The two most representative extended state estimation methods are the random hypersurface and random matrix methods. The random hypersurface method cannot accurately describe complex expansion states due to the limitation of function description capability. In addition, the star-convex random hypersurface method does not form a unified extended state estimation evaluation criterion, and is also an important factor for restricting the development of the random hypersurface method. As a mainstream extended state estimation research method, a random matrix method provides a convenient extended state estimation evaluation method while estimating a target extended state through strict mathematical modeling and probability statistical description, and has become a hotspot of extended state estimation research. The embodiment adopts a Bayesian estimation framework based on a random matrix extension state to realize the tracking of the resolvable group targets with unknown configuration.

A random matrix is defined as a matrix form of at least one random variable element. In the 30 s of the 20 th century, wishart proposed a random matrix concept and conducted intensive research on multidimensional random matrices, including defining random variable element joint distribution, matrix eigenvalue distribution, overall random matrix distribution form, and the like. In 1967, wigner first described physical phenomena such as sensor measurement variance with a random matrix. In 2008, koch first introduced a 2-dimensional random matrix form to estimate the target expansion state. The method models the expansion state of the target as an ellipse, and uses a 2-dimensional positive definite random matrix to represent the size and direction of the ellipse, as shown in the following formula:

wherein the first random matrix represents an ellipse having no direction, i.e., a direction parallel to the axis of the 2-dimensional cartesian coordinate system, and a and b represent the major and minor half axes of the ellipse, respectively, as shown in fig. 8 (a); the second random matrix is a general form of the first matrix, that is, the first matrix rotates in an arbitrary direction θ, the matrix direction coincides with the matrix eigenvector direction, and the long and short axes are obtained by matrix eigenvalue evolution, as shown in fig. 8 (b). The upper corner mark T here indicates a transpose operation.

In the method, a random matrix is described by inverse Weisaute distribution in a Bayesian estimation framework based on a random matrix expansion state. In addition, the target motion state and the metrology rate state are described by gaussian and gamma distributions, respectively. For clarity, the probability density function involved in the random matrix extended target tracking method is described below.

The gaussian probability density function is as follows:

wherein X here represents the variable of the function, a vector of d dimensions; m and P represent the mean and variance of the gaussian function, respectively. The upper corner mark T here indicates a transpose operation.

The inverse weisat probability density function is as follows:

wherein V and V represent the degree of freedom and matrix parameters of the inverse weisal distribution, respectively; where X represents the variation of the functionQuantity is a vector of d dimensions; Γ -shaped structure _d Representing a gamma function; etr represents an exponential function.

The weisauter probability density function is as follows:

wherein v is _w And V _w The degrees of freedom and matrix parameters of the Wittig distribution are respectively represented; in the Weisalt probability density function, n _N Represents n _N Normal samples of independent same distribution, n _W Represents n _W A normal sample is independently and uniformly distributed. Where X represents the variable of the function and is a vector of d dimensions.

The gamma probability density function is as follows:

Here, α and β represent a shape parameter and a scale parameter of the gamma distribution, respectively, and γ is a random variable.

The poisson probability density function is as follows:

where W and λ represent the random variable and the mean, respectively, of the Poisson distribution, the mean also being the variance.

(a) And carrying out state prediction on the ellipse expansion target.

When the extended target is measured randomly and distributed in a certain space, an ellipse can be generally used to approximately describe the extended state of the extended target, namely, the extended state of the extended target is tracked by the ellipse, wherein the range of the ellipse is marked to represent the space position occupied by the extended target. The ellipse-extending target state can be modeled as:

ξ＝(γ,x,X)

wherein, gamma, X and X respectively represent the measuring rate state and the motion of the expansion targetThe state and the extended state are both contained in the augmentation term ζ; gamma is a scalar subject to gamma distribution; vector x can be modeled as x= [ x, y, v, θ, ω] ^T ，[x,y] ^T Representing the two-dimensional coordinate position in the motion state of the target, v, θ, ω representing the speed, direction and steering rate of the target, respectively. In the following description, x is used to represent x, representing a motion state vector of an elliptical expansion target; using X for X, the extended state vector is represented.

The general elliptic expansion target setting obeys Gaussian distribution, and the probability density function of the elliptic expansion target state at the time k-1 obeys gamma Gao Sini Weisal (GGIW) distribution as follows:

Wherein, gamma _k-1 ,x _k-1 ,X _k-1 Respectively represent the measurement rate state, the motion state and the expansion state of the expansion target at the time k-1, Z ^k-1 Representing a measurement set from the start of tracking to the time k-1; the measurement rate state is subject to gamma distribution, alpha _k-1 、β _k-1 Gamma distribution parameters corresponding to the measurement rate state; the motion state of the target obeys Gaussian distribution, m _k-1 、P _k-1 Is a Gaussian distribution parameter corresponding to the motion state; the extended state of the target obeys the inverse Weisal distribution, v _k-1 、V _k-1 Is the inverse Weisal distribution parameter corresponding to the extended state. Because each target state is subject to a GGIW distribution, the described elliptical extended target tracking method is also referred to as the GGIW method.

Before giving the target prediction state, some necessary assumptions are given.

A.1: the measurement rate prediction state of the expansion target is irrelevant to the motion state and the expansion prediction state, and is subjected to a heuristic prediction form.

A.2: the motion prediction state of the expansion target is independent of the expansion state.

A.3: the extended predicted state of the extended target depends on the motion state.

The target predicted state probability density function at time k can be expressed as:

wherein f (gamma) _k |γ _k-1 )，f(x _k |x _k-1 ) And f (X) _k |X _k-1 ,x _k-1 ) Respectively representing a measurement rate state, a motion state and an expansion state prediction function; based on the above assumptions, there are the following probability functions:

Wherein F is _k-1 And Q _k-1 Transfer matrix and process noise matrix, delta, respectively representing target motion state _k Represents the degree of freedom of the Wishade distribution, M (x) _k-1 ) A directional rotation matrix function representing the steering rate, namely:

where q is the sampling interval, ω _k-1 The steering rate at time k-1 is shown.

The target predicted state probability density function at time k still has the form GGIW, i.e

/>

Wherein, the forgetting factor is 1/tau, and 1/tau is less than 1. Alpha _k|k-1 、β _k|k-1 Is the prediction parameter of the measurement rate state at the moment k; m is m _k|k-1 ,P _k|k-1 Is a prediction parameter of the motion state at the moment k; v _k|k-1 、V _k|k-1 Is the prediction parameter of the extended state at time k. In the method of the invention, the motion state and the extension state predict the parameter m due to the nonlinear modeling of the motion state _k|k-1 ,P _k|k-1 ,v _k|k-1 ,V _k|k-1 Calculation was performed using the extended kalman filter technique.

Measurement set w= { z composed of target multiple measurements received by sensor at k time _k,1 ,z _k,2 ,...,z _k,|W| W represents the number of measurements in the set W, z _k,j Representing the j-th measurement, the measurement likelihood function is:

wherein H is _k Representing the measurement matrix, R _k Represents measurement noise, eta represents a regulating parameter, and the larger eta indicates R _k The less the impact on the target state estimate,and->Representing the mean and variance vectors of the set W, respectively. Z is Z _k Representing predicted observation vector parameters.

The probability density function of the target posterior state at time k can be obtained:

wherein Z is ^k ＝Z _k ∪Z ^k-1 The measurement rate state is updated as:

α _k ＝α _k|k-1 +|W|,β _k ＝β _k|k-1 +1

the motion state is updated as follows:

m _k ＝m _k∣k-1 +K _k G _k

P _k ＝P _k∣k-1 -K _k H _k P _k∣k-1

wherein K is _k And G _k Is an intermediate parameter; here the superscript T denotes a transpose;representing the expected value of the extended state.

The extended state updates are:

v _k ＝v _k∣k-1 +|W|

wherein N is _k As an intermediate parameter, the upper corner mark T represents transposition; b (B) _k Representing covariance of the noise vector; where d is the motion state x _k Is defined in the vector dimension of (a).

The motion state of the extended target can be calculated by a probability density function of the posterior state of the elliptical extended target at the moment k.

(b) And carrying out state prediction on the non-elliptical expansion target.

When the spatial distribution of the plurality of measurements of the expansion target can approximately reflect the target shape, if only one ellipse is used to describe the expansion state, the information loss is large, and the target information contained in the expansion target measurement cannot be fully mined. At this time, a plurality of ellipses are used for describing the complex shape, so that more target shape details can be effectively captured. The present invention also uses the gamma Gao Sini Weisal distribution for the probability density function of non-elliptical extended target states. Here, emphasis is placed on a non-elliptical expansion target tracking framework based on a random matrix, and specific state prediction and updating formulas of each ellipse can refer to the description of state prediction of an elliptical expansion target.

The non-elliptical expansion target state may be modeled as:

wherein N is the number of used sub-ellipses, and the sub-ellipses are mutually independent; gamma ray ⁽ⁱ⁾ 、x ⁽ⁱ⁾ And X ⁽ⁱ⁾ The measurement rate state, the motion state and the expansion state of the ith sub-ellipse expansion target are respectively represented.

Before introducing non-elliptical expansion target filtering, some basic assumptions are explained.

A.4: the number of child ellipses N used is known and remains unchanged during the target tracking process.

A.5: the relationship between the sub-ellipses and the measurements is unknown, i.e., no information indicates that the measurements are explicitly associated with any sub-ellipses.

Based on the full probability theory and the Bayesian estimation framework, the posterior state probability density function of the ith sub-ellipse expansion target at the k moment is:

wherein,representing the total number of all possible associated events, if there are N sub-ellipses and N _k The measurement is commonA number of possible associated events; />Representing a first association event between the measurement and the child ellipse; />Representing associated eventsThe update state of the ith sub-ellipse; />Representing associated event->Probability of (2), i.e.)> Representing associated event->Is associated with the set of measurement components of the ith sub-ellipse. />Representing associationsEvent->The predicted state of the ith sub-ellipse. Given associated event +. > I.e. < ->Representing associated eventsThe ith sub-ellipse of (I) and its corresponding measurement set +.>Likelihood probability density functions between. The factor ellipse prediction state is independent of the associated event, available +.>

Based on Bayesian estimation theory, probabilityThe unfolding method comprises the following steps:

wherein,normalization constant for Bayesian estimation,>if the correlation event between the measurement and the sub-ellipse does not have any prior information, +.>Can be regarded as an average of the total number of events, i.eThe factors ellipses are mutually independent, and the likelihood function is measured>Can be obtained by multiplying likelihood functions between different sub-ellipses by using the edge integral theorem,/>The unfolding method comprises the following steps:

and calculating the motion state of each expansion target by using a probability density function of the posterior state of each sub-ellipse expansion target at the k moment.

And fifthly, a group tracking method based on a tolerance multi-model.

The invention aims to solve the problem of detection, tracking and identification of unmanned aerial vehicle cluster targets, which is a nonlinear estimation problem.

EKF is widely used to solve the problem of non-linear estimation. By linearizing a nonlinear system function, EKF is typically studied by gaussian process driven systems based on first order taylor expansion. However, even with optimal linearization, a non-optimal state estimate will be obtained since in the non-linear case the states and measurements are not co-gaussian distributions. Therefore, the EKF can obtain an approximate estimate only when the state deviation is small. As practical systems become more complex, EKFs are limited in many applications due to their reduced performance and even failure. The first-order taylor expansion can well approximate the original nonlinear function in a smaller expansion point range. That is, if a true value is found as an extension point representing the entire state space, it is not necessary to replace the true value with its best estimated value. This can greatly reduce the error between the estimated value and the true value. And even if the best estimator is used, the estimation error cannot be smaller than its Cramer-Rao lower bound. Moreover, since it is the number to be estimated accurately among all estimation problems, a true value cannot be obtained. Since the state space is continuous, it is computationally infeasible to enumerate all points probabilistically throughout the state space. Finding a better expansion point based on the first-order taylor expansion within the EKF framework, making it closer to the true value, and obtaining a more competitive nonlinear filter has become the main research direction. In EKF, the extension point determined by the best estimate calculated by the filter is considered to be the closest ideal extension point in the state space.

The multi-model extended kalman filter (MMEKF) can significantly reduce its large bias. The expansion points designed in the MMEKF algorithm obey gaussian distribution in the probability model process, and these probabilities represent approximately the entire state space by using a plurality of probability weighting methods. Compared to filters such as EKF, UKF and CKF, MMEKF shows a higher estimation accuracy for unreliable measurements.

From the point of view of decision making procedure, it is difficult to sayAnd->Is an ideal expansion point. That is, a hard decision is one point in the set of all probability possible decision choices. However, although the probability of an event occurring is small, this does not mean that it does not occur. Thus, soft decisions are decision processes that consider more choices in the case of low probability events. Compared with a direct two-stage strategy of single-model estimation after a hard decision process, the soft decision based on multi-model estimation realizes joint estimation and decision, and has the potential of obtaining a global optimal solution. When the probability mass function (pm) of each likelihood of selectionf) When counted, it has the main advantage of helping to reduce decision errors. Similarly, it can be used to find extension points and reduce errors.

(a) First, the design of the probability model in the MMEKF will be described.

Set E x _k |z ^k ]And cov (x) _k |z ^k ) Respectively x _k Posterior distribution p (x) _k |z ^k ) Mean and variance of x _k Representing the motion state of the group object at time k.

In a nonlinear system, given a distribution of random models, it can be assumed that the posterior probability distribution function (pdf) follows a gaussian distribution:

wherein,P _k|k-1 representing the predicted state vector and covariance at time k, respectively,/->P _k-1|k-1 The state vector and covariance, respectively, represent the posterior probability at time k-1.

In short, the probability model design problem is to approximate a continuous distribution with a discrete distribution, in other words, to find a discrete pmf to approximate this continuous gaussian distribution. Let p (x) _k-1 |z ^k-1 ) And p (x) _k |z ^k ) The approximation can be made with discrete pmf as follows:

in the same way as described above,

wherein, M points are selected as expansion points, M linear models are obtained, each model corresponds to a linear state model and a linear measurement model,representing the expansion point state of the ith linear model at time k-1,/for>And the expansion point state of the jth linear model at the k moment is represented. />Representing the probability of the ith linear model at time k-1,/for example>Representing the probability of the jth linear model at time k.

The linear state model of the i-th expansion point at the k moment is as follows:

/>

the linear measurement model of the j expansion point at the k moment is as follows:

wherein,representing a state transfer function, x, in a nonlinear system _k-1 And x _k State vectors respectively representing the k-1 time and the k time; w (w) _k-1 Gaussian white noise representing the linear state model at time k-1; v _k Gaussian white noise of a linear measurement model at the moment k is represented; />Representing a measurement matrix; intermediate parameters

It can be clearly seen that the original nonlinear filtering problem translates into a multi-model estimation problem. The method of the invention uses a multi-model estimation mode to solve the nonlinear estimation problem of the invention. Since the probability model is related to time k, the multi-model estimation problem is essentially a variable structure multi-model filtering problem.And->Can be obtained by the following method.

In order to minimize the distribution mismatch, the present invention sets the cumulative distribution function (cdf) F of the true continuous model c _c (x) As is known, x is a function variable. Let the cumulative distribution function of the discrete random variable be F _d (x) Two cumulative distribution functions F are defined by the following distance functions _c (x) And F _d (x) Distance d (F) _c ,F _d ) The following are provided:

in the above formula, R represents a real number set, and sup represents the maximum value of the function under the value range of x.

Distance function d (F) _c ,F _d ) The definition conforms to mathematical logic and satisfies positive qualitation and symmetry, there is the following quotation 1.

Lemma 1 (minimum model number): under the premise of tolerance epsilon of a given distance measure, F is solved by using constraint conditions _d (x)：

The minimum number of models m= [1/2 epsilon ] = minimum integer not less than 1/2 epsilon.

The appropriate tolerance epsilon is not always easy to produce. In some cases, the number M of models is predetermined directly from the available processing resources. Since the discrete distribution is characterized by the number of particles, the location of the particles and the probability at each particle can be determined by M.

Wherein m is _i Representing the ith particle, arg represents the satisfied by the calculation]X value of the expression in (a).

Thus, all models have the same load, as they all have the same probability of being effective, and each model covers an area with the same probability. It uses a minimum number of models and it is also completely consistent with the convention of assigning an equal initial probability to each model. As shown in fig. 9, the design of the model set is intuitively elucidated. M is obtained by _i May be located in a region of low probability density. For scalar continuous nonlinear functions, the corresponding pmf can be found using the criteria described above. At present due to no m _i And thus has no analytic form ofDirectly obtaining m by a method _i Is a value of (2). The m is indirectly obtained through a quantile function _i Is a value of (2).

The standard normal distribution has the following cumulative distribution function:/>

as with the process shown in fig. 9, if X is a univariate random variable and its distribution function is F, p e (0, 1), the p-minutes of X is obtained as follows:

i.e. the p quantile is the number x in the following formula:

P{X＜x}＝F(x)＝p

the statistical function "norm" in matlab was used to obtain m by the following two steps _i Values of (2)

Step A: a scalar pmf is generated that approximates the standard normal distribution described above using M mass points.

And (B) step (B): order theThen, norm v (p) _i ) Correspondingly get m _i ,i＝1,2…M。

The invention can also search the standard normal distribution minute table in advance by using the known particle number.

It should be noted that mass point m _i I=1, 2..m is actually the expansion point of the taylor expansion. ObtainingAndthe method of (2) is the same as in fig. 9. Since a plurality of probability models are used to compose a ladder box cdf, it can be approximately equivalent to the original cdf. If the number of probability models M is large enough, then consider multiple probabilitiesCumulative value of rate modelIt is the true value.

For the case where X is a vector, first, let X follow have an average value of Sum covariance matrix of Q _c Is a gaussian distribution of (c). Given n-dimensional Gaussian distribution->Wherein Q is _c Is symmetrical and positive. For Q _c After sull decomposition:

Q _c ＝UΣU′

wherein U is an orthogonal matrix; sigma is an upper triangular matrix consisting of diagonal jordan blocks. If y-N (0) _n×1 ,I _n×n ) The method comprises the following steps:

from above, through a linear transformationDiscrete approximations can be obtained one after the other>All M quality points can be combined and converted. When->And->When no longer obeying the normal distribution of the standard, it is possible to obtain +.>And->Is a value of (2).

(b) The multiple models extend the kalman filtering process.

One cycle of the multi-modal expansion kalman filter for a single sensor can be summarized as follows:

wherein,representing a state vector predicted by the ith linear model at the moment k; />Representing the prediction error of the ith linear model at time k,/for the linear model at time k>Covariance of prediction error for ith linear model; />Representing the predicted state vector of the ith linear model at time k, Q _k-1 Representing the noise covariance. />Can be represented by E x _k-1 ∣z ^k-1 ]Obtained. />

And calculating combined multi-model prediction by adopting a weighted addition method:

wherein,the state vector prediction covariance matrix of the j-th model is represented.

And (5) measurement and update:

Wherein,representing measurement vector covariance,/">Representing the Kalman gain, I representing the identity matrix; r is R _k Representing a measurement noise covariance matrix; />Measurement vector, z, representing k-time of jth linear model prediction _k Measurement vector representing time k +.>Representing the measurement matrix with a specific calculation parameter of +.> Represents the expansion point state at time k of the jth linear model,/->

Representing the measurement error, v, of the jth linear model _k Gaussian white noise representing the linear model at time k;the state vector, the state vector error and the covariance of the state vector after the j-th linear model k time update are respectively represented.

It should be noted that the time update and the measurement update are partially different from the EKF. Additional itemsAnd->The addition of each of the correction terms can be regarded as a correction term of the Kalman filter. The correction term is a key to reducing the estimation error of the multimode extended kalman filter.

Model probabilities based on Bayesian criteria are updated as follows:

/>

wherein,representing a function which obeys a Gaussian distribution, in the above formula +.>And->The probabilities of the jth linear model after and before updating at time k are respectively represented.

And finally, calculating and combining multi-model updating:

the obtained updated motion state of the unmanned plane cluster target And outputting the detection result as the detection result of the k moment.

From each local prediction density and the first two momentsAnd->It is also possible to generate M quality points and use them in a first-order taylor expansion of the measurement model. One cycle of this process is shown in fig. 10.

Other than the technical features described in the specification, all are known to those skilled in the art. Descriptions of well-known components and well-known techniques are omitted so as to not unnecessarily obscure the present application. The embodiments described in the above examples are not intended to represent all the embodiments consistent with the present application, and various modifications or variations may be made by those skilled in the art without the need for inventive effort on the basis of the technical solutions of the present application while remaining within the scope of the present application.

Claims (6)

1. The detection tracking identification method based on the unmanned aerial vehicle cluster target is characterized by comprising the following steps of:

firstly, according to different types of targets in a current unmanned aerial vehicle cluster, executing a corresponding cluster target detection tracking mode, and predicting the motion state of the unmanned aerial vehicle cluster target; the group target detection tracking mode comprises the following steps of 1 to 4:

step 1), when the targets in the group are indistinguishable, using a centroid group tracking algorithm to detect and track the targets in the group;

Step 2) when the targets in the group are distinguishable and are non-rigid targets, detecting and tracking the group targets by using a group target and multi-target information interaction tracking algorithm;

step 3) when the targets in the group are distinguishable and rigid targets, if the group configuration is known, using an extended group target tracking method A to detect and track the group targets; the extended group target tracking method A obtains the motion centroid of each group of tracks based on an extended Kalman filter;

step 4) when the targets in the group are distinguishable and rigid targets, if the configuration of the group is unknown, using an extended group target tracking method B to detect and track the group targets; the extended group target tracking method B carries out group target detection tracking by using a Bayesian estimation framework based on an extended state of a random matrix, wherein the random matrix is described by inverse Weisal distribution, gaussian distribution description is used for a target motion state, and gamma distribution description is used for a measurement rate state;

secondly, constructing a tolerance multi-model, as described in step 5; updating the predicted motion state of the unmanned aerial vehicle cluster target by using a tolerance multi-model, and taking the updated result as a final tracking and identifying result;

step 5) constructing a tolerance multi-model by using a multi-model extended Kalman filter, comprising the following steps:

The detection tracking recognition problem of the unmanned aerial vehicle cluster target is a nonlinear estimation problem, and a posterior probability distribution function of the motion state of the unmanned aerial vehicle cluster target is set to obey Gaussian distribution; setting M particles, wherein each particle uses a linear model to approximate the Gaussian distribution, and the probability accumulated value of the M linear models is approximate to a true value; m is a positive integer;

and calculating the corresponding expansion states of the predicted movement states of the unmanned aerial vehicle cluster targets by using M linear models, and carrying out weighted summation on the M expansion states to obtain updated movement states of the unmanned aerial vehicle cluster targets.

2. The method of claim 1, wherein in the step 3, the processing step of the extended group target tracking method a includes the steps of:

3.1 A step of prediction: calculating an estimated value of the current track centroid state by a prediction step of an extended Kalman filter according to the motion centroid of the track of the previous time group;

3.2 Gating and scoring steps: inputting all sensor measurement values at the current moment in the association process, screening out measurement falling into the current track wave gate through a gating function, scoring the screened measurement, if one measurement falls into a plurality of track wave gates at the same time, only reserving an optimal score for the measurement, and associating the measurement with the track with the optimal score; the gating function is obtained according to prior information of the group target formation, and represents the amount of information which exists in the current track and can accept the current uncertain state;

3.3 A) distribution step: performing clustering processing on the measurements which are not associated with any track, and distributing a new track for the measurements when the track generation condition is met;

3.4 Updating step): updating each track by an updating step of the extended Kalman filter, and calculating a group innovation amount and a group innovation covariance matrix according to the estimated amount and the associated measurement of the mass center state of each track to obtain a motion mass center of the group at the current moment;

wherein the group innovation covariance matrix

Wherein S is _apr (n) represents unmanned cluster state vector at time n, J _H Representing the Jacobian matrix, the upper corner T representing the transpose operation, P _apr (n) a priori estimate matrix of the state vector covariance matrix at time n, R _G Representing a measurement error covariance matrix, C _D Representing a track variance estimation matrix;

for each group track, calculating distance function for all measurements obtained at the update timeThe distance function represents the new information added for the current group track by the new measurement; limiting the amount of information accepted by chi-square detection requires that +.>d _max And (5) for presetting gating parameters, acquiring according to prior information of the group target formation.

3. The method of claim 1, wherein in the step 4, when describing the state of the extended object using an ellipse approximation, the processing step of predicting the state of the extended object using the extended group object tracking method B includes the steps of:

4.1 a) setting a probability density function of an elliptic expansion target state at the time k-1 to obey gamma Gao Sini Weisal distribution;

the method comprises the steps that an ellipse expansion target state is set, wherein the ellipse expansion target state comprises a measurement rate state, a motion state and an expansion state of an expansion target;

4.2 a) before making the state prediction, the following assumptions are made:

a.1: the measurement rate prediction state of the expansion target is irrelevant to the motion state and the expansion prediction state, and follows a heuristic prediction form;

a.2: the motion prediction state of the expansion target is irrelevant to the expansion state;

a.3: the extension prediction state of the extension target depends on the motion state;

4.3 a) performing state prediction, comprising:

let k moment sensor received measurement set w= { z _k,1 ,z _k,2 ,...,z _k,|W| W represents the number of measurements in the set W, z _k,j Represents the j-th measurement and the k-moment measurement Z _k Likelihood functions of (2) are as follows:

wherein, gamma _k ,x _k ,X _k Respectively representing the measurement rate state, the motion state and the expansion state of the expansion target at the moment k,as a poisson probability density function, < >>As a gaussian probability density function +.>As a Weisalt probability density function, H _k Representing the measurement matrix, R _k Representing measurementsNoise, eta represents the adjustment parameter, < >>And->Respectively represent the mean and variance vectors, m of the set W _k|k-1 The prediction parameter is the motion state at the moment k;

the probability density function of the posterior state of the extended target at time k is as follows:

wherein, by the moment k, the measurement set Z ^k ＝Z _k ∪Z ^k-1 ，Z ^k-1 A measurement set from the tracking start to the k-1 moment; measuring rate state gamma _k Obeying gamma distribution, alpha _k 、β _k Gamma distribution parameters of the measurement rate state at the moment k; state of motion x _k Obeying Gaussian distribution, m _k 、P _k Is a gaussian distribution parameter of the motion state; extended state X _k Obeying the inverse Weisalde distribution, v _k 、V _k Is the inverse Weisauter distribution parameter of the expansion state;

updating gamma distribution parameters of the measurement rate state: alpha _k ＝α _k|k-1 +|W|,β _k ＝β _k|k-1 +1；

Updating the Gaussian distribution parameters of the motion state: m is m _k ＝m _k∣k-1 +K _k G _k ,P _k ＝P _k∣k-1 -K _k H _k P _k∣k-1 ；

Intermediate parameter K _k And G _k The following are provided:

wherein, the upper corner mark T represents transposition and the upper corner mark-1 tableAn inversion matrix is shown;an expected value representing an extended state;

updating the inverse Weisal distribution parameters of the expansion state: v _k ＝v _k∣k-1 +|W|,

Intermediate parameter N _k The method comprises the following steps:

wherein B is _k Representing covariance of the noise vector; d is the vector dimension of the motion state; alpha _k|k-1 、β _k|k-1 Is the prediction parameter of the measurement rate state at the moment k; m is m _k|k-1 、P _k|k-1 Is a prediction parameter of the motion state at the moment k; v _k|k-1 、V _k|k-1 Is a prediction parameter of the expansion state at the moment k; the prediction parameters of the motion state and the expansion state are calculated by using an expansion Kalman filtering algorithm;

And calculating the motion state of the expansion target by using a probability density function of the posterior state of the expansion target at the moment k.

4. The method of claim 1, wherein in the step 4, when the extended target state is described using a plurality of ellipse approximations, the processing step of performing the state prediction using the extended group target tracking method B includes the steps of:

4.1 b) modeling a non-elliptical expansion target state into a set comprising N sub-elliptical expansion target states, wherein the sub-elliptical expansion target states are mutually independent, and each sub-elliptical expansion target state comprises a measurement rate state, a motion state and an expansion state of the expansion target; n is an integer greater than 1;

4.2 b) before making the state prediction, the following assumptions are made:

a.4: the number of sub-ellipses N is known and remains unchanged during the target tracking process;

a.5: the relationship between the sub-ellipses and the measurements is unknown;

4.3 b) performing state prediction, comprising:

based on the full probability theory and the Bayesian estimation framework, the probability density function of the posterior state of the ith sub-ellipse expansion target at the k moment is:

wherein,expanding the target state for the ith sub-ellipse, Z ^k Is a measurement set from the tracking start to the k moment; />To measure the total number of all possible association events with a child ellipse,/for example >n _k For measuring the number of->Representing the first associated event therein, +.>Representing associated event->Probability of (2); />Representing associated event->The update state of the ith sub ellipse expansion target is set down; z is Z _k Measuring a vector at the moment k; />Representing associated event->A set of measurements associated with the ith sub-ellipse;representing associated event->The prediction state of the ith sub ellipse expansion target is set down;

given an associated eventI.e. < ->The predicted state of the factor ellipse expansion target is irrelevant to the related event, there is +>

Based on Bayesian estimation theory, probabilityThe unfolding is as follows:

wherein,normalization constant for Bayesian estimation,>if the correlation event between the measurement and the sub-ellipse does not have any prior information, then +.>The factors ellipses are mutually independent, and the likelihood function is measured>Is obtained by multiplying likelihood functions among different sub-ellipses, and by utilizing an edge integral theorem,the unfolding is as follows:

and calculating the motion state of each expansion target by using a probability density function of the posterior state of each sub-ellipse expansion target at the k moment.

5. The method of claim 1 wherein in step 5, the number of particles M and the location of each particle are obtained as follows:

assuming that the nonlinear estimation problem is correspondingly expressed as a real continuous model, and knowing the cumulative distribution function F of the real continuous model _c (x) Let the cumulative distribution function of the discrete random variable be F _d (x) X is a function variable, and under the premise of given distance measurement tolerance epsilon, F is solved by using constraint conditions _d (x) The following are provided:

the required minimum model number M is a minimum integer not less than 1/2 epsilon;

let the ith particle be denoted as m _i The location of the particles and the probability at the particles are determined by M as follows:

wherein arg represents a value of x satisfying an expression in brackets;

indirect m acquisition by dividing digital function _i Obtaining m using statistical function norm in matlab _i The following are provided:

step A: generating a scalar probability mass function, and approximating a standard normal distribution by using M mass points;

and (B) step (B): order theThen, norm v (p) _i ) Correspondingly get m _i ,i＝1,2…M；

Or obtaining m by searching standard normal distribution fractional number table in advance with known particle number _i 。

6. The method according to claim 1 or 5, wherein in step 5, the expansion state of each linear model prediction is calculated first as follows:

let the motion state of unmanned aerial vehicle cluster target at k moment be x _k Solving the mean value Ex of posterior distribution of motion states at the first two moments _k-1 ∣z ^k-1 ]Calculating the expansion state of the ith linear model prediction at the k momentWherein x is _k-1 Representing the motion state at time k-1, z ^k-1 A measurement set from the tracking start to the k-1 moment;

calculating the estimated error of the ith linear model at time kCovariance of error->The following are provided:

wherein w is _k-1 Gaussian white noise, Q, representing a linear model at time k-1 _k-1 Representing the noise covariance at time k-1; intermediate parameters Representing the expansion point state of the ith linear model at the moment k-1, wherein the function f is a state transfer function under a nonlinear system; />P _k-1|k-1 The motion state and covariance of the posterior probability at the moment k-1 respectively; the upper corner mark T represents transposition;

then, a weighted addition method is adopted to calculate the prediction result of the combined multi-model as follows:

wherein,P _k|k-1 respectively representing state vectors and covariance of the combined M linear model predictions at k time points; />Representing the probability of the ith linear model at time k-1; />Representing a state vector prediction covariance of a j-th model;

the measurements of each linear model were updated as follows:

wherein,represents the j thMeasurement of k time, z, of a linear model prediction _k Measurement vector representing time k +.>Representing the measurement matrix, the intermediate parameter-> Represents the expansion point state at time k of the jth linear model,/->Representing the measurement error, v, of the jth linear model _k Gaussian white noise representing the linear model at time k; />Representing measurement vector covariance; r is R _k Representing a measurement noise covariance matrix; />Representing the kalman gain;

further obtaining the state vector calculated by the j-th linear model after updatingStatus vector error->Covariance of state vector->The following are provided:

wherein I represents an identity matrix;

updating probabilities of the linear model based on bayesian criteria, comprising: setting a functionProbability of jth linear model updated +.>

Wherein,and->Respectively representing the probability of the jth linear model after and before updating at the time k;

finally, the motion state of the unmanned aerial vehicle cluster target is updated by the combined multi-model to obtainThe following are provided:

to be obtainedAnd outputting the motion state of the unmanned aerial vehicle cluster target at the identified k moment.