Tracking of Midcourse Ballistic Target Group on Space-based Infrared <br/>Focal Plane using GM-CPHD Filter

Tracking of midcourse ballistic target group on space-based infrared focal plane plays a key role in the space-based early warning system. This paper proposes the Gaussian-mixture cardinalized probability hypothesis density (GM-CPHD) filter to solve this problem. The multi-target states and measurements on infrared focal plane are modeled by random finite set (RFS). The intensity function of RFS of multi-target states and the probability distribution of target number are jointly propagated by cardinalized probability hypothesis density (CPHD) recursion. Under the assumptions of linear Gaussian multi-target models, the Gaussian-mixture implementations of CPHD are presented, and the target number and the multi-target states on infrared focal plane are estimated. In order to enable track continuity, we propose 0-1 integer programming to associate the estimated states between frames. The simulation results show that the GM-CPHD filter can dramatically improve the accuracy of estimated target number and estimated target states compared with the Gaussian-mixture probability hypothesis density filter, and that the track continuity can be successfully achieved.

During the midcourse stage of a ballistic missile, many decoys are released to form closely spaced target group. For the awareness of threat situation, the early warning satellite continually stares this group by infrared sensor to acquire trajectories on infrared focal plane1-3 . However, the target number is time-varying and a lot of clutter appears in the images due to resident space objects (RSO), stars and cosmic rays. These factors present a challenge for tracking of target group on infrared focal plane.

Tracking of target group on infrared focal plane is a multitarget tracking (MTT) problem which involves the joint estimation of an unknown and time-varying number of targets as well as their individual states with uncertain data association. The joint probabilistic data association (JPDA)4 considers associations that survive gating and combines these associations in proportion to their likelihoods. But it can not cope with time-varying number of targets. Multiple hypothesis tracking (MHT)5 forms multiple data association hypotheses (new target, surviving target, false alarm). However, due to the number of hypotheses growing exponentially with the number of targets and clutter, MHT is very time consuming. Mahler has proposed probability hypothesis density (PHD) filter6 for MTT based on random finite set (RFS) theory. This filter models multitarget states and measurements using RFS, and recursively propagates the intensity function (Its integral in any region on state space is the expected number of targets in that region) of RFS of target states. The PHD filter operates on the single-target state space and avoids the sophisticated data association. The sequential Monte Carlo PHD (SMC-PHD) filter7-10 and the Gaussian-mixture PHD (GM-PHD) filter11-14 are two implementations of PHD filter. Mahler further proposes CPHD filter15 which jointly propagates the intensity function and the probability distribution of target number (called the cardinality distribution). Then, the accuracy of the estimation of target number can be improved using the maximum a posteriori (MAP) estimator. Vo16, et al. has proposed implementations of CPHD filter called Gaussian-mixture CPHD (GM-CPHD). Lin3, et al. has proposed SMC-PHD filter for tracking of target group on space-based infrared focal plane, but the track continuity is not achieved.

This paper proposes GM-CPHD filter for tracking of midcourse ballistic target group on space-based infrared focal plane. The multitarget states and measurements on infrared focal plane are modelled by RFS. The intensity function and the cardinality distribution are jointly propagated by CPHD recursion. Under the assumptions of linear Gaussian multitarget models, the Gaussian-mixture implementations of CPHD are presented. We propose 0-1 integer programming to associate the estimated states between frames for track continuity.

The single-target state at time k is a vector of position and velocity on two-dimensional infrared focal plane1-3, and follows a linear Gaussian motion model:                   (1)
where          $F=\left[\begin{array}{cccc}1& 0& T& 0\\ 0& 1& 0& T\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$, T denotes the sampling period, $G=[ T 2 /2 0 0 T 2 /2 T 0 0 T ] MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WHhbGaeyypa0ZaamWaaeaafaqabeabcaaaaeaacaWGubWaaWbaaSqa beaacaaIYaaaaOGaai4laiaaikdaaeaacaaIWaaabaGaaGimaaqaai aadsfadaahaaWcbeqaaiaaikdaaaGccaGGVaGaaGOmaaqaaiaadsfa aeaacaaIWaaabaGaaGimaaqaaiaadsfaaaaacaGLBbGaayzxaaaaaa@475C@$,and $δ k = ( δ x,k , δ y,k ) T MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WH0oWaaSbaaSqaaiaadUgaaeqaaOGaeyypa0Jaaiikaiabes7aKnaa BaaaleaacaWG4bGaaiilaiaadUgaaeqaaOGaaiilaiabes7aKnaaBa aaleaacaWG5bGaaiilaiaadUgaaeqaaOGaaiykamaaCaaaleqabaGa amivaaaaaaa@47DA@$ is the process noise. ${\delta }_{x,k}$ and ${\delta }_{y,k}$ are independent zero-mean Gaussian noise. The covariance of $G{\delta }_{k}$ is denoted by ${Q}_{k}$.

The single-target measurement ${z}_{k}={\left({z}_{x,k},{z}_{y,k}\right)}^{T}$ are the coordinates of two-dimensional focal plane at time K , and the measurement model is linear Gaussian:
${z}_{k}=H{x}_{k}+{\epsilon }_{k}$                      (2)

where $H=[ 1 0 0 0 0 1 0 0 ] MathType@MTEF@5@5@+= feaagGart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WHibGaeyypa0ZaamWaaeaafaqabeGaeaaaaeaacaaIXaaabaGaaGim aaqaaiaaicdaaeaacaaIWaaabaGaaGimaaqaaiaaigdaaeaacaaIWa aabaGaaGimaaaaaiaawUfacaGLDbaaaaa@421E@$, ${\epsilon }_{k}={\left({\epsilon }_{x,k},{\epsilon }_{y,k}\right)}^{T}$ is the measurement noise. ${\epsilon }_{x,k}$ and ${\epsilon }_{y,k}$ are independent zero-mean Gaussian noise. The covariance of ${\epsilon }_{k}$ is denoted by ${R}_{k}$.

The sets of target states and measurements at time are denoted by

$X k ={ x k,1 ,⋯, x k,| X k | } MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGybWaaSbaaSqaaiaadUgaaeqaaOGaeyypa0Jaai4EaiaahIhadaWg aaWcbaGaam4AaiaacYcacaaIXaaabeaakiaacYcacqWIVlctcaGGSa GaaCiEamaaBaaaleaacaWGRbGaaiilaiaacYhacaWGybWaaSbaaWqa aiaadUgaaeqaaSGaaiiFaaqabaGccaGG9baaaa@4B33@$ $Z k ={ z k,1 ,⋯, z k,| Z k | } MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGAbWaaSbaaSqaaiaadUgaaeqaaOGaeyypa0Jaai4EaiaahQhadaWg aaWcbaGaam4AaiaacYcacaaIXaaabeaakiaacYcacqWIVlctcaGGSa GaaCOEamaaBaaaleaacaWGRbGaaiilaiaacYhacaWGAbWaaSbaaWqa aiaadUgaaeqaaSGaaiiFaaqabaGccaGG9baaaa@4B3B@$
where $|{X}_{k}|$ $|{Z}_{k}|$ and are the cardinality of ${X}_{k}$ and ${Z}_{k}$ respectively.

Let ${S}_{k|k-1}\left(\varsigma \right)$ be the surviving RFS at time k that evolved from a target with state $\varsigma$ at time K - 1, and ${\Gamma }_{k}$ the RFS of spontaneous births at time K . The set of multi-target states can be written as (we do not consider target spawning)

$X k =[ ∪ ς∈ X k−1 S k|k−1 (ς) ]∪ Γ k MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGybWaaSbaaSqaaiaadUgaaeqaaOGaeyypa0ZaamWaaeaadaWfqaqa aiabgQIiidWcbaGaaCOWdiabgIGiolaadIfadaWgaaadbaGaam4Aai abgkHiTiaaigdaaeqaaaWcbeaakiaadofadaWgaaWcbaGaam4Aaiaa cYhacaWGRbGaeyOeI0IaaGymaaqabaGccaGGOaGaaCOWdiaacMcaai aawUfacaGLDbaacqGHQicYcqqHtoWrdaWgaaWcbaGaam4Aaaqabaaa aa@522C@$

Also, the set of multi-target measurements can be written as

$Z k =[ ∪ x∈ X k Θ k (x) ]∪ K k MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGAbWaaSbaaSqaaiaadUgaaeqaaOGaeyypa0ZaamWaaeaadaWfqaqa aiabgQIiidWcbaGaaCiEaiabgIGiolaadIfadaWgaaadbaGaam4Aaa qabaaaleqaaOGaeuiMde1aaSbaaSqaaiaadUgaaeqaaOGaaiikaiaa hIhacaGGPaaacaGLBbGaayzxaaGaeyOkIGSaam4samaaBaaaleaaca WGRbaabeaaaaa@4C5B@$

where ${\Theta }_{k}\left(x\right)$ denotes the RFS of measurements generated by single-target with the state x at time k , and ${K}_{k}$ denotes the RFS of clutter measurements at time k .

The PHD (also called intensity function) of RFS is a nonnegative function v with the property that for any closed subset S

$E[ |X∩S| ]= ∫ S v(x)dx MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGfbWaamWaaeaacaGG8bGaamiwaiabgMIihlaadofacaGG8baacaGL BbGaayzxaaGaeyypa0Zaa8qeaeaacaWG2bGaaiikaiaahIhacaGGPa GaamizaiaahIhaaSqaaiaadofaaeqaniabgUIiYdaaaa@49A4@$ where $|X\cap S|$ denotes the cardinality of X on the space S .

To improve the accuracy of the estimation of target number, the CPHD filter jointly propagates the PHD of target states and the cardinality distribution. The CPHD recursion consists of the following prediction and update steps.

Prediction Step: Given the posterior PHD ${v}_{k-1}$ and posterior cardinality distribution ${p}_{k-1}$ at time $k-1$ , the predicted cardinality distribution ${p}_{k|k-1}$ and predicted PHD ${v}_{k|k-1}$ can be calculated as $p k|k−1 (n)= ∑ j=0 n p Γ,k (n−j) ∑ l=j ∞ C j l 〈 p S,k , v k−1 〉 j 〈 1− p S,k , v k−1 〉 l−j 〈 1, v k−1 〉 l p k−1 (l) MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaaca WGWbWaaSbaaSqaaiaadUgacaGG8bGaam4AaiabgkHiTiaaigdaaeqa aOGaaiikaiaad6gacaGGPaGaeyypa0ZaaabCaeaacaWGWbWaaSbaaS qaaiabfo5ahjaacYcacaWGRbaabeaakiaacIcacaWGUbGaeyOeI0Ia amOAaiaacMcadaaeWbqaaiaadoeadaqhaaWcbaGaamOAaaqaaiaadY gaaaGcdaWcaaqaamaaamaabaGaamiCamaaBaaaleaacaWGtbGaaiil aiaadUgaaeqaaOGaaiilaiaadAhadaWgaaWcbaGaam4AaiabgkHiTi aaigdaaeqaaaGccaGLPmIaayPkJaWaaWbaaSqabeaacaWGQbaaaOWa aaWaaeaacaaIXaGaeyOeI0IaamiCamaaBaaaleaacaWGtbGaaiilai aadUgaaeqaaOGaaiilaiaadAhadaWgaaWcbaGaam4AaiabgkHiTiaa igdaaeqaaaGccaGLPmIaayPkJaWaaWbaaSqabeaacaWGSbGaeyOeI0 IaamOAaaaaaOqaamaaamaabaGaaGymaiaacYcacaWG2bWaaSbaaSqa aiaadUgacqGHsislcaaIXaaabeaaaOGaayzkJiaawQYiamaaCaaale qabaGaamiBaaaaaaaabaGaamiBaiabg2da9iaadQgaaeaacqGHEisP a0GaeyyeIuoakiaadchadaWgaaWcbaGaam4AaiabgkHiTiaaigdaae qaaOGaaiikaiaadYgacaGGPaaaleaacaWGQbGaeyypa0JaaGimaaqa aiaad6gaa0GaeyyeIuoaaaa@8110@$ $v k|k−1 (x)= ∫ p S,k (ζ) f k|k−1 (x|ζ) v k−1 ( ζ)dζ+ r k (x) MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WG2bWaaSbaaSqaaiaadUgacaGG8bGaam4AaiabgkHiTiaaigdaaeqa aOGaaiikaiaahIhacaGGPaGaeyypa0Zaa8qaaeaacaWGWbWaaSbaaS qaaiaadofacaGGSaGaam4AaaqabaGccaGGOaGaaCOTdiaacMcacaWG MbWaaSbaaSqaaiaadUgacaGG8bGaam4AaiabgkHiTiaaigdaaeqaaO GaaiikaiaahIhacaGG8bGaaCOTdiaacMcacaWG2bWaaSbaaSqaaiaa dUgacqGHsislcaaIXaaabeaakiaacIcaaSqabeqaniabgUIiYdGcca WH2oGaaiykaiaadsgacaWH2oGaey4kaSIaamOCamaaBaaaleaacaWG RbaabeaakiaacIcacaWH4bGaaiykaaaa@61F7@$

where ${p}_{\Gamma ,k}\left(·\right)$ is the cardinality distribution of births at time k , ${C}_{j}^{l}$ is the binomial coefficient ( $l!/j!\left(l-j\right)!$ ), ${p}_{S,k}\left(\zeta \right)$ is the probability of target existence at time k given state $\zeta$ at time , $k-1$, $〈·,·〉$ is the inner product defined between two real-valued functions $\alpha$ and $\beta$ by $〈\alpha ,\beta 〉=\int \alpha \left(x\right)\beta \left(x\right)dx$ (or $〈\alpha ,\beta 〉=\sum _{l=0}^{\infty }\alpha \left(l\right)\beta \left(l\right)$ when $\alpha$ and $\beta$ are real sequences), ${f}_{k|k-1}\left(·|\zeta \right)$ is the single-target transition density at time k given state $\zeta$ at time $k-1$ , and ${r}_{k}\left(·\right)$ is the PHD of spontaneous births at time k .

Update Step: Given the predicted PHD ${v}_{k|k-1}$ and predicted cardinality distribution ${p}_{k|k-1}$ at time k, the updated cardinality distribution ${p}_{k}$ and updated PHD ${v}_{k}$ can be calculated as

${p}_{k}\left(n\right)=\frac{{\Upsilon }_{k}^{0}\left[{v}_{k|k-1},{Z}_{k}\right]\left(n\right){p}_{k|k-1}\left(n\right)}{〈{\Upsilon }_{k}^{0}\left[{v}_{k|k-1},{Z}_{k}\right],{p}_{k|k-1}〉}$
where ${\Upsilon }_{k}^{u}\left[v,Z\right]\left(n\right)=\sum _{j=0}^{\mathrm{min}\left(|Z|,n\right)}\left(|Z|-j\right)!{p}_{K,k}\left(|Z|-j\right){P}_{j+u}^{n}\frac{{〈1-{p}_{D,k},v〉}^{n-\left(j+u\right)}}{{〈1,v〉}^{n}}{e}_{j}\left({\Xi }_{k}\left(v,Z\right)\right)$
${\psi }_{k,z}\left(x\right)=\frac{〈1,{k}_{k}〉}{{k}_{k}\left(z\right)}{g}_{k}\left(z|x\right){p}_{D,k}\left(x\right)$
${\Xi }_{k}\left(v,Z\right)=\left\{〈v,{\psi }_{k,z}〉:z\in Z\right\}$

${Z}_{k}$ is the measurement set at time k,${p}_{K,k}\left(·\right)$ is the cardinality distribution of clutter at time k,${P}_{j}^{n}$ is the permutation coefficient ($n!/\left(n-j\right)!$ ), ${p}_{D,k}\left(x\right)$ is the probability of target detection at time k given the state x ${e}_{j}\left(·\right)$ is the elementary symmetric function of order j defined for a finite set Z of real numbers by ${e}_{j}\left(Z\right)=\sum _{S\subseteq Z,|S|=j}\left(\prod _{\varsigma \in S}\varsigma \right)$ with ${e}_{0}\left(Z\right)=1$, ${k}_{k}\left(·\right)$ is the PHD of clutter measurements at time k , and ${g}_{k}\left(·|x\right)$ is the single-target measurement likelihood at time k given the state x .

The GM-CPHD recursion is possible under the following assumptions:
1) Single-target state and measurement follow linear Gaussian models: $f k|k−1 (x|ζ)=N(x;Fζ, Q k−1 ) MathType@MTEF@5@5@+= feaagGart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGMbWaaSbaaSqaaiaadUgacaGG8bGaam4AaiabgkHiTiaaigdaaeqa aOGaaiikaiaahIhacaGG8bGaaCOTdiaacMcacqGH9aqpcaWGobGaai ikaiaahIhacaGG7aGaaCOraiaahA7acaGGSaGaaCyuamaaBaaaleaa caWGRbGaeyOeI0IaaGymaaqabaGccaGGPaaaaa@4E0D@$ $g k (z|x)=N(z;Hx, R k ) MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGNbWaaSbaaSqaaiaadUgaaeqaaOGaaiikaiaahQhacaGG8bGaaCiE aiaacMcacqGH9aqpcaWGobGaaiikaiaahQhacaGG7aGaaCisaiaahI hacaGGSaGaaCOuamaaBaaaleaacaWGRbaabeaakiaacMcaaaa@4854@$
where $N\left(·;m,P\right)$ denotes a Gaussian density with mean m and covariance P . According to (1) and (2), these assumptions are satisfied.

