A Poisson multi-Bernoulli mixture filter for coexisting point and extended targets
Ángel F. García-Fernández, Jason L. Williams, Lennart Svensson, Yuxuan Xia
11 A Poisson multi-Bernoulli mixture filter forcoexisting point and extended targets ´Angel F. Garc´ıa-Fern´andez, Jason L. Williams, Lennart Svensson, Yuxuan Xia
Abstract —This paper proposes a Poisson multi-Bernoulli mix-ture (PMBM) filter for coexisting point and extended targets.The PMBM filter provides a recursion to compute the filteringposterior based on single-target predictions and updates. In thispaper, we first derive the PMBM filter update for a generalisedmeasurement model, which can include measurements originatedfrom point and extended targets. Second, we propose a single-target space that accommodates both point and extended targetsand derive the filtering recursion that propagates Gaussiandensities for single targets and gamma Gaussian inverse Wishartdensities for extended targets. As a computationally efficient ap-proximation of the PMBM filter, we also develop a Poisson multi-Bernoulli (PMB) filter for coexisting point and extended targets.The resulting filters are analysed via numerical simulations.
Index Terms —Multiple target filtering, point targets, extendedtargets.
I. I
NTRODUCTION
Multiple target filtering refers to the sequential estimation ofthe states of the current targets, which may appear, move anddisappear, given past and current noisy sensor measurements.This is a key component in many applications such as self-driving vehicles [1] and maritime navigation [2]. Multi-targetfiltering can be solved in a Bayesian framework by computingthe posterior density on the current set of targets, givenprobabilistic models for target births, dynamics and deaths,and for sensor measurements [3].If the target extent is small compared to the sensor res-olution, it is common to use point-target modelling. In thismodel, a target state typically contains kinematic information,such as position and velocity, and one target can generate atmost one measurement at each time step [4]. Conversely, ifthe target extent is large compared to the sensor resolution,a better choice is to use extended target modelling [5]. Here,the target state usually contains both kinematic informationand information on its extent, e.g., represented by an ellipse[6]. In addition, each extended target may generate more thanone measurement at each time step, represented via a Poissonpoint process (PPP) in the standard model [5], [7], [8].For both point and extended targets with Poisson birthmodel and the standard measurement models, the posteriordensity is a Poisson multi-Bernoulli mixture (PMBM), which
A. F. Garc´ıa-Fern´andez is with the Department of Electrical Engineer-ing and Electronics, University of Liverpool, Liverpool L69 3GJ, UnitedKingdom ([email protected]). J. L. Williams is withthe Commonwealth Scientific and Industrial Research Organization ([email protected]). L. Svensson and Y. Xia are with the Depart-ment of Electrical Engineering, Chalmers University of Technology, SE-41296 Gothenburg, Sweden (fi[email protected]). can be calculated by the corresponding PMBM filtering re-cursions [9]–[11]. The PMBM has a compact representationof global hypotheses, representing undetected targets via theintensity of a PPP and making use of probabilistic targetexistence in each global hypothesis. The PMBM recursion canalso handle a multi-Bernoulli birth model by setting the PPPintensity to zero and adding new Bernoulli components in theprediction step, resulting in the MBM filter [10], [12]. TheMBM filter can also be extended to consider multi-Bernoulliswith deterministic target existence, which we refer to as theMBM filter, at the expense of increasing the number ofglobal hypotheses [10, Sec. IV]. Both MBM and MBM filters can consider target states with labels, and the (labelled)MBM filtering recursion is analogous to the δ -generalisedlabelled multi-Bernoulli ( δ -GLMB) filter [13], [14].There are applications in which it is important to have moregeneral models than the standard point and extended targetmodels [3]. Specifically, there may be some targets that aresmall compared to the sensor resolution, while other targets arelarge, which implies that there are coexisting point/extendedtargets in the field of view. For example, in a self-drivingvehicle application, pedestrians may be modelled as pointtargets while other vehicles as extended targets. The distinctionbetween point and extended targets may also depend on thedistance, as sensor resolution is usually higher at short dis-tances. Therefore, it is of interest to develop multi-target filtersthat can handle coexisting point and extended targets. Theextended target measurement models in [14], [15] are generalenough to model measurements from coexisting point andextended target, but no target models or filter implementationsare presented.In this paper, we fill this gap and propose a PMBMfilter for coexisting point and extended targets. In order todo so, we first develop a PMBM filtering recursion for ageneralised measurement model, for which the standard pointand extended target measurement models are particular cases.As a result, the PMBM filter with the generalised measurementmodel can be used to address more general multi-target filter-ing problems than the standard point/extended target PMBMfilters [9]–[11]. For example, the generalised measurementmodel can also be used for diffuse multipath [16], extendedtargets composed of point-scatterers [5], and point targetswith stationary landmarks, modelled as extended targets. Theresulting PMBM recursion has a track-oriented form thatenables efficient implementation [17]. a r X i v : . [ s t a t . M E ] N ov ased on the developed PMBM filtering recursion, thesecond contribution is to derive a PMBM filter for coexistingpoint and extended targets. In this setting, each Bernoullicontains probabilistic information on target existence andtype, either point or extended target. The implementation isprovided for a linear Gaussian model for point targets [18]and a Gamma Gaussian Inverse Wishart (GGIW) model forextended targets [5], [11], [19]. Finally, we explain how aPMBM density in this context can be projected onto a Poissonmulti-Bernoulli (PMB) density [9]. Performing this projectionafter each update provides us with a PMB filter, which is afast approximation to the PMBM filter. Simulation results areprovided to analyse the performance of the filters.The rest of the paper is organised as follows. Section IIintroduces the problem formulation and an overview of thesolution. The update for the PMBM filter with generalisedmeasurement model is derived in Section III. Section IVexplains the PMBM filter for coexisting point and extendedtargets, and the PMB projection. Simulation results and con-clusions are given in Sections V and VI, respectively.II. P ROBLEM FORMULATION AND OVERVIEW OF THESOLUTION
This paper deals with multiple target tracking with bothpoint and extended target models. This section presents anoverview of a PMBM filter with a generalised measurementmodel that will be used to model coexisting point and extendedtargets in Section IV. We introduce the models in Section II-Aand the PMBM filter overview in Section II-B.
A. Models
A single target state x ∈ X , where X is a locally compact,Hausdorff and second-countable (LCHS) space [3], containsthe information of interest about the target, for example, itsposition, velocity and extent. The set of targets at time k isdenoted X k ∈ F ( X ) , where F ( X ) represents the set of finitesubsets of X .A main novelty in this paper is the development of a PMBMfilter with a generalised measurement model. Here, the set X k of targets at time step k , is observed through a set Z k ∈F ( R n x ) of noisy measurements, which consist of the unionof target-generated measurements and clutter, with the model: • Each target x ∈ X k generates an independent set Z ofmeasurements with density f ( Z | x ) . • Clutter is a PPP with intensity λ C ( · ) .It should be noted that the standard point and extendedmeasurement models [9], [11] can be recovered by suitablechoices of f ( Z | x ) .We also consider the standard dynamic model for targets.Given the set X k of targets at time step k , each target x ∈ X k survives with probability p S ( x ) and moves to a new state witha transition density g ( · | x ) , or dies with probability − p S ( x ) .At time step k , targets are born independently following aPoisson point process (PPP) with intensity λ Bk ( · ) . B. PMBM posterior
In this paper, we show that for the above-mentioned mea-surement and dynamic models, the density f k | k (cid:48) ( · ) of X k given the sequence of measurements ( Z , ..., Z k (cid:48) ) , where k (cid:48) ∈ { k − , k } , is a PMBM density. This section providesan overview of the PMBM posterior and its data associationhypotheses.The PMBM is of the form [9], [10] f k | k (cid:48) ( X k ) = (cid:88) Y (cid:93) W = X k f p k | k (cid:48) ( Y ) f mbm k | k (cid:48) ( W ) (1) f p k | k (cid:48) ( X k ) = e − (cid:82) λ k | k (cid:48) ( x ) dx (cid:89) x ∈ X k λ k | k (cid:48) ( x ) (2) f mbm k | k (cid:48) ( X k ) = (cid:88) a ∈A k | k (cid:48) w ak | k (cid:48) (cid:88) (cid:93) nk | k (cid:48) l =1 X l = X k n k | k (cid:48) (cid:89) i =1 f i,a i k | k (cid:48) (cid:0) X i (cid:1) (3)where λ k | k (cid:48) ( · ) is the intensity of the PPP f p k | k (cid:48) ( · ) , representingundetected targets, and f mbm k | k (cid:48) ( · ) is a multi-Bernoulli mixturerepresenting potential targets that have been detected at somepoint up to time step k (cid:48) . Symbol (cid:93) denotes the disjoint unionand the summation in (1) is taken over all mutually disjoint(and possibly empty) sets Y and W whose union is X k .In the PMBM posterior, there are n k | k (cid:48) Bernoulli compo-nents and for each Bernoulli there are h ik | k (cid:48) possible localhypotheses. By selecting a local hypothesis a i ∈ (cid:110) , ..., h ik | k (cid:48) (cid:111) for each Bernoulli, we obtain a global hypothesis a = (cid:0) a , ..., a n k | k (cid:48) (cid:1) ∈ A k | k (cid:48) , where A k | k (cid:48) is the set of global hy-potheses. Each global hypothesis represents a multi-Bernoullidistribution. The i -th Bernoulli component with local hypoth-esis a i has a density f i,a i k | k (cid:48) ( X ) = − r i,a i k | k (cid:48) X = ∅ r i,a i k | k (cid:48) p i,a i k | k (cid:48) ( x ) X = { x } (4)where r i,a i k | k (cid:48) is the probability of existence and p i,a i k | k (cid:48) ( · ) thesingle target density. The weight of global hypothesis a is w ak | k (cid:48) and meets w ak | k (cid:48) ∝ n k | k (cid:48) (cid:89) i =1 w i,a i k | k (cid:48) (5)where w i,a i k | k (cid:48) is the weight of the i -th Bernoulli with localhypothesis a i , and (cid:80) a ∈A k | k (cid:48) w ak | k (cid:48) = 1 .The set of feasible global hypotheses is defined as in theextended target case [11], [20]. We denote the measurementset at time step k as Z k = (cid:8) z k , ..., z m k k (cid:9) . We refer tomeasurement z jk using the pair ( k, j ) and the set of all suchmeasurement pairs up to (and including) time step k is denotedby M k . Then, a single target hypothesis a i for the i -thBernoulli component has a set of measurement pairs denotedas M i,a i k ⊆ M k . The set A k | k (cid:48) of all global hypotheses meets A k | k (cid:48) = (cid:110)(cid:0) a , ..., a n k | k (cid:48) (cid:1) : a i ∈ (cid:110) , ..., h ik | k (cid:48) (cid:111) ∀ i, n k | k (cid:48) (cid:91) i =1 M i,a i k = M k , M i,a i k ∩ M j,a j k = ∅ , ∀ i (cid:54) = j (cid:41) . hat is, all measurements must be assigned to a local hy-pothesis, and there cannot be more than one local hypothesiswith the same measurement. More than one measurementcan be associated to the same local hypothesis at the sametime step. Each global hypothesis therefore corresponds to aunique partition of M k [11, Sec. V], and the number of globalhypothesis is the Bell number of |M k | . At each time step, eachnon-empty subset of Z k generates a new Bernoulli component,corresponding to a potential target detected for the first timeor clutter. This implies that, at each time step, m k − newBernoulli components are generated.It should be noted that the prediction step of a PMBMdensity is closed-form for the standard dynamic models [9],[10], and is not affected by the choice of measurement model.Therefore, the next section focuses on the update step and weomit the details for prediction, which can be found in [9], [10].III. PMBM UPDATE FOR A GENERALISED MEASUREMENTMODEL
This section provides the PMBM filter update step with themeasurement model in Section II. We denote a Kronecker deltaas δ i [ · ] , with δ i [ u ] = 1 if u = i and δ i [ u ] = 0 , otherwise.Also, given two real-valued functions a ( · ) and b ( · ) on thetarget space, we denote their inner product as (cid:104) a, b (cid:105) = (cid:90) a ( x ) b ( x ) dx. (6) A. Update
The update of the predicted PMBM f k | k − ( · ) after observ-ing Z k is given in the following theorem. Theorem 1.
Assume the predicted density f k | k − ( · ) is aPMBM of the form (1) . Then, the updated density f k | k ( · ) with set Z k = (cid:8) z k , ..., z m k k (cid:9) is a PMBM with the followingparameters. The number of Bernoulli components is n k | k = n k | k − + 2 m k . The intensity of the PPP is λ k | k ( x ) = f ( ∅| x ) λ k | k − ( x ) . (7) For Bernoullis continuing from previous time steps i ∈ (cid:8) , ..., n k | k − (cid:9) , a new local hypothesis is included for eachprevious local hypothesis and either a misdetection or anupdate with a non-empty subset of Z k . The updated number oflocal hypotheses is h ik | k = 2 m k h ik | k − . For missed detectionhypotheses, i ∈ (cid:8) , ..., n k | k − (cid:9) , a i ∈ (cid:110) , ..., h ik | k − (cid:111) , weobtain M i,a i k = M i,a i k − (8) l i,a i , ∅ k | k = (cid:10) f i,a i k | k − , f ( ∅|· ) (cid:11) (9) w i,a i k | k = w i,a i k | k − (cid:104) − r i,a i k | k − + r i,a i k | k − l i,a i , ∅ k | k (cid:105) (10) r i,a i k | k = r i,a i k | k − l i,a i , ∅ k | k − r i,a i k | k − + r i,a i k | k − l i,a i , ∅ k | k (11) f i,a i k | k ( x ) = f ( ∅| x ) f i,a i k | k − ( x ) l i,a i , ∅ k | k . (12) Let Z k , ..., Z mk − k be the nonempty subsets of Z k . For aBernoulli i ∈ (cid:8) , ..., n k | k − (cid:9) with a single target hypothesis (cid:101) a i ∈ (cid:110) , ..., h ik | k − (cid:111) in the predicted density, the new localhypothesis generated by a set Z jk has a i = (cid:101) a i + h ik | k − j , r i,a i k | k = 1 , and M i,a i k = M i, (cid:101) a i k − ∪ (cid:110) ( k, p ) : z pk ∈ Z jk (cid:111) (13) l i,a i ,Z jk k | k = (cid:28) f i, (cid:101) a i k | k − , f (cid:16) Z jk |· (cid:17) (cid:29) w i,a i k | k = w i, (cid:101) a i k | k − l i,a i ,Z jk k | k (14) f i,a i k | k ( x ) = f (cid:16) Z jk | x (cid:17) f i, (cid:101) a i k | k − ( x ) l i,a i ,Z jk k | k . (15) For the new Bernoulli initiated by subset Z jk , whose index is i = n k | k − + j , we have two single target hypotheses ( h ik | k =2 ), one corresponding to a non-existent Bernoulli M i, k = ∅ , w i, k | k = 1 , r i, k | k = 0 (16) and the other M i, k = (cid:110) ( k, p ) : z pk ∈ Z jk (cid:111) (17) l Z jk k | k = (cid:28) λ k | k − , f (cid:16) Z jk |· (cid:17) (cid:29) (18) w i, k | k = δ (cid:104) | Z jk | (cid:105) (cid:89) z ∈ Z jk λ C ( z ) + l Z jk k | k (19) r i, k | k = l Z jk k | k w i,a i k | k (20) f i, k | k ( x ) = f (cid:16) Z jk | x (cid:17) λ k | k − ( x ) l Z jk k | k . (cid:3) (21)Theorem 1 is proved in Appendix A. We can see that theupdated PPP intensity in (7) corresponds to the predictedintensity multiplied by the probability of not receiving anymeasurements. This is expected as the PPP contains informa-tion on the undetected targets. Misdetection hypotheses lowerthe probability of existence of the Bernoullis via (9) and (11).If f ( ∅| x ) does not depend on x , the single-target densities ofmisdetection hypotheses remain unchanged, see (12).For the update of a previous Bernoulli component withsubset Z jk , the updated Bernoulli has a probability of existenceequal to one. Each non-empty subset Z jk ⊆ Z k creates anew Bernoulli component. If | Z jk | > , the existence prob-ability r i, k | k of the new Bernoulli component is one, whichimplies that, conditioned on the corresponding hypothesis,this Bernoulli represents an existing target. If | Z jk | = 1 , theexistence probability r i, k | k of the new Bernoulli componentdepends on the clutter intensity λ C ( · ) , as this Bernoulli maycorrespond to a target or to clutter. The higher λ C ( · ) , thelower the probability of existence of this potential target. . Relation to standard point/extended target models In the standard point target measurement model, a target x isdetected with probability p D ( x ) and, if detected, it generatesone measurement with density l ( ·| x ) . This model is obtainedby setting f ( Z | x ) = − p D ( x ) Z = ∅ p D ( x ) l ( z | x ) Z = { z } | Z | > . (22)If we use the above definitions of local and global hypothesesand Theorem 1 for point targets, many of the global hypothesescontain local hypotheses where more than one measurement isassociated to the same Bernoulli at the same time step. Sincethis is impossible according to (22), all these hypotheses wouldobtain weight zero. A more convenient way to handle pointtargets is to exclude these hypotheses from the set that weconsider, see [9].In the standard extended target model, a target x is detectedwith probability p D ( x ) and, if detected, it generates a PPPmeasurement with intensity γ ( x ) l ( ·| x ) , where l ( ·| x ) is asingle-measurement density and γ ( x ) is the expected numberof measurements. We can recover this model by setting f ( Z | x ) = (cid:40) − p D ( x ) + p D ( x ) e − γ ( x ) Z = ∅ p D ( x ) γ | Z | ( x ) e − γ ( x ) (cid:81) z ∈ Z l ( z | x ) | Z | > . (23)In this case, Theorem 1 becomes the standard extended-targetPMBM update in track-oriented form [11], [20]. C. Discussion
We have shown that the update of a PMBM prior with thegeneralised measurement model in Section II is also PMBM.The proposed measurement model contains the standard pointtarget and extended target measurement models as particularcases, and can be used for other types of measurementmodelling. For example, another important special case is thateach target could generate a union of independent Bernoullimeasurements, which can model extended targets that consistof reflection points [21], [22]. It can also model extendedtargets with binomially distributed target-generated measure-ments [23]. The considered measurement model also allows usto model coexisting point and extended targets, for example,modelling radar returns from vehicles (extended targets) andpedestrians (point targets), as will be explained in Section IV.The proposed PMBM update requires PPP clutter, which canbe relaxed in Bernoulli filters [24].We would also like to remark that we have presented theresults for PPP birth density, as we think this is generallythe most suitable birth process, due to the lower number ofgenerated hypotheses [10], [12]. Nevertheless, the presentedresults also hold for the following cases. For multi-Bernoullibirth, the above equations are valid, by setting the Poissonintensity equal to zero, and adding the Bernoulli componentsfor new born targets in the prediction step [10], [12]. In thiscase, the posterior is a multi-Bernoulli mixture (MBM), which can also be represented as
MBM [10, Sec. IV]. For multi-Bernoulli birth, one can also uniquely label each Bernoullicomponent, for which the labelled MBM recursion wouldcorrespond to the δ -GLMB filter recursion [14].IV. PMBM FILTER FOR COEXISTING POINT ANDEXTENDED TARGETS
This section presents the PMBM filter, and a track-orientedPMB filter, for coexisting point and extended targets. Thesingle target space for point targets is R n x , which representsthe kinematic state (e.g. position and velocity). We modelextended targets with the GGIW model [15], whose space is X e = R + × R n x × S d + , where R + represents the positive realnumbers and S d + the positive definite matrices of size d , whichis the dimension of the extent.The single target space for coexisting point/extended targetsis then X = R n x (cid:93)X e . If x ∈ X e , then x = ( γ, ξ, X ) , where γ represents the expected number of measurements per target, ξ is the kinematic state and X is the extent state that describesthe target’s size and shape. It should be noted that, though notnecessary, it is also possible to include a class variable in thetarget space to distinguish between point and extended targets,as in interacting multiple models [25], see Appendix B.We use a measurement model that corresponds to the stan-dard point and extended target measurement models dependingon the type of target we observe. That is, for x ∈ R n x , f ( Z | x ) is given by (22) with a probability p D ( x ) = p D of detection, l ( z | x ) = N ( z ; H x, R ) where H is the measurement matrix, R is the noise covariance matrix, and N ( · ; x, P ) is a Gaussiandensity with mean x and covariance P . For x ∈ X e , f ( Z | x ) is given by (23) with a probability p D ( x ) = p D of detection, γ ( x ) = γ , and l ( z | x ) = N ( z ; H ξ, X ) where H is themeasurement matrix.The rest of this section is organised as follows. Section IV-Apresents the considered single-target densities. The update andthe prediction are provided in Sections IV-B and IV-C. ThePMB approximation is addressed in Section IV-D. Practicalaspects are discussed in Section IV-E. A. Single-target densities
We develop a PMBM implementation in which we propa-gate a Gaussian for single target densities and a (factorised)GGIW density for extended target densities [5], [19], [26].The Gaussian density for x ∈ X with mean x i,a i , k | k (cid:48) andcovariance matrix P i,a i , k | k (cid:48) is N p (cid:16) x ; x i,a i , k | k (cid:48) , P i,a i , k | k (cid:48) (cid:17) = N (cid:16) x ; x i,a i , k | k (cid:48) , P i,a i , k | k (cid:48) (cid:17) (24)for x ∈ R n x and zero for x ∈ X e . Note that N p ( · ) is zeroevaluated at x ∈ X e , as it represents point targets.The Gamma density with parameters α > and β > isdenoted as G ( · ; α, β ) . The inverse Wishart density on matricesin S d + with v > d degrees of freedom and parameter matrix V ∈ S d + is denoted as IW (; v, V ) [27]. Then, the GGIWdensity for x ∈ X with parameters ζ i,a i k | k (cid:48) = (cid:16) α i,a i k | k (cid:48) , β i,a i k | k (cid:48) , x i,a i , k | k (cid:48) , P i,a i , k | k (cid:48) , v i,a i k | k (cid:48) , V i,a i k | k (cid:48) (cid:17) (25)s G e (cid:16) x ; ζ i,a i k | k (cid:48) (cid:17) = G (cid:16) γ ; α i,a i k | k (cid:48) , β i,a i k | k (cid:48) (cid:17) N (cid:16) ξ ; x i,a i , k | k (cid:48) , P i,a i , k | k (cid:48) (cid:17) × IW (cid:16) X ; v i,a i k | k (cid:48) , V i,a i k | k (cid:48) (cid:17) (26)for x ∈ X e and zero for x ∈ R n x .The single-target density of the i -th Bernoulli and localhypothesis a i is p i,a i k | k (cid:48) ( x ) = c i,a i k | k (cid:48) N p (cid:16) x ; x i,a i , k | k (cid:48) , P i,a i , k | k (cid:48) (cid:17) + (cid:16) − c i,a i k | k (cid:48) (cid:17) G e (cid:16) x ; ζ i,a i k | k (cid:48) (cid:17) (27)where c i,a i k | k (cid:48) and (cid:16) − c i,a i k | k (cid:48) (cid:17) are the probabilities that the targetis a point-target and and extended target, respectively. The PPPintensity is a mixture λ k | k (cid:48) ( x ) = n pk | k (cid:48) (cid:88) q =1 w p,qk | k (cid:48) N p (cid:16) x ; x p,q, k | k (cid:48) , P p,q, k | k (cid:48) (cid:17) + n ek | k (cid:48) (cid:88) q =1 w e,qk | k (cid:48) G e (cid:16) x ; ζ e,qk | k (cid:48) (cid:17) (28)where n pk | k (cid:48) is the number of components with point-targets,with weight w p,qk | k (cid:48) , mean x p,q, k | k (cid:48) and covariance P p,q, k | k (cid:48) , and n ek | k (cid:48) is the number of components with extended targets, withweight w e,qk | k (cid:48) and parameters ζ e,qk | k (cid:48) . It should be noted that (cid:80) n pk | k (cid:48) q =1 w p,qk | k (cid:48) and (cid:80) n ek | k (cid:48) q =1 w e,qk | k (cid:48) represent the expected numberof undetected point and extended targets, respectively. B. Update
We represent the update of a GGIW density with parameters ζ i,a i k | k − with a given measurement set Z jk as a function (cid:16) ζ e,qk | k , (cid:96) e,qk | k (cid:17) = u e (cid:16) ζ i,a i k | k − , Z jk (cid:17) where ζ e,qk | k is the updated GGIW and (cid:96) e,qk | k the marginallikelihood, see Appendix C. The Kalman filter update of aGaussian density with mean x i,a i , k | k − and covariance P i,a i , k | k − and measurement z is represented as (cid:16) x i,a i , k | k , P i,a i , k | k , (cid:96) i,a i , k | k (cid:17) = u p (cid:16) x i,a i , k | k − , P i,a i , k | k − , z (cid:17) where x i,a i , k | k and P i,a i , k | k are the updated mean and covariance,and (cid:96) i,a i , k | k is the marginal likelihood, see [18] for details.We apply Theorem 1 to obtain the specific parameters ofthe updated PMBM provided in the following lemma. Lemma 2.
The updated PMBM with a prior PMBM de-scribed by (1) , (27) and (28) , with measurement set Z k = (cid:8) z k , ..., z m k k (cid:9) has the structure in Theorem 1 with the follow-ing parameters. The number of PPP components is n pk | k = n pk | k − and n ek | k = 2 n ek | k − . For point targets, x p,q, k | k = x p,q, k | k − , P p,q, k | k = P p,q, k | k − , (29) w p,qk | k = (cid:0) − p D (cid:1) w p,qk | k − . (30) For extended targets and q ≤ n ek | k − , we have ζ e,qk | k = ζ e,qk | k − , w e,qk | k = (cid:0) − p D (cid:1) w e,qk | k − . (31) For q > n ek | k − , ˜ q = q − n ek | k − , (cid:16) ζ e,qk | k , (cid:96) e,qk | k (cid:17) = u e (cid:16) ζ e, ˜ qk | k − , ∅ (cid:17) (32) w e,qk | k = p D (cid:96) e,qk | k w e, ˜ qk | k . (33) For missed detection hypotheses of previous Bernoullis, p i,a i k | k ( x ) = c i,a i k | k N p (cid:16) x ; x i,a i , k | k , P i,a i , k | k (cid:17) + (cid:16) − c i,a i k | k (cid:17) × (cid:104) w G e (cid:16) x ; ζ i,a i , k | k (cid:17) + (1 − w ) G e (cid:16) x ; ζ i,a i , k | k (cid:17)(cid:105) (34) where x i,a i , k | k = x i,a i , k | k − , P i,a i , k | k = P i,a i , k | k − , ζ i,a i , k | k = ζ i,a i k | k − and (cid:16) ζ i,a i , k | k , (cid:96) i,a i , k | k (cid:17) = u e (cid:16) ζ i,a i k | k − , ∅ (cid:17) (35) l i,a i , ∅ k | k = c i,a i k | k − (cid:0) − p D (cid:1) + (cid:16) − c i,a i k | k − (cid:17) (cid:16) − p D + p D (cid:96) i,a i , k | k (cid:17) (36) w i,a i k | k = w i,a i k | k − (cid:104) − r i,a i k | k − + r i,a i k | k − l i,a i , ∅ k | k (cid:105) (37) r i,a i k | k = r i,a i k | k − l i,a i , ∅ k | k − r i,a i k | k − + r i,a i k | k − l i,a i , ∅ k | k (38) c i,a i k | k = (cid:0) − p D (cid:1) c i,a i k | k − l i,a i , ∅ k | k (39) w = 1 − p D − p D + p D (cid:96) i,a i , k | k . (40) The detection hypotheses of a previous Bernoulli with a subset Z jk , with (cid:12)(cid:12)(cid:12) Z jk (cid:12)(cid:12)(cid:12) = m jk , has r i,a i k | k = 1 , and w i,a i k | k = w i, (cid:101) a i k | k − r i, (cid:101) a i k | k − l i,a i ,Z jk k | k (41) (cid:16) ζ i,a i k | k , (cid:96) i,a i k | k (cid:17) = u e (cid:16) ζ i, (cid:101) a i k | k − , Z jk (cid:17) . (42) For m jk > , l i,a i ,Z jk k | k = p D (cid:96) i,a i k | k and c i,a i k | k = 0 , which impliesthat x i,a i , k | k and P i,a i , k | k are irrelevant. For m jk = 1 , Z jk = { z } ,we have (cid:16) x i,a i , k | k , P i,a i , k | k , (cid:96) i,a i , k | k (cid:17) = u p (cid:16) x i,a i , k | k − , P i,a i , k | k − , z (cid:17) (43) l i,a i ,Z jk k | k = c i,a i k | k − p D (cid:96) i,a i , k | k + (cid:16) − c i,a i k | k − (cid:17) p D (cid:96) i,a i k | k (44) c i,a i k | k = c i,a i k | k − p D (cid:96) i,a i , k | k l i,a i ,Z jk k | k . (45) For the new Bernoulli initiated by subset Z jk , the single targetdensity corresponding to existing Bernoulli is p i, k | k ( x ) = c i, k | k n pk | k − (cid:88) q =1 w q N p (cid:16) x ; x i, ,qk | k , P i, ,qk | k (cid:17) (cid:16) − c i, k | k (cid:17) n ek | k − (cid:88) q =1 w q G e (cid:16) x ; ζ i, ,qk | k (cid:17) (46) (cid:16) ζ i, ,qk | k , (cid:96) i, ,q ,k | k (cid:17) = u e (cid:16) ζ e,qk | k − , Z jk (cid:17) (47) w i, k | k = δ (cid:104) | Z jk | (cid:105) (cid:89) z ∈ Z jk λ C ( z ) + l Z jk k | k . (48) For m jk > , l Z jk k | k = p D n ek | k − (cid:88) q =1 w q (cid:96) i, ,q ,k | k , (49) r i, k | k = 1 and c i, k | k = 0 . For m jk = 1 , Z jk = { z } , we have (cid:16) x i, ,qk | k , P i, ,qk | k , (cid:96) i, ,q ,k | k (cid:17) = u p (cid:16) x p,q, k | k − , P p,q, k | k − , z (cid:17) (50) l Z jk k | k = p D n pk | k − (cid:88) q =1 w q (cid:96) i, ,q ,k | k + p D n ek | k − (cid:88) q =1 w q (cid:96) i, ,q ,k | k (51) r i, k | k = l Z jk k | k w i,a i k | k (52) c i, k | k = p D (cid:80) n pk | k − q =1 w q (cid:96) i, ,q ,k | k l Z jk k | k . (cid:3) (53)Lemma 2 is obtained by using Theorem 1 and the GGIWand Gaussian updates [11], [18]. We can see that the numberof components in the PPP corresponding to extended targetsdoubles in the update. This is due to the fact that the likelihoodfor misdetection for extended targets, see (23), has two terms − p D and p D e − γ , which create two updated PPP componentsfor each prior PPP component. For the same reason, inthe update of previous Bernoullis with a misdetection, theextended target updated density is a mixture of two GGIW,see (34). As only the Gamma distribution differs in the twoupdated GGIWs, we apply merging for Gamma densities [28]to obtain an updated single-target density of the form (27).For the detection of previous Bernoullis, the hypothesisrepresents with probability r i,a i k | k = 1 that there is target. If m jk > , the target is an extended target with probabilityone ( c i,a i k | k = 0 ). If m jk = 1 , the target may be a point or anextended target. For the new Bernoulli components, if m jk > ,the local hypotheses represent an existing extended target withprobability one. For m jk = 1 , the new Bernoulli may representclutter, a single target or an extended target. All possible clutterevents are accounted for in the hypotheses with m jk = 1 andso do not need to be duplicated in events with m jk > . We canalso see that the single target density (46) for new Bernoullicomponents is a mixture for both point and extended targets.To obtain an updated density as in (27), we perform mergingof the Gaussian and GGIW mixtures [28], [29]. C. Prediction
We consider that the probability of survival is a constant p S ( · ) = p S and linear/Gaussian dynamics for point targets.That is, for x ∈ R n x , we have g ( · | x ) = N ( · ; F x, Q ) (54)where F is the transition matrix and Q is the process noise co-variance matrix. For GGIW targets, there are several dynamicmodels [5], [6]. In the simulations, we use the one in [11].The target birth intensity is of the form λ Bk ( x ) = n b,pk (cid:88) q =1 w b,p,qk N p (cid:16) x ; x b,p,q, k , P b,p,q, k (cid:17) + n b,ek (cid:88) q =1 w b,e,qk G e (cid:16) x ; ζ b,e,qk (cid:17) . (55)We apply the PMBM prediction step [9], [10] to obtaina PMBM with the following parameters. Given a single-target filtering density p i,a i k − | k − ( · ) of the form (27), then thepredicted density is of the same form with c i,a i k | k − = c i,a i k − | k − and ζ i,a i k | k − = p e (cid:16) ζ i,a i k − | k − (cid:17) (56) (cid:16) x i,a i , k | k − , P i,a i , k | k − (cid:17) = p p (cid:16) x i,a i , k − | k − , P i,a i , k − | k − (cid:17) (57)where p p ( · ) and p e ( · ) denote the Kalman filter [18] and theextended target GGIW prediction [11, Tab. III], respectively.The predicted PPP is λ k | k − ( x ) = n pk − | k − (cid:88) q =1 p S w p,qk − | k − N p (cid:16) x ; x p,q, k | k − , P p,q, k | k − (cid:17) + n ek − | k − (cid:88) q =1 p S w e,qk − | k − G e (cid:16) x ; ζ e,qk | k − (cid:17) + λ Bk ( x ) where (cid:16) x p,q, k | k − , P p,q, k | k − (cid:17) = p p (cid:16) x p,q, k − | k − , P p,q, k − | k − (cid:17) and ζ e,qk | k − = p e (cid:16) ζ e,qk − | k − (cid:17) . D. PMB approximation
It is also useful to consider a PMB approximation to thePMBM (1) to develop a faster algorithm. If we perform thisapproximation after each update, we obtain the correspondingPMB filter [9], [30]. Given (1), the track-oriented PMBapproximation is f pmb k | k (cid:48) ( X k ) = (cid:88) Y (cid:93) W = X k f p k | k (cid:48) ( Y ) f mb k | k (cid:48) ( W ) (58) f mb k | k (cid:48) ( X k ) = (cid:88) (cid:93) nk | k (cid:48) l =1 X l = X k n k | k (cid:48) (cid:89) i =1 f ik | k (cid:48) (cid:0) X i (cid:1) (59)where f ik | k (cid:48) ( · ) is a Bernoulli density with probability r i ofexistence and single target density p i ( · ) such that r i = h i (cid:88) a i =1 w i,a i k | k (cid:48) r i,a i k | k (cid:48) (60) i ( x ) = (cid:80) h i a i =1 w i,a i k | k (cid:48) r i,a i k | k (cid:48) p i,a i k | k (cid:48) ( x ) r i (61) w i,a i k | k (cid:48) = (cid:88) b ∈A k | k (cid:48) : b i = a i w bk | k (cid:48) . (62)The PMB approximation (58)-(59) minimises the Kullback-Leibler divergence (KLD) on a single target space augmentedwith an auxiliary variable, which represents if the targetremains undetected or corresponds to the i -th Bernoulli com-ponent [31]. We can see that (61) is a mixture over all localhypotheses and that the PPP part of the PMBM (1) is notaffected by the PMB approximation.In the implementation for coexisting point-extended targets,we are interested in single target densities of the form (27). Byusing moment matching (KLD minimisation) for the mixturein p i ( · ) , we obtain the single-target density p i ( x ) = c i N p (cid:16) x ; x ik | k (cid:48) , P ik | k (cid:48) (cid:17) + (cid:0) − c i (cid:1) G e (cid:16) x ; ζ ik | k (cid:48) (cid:17) (63) c i = (cid:80) h i a i =1 w i,a i k | k (cid:48) r i,a i k | k (cid:48) c i,a i k | k (cid:48) r i (64) (cid:16) x ik | k (cid:48) , P ik | k (cid:48) (cid:17) = m G (cid:18) x i,a , k | k (cid:48) , P i,a , k | k (cid:48) , ..., , x i,a hi , k | k (cid:48) , P i,a hi , k | k (cid:48) ,β i,a G , ..., β i,a hi G (cid:19) (65) β i,a i G ∝ w i,a i k | k (cid:48) r i,a i k | k (cid:48) c i,a i k | k (cid:48) (66) ζ ik | k (cid:48) = m GG (cid:18) ζ i,a k | k (cid:48) , ..., ζ i,a hi k | k (cid:48) , β i,a GG , ..., β i,a hi GG (cid:19) (67) β i,a i GG ∝ w i,a i k | k (cid:48) r i,a i k | k (cid:48) (cid:16) − c i,a i k | k (cid:48) (cid:17) (68)where m G ( · ) is a function that obtains the meanand covariance of a Gaussian mixture with weights β i,a G , ..., β i,a hi G (normalised to sum to one) and mo-ments x i,a , k | k (cid:48) , P i,a , k | k (cid:48) , ..., , x i,a hi , k | k (cid:48) , P i,a hi , k | k (cid:48) [32]. The function m GG ( · ) obtains the GGIW parameters that minimise the KLDfrom a mixture with weights β i,a GG , ..., β i,a hi GG (normalised tosum to one) and parameters ζ i,a k | k (cid:48) , ..., ζ i,a hi k | k (cid:48) [28], [29]. E. Practical aspects
As in other multiple target filters with data associations, thenumber of global and local hypotheses increases unboundedlyin time. Therefore, in practice, it is necessary to performapproximations, with the objective of only propagating hy-potheses with relevant weights. In fact, due to the structureof the hypotheses of Theorem 1, the way to handle the dataassociation problem with coexisting point and extended targetsis quite similar to the extended target case [5], [11].In our implementation, the PMBM posterior is representedby a list of Bernoullis i ∈ (cid:8) , ..., n k | k (cid:48) (cid:9) , where each ofthem contains their local hypotheses with their parameters,a global hypothesis table, which contains indices to local hypotheses of each Bernoulli, and a vector with the globalhypotheses weights. To deal with the data association problemat each update, we first perform gating to obtain two sets ofmeasurements: 1) Measurements that are in the gate of at leastone previous Bernoulli and 2) Measurements that are only inthe gate of PPP birth components. Measurements that do notfall into these categories are discarded.Then, for measurements in group 1), we first generatepossible partitions of this set using the DBSCAN algorithmwith different distance thresholds [33], [34]. As there may berepeated subsets of measurements in the possible partitionsobtained by DBSCAN, we obtain the unique subsets ofmeasurements in the previous partitions. These subsets arethen used to generate the updated local hypotheses for previousBernoullis, see (8)-(15). For each previous global hypothesis,we run Murty’s algorithm [35] to obtain the global hypotheseswith highest weights.For measurements in group 2), which may corresponds tonew targets, we run the DBSCAN to obtain possible partitions.Each of these partitions in theory gives rise to differentglobal hypotheses corresponding to new born targets. Wesimplify this procedure by finding the partition with highestweight and only generating the Bernoulli components thatare generated by the sets in this partition [11]. These newBernoulli components are added to all the global hypotheses,whose weights remain unchanged.We also perform pruning of global hypotheses with lowweights, and pruning of Bernoulli components with low ex-istence probabilities [10], [12], [36]. A pseudocode of theresulting PMBM update is provided in Algorithm 1. The PMBfilter performs the same PMBM update and it is then followedby the PMB approximation, see Section IV-D. It is alsopossible to approximate the PMB marginal data associationprobabilities directly using belief propagation [37]. Algorithm 1
Pseudocode of the PMBM update - Perform gating to separate current measurements into the follow-ing disjoint categories: ◦
1. A set of measurements that are in the gate of at least oneprevious Bernoulli. ◦
2. A set of measurements that are only in the gate of PPPbirth components.- For measurements corresponding to 1: ◦ Run DBSCAN to generate possible partitions. ◦ Obtain unique subsets in the previous partitions. ◦ Generate new local hypotheses for previous Bernoullis. ◦ For each previous global hypothesis, run Murty’s algorithm toobtain updated global hypotheses.- Perform pruning of global hypotheses and Bernoulli components.- For measurements corresponding to 2 (new targets): ◦ Run DBSCAN to generate possible partitions. ◦ Find the partition with highest weight. ◦ Generate the new Bernoulli components for this partition. ◦ Add these Bernoullis to the global hypotheses.
V. S
IMULATIONS
In this section, we assess the PMBM and PMB filters forcoexisting point and extended targets via numerical simula-tions. In this section, we refer to these filters as point-extendedPMBM and PMB (PE-PMBM and PE-PMB) filters. The filtersare implemented with the following parameters: maximumumber of hypotheses , threshold for pruning the PPPweights − , threshold for pruning Bernoulli components − and threshold for pruning global hypotheses − .The DBSCAN algorithm [33] is run with distance thresholdsbetween 0.1 and 12, with a step size of 0.1, and minimumnumber of points to form a region equal to 1. Target stateestimation is performed by obtaining this global hypothesiswith highest weight and reporting the means of the Bernoulliswhose existence is above a threshold 0.5 (Estimator 1 in [10]).We have also implemented a point-extended MBM (PE-MBM)filter, see Section III-C.Extended target filters can in principle deal with point-targetdetections, as they do not place zero probability to this event.Therefore, we compare the proposed filters with extendedtarget PMBM and PMB filters, which we refer to as E-PMBMand E-PMB filters [11], [34]. We proceed to discuss the modelsand the simulations results. All the units in this section aregiven in the international system. A. Models
We consider a point target state [ p x , ˙ p x , p y , ˙ p y ] T , whichcontains position and velocity in a two-dimensional plane.Point targets move with a nearly-constant velocity model with F = I ⊗ (cid:18) τ (cid:19) , Q = qI ⊗ (cid:18) τ / τ / τ / τ (cid:19) where τ = 1 , q = 0 . , ⊗ denotes Kronecker product and I isan identity matrix of size 2. The probability of survival is p S =0 . . The extended target model is the GGIW model in SectionIV. Extended targets move with the previous nearly-constantvelocity model and their extent matrix and γ parameter remainconstant. The probability of survival is 0.99.The birth model is a PPP of the form (55). The PPPpoint target part has parameters n b,pk = 1 , w b,p, k = 0 . , x b,p, , k = [0 , , , T , P b,p, , k = diag([200 , , , ]) .The extended target part has: n b,ek = 1 , w b,e, k | k (cid:48) = 0 . , and ζ b,e,qk = (cid:16) , , x b,p, , k , P b,p, , k , , I (cid:17) . As the birth covariance matrix is large, new born targetsmay appear in a large area. The multi-Bernoulli birth modelfor the PE-MBM filter has a single Bernoulli with existenceprobability 0.06, point-target probability c = 1 / , point targetmean x b,p, , k and covariance P b,p, , k , and GGIW ζ b,e,qk .We measure the positions of the targets. For point targets,we have parameters: p D = 0 . , and H = (cid:18) (cid:19) , R = σ I where σ = 1 . For extended targets, the parameters are p D = p D , H = H . Clutter is uniformly distributed in thesurveillance area [ − , × [ − , with an averageof λ C = 8 false alarms per scan. We consider 100 timesteps and the set of trajectories shown in Figure 1, which hasbeen obtained by sampling from the dynamic process. The E-PMBM and E-PMB filters are recovered by setting the birthintensity for point targets to zero in the PE-PMBM and PE-PMB filters. -400 -300 -200 -100 0 100 200 300 400 x axis (m) -400-300-200-1000100200300400 y a x i s ( m )
115 10115 10115 10
Fig. 1: Scenario of the simulations. Two extended targets are born at timestep 1 and two point targets are born at time steps 5 and 10. The extendedtargets are alive at all time steps. The last time steps of the point targetsare 38 and 60. Targets are in close proximity at around time step 50.Target states at time of birth are marked with a red cross, and every 10time steps with a icross. The 3- σ ellipse for extended targets is shownevery 10 time steps. -30 -20 -10 0 10 20 x axis (m) -20-15-10-505101520 y a x i s ( m ) Fig. 2: Ground truth and estimated set of targets at time step 52 in anillustrative run. Measurements are shown as black crosses. Blue ellipsesrepresent the true extended targets, a green cross represents the true pointtarget. The red, dashed ellipses represent the estimated extended targetsand the pink cross represents the estimated point target. The three targetsare properly detected and classified.
B. Results
We first show the ground truth and the estimate of the setof targets at time step 52 in an illustrative run with the PE-PMBM filter in Figure 2. We can see that the each extendedtarget generates several measurements and are detected. Theellipses of the estimated targets are reasonably accurate. Thepoint target generates a single measurement at this time step,and it is also detected. Its estimate is close to its true state.We evaluate filter performance via Monte Carlo simulationwith 100 runs. We compute the error between the true set oftargets at each time and its estimate using the generalised op-timal subpattern assignment (GOSPA) metric with parameters
20 40 60 80 100
Time step R M S GO SPA e rr o r PE-PMBMPE-PMBPE-MBME-PMBME-PMB
Time step Lo c a li s a t i on e rr o r Time step F a l s e t a r ge t e rr o r Time step M i ss ed t a r ge t e rr o r Fig. 3: RMS-GOSPA error (m) for the position elements and its de-composition. The PE-PMBM and PE-PMB filters have very similarperformance. PE-MBM fails to detect one of the targets at time step1. The E-PMBM and E-PMB have a higher error at some time steps dueto missed point targets. Localisation errors are also higher at some timesteps, as point targets are estimated with some extent. α = 2 , p = 2 , c = 10 , and its decomposition into localisationerrors and costs for missed and false targets [38]. The basemetric for target states is the Gaussian Wasserstein distance,which measures error for position and extent [39]. In the basemetric, we consider a point target as an extended target withextent zero.The root mean square GOSPA (RMS-GOSPA) error againsttime and its decomposition are shown in Figure 3. We can seethat PE-PMBM and PE-PMB filter perform quite similarlyand outperform E-PMBM and E-PMB. PE-MBM performsquite similarly to PE-PMBM and PE-PMB but does not detectone of the targets at time step 1, as the birth model setsthe maximum number of new born targets to one. For PE-PMBM and PE-PMB, missed target errors are higher whennew targets are born. False target errors are higher when targetsdie and when targets get in close proximity. Localisationerrors are higher at the beginning of the simulation, and whentargets get in close proximity, as the data association problemis more complicated. ET-PMBM and ET-PMB also behavequite similarly and have more difficulty in detecting the pointtargets, so they show a higher missed target error at some timesteps. In addition, the localisation error is also higher at sometime steps, as point targets are estimated with a certain extent,which increases the error compared to the ground truth.The running times of the Matlab implementations (100 timesteps) on an Intel Core i5 laptop are: 56.4s (PE-PMBM),17.5s (PE-PMB), 64.5 (PE-MBM), 25.5s (E-PMBM) and 15.2(E-PMBM). The PMB filters are considerably faster thanPMBM/MBM, as they do not propagate a mixture throughthe filtering recursion. Only considering extended targets isalso faster, though it decreases performance.To provide more complete simulation results, we show theRMS-GOSPA errors, along with the GOSPA error decompo-sition, considering all time steps for different values of theprobability of detection and clutter rate in Table I. Due to space constraints, we do not show E-PMB, which behavesquite similarly to E-PMBM. In this table, “Tot.”, “Loc.”,“Fal.” and “Mis.” refer to total GOSPA, localisation, falsetarget and missed target costs, respectively. The filters withcoexisting point extended targets consistently provide moreaccurate results, especially due to a lower number of missedtargets. The PE-PMBM and PE-PMB filters provide quitesimilar results though the PE-PMB filter is slightly better.While the PE-PMBM filter provides the closed-form solutionto the filtering recursion, we apply approximations and ansuboptimal estimator, so the PE-PMB filter may work betterin some scenarios. Decreasing the probability of detection orincreasing the clutter rate, the GOSPA error for all filtersincreases, mainly due to a rise in missed target cost.VI. C ONCLUSIONS
We have derived the update of a PMBM filter with ameasurement model that can consider point and extendedtargets, and we have shown that the updated posterior is alsoa PMBM. We have also proposed an implementation of theresulting PMBM recursion to consider coexisting point andextended targets. In order to do so, we first set the suitablesingle-target space and single target densities, which are basedon Gaussian densities for single targets, and GGIW densitiesfor extended targets. Finally, based on the previous results, wehave explained how to obtain a computationally-lighter PMBfilter for coexisting point and extended targets.We think there are many lines of future work. In manyapplications, there are coexisting point and extended targets,and one can perform research into tailored measurement andtarget models for each application. Another line of future workis to extend the above results to consider PMBMs on setsof trajectories, with coexisting point and extended targets, toprovide full trajectory information [20], [40], [41].R
EFERENCES[1] J. Choi, S. Ulbrich, B. Lichte, and M. Maurer, “Multi-target trackingusing a 3D-lidar sensor for autonomous vehicles,” in , 2013, pp. 881–886.[2] E. F. Brekke, E. F. Wilthil, B.-O. H. Eriksen, D. K. M. Kufoalor, Ø. K.Helgesen, I. B. Hagen, M. Breivik, and T. A. Johansen, “The Autoseaproject: Developing closed-loop target tracking and collision avoidancesystems,”
Journal of Physics: Conference Series , vol. 1357, pp. 1–12,Oct. 2019.[3] R. P. S. Mahler,
Advances in Statistical Multisource-Multitarget Infor-mation Fusion . Artech House, 2014.[4] S. S. Blackman, “Multiple hypothesis tracking for multiple targettracking,”
IEEE Aerospace and Electronic Systems Magazine , vol. 19,no. 1, pp. 5–18, Jan. 2004.[5] K. Granstr¨om, M. Baum, and S. Reuter, “Extended object track-ing:introduction, overview, and applications,”
Journal of Advances inInformation Fusion , vol. 12, no. 2, pp. 139–174, Dec. 2017.[6] J. W. Koch, “Bayesian approach to extended object and cluster trackingusing random matrices,”
IEEE Transactions on Aerospace and ElectronicSystems , vol. 44, no. 3, pp. 1042–1059, Jul. 2008.[7] K. Gilholm, S. Godsill, S. Maskell, and D. Salmond, “Poisson modelsfor extended and group tracking,” in
Proc. SPIE 5913, Signal and DataProcessing of Small Targets , vol. 5913, 2005, pp. 1–12.[8] X. Tang, M. Li, R. Tharmarasa, and T. Kirubarajan, “Seamless trackingof apparent point and extended targets using Gaussian process PMHT,”
IEEE Transactions on Signal Processing , vol. 67, no. 18, pp. 4825–4838,Sep. 2019.ABLE I: RMS-GOSPA errors and their decompositions for the filters and different parametersPE-PMBM PE-PMB PE-MBM E-PMBM p D = p D λ c Tot. Loc. Fal. Mis. Tot. Loc. Fal. Mis. Tot. Loc. Fal. Mis. Tot. Loc. Fal. Mis. .
95 8 .
95 16 .
85 8 .
85 16
IEEE Transactions onAerospace and Electronic Systems , vol. 51, no. 3, pp. 1664–1687, July2015.[10] A. F. Garc´ıa-Fern´andez, J. L. Williams, K. Granstr¨om, and L. Svensson,“Poisson multi-Bernoulli mixture filter: direct derivation and imple-mentation,”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 54, no. 4, pp. 1883–1901, Aug. 2018.[11] K. Granstr¨om, M. Fatemi, and L. Svensson, “Poisson multi-Bernoullimixture conjugate prior for multiple extended target filtering,”
IEEETransactions on Aerospace and Electronic Systems , vol. 56, no. 1, pp.208–225, Feb. 2020.[12] A. F. Garc´ıa-Fern´andez, Y. Xia, K. Granstr¨om, L. Svensson, and J. L.Williams, “Gaussian implementation of the multi-Bernoulli mixture fil-ter,” in
Proceedings of the 22nd International Conference on InformationFusion , 2019.[13] B. T. Vo and B. N. Vo, “Labeled random finite sets and multi-objectconjugate priors,”
IEEE Transactions on Signal Processing , vol. 61,no. 13, pp. 3460–3475, July 2013.[14] M. Beard, S. Reuter, K. Granstr¨om, B. Vo, B. Vo, and A. Scheel,“Multiple extended target tracking with labeled random finite sets,”
IEEETransactions on Signal Processing , vol. 64, no. 7, pp. 1638–1653, 2016.[15] C. Lundquist, K. Granstr¨om, and U. Orguner, “An extended target CPHDfilter and a gamma Gaussian inverse Wishart implementation,”
IEEEJournal of Selected Topics in Signal Processing , vol. 7, no. 3, pp. 472–483, June 2013.[16] Y. Ge, F. Wen, H. Kim, M. Zhu, S. Kim, L. Svensson, and H. Wymeer-sch, “5G SLAM using the clustering and assignment approach withdiffuse multipath,”
Sensors , vol. 20, 4656.[17] T. Kurien, “Issues in the design of practical multitarget tracking al-gorithms,” in
Multitarget-Multisensor Tracking: Advanced Applications ,Y. Bar-Shalom, Ed. Artech House, 1990.[18] S. S¨arkk¨a,
Bayesian Filtering and Smoothing . Cambridge UniversityPress, 2013.[19] K. Granstr¨om, A. Natale, P. Braca, G. Ludeno, and F. Serafino, “GammaGaussian inverse Wishart probability hypothesis density for extendedtarget tracking using X-band marine radar data,”
IEEE Transactions onGeoscience and Remote Sensing , vol. 53, no. 12, pp. 6617–6631, Dec.2015.[20] Y. Xia, K. Granstr¨om, L. Svensson, A. F. Garc´ıa-Fern´andez, and J. L.Williams, “Extended target Poisson multi-Bernoulli mixture trackersbased on sets of trajectories,” in
Proceedings of the 22nd InternationalConference on Information Fusion , 2019.[21] T. J. Broida, S. Chandrashekhar, and R. Chellappa, “Recursive 3-D mo-tion estimation from a monocular image sequence,”
IEEE Transactionson Aerospace and Electronic Systems , vol. 26, no. 4, pp. 639–656, Jul.1990.[22] L. Hammarstrand, L. Svensson, F. Sandblom, and J. Sorstedt, “Extendedobject tracking using a radar resolution model,”
IEEE Transactions onAerospace and Electronic Systems , vol. 48, no. 3, pp. 2371–2386, Jul.2012.[23] B. Ristic and J. Sherrah, “Bernoulli filter for joint detection and trackingof an extended object in clutter,”
IET Radar, Sonar Navigation , vol. 7,no. 1, pp. 26–35, Jan. 2013.[24] X. Shen, Z. Song, H. Fan, and Q. Fu, “General Bernoulli filter for arbi-trary clutter and target measurement processes,”
IEEE Signal ProcessingLetters , vol. 25, no. 10, pp. 1525–1529, Oct. 2018.[25] E. Mazor, A. Averbuch, Y. Bar-Shalom, and J. Dayan, “Interactingmultiple model methods in target tracking: a survey,”
IEEE Transactionson Aerospace and Electronic Systems , vol. 34, no. 1, pp. 103–123, Jan.1998.[26] M. Feldmann, D. Fr¨anken, and W. Koch, “Tracking of extended objectsand group targets using random matrices,”
IEEE Transactions on SignalProcessing , vol. 59, no. 4, pp. 1409–1420, Apr. 2011.[27] A. K. Gupta and D. K. Nagar,
Matrix Variate Distributions . Chapman& Hall, 1999. [28] K. Granstr¨om and U. Orguner, “Estimation and maintenance of measure-ment rates for multiple extended target tracking,” in , 2012, pp. 2170–2176.[29] ——, “On the reduction of Gaussian inverse Wishart mixtures,” in , 2012, pp. 2162–2169.[30] J. L. Williams, “An efficient, variational approximation of the best fittingmulti-Bernoulli filter,”
IEEE Transactions on Signal Processing , vol. 63,no. 1, pp. 258–273, Jan. 2015.[31] A. F. Garc´ıa-Fern´andez, L. Svensson, J. L. Williams, Y. Xia, andK. Granstr¨om, “Trajectory Poisson multi-Bernoulli filters,”
IEEE Trans-actions on Signal Processing , vol. 68, pp. 4933–4945, 2020.[32] C. M. Bishop,
Pattern Recognition and Machine Learning . Springer,2006.[33] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “A density-based algorithmfor discovering clusters in large spatial datasets with noise,” in ,1996, pp. 226–231.[34] Y. Xia, K. Granstr¨om, L. Svensson, and M. Fatemi, “Extendedtarget Poisson multi-Bernoulli filter,” 2018. [Online]. Available:https://arxiv.org/abs/1801.01353[35] K. G. Murty, “An algorithm for ranking all the assignments in orderof increasing cost.”
Operations Research , vol. 16, no. 3, pp. 682–687,1968.[36] A. F. Garc´ıa-Fern´andez and S. Maskell, “Continuous-discrete multipletarget filtering: PMBM, PHD and CPHD filter implementations,”
IEEETransactions on Signal Processing , vol. 68, pp. 1300–1314, 2020.[37] F. Meyer and J. L. Williams, “Scalable detection and tracking ofextended objects,” in
IEEE International Conference on Acoustics,Speech and Signal Processing , 2020, pp. 8916–8920.[38] A. S. Rahmathullah, A. F. Garc´ıa-Fern´andez, and L. Svensson, “Gen-eralized optimal sub-pattern assignment metric,” in , 2017, pp. 1–8.[39] S. Yang, M. Baum, and K. Granstr¨om, “Metrics for performanceevaluation of elliptic extended object tracking methods,” in
IEEEInternational Conference on Multisensor Fusion and Integration forIntelligent Systems , 2016, pp. 523–528.[40] A. F. Garc´ıa-Fern´andez, L. Svensson, and M. R. Morelande, “Multipletarget tracking based on sets of trajectories,”
IEEE Transactions onAerospace and Electronic Systems , vol. 56, no. 3, pp. 1685–1707, Jun.2020.[41] Y. Xia, K. Granstr¨om, L. Svensson, A. F. Garc´ıa-Fern´andez, and J. L.Wlliams, “Multi-scan implementation of the trajectory Poisson multi-Bernoulli mixture filter,”
Journal of Advances in Information Fusion ,vol. 14, no. 2, pp. 213–235, Dec. 2019. upplementary material: A Poissonmulti-Bernoulli mixture filter for co-existing point and extended targets A PPENDIX
AIn this appendix, we prove Theorem 1, which provides theupdate step, making use of probability generating functionals(PGFLs). A PGFL is an alternative representation of a multi-object density, in the same way as Fourier and z -transformsare for signals defined in the time domain.For a multi-object density f ( · ) , its PGFL G f [ · ] is given bythe set integral [3] G f [ h ] = (cid:90) h X f ( X ) δX (69)where h ( · ) is a unitless function of state space, and h X = (cid:81) x ∈ X h ( x ) , h ∅ = 1 . The test function for PGFLs related todensities defined for targets and measurements are denoted as h ( · ) and g ( · ) , respectively.Given the PGFL G f [ · ] , we can recover its multi-objectdensity f ( · ) by the set derivative [3] f ( X ) = δδX G f [ h ] (cid:12)(cid:12)(cid:12)(cid:12) h =0 . (70) A. PGFLs of targets and measurements
The density (1) in PGFL form is represented as [9] G k | k (cid:48) [ h ] = G p k | k (cid:48) [ h ] · G mbm k | k (cid:48) [ h ] (71) G p k | k (cid:48) [ h ] = exp (cid:0) (cid:104) λ k | k (cid:48) , h − (cid:105) (cid:1) ∝ exp (cid:0) (cid:104) λ k | k (cid:48) , h (cid:105) (cid:1) (72) G mbm k | k (cid:48) [ h ] = (cid:88) a ∈A k (cid:48)| k w ak | k (cid:48) n k (cid:48)| k (cid:89) i =1 G i,a i k | k (cid:48) [ h ] (73) ∝ (cid:88) a ∈A k (cid:48)| k n k | k (cid:48) (cid:89) i =1 (cid:104) w i,a i k | k (cid:48) G i,a i k | k (cid:48) [ h ] (cid:105) (74)hwhere G i,a i k | k (cid:48) [ h ] = 1 − r i,a i k | k (cid:48) + r i,a i k | k (cid:48) (cid:104) f i,a i k | k (cid:48) , h (cid:105) . (75)Given the multi-target state X , measurements from eachtarget are independent, and there is also independent PPPclutter. Therefore, the PGFL G Z [ g | X ] of the measurementsgiven X is the product of PGFL G Z [ g | X ] = exp (cid:0) (cid:104) λ C , g − (cid:105) (cid:1) (cid:89) x ∈ X G [ g | x ] (76)where G [ g | x ] is the PGFL of f ( Z | x ) . B. Joint PGFL of targets and measurements
The joint PGFL of measurements and targets is [3], [9] F [ g, h ] = (cid:90) (cid:90) g Z k h X k f ( Z k | X k ) f k | k − ( X k ) δZ k δX k (77) = (cid:90) G [ g | X k ] h X k f k | k − ( X k ) δX k (78) = exp (cid:0) (cid:104) λ C , g − (cid:105) (cid:1) G k | k − [ hG Z [ g |· ]] (79) ∝ exp (cid:0) (cid:104) λ C , g (cid:105) + (cid:104) λ k | k − , hG [ g |· ] (cid:105) (cid:1) × (cid:88) a ∈A k | k − n k | k − (cid:89) i =1 (cid:104) w i,a i k | k − G i,a i k | k − (cid:2) hG [ g |· ] (cid:3)(cid:105) . (80)We denote the first line of (80) as F [ g, h ] = exp (cid:0) (cid:104) λ C , g (cid:105) + (cid:10) λ k | k − , hG [ g |· ] (cid:11)(cid:1) (81)which represents the joint PGFL of measurements (includingfalse alarms) and targets in the PPP, up to a proportionalityconstant. We also denote F i,a i [ g, h ] = G i,a i k | k − (cid:2) hG [ g |· ] (cid:3) (82) = 1 − r i,a i k | k − + r i,a i k | k − (cid:28) f i,a i k | k − , hG [ g |· ] (cid:29) (83)which represents the joint PGFL of measurements (not includ-ing false alarms) and the i -th potential target. Then, using (5),we can write (80) as F [ g, h ] ∝ F [ g, h ] (cid:88) a ∈A k | k − n k | k − (cid:89) i =1 (cid:104) w i,ak | k − F i,a i [ g, h ] (cid:105) . (84) C. Updated PGFL
We calculate the updated density f k | k ( · ) via its PGFL G k | k [ h ] , which is given by the set derivative of F [ g, h ] w.r.t. Z k evaluated at g = 0 [3, Sec. 5.8] [9, Eq. (25)] G k | k [ h ] ∝ δδZ k F [ g, h ] (cid:12)(cid:12)(cid:12)(cid:12) g =0 . (85)Applying the product rule [9, Eq. (31)], we obtain G k | k [ h ] ∝ (cid:88) W (cid:93)···(cid:93) W nk | k − = Z k δδW F [ g, h ] × (cid:88) a ∈A k | k − n k | k − (cid:89) i =1 δδW i (cid:16) w i,a i k | k − F i,a i [ g, h ] (cid:17) (cid:12)(cid:12)(cid:12)(cid:12) g =0 . (86)The sum in (86) is over all decompositions of the measurementset Z into ( n k | k − + 1) subsets, where one subset, W ,represents measurements which are either false alarms, orcorrespond to a target represented by the PPP component (i.e.,a target which has never been detected so far), and subset W i , i > , represents measurements assigned to the i -th Bernoulli.We develop the required set derivatives over the followinglemmas, starting with the Bernoulli component F i [ g, h ] inSection A-C1, and then moving on to the update of the PPP, F [ g, h ] , in Section A-C2.
1) Bernoulli update:
We calculate the set derivative of F i,a i [ g, h ] . For W i (cid:54) = ∅ , we have δδW i F i,a i [ g, h ] = r i,a i k | k − (cid:68) f i,a i k | k − , h δδW i G [ g |· ] (cid:69) (87) δδW i F i,a i [ g, h ] (cid:12)(cid:12)(cid:12)(cid:12) g =0 = r i,a i k | k − (cid:68) f i,a i k | k − , hf ( W i |· ) (cid:69) (88)here we have applied (70). For W i = ∅ , the set derivativedoes not change F i,a i [ g, h ] .Equation (88) and F i,a i [0 , h ] are the PGFL of a weightedBernoulli component [9, Lem. 2] with parameters given in thefollowing lemma. Lemma 3.
The update of the weighted PGFL component w i,a i k | k − G i,a i k | k − [ h ] (weighted Bernoulli) with measurement set W i , i.e., w i,a i ,W i k | k G i,a i ,W i k | k [ h ] = δδW i (cid:16) w i,a i k | k − F i,a i [ g, h ] (cid:17) (cid:12)(cid:12) g =0 is the PGFL of a weighted Bernoulli distribution, i.e., adistribution of the form [9, Lem. 2] f i,a i ,W i k | k ( X ) = w i,a i ,W i k | k × − r i,a i ,W i k | k X = ∅ r i,a i ,W i k | k f i,a i ,W i k | k ( x ) X = { x } | X | > (89) where for W i = ∅ , w i,a i ,W i k | k = w i,a i k | k − (cid:104) − r i,a i k | k − + r i,a i k | k − l i,a i ,W i k | k (cid:105) (90) l i,a i ,W i k | k = (cid:10) f i,a i k | k − , (cid:10) f i,a i k | k − , f ( ∅|· ) (cid:11) (91) r i,a i ,W i k | k = r i,a i k | k − l i,a i ,W i k | k − r i,a i k | k − + r i,a i k | k − l i,a i ,W i k | k (92) f i,a i ,W i k | k ( x ) = f i,a i k | k − ( x ) f ( ∅| x ) l i,a i ,W i k | k . (93) For | W i | ≥ , w i,a i ,W i k | k = w i,a i k | k − r i,a i k | k − l i,a i ,W i k | k (94) l i,a i ,W i k | k = (cid:10) f i,a i k | k − , f ( W i |· ) (cid:11) (95) r i,a i ,W i k | k = 1 (96) f i,a i ,W i k | k ( x ) = f i,a i k | k − ( x ) f ( W i | x ) l i,a i ,W i k | k . (97)This lemma therefore proves how to update a previousBernoulli with a misdetection or a detection hypothesis inTheorem 1.
2) PPP update:
We now turn to calculating the update ofthe PPP in (86) via the set derivative of F [ g, h ] , see (81). Lemma 4.
The set derivative of F [ g, h ] is δδW F [ g, h ] = F [ g, h ] (cid:88) P ∠ W (cid:89) V ∈ P d V [ g, h ] (98) where d V [ g, h ] = δδV (cid:0)(cid:10) λ C , g (cid:105) + (cid:10) λ k | k − , hG [ g |· ] (cid:11)(cid:1) (99) and (cid:80) P ∠ W denotes the sum over all partitions P of W . (cid:3) The proof of Lemma 4 is in Section A-D.Following (86) , we evaluate the first factor of (98), F [ g, h ] ,at g = 0 to obtain F [0 , h ] = exp (cid:0)(cid:10) λ k | k − , hf ( ∅|· ) (cid:11)(cid:1) . (100) This is proportional to the PGFL of a PPP with intensity λ k | k ( x ) = f ( ∅| x ) λ k | k − ( x ) , which proves (7).We now need to compute the set derivatives in (99), evaluatethem at g = 0 and compute the corresponding multi-objectdensities. For V = { v } (set with a single element), we have d { v } [ g, h ] = λ C ( v ) + (cid:10) λ k | k − , h δδ { v } G [ g |· ] (cid:11) (101) d { v } [0 , h ] = λ C ( v ) + (cid:10) λ k | k − , hf ( { v } |· ) (cid:11) (102)where we have applied the linear rule [3].For | V | > , we have d V [ g, h ] = (cid:10) λ k | k − , h δδV G [ g |· ] (cid:11) (103) d V [0 , h ] = (cid:10) λ k | k − , hf ( V |· ) (cid:11) . (104)where we have applied that the set derivative of a constant iszero. This is why the term λ C ( v ) is not present in (103).The following lemma provides the form of the multi-objectdensities whose PGFL is d V [0 , h ] in (102) and (104). Lemma 5.
The update of the PGFL of the PPP prior withmeasurement subset Vw Vk | k G Vk | k [ h ] = d V [ g, h ] (cid:12)(cid:12)(cid:12) g =0 are PGFLs of weighted Bernoulli distributions with the form[9, Lem. 2] f Vk | k ( X ) = w Vk | k × − r Vk | k X = ∅ r Vk | k f Vk | k ( x ) X = { x } | X | > (105) where w Vk | k = (cid:34) δ [ | V | ] (cid:89) z ∈ V λ C ( z ) (cid:35) + l Vk | k (106) l Vk | k = (cid:28) λ k | k − , f ( V |· ) (cid:29) (107) r Vk | k = l Vk | k w Vk | k (108) f Vk | k ( x ) = f ( V | x ) λ k | k − ( x ) l Vk | k . (cid:3) (109)Therefore, the PGFL of the updated PPP in (98) correspondsto the union of a PPP for undetected targets, with intensity λ k | k ( x ) = f ( ∅| x ) λ k | k − ( x ) and, a multi-Bernoulli mixturewhere each term in the mixture is a partition of W and eachBernoulli component has a weight and density provided inLemma 98. This concludes the proof of Theorem 1.It should be noted that to perform the PMBM updatewe first take all possible sets W (cid:93) · · · (cid:93) W n k | k − = Z k ,which represents subsets of Z k associated to the PPP ( W )or the previous Bernoullis ( W i , i > ). Then, we take allpossible partitions P of W , P ∠ W , to generate the newBernoulli components. A compact way to represent thesedecompositions is to take all possible subsets of Z k to generatethe new Bernoulli components and represent the possible dataassociations to previous Bernoulli components, as in done inTheorem 1. . Set derivative of F [ g, h ] We prove Lemma 4 by induction. The set derivative of F [ g, h ] , see (81), with respect to a set with | W | = 1 isstraightforward, as there is a single partitioning of a oneelement set. For induction, we assume that the lemma holdsup to a given size | W | , and we show that it holds for ˜ W = W ∪ { z } : δδ ˜ W F [ g, h ]= δδ { z } δδW F [ g, h ] (110) = δδ { z } (cid:32) F [ g, h ] (cid:88) P ∠ W (cid:89) V ∈ P d V [ g, h ] (cid:33) (111) = (cid:18) δδ { z } F [ g, h ] (cid:19) (cid:88) P ∠ W (cid:89) V ∈ P d V [ g, h ]+ F [ g, h ] (cid:88) P ∠ W (cid:32) δδ { z } (cid:89) V ∈ P d V [ g, h ] (cid:33) (112) = F [ g, h ] d { z } [ g, h ] (cid:88) P ∠ W (cid:89) V ∈ P d V [ g, h ]+ F [ g, h ] (cid:88) P ∠ W (cid:88) V ∈ P (cid:18) δδ { z } d V [ g, h ] (cid:19) (cid:89) V (cid:48) ∈ P \{ V } d V (cid:48) [ g, h ] . (113)Each step in the previous derivation results from the productrule [3].From (99), we have ∂∂ { z } d V [ g, h ] = d V ∪{ z } [ g, h ] . In addi-tion, each partitioning of ˜ W consists of either a partitioningof W with an additional single element subset { z } ; or apartitioning of W , adding element z to one of the existingsubsets [3, App. D.2]. Since the top line in (113) handles theformer case and the bottom line handles the latter case, we findthat (113) is equivalent to F [ g, h ] (cid:80) P ∠ ˜ W (cid:81) V ∈ P d V [ g, h ] ,which proves Lemma 4.A PPENDIX
BIn this appendix, we explicitly relate the space of coexistingpoint extended targets, X = R n x (cid:93)X e , in Section IV to spacesused in interacting multiple models [25], which usually includea class variable to distinguish between different models. Given x ∈ R n x (cid:93) X e , we know if x represents a point target or anextended target as R n x and X e are disjoint. Therefore, it isnot necessary to extend the single target space with a classvariable to distinguish both types of targets.Nevertheless, it is possible to add a class variable c suchthat the single target state becomes y = ( c, x ) , where c = 0 for point targets and c = 1 to extended targets. In this case,the single target space is ( { } × R n x ) (cid:93) ( { } × X e ) and thePMBM filtering recursions remain unchanged.A PPENDIX
CFor completeness, in this appendix, we provide the (approxi-mate) single extended target update for factorised GGIW priors[11], [26]. The resulting expressions are provided in Table II.The update for the parameters of the Gamma distribution isexact due to the Poisson-Gamma conjugacy.
TABLE II: Update and marginal likelihood of a GGIW density
Input:
Prior GGIW parameters ζ + = ( α + , β + , x + , P + , v + , V + ) , set W ofmeasurements. Output: ( ζ, (cid:96) ) = u e ( ζ + , W ) , where ζ are the updated GGIW parametersand (cid:96) the marginal likelihood evaluated at W .If | W | > ζ = α = α + + | W | β = β + + 1 x = x + + KεP = P + − KHP + v = v + + | W | V = V + + N + Z where z = 1 | W | (cid:88) z ∈ W zZ = (cid:88) z ∈ W ( z − z ) ( z − z ) T ˆ X = V + ( v + − d − − ε = z − Hx + S = HP + H T + ˆ X | W | K = P + H T S − N = ˆ X / S − / εε T S − T/ ˆ X T/ (cid:96) = (cid:16) π | W | | W | (cid:17) − d/ | V + | v + − d − Γ d (cid:16) v − d − (cid:17) (cid:12)(cid:12)(cid:12) ˆ X (cid:12)(cid:12)(cid:12) / Γ ( α ) ( β + ) α + | V | v − d − Γ d (cid:16) v + − d − (cid:17) | S | / Γ ( α + ) ( β ) α . If | W | = 0 ζ = ( α + , β + + 1 , x + , P + , v + , V + ) (cid:96) = (cid:18) β + β + + 1 (cid:19) α ++