CN104156984A - PHD (Probability Hypothesis Density) method for multi-target tracking in uneven clutter environment - Google Patents

Abstract

The embodiment of the invention provides a PHD (Probability Hypothesis Density) method for multi-target tracking in an uneven clutter environment, and relates to the field of the intelligent information processing. The method can improve the accuracy of the calculated result under the condition that the calculating amount is reduced. The method comprises the following steps: if the surveillance region is a sparse clutter region, the standard PHD can be directly predicted and updated to obtain estimated target amount and target state; if the surveillance region is an intensive clutter region, convex hulls are searched to determine a clutter region; when the PHD is predicted gradually at different moment, return waves not contained in the clutter region are selected for calculation ; when the PHD is updated, the return waves of which the real target measurement exactly drops into the clutter region are added; Gauss components in the posteriori intensity are cut and merged to obtain the target intensity; then, the weight sum of all the Gauss components in the target intensity are calculated to obtain estimated target amount in the surveillance region; the Gauss components with weight greater than Tau in the target intensity are extracted to be in the estimated target state.

Description

A kind of probability hypothesis density method of multiple target tracking under inhomogeneous clutter environment

Technical field

The present invention relates to Intelligent Information Processing field, relate in particular to the probability hypothesis density method of multiple target tracking under a kind of inhomogeneous clutter environment.

Background technology

Multiple target tracking be research from target measure and clutter the method for estimating target number and each dbjective state.The method of traditional processing multiple target tracking problem is mostly based on data correlation, need to set up the corresponding relation of measurement and target, as nearest neighbor method (NN), Joint Probabilistic Data Association algorithm (JPDA), many Hypothesis Tracking Algorithms (MHT).In recent years, obtain very large concern based on the theoretical multiple target tracking algorithm of random finite set (Random Finite Set, RFS).Probability hypothesis density (the Probability Hypothesis Density that Mahler proposes, PHD) filtering algorithm is exactly a kind of under random finite set framework, compare complete multi-objective Bayesian wave filter, computational complexity obtains the algorithm effectively reducing.

In traditional PHD wave filter, conventionally hypothesis clutter in observable scope is equally distributed, but the space distribution of clutter is inhomogeneous in observation area in actual environment, clutter district can be divided into sparse clutter district and dense clutter district, in dense clutter district, the number of false-alarm is far longer than the measurement of real goal, and carrying out in this case multiple target tracking is a difficult problem.Therefore,, in the time of clutter skewness, the tracking performance of traditional PHD wave filter can sharply decline.For this problem, in prior art, having a kind of method is in the process of PHD particle filter, by setting thresholding, by the observed reading filtering being positioned at outside thresholding, only utilizes the observed reading that is positioned at thresholding to carry out the renewal of particle weight.Also having a kind of method is to propose a kind of application target initial sum to maintain the auxiliary SMC-PHD filtering method of rule, each particle is except state vector and weights, two auxiliary variables are also increased, one is used for representing the living age of particle, and one is used for representing that nearest N cycle particle is subject to more news of observed reading.

Above-mentioned these methods, it is all the PHD method for target following under dense clutter environment, apply these methods and have at existing sparse clutter under the complicated clutter environment of dense clutter, will cause in sparse clutter district and dense clutter district calculated amount all very large, can not realize quick computing, but originally need not be so large in sparse clutter district calculated amount, apply these methods and cause on the contrary calculated amount large, and the result obtaining is also not necessarily desirable.

Summary of the invention

Embodiments of the invention provide the probability hypothesis density method of multiple target tracking under a kind of inhomogeneous clutter environment, can, in the situation that calculated amount reduces, improve the accuracy of result of calculation.

For achieving the above object, embodiments of the invention adopt following technical scheme:

A probability hypothesis density method for multiple target tracking under inhomogeneous clutter environment, comprising:

Whether the number of echoes that judges n frame accumulation in monitor area is greater than threshold value ε, wherein, described n, ε is default value;

If not, directly carry out standard P HD prediction and PHD and upgrade, obtain number of targets and the dbjective state estimated;

If so, find convex closure, determine clutter district;

Initialization, the posteriority intensity of acquisition initial time, described initial time was 0 moment;

For the Echo Rating that the k moment records, select the echo that is not included in described clutter district to carry out PHD prediction and calculation, according to the posteriority intensity in k-1 moment, obtain the predicted intensity in k moment, wherein, described k is more than or equal to 1 integer;