2) The detection probability and the survival probability are constant over the entire observed area: $P D,k (x)= P D,k MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGqbWaaSbaaSqaaiaadseacaGGSaGaam4AaaqabaGccaGGOaGaaCiE aiaacMcacqGH9aqpcaWGqbWaaSbaaSqaaiaadseacaGGSaGaam4Aaa qabaaaaa@42AC@$ $P S,k (x)= P S,k MathType@MTEF@5@5@+= feaagGart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXafv3ySLgzGmvETj2BSbqefm0B1jxALjhiov2D aebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaq Fr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qq Q8frFve9Fve9Ff0dmeaabaqaciGacaGaaeWabaWaaeaaeaaakeaaca WGqbWaaSbaaSqaaiaadofacaGGSaGaam4AaaqabaGccaGGOaGaaCiE aiaacMcacqGH9aqpcaWGqbWaaSbaaSqaaiaadofacaGGSaGaam4Aaa qabaaaaa@42CA@$

3) The birth intensity can be considered as Gaussian mixture:
${r}_{k}\left(x\right)=\sum _{i=1}^{{J}_{r,k}}{w}_{r,k}^{\left(i\right)}N\left(x;{m}_{r,k}^{\left(i\right)},{P}_{r,k}^{\left(i\right)}\right)$

where ${w}_{r,k}^{\left(i\right)}$ , ${m}_{r,k}^{\left(i\right)}$ , and ${P}_{r,k}^{\left(i\right)}$ are the weight, mean, and covariance of the birth Gaussians, and ${J}_{r,k}$ is their number.

The GM-CPHD recursion consists of the following prediction and update steps.

Prediction Step: At time $k-1$ , the posterior PHD ${v}_{k-1}$ and posterior cardinality distribution ${p}_{k-1}$ are given, and ${v}_{k-1}$ is a Gaussian mixture of the form

${v}_{k-1}\left(x\right)=\sum _{i=1}^{{J}_{k-1}}{w}_{k-1}^{\left(i\right)}N\left(x;{m}_{k-1}^{\left(i\right)},{P}_{k-1}^{\left(i\right)}\right)$

The predicted cardinality distribution ${p}_{k|k-1}$ and predicted PHD ${v}_{k|k-1}$ can be calculated as

${p}_{k|k-1}\left(n\right)=\sum _{j=0}^{n}{p}_{\Gamma ,k}\left(n-j\right)\sum _{l=j}^{\infty }{C}_{j}^{l}{p}_{k-1}\left(l\right){p}_{S,k}^{j}{\left(1-{p}_{S,k}\right)}^{l-j}$
${v}_{k|k-1}\left(x\right)={p}_{S,k}\sum _{j=1}^{{J}_{k-1}}{w}_{k-1}^{\left(j\right)}N\left(x;{m}_{S,k|k-1}^{\left(j\right)},{P}_{S,k|k-1}^{\left(j\right)}\right)+{r}_{k}\left(x\right)$

