CN110780269B - Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity - Google Patents
Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity Download PDFInfo
- Publication number
- CN110780269B CN110780269B CN201910949144.4A CN201910949144A CN110780269B CN 110780269 B CN110780269 B CN 110780269B CN 201910949144 A CN201910949144 A CN 201910949144A CN 110780269 B CN110780269 B CN 110780269B
- Authority
- CN
- China
- Prior art keywords
- target
- gaussian
- new
- component
- label
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 238000005259 measurement Methods 0.000 claims abstract description 57
- 238000000605 extraction Methods 0.000 claims abstract description 24
- 230000004083 survival effect Effects 0.000 claims abstract description 20
- 238000012545 processing Methods 0.000 claims abstract description 11
- 238000002372 labelling Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 60
- 230000001052 transient effect Effects 0.000 claims description 49
- 238000001514 detection method Methods 0.000 claims description 27
- 238000007781 pre-processing Methods 0.000 claims description 7
- 238000013138 pruning Methods 0.000 claims description 6
- 239000000203 mixture Substances 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 241000764238 Isis Species 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 230000003252 repetitive effect Effects 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 2
- 230000003044 adaptive effect Effects 0.000 claims 4
- 238000005516 engineering process Methods 0.000 abstract description 6
- 238000002474 experimental method Methods 0.000 description 16
- 238000004088 simulation Methods 0.000 description 9
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 4
- 238000012986 modification Methods 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 239000013589 supplement Substances 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 238000012790 confirmation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000009469 supplementation Effects 0.000 description 1
- 230000001502 supplementing effect Effects 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
- G01S7/418—Theoretical aspects
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/66—Radar-tracking systems; Analogous systems
- G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
- G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
- G01S13/726—Multiple target tracking
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
The invention discloses an explicit multi-target tracking method based on a GM-PHD filter under self-adaptive newness intensity, which has the core technology that Gaussian components are divided into three categories by labeling, and different categories of Gaussian components adopt different state extraction methods to shield clutter in a survival target prior density area or interference of other target measurement; meanwhile, the wave gate method is used for obtaining the prior information of the new target, so that the multi-target explicit track under the condition that the new target is unknown a priori is obtained, and additional association processing is not needed to be carried out on the multi-target state estimation. Compared with a basic GM-PHD filter known by a new object in a priori manner, the precision is slightly good, and the real-time performance is slightly low; compared with the GLMB filter known by the new target a priori, the precision is low, and the real-time performance is excellent. In a multi-target tracking scene with unknown new targets in a priori manner, which requires both real-time performance and precision, the filter designed by the invention is an ideal choice.
Description
Technical Field
The invention relates to an explicit multi-target tracking method based on a Gaussian mixture probability hypothesis density (GM-PHD) filter under a new target priori unknown scene, and belongs to the technical field of signal processing.
Background
The multi-target jointly estimates the number of a plurality of targets and the states of the targets from a series of measured data on line, and is widely applied to the military and civil fields, such as radar multi-target tracking. With the continuous improvement of the demand of people on the radar function, the application scene of the radar is increasingly complex. Such as signal-to-noise ratio, low signal-to-noise, and dense targets, which increases the false alarm probability and reduces the detection probability of the target, thereby directly affecting the accuracy of target state-measurement association and reducing the accuracy of state estimation. The related problems of shielding various interferences in multi-target tracking and realizing accurate multi-target state estimation and flight path become a hot spot of research in the field at home and abroad. Currently, joint probabilistic data correlation, multi-hypothesis correlation, and Random Finite Set (RFS) are most applied.
In recent years, RFS has gained wide attention in multi-target tracking applications, since it avoids traditional data correlation. Based on RFS, various approximate implementations of bayesian multi-objective filters have been proposed, mainly including Probability Hypothesis Density (PHD) filters, potential probability hypothesis density filters, and multibarry filters. Although these filters do not provide multi-target tracks, they are still widely used.
The PHD filter iterates only the first moment of the multi-target a posteriori probability. From the practical application perspective, the PHD filter is particularly suitable for scenes with high real-time requirements. But it does not provide explicit multi-target tracks. And, the basic PHD filter is based on the assumption that the new prior object is known a priori. For PHD filters that cannot provide explicit multi-target tracks, the methods proposed in recent years require additional track management. When targets are very close to each other or in dense clutter or missing detection, the methods have difficulty in giving correct multi-target tracks. In recent two years, explicit multi-target tracking algorithms have been proposed in the literature based on Sequential Monte Carlo (SMC) approximate implementation of a PHD filter and Gaussian Mixture (GM) approximate implementation of a PHD filter, respectively, but they are based on the assumption that the target is known a priori. This makes it applicable only to some special scenarios, such as airport. Although some solutions are proposed for the problem that the new targets are not known a priori, the explicit track management problem of multi-target tracking is not solved at the same time.
In recent years, labeled RFS-based labeled multi-bernoulli force (GLMB) filters have been proposed that can significantly improve the accuracy of multi-target state extraction and provide explicit multi-target trajectories. However, such filters fuse data correlation techniques, so that the computational burden is increased. They are difficult to apply in real-time demanding scenarios. Furthermore, such filters are also based on the assumption that the new object is known a priori.
Therefore, a filter which is suitable for a new object in an unknown scene a priori, relatively small in calculation amount and capable of providing an explicit multi-target track has to be designed.
Disclosure of Invention
Aiming at the technical problem that multi-target track management under a new target priori unknown scene is not solved by a multi-target tracking technology based on a GM-PHD filter, the invention aims to provide an explicit multi-target tracking algorithm based on the GM-PHD filter under the new target priori unknown scene, and the core technology of the explicit multi-target tracking algorithm is that Gaussian components are divided into three categories by labeling, and different categories of Gaussian components adopt different state extraction methods to shield clutter in a survival target priori density region or interference of other target measurement; meanwhile, the wave gate method is used for obtaining the prior information of the new target, so that the multi-target explicit track under the condition that the new target is unknown a priori is obtained, and additional association processing is not needed to be carried out on the multi-target state estimation. Compared with a basic GM-PHD filter known by a new object in a priori manner, the precision is slightly good, and the real-time performance is slightly low; compared with the GLMB filter known by the new target a priori, the precision is low, and the real-time performance is excellent. In a multi-target tracking scene with unknown new targets in a priori manner, which requires both real-time performance and precision, the filter designed by the invention is an ideal choice.
The invention adopts the following technical scheme for solving the technical problems:
the invention provides an explicit multi-target tracking method based on a GM-PHD filter under self-adaptive new generation strength, which is based on the assumption that new generation targets are unknown a priori and provides two one-to-one principles suitable for multi-target state extraction of the GM-PHD filter: the non-new gaussian component of a cluster of the same label corresponds to only one measurement, and the measurement also corresponds to only the cluster of gaussian components.
The method comprises the following steps:
When k is 0, determining the label r of the Gaussian componentmaxEach new determined track is built by taking 0 as a base number, rmaxAdding 1; the label of the second order Gaussian component is given by rmax2=Vun2As the base number, r is added every time a second-level transient track is createdmax2 Adding 1; labeling of first order transient Gaussian components with rmax1=VunIs radix, every newEstablishing a first-level transient track rmax1Adding 1;
When k is 1, there is no prior information about the target, and the gaussian component prediction step is not performed;
when k is 2, performing a Gaussian component prediction step, wherein only the prediction of the new multi-target information is included;
when k is larger than or equal to 3, a Gaussian component prediction step is executed, wherein the Gaussian component prediction step comprises the prediction of survival multi-target information and the prediction of new multi-target information;
and predicting survival multi-target information: the posterior intensity D of the surviving multiple targets is obtained at the known k-1 momentk-1|k-1(x) I.e. a priori information about surviving objects is known at time kIs composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,denotes the Gaussian density, x is the state of a single target, x ═ xp xv yp yv]T,[xp yp]TIndicates the position of a single target, [ x ]v yv]TRepresenting the speed of a single target;andthe weight, the mean, the covariance matrix and the label of the j0 th component at the moment k-1 are respectively; predicting survival multi-target prior information at the moment k to obtain survival multi-target prediction strength D at the moment kS,k|k-1(x),Is composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,andthe weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k, wherein p isS,kProbability of survival of a single target, Fk-1And Qk-1Respectively a single-target state transition matrix and a process noise covariance matrix;
and predicting the new multi-target information: the new multi-target intensity gamma at the k moment is obtained at the k-1 momentk(x),Is composed of Jγ,kParameter set of gaussian componentsApproximation of, whereinAndrespectively is the weight, the mean, the covariance matrix and the label of the ith new target component; predicting the new multi-target information at the k moment to obtain new multi-target prediction strength D at the k momentγ,k|k-1(x),Is composed of Jγ,kParameter set of gaussian componentsAnd (ii) an approximation, wherein,andrespectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
mergingAndobtaining a Gaussian component parameter setJk|k-1=Jk-1+Jγ,k,Andrespectively obtaining the weight, the mean, the covariance matrix and the label of the jth prediction component, and obtaining the multi-target prediction strength D at the moment kk|k-1(x),
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radarAnd the prediction Gaussian component parameter set obtained in step 2Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtainedWherein,is the xy-plane ze1 measurements obtained by radar, is the ze number effective measurement obtained by adopting a wave gate method, and | is the aggregate potential;
(3b) execution of Zk,resi=Zk-Zk,efObtaining a residual measurement setWhen k is 1, no gaussian component parameter set is predicted, then Zk,efIs empty, Zk,resi=Zk;
Step 4, updating Gaussian components
When k is 1, there is no prediction gaussian component parameter set, Zk,efIf the value is null, the step of updating the Gaussian component is not executed;
when k is more than or equal to 2, because the gaussian component prediction step is executed, the gaussian component updating step exists, specifically:
according to the obtained Zk,efAndupdating the multiple posterior intensities D at the time kk|k(x) WhereinpD,kThe detection probability of a single target at the moment k is (1-p)D,k)Dk|k-1(x) Is composed of Jk|k-1Parameter set for updating Gaussian componentApproximation;wherein,andrespectively according to the effective measurement zk,efAnd the weight, mean, covariance matrix, and label of the updated gaussian component obtained from the jth predicted component,λ is the average number of clutter obtained by the radar at time k, c (z)k) Arbitrary measurements z obtained for radarkThe probability density of the spatial distribution of (a),Hkfor the measurement matrix, RkMeasuring a noise covariance matrix; i is an identity matrix and is a matrix of the identity,obtaining a multi-target posterior intensity D at an approximate k momentk|k(x),Is composed of Jk|k-1(1+|Zk,ef|) parameter sets for updating gaussian componentsIn the approximation that the difference between the first and second values,andthe weight, mean, covariance matrix and label of the j1 th updated gaussian component respectively;
step 5, multi-target state extraction
When k is 1, no Gaussian component is updated, and the multi-target state extraction step is not executed;
when k is larger than or equal to 2, executing a multi-target state extraction step, and extracting a plurality of single-target state estimations with identity identifications from the updated Gaussian components, wherein the method specifically comprises the following steps:
(5a) setting four query matrixes Uw、Um、UP、UlAre respectively Jk|k-1(1+|Zk,efI) matrix of dimension, J obtained in step 4k|k-1(1+|Zk,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in Uw、Um、UPAnd UlIn (1), wherein,are stored to U in sequencewThe first column of (a) is,are stored to U in sequencemThe first column of (a) is,are stored to U in sequencePThe first column of (a) is,are stored to U in sequencelThe first column of (1); press ze ═ 1., | Zk,efL, stored in sequence byUpdated parameters of Gaussian component to Uw、Um、UPAnd UlIn the specification:are stored to U in sequencewThe ze +1 th column of (a),are stored to U in sequencemThe ze +1 th column of (a),are stored to U in sequencePThe ze +1 th column of (a),are stored to U in sequencelZe +1 column (b);
(5d) Get UwMaximum value of all elements from the second column to the last columnIf this maximum value is presentLess than basic weight threshold w of state extraction1Directly jumping to the step (5 h); otherwise, according to the maximum valueAt a different position, from UlTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label setIs provided withThe system is used for storing a part of column number sets of effective measurement in the query matrix;
(5e) push buttonFor all labels in turn areThe processing of the gaussian component specifically includes:
(5e.1) the weight is givenThe label isIs at UlIn (1), the row number and column number are denoted as rn(i1)And cn(i1)Is recorded in UmElements in corresponding positions inIs m(i1),UlWherein all tags areRow number set ofRecord as
(5e.2) judging the tagAccording to the attribute, the corresponding operation is carried out on the attribute, and the following three conditions exist:
A. when in useThenA label for a new component, wherein Vnew>Vun2The method comprises the following specific operation steps:
firstly, whenEstablishing a two-level transient track, wherein w2For the second-level transient weight threshold, the specific operation is as follows:
rmax2=rmax2+1
wherein r ismax2=rmax2+1, updating the label of the secondary transient track;note the bookMiddle (rn)(i1)The elements of a row in all columns areWill be provided withAll elements in (1) are changed to rmax2(ii) a Note the bookMiddle rn(i1)The elements of a row in all columns areNote UwMiddle rn(i1)Rows and cn(i1)The elements of the column areFirstly, the method is carried outAll elements in (1) are set to zero and thenMiddle rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoRespectively counting the transient tracks rmax2Secondary transient set stored at time kAnd first order transient setPerforming the following steps;
rmax1=rmax1+1
wherein r ismax1=rmax1+1, updating the label of the primary transient track; will be provided withAll elements in (1) are changed to rmax1(ii) a Firstly, the method is carried outAll elements in (1) are set to zero and thenMiddle rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoThe number r of transient tracksmax1First order transient set stored at time kPerforming the following steps;
③ the number cn of the row(i1)Is stored toIn CI, i.e. CI ═ CI cn(i1)](ii) a Changes are made toIs rn(i1)I.e. by
B. When in useThenThe method is a label of a transient component, and comprises the following specific steps:
firstly, whenThen, further on,for labels of the two types of transient components, the following steps are sequentially executed:
secondly, if k is not less than 3 andfirstly, a new determined track r is establishedmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)];
② whenThen, further on,for a class of transient classified tags, the following steps are performed in sequence:
Secondly, if k is more than or equal to 4,And isFirstly, a new determined track r is establishedmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)];
C. When in useThenTo determine the label of the component, from labelsThe method comprises the following specific steps of:
let w3Is a weight threshold value of state preprocessing and satisfies w1≤w3<w2;
Wherein,is at UwMiddle labelAll the rows correspond to the row number of the maximum value in all the elements of the first column; | | b1-b2| | represents the euclidean distance between b1 and b 2; note UmTo (1) aThe elements in row number 1 areKnown zero mean metrology noise covariance matrix Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]T,d(1)=||σ-0||;
Otherwise, executing the step II;
② will m(i1)As a single target state estimate, store it to a set of state estimatesIn, i.e.Will be provided withAnd k are stored separately to the track tag setAnd set of occurrence momentsIn, i.e.Andat the same time, adding cn(i1)Store to CI, i.e. CI ═ CI cn(i1)];
(iii) correction of weight of Gaussian component
Firstly, calculate m one by one(i1)Euclidean distance to the same tag component:
and correcting the weight according to the conditions:
Or
Or
Wherein, UwAnd UmTo (1) aThe elements in the row at the ze1 column are respectively marked asAndattenuation factor a1 ═ 1/3;
(5e.3) based on two "one-to-one" principles, U is dividedwThe middle row has the number ofIs set to zero, i.e. is performed(5f) Is executed completelyAfter circulation, the elements of CI are used as a column number set, and then the elements are listed in U according to two 'one-to-one' principleswAll the elements of any row of (1) are set to zero;
(5g) returning to the step (5 d);
(5h) respectively to be provided withAndj in (1)k|k-1(1+|Zk,efI) conversion of elements column by column to 1 × (J)k|k-1(1+|Zk,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
(5i) Connecting state estimates with the same identity at different moments to obtain an explicit multi-target track;
(6a) Cut out a succession of nnostThe method comprises the following specific steps of:
(6a.1) takingIn non-repetitive tab sets, from which to deleteThe label set of the component of which the state estimation is not obtained at the moment k is obtained
(6a.2) pressingSequentially discriminatingIn thatk1=k-nnotrThe number of occurrences in k-1: if it is notAt each time k1If all appear in (1), then execute
(6b) All the weights are larger than a small weight threshold value wetThe sequence number of the Gaussian component of (1) is recorded in I1, i.e.
(6c) Let Jk|k=|I1|,By analogy, obtainWherein,andthe weight, the mean, the covariance matrix and the label of the j3 th Gaussian component respectively;
step 7, obtaining parameter sets of Gaussian components based on the pruning in step 6Performing Gaussian component combination to makeIs a set of sequence numbers of the Gaussian components, andthe method comprises the following specific steps:
(7b) Taking andthe sequence number set of the co-tagged and combinable components is I3,wherein d isThIs a merged distance threshold;
(7c) all gaussian component parameters in the sequence number set I3 are merged, obtaining a new Gaussian component;
(7d) i2 ═ I2-I3, and return to step (7 a);
Step 8, generating new target intensity
According to Zk,resiPress ze2 ═ 1., | Zk,resiObtaining the parameters of 2 th new components according to the following operations:
obtaining the prior intensity gamma of the new target at the k +1 momentk+1(x),Jγ,k+1=|Zk,resi|,i=ze2;γk+1(x) Is composed of Jγ,k+1Parameter set of gaussian componentsApproximation;
and 9, returning to the step 2 when k is equal to k + 1.
Further, the specific steps in step 3a are:
(3a.1) known zero mean of the measured noise covariance matrix Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]TBased on this, calculating a gate threshold d (a) ═ a | | | σ -0|, where a is a confidence coefficient;
(3a.3)For the latest measurement set obtained at time k, ze1 ═ 1kL, sequentially judging eachIs determined by the attribute ofThen, press J ═ 1k|k-1One by one calculateTo each gaussian component meanIs a distance ofIf d is(ze1,j)D (a) or less, then judgingNearby Gaussian component, seeFor effective measurement, add it to the effective measurement set Zk,efThen jump out of the loop and execute the pairAnd (4) judging.
Further, step 1, Vun<Vun2。
Further, Vun=500,Vun2=1000。
Further, step 7 further includes step (7f), specifically:
the maximum Gaussian component number of the iteration to the k +1 moment does not exceed Jmax1When is coming into contact withWhen it comes toBefore taking Jmax1The Gaussian component with the maximum weight is recorded asJk=Jmax1;When in use, willRelabeling toThereby obtaining the posterior intensity of survival multiple targets
Compared with the prior art, the invention adopting the technical scheme has the following technical effects:
(1) obtaining effective measurement for updating and predicting Gaussian components by using a wave gate technology to obtain residual measurement; based on the residual measurement, obtaining prior information of a new target at the next moment;
(2) labels of the Gaussian components are divided into three categories, and the transient components are further divided into two stages, so that the accuracy of establishing the new target track is improved;
(3) based on the labels of Gaussian components, two one-to-one principles and corresponding processing schemes suitable for multi-target state extraction of a GM-PHD filter under the scene of unknown strength of a new target are provided, short-distance clutter or interference of other target measurement can be shielded, state estimation of each single target is obtained synchronously with the labels, and explicit multi-target flight paths can be obtained without additional correlation processing.
Drawings
FIG. 1 is a diagram of a lookup matrix structure according to the present invention;
FIG. 2 is a flow diagram of multi-target state extraction of the present invention;
FIG. 3 is a flight path of each target in the x-axis and the y-axis obtained in a single experiment with a clutter number of 20 and a detection probability of 0.92;
FIG. 4 is a plot of the x-axis trajectory over time for each target obtained in a single experiment with a clutter number of 20 and a detection probability of 0.92;
FIG. 5 is a y-axis trajectory of each target over time obtained in a single experiment with a clutter number of 20 and a detection probability of 0.92;
FIG. 6 is the average OSPA distance for each filter obtained from 100 experiments with a clutter number of 20 and a detection probability of 0.92;
FIG. 7 shows the average calculation time of each filter obtained in 100 experiments when the number of clutter is 20 and the detection probability is 0.92;
FIG. 8 is a plot of the x-axis and y-axis trajectories of each target obtained in a single experiment with a clutter number of 50 and a detection probability of 0.98;
FIG. 9 is a plot of the x-axis trajectory over time for each target obtained in a single experiment with a clutter number of 50 and a detection probability of 0.98;
FIG. 10 is a y-axis trajectory of each target over time obtained in a single experiment with a clutter number of 50 and a detection probability of 0.98;
FIG. 11 is the average OSPA distance for each filter obtained from 100 experiments with a clutter number of 50 and a detection probability of 0.98;
fig. 12 shows the average calculation time of each filter obtained by 100 experiments when the number of clutter is 50 and the detection probability is 0.98.
Detailed Description
The technical scheme of the invention is further explained in detail by combining the attached drawings:
aiming at the technical problem that multi-target track management under the condition that a new target is unknown a priori is not solved by a multi-target tracking technology based on a GM-PHD filter, the invention aims to provide an explicit multi-target tracking method based on a Gaussian mixture probability hypothesis density filter under the condition that the new target is unknown a priori, and the core technology of the explicit multi-target tracking method is to remove clutter by using a wave gate method and obtain the prior information of the new target based on the clutter; labeling all Gaussian components, and respectively storing the four types of parameters of all Gaussian components to a search matrix according to corresponding sequences; by utilizing the search matrix, the interference of the clutter in the high prior density area of the target is shielded by judging the labels and the weight of the Gaussian component, the multi-target state extraction is converted into a plurality of single-target state estimations which can provide unique unchanged identity labels, and therefore the explicit multi-target tracking without additional associated programs is obtained. Compared with a basic GM-PHD filter under a new priori known scene, the precision is slightly superior to that of the GM-PHD filter, and the real-time performance is slightly inferior to that of the GM-PHD filter; compared with the GLMB filter under the new-born priori known scene, the precision is inferior to the GLMB filter, and the real-time performance is far superior to the GLMB filter. The method is easy to realize in engineering, and is suitable for a multi-target tracking scene which has high real-time requirement, is unknown in new born prior and needs to provide an explicit track.
In order to achieve the technical purpose, the invention is realized by adopting the following technical scheme.
The invention is suitable for a scene with unknown newborn target priori, and provides two one-to-one principles suitable for multi-target state extraction of a GM-PHD filter: the non-new gaussian component of a cluster of the same label corresponds to only one measurement, and the measurement also corresponds to only the cluster of gaussian components.
Based on two 'one-to-one' principles of the proposed multi-target state extraction, the invention designs an explicit multi-target tracking algorithm based on a GM-PHD filter under a new target priori unknown scene, which comprises the following steps:
When k is 0, determining the label r of the Gaussian componentmaxEach new determined track is built by taking 0 as a base number, rmaxAdding 1; secondary transient gaussian component with rmax2=Vun2As the base number, r is added every time a second-level transient track is createdmax2Adding 1; first order transient Gaussian component with rmax1=VunAs a base number, r is a transient state track every time a new one is builtmax1Adding 1; wherein, Vun2And VunCan take different values in different scenes, and Vun<Vun2In this patent, V is specifiedun=500,Vun2=1000。
When k is 1, there is no prior information about the target, and the gaussian component prediction step is not performed;
when k is 2, performing a Gaussian component prediction step, wherein only the prediction of the new multi-target information is included;
when k is larger than or equal to 3, a Gaussian component prediction step is executed, wherein the Gaussian component prediction step comprises the prediction of survival multi-target information and the prediction of new multi-target information;
and (4) survival multi-target information prediction. The posterior intensity D of the surviving multiple targets is obtained at the known k-1 momentk-1|k-1(x) I.e. a priori information about surviving objects is known at time kIs composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,denotes the Gaussian density, x is the state of a single target, x ═ xp xv yp yv]T,[xp yp]TIndicates the position of a single target, [ x ]v yv]TRepresenting the speed of a single target; andthe weight, mean, covariance matrix, and label for the j0 th component at time k-1, respectively. Predicting survival multi-target prior information at the moment k to obtain survival multi-target prediction strength D at the moment kS,k|k-1(x),Is composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,andthe weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k, wherein p isS,kProbability of survival of a single target, Fk-1And Qk-1Respectively, a single target state transition matrix and a process noise covariance matrix.
And (4) predicting the new multi-target information. The new multi-target intensity gamma at the k moment is obtained at the k-1 momentk(x),Is composed of Jγ,kParameter set of gaussian componentsApproximation of, whereinAndrespectively is the weight, the mean, the covariance matrix and the label of the ith new target component; predicting the new multi-target information at the k moment to obtain new multi-target prediction strength D at the k momentγ,k|k-1(x),Is composed of Jγ,kParameter set of gaussian componentsAnd (ii) an approximation, wherein,andrespectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
mergingAndobtaining a Gaussian component parameter setJk|k-1=Jk-1+Jγ,k,Andrespectively obtaining the weight, the mean, the covariance matrix and the label of the jth prediction component, and obtaining the multi-target prediction strength D at the moment kk|k-1(x),
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radarAnd the prediction Gaussian component parameter set obtained in step 2Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtainedWherein,is the xy-plane ze1 measurements obtained by radar,is the ze number effective measurement obtained by adopting a wave gate method, and | is the aggregate potential.
The method specifically comprises the following steps:
(3a.1) known zero mean of the measured noise covariance matrix Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]TBased on this, the gate threshold d (a) is calculated as a | | | σ -0 |. Wherein a is a confidence coefficient.
(3a.3)For the latest measurement set obtained at time k, ze1 ═ 1kL, sequentially judging eachThe attribute of (2). Judgment ofThen, press J ═ 1k|k-1One by one calculateTo each gaussian component meanIs a distance ofIf d is(ze1,j)D (a) or less, then judgingNearby Gaussian component, seeFor effective measurement, add it to the effective measurement set Zk,efThen jump out of the loop and execute the pairAnd (4) judging.
(3b) Execution of Zk,resi=Zk-Zk,efObtaining a residual measurement setWhen k is 1, no gaussian component parameter set is predicted, then Zk,efIs empty, Zk,resi=Zk。
Step 4, updating Gaussian components
According to the obtained effective measurement set Zk,efAnd a set of Gaussian components for the prediction of time kUpdating the multiple posterior intensities D at the time kk|k(x) In that respect Known as pD,kThe detection probability of a single target at the time k is specifically as follows:
wherein,
wherein (1-p)D,k)Dk|k-1(x) From a set of Gaussian component parametersAnd (4) approximation. DD,k(x;zk,ef) Is composed of Jk|k-1Parameter set of gaussian componentsIn the approximation that the difference between the first and second values,andrespectively according to the effective measurement zk,efAnd the weight, the mean, the covariance matrix and the label of the updated Gaussian component obtained by the jth predicted Gaussian component are called the jth predicted Gaussian component as the parent component of the updated Gaussian component. λ is the average number of clutter obtained by the radar at time k, c (z)k) Arbitrary measurements z obtained for radarkProbability density of spatial distribution of (1), HkFor the measurement matrix, RkTo measure the noise covariance matrix, I is the identity matrix.
Accordingly, the multi-target posterior intensity at the k time is determined byParameter set of gaussian componentsAnd (4) approximation. Wherein,andthe weight, mean, covariance matrix, and label of the j1 th updated gaussian component, respectively.
Step 5, multi-target state extraction
Extracting a plurality of single-target state estimates with identity identifiers from the updated gaussian component, as shown in fig. 1, specifically includes:
(5a) setting four query matrixes Uw、Um、UP、UlIs Jk|k-1(1+|Zk,efI) matrix of dimension, J obtained in step 4k|k-1(1+|Zk,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in Uw、Um、UPAnd UlIn (1). FIG. 2 shows Uw、Um、UPAnd UlThe structure of the matrix. Wherein,are stored to U in sequencewThe first column of (a) is,are stored to U in sequencemThe first column of (a) is,are stored to U in sequencePThe first column of (a) is,are stored to U in sequencelThe first column of (2). Press ze ═ 1., | Zk,efL, stored in sequence byUpdated parameters of Gaussian component to Uw、Um、UPAnd UlIn the specification:are stored to U in sequencewThe ze +1 th column of (a),are stored to U in sequencemThe ze +1 th column of (a),are stored to U in sequencePThe ze +1 th column of (a),are stored to U in sequencelZe +1 column (b). Obviously, UPThe elements in each row are the same, UlThe elements of each row in the series are the same.
Storing the updated Gaussian component to Uw、Um、UPAnd UlIn this way, the parent component of each update component and the corresponding measurement are conveniently searched.
(5b) Creating four temporary matricesAndparameters for recording various types of modifications in the state extraction process:
(5c) set of state estimates for time kTrack label setSet of occurrence momentsEffectively measured at UwColumn number set in (1)
(5d) Memory matrix UwThe elements of all rows in the second to last columns areTaking the maximum value from the above, and recording the maximum value asUwThere may be more than oneWhen in useAnd (5h) jumping to the step. Otherwise, whenAccording toAt UwRow and column numbers at different positions in the slave UlTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label setWherein, w1Taking w as a basic weight threshold value extracted for the state according to a plurality of experimental simulation results1=0.02。
(5e) Push buttonFor all labels in turn areThe processing of the gaussian component specifically includes:
(5e.1) the weight is givenThe label isIs at UlIn (1), the row number and column number are denoted as rn(i1)And cn(i1)Is recorded in UmElements in corresponding positions inIs m(i1),UlWherein all tags areIs recorded as
(5e.2) judging the tagAccording to the attribute, the corresponding operation is carried out on the attribute. There are three cases:
A. when in useThenIs a label for the nascent component. VnewThe value of (A) can be different in different scenes, is far greater than the total number of tracks possibly appearing in the scenes, and satisfies Vnew>Vun2In this patent, V is specified through multiple simulation experimentsnew2000. The method comprises the following specific operation steps:
firstly, whenAnd establishing a secondary transient state track. Wherein, w2Is a secondary transient weight threshold, and is designated as w in the patent according to the results of multiple simulation experiments20.5. The specific operation is as follows:
rmax2=rmax2+1
wherein r ismax2=rmax2+1, updating the label of the secondary transient track; note the bookMiddle (rn)(i1)The elements of a row in all columns areWill be provided withAll elements in (1) are changed to rmax2(ii) a Note the bookMiddle rn(i1)The elements of a row in all columns areNote UwMiddle rn(i1)Rows and cn(i1)The elements of the column areFirstly, the method is carried outAll elements in (1) are set to zero and thenMiddle rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoRespectively counting the transient tracks rmax2Secondary transient set stored at time kAnd first order transient setIn (1).
rmax1=rmax1+1
wherein r ismax1=rmax1+1, updating the label of the primary transient track; will be provided withAll elements in (1) are changed to rmax1(ii) a Firstly, the method is carried outAll elements in (1) are set to zero and thenMiddle rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoThe number r of transient tracksmax1First order transient set stored at time kIn (1).
③ the number cn of the row(i1)Stored into CI, i.e. CI ═ CI cn(i1)](ii) a Changes are made toIs rn(i1)I.e. by
B. When in useThenThe method is a label of a transient component, and comprises the following specific steps:
secondly, if k is not less than 3 andfirstly, a new confirmation is establishedDetermining a flight path rmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)]。
Secondly, if k is more than or equal to 4,And isFirstly, a new determined track r is establishedmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)]。
C. When in useThenTo determine the label of the component. Slave signIs signed asThe method comprises the following specific steps of:
first, in order to reduce measurement noise, detection uncertainty, and interference caused by clutter, state estimation preprocessing is required. Let w3For the weight threshold of state preprocessing, it should satisfy w1≤w3<w2In this patent, w is specified through multiple simulation experiments3=0.2。
Wherein,is at UwMiddle labelAll the rows correspond to the row number of the maximum value in all the elements of the first column; b1-b2| | | represents the euclidean distance between b1 and b 2; note UmTo (1) aThe elements in row number 1 areKnown as zero meanMeasured noise covariance matrix of values Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]T,d(1)=||σ-0||。
② will m(i1)As a single target state estimate, store it to a set of state estimatesIn, i.e.Will be provided withAnd k are stored separately to the track tag setAnd set of occurrence momentsIn, i.e.Andat the same time, adding cn(i1)Store to CI, i.e. CI ═ CI cn(i1)]。
And correcting the weight of the Gaussian component. The larger the value w is, the more concentrated the posterior distribution of the single target is, and the smaller the radius of the effective area is; the smaller the value, the more dispersed the posterior distribution of the single target, and the larger the effective area radius should be.
For this purpose, m is first calculated one by one(i1)Distance d to the same label component(i2,ze1),
Then, the weight value is corrected according to the condition,
Or
Or
Wherein, UwAnd UmTo (1) aThe elements in the row at the ze1 column are respectively marked asAndin the effective area, the weight of the Gaussian component is unchanged; the gaussian component weights outside the region are attenuated by a factor of a 1. In the invention, a1 is 1/3, which is used for weakening the weight of the component outside the effective region and reducing the influence of various types of interference on the posterior information.
(5e.3) extracting U according to two 'one-to-one' principles of multi-target state extractionwThe middle row has the number ofIs set to zero, i.e. is performed
(5f) Is executed completelyAfter circulation, the elements of the CI are taken as a column number set, and the elements are listed in the U according to two 'one-to-one' principles of multi-target state extractionwAll of the elements of any row of (c) are zeroed.
(5g) And (5) returning to the step (5 d).
(5h) Respectively to be provided withAndj in (1)k|k-1(1+|Zk,efI) conversion of elements column by column to 1 × (J)k|k-1(1+|Zk,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
(5i) And connecting the state estimations with the same identity at different moments to obtain the explicit multi-target track.
When k is 1, the Gaussian component is not updated, and the step of multi-target state extraction is omitted;
and 6, pruning the Gaussian components. The method comprises the following specific steps:
(6a) cut out a succession of nnostThe next component without state estimation, n is specified in this patent nost6. The method comprises the following specific steps:
(6a.1) takingIn non-repetitive tab sets, from which to deleteThe tag contained in, i.e. executingTag set for components for which state estimates are not obtained at time k is obtained
(6a.2) pressingSequentially discriminatingIn thatk1=k-nnotrThe number of occurrences in k-1: if it is notAt each time k1If all appear in (1), then execute
(6b) All the weights are larger than a small weight threshold value wetThe number of the Gaussian component of (2) is recorded in I1, specified in this patentwet=10-5I.e. by
(6c) Let Jk|k=|I1|,By analogy, obtainWherein,andthe weight, mean, covariance matrix, and label of the j3 th gaussian component, respectively.
Step 7, parameter set based on Gaussian component obtained after pruningAnd carrying out Gaussian component combination. When multiple targets are adjacent or crossed, the posterior information among the targets is easy to generate interference during combination, and the Gaussian components of the targets with relatively small weights caused by missing detection or measurement noise are easy to be combined by the Gaussian components of the targets with larger weights, so that part of real information of the targets with small weights is lost. Thus, only components that are tagged identically can be merged. Order toIs a set of sequence numbers of the Gaussian components, andthe method comprises the following specific steps:
(7b) Taking andthe sequence number set of the components which are identical in label and can be combined is I3
Wherein d isThFor the combined distance threshold, d is specified in this patentTh=4。
(7c) All the gaussian component parameters in the sequence number set I3 are merged to obtain a new gaussian component:
(7d) i2 ═ I2-I3, and return to step (7 a).
(7f) The maximum Gaussian component number of the iteration to the k +1 moment does not exceed Jmax1In this patent, J is designatedmax1=150。
When in useWhen it comes toBefore taking Jmax1The Gaussian component with the maximum weight is recorded asJk=Jmax1;
And 8, generating new target intensity. According to Zk,resiPress ze2 ═ 1., | Zk,resiObtaining the parameters of 2 th new components according to the following operations:
thus, the prior intensity γ of the new object at the time k +1 is obtainedk+1(x),Jγ,k+1=|Zk,resi|,i=ze2;γk+1(x) Is composed of Jγ,k+1Parameter set of gaussian componentsAnd (4) approximation.
And 9, returning to the step 2 when k is equal to k + 1.
The explicit multi-target tracking method based on the GM-PHD filter under the self-adaptive new-generation strength is finished.
Examples
The effect of the invention is further verified and explained by the following simulation implementation.
1. The experimental conditions are as follows:
in a two-dimensional scene [ -10001000 ] mx[ -10001000 ] m, the equation of motion for each object is:
wherein x isk=[xp xv yp yv]T,[xp yp]TIs the position of the target at time k, [ x ]v yv]TIs the velocity of the target at time k. Δ ═ 1s for the sampling interval, σω=5m/s2. The object can appear or disappear at any time in the scene, and the survival probability pS,k0.99. In a GLMB filter under the scene that the new object is known a priori, a new object model is a labeled multi-Bernoulli parameterWherein PB=diag([10 10 10 10]T)2. In a basic GM-PHD filter under a scene in which the new object is known a priori, the intensity function of the new object appearance isWherein,Pγ,k=PB,i=1,2,3,4。
the measurement equation of the target is
Wherein upsilon isx,kAnd upsilony,kIs independent zero mean Gaussian white noise, and has mean square error of sigmax=10m、σy10 m. The noise is uniformly distributed in the monitoring area [ -10001000 [)]m×[-1000 1000]m, the number of clutter per frame is λ. In the basic GM-PHD filter, the clutter intensity is k ═ λ/20002. The filter provided by the invention is compared with a basic GM-PHD filter under a new object priori known scene. Because the filter provided by the invention uses the wave gate method to remove clutter before updating the Gaussian components, for fairness, the compared GM-PHD filter also uses the wave gate method to remove clutter, and the filter measures the gate probability p based on all predicted Gaussian components00.999. The maximum number of Gaussian components is JmaxTruncate threshold w to 100et=10-5Merging the thresholds dTh=4。
The optimal sub-pattern assignment (OSPA) distance is used to evaluate multi-target tracking performance, c is 100 and p is 1. Based on the same target trajectory, 100 MC experiments were run to obtain average performance, with measurements in each experiment being independent of each other.
2. Emulated content
3. And (3) simulation result analysis:
as can be seen from fig. 3 to 5 and fig. 8 to 10, the filter designed by the present invention can implement complex multi-target tracking under dense clutter, and give a correct multi-target track. In fig. 3-5, there is a pseudo-track, which occurs at k 59; in fig. 8 to 10, there are four pseudo-tracks, which occur at k-24, k-28, k-44 and k-85, respectively. However, the existence of these pseudo tracks does not affect the continuity of tracking of each determined track, and disappears soon.
In fig. 4, 5, 9 and 10, at the time steps where the new target appears, k is 1, k is 20, k is 40, k is 60 and k is 80, since the new target is unknown a priori, three or more time steps are required to determine a new target track, and therefore, compared with the other two filters in the scene where the new target is known a priori, the filter designed by the present invention shows the maximum OSPA distance, i.e. the maximum error, in fig. 6 and 11 at these times and three or more time steps thereafter.
However, followThe OSPA distance of the filter designed by the invention tends to be stable after the track is determined to be established. In FIG. 6, when the detection probability is relatively low, pD,kWhen the OSPA distance is equal to 0.92, the OSPA distance of the method is gradually equal to a GM-PHD filter based on a new target priori known scene until the OSPA distance is lower than the GM-PHD filter based on the new target priori known scene, namely the tracking performance is superior to the GM-PHD filter, and the contribution of clutter at a target short distance can be effectively shielded by the multi-target state extraction method designed by the invention; since the filter designed by the invention does not supplement the state estimation of each target in the process of missing detection, it can be seen in fig. 4 and 5 that there is a breakpoint in each determined track, which results in that the OSPA distance of the invention is smaller than that of the GLMB filter in the scene that the new target is known a priori, i.e. the tracking performance is inferior thereto, however, this does not affect the determination and tracking of each target track, see fig. 3. On the other hand, the state estimation of each target in the missed detection is supplemented by the GLMB filter in the new target priori known scene, so that when a target disappears, for example, when k is 71-73 in fig. 6, the state estimation of the disappeared target is still supplemented, and the OSPA distance of the filter is obviously larger than that of the other two filters. The filter designed by the invention does not have the phenomenon.
In FIG. 11, when the detection probability is relatively high, pD,kWhen the track is determined to be 0.98, the OSPA distance of the invention is smaller than the basic GM-PHD filter under the scene that the new target is known a priori, namely the tracking performance is better than the OSPA distance; likewise, it can be seen in fig. 9 and 10 that there is a breakpoint for each determined track, resulting in the OSPA distance of the present invention being smaller than the GLMB filter in the scenario where the new object is known a priori, i.e. the tracking performance is inferior thereto, however, this does not affect the determination and tracking of each target track as well, see fig. 8. In fig. 11, similarly, when k is 71-73, the OSPA distance of the GLMB filter in the new object a priori known scene is significantly larger than that of the other two filters. The filter designed by the invention does not have the phenomenon.
Comparing fig. 6 and fig. 11, it can be seen that the OSPA distance of the filter designed by the present invention is smaller than the OSPA distance of the filter designed by the present invention when the number of spurs λ is 50, that is, the tracking performance is rather improved because the corresponding detection probability is higher when λ is 50Rate pD,kCorresponding detection probability p when λ is greater than 20 at 0.98D,k0.92. Obviously, the filter designed by the invention does not supplement the state estimation of each target in the missing detection, and the OSPA distance of each target is greatly influenced. However, artificial supplementation will again cause the phenomenon of blindly supplementing the state estimation when the target disappears. Therefore, as long as the tracking of each target and the determination of the flight path are not influenced, the state estimation of each target in the missing detection can not be supplemented.
Comparing fig. 7 and fig. 12, the calculation time of the filter designed by the invention is slightly higher than that of the basic GM-PHD filter in the new-object a priori known scene, and is far lower than that of the GLMB filter in the new-object a priori known scene.
In conclusion, the simulation experiment verifies the correctness, the effectiveness and the reliability of the method. Therefore, when the real-time performance requirement is high, the new targets are not known a priori, and the tracks of all the targets are required to be provided, the filter designed by the invention is an ideal choice.
It will be apparent to those skilled in the art that various changes and modifications may be made in the present invention without departing from the spirit and scope of the invention; thus, if such modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include such modifications and variations.
Claims (5)
1. The explicit multi-target tracking method based on the GM-PHD filter under the self-adaptive new-generation intensity is characterized in that two one-to-one principles suitable for multi-target state extraction of the GM-PHD filter with Gaussian mixture probability hypothesis density are provided based on the hypothesis that new-generation targets are unknown a priori: the non-new Gaussian component of the same label of a certain cluster only corresponds to one measurement, and the measurement only corresponds to the Gaussian component of the cluster;
the method comprises the following steps:
step 1, initializing the number of various tracks
When k is 0, determining the label r of the Gaussian componentmaxEach new determined track is built by taking 0 as a base number, rmaxAdding 1; the label of the second order Gaussian component is given by rmax2=Vun2As the base number, r is added every time a second-level transient track is createdmax2Adding 1; labeling of first order transient Gaussian components with rmax1=VunAs a base number, r is a transient state track every time a new one is builtmax1Adding 1;
step 2, Gaussian component prediction
When k is 1, there is no prior information about the target, and the gaussian component prediction step is not performed;
when k is 2, performing a Gaussian component prediction step, wherein only the prediction of the new multi-target information is included;
when k is larger than or equal to 3, a Gaussian component prediction step is executed, wherein the Gaussian component prediction step comprises the prediction of survival multi-target information and the prediction of new multi-target information;
and predicting survival multi-target information: the posterior intensity D of the surviving multiple targets is obtained at the known k-1 momentk-1|k-1(x) I.e. a priori information about surviving objects is known at time kIs composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,denotes the Gaussian density, x is the state of a single target, x ═ xp xv yp yv]TWherein [ x ]p yp]TIndicates the position of a single target, [ x ]v yv]TRepresenting the speed of a single target;andthe weight, the mean, the covariance matrix and the label of the j0 th component at the moment k-1 are respectively; multi-objective prior living at time kPredicting the information to obtain survival multi-target prediction strength D at the moment kS,k|k-1(x),Is composed of Jk-1Parameter set of gaussian componentsAnd (ii) an approximation, wherein,andthe weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k, wherein p isS,kProbability of survival of a single target, Fk-1And Qk-1Respectively a single-target state transition matrix and a process noise covariance matrix;
and predicting the new multi-target information: the new multi-target intensity gamma at the k moment is obtained at the k-1 momentk(x),Is composed of Jγ,kParameter set of gaussian componentsApproximation of, whereinAndrespectively is the weight, the mean, the covariance matrix and the label of the ith new target component; predicting the new multi-target information at the k moment to obtain new multi-target prediction strength D at the k momentγ,k|k-1(x),Is composed of Jγ,kParameter set of gaussian componentsAnd (ii) an approximation, wherein,andrespectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
mergingAndobtaining a Gaussian component parameter setJk|k-1=Jk-1+Jγ,k,Andrespectively obtaining the weight, the mean, the covariance matrix and the label of the jth prediction component, and obtaining the multi-target prediction strength D at the moment kk|k-1(x),
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radarAnd the prediction Gaussian component parameter set obtained in step 2Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtainedWherein,is the xy-plane ze1 measurements obtained by radar, is the ze number effective measurement obtained by adopting a wave gate method, and | is the aggregate potential;
(3b) execution of Zk,resi=Zk-Zk,efObtaining a residual measurement setWhen k is 1, no gaussian component parameter set is predicted, then Zk,efIs empty, Zk,resi=Zk;
Step 4, updating Gaussian components
When k is 1, there is no prediction gaussian component parameter set, Zk,efIf the value is null, the step of updating the Gaussian component is not executed;
when k is more than or equal to 2, because the gaussian component prediction step is executed, the gaussian component updating step exists, specifically:
according to the obtained Zk,efAndupdating the multiple posterior intensities D at the time kk|k(x) WhereinpD,kThe detection probability of a single target at the moment k is (1-p)D,k)Dk|k-1(x) Is composed of Jk|k-1Parameter set for updating Gaussian componentApproximation;wherein,andrespectively according to the effective measurement zk,efAnd the weight, mean, covariance matrix, and label of the updated gaussian component obtained from the jth predicted component,λ is the average number of clutter obtained by the radar at time k, c (z)k) Arbitrary measurements z obtained for radarkThe probability density of the spatial distribution of (a),Hkfor the measurement matrix, RkMeasuring a noise covariance matrix; i is an identity matrix and is a matrix of the identity,obtaining a multi-target posterior intensity D at an approximate k momentk|k(x),Is composed of Jk|k-1(1+|Zk,ef|) parameter sets for updating gaussian componentsIn the approximation that the difference between the first and second values,andthe weight, mean, covariance matrix and label of the j1 th updated gaussian component respectively;
step 5, multi-target state extraction
When k is 1, no Gaussian component is updated, and the multi-target state extraction step is not executed;
when k is larger than or equal to 2, executing a multi-target state extraction step, and extracting a plurality of single-target state estimations with identity identifications from the updated Gaussian components, wherein the method specifically comprises the following steps:
(5a) setting four query matrixes Uw、Um、UP、UlAre respectively Jk|k-1(1+|Zk,efI) matrix of dimension, J obtained in step 4k|k-1(1+|Zk,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in Uw、Um、UPAnd UlIn (1), wherein,are stored to U in sequencewThe first column of (a) is,are stored to U in sequencemThe first column of (a) is,are stored to U in sequencePThe first column of (a) is,are stored to U in sequencelThe first column of (1); press ze ═ 1., | Zk,efL, stored in sequence byUpdated parameters of Gaussian component to Uw、Um、UPAnd UlIn the specification:are stored to U in sequencewThe ze +1 th column of (a),are stored to U in sequencemThe ze +1 th column of (a),are stored to U in sequencePThe ze +1 th column of (a),are stored to U in sequencelZe +1 column (b);
(5d) Get UwMaximum value of all elements from the second column to the last columnIf this maximum value is presentLess than basic weight threshold w of state extraction1Directly jumping to the step (5 h); otherwise, according to the maximum valueAt a different position, from UlTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label setIs provided withThe system is used for storing a part of column number sets of effective measurement in the query matrix;
(5e) push buttonFor all labels in turn areThe processing of the gaussian component specifically includes:
(5e.1) the weight is givenThe label isIs at UeIn (1), the row number and column number are denoted as rn(i1)And cn(i1)Is recorded in UmElements in corresponding positions inIs m(i1),UeWherein all tags areIs recorded as
(5e.2) judging the tagAccording to the attribute, the corresponding operation is carried out on the attribute, and the following three conditions exist:
A. when in useThenA label for a new component, wherein Vnew>Vun2The method comprises the following specific operation steps:
firstly, whenEstablishing a two-level transient track, wherein w2For the second-level transient weight threshold, the specific operation is as follows:
rmax2=rmax2+1
wherein r ismax2=rmax2+1, updating the label of the secondary transient track; note the bookMiddle (rn)(i1)The elements of a row in all columns areWill be provided withAll elements in (1) are changed to rmax2(ii) a Note the bookMiddle rn(i1)The elements of a row in all columns areNote UwMiddle rn(i1)Rows and cn(i1)The elements of the column areFirstly, the method is carried outAll elements in (1) are set to zero and thenMiddle rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoRespectively counting the transient tracks rmax2Secondary transient set stored at time kAnd first order transient setPerforming the following steps;
rmax1=rmax1+1
wherein r ismax1=rmax1+1, updating the label of the primary transient track; will be provided withAll elements in (1) are changed to rmax1(ii) a Firstly, the method is carried outAll elements in (1) are set to zero and thenMiddle in rn(i1)Rows and cn(i1)Elements of a columnInstead, it is changed intoThe number r of transient tracksmax1First order transient set stored at time kPerforming the following steps;
③ the number cn of the row(i1)Stored into CI, i.e. CI ═ CI cn(i1)](ii) a Changes are made toIs rn(i1)I.e. by
B. When in useThenThe method is a label of a transient component, and comprises the following specific steps:
firstly, whenThen, further on,for labels of the two types of transient components, the following steps are sequentially executed:
secondly, if k is not less than 3 andfirstly, a new determined track r is establishedmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)];
② whenThen, further on,for a class of transient classified tags, the following steps are performed in sequence:
Secondly, if k is more than or equal to 4,And isFirstly, a new determined track r is establishedmax=rmax+1, then get single target state estimateAnd modifying the corresponding parametersFinally updating CI ═ CI cn(i1)];
C. When in useThenTo determine the label of the component, from labelsThe method comprises the following specific steps of:
let w3Is a weight threshold value of state preprocessing and satisfies w1≤w3<w2;
Wherein,is at UwMiddle labelAll the rows correspond to the row number of the maximum value in all the elements of the first column; | | b1-b2| | represents the euclidean distance between b1 and b 2; note UmTo (1) aThe elements in row number 1 areKnown zero mean metrology noise covariance matrix Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]T,d(1)=||σ-0||;
Otherwise, executing the step II;
② will m(i1)As a single target state estimate, store it to a set of state estimatesIn, i.e.Will be provided withAnd k are stored separately to the track tag setAnd set of occurrence momentsIn, i.e.Andat the same time, adding cn(i1)Store to CI, i.e. CI ═ CI cn(i1)];
(iii) correction of weight of Gaussian component
Firstly, calculate m one by one(i1)Euclidean distance to the same tag component:
and correcting the weight according to the conditions:
wherein, UwAnd UmTo (1) aThe elements in the row at the ze1 column are respectively marked asAndattenuation factor a1 ═ 1/3;
(5e.3) based on two "one-to-one" principles, U is dividedwThe middle row has the number ofIs set to zero, i.e. is performed
(5f) Is executed completelyAfter circulation, the elements of CI are used as a column number set, and then the elements are listed in U according to two 'one-to-one' principleswAll the elements of any row of (1) are set to zero;
(5g) returning to the step (5 d);
(5h) respectively to be provided withAndj in (1)k|k-1(1+|Zk,efI) conversion of elements column by column to 1 × (J)k|k-1(1+|Zk,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
(5i) Connecting state estimates with the same identity at different moments to obtain an explicit multi-target track;
step 6, pruning with Gaussian components
(6a) Cut out a succession of nnostThe method comprises the following specific steps of:
(6a.1) takingIn non-repetitive tab sets, from which to deleteThe label set of the component of which the state estimation is not obtained at the moment k is obtained
(6a.2) pressingSequentially discriminatingIn thatk1=k-nnotrThe number of occurrences in k-1: if it is notAt each time k1If all appear in (1), then execute(6b) All the weights are larger than a small weight threshold value wetThe sequence number of the Gaussian component of (1) is recorded in I1, i.e.
(6c) Let Jk|k=|I1|,By analogy, obtainWherein,andthe weight, the mean, the covariance matrix and the label of the j3 th Gaussian component respectively;
step 7, obtaining parameter sets of Gaussian components based on the pruning in step 6Performing Gaussian component combination to makeIs a set of sequence numbers of the Gaussian components, andthe method comprises the following specific steps:
(7b) Taking andorder of co-tagged, combinable componentsThe set of numbers is I3,wherein d isThIs a merged distance threshold;
(7c) all gaussian component parameters in the sequence number set I3 are merged, obtaining a new Gaussian component;
(7d) i2 ═ I2-I3, and return to step (7 a);
Step 8, generating new target intensity
According to Zk,resiPress ze2 ═ 1., | Zk,resiObtaining the parameters of 2 th new components according to the following operations:
obtaining the prior intensity gamma of the new target at the k +1 momentk+1(x),Jγ,k+1=|Zk,resi|,i=ze2;γk+1(x) Is composed of Jγ,k+1Parameter set of gaussian componentsApproximation;
and 9, returning to the step 2 when k is equal to k + 1.
2. The explicit multi-target tracking method based on GM-PHD filter under adaptive new birth intensity as claimed in claim 1, wherein the specific steps in step 3a are:
(3a.1) known zero mean of the measured noise covariance matrix Andthe measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ]x σy]TBased on this, calculating a gate threshold d (a) ═ a | | | σ -0|, where a is a confidence coefficient;
(3a.3)For the latest measurement set obtained at time k, ze1 ═ 1kL, sequentially judging eachIs determined by the attribute ofThen, press J ═ 1k|k-1One by one calculateTo each gaussian component meanIs a distance ofIf d is(ze1,j)D (a) or less, then judgingNearby Gaussian component, seeFor effective measurement, add it to the effective measurement set Zk,efThen jump out of the loop and execute the pairAnd (4) judging.
3. The explicit multi-target tracking method based on GM-PHD filter under adaptive newness intensity of claim 1, where V in step 1un<Vun2。
4. The explicit multi-target tracking method under adaptive newness strength based on GM-PHD filter of claim 3 where V isun=500,Vun2=1000。
5. The explicit multi-target tracking method based on GM-PHD filter under adaptive newness strength according to claim 1, wherein step 7 further comprises step (7f), specifically:
the maximum Gaussian component number of the iteration to the k +1 moment does not exceed Jmax1When is coming into contact withWhen it comes toBefore taking Jmax1The Gaussian component with the maximum weight is recorded as When in use, willRelabeling toThereby obtaining the posterior intensity of survival multiple targets
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910949144.4A CN110780269B (en) | 2019-10-08 | 2019-10-08 | Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910949144.4A CN110780269B (en) | 2019-10-08 | 2019-10-08 | Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110780269A CN110780269A (en) | 2020-02-11 |
CN110780269B true CN110780269B (en) | 2022-03-25 |
Family
ID=69385418
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910949144.4A Active CN110780269B (en) | 2019-10-08 | 2019-10-08 | Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110780269B (en) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111562571B (en) * | 2020-05-28 | 2022-04-29 | 江南大学 | Maneuvering multi-target tracking and track maintaining method for unknown new-born strength |
CN111856442A (en) * | 2020-07-03 | 2020-10-30 | 哈尔滨工程大学 | Multi-target tracking method for self-adaptively estimating strength of newborn target based on measured value driving |
CN111811515B (en) * | 2020-07-03 | 2022-10-04 | 浙江工业大学 | Multi-target track extraction method based on Gaussian mixture probability hypothesis density filter |
CN112328959B (en) * | 2020-10-14 | 2023-04-11 | 哈尔滨工程大学 | Multi-target tracking method based on adaptive extended Kalman probability hypothesis density filter |
CN112688667B (en) * | 2020-12-22 | 2022-10-11 | 中国人民解放军63921部队 | Design method of GM-PHD filter |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH1138130A (en) * | 1997-07-23 | 1999-02-12 | Tech Res & Dev Inst Of Japan Def Agency | Multi-target tracking device |
CN105761276A (en) * | 2015-12-15 | 2016-07-13 | 江南大学 | Iteration RANSAC-based adaptive birth target intensity estimation GM-PHD multi-target tracking algorithm |
CN105844217A (en) * | 2016-03-11 | 2016-08-10 | 南京航空航天大学 | Multi-target tracking method based on measure-driven target birth intensity PHD (MDTBI-PHD) |
CN107526070A (en) * | 2017-10-18 | 2017-12-29 | 中国航空无线电电子研究所 | The multipath fusion multiple target tracking algorithm of sky-wave OTH radar |
CN107797106A (en) * | 2017-05-08 | 2018-03-13 | 南京航空航天大学 | A kind of PHD multiple target tracking smooth filtering methods of the unknown clutter estimations of acceleration EM |
-
2019
- 2019-10-08 CN CN201910949144.4A patent/CN110780269B/en active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH1138130A (en) * | 1997-07-23 | 1999-02-12 | Tech Res & Dev Inst Of Japan Def Agency | Multi-target tracking device |
CN105761276A (en) * | 2015-12-15 | 2016-07-13 | 江南大学 | Iteration RANSAC-based adaptive birth target intensity estimation GM-PHD multi-target tracking algorithm |
CN105844217A (en) * | 2016-03-11 | 2016-08-10 | 南京航空航天大学 | Multi-target tracking method based on measure-driven target birth intensity PHD (MDTBI-PHD) |
CN107797106A (en) * | 2017-05-08 | 2018-03-13 | 南京航空航天大学 | A kind of PHD multiple target tracking smooth filtering methods of the unknown clutter estimations of acceleration EM |
CN107526070A (en) * | 2017-10-18 | 2017-12-29 | 中国航空无线电电子研究所 | The multipath fusion multiple target tracking algorithm of sky-wave OTH radar |
Non-Patent Citations (1)
Title |
---|
高斯混合概率假设密度滤波器在多目标跟踪中的应用;吕学斌 等;《计算机学报》;20120229;第35卷(第2期);第397-404页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110780269A (en) | 2020-02-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110780269B (en) | Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity | |
CN107424171B (en) | Block-based anti-occlusion target tracking method | |
CN106407677B (en) | A kind of multi-object tracking method in the case of missing measurement | |
CN110967690B (en) | Multi-target tracking method based on multiple Bernoulli distributed multiple sensors | |
Xia et al. | Extended target Poisson multi-Bernoulli mixture trackers based on sets of trajectories | |
CN104133983B (en) | Multiple-target Bernoulli distribution type fusion method based on weighting negative first-order RD (Renyi Divergence) | |
CN110376581B (en) | Explicit multi-target tracking method based on Gaussian mixture probability hypothesis density filter | |
CN105182291A (en) | Multi-target tracking method for PHD smoother adaptive to target nascent strength | |
CN106934324A (en) | Based on the radar data correlating methods for simplifying many hypothesis algorithms | |
CN105761276A (en) | Iteration RANSAC-based adaptive birth target intensity estimation GM-PHD multi-target tracking algorithm | |
CN105975772B (en) | Tracking before multi-target detection based on probability hypothesis density filtering | |
CN111562569A (en) | Weighted group sparse constraint-based multi-target constant false alarm detection method under Weibull background | |
CN108717702B (en) | Probabilistic hypothesis density filtering smoothing method based on segmented RTS | |
CN115204212A (en) | Multi-target tracking method based on STM-PMBM filtering algorithm | |
CN113296089B (en) | LMB density fusion method and device for multi-early-warning-machine target tracking system | |
CN108153519B (en) | Universal design framework for target intelligent tracking method | |
CN116299195A (en) | Radar signal processing method based on TOA sequence relativity | |
CN115657008A (en) | Multi-target tracking method and device for airborne terahertz radar | |
CN104850856A (en) | Multi-extension target tracking method for affinity propagation cluster observation | |
CN115097437A (en) | Underwater target tracking track approaching and crossing solution method based on label multi-Bernoulli pre-detection tracking algorithm | |
CN114037017A (en) | Data fusion method based on error distribution fitting | |
Bai et al. | RCS Information Aided PMBM Filter | |
Yu et al. | A hypergraph matching labeled multi-Bernoulli filter for group targets tracking | |
CN110880012A (en) | Correlation method for frequency information of agile radar radiation source between multiple scout platforms | |
CN111811515A (en) | Multi-target track extraction method based on Gaussian mixture probability hypothesis density filter |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |