CN110780269A - 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 PDF

Info

Publication number
CN110780269A
CN110780269A CN201910949144.4A CN201910949144A CN110780269A CN 110780269 A CN110780269 A CN 110780269A CN 201910949144 A CN201910949144 A CN 201910949144A CN 110780269 A CN110780269 A CN 110780269A
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.)
Granted
Application number
CN201910949144.4A
Other languages
Chinese (zh)
Other versions
CN110780269B (en
Inventor
高乙月
蒋德富
蒋康辉
韩艳
王松
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN201910949144.4A priority Critical patent/CN110780269B/en
Publication of CN110780269A publication Critical patent/CN110780269A/en
Application granted granted Critical
Publication of CN110780269B publication Critical patent/CN110780269B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

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

Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity
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:
step 1, initializing the number of various tracks
When k is 0, determining the label r of the Gaussian component maxEach new determined track is built by taking 0 as a base number, r maxAdding 1; the label of the second order Gaussian component is given by r max2＝V un2As the base number, r is added every time a second-level transient track is created max2 Adding 1; labeling of first order transient Gaussian components with r max1＝V unAs a base number, r is a transient state track every time a new one is built max1 Adding 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 moment k-1|k-1(x) I.e. a priori information about surviving objects is known at time k
Figure BDA0002225001040000021
Is composed of J k-1Parameter set of gaussian components
Figure BDA0002225001040000022
And (ii) an approximation, wherein,
Figure BDA0002225001040000023
representing a Gaussian density, x being a single targetX ═ x px vy py v] T，[x py p] TIndicates the position of a single target, [ x ] vy v] TRepresenting the speed of a single target; and
Figure BDA0002225001040000032
the 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 k S,k|k-1(x)，
Figure BDA0002225001040000033
Is composed of J k-1Parameter set of gaussian components
Figure BDA0002225001040000034
And (ii) an approximation, wherein,
Figure BDA0002225001040000035
and
Figure BDA0002225001040000036
the weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k,
Figure BDA0002225001040000038
wherein p is S,kProbability of survival of a single target, F k-1And Q k-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 moment k(x)，
Figure BDA0002225001040000039
Is composed of J γ,kParameter set of gaussian components
Figure BDA00022250010400000310
Approximation of, wherein
Figure BDA00022250010400000311
And
Figure BDA00022250010400000312
respectively 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)，
Figure BDA00022250010400000313
Is composed of J γ,kParameter set of gaussian components
Figure BDA00022250010400000330
And (ii) an approximation, wherein,
Figure BDA00022250010400000315
and respectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
Figure BDA00022250010400000317
Figure BDA00022250010400000318
merging And obtaining a Gaussian component parameter set
Figure BDA00022250010400000321
J k|k-1＝J k-1+J γ，k，
Figure BDA00022250010400000323
And respectively 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 k k|k-1(x)，
Figure BDA00022250010400000325
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radar
Figure BDA00022250010400000326
And the prediction Gaussian component parameter set obtained in step 2
Figure BDA00022250010400000327
Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtained
Figure BDA00022250010400000328
Wherein the content of the first and second substances,
Figure BDA00022250010400000329
is the xy-plane ze1 measurements obtained by radar,
Figure BDA0002225001040000041
Figure BDA0002225001040000042
is the ze number effective measurement obtained by adopting a wave gate method, and | is the aggregate potential;
(3b) execution of Z k,resi＝Z k-Z k,efObtaining a residual measurement set
Figure BDA0002225001040000043
When k is 1, no gaussian component parameter set is predicted, then Z k,efIs empty, Z k,resi＝Z k；
Step 4, updating Gaussian components
When k is 1, there is no prediction gaussian component parameter set, Z k,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 Z k,efAnd updating the multiple posterior intensities D at the time k k|k(x) Wherein
Figure BDA0002225001040000045
p D,kThe detection probability of a single target at the moment k is (1-p) D,k)D k|k-1(x) Is composed of J k|k-1Parameter set for updating Gaussian component
Figure BDA00022250010400000419
Approximation; wherein the content of the first and second substances,
Figure BDA0002225001040000048
and
Figure BDA0002225001040000049
respectively according to the effective measurement z k,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 radar kThe probability density of the spatial distribution of (a),
Figure BDA00022250010400000411
H kfor the measurement matrix, R kMeasuring a noise covariance matrix;
Figure BDA00022250010400000412
i is an identity matrix and is a matrix of the identity,
Figure BDA00022250010400000420
obtaining a multi-target posterior intensity D at an approximate k moment k|k(x)，
Figure BDA00022250010400000415
Is composed of J k|k-1(1+|Z k,ef|) parameter sets for updating gaussian components In the approximation that the difference between the first and second values, and
Figure BDA00022250010400000418
the 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 U w、U m、U P、U lAre respectively J k|k-1(1+|Z k,efI) matrix of dimension, J obtained in step 4 k|k-1(1+|Z k,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in U w、U m、U PAnd U lIn (1), wherein,
Figure BDA0002225001040000051
are stored to U in sequence wThe first column of (a) is,
Figure BDA0002225001040000052
are stored to U in sequence mThe first column of (a) is, are stored to U in sequence PThe first column of (a) is,
Figure BDA0002225001040000054
are stored to U in sequence lThe first column of (1); press ze ═ 1., | Z k,efL, stored in sequence by
Figure BDA0002225001040000055
Updated parameters of Gaussian component to U w、U m、U PAnd U lIn the specification:
Figure BDA0002225001040000056
are stored to U in sequence wThe ze +1 th column of (a), are stored to U in sequence mThe ze +1 th column of (a),
Figure BDA0002225001040000058
are stored to U in sequence PThe ze +1 th column of (a),
Figure BDA0002225001040000059
are stored to U in sequence lZe +1 column (b);
(5b) order to
Figure BDA00022250010400000510
(5c) Let state estimate set at time k
Figure BDA00022250010400000511
Track label set Set of occurrence moments
Figure BDA00022250010400000513
(5d) Get U wMaximum value of all elements from the second column to the last column If this maximum value is present
Figure BDA00022250010400000515
Less than basic weight threshold w of state extraction 1Directly jumping to the step (5 h); otherwise, according to the maximum value
Figure BDA00022250010400000516
At a different position, from U lTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label set
Figure BDA00022250010400000517
Is provided with The system is used for storing a part of column number sets of effective measurement in the query matrix;
(5e) push button
Figure BDA00022250010400000519
For all labels in turn are
Figure BDA00022250010400000520
The processing of the gaussian component specifically includes:
(5e.1) the weight is given
Figure BDA00022250010400000521
The label is Is at U lIn (1), the row number and column number are denoted as rn (i1)And cn (i1)Is recorded in U mElements in corresponding positions in
Figure BDA00022250010400000523
Is m (i1)，U lWherein all tags are
Figure BDA00022250010400000524
Is recorded as
Figure BDA00022250010400000525
(5e.2) judging the tag
Figure BDA00022250010400000526
According to the attribute, the corresponding operation is carried out on the attribute, and the following three conditions exist:
A. when in use Then
Figure BDA00022250010400000528
A label for a new component, wherein V new＞V un2The method comprises the following specific operation steps:
① when
Figure BDA00022250010400000529
Establishing a two-level transient track, wherein w 2For the second-level transient weight threshold, the specific operation is as follows:
r max2＝r max2+1
Figure BDA00022250010400000531
Figure BDA00022250010400000533
wherein r is max2＝r max2+1, updating the label of the secondary transient track; note the book
Figure BDA00022250010400000534
Middle (rn) (i1)The elements of a row in all columns are Will be provided with
Figure BDA0002225001040000062
All elements in (1) are changed to r max2(ii) a Note the book Middle rn (i1)The elements of a row in all columns are Note U wMiddle rn (i1)Rows and cn (i1)The elements of the column are
Figure BDA0002225001040000065
Firstly, the method is carried out
Figure BDA0002225001040000066
All elements in (1) are set to zero and then Middle rn (i1)Rows and cn (i1)Elements of a column
Figure BDA0002225001040000068
Instead, it is changed into
Figure BDA0002225001040000069
Respectively counting the transient tracks r max2Secondary transient set stored at time k And first order transient set
Figure BDA00022250010400000611
Performing the following steps;
② when
Figure BDA00022250010400000612
Establishing a first-level transient track, and specifically operating as follows:
r max1＝r max1+1
Figure BDA00022250010400000613
Figure BDA00022250010400000615
wherein r is max1＝r max1+1, updating the label of the primary transient track; will be provided with
Figure BDA00022250010400000616
All elements in (1) are changed to r max1(ii) a Firstly, the method is carried out
Figure BDA00022250010400000617
All elements in (1) are set to zero and then
Figure BDA00022250010400000618
Middle rn (i1)Rows and columnscn (i1)Elements of a column
Figure BDA00022250010400000619
Instead, it is changed into
Figure BDA00022250010400000620
The number r of transient tracks max1First order transient set stored at time k
Figure BDA00022250010400000621
Performing the following steps;
③ column number cn (i1)Stored into CI, i.e. CI ═ CI cn (i1)](ii) a Changes are made to
Figure BDA00022250010400000622
Is rn (i1)I.e. by
Figure BDA00022250010400000623
B. When in use Then
Figure BDA00022250010400000625
The method is a label of a transient component, and comprises the following specific steps:
① when
Figure BDA00022250010400000626
Then, further on,
Figure BDA00022250010400000627
for labels of the two types of transient components, the following steps are sequentially executed:
first, it executes
Figure BDA00022250010400000628
And
Figure BDA00022250010400000629
secondly, if k is not less than 3 and
Figure BDA00022250010400000630
firstly, a new determined track r is established max＝r max+1, then get single target state estimate
Figure BDA00022250010400000631
And modifying the corresponding parameters
Figure BDA00022250010400000632
Finally updating CI ═ CI cn (i1)]；
② when Then, further on,
Figure BDA00022250010400000634
for a class of transient classified tags, the following steps are performed in sequence:
first, it executes
Figure BDA00022250010400000635
Secondly, if k is more than or equal to 4, And is
Figure BDA00022250010400000637
Firstly, a new determined track r is established max＝r max+1, then get single target state estimate
Figure BDA00022250010400000638
And modifying the corresponding parameters
Figure BDA00022250010400000639
Finally updating CI ═ CI cn (i1)]；
C. When in use
Figure BDA0002225001040000071
Then
Figure BDA0002225001040000072
To determine the label of the component, from labels The method comprises the following specific steps of:
① order w 3Is a weight threshold value of state preprocessing and satisfies w 1≤w 3＜w 2；
If it is not
Figure BDA0002225001040000074
The following state estimation preprocessing steps are performed in sequence:
Figure BDA0002225001040000075
Figure BDA0002225001040000076
if it is not
Figure BDA0002225001040000077
Wherein the content of the first and second substances,
Figure BDA0002225001040000078
is at U wMiddle label
Figure BDA0002225001040000079
All 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 U mTo (1) a
Figure BDA00022250010400000710
The elements in row number 1 are Known zero mean metrology noise covariance matrix
Figure BDA00022250010400000712
Figure BDA00022250010400000713
And the measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ] xσ y] T，d(1)＝||σ-0||；
Otherwise, go to step ②;
② mixing m (i1)As a single target state estimate, store it to a set of state estimates
Figure BDA00022250010400000715
In, i.e.
Figure BDA00022250010400000716
Will be provided with
Figure BDA00022250010400000717
And k are stored separately to the track tag set
Figure BDA00022250010400000718
And set of occurrence moments
Figure BDA00022250010400000719
In, i.e.
Figure BDA00022250010400000720
And at the same time, adding cn (i1)Store to CI, i.e. CI ═ CI cn (i1)]；
③ weight correction of Gaussian component
Firstly, calculate m one by one (i1)Euclidean distance to the same tag component:
Figure BDA00022250010400000722
and correcting the weight according to the conditions:
Figure BDA00022250010400000723
and (d) (i2,ze1)＞d(3)),
Or
Figure BDA00022250010400000724
And (d) (i2,ze1)＞d(4)),
Figure BDA00022250010400000725
if or
Figure BDA00022250010400000726
And (d) (i2,ze1)＞d(5)),
Or
Figure BDA00022250010400000727
And (d) (i2,ze1)＞d(6))
Wherein, U wAnd U mTo (1) a The elements in the row at the ze1 column are respectively marked as
Figure BDA0002225001040000082
And
Figure BDA0002225001040000083
attenuation factor a1 ═ 1/3;
(5e.3) based on two "one-to-one" principles, U is divided wThe middle row has the number of
Figure BDA0002225001040000084
Is set to zero, i.e. is performed (5f) Is executed completely
Figure BDA0002225001040000086
After 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' principles wAll the elements of any row of (1) are set to zero;
(5g) returning to the step (5 d);
(5h) respectively to be provided with
Figure BDA0002225001040000087
And
Figure BDA0002225001040000088
j in (1) k|k-1(1+|Z k,efI) conversion of elements column by column to 1 × (J) k|k-1(1+|Z k,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
Figure BDA0002225001040000089
(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 n nostThe method comprises the following specific steps of:
(6a.1) taking
Figure BDA00022250010400000810
In non-repetitive tab sets, from which to delete The label set of the component of which the state estimation is not obtained at the moment k is obtained
Figure BDA00022250010400000812
(6a.2) pressing
Figure BDA00022250010400000813
Sequentially discriminating In that
Figure BDA00022250010400000815
k1＝k-n notrThe number of occurrences in k-1: if it is not
Figure BDA00022250010400000816
At each time k1
Figure BDA00022250010400000817
If all appear in (1), then execute
Figure BDA00022250010400000818
(6b) All the weights are larger than a small weight threshold value w etThe sequence number of the Gaussian component of (1) is recorded in I1, i.e.
Figure BDA00022250010400000819
(6c) Let J k|k＝|I1|，
Figure BDA00022250010400000820
By analogy, obtain
Figure BDA00022250010400000821
Wherein the content of the first and second substances,
Figure BDA00022250010400000822
and
Figure BDA00022250010400000823
the 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 6
Figure BDA00022250010400000824
Performing Gaussian component combination to make
Figure BDA00022250010400000825
Is a set of sequence numbers of the Gaussian components, and the method comprises the following specific steps:
(7a) when I2 is empty, go to step (7e), otherwise,
Figure BDA00022250010400000827
obtaining a label
Figure BDA00022250010400000828
(7b) Taking and
Figure BDA00022250010400000829
the sequence number set of the co-tagged and combinable components is I3, wherein d is ThIs a merged distance threshold;
(7c) all gaussian component parameters in the sequence number set I3 are merged,
Figure BDA0002225001040000093
Figure BDA0002225001040000094
obtaining a new Gaussian component;
(7d) i2 ═ I2-I3, and return to step (7 a);
(7e) obtaining a combined Gaussian component set
Figure BDA0002225001040000095
Step 8, generating new target intensity
According to Z k,resiPress ze2 ═ 1., | Z k,resiL in turn according toThe following operations yield the parameters for the ze2 th new component:
Figure BDA0002225001040000096
Figure BDA0002225001040000097
Figure BDA0002225001040000098
Figure BDA0002225001040000099
obtaining the prior intensity gamma of the new target at the k +1 moment k+1(x)，
Figure BDA00022250010400000910
J γ,k+1＝|Z k,resi|，i＝ze2；γ k+1(x) Is composed of J γ,k+1Parameter set of gaussian components
Figure BDA00022250010400000911
Approximation;
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
Figure BDA00022250010400000912
And
Figure BDA00022250010400000914
the 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.2) setting up effective measurement set
Figure BDA00022250010400000915
(3a.3)
Figure BDA00022250010400000916
For the latest measurement set obtained at time k, ze1 ═ 1 kL, sequentially judging each
Figure BDA00022250010400000917
Is determined by the attribute of
Figure BDA00022250010400000918
Then, press J ═ 1 k|k-1One by one calculate
Figure BDA00022250010400000919
To each gaussian component mean Is a distance of
Figure BDA00022250010400000921
If d is (ze1,j)D (a) or less, then judging
Figure BDA00022250010400000922
Nearby Gaussian component, see
Figure BDA00022250010400000923
For effective measurement, add it to the effective measurement set Z k,efThen jump out of the loop and execute the pair
Figure BDA0002225001040000101
And (4) judging.
Further, step 1, V un＜V un2。
Further, V un＝500，V un2＝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 J max1When is coming into contact with
Figure BDA0002225001040000102
When it comes to Before taking J max1The Gaussian component with the maximum weight is recorded as
Figure BDA0002225001040000104
J k＝J max1； When in use, will
Figure BDA0002225001040000106
Relabeling to
Figure BDA0002225001040000107
Thereby obtaining the posterior intensity of survival multiple targets
Figure BDA0002225001040000108
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:
step 1, initializing the number of various tracks
When k is 0, determining the label r of the Gaussian component maxEach new determined track is built by taking 0 as a base number, r maxAdding 1; secondary transient gaussian component with r max2＝V un2As the base number, r is added every time a second-level transient track is created max2Adding 1; first order transient gaussian componentIs expressed by r max1＝V unAs a base number, r is a transient state track every time a new one is built max1Adding 1; wherein, V un2And V unCan take different values in different scenes, and V un＜V un2In this patent, V is specified un＝500，V un2＝1000。
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 (4) survival multi-target information prediction. The posterior intensity D of the surviving multiple targets is obtained at the known k-1 moment k-1|k-1(x) I.e. a priori information about surviving objects is known at time k Is composed of J k-1Parameter set of gaussian components
Figure BDA0002225001040000122
And (ii) an approximation, wherein,
Figure BDA0002225001040000123
denotes the Gaussian density, x is the state of a single target, x ═ x px vy py v] T，[x py p] TIndicates the position of a single target, [ x ] vy v] TRepresenting the speed of a single target;
Figure BDA0002225001040000124
Figure BDA0002225001040000125
and
Figure BDA0002225001040000126
the 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 k S,k|k-1(x)，
Figure BDA0002225001040000127
Is composed of J k-1Parameter set of gaussian components
Figure BDA0002225001040000128
And (ii) an approximation, wherein,
Figure BDA0002225001040000129
and
Figure BDA00022250010400001210
the weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k,
Figure BDA00022250010400001211
Figure BDA00022250010400001212
wherein p is S,kProbability of survival of a single target, F k-1And Q k-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 moment k(x)，
Figure BDA00022250010400001213
Is composed of J γ,kParameter set of gaussian components
Figure BDA00022250010400001214
Approximation of, wherein And
Figure BDA00022250010400001216
respectively 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)，
Figure BDA00022250010400001217
Is composed of J γ,kParameter set of gaussian components
Figure BDA00022250010400001218
And (ii) an approximation, wherein,
Figure BDA00022250010400001219
and
Figure BDA00022250010400001220
respectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
Figure BDA00022250010400001221
Figure BDA00022250010400001222
merging
Figure BDA00022250010400001223
And
Figure BDA00022250010400001224
obtaining a Gaussian component parameter set J k|k-1＝J k-1+J γ,k，
Figure BDA00022250010400001226
And weight, mean, covariance matrix, and label for the jth prediction component, respectivelyThen obtaining the multi-target prediction strength D at the k moment k|k-1(x)，
Figure BDA0002225001040000131
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radar
Figure BDA0002225001040000132
And the prediction Gaussian component parameter set obtained in step 2
Figure BDA0002225001040000133
Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtained
Figure BDA0002225001040000134
Wherein the content of the first and second substances,
Figure BDA0002225001040000135
is the xy-plane ze1 measurements obtained by radar,
Figure BDA0002225001040000136
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
Figure BDA0002225001040000137
Figure BDA0002225001040000138
And
Figure BDA0002225001040000139
the 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.2) setting up effective measurement set
Figure BDA00022250010400001310
(3a.3)
Figure BDA00022250010400001311
For the latest measurement set obtained at time k, ze1 ═ 1 kL, sequentially judging each
Figure BDA00022250010400001312
The attribute of (2). Judgment of
Figure BDA00022250010400001313
Then, press J ═ 1 k|k-1One by one calculate
Figure BDA00022250010400001314
To each gaussian component mean
Figure BDA00022250010400001315
Is a distance of If d is (ze1,j)D (a) or less, then judging Nearby Gaussian component, see
Figure BDA00022250010400001318
For effective measurement, add it to the effective measurement set Z k,efThen jump out of the loop and execute the pair
Figure BDA00022250010400001319
And (4) judging.
(3b) Execution of Z k,resi＝Z k-Z k,efObtaining a residual measurement set
Figure BDA00022250010400001320
When k is 1, no gaussian component parameter set is predicted, then Z k,efIs empty，Z k,resi＝Z k。
Step 4, updating Gaussian components
According to the obtained effective measurement set Z k,efAnd a set of Gaussian components for the prediction of time k
Figure BDA00022250010400001321
Updating the multiple posterior intensities D at the time k k|k(x) In that respect Known as p D,kThe detection probability of a single target at the time k is specifically as follows:
Figure BDA00022250010400001322
wherein the content of the first and second substances,
Figure BDA00022250010400001323
Figure BDA0002225001040000141
Figure BDA0002225001040000142
Figure BDA0002225001040000143
Figure BDA0002225001040000144
Figure BDA0002225001040000145
Figure BDA0002225001040000146
wherein (1-p) D,k)D k|k-1(x) From a set of Gaussian component parameters
Figure BDA0002225001040000147
And (4) approximation. D D,k(x；z k,ef) Is composed of J k|k-1Parameter set of gaussian components In the approximation that the difference between the first and second values,
Figure BDA0002225001040000149
and
Figure BDA00022250010400001410
respectively according to the effective measurement z k,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 radar kProbability density of spatial distribution of (1), H kFor the measurement matrix, R kTo measure the noise covariance matrix, I is the identity matrix.
Accordingly, the multi-target posterior intensity at the k time is determined by
Figure BDA00022250010400001411
Parameter set of gaussian components And (4) approximation. Wherein the content of the first and second substances,
Figure BDA00022250010400001413
and the 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 U w、U m、U P、U lIs J k|k-1(1+|Z k,efI) matrix of dimension, J obtained in step 4 k|k-1(1+|Z k,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in U w、U m、U PAnd U lIn (1). FIG. 2 shows U w、U m、U PAnd U lThe structure of the matrix. Wherein the content of the first and second substances,
Figure BDA00022250010400001415
are stored to U in sequence wThe first column of (a) is,
Figure BDA00022250010400001416
are stored to U in sequence mThe first column of (a) is,
Figure BDA00022250010400001417
are stored to U in sequence PThe first column of (a) is,
Figure BDA00022250010400001418
are stored to U in sequence lThe first column of (2). Press ze ═ 1., | Z k,efL, stored in sequence by
Figure BDA00022250010400001419
Updated parameters of Gaussian component to U w、U m、U PAnd U lIn the specification:
Figure BDA00022250010400001420
are stored to U in sequence wThe ze +1 th column of (a), are stored to U in sequence mThe ze +1 th column of (a),
Figure BDA0002225001040000152
are stored to U in sequence PThe ze +1 th column of (a), are stored to U in sequence lZe +1 column (b). Obviously, U PThe elements in each row are the same, U lThe elements of each row in the series are the same.
Storing the updated Gaussian component to U w、U m、U PAnd U lIn this way, the parent component of each update component and the corresponding measurement are conveniently searched.
(5b) Creating four temporary matrices
Figure BDA0002225001040000154
And parameters for recording various types of modifications in the state extraction process:
Figure BDA0002225001040000156
(5c) set of state estimates for time k
Figure BDA0002225001040000157
Track label set
Figure BDA0002225001040000158
Set of occurrence moments Effectively measured at U wColumn number set in (1)
(5d) Memory matrix U wThe elements of all rows in the second to last columns are
Figure BDA00022250010400001511
Taking the maximum value from the above, and recording the maximum value as
Figure BDA00022250010400001512
U wThere may be more than one
Figure BDA00022250010400001513
When in use
Figure BDA00022250010400001514
And (5h) jumping to the step. Otherwise, when
Figure BDA00022250010400001515
According to
Figure BDA00022250010400001516
At U wRow and column numbers at different positions in the slave U lTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label set
Figure BDA00022250010400001517
Wherein, w 1Taking w as a basic weight threshold value extracted for the state according to a plurality of experimental simulation results 1＝0.02。
(5e) Push button
Figure BDA00022250010400001518
For all labels in turn are
Figure BDA00022250010400001529
The processing of the gaussian component specifically includes:
(5e.1) the weight is given
Figure BDA00022250010400001520
The label is
Figure BDA00022250010400001521
Is at U lIn (1), the row number and column number are denoted as rn (i1)And cn (i1)Is recorded in U mElements in corresponding positions in
Figure BDA00022250010400001522
Is m (i1)，U lWherein all tags are
Figure BDA00022250010400001523
Is recorded as
(5e.2) judging the tag
Figure BDA00022250010400001525
According to the attribute, the corresponding operation is carried out on the attribute. There are three cases:
A. when in use
Figure BDA00022250010400001526
Then
Figure BDA00022250010400001527
Is a label for the nascent component. V newThe 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 V new＞V un2In this patent, V is specified through multiple simulation experiments new2000. The method comprises the following specific operation steps:
① when
Figure BDA00022250010400001528
And establishing a secondary transient state track. Wherein, w 2Is a secondary transient weight threshold, and is designated as w in the patent according to the results of multiple simulation experiments 20.5. The specific operation is as follows:
r max2＝r max2+1
Figure BDA0002225001040000161
Figure BDA0002225001040000162
wherein r is max2＝r max2+1, updating the label of the secondary transient track; note the book
Figure BDA0002225001040000165
Middle (rn) (i1)The elements of a row in all columns are
Figure BDA0002225001040000166
Will be provided with
Figure BDA0002225001040000167
All elements in (1) are changed to r max2(ii) a Note the book
Figure BDA0002225001040000168
Middle rn (i1)The elements of a row in all columns are
Figure BDA0002225001040000169
Note U wMiddle rn (i1)Rows and cn (i1)The elements of the column are
Figure BDA00022250010400001610
Firstly, the method is carried out
Figure BDA00022250010400001638
All elements in (1) are set to zero and then
Figure BDA00022250010400001639
Middle rn (i1)Rows and cn (i1)Elements of a column
Figure BDA00022250010400001637
Instead, it is changed into
Figure BDA00022250010400001614
Respectively counting the transient tracks r max2Secondary transient set stored at time k
Figure BDA00022250010400001615
And first order transient set In (1).
② when And establishing a first-level transient track. The specific operation is as follows:
r max1＝r max1+1
Figure BDA00022250010400001618
Figure BDA00022250010400001619
Figure BDA00022250010400001620
wherein r is max1＝r max1+1, updating the label of the primary transient track; will be provided with
Figure BDA00022250010400001621
All elements in (1) are changed to r max1(ii) a Firstly, the method is carried out
Figure BDA00022250010400001622
All elements in (1) are set to zero and then
Figure BDA00022250010400001623
Middle rn (i1)Rows and cn (i1)Elements of a column
Figure BDA00022250010400001624
Instead, it is changed into
Figure BDA00022250010400001625
The number r of transient tracks max1First order transient set stored at time k
Figure BDA00022250010400001626
In (1).
③ column number cn (i1)Stored into CI, i.e. CI ═ CI cn (i1)](ii) a Changes are made to
Figure BDA00022250010400001627
Is rn (i1)I.e. by
Figure BDA00022250010400001628
B. When in use
Figure BDA00022250010400001629
Then
Figure BDA00022250010400001630
The method is a label of a transient component, and comprises the following specific steps:
① when Then, further on,
Figure BDA00022250010400001632
labels for transient components of two classes.
First, it executes
Figure BDA00022250010400001633
And
Figure BDA00022250010400001634
secondly, if k is not less than 3 and
Figure BDA00022250010400001635
firstly, a new determined track r is established max＝r max+1, then get single target state estimate And modifying the corresponding parameters
Figure BDA0002225001040000171
Finally updating CI ═ CI cn (i1)]。
② when
Figure BDA0002225001040000172
Then, further on,
Figure BDA0002225001040000173
is a class of temporally classified labels.
First, it executes
Figure BDA0002225001040000174
Secondly, if k is more than or equal to 4,
Figure BDA0002225001040000175
And is
Figure BDA0002225001040000176
Firstly, a new determined track r is established max＝r max+1, then get single target state estimate
Figure BDA0002225001040000177
And modifying the corresponding parameters
Figure BDA0002225001040000178
Finally updating CI ═ CI cn (i1)]。
C. When in use
Figure BDA0002225001040000179
Then
Figure BDA00022250010400001710
To determine the label of the component. From the label is
Figure BDA00022250010400001711
Extracting the state estimation of the single target from the components, and the specific stepsThe method comprises the following steps:
① to reduce measurement noise, detection uncertainty, and clutter interference, it is necessary to perform state estimation preprocessing 3For the weight threshold of state preprocessing, it should satisfy w 1≤w 3＜w 2In this patent, w is specified through multiple simulation experiments 3＝0.2。
If it is not
Figure BDA00022250010400001730
The following state estimation preprocessing steps are performed in sequence:
Figure BDA00022250010400001713
if it is not
Figure BDA00022250010400001714
Wherein the content of the first and second substances,
Figure BDA00022250010400001715
is at U wMiddle label
Figure BDA00022250010400001716
All 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 U mTo (1) a
Figure BDA00022250010400001717
The elements in row number 1 are Known zero mean metrology noise covariance matrix
Figure BDA00022250010400001719
Figure BDA00022250010400001720
And
Figure BDA00022250010400001721
the measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ] xσ y] T，d(1)＝||σ-0||。
If it is not
Figure BDA00022250010400001722
The following steps are performed directly.
② mixing m (i1)As a single target state estimate, store it to a set of state estimates
Figure BDA00022250010400001723
In, i.e.
Figure BDA00022250010400001724
Will be provided with
Figure BDA00022250010400001725
And k are stored separately to the track tag set
Figure BDA00022250010400001726
And set of occurrence moments
Figure BDA00022250010400001727
In, i.e.
Figure BDA00022250010400001728
And at the same time, adding cn (i1)Store to CI, i.e. CI ═ CI cn (i1)]。
③ Gaussian component weight correction, the larger the value of w, the more concentrated the posterior distribution of the single target and the smaller the effective area radius, and the smaller the value, the more dispersed the posterior distribution of the single target and the larger the effective area radius.
For this purpose, m is first calculated one by one (i1)Distance d to the same label component (i2,ze1)，
Figure BDA0002225001040000181
Then, the weight value is corrected according to the condition,
Figure BDA0002225001040000182
and (d) (i2,ze1)＞d(3)),
Or
Figure BDA0002225001040000183
And (d) (i2,ze1)＞d(4)),
Figure BDA0002225001040000184
if or
Figure BDA0002225001040000185
And (d) (i2,ze1)＞d(5)),
Or
Figure BDA00022250010400001819
And (d) (i2,ze1)＞d(6)).
Wherein, U wAnd U mTo (1) a
Figure BDA00022250010400001820
The elements in the row at the ze1 column are respectively marked as
Figure BDA0002225001040000188
And
Figure BDA0002225001040000189
in 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 extraction wThe middle row has the number of
Figure BDA00022250010400001810
Is set to zero, i.e. is performed
(5f) Is executed completely After 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 extraction wAll of the elements of any row of (c) are zeroed.
(5g) And (5) returning to the step (5 d).
(5h) Respectively to be provided with
Figure BDA00022250010400001813
And
Figure BDA00022250010400001814
j in (1) k|k-1(1+|Z k,efI) conversion of elements column by column to 1 × (J) k|k-1(1+|Z k,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
Figure BDA00022250010400001815
(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 n nostThe next component without state estimation, n is specified in this patent nost6. The method comprises the following specific steps:
(6a.1) taking
Figure BDA00022250010400001816
In non-repetitive tab sets, from which to delete
Figure BDA00022250010400001817
The tag contained in, i.e. executing
Figure BDA00022250010400001818
Tag set for components for which state estimates are not obtained at time k is obtained
Figure BDA0002225001040000191
(6a.2) pressing
Figure BDA0002225001040000192
Sequentially discriminating
Figure BDA0002225001040000193
In that
Figure BDA0002225001040000194
k1＝k-n notrThe number of occurrences in k-1: if it is not
Figure BDA0002225001040000195
At each time k1 If all appear in (1), then execute
Figure BDA0002225001040000197
(6b) All the weights are larger than a small weight threshold value w etThe number of the gaussian component of (a) is recorded in I1, where w is specified in this patent et＝10 -5I.e. by
Figure BDA0002225001040000198
(6c) Let J k|k＝|I1|，
Figure BDA0002225001040000199
By analogy, obtain
Figure BDA00022250010400001910
Wherein the content of the first and second substances,
Figure BDA00022250010400001911
and
Figure BDA00022250010400001912
the weight, mean, covariance matrix, and label of the j3 th gaussian component, respectively.
Step 7, parameter set based on Gaussian component obtained after pruning And 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 to
Figure BDA00022250010400001914
Is a set of sequence numbers of the Gaussian components, and
Figure BDA00022250010400001915
the method comprises the following specific steps:
(7a) when I2 is empty, go to step (7e), otherwise,
Figure BDA00022250010400001916
obtaining a label
Figure BDA00022250010400001917
(7b) Taking and
Figure BDA00022250010400001918
the sequence number set of the components which are identical in label and can be combined is I3
Figure BDA00022250010400001919
Wherein d is ThFor the combined distance threshold, d is specified in this patent Th＝4。
(7c) All the gaussian component parameters in the sequence number set I3 are merged to obtain a new gaussian component:
Figure BDA00022250010400001920
Figure BDA00022250010400001921
Figure BDA00022250010400001922
Figure BDA00022250010400001923
(7d) i2 ═ I2-I3, and return to step (7 a).
(7e) Obtaining a combined Gaussian component set
(7f) The maximum Gaussian component number of the iteration to the k +1 moment does not exceed J max1In this patent, J is designated max1＝150。
When in use
Figure BDA0002225001040000202
When it comes to
Figure BDA0002225001040000203
Before taking J max1The Gaussian component with the maximum weight is recorded as
Figure BDA0002225001040000204
J k＝J max1；
Figure BDA0002225001040000205
When in use, will
Figure BDA0002225001040000206
Relabeling to
Figure BDA0002225001040000207
Thereby obtaining the posterior intensity of survival multiple targets
And 8, generating new target intensity. According to Z k,resiPress ze2 ═ 1., | Z k,resiObtaining the parameters of 2 th new components according to the following operations:
Figure BDA00022250010400002010
thus, the prior intensity γ of the new object at the time k +1 is obtained k+1(x)，
Figure BDA00022250010400002013
J γ,k+1＝|Z k,resi|，i＝ze2；γ k+1(x) Is composed of J γ,k+1Parameter set of gaussian components
Figure BDA00022250010400002014
And (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:
Figure BDA00022250010400002015
wherein x is k＝[x px vy py v] T，[x py p] TIs the position of the target at time k, [ x ] vy v] TIs the velocity of the target at time k. Δ ═ 1s for the sampling interval, σ ω＝5m/s 2. The object can appear or disappear at any time in the scene, and the survival probability p S,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 parameter
Figure BDA0002225001040000211
Wherein
Figure BDA0002225001040000212
Figure BDA0002225001040000213
Figure BDA0002225001040000214
P B＝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 is
Figure BDA0002225001040000215
Wherein the content of the first and second substances,
Figure BDA0002225001040000216
P γ,k＝P B，
Figure BDA0002225001040000217
i＝1,2,3,4。
the measurement equation of the target is
Wherein upsilon is x,kAnd upsilon y,kIs independent zero mean Gaussian white noise, and has mean square error of sigma x＝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 ═ λ/2000 2. 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 components 00.999. The maximum number of Gaussian components is J maxTruncate threshold w to 100 et＝10 -5Merging the thresholds d Th＝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
Simulation 1, clutter number λ 20, detection probability p D,kFig. 3 to 5 show the state estimation and the track of the method of the present invention, which are determined in a single experiment, when the value is 0.92. In fig. 3 to 5, the marks with the same gray scale and different marks belong to the same track. FIG. 6 is the average OSPA distance of 100 experiments according to the method of the present invention, and FIG. 7 is the corresponding calculationTime.
Simulation 2, clutter number λ 50, detection probability p D,kFig. 8 to 10 show the state estimation and the track of the method of the present invention, which are determined in a single experiment, when the value is 0.98. In fig. 8 to 10, the marks with the same gray scale and different marks belong to the same track. FIG. 11 is the average OSPA distance of the proposed method over 100 experiments, and FIG. 12 is the corresponding calculated time.
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, with the establishment of the determined track, the OSPA distance of the filter designed by the invention tends to be smooth. In FIG. 6, when the detection probability is relatively low, p D,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 during missing detection, it can be seen in fig. 4 and 5 that each determined track has a breakpoint, so that the OSPA distance of the invention is smaller than that of the new target in the scene known a prioriThe GLMB filter, i.e. the tracking performance, is inferior to it, however this does not affect the determination and tracking of the respective 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, p D,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 p is 50 D,kCorresponding detection probability p when λ is greater than 20 at 0.98 D,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 component maxEach new determined track is built by taking 0 as a base number, r maxAdding 1; the label of the second order Gaussian component is given by r max2＝V un2As the base number, r is added every time a second-level transient track is created max2Adding 1; labeling of first order transient Gaussian components with r max1＝V unAs a base number, r is a transient state track every time a new one is built max1Adding 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 moment k-1|k-1(x) I.e. a priori information about surviving objects is known at time k Is composed of J k-1Parameter set of gaussian components And (ii) an approximation, wherein,
Figure FDA0002225001030000013
denotes the Gaussian density, x is the state of a single target, x ═ x px vy py v] TWherein [ x ] py p] TIndicates the position of a single target, [ x ] vy v] TRepresenting the speed of a single target;
Figure FDA0002225001030000014
and
Figure FDA0002225001030000015
the 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 k S,k|k-1(x)，
Figure FDA0002225001030000016
Is composed of J k-1Parameter set of gaussian components
Figure FDA0002225001030000017
And (ii) an approximation, wherein,
Figure FDA0002225001030000018
and
Figure FDA0002225001030000019
the weight, mean, covariance matrix, and label of the j0 th surviving prediction component at time k,
Figure FDA00022250010300000110
Figure FDA00022250010300000111
wherein p is S,kProbability of survival of a single target, F k-1And Q k-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 moment k(x)， Is composed of J γ,kParameter set of gaussian components
Figure FDA0002225001030000022
Approximation of, wherein
Figure FDA0002225001030000023
And
Figure FDA0002225001030000024
respectively 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)，
Figure FDA0002225001030000025
Is composed of J γ,kParameter set of gaussian components
Figure FDA0002225001030000026
And (ii) an approximation, wherein,
Figure FDA0002225001030000027
and
Figure FDA0002225001030000028
respectively the weight, mean, covariance matrix and label of the ith new prediction component at the moment k,
Figure FDA0002225001030000029
Figure FDA00022250010300000210
merging And
Figure FDA00022250010300000212
obtaining a Gaussian component parameter set
Figure FDA00022250010300000213
J k|k-1＝J k-1+J γ,k，
Figure FDA00022250010300000214
And
Figure FDA00022250010300000215
respectively 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 k k|k-1(x)，
Figure FDA00022250010300000216
Step 3, eliminating clutter
(3a) Measurement set of k time obtained from radar
Figure FDA00022250010300000217
And the prediction height obtained in step 2Parameter set of gaussian component
Figure FDA00022250010300000218
Clutter is removed based on a wave gate method suitable for a GM-PHD filter, and an effective measurement set at the moment k is obtained
Figure FDA00022250010300000219
Wherein the content of the first and second substances,
Figure FDA00022250010300000220
is the xy-plane ze1 measurements obtained by radar,
Figure FDA00022250010300000221
Figure FDA00022250010300000222
is the ze number effective measurement obtained by adopting a wave gate method, and | is the aggregate potential;
(3b) execution of Z k,resi＝Z k-Z k,efObtaining a residual measurement set
Figure FDA00022250010300000223
When k is 1, no gaussian component parameter set is predicted, then Z k,efIs empty, Z k,resi＝Z k；
Step 4, updating Gaussian components
When k is 1, there is no prediction gaussian component parameter set, Z k,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 Z k,efAnd updating the multiple posterior intensities D at the time k k|k(x) Wherein
Figure FDA00022250010300000225
p D,kThe detection probability of a single target at the moment k is (1-p) D,k)D k|k-1(x) Is composed of J k|k-1Parameter set for updating Gaussian component
Figure FDA0002225001030000031
Approximation;
Figure FDA0002225001030000032
wherein the content of the first and second substances,
Figure FDA0002225001030000033
and
Figure FDA0002225001030000034
respectively according to the effective measurement z k,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 radar kThe probability density of the spatial distribution of (a),
Figure FDA0002225001030000036
H kfor the measurement matrix, R kMeasuring a noise covariance matrix;
Figure FDA0002225001030000037
Figure FDA0002225001030000038
i is an identity matrix and is a matrix of the identity,
Figure FDA0002225001030000039
obtaining a multi-target posterior intensity D at an approximate k moment k|k(x)，
Figure FDA00022250010300000310
Is composed of J k|k-1(1+|Z k,ef|) parameter sets for updating gaussian components
Figure FDA00022250010300000311
In the approximation that the difference between the first and second values,
Figure FDA00022250010300000312
and the 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 U w、U m、U P、U lAre respectively J k|k-1(1+|Z k,efI) matrix of dimension, J obtained in step 4 k|k-1(1+|Z k,efI) weight, mean, covariance matrix and label of updated Gaussian components are respectively stored in U w、U m、U PAnd U lIn (1), wherein,
Figure FDA00022250010300000314
are stored to U in sequence wThe first column of (a) is,
Figure FDA00022250010300000315
are stored to U in sequence mThe first column of (a) is,
Figure FDA00022250010300000316
are stored to U in sequence PThe first column of (a) is,
Figure FDA00022250010300000317
are stored to U in sequence lThe first column of (1); press ze ═ 1., | Z k,efL, stored in sequence by Updated parameters of Gaussian component to U w、U m、U PAnd U lIn the specification:
Figure FDA00022250010300000319
are stored to U in sequence wThe ze +1 th column of (a),
Figure FDA00022250010300000320
are stored to U in sequence mThe ze +1 th column of (a),
Figure FDA00022250010300000321
are stored to U in sequence PThe ze +1 th column of (a),
Figure FDA00022250010300000322
are stored to U in sequence lZe +1 column (b);
(5b) order to
Figure FDA00022250010300000323
(5c) Let state estimate set at time k
Figure FDA0002225001030000041
Track label set
Figure FDA0002225001030000042
Set of occurrence moments
Figure FDA0002225001030000043
(5d) Get U wMaximum value of all elements from the second column to the last column If this maximum value is present
Figure FDA0002225001030000045
Less than basic weight threshold w of state extraction 1Directly jumping to the step (5 h); otherwise, according to the maximum value
Figure FDA0002225001030000046
At a different position, from U lTaking the labels at different positions, processing the labels to remove repeated elements to obtain a non-repeated label set
Figure FDA0002225001030000047
Is provided with
Figure FDA0002225001030000048
The system is used for storing a part of column number sets of effective measurement in the query matrix;
(5e) push button For all labels in turn are
Figure FDA00022250010300000410
The processing of the gaussian component specifically includes:
(5e.1) the weight is given
Figure FDA00022250010300000411
The label is
Figure FDA00022250010300000412
Is at U eIn (1), the row number and column number are denoted as rn (i1)And cn (i1)Is recorded in U mElements in corresponding positions in
Figure FDA00022250010300000413
Is m (i1)，U eWherein all tags are
Figure FDA00022250010300000414
Is recorded as
Figure FDA00022250010300000415
(5e.2) judging the tag
Figure FDA00022250010300000416
According to the attribute, the corresponding operation is carried out on the attribute, and the following three conditions exist:
A. when in use
Figure FDA00022250010300000417
Then
Figure FDA00022250010300000418
A label for a new component, wherein V new＞V un2The method comprises the following specific operation steps:
① when
Figure FDA00022250010300000419
Establishing a two-level transient track, wherein w 2For the second-level transient weight threshold, the specific operation is as follows:
r max2＝r max2+1
Figure FDA00022250010300000421
Figure FDA00022250010300000422
Figure FDA00022250010300000423
wherein r is max2＝r max2+1, updating the label of the secondary transient track; note the book Middle (rn) (i1)The elements of a row in all columns are Will be provided with
Figure FDA00022250010300000426
All elements in (1) are changed to r max2(ii) a Note the book
Figure FDA00022250010300000427
Middle rn (i1)The elements of a row in all columns are
Figure FDA00022250010300000428
Note U wMiddle rn (i1)Rows and cn (i1)The elements of the column are
Figure FDA00022250010300000429
Firstly, the method is carried out
Figure FDA00022250010300000430
All elements in (1) are set to zero and then
Figure FDA00022250010300000431
Middle rn (i1)Rows and cn (i1)Elements of a column
Figure FDA00022250010300000432
Instead, it is changed into
Figure FDA00022250010300000433
Respectively counting the transient tracks r max2Secondary transient set stored at time k
Figure FDA00022250010300000434
And first order transient set
Figure FDA00022250010300000435
Performing the following steps;
② when
Figure FDA00022250010300000436
Establishing a first-level transient track, and specifically operating as follows:
r max1＝r max1+1
Figure FDA00022250010300000437
Figure FDA00022250010300000438
Figure FDA00022250010300000439
wherein r is max1＝r max1+1, updating the label of the primary transient track; will be provided with
Figure FDA0002225001030000051
All elements in (1) are changed to r max1(ii) a Firstly, the method is carried out
Figure FDA0002225001030000052
All elements in (1) are set to zero and then
Figure FDA0002225001030000053
Middle in rn (i1)Rows and cn (i1)Elements of a column Instead, it is changed into
Figure FDA0002225001030000055
The number r of transient tracks max1First order transient set stored at time k
Figure FDA0002225001030000056
Performing the following steps;
③ column number cn (i1)Stored into CI, i.e. CI ═ CI cn (i1)](ii) a Changes are made to Is rn (i1)I.e. by
Figure FDA0002225001030000058
B. When in use
Figure FDA0002225001030000059
Then
Figure FDA00022250010300000510
The method is a label of a transient component, and comprises the following specific steps:
① when
Figure FDA00022250010300000511
Then, further on,
Figure FDA00022250010300000512
for labels of the two types of transient components, the following steps are sequentially executed:
first, it executes
Figure FDA00022250010300000513
And
Figure FDA00022250010300000514
secondly, if k is not less than 3 and
Figure FDA00022250010300000515
firstly, a new determined track r is established max＝r max+1, then get single target state estimate
Figure FDA00022250010300000516
And modifying the corresponding parameters
Figure FDA00022250010300000517
Finally updating CI ═ CI cn (i1)]；
② when
Figure FDA00022250010300000518
Then, further on,
Figure FDA00022250010300000519
for a class of transient classified tags, the following steps are performed in sequence:
first, it executes
Figure FDA00022250010300000520
Secondly, if k is more than or equal to 4, And is
Figure FDA00022250010300000522
Firstly, a new determined track r is established max＝r max+1, then get single target state estimate And modifying the corresponding parameters
Figure FDA00022250010300000524
Finally updating CI ═ CI cn (i1)]；
C. When in use
Figure FDA00022250010300000525
Then To determine the label of the component, from labels
Figure FDA00022250010300000527
The method comprises the following specific steps of:
① order w 3Is a weight threshold value of state preprocessing and satisfies w 1≤w 3＜w 2；
If it is not
Figure FDA00022250010300000528
The following state estimation preprocessing steps are performed in sequence:
Figure FDA00022250010300000529
if it is not
Figure FDA00022250010300000531
Wherein the content of the first and second substances,
Figure FDA00022250010300000532
is at U wMiddle label
Figure FDA00022250010300000533
All 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 U mTo (1) a
Figure FDA00022250010300000534
The elements in row number 1 are
Figure FDA00022250010300000535
Known zero mean metrology noise covariance matrix
Figure FDA0002225001030000061
Figure FDA0002225001030000062
And the measured variances of the x-axis and the y-axis are respectively, and the sigma is [ sigma ] xσ y] T，d(1)＝||σ-0||；
Otherwise, go to step ②;
② mixing m (i1)As a single target state estimate, store it to a set of state estimates
Figure FDA0002225001030000064
In, i.e.
Figure FDA0002225001030000065
Will be provided with And k are stored separately to the track tag set
Figure FDA0002225001030000067
And set of occurrence moments
Figure FDA0002225001030000068
In, i.e. And at the same time, adding cn (i1)Store to CI, i.e. CI ═ CI cn (i1)]；
③ weight correction of Gaussian component
Firstly, calculate m one by one (i1)Euclidean distance to the same tag component:
Figure FDA00022250010300000611
and correcting the weight according to the conditions:
Figure FDA00022250010300000612
if
Figure FDA00022250010300000613
wherein, U wAnd U mTo (1) a
Figure FDA00022250010300000614
The elements in the row at the ze1 column are respectively marked as
Figure FDA00022250010300000615
And
Figure FDA00022250010300000616
attenuation factor a1 ═ 1/3;
(5e.3) based on two "one-to-one" principles, U is divided wThe middle row has the number of
Figure FDA00022250010300000617
Is set to zero, i.e. is performed
Figure FDA00022250010300000618
(5f) Is executed completely After 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' principles wAll the elements of any row of (1) are set to zero;
(5g) returning to the step (5 d);
(5h) respectively to be provided with
Figure FDA00022250010300000620
And
Figure FDA00022250010300000621
j in (1) k|k-1(1+|Z k,efI) conversion of elements column by column to 1 × (J) k|k-1(1+|Z k,ef|))) dimensional array, resulting in a modified parameter set that updates the gaussian component
Figure FDA00022250010300000622
(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 n nostThe method comprises the following specific steps of:
(6a.1) taking
Figure FDA0002225001030000071
In non-repetitive tab sets, from which to delete
Figure FDA0002225001030000072
The label set of the component of which the state estimation is not obtained at the moment k is obtained
Figure FDA0002225001030000073
(6a.2) pressing
Figure FDA0002225001030000074
Sequentially discriminating
Figure FDA0002225001030000075
In that
Figure FDA0002225001030000076
k1＝k-n notrThe number of occurrences in k-1: if it is not
Figure FDA0002225001030000077
At each time k1
Figure FDA0002225001030000078
If all appear in (1), then execute
Figure FDA0002225001030000079
(6b) All the weights are larger than a small weight threshold value w etThe sequence number of the Gaussian component of (1) is recorded in I1, i.e.
Figure FDA00022250010300000710
(6c) Let J k|k＝|I1|，
Figure FDA00022250010300000711
By analogy, obtain
Figure FDA00022250010300000712
Wherein the content of the first and second substances, and
Figure FDA00022250010300000714
the 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 6
Figure FDA00022250010300000715
Performing Gaussian component combination to make
Figure FDA00022250010300000716
Is a set of sequence numbers of the Gaussian components, and
Figure FDA00022250010300000717
the method comprises the following specific steps:
(7a) when I2 is empty, go to step (7e), otherwise,
Figure FDA00022250010300000718
obtaining a label
Figure FDA00022250010300000719
(7b) Taking and
Figure FDA00022250010300000720
the sequence number set of the co-tagged and combinable components is I3,
Figure FDA00022250010300000721
wherein d is ThIs a merged distance threshold;
(7c) all gaussian component parameters in the sequence number set I3 are merged,
Figure FDA00022250010300000722
obtaining a new Gaussian component;
(7d) i2 ═ I2-I3, and return to step (7 a);
(7e) obtaining a combined Gaussian component set
Figure FDA00022250010300000725
Step 8, generating new target intensity
According to Z k,resiPress ze2 ═ 1., | Z k,resiObtaining the parameters of 2 th new components according to the following operations:
Figure FDA0002225001030000081
Figure FDA0002225001030000082
Figure FDA0002225001030000083
Figure FDA0002225001030000084
obtaining the prior intensity gamma of the new target at the k +1 moment k+1(x)，
Figure FDA0002225001030000085
J γ,k+1＝|Z k,resi|，i＝ze2；γ k+1(x) Is composed of J γ,k+1Parameter set of gaussian components
Figure FDA0002225001030000086
Approximation;
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
Figure FDA0002225001030000087
Figure FDA0002225001030000088
And
Figure FDA0002225001030000089
the 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.2) setting up effective measurement set
(3a.3)
Figure FDA00022250010300000811
For the latest measurement set obtained at time k, ze1 ═ 1 kL, sequentially judging each
Figure FDA00022250010300000812
Is determined by the attribute of
Figure FDA00022250010300000813
Then, press J ═ 1 k|k-1One by one calculate
Figure FDA00022250010300000814
To each gaussian component mean
Figure FDA00022250010300000815
Is a distance of
Figure FDA00022250010300000816
If d is (ze1,j)D (a) or less, then judging
Figure FDA00022250010300000817
Nearby Gaussian component, see
Figure FDA00022250010300000818
For effective measurement, add it to the effective measurement set Z k,efThen jump out of the loop and execute the pair And (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 1 un＜V un2。
4. The explicit multi-target tracking method under adaptive newness strength based on GM-PHD filter of claim 3 where V is un＝500，V un2＝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 J max1When is coming into contact with When it comes to
Figure FDA00022250010300000821
Before taking J max1The Gaussian component with the maximum weight is recorded as
Figure FDA0002225001030000091
Figure FDA0002225001030000092
When in use, will
Figure FDA0002225001030000093
Relabeling to
Figure FDA0002225001030000094
Thereby obtaining the posterior intensity of survival multiple targets
CN201910949144.4A 2019-10-08 2019-10-08 Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity Active CN110780269B (en)

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 true CN110780269A (en) 2020-02-11
CN110780269B 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)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111562571A (en) * 2020-05-28 2020-08-21 江南大学 Maneuvering multi-target tracking and track maintaining method for unknown new-born strength
CN111665495A (en) * 2020-06-16 2020-09-15 苏州慧至智能科技有限公司 VSMM-GMPLD-based multi-target tracking method
CN111811515A (en) * 2020-07-03 2020-10-23 浙江工业大学 Multi-target track extraction method based on Gaussian mixture probability hypothesis density filter
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
CN112328959A (en) * 2020-10-14 2021-02-05 哈尔滨工程大学 Multi-target tracking method based on adaptive extended Kalman probability hypothesis density filter
CN112688667A (en) * 2020-12-22 2021-04-20 中国人民解放军63921部队 Design method of GM-PHD filter

Citations (5)

* Cited by examiner, † Cited by third party
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

Patent Citations (5)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
Title
吕学斌 等: "高斯混合概率假设密度滤波器在多目标跟踪中的应用", 《计算机学报》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111562571A (en) * 2020-05-28 2020-08-21 江南大学 Maneuvering multi-target tracking and track maintaining method for unknown new-born strength
CN111562571B (en) * 2020-05-28 2022-04-29 江南大学 Maneuvering multi-target tracking and track maintaining method for unknown new-born strength
CN111665495A (en) * 2020-06-16 2020-09-15 苏州慧至智能科技有限公司 VSMM-GMPLD-based multi-target tracking method
CN111811515A (en) * 2020-07-03 2020-10-23 浙江工业大学 Multi-target track extraction method based on Gaussian mixture probability hypothesis density filter
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
CN112328959A (en) * 2020-10-14 2021-02-05 哈尔滨工程大学 Multi-target tracking method based on adaptive extended Kalman probability hypothesis density filter
CN112688667A (en) * 2020-12-22 2021-04-20 中国人民解放军63921部队 Design method of GM-PHD filter
CN112688667B (en) * 2020-12-22 2022-10-11 中国人民解放军63921部队 Design method of GM-PHD filter

Also Published As

Publication number Publication date
CN110780269B (en) 2022-03-25