where           ${m}_{S,k|k-1}^{\left(j\right)}=F{m}_{k-1}^{\left(j\right)}$
${P}_{S,k|k-1}^{\left(j\right)}={Q}_{k-1}+F{P}_{k-1}^{\left(j\right)}{F}^{T}$

Update Step: At time k, the predicted PHD ${v}_{k|k-1}$ and predicted cardinality distribution ${p}_{k|k-1}$ are given, and ${v}_{k|k-1}$ is a Gaussian mixture of the form

${v}_{k|k-1}\left(x\right)=\sum _{i=1}^{{J}_{k|k-1}}{w}_{k|k-1}^{\left(i\right)}N\left(x;{m}_{k|k-1}^{\left(i\right)},{P}_{k|k-1}^{\left(i\right)}\right)$

The updated cardinality distribution ${p}_{k}$ and updated PHD ${v}_{k}$ can be calculated as

${p}_{k}\left(n\right)=\frac{{\Psi }_{k}^{0}\left[{w}_{k|k-1},{Z}_{k}\right]\left(n\right){p}_{k|k-1}\left(n\right)}{〈{\Psi }_{k}^{0}\left[{w}_{k|k-1},{Z}_{k}\right],{p}_{k|k-1}〉}$
${v}_{k}\left(x\right)=\frac{〈{\Psi }_{k}^{1}\left[{w}_{k|k-1},{Z}_{k}\right],{p}_{k|k-1}〉}{〈{\Psi }_{k}^{0}\left[{w}_{k|k-1},{Z}_{k}\right],{p}_{k|k-1}〉}\left(1-{p}_{D,k}\right){v}_{k|k-1}\left(x\right)+\sum _{z\in {Z}_{k}}\sum _{j=1}^{{J}_{k|k-1}}{w}_{k}^{\left(j\right)}\left(z\right)N\left(x;{m}_{k}^{\left(j\right)}\left(z\right),{P}_{k}^{\left(j\right)}\right)$

where

${\Psi }_{k}^{u}\left[w,Z\right]\left(n\right)=\sum _{j=0}^{\mathrm{min}\left(|Z|,n\right)}\left(|Z|-j\right)!{p}_{K,k}\left(|Z|-j\right){P}_{j+u}^{n}\frac{{\left(1-{p}_{D,k}\right)}^{n-\left(j+u\right)}}{{〈1,w〉}^{j+u}}{e}_{j}\left({\Lambda }_{k}\left(w,Z\right)\right)$
${\Lambda }_{k}\left(w,Z\right)=\left\{\frac{〈1,{k}_{k}〉}{{k}_{k}\left(z\right)}{p}_{D,k}{w}^{T}{q}_{k}\left(z\right):z\in Z\right\}$
${w}_{k|k-1}={\left[{w}_{k|k-1}^{\left(1\right)},\cdots ,{w}_{k|k-1}^{\left({J}_{k|k-1}\right)}\right]}^{T}$
${q}_{k}\left(z\right)={\left[{q}_{k}^{\left(1\right)}\left(z\right),\cdots ,{q}_{k}^{\left({J}_{k|k-1}\right)}\left(z\right)\right]}^{T}$
${q}_{k}^{\left(j\right)}\left(z\right)=N\left(z;{\eta }_{k|k-1}^{\left(j\right)},{S}_{k|k-1}^{\left(j\right)}\right)$
${\eta }_{k|k-1}^{\left(j\right)}=H{m}_{k|k-1}^{\left(j\right)}$
${S}_{k|k-1}^{\left(j\right)}=H{P}_{k|k-1}^{\left(j\right)}{H}^{T}+{R}_{k}$
${w}_{k}^{\left(j\right)}\left(z\right)={p}_{D,k}{w}_{k|k-1}^{\left(j\right)}{q}_{k}^{\left(j\right)}\left(z\right)\frac{〈{\Psi }_{k}^{1}\left[{w}_{k|k-1},{Z}_{k}\\left\{z\right\}\right],{p}_{k|k-1}〉〈1,{\kappa }_{k}〉}{〈{\Psi }_{k}^{0}\left[{w}_{k|k-1},{Z}_{k}\right],{p}_{k|k-1}〉{\kappa }_{k}\left(z\right)}$
${m}_{k}^{\left(j\right)}\left(z\right)={m}_{k|k-1}^{\left(j\right)}+{K}_{k}^{\left(j\right)}\left(z-{\eta }_{k|k-1}^{\left(j\right)}\right)$
${P}_{k}^{\left(j\right)}=\left(I-{K}_{k}^{\left(j\right)}H\right){P}_{k|k-1}^{\left(j\right)}$
${K}_{k}^{\left(j\right)}={P}_{k|k-1}^{\left(j\right)}{H}^{T}{\left[{S}_{k|k-1}^{\left(j\right)}\right]}^{-1}$

