Improved Bearings-Only Multi-Target Tracking with GM-PHD Filtering

18 Pages·2016·0.93 MB·English

sensors
Article
Improved Bearings-Only Multi-Target Tracking with
GM-PHD Filtering
QianZhangandTaekLyulSong∗
DepartmentofElectronicSystemsEngineering,HanyangUniversity,Ansan,Gyeonggi-do15588,Korea;
[email protected]
*Correspondence:[email protected];Tel.:[phone]
AcademicEditor:VittorioM.N.Passaro
Received:25July2016;Accepted:7September2016;Published:10September2016
Abstract: In this paper, an improved nonlinear Gaussian mixture probability hypothesis density
(GM-PHD) filter is proposed to address bearings-only measurements in multi-target tracking.
Theproposedmethod,calledtheGaussianmixturemeasurements-probabilityhypothesisdensity
(GMM-PHD)filter,notonlyapproximatestheposteriorintensityusingaGaussianmixture,butalso
modelsthelikelihoodfunctionwithaGaussianmixtureinsteadofasingleGaussiandistribution.
Besides,thetargetbirthmodeloftheGMM-PHDfilterisassumedtobepartiallyuniforminsteadof
aGaussianmixture. SimulationresultsshowthattheproposedfilteroutperformstheGM-PHDfilter
embeddedwiththeextendedKalmanfilter(EKF)andtheunscentedKalmanfilter(UKF).
Keywords: nonlinear estimation; bearings-only measurement; multi-target tracking; Gaussian
mixturemeasurements;passivesensor
1. Introduction
Bearings-onlymulti-targettracking(MTT)[1–3]withclutterandmisseddetectionsisachallenging
nonlinear problem. The filter for such a problem should not only deal with the nonlinearity of
measurements,butalsosolvethemeasurementoriginuncertaintyofmultipletargets. Bearings-only
measurements [4] are often obtained by passive sensors, such as sonars, and contain a low level
of information about detected targets, leading to the observability problem [5]. For bearings-only
single-target tracking, the observable level can be evaluated by the Cramer–Rao lower bound
(CRLB) [6,7]. However, it is difficult to evaluate the level of observability in a theoretical context
forthebearings-onlymulti-targettrackingwithclutterandmisseddetections[8].Inordertoovercome
the non-observable problem, the sensor needs to outmaneuver the targets [5,9]. As the origin of
a measurement is not known for target tracking in clutter, the tracking filter needs to have a
proper track-to-measurement-association method [1] or something equivalent, such as symmetric
measurementequations[10]andrandomfinitesets[11].
Fortrackingtargetsusingnonlinearmeasurements,themostpopularmethodistheextended
Kalman filter (EKF) [12]. It linearizes the nonlinear measurement function around the predicted
target state under the assumption that the true target state is close enough to the predicted state.
Then, the measurement update of the filter can be performed using the Kalman filter (KF) [12],
whichhasagoodperformanceforlinearsystems. TheunscentedKalmanfilter(UKF)[13]appliesthe
unscentedtransforminsteadoflinearizationtothenonlinearmeasurementfunction. Theprobability
density functions (pdfs) are represented and propagated using several sigma points. The particle
filter(PF)[14]representsthenonlinearpdfsbyconsideringtheamountofparticlesandcanobtain
betterperformancesthantheEKFandtheUKF.However,particlefiltershaveaheavycomputational
load. Apartfromthesenonlinearfilters,manyotherfiltersexistfornonlinearsystems,suchasthe
Sensors2016,16,1469;doi:10.3390/s16091469 www.mdpi.com/journal/sensors
Sensors2016,16,1469 2of18
cubatureKalmanfilter(CKF)[15],theshiftedRayleighfilter(SRF)[16]andtheensembleKalmanfilter
(EnKF)[17].
TohandlethemeasurementoriginuncertaintyprobleminMTTwithclutters,manytechniques
havebeendeveloped. Manyofthesemethodsbelongtothefollowingcategories: jointprobabilitydata
association(JPDA)[18],multiplehypothesistracking(MHT)[19,20]andrandomfiniteset(RFS)[11]
basedmethods. TheJPDAfilterscomefromcombingtheprobabilitydataassociation(PDA)[1]filter
withjointeventsandcanonlytrackfixedandknownnumbersoftargets. Toaccommodatevariedand
unknownnumbersoftargets,thejointintegratedPDA(JIPDA)[21]filter,whichisabletoestimatethe
probabilityoftargetexistence,isproposed. Toimprovethetrackingaccuracy,amulti-scanmulti-target
trackingalgorithm,thejointintegratedtracksplitting(JITS),wasdevelopedin[22]. AstheJIPDAand
theJITSsufferfromheavycomputationalloadwhentrackingtargetsinmutualproximity,thelinear
multi-target IPDA (LMIPDA) and the linear multi-target ITS (LMITS) were proposed in [22,23],
respectively. Besides, the iterative JIPDA (iJIPDA) [24] tracker is also a computationally-efficient
algorithm for the MTT. The MHT filters attempt to maintain and evaluate a set of measurement
hypotheseswithhightrackscores. TherearemanyversionsoftheMHTfilter,andmostofthemcan
begroupedintotwoclasses: thetrack-orientedclass[25]andthemeasurement-orientedclass[19].
TheJPDAandMHTapproachesareformulatedviadataassociation,whichtakesalargeproportionof
computationalresources,whiletheRFSapproachisanemergingparadigmandisestablishedwithout
dataassociation. TheRFS-basedfilters[26–28]treatmulti-targetstatesandmeasurementsasthestate
finitesetandthemeasurementfiniteset,respectively. TheRFS-basedfilterstrytoestimatethetarget
setbasedonthemeasurementset. Inthisway,theRFS-basedfiltersareperformedinacomputational
efficientwayforthemulti-targettracking.
AmongRFS-basedfilters,theprobabilityhypothesisdensity(PHD)[26]filterisafirstmoment
approximationtothemulti-targetpredictedandposteriordensities. Itpropagatesthetargetintensities
withoutconsideringdataassociationsbetweentargetsandmeasurements. Asthereisnocloseformto
thePHDfilter,asequentialMonteCarlo(SMC)orGaussianmixture(GM)techniquewasimplemented,
whichresultedintheSMC-PHD[29,30]filterortheGM-PHD[31]filter. Forbearings-onlyMTT,itis
difficult to apply the SMC-PHD filter because of the problem of extracting estimated target states.
Forthisreason,theGM-PHDfilterisconsideredinthispaper. TheGM-PHDfilterassociatedwiththe
EKFandtheUKFaretermedasGM-PHD-EKFandGM-PHD-UKF,respectively.
IntheGM-PHD-EKFandtheGM-PHD-UKF,thepredictedintensityismodeledbyaGaussian
mixture, and the likelihood function [12] is approximated using a single Gaussian distribution.
The updated intensity of each filter is obtained after the predicted intensity is updated by the
measurements.Asbearings-onlymeasurementssufferfromseverenonlinearity,thelikelihoodfunction
approximatedbyasingleGaussianisnotaccurateenough. Furthermore,theaccuracyoftheupdated
intensitycannotgetenoughimprovementsafterthepredictedintensityisupdatedbytheinaccurate
likelihoodfunctionofthemeasurements. Inthispaper,theGaussianmixturemeasurements-PHD
(GMM-PHD) filter is proposed to address bearings-only MTT with clutter and missed detections.
To improve the tracking accuracy, the target intensities are approximated by Gaussian mixtures,
andthelikelihoodfunctionisalsomodeledbyaGaussianmixtureintheGMM-PHDfilter. Inthis
way,theupdatedintensityoftheGMM-PHDismuchmoreaccuratethanthoseoftheGM-PHD-EKF
andtheGM-PHD-UKFafterthepredictedintensityisupdatedbyamoreaccuratelikelihoodfunction,
which is modeled by Gaussian mixtures. The proposed filter is a nonlinear MTT algorithm that
can address not only bearings-only measurements, but also other nonlinear measurements. For
bearings-onlyMTT,thetargetsmayappearatanyplaceinthemeasurementspace. Inordertoreduce
thenumberofselectedparametersandavoidtheundesirableeffectofpoorparameterselection[8],
thederivationoftheGMM-PHDfilterisprocessedwithapartially-uniformtargetbirthmodel[8,32],
but not a Gaussian mixture birth model. In the simulation experiment, the performance of the
GMM-PHD is compared with the GM-PHD-EKF and the GM-PHD-UKF in terms of the optimal
Sensors2016,16,1469 3of18
subpattern assignment (OSPA) [33] distance, the OSPA localization, the OSPA cardinality and the
averageCPUtime.
Therestofthepaperisarrangedasfollows. Thebearings-onlyMTTproblemispresentedin
Section 2. The proposed GMM-PHD filter is derived in Section 3. A simulation study is given in
Section4. Finally,aconclusionisproposedinSection5.
2. TheBearings-OnlyMTTProblem
Inthispaper,bearings-onlyMTTwithonepassivesensoronamaneuveringplatformisstudied.
Thetargetsobeythecontinuouswhitenoiseacceleration(CWNA)motionmodel[1,12]. Thetarget
motionmodelandthesensormeasurementmodelinthe2DCartesiancoordinatescaseareconsidered
inthissection. Forthetargett,thetargetstateet withposition(xt,yt)andvelocity(x˙t,y˙t)attimekis
k k k k k
expressedas:
et = [xt,yt,x˙t,y˙t](cid:48). (1)
k k k k k
AnRFSX oftargetstatescancontainanynumberN oftargetsattimekandisgivenby:
k k
X = {et ,...,et }. (2)
k k,1 k,Nk
Let χ denote the single target state space and F(χ) denote the set of all finite subsets of χ,
i.e.,et ,...,et ∈ χandX ∈ F(χ).
k,1 k,Nk k
As the target t is assumed to follow the CWNA motion model, its dynamic model can be
expressedas:
etk = Φetk−1+νk−1, (3)
wherethestatepropagationmatrixΦistime-invariant,
(cid:34) (cid:35)
1 T
Φ = ⊗I , (4)
2
0 1
νk−1isasequenceofzeromean,whiteGaussianprocessnoiseswithcovariance:
(cid:34) (cid:35)
T4/4 T3/2
Q = q ⊗I , (5)
k−1 T3/2 T2 2
Tisthesamplingtime,I isthe2×2identitymatrixandqisthepowerspectraldensity(PSD)[1].
2
Thesensorstatees isgivenas:
k
es = [xs,ys,x˙s,y˙s](cid:48). (6)
k k k k k
Each target t can be detected with the probability P at time k. The sensor can generate a
D,k
measurementwhenthetargetisdetected. Weuseθt todenotethetargetmeasurementattimekfor
k
targett. Themeasurementfromtheradaris:
∆
w = θt = h(et,es)+(cid:118) , (7)
k k k k k
where:
(cid:32) (cid:33)
xt −xs
h(et,es) =tan−1 k k , (8)
k k yt −ys
k k
(cid:118) iszeromean,whiteGaussianmeasurementnoisewithcovariance:
k
R = σ2, (9)
k θ
Sensors2016,16,1469 4of18
thatitisuncorrelatedwithν . Therefore,eachtargettcangenerateameasurementRFSΘ (et),which
k k k
canbeeither{θt}(thetargettisdetectable)or∅(thetargettisnotdetected). Inaddition,thesensor
k
alsoproducesfalsemeasurements,whichformanRFSK ateachtimek. Then,theMTTmeasurement
k
RFSW attimekcanbeexpressedas:
k
 
Wk =Kk∪ (cid:91) Θk(e). (10)
e∈Xk
3. TheProposedGMM-PHDFilter
InRFS-basedmethods,theBayesianrecursionofmulti-targetposteriordensityispropagatedin
timeas:
(cid:90)
pk|k−1(Xk|W1:k−1)= fk|k−1(Xk|X)pk−1(X|W1:k−1)µs(dX), (11)
g (W |X )p (X |W )
pk(Xk|W1:k) = (cid:82) gk(kWk|Xk)pkk|k−k1|k(−X1|Wk1:k−11):k−µs1(dX), (12)
where f (X |X)andg (W |X )denotethemulti-targettransitiondensityandlikelihoodfunction,
k|k−1 k k k k
respectively.Here,p (X |W )denotesthemulti-targetposterior,densityandµ (dX)isanappropriate
k k 1:k s
referencemeasureonF(χ)[31].
Toreducethecomputationalintractability,afirstmomentapproximation(thePHDfilter)ofthe
recursionisproposedin[26]. ThegeneralformofthePHDfilterrecursion(withouttargetspawning)
isgivenby:
(cid:90)
vk|k−1(e)= PS,k(η) fk|k−1(e|η)vk−1(η)dη+γk(e), (13)
P (e)g (w|e)v (e)
vk|k(e) = [1−PD,k(e)]vk|k−1(e)+ ∑ κ (w)+D(cid:82),kP (ke)g (w|ke|k)−v1 (e)de, (14)
w∈Wk k D,k k k|k−1
where v (e) and v (e) denote the predicted intensity from time k−1 to k and the updated
k|k−1 k|k
intensity,respectively,P (η)isthesurvivalprobability,meaningthetargetstillsurvivesattimek,
S,k
f (e|η)representsthesingletargettransitiondensityfromtimek−1tok,γ (e)denotestheprior
k|k−1 k
intensityofspontaneoustargetbirthsattimek,g (w|e)isthesingletargetmeasurementlikelihood
k
functionandκ (w)istheclutterintensity.
k
3.1. IntensityPrediction
According to the Gaussian mixture assumption, the posterior intensity at time k−1 can be
expressedas:
vk−1(e) = J∑k−1ωk(i−)1N (cid:16)e;eˆ(ki−)1,J(ki−)1(cid:17). (15)
i=1
Letθ,r,c,andsdenotethebearing,range,courseandspeedinpolarcoordinates,respectively.
ThefunctionΨ(e;γ (θ,r,c,s))representsthetransformationofthedensityγ (θ,r,c,s)frompolar
k k
coordinatestoCartesiancoordinates. Then,thepredictedintensityattimekisgivenby:
v (e) = v (e)+Ψ(e;γ (θ,r,c,s)), (16)
k|k−1 S,k|k−1 k
where:
v (e) = P J∑k−1ω(i) N (cid:16)e; eˆ(i) ,J(i) (cid:17), (17)
S,k|k−1 S,k k−1 k|k−1 k|k−1
i=1
(cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17)
γ (θ,r,c,s) = ωbU(θ;H )N r;r¯,σ2 N c;θ−π,σ2 N s;s¯,σ2 . (18)
k k θ r c s
Sensors2016,16,1469 5of18
In Equation (17), the survived probability P (e) is assumed to be independent of the target
S,k
stateandisgivenasP . InEquation(18),ωb istheexpectednumberoftargetsappearingattimek,
S,k k
U(θ;H )istheuniformdistributioninθovertheregionH ,r¯isthepriormeanoftherangeandσ2
θ θ r
isitspriorvariance. Similarly,s¯isthepriormeanofthespeedassociatedwithitspriorvarianceσ2,
s
and σ2 is the prior course variance. The predicted estimates eˆ(i) and J(i) are obtained via the
c k|k−1 k|k−1
Kalmanprediction:
eˆ(i) = Φeˆ(i) , (19)
k|k−1 k−1
J(i) = ΦJ(i) Φ(cid:48) +Q . (20)
k|k−1 k|k−1 k−1
Thetargetstateisaugmentedbyabinaryvariableβ,sowemaydistinguishbetweenthesurviving
componentsandthebirthcomponents:

 Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17),β =0
v (e,β) = k|k−1 k|k−1 k|k−1 (21)
k|k−1 i=1
 Ψ(e;γ (θ,r,c,s)),β =1
k
whereω(i) = P ω(i) . Asnotedin[34],thesurvivingandbirthcomponentsareseparatedtoavoid
k|k−1 S,k k−1
biasingthecardinalityestimates.
3.2. GMMLikelihoodApproximation
In GM-PHD filters, the likelihood function g (w|e) is always modeled as a single Gaussian
k
distribution. However,itisapproximatedbyaGaussianmixtureintheproposedGMM-PHDfilter.
Though the measurement noise is assumed to be Gaussian, the measurement uncertainty is
non-Gaussian (non-ellipse) in Cartesian coordinates. The GMM measurement presentation first
divides the non-elliptical measurement uncertainty into several segments. Then, each segment is
modeledbyoneGaussiandistribution,andthewholemeasurementuncertaintycanbeapproximated
byaGaussianmixture.
Suppose the measurement uncertainty in Cartesian coordinates is determined by the range
interval [r ,r ] and the measurement θt with standard deviation σ . The range interval is
k,min k,max k θ
dividedinto A subintervals,givenby[35,36]:
k
r
k,a+1 = τ ; a =1,...,A , (22)
r k k
k,a
where:
(cid:18)r (cid:19)1/Ak
τ = k,max . (23)
k r
k,min
Then, each segment a can be determined by rk,a,rk,a+1,θkt −σθ,θkt +σθ in polar coordinates.
Letr¯ = (rk,a+rk,a+1)/2and∆r = (rk,a+1−rk,a)/2,thenthesegmentaisapproximatedbyaGaussian
distributionwithmeanwˆ andcovarianceR inCartesiancoordinates:
k,a k,a
(cid:34) xs (cid:35) (cid:34) sin(cid:0)θt(cid:1) (cid:35)
wˆk,a = ysk +r¯ cos(cid:0)θkt(cid:1) , (24)
k k
(cid:34) (cid:35)
R = φ (∆r)2 0 φ(cid:48), (25)
k,a 0 r¯2σ2
θ
where:
(cid:34) sin(cid:0)θt(cid:1) −cos(cid:0)θt(cid:1) (cid:35)
φ = cos(cid:0)θkt(cid:1) sin(cid:0)θt(cid:1)k . (26)
k k
Sensors2016,16,1469 6of18
Obviously,theareaofeachsegmentisdifferent,andtheweightofthesegmentaisproportional
toitsarea,givenby[35]:
(cid:112)
det(R )
λ = k,a , (27)
k,a Ak (cid:112)
∑ det(R )
k,a
a=1
and
∑Ak
λ =1. (28)
k,a
a=1
Then,themeasurementlikelihoodfunctioninCartesiancoordinate(β =0)isapproximatedasa
Gaussianmixture:
∑Ak
g (w|e,β =0) ≈C λ N(wˆ ;He ,R ), (29)
k k k,a k,a k k,a
a=1
(cid:34) (cid:35)
1 0 0 0
where,H= istheobservationmatrixandtheconstantC iscalculatedas[37]:
k
0 1 0 0
C = (cid:90) rk,maxrdr = rk2,max−rk2,min. (30)
k
rk,min 2
SincethemeasurementnoiseismodeledasGaussian,thelikelihoodfunctioninpolarcoordinates
(β =1)isexpressedas:
(cid:16) (cid:17)
g (w|e,β =1) = N θt;θ,σ2 (31)
k k θ
Figure1givesanexampleoftheGMMpresentation. Inthefigure,themeasurementuncertainty
oftheGMM-PHDisapproximatedbysixmeasurementcomponents(shownasthesolidellipses),and
thatoftheGM-PHD-EKFismodeledbyasingleGaussiandistribution(shownasthedashedellipse).
Actually,themeasurementuncertaintyoftheGM-PHD-UKFisalsoapproximatedbyoneGaussian
distribution. Thetruetargetisdisplayedbyacross. Itisobviousthatthemeasurementlikelihood
functionapproximatedusingGaussianmixturesintheGMM-PHDismuchmoreaccuratethanthe
approximationwithoneGaussianintheGM-PHD-EKF.Thus,theGMM-PHDcanperformwitha
muchhighertrackingaccuracythantheGM-PHD-EKFandtheGM-PHD-UKF.
r
k,max
s
q
r
k,min
GM-PHD-EKF
GMM-PHD
sensor
Figure1.AnexampleofGMMpresentation.
Sensors2016,16,1469 7of18
3.3. IntensityUpdate
New targets are assumed to be always detected at their time of birth. Furthermore, we also
assumethatthedetectionprobabilityforasurvivingtargetisindependentoftheirstateand:
(cid:40)
P ,β =0
P (e,β) = D,k . (32)
D,k 1,β =1
AccordingtoEquation(14),theposteriorintensityattimekisgivenby:
v (e,β =0)
k|k
= [1−P ]Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17)
D,k k|k−1 k|k−1 k|k−1
i=1
P C A∑k λ N(wˆ ;He ,R )ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17)
+ ∑ Jk∑−1 D,k ka=1 k,a k,a k k,a k|k−1 k|k−1 k|k−1
w∈Wk i=1 κk(w)+(cid:82) PD,k(e)gk(w|e)vk|k−1(e)de (33)
= [1−P ]Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17)
D,k k|k−1 k|k−1 k|k−1
i=1 (cid:16) (cid:17)
+ ∑ Jk∑−1 A∑k PD,kCkλk,aN (wˆk(cid:82),a;Hek,Rk,a)ωk(i|)k−1N e;e(ki|)k−1,J(ki|)k−1 ,
w∈Wk i=1a=1 κk(w)+ PD,k(e)gk(w|e)·vk|k−1(e)de
and:
v (e,β =1)
k|k
g (w|e)Ψ(e;γ (θ,r,c,s))
= ∑ k k
(cid:82)
w∈Wk κk(w)+ PD,k(e)gk(w|e)vk|k−1(e)de (34)
Ψ(e;g (w|e)γ (θ,r,c,s))
= ∑ (cid:82) k k .
w∈Wk κk(w)+ PD,k(e)gk(w|e)vk|k−1(e)de
In the first step of Equation (34), the likelihood function g (w|e) is expressed by a Gaussian
k
distribution Equation (31), but not a Gaussian mixture Equation (29), as the target birth model
γ (θ,r,c,s) is given in polar coordinates (uniform across the bearing space and Gaussian in the
k
rangeandvelocity). Asbothofthemareexpressedinthesamepolarcoordinates,wehave:
g (w|e)Ψ(e;γ (θ,r,c,s)) = Ψ(e;g (w|e)γ (θ,r,c,s)). (35)
k k k k
Intheaboveequation,
ϕ (θ,r,c,s)(cid:44) g (w|e)γ (θ,r,c,s)
k k k
= N(cid:0)θt;θ,σ2(cid:1)ωbU(θ;H )N (cid:0)r;r¯,σ2(cid:1)N (cid:0)c;θ−π,σ2(cid:1)
k θ k θ r c
(cid:0) (cid:1)
·N s;s¯,σ2
s
= ωb1H(θ)N(cid:0)θ;θt,σ2(cid:1)N (cid:0)r;r¯,σ2(cid:1)N (cid:0)c;θ−π,σ2(cid:1) (36)
k VH (cid:0) k(cid:1) θ r c
·N s;s¯,σ2
s
≈ ωkbN(cid:0)θ;θt,σ2(cid:1)N (cid:0)r;r¯,σ2(cid:1)N (cid:0)c;θ−π,σ2(cid:1)N (cid:0)s;s¯,σ2(cid:1),
2π k θ r c s
where1H(θ)istheindicatorfunctionofthebearingspaceregionHθ andVH isthevolumeofHθ,and
weusetheapproximation:
(cid:16) (cid:17) (cid:16) (cid:17)
1H(θ)N θ;θkt,σθ2 ≈ N θ;θkt,σθ2 (37)
Sensors2016,16,1469 8of18
which is reasonable in practice by assuming that σ is small compared to the region H .
θ θ
The transformation function Ψ(e;ϕ (θ,r,c,s)) transforms the function ϕ (θ,r,c,s) from polar
k k
coordinatestoCartesiancoordinates,givenby:
Ψ(e;ϕ (θ,r,c,s)) ≈ ωkb ∑Ak λ N (cid:0)e;e˜ (cid:0)θt(cid:1),J˜ (cid:0)θt(cid:1)(cid:1), (38)
k 2π k,a k k k k
a=1
wheretheweightλ isgivenbyEquation(27),thepositionpartofthemeane˜ (cid:0)θt(cid:1)andthecovariance
k,a k k
J˜ (cid:0)θt(cid:1)aregivenbyEquations(24)and(25),respectively,andthevelocitypartiscalculatedbythe
k k
approximation[38,39]:
 xs +r¯sin(cid:0)θt(cid:1) 
k k
e˜k(cid:0)θkt(cid:1) =  ys¯sks+inr(cid:0)¯θcots−(cid:0)πθkt(cid:1)(cid:1) , (39)
 k 
s¯cos(cid:0)θt −π(cid:1)
k
 
P P 0 0
xx xy
J˜ (cid:0)θt(cid:1) =  Pyx Pyy 0 0 , (40)
k k  0 0 P P 
 x˙x˙ x˙y˙ 
0 0 P P
y˙x˙ y˙y˙
where:
P = (∆r)2sin2(cid:0)θt(cid:1)+r¯2σ2cos2(cid:0)θt(cid:1), (41)
xx k θ k
P = (∆r)2cos2(cid:0)θt(cid:1)+r¯2σ2sin2(cid:0)θt(cid:1), (42)
yy k θ k
1 (cid:104) (cid:105)
P = P = sin2θt (∆r)2−r¯2σ2 , (43)
xy yx 2 k θ
P = σ2sin2(cid:0)θt −π(cid:1)+σ2s¯2cos2(cid:0)θt −π(cid:1), (44)
x˙x˙ s k c k
P = σ2cos2(cid:0)θt −π(cid:1)+σ2s¯2sin2(cid:0)θt −π(cid:1), (45)
x˙x˙ s k c k
P = P = 1sin(cid:0)2(cid:0)θt −π(cid:1)(cid:1)(cid:104)σ2−σ2s¯2(cid:105). (46)
x˙y˙ y˙x˙ 2 k s c
InthedenominatorofEquations(33)and(34),wehavethefollowing:
(cid:82)
P (e)g (w|e)v (e)de
D,k k k|k−1
(cid:82)(cid:82)
= P (e,β)g (w|e,β)v (e,β)dedβ
D,k k k|k−1
(cid:82) 1
= ∑ P (e,β)g (w|e,β)v (e,β)de
D,k k k|k−1
β=0
= (cid:82) P C A∑k λ N (wˆ ;He ,R )Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17)de
D,k k k,a k,a k k,a k|k−1 k|k−1 k|k−1
a=1 i=1
+(cid:82) Ψ(e;g (w|e)γ (θ,r,c,s))de (47)
k k
= P C A∑k Jk∑−1λ ω(i) N (cid:16)wˆ ;Heˆ(i) ,S(i) (cid:17)(cid:82) N (cid:16)e;eˆ(i,a),J(i,a)(cid:17)de
D,k k k,a k|k−1 k,a k|k−1 k|k−1 k|k k|k
a=1i=1
+(cid:82) ωkb A∑k λ N (cid:0)e;e˜ (cid:0)θt(cid:1),J˜ (cid:0)θt(cid:1)(cid:1)de
2π k,a k k k k
a=1
= P C A∑k Jk∑−1λ ω(i) q(i)(wˆ )+ωkb,
D,k k k,a k|k−1 k k,a 2π
a=1i=1
(cid:16) (cid:17)
whereq(i)(wˆ ) = N wˆ ;Heˆ(i) ,S(i) isthelikelihoodofmeasurementcomponentwˆ against
k k,a k,a k|k−1 k|k−1 k,a
thecomponentiwiththeinnovationcovariance:
S(i) =HJ(i) H(cid:48) +R , (48)
k|k−1 k|k−1 k,a
Sensors2016,16,1469 9of18
andtheupdatedtargetstatesandcorrespondingcovariancearegivenby:
(cid:16) (cid:17)
eˆ(i,a) =eˆ(i) +K(i) wˆ −Heˆ(i) , (49)
k|k k|k−1 k k,a k|k−1
J(i,a) =J(i) −K(i)HJ(i) , (50)
k|k k|k−1 k k|k−1
withtheKalmangain:
K(i) =J(i) H(cid:48)(cid:16)S(i) (cid:17)−1. (51)
k k|k−1 k|k−1
Thefollowingapproximationfortheposteriorintensitycanbeobtained:
v (e,β =0)
k|k
= [1−P ]Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17)
D,k k|k−1 k|k−1 k|k−1
i=1
+ ∑ Jk∑−1 A∑k PD,kCkλk,aN (wˆk,a;Hek,Rk,a)ωk(i|)k−1N (cid:16)e;eˆ(ki|)k−1,J(ki|)k−1(cid:17) (52)
w∈Wk i=1a=1 κ (w)+P C A∑k Jk∑−1λ ω(i) q(i)(wˆ )+ωkb
k D,k k k,a k|k−1 k k,a 2π
a=1i=1
= Jk∑−1ω(i) N (cid:16)e;eˆ(i) ,J(i) (cid:17) + ∑ Jk∑−1 A∑k ω(i,a)N (cid:16)e;eˆ(i,a),J(i,a)(cid:17),
m,k k|k−1 k|k−1 s,k k|k k|k
i=1 w∈Wk i=1a=1
and:
v (e,β =1)
k|k
= ∑ A∑k ω2πkbλk,aN(e;e˜k(θkt),J˜k(θkt))
w∈Wka=1 κ (w)+P C A∑k Jk∑−1λ ω(i) q(i)(wˆ )+ωkb (53)
k D,k k k,a k|k−1 k k,a 2π
a=1i=1
= ∑ A∑k ω N (cid:0)e;e˜ (cid:0)θt(cid:1),J˜ (cid:0)θt(cid:1)(cid:1),
b,k k k k k
w∈Wka=1
where:
ω(i) = [1−P ]ω(i) , (54)
m,k D,k k|k−1
P C λ ω(i) q(i)(wˆ )
ω(i,a) = D,k k k,a k|k−1 k k,a , (55)
s,k κ (w)+P C A∑k Jk∑−1λ ω(i) q(i)(wˆ )+ωkb
k D,k k k,a k|k−1 k k,a 2π
a=1i=1
(cid:46)
λ ωb 2π
ω(a) = k,a k . (56)
b,k κ (w)+P C A∑k Jk∑−1λ ω(i) q(i)(wˆ )+ωkb
k D,k k k,a k|k−1 k k,a 2π
a=1i=1
3.4. ComponentManagementandStateExtraction
Just like the GM-PHD filter, the number of Gaussian components of the updated intensity
exponentially increases in time. To reduce the computational load, same techniques (component
merging and pruning) in [27,31] are also implemented in the GMM-PHD filter. If the weight of a
Gaussiancomponentislowerthanthepresetthreshold,itwillbediscarded.Whensomeofcomponents
arecloseenough,theywillbemergedintooneGaussiancomponent. IfthetotalnumberofGaussian
componentsisbiggerthanthemaximumvalue M ,only M componentswithhighprobabilitywillbe
k k
retained. Moredetailsaboutcomponentmanagementcanbefoundin[31]. Thestateextractionisalso
Sensors2016,16,1469 10of18
thesameasthemethodin[31]. Iftheweightsarebiggerthanathreshold,thecorrespondingstatesare
extractedastheoutputsofthefilter.
4. SimulationExperiments
ToillustratetheimprovementsoftheproposedGMM-PHDfilter,wecompareitwiththeGM-PHD-EKF
and the GM-PHD-UKF, which also use a partially-uniform birth model. The GM-PHD-EKF and
GM-PHD-UKFwereintroducedin[8]forachallengingbearings-onlymulti-targetscenario.
4.1. SimulationScenarios
4.1.1. Case1
The sensor is located on a maneuvering platform with the initial position [−4200m,3500m].
The sensor platform moves at a speed of five knots and changes course two times to ensure the
observabilityofthebearings-onlytargettracking. Notethatthesensorshouldchangecourseagain
totrackthetargetsthatappearafterthesecondsensormaneuvering. Theinitialcourseofthesensor
is 220◦and changes to 60◦ from 14 min to 17 min. The second course change occurs from 31 min
to 34 min, with the course changed from 60◦ to 220◦. Note that the clockwise rotation from the
positive Y-axis is considered to be positive. A time-varying number of targets is contained in the
scenario,andthemotionprofileoftargetsispresentedinTable1. ThestartingpositionsofTargets
#1and#2are[−8000m,−2500m]and[−3000m,−6500m],respectively.Similarly,theinitialpositions
ofTargets#3and#4are[4100m,−6100m]and[4200m,−2200m],respectively. Finally,thetravelof
Target#5startsat[6300m,4000m]. Anillustrationofthesensorandthetargetswithoutprocessnoise
andclutterispresentedinFigure2.
Table1.Motionprofileoftargets.
Target SurvivalTime(s) Course(degree) Speed(knots)
#1 [0, 2400] 95 8
#2 [300, 3000] 20 7
#3 [500, 3000] 280 8
#4 [0, 3000] 275 7
#5 [0, 2700] 215 10

4000
2000
0
Y (m)−2000
−4000
−6000
sensor
start points
−8000
end points

−8000 −6000 −4000 −2000 0 2000 4000 6000
X (m)
Figure2.Thegeometryofsensorandtargets.