Similar Documents

Publication Publication Date Title
CN110780269B (en) Explicit multi-target tracking method based on GM-PHD filter under self-adaptive new growth intensity
CN110376581B (en) Explicit multi-target tracking method based on Gaussian mixture probability hypothesis density filter
CN107424171B (en) Block-based anti-occlusion target tracking method
CN114972418B (en) Maneuvering multi-target tracking method based on combination of kernel adaptive filtering and YOLOX detection
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
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
CN103985120A (en) Remote sensing image multi-objective association method
CN116128932B (en) Multi-target tracking method
CN110688940A (en) Rapid face tracking method based on face detection
CN105761276A (en) Iteration RANSAC-based adaptive birth target intensity estimation GM-PHD multi-target tracking algorithm
CN108717702B (en) Probabilistic hypothesis density filtering smoothing method based on segmented RTS
CN111562569A (en) Weighted group sparse constraint-based multi-target constant false alarm detection method under Weibull background
CN115204212A (en) Multi-target tracking method based on STM-PMBM filtering algorithm
CN108153519B (en) Universal design framework for target intelligent tracking method
CN116299195A (en) Radar signal processing method based on TOA sequence relativity
CN113296089B (en) LMB density fusion method and device for multi-early-warning-machine target tracking system
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
CN110880012B (en) Inter-pulse agile radar radiation source frequency information correlation method for multi-reconnaissance platform
Su et al. Real-time multiple object tracking based on optical flow
Yu et al. A hypergraph matching labeled multi-Bernoulli filter for group targets tracking
Zhu et al. The modified probability hypothesis density filter with adaptive birth intensity estimation for multi-target tracking in low detection probability
Gao et al. An explicit track continuity algorithm for the GM-PHD 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