According to GM-CPHD recursion, the number of Gaussian components increases dramatically with time. To reduce the computation time, the 'pruning' and 'merging' procedure11-13 can be used for GM-CPHD. The components with negligible weights can be discarded and the closed spaced components are merged into one as they are more efficiently approximated by a single Gaussian term.
The number of targets can be estimated using MAP estimator

${\stackrel{^}{N}}_{k}=\mathrm{arg}\mathrm{max}{p}_{k}\left(·\right)$

and the state estimates can be extracted by picking the means of ${\stackrel{^}{N}}_{k}$ Gaussian terms (from posterior PHD Vk ) with the largest weights.

Although the GM-CPHD filter can estimate the multi-target states at each time step, the track of individual targets is not produced. Here, we construct 0-1 integer programming to associate the estimated states between frames for track continuity.

Let $\left\{{\stackrel{^}{x}}_{k-1,1},\cdots ,{\stackrel{^}{x}}_{k-1,\stackrel{^}{N}\left(k-1\right)}\right\}$ be the estimated states at time $k-1$ , and $\left\{{P}_{k-1}^{\left(1\right)},\cdots ,{P}_{k-1}^{\left(\stackrel{^}{N}\left(k-1\right)\right)}\right\}$ the covariance matrixes of corresponding Gaussian terms of posterior PHD ${v}_{k-1}$ . In the same way, $\left\{{\stackrel{^}{x}}_{k,1},\cdots ,{\stackrel{^}{x}}_{k,\stackrel{^}{N}\left(k\right)}\right\}$ and $\left\{{P}_{k}^{\left(1\right)},\cdots ,{P}_{k}^{\left(\stackrel{^}{N}\left(k\right)\right)}\right\}$ are defined for time k . For and , we define scalar quantity ${d}_{ij}$ as follows

${d}_{ij}={\left({\stackrel{^}{x}}_{k,i}-F{\stackrel{^}{x}}_{k-1,j}\right)}^{T}{\left({P}_{k}^{\left(i\right)}+F{P}_{k-1}^{\left(j\right)}{F}^{T}\right)}^{-1}\left({\stackrel{^}{x}}_{k,i}-F{\stackrel{^}{x}}_{k-1,j}\right)$

If $\stackrel{^}{N}\left(k-1\right)\le \stackrel{^}{N}\left(k\right)$ , new targets may appear at time k . We construct the following 0-1 integer programming

$\mathrm{min}\sum _{i=1}^{\stackrel{^}{N}\left(k\right)}\sum _{j=1}^{\stackrel{^}{N}\left(k-1\right)}{y}_{ij}{d}_{ij}$

where ${y}_{ij}$ is a binary variable. If ${y}_{ij}=1$ , ${\stackrel{^}{x}}_{k,i}$ and ${\stackrel{^}{x}}_{k-1,j}$ are associated. If ${y}_{ij}=0$ , and ${\stackrel{^}{x}}_{k-1,j}$ belong to different targets.

We solve this 0-1 integer programming using the branch-and-bound algorithm17. If the solutions $\left\{{y}_{ij}\right\}$ satisfy $\sum _{j=1}^{\stackrel{^}{N}\left(k-1\right)}{y}_{ij}=0$ , we set ${\stackrel{^}{x}}_{k,i}$ as the initial state of a spontaneous new target at time k . If $\stackrel{^}{N}\left(k-1\right)>\stackrel{^}{N}\left(k\right)$ , some targets disappear at time k . We construct the following 0-1 integer programming

$\mathrm{min}\sum _{i=1}^{\stackrel{^}{N}\left(k\right)}\sum _{j=1}^{\stackrel{^}{N}\left(k-1\right)}{y}_{ij}{d}_{ij}$

If the solutions $\left\{{y}_{ij}\right\}$ satisfy $\sum _{i=1}^{\stackrel{^}{N}\left(k\right)}{y}_{ij}=0$ , the target with the state ${\stackrel{^}{x}}_{k-1,j}$ disappears at time k. We end the track of this target.