According to the predicted intensity in described k moment, with upgrade echo carry out PHD renewal, obtain the posteriority intensity in k moment, wherein, the echo of described renewal comprises the echo that is not included in described clutter district described in the Echo Rating that the k moment records, and the measurement of k moment real goal just falls into the echo in described clutter district; Described posteriority intensity is Gaussian Mixture form, comprises weights, average and the covariance of each gauss component;

Weights in described posteriority intensity are less than to the gauss component filtering of T, the gauss component that the distance between gauss component is less than to U merges, and obtains target strength; Wherein, U≤10, T≤10 ^-5;

Calculate all gauss components in described target strength weights and, obtain the estimating target number in described monitor area; Extract the gauss component that weights in described target strength are greater than τ, as the dbjective state of estimating, wherein, described τ >=0.5.

Optionally, described searching convex closure, determines clutter district, comprising:

Approximate (AP) clustering algorithm of propagating of application carries out cluster to the echo of the n frame accumulation in described monitor area, obtains cluster result;

Count and be greater than the region of η for cluster in described cluster result, find convex closure by Graham scan method, determine clutter district, described η determines according to simulating scenes.

The probability hypothesis density method of multiple target tracking under the inhomogeneous clutter environment that technique scheme provides, utilize AP clustering algorithm, first the echo of multiple frame cumulation in monitor area is carried out to cluster, determine clutter district with finding the method for convex closure, and then by observing the moment carry out PHD prediction and renewal, in the time that PHD predicts, without the echo in clutter district, but in the time carrying out PHD renewal, the echo that the measurement of the real goal of discovery just need to be fallen into clutter district takes out and carries out PHD renewal, the high precision that can ensure like this algorithm, makes estimated result more accurate.In addition, in the time that being greater than certain thresholding, the number of echoes in monitor area just carries out the computing in cluster and definite clutter district, otherwise not doing this computing directly carries out the PHD prediction of standard and upgrades, computation process is simple like this, for inhomogeneous clutter, district can realize quick computing, thereby realizes the lifting of PHD method tracking performance under inhomogeneous clutter environment.

Brief description of the drawings

The schematic flow sheet of the probability hypothesis density method of multiple target tracking under a kind of inhomogeneous clutter environment that Fig. 1 provides for the embodiment of the present invention;

Between a kind of AP clustering algorithm data point that Fig. 2 provides for the embodiment of the present invention, disappear and transmit schematic diagram;

A kind of Graham of application scan algorithm that Fig. 3 provides for the embodiment of the present invention is found the schematic diagram of convex closure;

Under a kind of inhomogeneous clutter environment that Fig. 4 provides for the embodiment of the present invention, carry out the simulating scenes figure of multiple target tracking;

A kind of method emulation provided by the invention state estimation figure out that applies that Fig. 5 provides for the embodiment of the present invention;

A kind of method emulation provided by the invention number of targets drawing for estimate out of applying that Fig. 6 provides for the embodiment of the present invention;

A kind of Gaussian Mixture PHD method emulation state estimation figure out that applies that Fig. 7 provides for the embodiment of the present invention;

A kind of Gaussian Mixture PHD method emulation number of targets drawing for estimate out of applying that Fig. 8 provides for the embodiment of the present invention.

Embodiment

Below in conjunction with the accompanying drawing in the embodiment of the present invention, the technical scheme in the embodiment of the present invention is clearly and completely described, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiment.Based on the embodiment in the present invention, those of ordinary skill in the art, not making the every other embodiment obtaining under creative work prerequisite, belong to the scope of protection of the invention.

The embodiment of the present invention provides the probability hypothesis density method of multiple target tracking under a kind of inhomogeneous clutter environment, as shown in Figure 1, said method comprising the steps of:

101, judge whether the number of echoes that in monitor area, n frame is accumulated is greater than ε.

Wherein, described n, ε is default value; Described ε is preset by user according to the number of echoes in monitor area, and described n is also default by user oneself.

If the number of echoes of n frame accumulation is greater than ε in monitor area, show that described monitor area is dense clutter region, carries out step 103-107; If the number of echoes of n frame accumulation is less than ε in monitor area, show that described monitor area is sparse clutter region, carry out step 102.

102, directly carry out standard P HD prediction and PHD and upgrade, obtain number of targets and the dbjective state estimated.

Here described standard P HD prediction and PHD upgrades, and is traditional PHD prediction of the prior art and PHD and upgrades, and those skilled in the art have a clear understanding of its computation process, therefore no longer repeat at this.

103, find convex closure, determine clutter district.

Optionally, step 103 comprises the following steps:

A, application AP clustering algorithm carry out cluster to the echo of the n frame accumulation in described monitor area, obtain cluster result.

AP (Affinity Propagation, approximate propagation) clustering algorithm is to carry out cluster according to the similarity between N data point, these similarities can be symmetrical, i.e. two data points similarity the same (as Euclidean distance) between mutually; Can be asymmetric, two data points similarity between be mutually etc. yet.The similarity matrix S (wherein N is for there being N data point) of these similarity compositions N × N.AP algorithm does not need to specify in advance clusters number, and it,, using all data points all as potential cluster centre, is referred to as exemplar on the contrary.

The main thought of AP clustering algorithm is: with the numerical value S on the diagonal line of s-matrix (k ', k ') as k ' some judgment criteria that can become cluster centre, this means that this value is larger, the possibility that this point becomes cluster centre is also just larger, and this value is with reference to degree P.The quantity of cluster is subject to the impact with reference to degree P, if think that each data point is likely as cluster centre, P just should get identical value so.If get the average of similarity of input as the value of P, it is medium obtaining number of clusters.If the minimum value of getting, obtains the less cluster of class number.In AP clustering algorithm, transmit the message of two types, i.e. Responsibility and Availability.R (i ', k ') represents to send to from an i ' the numerical value message of candidate's cluster centre k ', reflects at k ' and whether is suitable as the cluster centre of i '.Whether A (i ', k ') represents to send to from candidate's cluster centre k ' the numerical value message of i ', reflect at i ' and select k ' as its cluster centre.R (i ', k ') and A (i ', k ') stronger, k ' some possibility as cluster centre is just larger, and i ' is under the jurisdiction of taking k ' some possibility as the cluster of cluster centre also larger.AP clustering algorithm is constantly updated Attraction Degree and the degree of membership value of each point by iterative process, until produce M high-quality exemplar, remaining data point is assigned in corresponding cluster simultaneously.

Here it should be noted that:

1, exemplar: refer to cluster centre.

2, the similarity of similarity: data point i ' and some j ' is designated as S (i ', j ').It is the similarity of giving directions the cluster centre of j ' conduct point i '.Generally calculate by Euclidean distance.In this algorithm, a little with point similarity be all taken as negative value.Like this, the larger explanation point of similarity value is nearer with the distance of putting, and is convenient to subsequent calculations.

3, the reference degree of preference: data point i ' is called P (i ') or S (i ', i ').To give directions the reference degree of i ' as cluster centre.Generally get the intermediate value of S similarity value.

4, Responsibility:R (i ', k ') is used for describing some k ' and is suitable as the degree of the cluster centre of data point i '.

Availability:A (i ', k ') is used for describing the appropriateness of some i ' selected element k ' as its cluster centre.

As shown in Figure 2, be Responsibility (being designated as R in Fig. 2) and Availability (being designated as A in Fig. 2) relation between the two.

The computing formula of R (i ', k '), A (i ', k ') and R (k ', k ') is:

R(i′，k′)＝S(i′，k′)-max{A(i′，j′)+S(i′，j′)}(j′∈{1，2，…，N}&j′≠k′) (1)

R(k′，k′)＝P(k′)-max{A(k′，j′)+S(k′，j′)}(j′∈{1，2，…，N}&j′≠k′) (3)