In this section, we demonstrate the performance of tracking ballistic target group on space-based infrared focal plane using GM-CPHD filter. The ballistic target group consists of one warhead and eight decoys. An early warning satellite continually tracks this group for duration of 800s. Due to the finite resolution of infrared sensor, four decoys appear on infrared focal plane at 87s, 146s, 194s and 665s respectively, two decoys synchronously appear at 278s and other two decoys synchronously appear at 547s. All the targets do not disappear. The infrared focal plane is the square $\chi =\left[-64,64\right]×\left[-64,64\right]$ (in pixel). The sampling period $T=1$ s. The standard deviation of the measurement noise is 1pixel. The detection probability ${P}_{D,k}$ means that each target is either detected with probability ${P}_{D,k}$ or missed with probability $1-{P}_{D,k}$ at time k . We set ${P}_{D,k}=0.99$ for this simulation experiment. The survival probability ${P}_{S,k}$ means that each target at time k - 1 either continues to exist at time k with probability ${P}_{S,k}$ or dies with the probability $1-{P}_{S,k}$ . Since there are no target deaths, the survival probability is ${P}_{S,k}=1$ . The clutter is modeled as a Poisson RFS with intensity function $k\left(z\right)={\lambda }_{c}{f}_{c}\left(z\right)$ , where ${f}_{c}\left(z\right)$ represents the uniform probability density over X, and ${\lambda }_{c}=3$ is the average number of clutter per frame.
We use Wasserstein distance (WD)18 to evaluate the accuracy of multi-target state estimates. Let $X=\left\{{x}_{1},\cdots ,{x}_{|X|}\right\}$ and $\stackrel{^}{X}=\left\{{\stackrel{^}{x}}_{1},\cdots ,{\stackrel{^}{x}}_{|\stackrel{^}{X}|}\right\}$ be the true and estimated multi-target states. The WD between $X$ and $\stackrel{^}{X}$ is defined by

$d\left(X,\stackrel{^}{X}\right)=\underset{C}{\mathrm{min}}{\left[\sum _{i=1}^{|\stackrel{^}{X}|}\sum _{j=1}^{|X|}{C}_{ij}|{\stackrel{^}{x}}_{i}-{x}_{j}{|}^{2}\right]}^{1/2}$

where the minimum is taken over the set of all transportation matrices $C=\left\{{C}_{ij}\right\}$ , and each entry of the matrix satisfies ${C}_{ij}\ge 0$ , $\sum _{i=1}^{|\stackrel{^}{X}|}{C}_{ij}=1/|X|$ , and $\sum _{j=1}^{|X|}{C}_{ij}=1/|\stackrel{^}{X}|$ .

Fig. 1. shows true tracks together with measurements for the duration of 800s on the focal plane. The warhead stays at the origin, and the decoys appear around the origin and move radially outwards.
Fig. 2 shows that the estimated tracks for individual targets on focal plane. We can see that the GM-CPHD filter can efficiently eliminates the clutter, and that the estimates of states approximate to the real states. Further more, the track continuity for each target is successfully achieved.

Figure 1.Measurements and true tracks on focal plane.

Figure 2. Estimated tracks for individual targets on focal plane.

1000 Monte Carlo runs are implemented on the same target trajectories but with independently generated measurements for each trial. Fig. 3. and Fig. 4 show that the mean of estimated target number versus time for GM-PHD and GM-CPHD filters respectively. We can see that both filters are unbiased estimation for target number. The standard deviation of estimated target number versus time for GM-PHD and GM-CPHD filters is shown in Fig. 5 The average variance of estimated target number of GM-CPHD filter is reduced by 75%, compared with the GM-PHD filter. It can be seen that the GM-CPHD filter is much more reliable and accurate for the estimation of target number than GM-PHD filter. The plots also show that the accuracy of estimation of target number dramatically drops (exhibit high peaks in Fig. 5.) for both filters when the target number changes.

Figure 3. Mean of estimated target number versus time for GM-PHD filter.

Figure 4. Mean of estimated target number versus time for GM-CPHD filter.

Figure 6 shows that the WD versus time for GM-PHD and GM-CPHD filters. When the target number changes, the WD exhibits high peaks for both filters. The WD penalizes errors in both the target localization and the target number estimates. The GM-CPHD and GM-PHD filters have the similar error of target number estimates when the target number changes (see Fig. 5). The WD peaks of GM-CPHD filter are higher than GM-PHD filter since the GM-CPHD filter has more error of target localization at these time steps. When the target number is steady, the WD of GM-CPHD filter is about 0.5pixel (much smaller than GM-PHD filter). The average WD of GM-CPHD filters is reduced by 43.39%, compared with GM-PHD filter. It can be seen that the GM-CPHD filter acquires much more accurate estimation of multi-target states than GM-PHD filter. From Fig. 6, we can see that GM-PHD filter has an increasing trend of WD. This phenomenon is due to the dramatic increase of cardinality error as time progressing (see Fig. 5).

Figure 5. Comparison of standard deviation of estimated target number versus time between GM-PHD and GM-CPHD

Figure 6. Comparison of WD versus time between GM-PHD and GM-CPHD

This paper proposes GM-CPHD filter for tracking of midcourse ballistic target group on space-based infrared focal plane. In order to improve the accuracy of the estimation of target number, the GM-CPHD filter jointly propagates the intensity function of target states and the cardinality distribution. 0-1 integer programming is constructed to associate the estimated states between frames for track continuity. The simulation results show that the target number estimation in GM-CPHD filter is unbiased, and that the variance of estimated target number and the errors of estimated states in GM-CPHD filter are reduced by 75% and 43.39% respectively, compared with the GM-PHD filter. Moreover, the simulation suggests that the track continuity for each target is successfully achieved.

1. Korn, J.; Holtz, H. & Farber, M.S. Trajectory estimation of closely spaced objects (CSO) using infrared focal plane data of an STSS (Space Tracking and Surveillance System) platform. SPIE, 2004, 5428, 387-399.

2. Rago, C. & Landau, H. Stereo spatial super-resolution technique for multiple reentry vehicles. In Proceedings of 2004 IEEE Aerospace Conference, Big Sky, Montana, USA, March 6-13, pp. 1834-1841.

3. Lin, L.; Xu, D.; Sheng, W.; An, W. & Xu, H. Tracking of midcourse ballistic target group with space-based infrared FPA based on random finite set. J. Infrared Millim. W., 2010, 29(6), 465-470(Chinese).

4. Fortmann, T.E.; Bar-Shalom, Y. & Scheffe, M. Sonar tracking of multiple targets using joint probabilistic data association. IEEE J. Oceanic Eng., 1983, 8(3), 173-184.

5. Reid, D.B. An algorithm for tracking multiple targets. IEEE T. Automat. Contr., 1979, 24(6), 843-854.

6. Mahler, R. Multitarget Bayes filtering via first-order multitarget moments. IEEE T. Aero. Elec. Sys., 2003, 39(4), 1152-1178.

7. Vo, B.-N.; Singh, S. & Doucet, A. Sequential Monte Carlo methods for multi-target filtering with random finite sets. IEEE T. Aero. Elec. Sys., 2005, 41(4), 1224-1245.

8. Johansen, A.M.; Singh, S.S.; Doucet, A. & Vo, B.-N. Convergence of the SMC implementation of the PHD filter. Methodol. Comput. Appl. Probab., 2006, 8(2), 265-291.

9. Panta, K.; Vo, B.-N. & Singh, S. Novel data association schemes for the probability hypothesis density filter. IEEE T. Aero. Elec. Sys., 2007, 43(2), 556-570.

10. Clark, D.E. & Bell, J. Multi-target state estimation and track continuity for the particle PHD filter. IEEE T. Aero. Elec. Sys., 2007, 43(4), 1441-1453.

11. Vo, B.-N. & Ma, W.-K. The Gaussian mixture probability hypothesis density filter. IEEE T. Signal Proces., 2006, 54(11), 4091-4104.

12. Clark, D. & Vo, B.-N. Convergence analysis of the Gaussian mixture PHD filter. IEEE T. Signal Proces., 2007, 55(4), 1204-1212.

13. Panta, K.; Clark, D.E. & Vo, B.-N. Data association and track management for the Gaussian mixture probability hypothesis density filter. IEEE T. Aero. Elec. Sys., 2009, 45(3), 1003-1016.

14. Pasha, S.A.; Vo, B.-N.; Tuan, H.D. & Ma, W.-K. A Gaussian mixture PHD filter for jump Markov system models. IEEE T. Aero. Elec. Sys., 2009, 45(3), 919-936.

15. Mahler, R. PHD filters of higher order in target number. IEEE T. Aero. Elec. Sys., 2007, 43(4), 1523-1543.

16. Vo, B.-T.; Vo, B.-N. & Cantoni, A. Analytic implementations of the cardinalized probability hypothesis density filter. IEEE T. Signal Proces., 2007, 55(7): 3553-3567.

17. Hillier, F.S. & Lieberman, G.J. Introduction to operations research. McGraw-Hill Press, New York, US, 2001, pp.604-16.

18. Hoffman, J.R. & Mahler R. Multitarget miss distance via optimal assignment. IEEE Trans. Syst. Man CY. A, 2004,34(3), 327-36.

 Mr Dong Li received BS and MS degrees from National University of Defense Technology in 2006 and 2008 respectively. His research interests are target tracking and signal processing. Mr Dong-Yun Yi received PhD degree from National University of Defense Technology in 2003. Since then he has been a professor at the Department of Mathematics and System Science, National University of Defense Technology. His main research interests are signal processing, data fusion and dynamic system analysis.