Can find out above, in the time that P (k ') makes more greatly R (k ', k ') larger, A (i ', k ') also larger, thus class represents that the possibility of the final cluster centre of k ' conduct is larger; Equally, in the time that more P (i ') are larger, more classes represent the cluster centre centered by tending to into.Therefore, increase or reduce P and can increase or reduce the clusters number of AP output.

5, Damping factor (ratio of damping) λ: be mainly converging action.Each iteration, Attraction Degree R _{i '}with degree of membership A _{i '}will with last R _{i '-1}and A _{i '-1}be weighted renewal.Formula is as follows:

R _i′＝(1-λ)*R _i′+λ*R _i′-1 (4)

A _i′＝(1-λ)*A _i′+λ*A _i′-1 (5)

Wherein, λ ∈ [0.5,1).

The concrete steps of AP clustering algorithm are as follows:

(1) calculate N the similarity value between point, value is placed in s-matrix, then chooses P value (generally getting the intermediate value of S).

(2) maximum iteration time is set, after iterative process starts, calculate R value and A value each time, according to R (k ', k ')+A (k ', k ') value determine whether cluster centre, when iterations exceed maximal value or in the time that the continuous how many times iteration of cluster centre does not change stop calculate, obtain cluster result.

B, count and be greater than the region of η for cluster in described cluster result, find convex closure by Graham scan method, determine clutter district.

Described η is preset value, can determine according to simulating scenes.

Wherein, Graham scan algorithm is the algorithm of the convex closure of one group of plane finite point of a kind of calculating.Its concrete steps that solve convex closure are as follows:

(1) choose 1 H of y coordinate minimum in a little, be used as basic point.If having the y coordinate of multiple points is all minimum value, choose a bit of x coordinate minimum.The point that coordinate is identical should be got rid of.Then the vectorial < H forming according to other each point P and basic point, the angle of P > and x axle sorts, and angle scans from large to small clockwise, otherwise scans counterclockwise.In realization, without trying to achieve angle, only need obtain vectorial mould according to vectorial inner product formula.Taking Fig. 3 as example, basic point is H, after sorting from small to large, is followed successively by H, K, C, D, L, F, G, E, I, B, A, J according to angle.Scan counterclockwise below.

(2) line segment < H, K > mono-fixes on convex closure, then adds C.Suppose line segment < H, C > also on convex closure because with regard to H, K, 3 of C, their convex closure be exactly thus 3 institutes form.But while next adding D, can find, line segment < K, D > just can be on convex closure, so by line segment < K, C > gets rid of, C point can not be convex closure.When adding when some, must consider whether line segment above there will be on convex closure.

(3), from basic point, the sense of rotation of the line segment that on convex closure, every is faced mutually should be consistent, and with scanning opposite direction.If find that the point newly adding changes the sense of rotation of new line segment and upper line segment, can judge upper a bit inevitable not on convex closure.While realization, availability vector cross product judges, establishing the point newly adding is P _n+1, upper is some P _n, then be a bit above P _n-1.While scanning clockwise, if vectorial < P _n-1, P _n> and < P _n, P _n+1the cross product of > is for just (scanning and determine whether to bear counterclockwise), by a upper point deletion.Delete procedure need to be recalled, and the point that all cross product symbols are contrary before is all deleted, and then will newly put and add convex closure.In Fig. 3, while adding K point, due to line segment < H, K > is with respect to < H, and C >, for turning clockwise, so C point is not on convex closure, should delete, and retains K point.Then add D point, due to line segment < K, the relative < H of D >, K > is for being rotated counterclockwise, therefore D point retains.Scan according to above-mentioned steps, until point concentrates all points all to travel through, obtain the convex closure being formed by a H, J, B, G, D and K line as shown in Figure 3, determine clutter district.

Certainly, find convex closure, determine that can also there be additive method of the prior art in clutter district, here no longer describe in detail.

1041, initialization, the posteriority intensity of acquisition initial time, described initial time was 0 moment.

Carry out PHD when prediction entering simulating scenes, should hypothetical target motion model and observation model all meet linear Gauss's condition; Survival probability and the detection probability of hypothetical target are separate, and irrelevant with dbjective state; The intensity of supposing newborn target random set and derivative goal random set is Gaussian Mixture form; Use J ₀the weighted sum of individual gauss component is carried out initialization.

Hypothetical target motion model and observation model all meet linear Gauss's condition:

Target movement model: f _k|k-1(x| ζ)=N (x; F _k-1ζ, Q _k-1) (6)

Observation model: g _k(z|x)=N (z; H _kx, R _k) (7)

Wherein, N (; M, P) represent that its density average is m, the Gaussian distribution that variance is P; F _k-1for state-transition matrix, Q _k-1for system noise covariance matrix; H _kfor observing matrix, R _kfor measuring noise covariance matrix.

Survival probability and the detection probability of hypothetical target are separate, and irrelevant with dbjective state:

Survival probability: p _{s, k}(x)=p _{s, k}(8)

Detection probability: p _{d, k}(x)=p _{d, k}(9)

The intensity of supposing newborn target random set and derivative goal random set is Gaussian Mixture form:

The intensity of newborn target random set: ${\mathrm{\&gamma;}}_k\left(x\right)=\underset{i=1}{\overset{J_{\mathrm{\&gamma;},k}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{\&gamma;},k}^{\left(i\right)}N(x;m_{\mathrm{\&gamma;},k}^{\left(i\right)},P_{\mathrm{\&gamma;},k}^{\left(i\right)})---\left(10\right)$

The intensity of derivative goal random set:

${\mathrm{\&beta;}}_{k|k-1}\left(x\right|\mathrm{\&zeta;})=\underset{j=1}{\overset{J_{\mathrm{\&beta;},k}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{\&beta;},k}^{\left(j\right)}N(x;F_{\mathrm{\&beta;},k-1}^{\left(j\right)}\mathrm{\&zeta;}+d_{\mathrm{\&beta;},k-1}^{\left(j\right)},Q_{\mathrm{\&beta;},k-1}^{\left(j\right)})---\left(11\right)$

Wherein, J _{γ, k}, i=1,2 ..., J _{γ, k}determined the intensity of newborn target random set etc. parameter, for weights, average and the covariance of i gauss component of newborn target random set intensity, J _{γ, k}it is the number of newborn target gauss component of k moment.In like manner: J _{β, k}, j=1,2 ..., J _{β, k}determined the intensity of the target random set being derived by target ζ etc. parameter, for weights, average and the covariance of i gauss component of newborn target random set intensity, J _{β, k}it is the number of newborn target gauss component of k moment.Here the value in the time that initial time was 0 moment of described parameters is all known.Value in the time that the value in follow-up moment can be 0 moment according to initial time is derived and is calculated.

Use J ₀the weighted sum of individual gauss component is carried out initialization, obtains the posteriority intensity of initial time:

$v_0=\underset{i=1}{\overset{J_0}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_0^{\left(i\right)}N(x;m_0^{\left(i\right)},P_0^{\left(i\right)})---\left(12\right)$

Wherein N (; M, P) be that average is m, the Gaussian distribution that variance is P.Its weights and be expect initial target number: other correlation parameters J in formula (12) ₀, apply the present invention while carrying out emulation user, set own setting by user according to concrete simulated environment.

1042, the Echo Rating obtaining for the k moment, selects the echo that is not included in described clutter district to carry out PHD prediction and calculation, according to the posteriority intensity in k-1 moment, obtains the predicted intensity in k moment.

K moment, k-1 moment in the embodiment of the present invention are the observation moment, and wherein the k-1 moment is the last observation moment in k moment.The k-1=0 moment is initial time, now, and k=1.

When carry out PHD prediction by the moment, for the Echo Rating that each observation moment obtains, if wherein there be the clutter district of echo in step 103, remove those echoes, then carry out PHD prediction.

The posteriority intensity of supposing the k-1 moment has following mixed Gaussian form:

$v_{k-1}\left(x\right)=\underset{i=1}{\overset{J_{k-1}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{\left(i\right)}N(x;m_{k-1}^{\left(i\right)},P_{k-1}^{\left(i\right)})---\left(13\right)$

In the time of k-1=0, formula (13) can use formula (12) to represent:

$v_{k-1}\left(x\right)=v_0=\underset{i=1}{\overset{J_0}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_0^{\left(i\right)}N(x;m_0^{\left(i\right)},P_0^{\left(i\right)})---\left(13a\right)$

The Echo Rating obtaining for the k moment, selects the echo that is not included in described clutter district to carry out PHD prediction and calculation, according to the posteriority intensity in k-1 moment, obtains the predicted intensity in k moment, and the predicted intensity in k moment is also Gaussian Mixture form, provides as follows:

v _k|k-1(x)＝v _S，k|k-1(x)+v _β，k|k-1(x)+γ _k(x) (14)

Wherein, γ _k(x) in formula (10), provide.

$v_{S,k|k-1}\left(x\right)=p_{S,k}\underset{j=1}{\overset{J_{k-1}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{\left(j\right)}N(x;m_{S,k|k-1}^{\left(j\right)},P_{S,k|k-1}^{\left(j\right)})---\left(15\right)$

$m_{S,k|k-1}^{\left(j\right)}=F_{k-1}{m}_{k-1}^{\left(j\right)}---\left(16\right)$

$P_{S,k|k-1}^{\left(j\right)}=Q_{k-1}+F_{k-1}{P}_{k-1}^{\left(j\right)}{F}_{k-1}^T---\left(17\right)$

$v_{\mathrm{\&beta;},k|k-1}\left(x\right)=\underset{j=1}{\overset{J_{k-1}}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{J_{\mathrm{\&beta;},k}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k-1}^{\left(j\right)}{\mathrm{\&omega;}}_{\mathrm{\&beta;},k}^{\left(l\right)}N(x;m_{\mathrm{\&beta;},k|k-1}^{(j,l)},P_{\mathrm{\&beta;},k|k-1}^{(j,l)})---\left(18\right)$

$m_{\mathrm{\&beta;},k|k-1}^{(j,l)}=F_{\mathrm{\&beta;},k-1}^{\left(l\right)}{m}_{k|k-1}^{\left(j\right)}+d_{\mathrm{\&beta;},k-1}^{\left(l\right)}---\left(19\right)$

$P_{\mathrm{\&beta;},k|k-1}^{(j,l)}=Q_{\mathrm{\&beta;},k-1}^{\left(l\right)}+F_{\mathrm{\&beta;},k-1}^{\left(l\right)}{P}_{\mathrm{\&beta;},k-1}^{\left(j\right)}{\left(F_{\mathrm{\&beta;},k-1}^{\left(l\right)}\right)}^T---\left(20\right)$

By the parameter containing in formula (13) and J _k-1can calculate formula (10) and formula (15)-(20), then just can calculate the predicted intensity v in k moment in formula (14) by formula (10) and formula (15)-(20) _k|k-1(x).

Here, can be by by formula (10), and (15)-(20), the predicted intensity that calculates the k moment is designated as following Gaussian Mixture form:

$v_{k|k-1}\left(x\right)=\underset{i=1}{\overset{J_{k|k-1}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k|k-1}^{\left(i\right)}N(x;m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})---\left(21\right)$

105, according to the predicted intensity in described k moment, with upgrade echo carry out PHD renewal, obtain the posteriority intensity in k moment, wherein, the echo of described renewal comprises the echo that is not included in described clutter district described in the Echo Rating of k moment gained, and the measurement of k moment real goal just falls into the echo in described clutter district.

Described posteriority intensity is Gaussian Mixture form, comprises weights, average and the covariance of each gauss component.

Carrying out PHD while upgrading, need to carry out self-adaptation judgement to the echo in clutter district, in clutter district, be all clutter if determine, and without the measurement of real goal, the echo in clutter district all need not; If but the measurement of finding real goal just falls into clutter district, its taking-up is put into PHD and upgrades the echo of using, change the size in clutter district.Judge whether that the method that falls into clutter district is: calculate and measure by dbjective state, by the echo ratio pair in the measurement drawing and clutter district, measure and fall into clutter district if find, will measure taking-up, put into PHD and upgrade the echo that needs use; Otherwise, while carrying out PHD renewal, use the echo that does not comprise clutter district.

The predicted intensity that can be obtained the k moment by step 1042 is following Gaussian Mixture form:

$v_{k|k-1}\left(x\right)=\underset{i=1}{\overset{J_{k|k-1}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{k|k-1}^{\left(i\right)}N(x;m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})---\left(21\right)$

The posteriority intensity in k moment is also Gaussian Mixture form:

$v_k\left(x\right)=(1-p_{D,k})v_{k|k-1}\left(x\right)+\underset{z\&Element;Z_k}{\mathrm{\&Sigma;}}{v}_{D,k}(x;z)---\left(22\right)$

Wherein,

$v_{D,k}(x;z)=\underset{j=1}{\overset{J_{k|k-1}}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^{\left(j\right)}\left(z\right)N(x;m_{k|k}^{\left(j\right)}\left(z\right),P_{k|k}^{\left(j\right)})---\left(23\right)$

${\mathrm{\&omega;}}_k^{\left(j\right)}\left(z\right)=\frac{p_{D,k}{\mathrm{\&omega;}}_{k|k-1}^{\left(j\right)}{q}_k^{\left(j\right)}\left(z\right)}{k_k\left(z\right)+p_{D,k}{\mathrm{\&Sigma;}}_{l=1}^{J_{k|k-1}}{{\mathrm{\&omega;}}_{k|k-1}^{\left(l\right)}q}_k^{\left(l\right)}\left(z\right)}---\left(24\right)$

$m_{k|k}^{\left(j\right)}\left(z\right)=m_{k|k-1}^{\left(j\right)}+K_k^{\left(j\right)}(z-H_k{m}_{k|k-1}^{\left(j\right)})---\left(25\right)$

$P_{k|k}^{\left(j\right)}=[I-K_k^{\left(j\right)}{H}_k]P_{k|k-1}^{\left(j\right)}---\left(26\right)$

$K_k^{\left(j\right)}=P_{k|k-1}^{\left(j\right)}{H}_k^T{(H_k{P}_{k|k-1}^{\left(j\right)}{H}_k^T+R_k)}^{-1}---\left(27\right)$

By the parameter containing in formula (21) and J _k|k-1can calculate formula (23)-(27), then just can calculate the posteriority intensity v in k moment in formula (22) by formula (23)-(27) _k(x).

106, weights in described posteriority intensity are less than to the gauss component filtering of T, the gauss component that the distance between gauss component is less than to U merges, and obtains target strength; Wherein, U≤10, T≤10 ^-5.

Set threshold value T (T≤10 ^-5), in the time that the weights of some gauss components are less than T, these gauss components being filtered out, concrete formula can be expressed as follows:

${\stackrel{-}{v}}_k\left(x\right)=\frac{{\mathrm{\&Sigma;}}_{l=1}^{J_k}{\mathrm{\&omega;}}_k^{\left(l\right)}}{{\mathrm{\&Sigma;}}_{j=N_p+1}^{J_k}{\mathrm{\&omega;}}_k^{\left(j\right)}}\underset{i=N_p+1}{\overset{J_k}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_k^{\left(i\right)}N(x;m_k^{\left(i\right)},P_k^{\left(i\right)})---\left(28\right)$

Wherein, i=N _p+ 1 ..., J _kthe weights that are greater than threshold value.

Set a threshold value U (U < 10), in the time that the distance between some gauss components is less than U, these gauss components are merged.Provide a maximum gauss component number J who allows _max.

Make l=0, $I=\{i=\mathrm{1,2},...,J_k|{\mathrm{\&omega;}}_k^{\left(i\right)}>T\}.$

Repeat

l：＝l+1 (29)

$j=\arg \underset{i\&Element;I}{\max }{\mathrm{\&omega;}}_k^{\left(i\right)}---\left(30\right)$

$L:=\{i\&Element;I|{(m_k^{\left(i\right)}-m_k^{\left(j\right)})}^T{\left(P_k^{\left(i\right)}\right)}^{-1}(m_k^{\left(i\right)}-m_k^{\left(j\right)})\&le;U\}---\left(31\right)$

${\widetilde{\mathrm{\&omega;}}}_k^{\left(l\right)}=\underset{i\&Element;L}{\mathrm{\&Sigma;}}{\widetilde{\mathrm{\&omega;}}}_k^{\left(i\right)}---\left(32\right)$

${\tilde{m}}_k^{\left(l\right)}=\frac{1}{{\widetilde{\mathrm{\&omega;}}}_k^{\left(l\right)}}\underset{i\&Element;L}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_k^{\left(i\right)}{x}_k^{\left(i\right)}---\left(33\right)$

${\tilde{P}}_k^{\left(l\right)}=\frac{1}{{\widetilde{\mathrm{\&omega;}}}_k^{\left(l\right)}}\underset{i\&Element;L}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_k^{\left(i\right)}(P_k^{\left(i\right)}+({\tilde{m}}_k^{\left(l\right)}-m_k^{\left(i\right)}){({\tilde{m}}_k^{\left(l\right)}-m_k^{\left(i\right)})}^T)---\left(34\right)$

I：＝I\L (35)

Until

If l > is J _max, use J _maxthe gauss component that individual weights are larger replaces obtain be the gauss component after merging.

Its concrete merging process those skilled in the art have a clear understanding of, and are not described in detail in this.

The described target strength of carrying out after Gauss's merging is also Gaussian Mixture form, can be as follows:

$\tilde{v}\left(x\right)=\underset{i=1}{\overset{J_{\max }}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&omega;}}}_k^{\left(i\right)}N(x;{\tilde{m}}_k^{\left(i\right)},{\tilde{P}}_k^{\left(i\right)})---\left(36\right)$

107, calculate all gauss components in described target strength weights and, obtain the estimating target number in described monitor area; Extract the gauss component that weights in described target strength are greater than τ, as the dbjective state of estimating, wherein, described τ >=0.5.

Obtain estimating target number: weights by mixed Gaussian item and obtain the number of targets in monitor area,

${\hat{n}}_k=\underset{i=1}{\overset{J_{\max }}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&omega;}}}_k^{\left(i\right)}---\left(37\right)$

Obtain dbjective state: extract the gauss component that weights are greater than τ (τ >=0.5),

${\hat{X}}_k=\{{\tilde{m}}_k^{\left(i\right)}:{\widetilde{\mathrm{\&omega;}}}_k^{\left(i\right)}>\mathrm{\&tau;}\}---\left(38\right)$

In order to verify the present invention program, provide below towards the adaptive probability assumed density method simulation result of multiple target tracking under inhomogeneous clutter environment.

Simulating scenes: the scene of considering to have in monitor area under inhomogeneous clutter environment three targets.

Target 1 occurs in the k=1 moment, and original state is [25m, 4m/s, 10m, 2m/s] ^t, disappear in the k=38 moment; Target 2 occurs in the k=10 moment, and original state is [550m ,-15m/s, 100m, 4m/s] ^t, disappear in the k=60 moment; Target 3 occurs in the k=22 moment, and original state is [270m, 18m/s ,-190m, 6m/s] ^t, disappear in the k=51 moment.

Target state equation is:

X _k+1＝F _kX _k+v _k (39)

Wherein, state vector X _k=[x _{1, k}, x _{2, k}, x _{3, k}, x _{4, k}] ^t, [x _{1, k}, x _{3, k}] ^tthe target location in k moment, [x _{2, k}, x _{4, k}] ^tit is the target velocity in k moment.Sampling interval T=1s.State-transition matrix F _k, process noise covariance be respectively $F_k=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],Q=\left[\begin{array}{cccc}\frac{T^4}{4}& 0& \frac{T^3}{2}& 0\\ 0& T^2& 0& \frac{T^3}{2}\\ \frac{T^3}{2}& 0& \frac{T^4}{4}& 0\\ 0& \frac{T^3}{2}& 0& T^2\end{array}\right].$ In formula, process noise is the white Gaussian noise of zero-mean, and separate with measurement noise sequence.Target survival probability is p _s=0.975.

Measurement equation is:

Z _k＝H _kX _k+ω _k (40)

In formula, measurement noise is the white noise of zero-mean.Measurement matrix and measurement noise covariance are respectively $H_k=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],R=\left[\begin{array}{cc}100& 0\\ 0& 100\end{array}\right].$ Acquisition probability is p _d=0.96.

The space distribution of clutter is heterogeneous, very intensive at certain several regions clutter.Whole supervision spatial domain is [400m, 600m] × [400m, 800m], and in whole region, the average clutter number in each observation moment is 5.Dense clutter region is: in [200m ,-150m] × [200m ,-150m] region, the average clutter number in each observation moment is 20; In [50m, 100m] × [0m, 100m] region, the average clutter number in each observation moment is 20; In [300m, 350m] × [100m, 150m] region, the average clutter number in each observation moment is 20.

Simulating scenes is shown in Fig. 4.

When emulation, n=5, ε=200, η=80.

In emulation experiment, and the inventive method and traditional Gaussian Mixture PHD (referred to as: GM-PHD) method (a kind of implementation of standard P HD) contrasts.Simulation result of the present invention is shown in Fig. 5 and Fig. 6, and the simulation result of Gaussian Mixture PHD is shown in Fig. 7 and Fig. 8.

Contrast two kinds of methods, can find out, under inhomogeneous clutter environment, the in the situation that of having again dense clutter district in existing sparse clutter district, the present invention's (Fig. 5 and Fig. 6 show) can better estimate dbjective state and the number of targets in each moment, as can be seen from Figure, compares Gaussian Mixture PHD, the present invention estimates that the precision of each moment dbjective state and number of targets increases significantly, and has reached the object that promotes performance of target tracking.

The above; be only the specific embodiment of the present invention, but protection scope of the present invention is not limited to this, any be familiar with those skilled in the art the present invention disclose technical scope in; can expect easily changing or replacing, within all should being encompassed in protection scope of the present invention.Therefore, protection scope of the present invention should described be as the criterion with the protection domain of claim.

Claims (2)

1. a probability hypothesis density method for multiple target tracking under inhomogeneous clutter environment, is characterized in that, comprising:

Whether the number of echoes that judges n frame accumulation in monitor area is greater than threshold value ε, wherein, described n, ε is default value;

If not, directly carry out standard P HD prediction and PHD and upgrade, obtain number of targets and the dbjective state estimated;

If so, find convex closure, determine clutter district;

Initialization, the posteriority intensity of acquisition initial time, described initial time was 0 moment;

For the Echo Rating that the k moment records, select the echo that is not included in described clutter district to carry out PHD prediction and calculation, according to the posteriority intensity in k-1 moment, obtain the predicted intensity in k moment, wherein, described k is more than or equal to 1 integer;

According to the predicted intensity in described k moment, with upgrade echo carry out PHD renewal, obtain the posteriority intensity in k moment, wherein, the echo of described renewal comprises the echo that is not included in described clutter district described in the Echo Rating that the k moment records, and the measurement of k moment real goal just falls into the echo in described clutter district; Described posteriority intensity is Gaussian Mixture form, comprises weights, average and the covariance of each gauss component;

Weights in described posteriority intensity are less than to the gauss component filtering of T, the gauss component that the distance between gauss component is less than to U merges, and obtains target strength; Wherein, U≤10, T≤10 ^-5;

Calculate all gauss components in described target strength weights and, obtain the estimating target number in described monitor area; Extract the gauss component that weights in described target strength are greater than τ, as the dbjective state of estimating, wherein, described τ >=0.5.

2. method according to claim 1, is characterized in that, described searching convex closure, determines clutter district, comprising:

Approximate (AP) clustering algorithm of propagating of application carries out cluster to the echo of the n frame accumulation in described monitor area, obtains cluster result;

Count and be greater than the region of η for cluster in described cluster result, find convex closure by Graham scan method, determine clutter district, described η determines according to simulating scenes.