Decentralized Poisson Multi-Bernoulli Filtering for Vehicle Tracking

1 Decentralized Poisson Multi-Bernoulli Filtering for V ehicle T racking Markus Fr ¨ ohle, Karl Granstr ¨ om, Henk W ymeersch Abstract —A decentralized Poisson multi-Ber noulli filter is proposed to track multiple vehicles using multiple high-r esolution sensors. Independent filters estimate the vehicles’ presence, state, and shape using a Gaussian process extent model; a decentralized filter is realized thr ough fusion of the filters posterior densities. An efficient implementation is achieved by parametric state repr esentation, utilization of single hypothesis tracks, and fusion of vehicle inf ormation based on a fusion mapping. Numerical results demonstrate the performance. Index T erms —Gaussian processes, multitarget tracking, poste- rior fusion, tar get extent. I . I N T RO D U C T I O N Multitarget tracking (MTT), i.e., tracking of independently moving tar gets, is important for surveillance and safety ap- plications [1], [2]. T raditionally , it has been developed for surveillance of the sk y using ground-to-air radar sensors. An MTT filter allo ws to incorporate the peculiarities of those kind of sensors: false alarm measurements due to clutter; missed detections; unknown measurement-to-target correspondence; and target appearance and disappearance; which are all chal- lenges that arise for radar -like sensors. In many typical MTT scenarios, the sensor resolution is low with respect to (w .r .t.) the target size, and a reasonable assumption is to model the targets as points having a kinematic state (e.g., position and velocity). An elegant way to track multiple targets is via the Poisson multi-Bernoulli mixture (PMBM) filter [3], [4], which preserves a PMBM form during prediction and update steps. W ith the av ailability of high resolution sensors, the point target assumption does not hold anymore, see, e.g., [5]. For instance, a high resolution Lidar sensor can obtain in one scan multiple detections from a single v ehicle [6] or cyclist [7]. This is because the sensor resolution is high w .r .t. the tar get (here, a vehicle) size. In such an application scenario, the vehicle extent needs to be modeled (and estimated) as well in the MTT filter , leading to an extended target tracking (ETT) filter . In ETT, extended targets (ETs) giv e rise to possibly multiple noisy detections, the vehicle extent (shape and size) is a priori unknown and may vary over time, and the objectiv e is to estimate the ET’ s kinematic state as well as its extent [5]. M. Fr ¨ ohle was with the Department of Electrical Engineering, Chalmers Univ ersity of T echnology , Gothenburg and is now with Zenuity AB, Sweden. E-mail: markus.frohle@zenuity.com . K. Granstr ¨ om, and H. W ymeer- sch are with the Department of Electrical Engineering, Chalmers Univ er- sity of T echnology , Gothenburg, Sweden. E-mail: { karl.granstrom, henkw } @chalmers.se . This work was supported, in part, by the EU- H2020 project HIGHTS (High Precision Positioning for Cooperative ITS Applications) under grant no. MG-3.5a-2014-636537 and COPPLAR (campus shuttle cooperativ e perception and planning platform) project funded under grant no. 2015-04849 from Vinno va. A common model for the ET extent is based on Gaussian process (GP) modelling [8], [9]. Neither MTT nor ETT filters are bound to single sensors. When multiple sensors are used, it becomes necessary to fuse the information from the different sensors. Centralized fusion consists of transmitting all measurements to a central process- ing unit. Decentralized fusion consists of local processing (i.e., perform recursive state-space estimation) at each sensor unit, and transmission of the tracking results. Such an approach was proposed in [10], for fusing two multi-Bernoulli (MB) densities, though it did not account for a Poisson point process (PPP) component. In this paper , we utilize and extend the shape model from [8], [9], and integrate it into an ETT filter, which allows tracking of multiple ETs. The resulting ETT filter’ s multiobject posterior is of the so-called Poisson multi-Bernoulli (PMB) form [11]. Furthermore, we propose a nov el fusion strategy that performs fusion separately for the PPP and MB parts of the PMB density . The implementation of the PMB filter yields a tracking filter with low computational cost 1 locally at each sensor , and globally lo w computational cost through the introduction of a fusion map based on the Kullback-Leibler div ergence (KLD) between target tracks. Simulation results demonstrate the performance of the proposed independent ETT filter as well as of the decentralized ETT filtering approach. The main contributions of this paper are: • W e apply a state-of-the-art ETT filter [3], [4] and propose a novel distributed fusion strategy based on the Kullback- Leibler av erage (KLA) by fusion of the filters’ posterior multiobject densities of PMB form; and • W e enable lo w complexity in distributed fusion through introduction of a fusion map based on KLD between target tracks. • W e extend, in Section VII-A2, the GP model from [8], [9] for the ET shape description to a multisensor scenario incorporating the sensors’ state (position, orientation); The remainder of this paper is organized as follows; Sec- tion II presents relev ant related work to this paper’ s work, Section III gives some background knowledge on random finite sets (RFSs), and Section IV introduces the system model and the problem formulation. Section V details the proposed ETT filter , Section VI presents the decentralized posterior fusion approach using independent ETT filters. Simulation results are gi ven in Section VII, and conclusions are drawn in Section VIII. 1 As measured by the a verage time it takes for one cycle of prediction and update. 2 I I . R E L A T E D W O R K Here, we discuss related work rele vant for this paper’ s work co vering extent modeling for tracking an ET, ETT, and information fusion with multiple sensors. A. Extent modelling Different models for the target extent (shape and size) exist in ETT, which may be classed according to complexity , ranging from assuming a specific geometric shape [6], [7], [12] with, e.g., unkno wn translation and rotation, to models that describe general shapes, see, e.g., [13]–[15]. T ypically , more complex models provide a richer shape description. T wo popular models are the random matrix approach [12], where the target shape is described by an ellipsoid; and a GP based approach [8], [9], where a star con vex tar get shape is described by a GP. See [5] for an extensi v e overvie w of w orks on ETT. B. T rac king multiple extended tar gets For tracking multiple extended targets, sev eral different filters have been dev eloped: Probability Hypothesis Density (PHD) filter [16]–[19]; Cardinalized PHD (CPHD) filter [20]; δ -generalized labelled multi-Bernoulli ( δ -GLMB) filter [21]; and PMBM filter [3], [4]. Multiple extended target tracking filters have been applied to tracking of different tar get types (e.g., for cars [22]–[24]). The δ -GLMB filter and PMBM filter are so called multi- object conjugate priors, meaning that if we start with the conjugate density form ( δ -GLMB or PMBM), then all sub- sequent predicted and updated densities will be of the same form. Based on the conjugate priors, computationally cheaper , approximate filters have been presented. The labeled multi- Bernoulli (LMB) filter is an approximation of the δ -GLMB filter , see [21]. The PMB filter is an approximation of the PMBM filter , see [11], [25]. The PMBM conjugate prior [26] was originally developed for point targets; a PMBM conjugate prior for extended targets was presented in [3], [4]. In several simulation studies it has been shown that, compared to tracking filters b uilt upon la- belled RFS, the PMBM conjugate prior has good performance for tracking the set of present target states, for both point targets [27]–[29] and extended tar gets [3], [4], [30]. The PMBM conjugate priors for point targets and e xtended targets have been sho wn to be versatile, and ha ve been used with data from Lidars [31]–[34], radars [32], [33], and cameras [33], [35]. They hav e been successfully applied not only to tracking of moving targets, but also mapping of stationary objects [36], as well as joint tracking and sensor localisation [37]. Thus, it is well motiv ated to use the PMBM filter in this work; specifically , we work with the computationally cheaper PMB filter . Dev eloping similar decentralized tracking algorithms for the δ -GLMB or LMB filter [21] is a topic for future work. C. T rac king using multiple sensors Different approaches to incorporate measurements from multiple sensors which were acquired within one scan exist, where the simplest may be seen as performing multiple update steps (e.g., Kalman update or likewise) by augmenting the measurement model to incorporate all sensor models, see, e.g., [38]. When sensors are geographically separated, measurements from all sensors need to be transmitted to the central pro- cessing unit where the filter is run. In the absence of such a central unit or simply due to the limited capacity of the communication channel, one can perform filtering already at the sensor and share only target track information, e.g., the parameters of a known probability density function (PDF) family . Through incorporation of target track information ob- tained by independent filters, a decentralized filtering approach (i.e., information fusion) can be realized without the need for additional prior information. Such methods need to ensure that (unknown) common information obtained by the independent tracking filters is not double counted, e.g., the prior target density [39]. Depending on the type of posterior multiobject density in MTT/ETT filtering, several sub-optimal information fusion strategies have been developed based on, e.g, covariance inter- section (CI) for Gaussian densities. This replaces the product form of Bayes’ rule with the KLA 2 . CI (and consequently KLA) is a method to fuse information with unknown priors in a robust way in the sense that the fused posterior is con- servati v e and nev er ov erconfident about the estimates and thus implicitly sub-optimal [39], [40]. Examples of CI/KLA include fusion of Bernoulli and independent identically distributed (IID) cluster processes posteriors [41], fusion for Bernoulli filters [42], fusion for PHD and CPHD filters [43]–[48], and fusion for LMB filters [10], [47], [49]. I I I . B AC K G RO U N D O N R A N D O M F I N I T E S E T S T wo types of RFSs relev ant for this work deserve special attention: Bernoulli RFS, and a PPP. They can be extended to MB RFS, multi-Bernoulli mixture (MBM) RFS and combined in a PMBM RFS. These are described belo w , for more details on RFSs the reader is referred to [50]. A. Common RFS densities Bernoulli RFS: A Bernoulli RFS X has a multiobject density [50] f ( X ) =      1 − r , X = ∅ , r f ( x ) , X = { x } , 0 , | X | ≥ 2 , (1) where r ∈ [0 , 1] denotes the probability that a target exists and if it exists f ( x ) is its PDF. Multi-Bernoulli RFS: An MB RFS is the disjoint union of independent Bernoulli RFSs indexed by i . It is fully parametrized by { r i , f i ( x ) } i ∈ I , where I is its index set. For X = { x 1 , . . . , x n } , the multiobject density can be written as f ( X ) = X ] i ∈ I X i = X Y i ∈ I f i ( X i ) (2) 2 In some literature this is known as exponential mixture density (EMD). 3 for | X | ≤ | I | , and f ( X ) = 0 otherwise. The notation, X 1 ] X 2 = X means X 1 ∪ X 2 = X and X 1 ∩ X 2 = ∅ . Multi-Bernoulli Mixtur e RFS: The multiobject density of an MBM is the normalized, weighted sum of multiobject densities of MBs, which can be stated as [3] f ( X ) = X j ∈ J w j X ] i ∈ I j X i = X Y i ∈ I j f j,i ( X i ) . (3) The MBM multiobject density is parametrized by { w j,i , { r j,i , f j,i ( x ) } i ∈ I j } j ∈ J , where w j is the weight of MB j , and J is the index set of the MBs in the MBM. An MB is therefore a special case of an MBM with | J | = 1 . PPP: A PPP is a type of RFS, where the cardinality follows a Poisson distribution, and its elements are IID. It is parametrized by the intensity function D ( x ) = λf ( x ) , where λ > 0 is the Poisson rate and f ( x ) is a PDF on the single element state x . The multiobject density of a PPP is [3], [50] f ( X ) = e − λ n Y i =1 λf ( x i ) . (4) PMBM: An PMBM RFS is the disjoint set union of an PPP and an MBM having multiobject density [3] f ( X ) = X X u ] X d = X f u ( X u ) f d ( X d ) , (5) where f u ( · ) has multiobject density (4), and f d ( X d ) has multiobject density (3). For | J | = 1 , (5) is called a PMB distribution. B. RFS Bayesian F ilter Similar to the random vector (R V) case, an RFS based filter can be described, conceptually at least, within the Bayesian filtering frame work by performing a prediction step using the motion model [2, Ch. 14] f + ( X ) = Z f ( X | X 0 ) f − ( X 0 ) δ X 0 , (6) where f − ( X 0 ) is the prior RFS density , f ( X | X 0 ) is the multiobject process model, and a Bayesian update step f ( X | Z ) ∝ ` ( Z | X ) f + ( X ) . (7) Here, f + ( X ) is the predicted RFS density , and ` ( Z | X ) is the RFS measurement likelihood for measurement set Z . A typical way to estimate the set states from a Bernoulli process with RFS density f ( X ) is by comparing the proba- bility of existence r to an existence threshold r th . For r > r th , the target is said to e xist and has PDF f ( x ) . Its state can then be estimated by the mean ˆ x = R x f ( x )d x . See, e.g., [28] for a elongated discussion on multiobject estimation. I V . P R O B L E M F O R M U L A T I O N A N D S Y S T E M M O D E L Here, we present first the problem formulation of this paper followed by the ET state and transition model and the ET set measurement likelihood function. Inter - sensor communication Fig. 1. Scenario with three ETs observed by two sensors with partially overlapping FoVs. A. Pr oblem F ormulation W e consider a scenario with N sens sensing systems (each composed of sensors plus filter), each collecting measure- ments using its local sensors, with the aim to jointly surveil an en vironment E where vehicles pass (c.f. Fig. 1). Our goals are (i) to derive a low-comple xity PMB-ETT filter for each sensing system s = 1 , 2 , . . . , N sens , computing in ev ery time step k the posterior density f s ( X k | Z s,k :1 ) of the ETs, using only its own sensors with measurement set Z s,k :1 ; and (ii) to deri ve a decentralized method to combine posterior information of ETs obtained by N sens independent ETT filters in order to obtain a global posterior density ¯ f w ( X k | Z 1 ,k :1 , Z 2 ,k :1 , . . . , Z N sens ,k :1 ) (i.e., using only the posterior densities of each ETT filter). B. ET State and T ransition Model A standard motion model for the ET is assumed, where ET motion follows IID Markov processes with single ET transition PDF f k +1 | k ( x k +1 | x k ) , where x k denotes the state at time k . T argets arriv e according to a non-homogeneous PPP with intensity D b ( x k ) , and depart according to IID Markov processes, where the surviv al probability in x k is p S ( x k ) . Additionally to this, x k comprises an unknown Poisson rate γ k describing the a verage number of measurements generated by the ET and the ET spatial state y k (this will become clear in Sec. VII-A where a detailed description of the state and transition density is provided). C. ET Set Measur ement Likelihood Function In one scan, a sensor receives a set of measur ements Z consisting of target-generated measurements z and clutter, where the ETs are independently detected with state-dependent probability of detection p D ( x ) , which depends on the sensor FoV. Clutter is modeled by a PPP with intensity κ ( z ) = λc ( z ) with mean λ and spatial distribution c ( z ) . T ar get- generated measurements are modeled by a PPP with intensity γ ( x ) f ( z | x ) , where both the Poisson measurement rate γ ( x ) and the single measurement likelihood f ( z | x ) are state de- 4 pendent. The measurement likelihood for ETs { x 1 , . . . , x n } and measurement set Z is ` ( Z |{ x 1 , . . . , x n } ) = e − λ X Z c ] Z 1 ] ... Z n = Z [ c ( · )] Z c n Y i =1 ` Z i ( x i ) , (8) where [ c ( · )] Z c is shorthand for Q z ∈ Z c c ( z ) , [ c ( · )] ∅ = 1 by definition, and [3] ` Z ( x ) = ( p D ( x ) e − γ ( x ) Q z ∈ Z γ ( x ) f ( z | x ) , | Z | > 0 , (1 − p D ( x )) + p D ( x ) e − γ ( x ) , Z = ∅ . (9) Note that (i) equation (8) in volv es potentially multiple ETs, leading to a data association (D A) problem; (ii) a single ET can generate multiple measurements. V . I N D E P E N D E N T P M B - E T T F I LTE R In this section, we briefly describe the processing performed by each PMB-ETT (PMB-ETT) filter from [3]. W e omit the time index k for brevity . The PMB model is a combination of a PPP describing the distrib ution of unkno wn tar gets, i.e., targets which are hypothesized to e xist, but have not yet been detected; and a MB which describes targets that hav e been detected at least once. The target set can therefore be split into two disjoint subsets X = X u ] X d corresponding to the unknown target set X u (with PPP intensity D u ( x ) , modelled as a non-normalized mixture density with mixture components located in the region of interest, denoted E ) and the detected target set X d with density f d ( X d ) and index set I . Hence, the PMB density is fully described by a MB component described by { r i , f i ( x ) } i ∈ I and a PPP component D u ( x ) . 1) PMB-ETT F ilter Pr ediction: The predicted density is a PMB density with parameters [3, Sec. IV], [4] D u + = D b ( x ) +  D u , p S f k +1 | k  , (10) r i + =  f i , p S  r i , (11) f i + =  f i , p S f k +1 | k  h f i , p S i , (12) where h g , h i = R g ( x ) f ( x )d x denotes the inner product. The proof of the prediction step can be found in, e.g., [11]. 2) PMB-ETT F ilter Update: W e introduce the set of valid D As A = P ( M ∪ I ) , where M is the index set for Z . Here, A ∈ A is a partition of M ∪ I into non-empty disjoint subsets C ∈ A (called index cells), with the constraint that for each C : | C ∩ I | ≤ 1 (i.e., measurements can only be associated with a single target). When | C ∩ I | = 1 , let the entry in C ∩ I be denoted by i C and let C C = ∪ m ∈ C ∩ M z m contain the associated measurements. Giv en the predicted prior PMB density with parameters (10), (11), (12), and a set of measurements Z ; the updated density is a PMBM density [3, Sec. IV], [4] f ( X | Z ) = X X u ] X d = X f u ( X u ) X A ∈A w A f d A ( X d ) , (13) f u ( X u ) = e −h D u , 1 i Y x ∈ X u D u ( x ) , (14) f d A ( X d ) = X ] C ∈ A X C = X d Y C ∈ A f C ( X C ) , (15) D u ( x ) = q D ( x ) D u + ( x ) , (16) where f C ( X C ) is a Bernoulli density , with existence probabil- ity and spatial distribution provided in Appendix A, together with the weights w A of each DA hypothesis. Abov e, q D ( x ) denotes the probability that the target x is not detected and is defined as q D ( x ) = 1 − p D ( x ) + p D ( x ) e − γ ( x ) . (17) T o reduce computational complexity , we use standard methods and truncate the space of possible partitions by clustering measurements and consider DA w .r .t. different clusters [17], [18], [20]. Finally , the PMBM in (13) is con verted to a PMB, which was detailed in [51] for point targets and in [25] for ETs. V I . D E C E N T R A L I Z E D P O S T E R I O R F U S I O N Here, we present the decentralized approach to rob ust fusion of posterior densities f s ( X | Z s ) computed by independent PMB-ETT filters s with unkno wn prior densities. Fusion can be performed after ev ery update step of the filters or based on a lower rate depending on the application and communication capabilities. A. Robust P osterior Fusion: Principle Robust posterior fusion can be achieved by minimizing the KLA between RFS densities f ( · ) and f s ( · ) for s = 1 , . . . , N sens with respect to f ( · ) . The KLA is defined as [46] ¯ f ω = arg inf f N sens X s =1 ω s D ( f k f s ) , (18) for any combinations of weights ω s ∈ [0 , 1] : P s ω s = 1 . W e hav e introduced D ( f k f s ) as the KLD between RFS densities f ( X | Z 1: N sens ) and f s ( X | Z s ) defined as [46] D ( f k f s ) = Z f ( X | Z 1: N sens ) log f ( X | Z 1: N sens ) f s ( X | Z s ) δ X . (19) The fused posterior (18) is robust in the sense that it is conservati v e and never overconfident w .r .t. the true target uncertainty [39]. Problem (18) was shown to have closed-form solution [46] ¯ f ω ( X | Z 1: N sens ) = Q N sens s =1 f s ( X | Z s ) ω s R Q N sens s 0 =1 f s 0 ( X | Z s 0 ) ω s 0 δ X . (20) Note that (20) is a generalization of the Uhlmann-Julier cov ariance intersection method (c.f. [39]) for posterior RFS densities with unknown priors [40]. The weights ω i can be chosen such that (20) is as peak y as possible [40]. 5 B. Robust PMB P osterior Fusion The posterior density of each PMB-ETT filter is a PMB, therefore, to be able to utilize the fused posterior density as prior for the next time step in each PMB-ETT filter , the fused posterior should be of the PMB form, i.e., ¯ f ω ( X | Z 1: N sens ) = X X u ] X d = X ¯ f u ω ( X u | Z 1: N sens ) × ¯ f d ω ( X d | Z 1: N sens ) (21) comprised of an PPP ¯ f u ω ( X u | Z 1: N sens ) modeling the un- known tar gets and an MB ¯ f d ω ( X d | Z 1: N sens ) modeling detected targets. Note, for brevity we avoid the conditioning on the measurement set in the remainder of this section. A challenge in the fusion step is the fact that sensors do not hav e the same FoVs (c.f. Fig. 1). In the process of fusion, any target that is in the FoV of one sensor and has been detected, but outside the FoV of another sensor, must be treated carefully . In this situation, prior to the fusion, for the first sensor the target corresponds to a detected target represented by a Bernoulli density , whereas for the second sensor the target corresponds to an unknown target, which is represented by the PPP. Because of this, in the closed form KLA solution (20), we must fuse a Bernoulli with a part of the PPP intensity . From this it follo ws that if we solve the KLA (18) separately for the PPP and the MB then we will obtain incorrect results. T o enable a v alid fusion, we propose the approach outlined in the following sub-sections. Consider PMB densities with PPP intensity λ s ( x ) and N s Bernoullis with parameters r i s and f i s ( x ) indexed i = 1 , . . . , N s . Fusion of the PMB densities is simplified if all hav e the same number of Bernoullis, howe ver , in the general case one cannot assume this. Instead we rely on a result from [50, Sec. 4.3.1]: any PPP with intensity λ s ( x ) can be divided into multiple independent PPPs with intensities λ j s ( x ) , where P j λ j s ( x ) = λ s ( x ) . Based on this result, for each sensor we divide λ s ( x ) into M s parts such that N s + M s = K for all sensors s , where parameter K is determined depending on the number of components of the MBs. A robust choice is K = P s N s , so that each Bernoulli component in one sensor can be assigned to any combination of PPP or Bernoullis in the other sensors. The FoV can be taken into account in order to reduce K . 3 If follows from the division into K parts that the PMB densities can be expressed as follows, f s ( X ) = X ] K i =1 X i = X K Y i =1 f i s ( X i ) (22) where f i s ( X i ) is Bernoulli for i ∈ { 1 , . . . , N s } and PPP with intensity λ i s ( x ) for i ∈ { N s + 1 , . . . , N s + M s } . Note that the sum in (22) is implicit and nev er has to be computed; it is sufficient to represent the parameters of the densities f i s ( X i ) . 3 In particular, we partition the deployment region into 2 N sens sub-regions, each determined by a subset of sensors that have each partition in the FoV. For each sub-region l ∈ { 1 , . . . , 2 N sens } , sensor s has N s,l Bernoulli components (with N s = P l N s,l ), we determine K l = max N s,l and K = P l K l . W e are seeking the fused density ¯ f ( X ) that minimizes the KLA (18). Assume that ¯ f ( X ) is of the format (22), i.e., ¯ f ( X ) = X ] K i =1 X i = X K Y i =1 ¯ f i ( X i ) . (23) In [51] it is sho wn that, under this assumption, the minimisa- tion problem (18) can be solved approximately by minimizing an upper bound, N sens X s =1 ω s " K X i =1 D  ¯ f i ( X ) || f π s ( i ) s ( X )  # (24) where π s ∈ Π (1: K ) for all s , and Π (1: K ) is the set of all permutations of the integers { 1 , . . . , K } . For π ∈ Π and i ∈ { 1 , . . . , K } we have π ( i ) ∈ { 1 , . . . , K } and π ( i ) 6 = π ( i 0 ) for i 6 = i 0 . The fusion results will depend on ho w the permutations π s are chosen, which we discuss in Section VI-D. Based on (24), we compute the K components of the fused density ¯ f ( X ) as ¯ f i ( X ) = Q N sens s =1  f π s ( i ) s ( X )  ω s R Q N sens s =1  f π s ( i ) s ( X )  ω s δ X . (25) In summary , our proposed approach to finding the (approx- imately) optimal fused PMB density consists of the following steps: 1) Depending on how the FoVs ov erlap, we select how to divide the PPP intensities such that the PMB densities all hav e K parts, and we find permutations π s such that the fused densities f π s ( i ) s ( X ) have a high degree of similarity in the sense of the KL-div ergence. Note that this entails determining K , λ j s ( x ) , and π s . 2) Fuse the matched PPPs and Bernoullis, see details below in Section VI-C. 3) After the fusion of the K parts of ¯ f ( X ) , we add all the PPP intensities such that a single PPP is obtained. 4) Lastly , we recycle any Bernoullis in ¯ f ( X ) that hav e very low probability of existence. C. Fusion of Bernoullis and PPPs In this subsection, we provide expressions for the fusion of Bernoullis, fusion of PPPs, as well as the fusion of both Bernoullis and PPPs. Lastly , expressions for fusion of Gaussian densities are provided. 1) Fusion of PPPs: Let I be an index set for PPP densities f i ( X ) with intensities λ i ( x ) = µ i f i ( x ) . Fusion of densities f i ( X ) , i ∈ I with fusion weights w i , P i ∈ I w i = 1 , yields a PPP density with intensity [41] ¯ λ ( x ) = Y i ∈ I ( λ i ( x )) w i = ¯ µ ¯ f ( x ) , (26a) ¯ µ = C Y i ∈ I µ w i i , (26b) ¯ f ( x ) = Q i ∈ I ( f i ( x )) w i C , (26c) C = Z Y i ∈ I ( f i ( x )) w i d x. (26d) 6 2) Fusion of Bernoulli RFSs: Consider Bernoulli densities f i ( X ) with probability of existence r i and state PDF f i ( x ) . Fusion of densities f i ( X ) , i ∈ I , | I | ≥ 1 , with fusion weights w i , P i ∈ I w i = 1 , yields a Bernoulli density with parameters ¯ r = C Q i ∈ I r w i i Q i ∈ I (1 − r i ) w i + C Q i ∈ I r w i i , (27a) ¯ f ( x ) = Q i ∈ I ( f i ( x )) w i C , (27b) C = Z Y i ∈ I ( f i ( x )) w i d x. (27c) 3) Fusion of Bernoulli RFSs and PPPs: Let I , | I | ≥ 1 , be an index set for PPP densities f i ( X ) with intensities λ i ( x ) = µ i f i ( x ) , and let J , | J | ≥ 1 , be an index set for Bernoulli densities f j ( X ) with probabilities of existence r j and state densities f j ( x ) . Fusion of PPPs and Bernoullis with fusion weights w i and w j , P i ∈ I w i + P j ∈ J w j = 1 , yields (without approximation) a Bernoulli density with parameters ¯ r = C Q j ∈ J r w j j Q i ∈ I µ w i i Q j ∈ J (1 − r j ) w j + C Q j ∈ J r w j j Q i ∈ I µ w i i , (28a) ¯ f ( x ) = Q i ∈ I ( f i ( x )) w i Q j ∈ J ( f j ( x )) w j C , (28b) C = Z Y i ∈ I ( f i ( x )) w i Y j ∈ J ( f j ( x )) w j d x. (28c) The fusion in (28) holds for any PPP intensities λ i ( x ) , howe v er , given that a Bernoulli RFS represents zero or one object, the fusion results will be more accurate for µ i < 1 . 4) Fusion of Gaussian densities: When the Bernoulli state densities, and/or the PPP intensities are Gaussian, the fused densities ¯ f ( · ) and the normalizing constants C can be com- puted exactly , see, e.g., [46, Eqn. 36] and [52]. Let f i ( x ) = N ( x ; m i , P i ) for i ∈ I . Fusion of the densities, with weights w i , P i ∈ I w i = 1 , yields a Gaussian density , ¯ f ( x ) = Q i ∈ I ( f i ( x )) w i C = N  x ; ¯ m, ¯ P  , (29a) where ¯ m = ¯ P X i ∈ I w i ( P i ) − 1 m i ! , (29b) ¯ P = X i ∈ I w i ( P i ) − 1 ! − 1 , (29c) C = Z Y i ∈ I ( f i ( x )) w i d x (29d) = det(2 π ¯ P ) 1 / 2 Q i ∈ I det(2 π P i ) w i / 2 exp " (29e) 1 2 ¯ m T ¯ P − 1 ¯ m − X i ∈ I w i m T i P − 1 i m i ! # . D. Complexity Reduction From (24), we see that we need to fuse components of the PMBs ov er all permutations π s ∈ Π (1: K ) . The KLD is lowest when the RFS densities are identical. In our application, this means when the single-target posterior densities represent the same tar get. T o reduce computational comple xity tow ards finding the fused posterior density minimizing the KLA, we propose to use the optimal permutation, denoted the best possible fusion map . It is defined and found as follows. Let there be two RFSs densities with index sets I 1 and I s respectiv ely . W e define the best fusion map π ∗ s ∈ Π (1: K ) as the solution of the optimal assignment problem [53] minimize a | I 1 | X n =1 | I s | X m =1 a n,m C n,m (30) sub ject to | I s | X m =1 a n,m = 1 , ∀ n, | I 1 | X n =1 a n,m = 1 , ∀ m, a n,m ∈ { 0 , 1 } , where C n,m denotes the cost for assigning (mapping) com- ponent n in f 1 ( X ) with component m in f s ( X ) . T o solve (30), we can use, e.g., the Munkres algorithm [54], [55]. For N sens > 2 the best fusion map is found by sequentially solving (30) for s = 2 , . . . , N sens . W e define the cost metric in terms of the KLD between the PDFs of the components in f 1 ( X ) and f s ( X ) . For a component n ∈ I 1 with PDF f 1 ,n ( y ) , and similarly for m ∈ I s , the cost is defined C n,m = 1 2 [ D ( f 1 ,n k f s,m ) + D ( f s,m k f 1 ,n )] , (31) which admits a closed-form expression for Gaussian PDFs. Using the procedure above, we find the best fusion map and use it to solve (24). V I I . N U M E R I C A L R E S U L T S W e present first the target extent model used for the simulations. This is followed by the simulation setup, the used performance metrics, and a discussion of the obtained results using the independent PMB-ETT filter with and without posterior fusion. A. Single ET Model The specific target extent model is based on [8], here extended to include the sensor state (including position and orientation) in the measurement model. In the follo wing, we describe the ET state and measurement model, and how this leads to extended Kalman filter (EKF) prediction and update equations which are utilized in the PMB-ETT filter . 7 − 100 − 50 0 50 100 − 50 0 50 100 sensor position sensor orientation sensor FoV Fig. 2. Simulation scenario with two ETs observed by two sensors with overlapping FoVs. The ETs state is plot e very 20 time steps. 1) ET State and Motion Model: The augmented state x k = [ γ k , y T k ] T of a single ET at time k comprises an unknown Poisson rate γ k of number of measurements generated by the ET and the ET spatial state y k . The state has prior 4 PDF f ( x k ) = G ( γ k ; α k , β k ) N ( y k ; ˆ y k , P k ) , (32) where the gamma distribution with parameters α k and β k is a conjugate prior for the rate γ k , and the Gaussian distribution with mean ˆ y k and covariance matrix P k describes the a priori knowledge regarding the ET spatial state y k . The ET spatial state (which includes the ET center, orientation, and extent) and motion model are based on [8] and detailed in Appendix B. Since the measurement rate is independent of the sensor state, decentralized fusion of Sec. VI is only applied to the spatial state y k , while each sensor maintains a local density of the rate γ k (i.e., av erage number of measurements obtained from target with specific sensor) of each target. 2) ET Measur ement Model: In this section, the time index k will be omitted for the sake of brevity . A sensor s located at p s ∈ R 2 with orientation α s observes the ET contour in its local coordinate frame. W e distinguish between three coordinate frames: a quantity in the sensor coor dinate frame is indicated by superscript S , in the ET coor dinate frame by superscript L , and in the global coordinate frame by superscript G . A measurement z S ∈ R 2 is thus z S = h ( y ) + w , (33) where w ∼ N (0 , R ) with measurement noise cov ariance R and h ( y ) = y S c + e ( θ S ) T f ( θ L ) , (34) where θ S = ∠ ( z S − y S c ) , e ( θ S ) is a unit vector in direction θ S , and θ L = θ S − ψ S in which y S c and ψ S are the target 4 Note that in the augmented state PDF (32) the measurement rate γ k and the ET state (including its extent) y k are considered independent, which is a common assumption in ETT (see e.g. [56]). Due to the sensor-to- target geometry (e.g., for a Lidar its angular resolution and the sensor-target distance), the estimated measurement rate of the ET is implicitly sensor dependent, i.e., a different sensor configuration yields a different estimate for γ k . The explicit modeling of the dependence of γ k on the sensor state and y k is not considered here. location and orientation in the sensor frame of reference. 5 Here, f ( θ L ) ≥ 0 is the extent of the target along local angle θ L . W e now express the observation in the global coordinate frame. It can be shown that z S ∼ N ( ˜ h ( y ) , ˜ R ) with ˜ h ( y ) = R ( α s )( y c − p s ) + e T ( θ S ) H f ( θ L ) y f , (35) ˜ R = R + k f e ( θ S ) e T ( θ S ) , (36) where R ( α s ) is a rotation matrix, H f ( θ L ) , and k f are defined in Appendix C. 3) ET Extended Kalman F ilter Equations: Prediction and update equations of an EKF filter using the ET motion and measurement models are given in Appendix D. The predic- tion and update equations for the measurement rate and the predicted likelihood are also stated. All these are utilized in the prediction and update step of the PMB-ETT filter (c.f. Section V). B. Setup If not stated otherwise, there are two ETs present in the scene and their visibility and number of measurements pro- duced per scan depends on the sensor F oV and its configura- tion. W e use a rectangular target of length 5 m and width 3 m to model ETs representing vehicles. Furthermore, we use Lidar type sensors with the follo wing simplified sensor models. Sensor 1 is located at p S 1 = [ − 115 , 120] T , with orientation α S 1 = − 45 ◦ , opening angle of 80 ◦ , angular resolution of 0 . 15 ◦ , and maximum range of 300 m . W e generate a measurement when a ray from the sensor hits an ET and add noise with covariance matrix R S 1 = 0 . 5 I 2 . Sensor 2 is located at p S 2 = [ − 105 , − 80] T , with orientation α S 2 = 45 ◦ , opening angle of 90 ◦ , angular resolution of 0 . 15 ◦ , maximum range of 300 m , and R S 2 = 0 . 02 I 2 . Each sensor produces clutter measurements with rate λ = 2 and uniform spatial distribution c ( z ) = U [ − 200 , 200] 2 . In the simulation, each independent PMB-ETT filter has only access to measurements from one sensor (denoted indep. filter). The fusion filters are independent PMB-ETT filters (denoted fusion filter), but perform posterior fusion according to Sec. VI. W e approximate the PMBM posterior by a PMB which consists of the MB in the MBM that has highest weight. If not stated otherwise, posterior fusion is applied in ev ery time step. W e set the weights ω s = 1 / N sens , resulting in a more conservati v e estimate than achie vable. The fused posterior is then used as the prior for the next filter iteration. For com- parison, an PMB-ETT filter which incorporates measurements from all sensors is used (denoted centralized filter). There, measurements from each sensor are incorporated separately by applying multiple sequential PMB-ETT filter update steps. For all filter v ariants, spatially close measurements are clus- tered into measurement cells using the DBSCAN clustering algorithm [57], where we set the maximum radius for the neighborhood to 4 m and the minimum number of points for a core point to 4 . The simulation scenario is outlined in Fig. 2. The hyper parameters of the GP (see Appendix B) 5 Note that here the unkno wn angle θ S is replaced by a point estimate, which is a simple but inaccurate approach and can be seen as a greedy association model [5]. It is also the approach taken by [8]. 8 − 100 − 80 − 60 − 40 − 40 − 20 0 x y (a) Independent filter 1 − 100 − 80 − 60 − 40 − 40 − 20 0 x y (b) Independent filter 2 − 100 − 80 − 60 − 40 − 40 − 20 0 x y (c) Centralized filter − 100 − 80 − 60 − 40 − 40 − 20 0 x y (d) Fusion filter 1 and 2 Fig. 3. Estimated shape of one ET for different PMB-ETT filters. The true ET center and extent is plotted (green dot, blue solid line), as well as the estimated ones (red dot, red solid line for mean extent, red dotted line for one standard deviation). are l 2 = π / 8 , σ 2 f = 2 , σ 2 r = 2 , and 20 support points are used to track the tar get e xtent. Note that the dimension of the ET state is 26 ( xy -position, orientation, xy -velocity , angular velocity , target extent support points). T arget motion model and its model in the filter is identical with sampling time T = 0 . 5 s , ¯ F k =  1 T 0 1  ⊗ I 3 , (37) ¯ W k = " T 3 3 T 2 2 T 2 2 T # ⊗ diag([0 . 01 , 0 . 01 , 0 . 001]) , (38) β = 0 . 001 in F f (c.f. (51)), and the for getting factor is set to 1 η = 1 1 . 11 . In the filters, the birth intensity has rate λ b = 1 10 and for the spatial distribution we use a single Gaussian centered at location x = [0 , 100] T with cov ariance matrix P = 30 I 2 . The probability of ET surviv al is p S = 0 . 999 , the probability of detection is p D = 0 . 99 . C. P erformance Metrics Multiple target tracking performance is measured by three errors: estimation error for localized targets, number of missed targets, and number of false targets, see, e.g., [58, Sec. 13.6]. The generalized optimal subpattern assignment (GOSP A) met- ric [59], [60] measures all three errors, hence we use it for performance e valuation. 6 Performance of the estimated num- 6 In tracking literature, the OSP A metric [61] is often used. Ho wev er, recent work [62] has shown that the OSP A metric is susceptible to “spooky action at a distance”, which is undesirable. Hence, we do not use the OSP A metric in this paper . ber of targets as well as their center location is thus assessed as follows. Let sets ˆ X = { x 1 , . . . , x n } and Y = { y 1 , . . . , y m } be finite subsets of R N , where without loss of generality n ≤ m . Then [59], [60] d ( c,α,p ) GOSP A = min I n ∈F n ( { 1 ,...,n } ) n X i =1 d ( c ) ( x i , y I n ( i ) ) p + c p α ( m − n )  1 p , (39) where we set the power parameter p = 2 , cut-off distance c = 20 , α = 2 , and d ( c ) ( x , y ) = min( k x − y k 2 , c ) , i.e., the minimum of the Euclidean distance and value c . T o obtain ˆ X , we estimate the detected ETs from the (fused) posterior through comparison of the probability of existence of each Bernoulli component against the threshold r th = 0 . 5 . The set Y contains the true ETs. Note that we only use the mean position of the target center (c.f. Sec. III-B). Performance of the target extent estimation is assessed with the intersection o ver union (IOU) of the true target shape (in the global coordinate frame) and the estimated shape. Let A k be the true ET area in xy -dimension at time step k , and ˆ A k its estimate. Then, the IOU is defined as, see, e.g., [8], IOU( A k , ˆ A k ) = area( A k ∩ ˆ A k ) area( A k ∪ ˆ A k ) . (40) Note that the IOU is, by definition, alw ays between zero for non-ov erlapping target shapes and one when they fully ov er- lap. Thus a well performing ETT filter will yield a high IOU for e very ET. For both GOSP A and IOU average performance results were obtained by a veraging ov er 50 Monte-Carlo runs. D. Discussion of Results Here, we first discuss the ETT filter performance with decentralized posterior fusion in terms of GOSP A and IOU performance metric. This is followed by performing posterior fusion at a lower rate. After that, we inv estigate the case when more than two ETT filters are used for posterior fusion. 1) F ilter P erformance with P osterior Fusion: In Fig. 3, the true and estimated target state is plotted for ET 2 using the different filter variants. W e observe that independent filter 1 (Fig. 3a) estimates the target state with a clear position error (in the positiv e y -direction). This filter uses measurements from sensor 1, where most of the measurements provide information only from the ET’ s top edge due to the horizontal mov ement of the tar get and the sensor pose. In the measure- ment model (33), occlusions caused by the ET itself are not modeled. Due to the simple data association that is used, the receiv ed measurements can then be explained by an ET whose target contour is in proximity of the measurements and the target center is placed north of it. In contrast, independent filter 2 and the centralized filter estimate the target center closer to the true position (Fig. 3b and Fig. 3c). The former filter overestimates the target size in the direction where no measurement is provided, whereas the latter filter utilizes measurements from both sensors and can therefore estimate the target size accurately . The fusion filter utilizes information from both sensors through posterior fusion. 9 0 20 40 60 80 100 10 15 20 time step k measurement rate ET1 (S1) ET1 (S2) ET1 (fusion filter 1) ET1 (indep. filter 1) ET1 (indep. filter 2) ET1 (central. filter) (a) ET 1 0 20 40 60 80 100 15 20 25 30 time step k measurement rate ET2 (S1) ET2 (S2) ET2 (fusion filter 1) ET2 (indep. filter 1) ET2 (indep. filter 2) ET2 (central. filter) (b) ET 2 Fig. 4. The true and estimated measurement rate for each ET and sensor/filter is plotted over time. In Fig. 4, the true and estimated measurement rate γ is plotted over time. W e can observe that the true measurement rate of the ETs is varying o ver time. Although only a simple process model for the measurement rate is used in the filters, they are able to correctly track the rate. This is true for all filter variants except the centralized filter . This filter performs two sequential filter update steps using measurement from dif ferent sensors. Since the number of measurements for an ET are different for each sensor , it follows that after centralized fusion the filter estimates the measurement rate of the ET as the av erage of the two. In Fig. 5, the av erage GOSP A is plotted over time for the simulation scenario illustrated in Fig. 2. The centralized filter has best performance followed by the fusion filter . Independent filter 2 has superior performance compared to the independent filter 1. Note that independent filter 1 is provided with measurements from sensor 1, which has a higher measurement noise compared to sensor 2. In T able I, the a verage IOU is stated for the dif ferent ETs and filter v ariants. W e observe that the IOU for ET 2 is low with independent filter 1 due to the misplaced target center, and with independent filter 2 due to the ov erestimation of the target size. The fusion and the centralized filter show comparable performance. T ABLE I A VE R AG E I O U Filter ET 1 ET 2 Independent 1 0.65 0.25 Independent 2 0.71 0.43 Fusion 0.72 0.75 Centralized 0.70 0.85 2) Low Rate P osterior Fusion: In a real system, it may not be feasible to perform posterior fusion after every filter update step. This can occur when the computers on which the filters run are geographically separated and need to communicate ov er the wireless channel. Therefore, it is worthwhile to in vestigate the filter performance when posterior fusion is performed only ev ery N time steps. In T able II, the average GOSP A is stated for different values of N . W e see that with increasing N the performance of the fusion filters deteriorates, since the information transfer (through fusion) between the filters over time is too low . Fusion filter 1 has w orse performance compared to fusion filter 2 for N ≥ 15 , since it is equipped with the low quality sensor 1. T ABLE II A VE R AG E G O S P A Filter N = 1 N = 15 N = 30 N = 50 Fusion 1 1.12 1.65 2.20 2.76 Fusion 2 1.12 1.80 2.03 2.10 Independent 1 3.20 3.20 3.20 3.20 Independent 2 2.17 2.17 2.17 2.17 Centralized 0.89 0.89 0.89 0.89 3) P osterior Fusion with more filters: In Sec. VI, we pro- posed a procedure to fuse PMB posteriors when there are more than two independent PMB-ETT filters. W e implement this sequentially , where first posteriors from two filters are fused. The outcome is then used to fuse with a not yet fused posterior from one of the remaining filters. This process is repeated until all filter posteriors hav e been incorporated. W e placed four sensors at p S 1 = [ − 150 , − 80] T , p S 2 = [ − 150 , − 50] T , p S 3 = [ − 150 , − 20] T , and p S 4 = [ − 150 , 10] T all with ov erlapping sensor FoVs towards the ETs. The remaining sensor parameters are the same as for sensor 1 used in the previous simulations. In T able III, the average GOSP A, as well as the av erage IOU per ET are stated for a single independent PMB-ETT filter (no fusion), two filters with posterior fusion performed after ev ery filter update, and four filters with posterior fusion. With posterior fusion, performance increases, visible by a decrease of the GOSP A value. Also the a verage IOU increases for all ETs with posterior fusion. Furthermore, fusion of two filter posteriors and four filter posteriors show similar performance in this scenario. T ABLE III M U L T I PL E P O S T ER I O R F U SI O N No. Posteriors Fused GOSP A IOU (ET 1) IOU (ET 2) (no fusion) 1.22 0.61 0.62 2 0.90 0.66 0.67 4 0.86 0.68 0.66 10 0 20 40 60 80 100 0 1 2 3 4 time step k avg. GOSP A Fusion filter 1 Fusion filter 2 Indep. filter 1 Indep. filter 2 Central. filter Fig. 5. The average GOSP A value is plotted o ver time. Posterior fusion is performed in e very time step. V I I I . C O N C L U S I O N S W e proposed a low-complexity decentralized ETT filter that is capable of estimating the presence, state, and shape of ETs accurately . It operates by combining the multibobject densities of PMB form computed by independent ETT filters. Fusion is performed in the minimum KLA sense yielding a fused posterior which is conserv ativ e but ne ver ov erconfident about the estimated states. A lo w-complexity implementation is highlighted, which permits the use of an optimal fusion map- ping between pairs of sensors. The fusion map was identified as the solution of a linear optimal assignment problem based on a cost matrix comprised of the symmetric KLD between target state estimates. In the simulation results, we observed ho w the independent PMB-ETT filters together with the GP target extent model are capable of estimating the state and shape of the present ETs. Furthermore, we observed that fusion of the filters posterior permits a holistic view of the surveillance region spanned by all sensors combined. This resulted in a reduced state estimation error quantified by the GOSP A distance metric, and for the ET shape estimation by an increased area overlap quantified by the IOU value. A P P E N D I X A E T - P M B P A R A M E T E R S Follo wing [3, Sec. IV], we find that w A = Q C ∈ A L C P A ∈A Q C ∈ A L C , (41) L C =          κ C C +  D u + , ` C C  , if C ∩ I = ∅ , | C C | = 1 ,  D u + , ` C C  , if C ∩ I = ∅ , | C C | > 1 , 1 − r i C + + r i C +  f i C + , q D  , if C ∩ I 6 = ∅ , C C = ∅ , r i C +  f i C + , ` C C  , if C ∩ I 6 = ∅ , C C 6 = ∅ . (42) The density f C ( X C ) in (15) is a Bernoulli density with parameters r C =                h D u + ,` C C i κ C C + h D u + ,` C C i , if C ∩ I = ∅ , | C C | = 1 , 1 , if C ∩ I = ∅ , | C C | > 1 , r i C + h f i c + ,q D i 1 − r i C + + r i C + D f i C + ,q D E , if C ∩ I 6 = ∅ , C C = ∅ , 1 , if C ∩ I 6 = ∅ , C C 6 = ∅ , (43) f C ( x ) =              ` C C D u + ( x ) h D u + ,` C C i , if C ∩ I = ∅ , q D ( x ) f i C + ( x ) h f i C ,q D i , if C ∩ I 6 = ∅ , C C = ∅ , ` C C ( x ) f i C + ( x ) D f i C + ,` C C E , if C ∩ I 6 = ∅ , C C 6 = ∅ . (44) A P P E N D I X B E T S P A T I A L S TA T E A N D M O T I O N M O D E L A. Spatial State The ET spatial state including target extent is giv en by y = [( ¯ y ) T , ( y f ) T ] T , (45) where ¯ y = [( y c ) T , ψ , y ∗ T ] T (46) comprising the ET center y c , the ET orientation ψ , and any additional quantities (e.g, velocity) in y ∗ . The variable y f models the target extent, follo wing [8]: Let u denote the local angle w .r .t. the ET orientation and y f i denotes the unknown target extent along input (angle) u i , for a fix ed and finite set of N angles. Then y f = [ y f 1 , . . . , y f N ] T . (47) This vector is modeled as a zero mean GP [63], [64] y f ∼ G P ( 0 , K ( u , u )) (48) with covariance matrix [ K ( u , u )] i,j = k ( u i , u j ) , in which k ( · , · ) is a periodic kernel function. The input of the GP is u = [ u 1 , . . . , u N ] T , and the output is y f . W e utilize the periodic kernel function proposed in [8] k ( u, u 0 ) = σ 2 f exp   − 2 sin 2  | u − u 0 | 2  l 2   + σ 2 r , (49) where σ f , l and σ r are the (known) model hyper-parameters. This function is 2 π periodic, i.e., k ( u + 2 π , u 0 ) = k ( u, u 0 ) , and, thanks to σ r , star con vex object shapes of dif ferent sizes can be described. See [8], [9] for further details and dif ferent choices for the kernel function to describe the extent of an ET with the help of a GP. 11 B. Motion Model Extending Sec. IV -B, the ET follo ws the linear dynamic model ¯ y k +1 = ¯ F k ¯ y k + ¯ w k , (50) where ¯ F k denotes the state transition matrix, and ¯ w k ∼ N (0 , ¯ W k ) with process noise cov ariance ¯ W k . The motion model of the ET contour is [8] y f k +1 = F f k y f k + w f k , (51) where F f = e − β T I N with I N denoting the identity matrix of dimension N , and w f k ∼ N (0 , W f ) with W f = (1 − e − 2 β T ) K ( u f , u f ) . Here, β ≥ 0 denotes the forgetting factor allowing to accommodate targets with extents that change slowly , and T is the sampling time. According to (32), the measurement rate of the ET is assumed independent of the ETs’ spatial state. T o allow the measurement rate to change ov er time an exponential forgetting factor 1 /η is used in the motion model, where the predicted rate is gi ven by the motion of the gamma distrib ution parameters with [3], [56] α k | k − 1 = α k − 1 /η , (52) β k | k − 1 = β k − 1 /η . (53) A P P E N D I X C P RO O F O F ( 4 9 ) – ( 5 0 ) Due to the GP model, we follo w [8] and express f ( θ L ) as f ( θ L ) = H f ( θ L ) y f + e f , (54) H f ( θ L ) = K ( θ L , u f )[ K ( u f , u f )] − 1 , (55) in which e f ∼ N (0 , k f ) with k f = k ( θ L , θ L ) − H f ( θ L ) K ( θ L , u f ) T . (56) Separating signal and noise contribution, and expressing the local states in global frame of reference through y S c = R ( α s )( y c − p s ) (57) ψ S = ψ − α s , (58) where R ( α ) denotes the rotation matrix R ( α ) =  cos( α ) − sin( α ) sin( α ) cos( α )  , (59) we find that with the GP contour model measurement z S in (33) has mean and covariance given by (35)–(36). A P P E N D I X D E T P R E D I C T I O N A N D U P D A T E S T E P S Here, we first describe the EKF prediction and update steps for the ET’ s spatial state. This is follo wed by the update step of the ET measurement rate, and lastly the predicted likelihood utilized in the update step of the PMB-ETT filter for the ET state model of Section IV. A. EKF prediction and update equations W ith the linear ET motion model the standard EKF predic- tion step with initial state y 0 ∼ N ( ˆ y 0 , P 0 ) is [38] ˆ y k | k − 1 = F k ˆ y k − 1 , (60) P k | k − 1 = F k P k − 1 F T k + W k , (61) where F k = blkdiag( ¯ F k , F f k ) , and W k = blkdiag( ¯ W k , W f k ) . W e now extend the EKF update steps deri ved in [8] to in- corporate the (known) sensor state (position p s and orientation α s ). The standard EKF measurement update equations for a detection z k are [38] H k = d d y k ˜ h ( y ) | y = ˆ y k | k − 1 , (62) S k = H k P k | k − 1 H T k + R k , (63) K k = P k | k − 1 H T k S − 1 k , (64) ˆ y k = ˆ y k | k − 1 + K k ( z k − ˜ h ( ˆ y k | k − 1 )) , (65) P k = P k | k − 1 − K k H k P k | k − 1 , (66) where ˜ h ( · ) was defined in (35) and (36). T o linearize the measurement function, we need to compute [8] H k = h d ˜ h ( y ) d y c , d ˜ h ( y ) d ψ , d ˜ h ( y ) d y ∗ , d ˜ h ( y ) d y f i , (67) where in our case d ˜ h ( y ) d y ∗ = 0 . W e get d ˜ h ( y ) d y f = e ( θ S ) H f ( θ L ) , (68) d ˜ h ( y ) d ψ = e ( θ S ) d d ψ H f ( θ L ) y f , (69) d d ψ H f ( θ L ) = − ∂ ∂ u H f ( u )     u = θ L , (70) where [8] d H f ( u ) d u = d d u K ( u, u f )[ K ( u f , u f )] − 1 , (71) d K ( u, u f ) d u = d d u [ k ( u, u f 1 ) , . . . , k ( u, u f N )] , (72) d k ( u, u f i ) d u = − 1 l 2 sin( u − u f i ) k ( u, u f i ) . (73) Further , d ˜ h ( y ) d y c = R ( α s ) + d d u e ( u )     u = y c ( H f ( θ L ) y f ) T + e ( y c ) ×    d d u H f ( u )     u = θ L  T d d w θ L ( w )     w = y c ! T y f   T , (74) where [8] d d u e ( u ) = ( z S − u )( z S − u ) T k z S − u k 3 − 1 k z S − u k I 2 , (75) d d w θ L ( w ) = 1 k z S − y S c k 2 × h z S Y − y S c Y , − ( z S X − y S c X ) i R ( α s ) . (76) 12 Here, the superscript b X and b Y correspond to the first and second dimension of the vector b . Furthermore, we wrote e ( y c ) for e ( θ S ) to indicate the state dependency . T o update the ET spatial state by a set of detections W = { z k,l } n k l =1 , we augment the measurement vector z k = [ z T k, 1 , . . . , z T k,n k ] T , (77) where R k = diag( R k, 1 , . . . , R k,n k ) , (78) ˜ h k ( y k ) = [ ˜ h k, 1 ( y k ) T , . . . , ˜ h k,n k ( y k ) T ] T . (79) B. ET measurement r ate The predicted ET measurement rate γ k | k − 1 has parameters α k | k − 1 and β k | k − 1 (c.f. (32), (53)), which are updated for a set of detections W by [3] α k = α k | k − 1 + | W | , (80) β k = β k | k − 1 + 1 . (81) C. Predicted likelihood for ET -PMB filter The predicted likelihood, used in the update step of the PMB-ETT filter (c.f. Sec. V -2 and Appendix A), for a set of detections W for a single ET is ` W = Γ( α k ) β α k | k − 1 k | k − 1 Γ( α k | k − 1 ) β α k k | W | Y l =1 N ( z k,l − ˜ h ( ˆ y k | k − 1 ,l ) , S k,l ) , (82) where Γ( · ) denotes the gamma function. R E F E R E N C E S [1] Y . Bar -Shalom and X.-R. Li, Multitar get-multisensor tracking: princi- ples and techniques . YBs London, UK:, 1995, v ol. 19. [2] R. P . Mahler , Statistical multisource-multitar get information fusion . Artech House, Inc., 2007. [3] K. Granstr ¨ om, M. Fatemi, and L. Svensson, “Poisson multi- Bernoulli conjugate prior for multiple extended object estimation, ” IEEE T ransactions on Aer ospace and Electronic Systems , DOI 10.1109/T AES.2019.2920220. [4] K. Granstr ¨ om, M. Fatemi, and L. Svensson, “Gamma Gaussian inv erse- Wishart Poisson multi-Bernoulli filter for extended target tracking, ” in 19th International Confer ence on Information Fusion (FUSION) . IEEE, 2016, pp. 893–900. [5] K. Granstr ¨ om, M. Baum, and S. Reuter , “Extended Object Tracking: Introduction, Overview , and Applications, ” JAIF J ournal of Advances in Information Fusion , vol. 12, no. 2, pp. 139–174, Dec. 2017. [6] K. Granstr ¨ om, S. Reuter, D. Meissner , and A. Scheel, “A multiple model PHD approach to tracking of cars under an assumed rectangular shape, ” in International Confer ence on Information Fusion , Salamanca, Spain, Jul. 2014. [7] K. Granstr ¨ om and C. Lundquist, “On the Use of Multiple Measurement Models for Extended T ar get T racking, ” in International Confer ence on Information Fusion , Istanbul, T urke y , Jul. 2013. [8] N. W ahlstr ¨ om and E. ¨ Ozkan, “Extended tar get tracking using Gaussian processes, ” IEEE T ransactions on Signal Pr ocessing , vol. 63, no. 16, pp. 4165–4178, 2015. [9] T . Hirscher, A. Scheel, S. Reuter, and K. Dietmayer, “Multiple extended object tracking using Gaussian processes, ” in Information Fusion (FU- SION), 2016 19th International Conference on . IEEE, 2016, pp. 868– 875. [10] B. W ang, W . Y i, R. Hoseinnezhad, S. Li, L. Kong, and X. Y ang, “Distributed fusion with multi-Bernoulli filter based on generalized co- variance intersection, ” IEEE Tr ansactions on Signal Pr ocessing , vol. 65, no. 1, pp. 242–255, 2017. [11] J. L. Williams, “Marginal multi-bernoulli filters: RFS deriv ation of MHT , JIPD A, and association-based MeMBer, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 51, no. 3, pp. 1664–1687, Jul. 2015. [12] J. W . K och, “Bayesian approach to extended object and cluster tracking using random matrices, ” IEEE T ransactions on Aerospace and Electr onic Systems , vol. 44, no. 3, 2008. [13] M. Baum and U. Hanebeck, “Extended object tracking with random hypersurface models, ” IEEE T ransactions on Aerospace and Electr onic Systems , vol. 50, no. 1, pp. 149–159, Jan. 2013. [14] W . Aftab, A. De Freitas, M. Arvaneh, and L. Mihaylova, “ A Gaussian process approach for extended object tracking with random shapes and for dealing with intractable likelihoods, ” in Digital Signal Pr ocessing (DSP), 2017 22nd International Confer ence on . IEEE, 2017, pp. 1–5. [15] E. ¨ Ozkan, N. W ahlstr ¨ om, and S. J. Godsill, “Rao-Blackwellised particle filter for star-con vex extended target tracking models, ” in 19th Interna- tional Conference on Information Fusion (FUSION) . IEEE, 2016, pp. 1193–1199. [16] R. Mahler, “PHD filters for nonstandard targets, I: Extended targets, ” in Proceedings of the International Conference on Information Fusion , Seattle, W A, USA, Jul. 2009, pp. 915–921. [17] K. Granstr ¨ om, C. Lundquist, and O. Orguner , “Extended target tracking using a Gaussian-mixture PHD filter , ” IEEE T ransactions on Aer ospace and Electronic Systems , vol. 48, no. 4, pp. 3268–3286, 2012. [18] K. Granstr ¨ om and U. Orguner , “A PHD filter for tracking multiple extended targets using random matrices, ” IEEE T ransactions on Signal Pr ocessing , vol. 60, no. 11, pp. 5657–5671, Nov . 2012. [19] A. Swain and D. Clark, “The PHD filter for extended tar get tracking with estimable shape parameters of varying size, ” in Pr oceedings of the International Conference on Information Fusion , Singapore, Jul. 2012. [20] C. Lundquist, K. Granstr ¨ om, and U. Or guner, “ An e xtended target CPHD filter and a gamma Gaussian in verse Wishart implementation, ” IEEE Journal of Selected T opics in Signal Pr ocessing, Special Issue on Multi- tar get T racking , vol. 7, no. 3, pp. 472–483, Jun. 2013. [21] M. Beard, S. Reuter, K. Granstr ¨ om, B.-T . V o, B.-N. V o, and A. Scheel, “Multiple extended target tracking with labelled random finite sets, ” IEEE T ransactions on Signal Pr ocessing , vol. 64, no. 7, pp. 1638–1653, Apr . 2016. [22] A. Scheel, S. Reuter, and K. Dietmayer , “Using separable likelihoods for laser-based vehicle tracking with a labeled multi-bernoulli filter , ” in Information Fusion (FUSION), 2016 19th International Confer ence on . IEEE, 2016, pp. 1200–1207. [23] ——, “V ehicle tracking using extended object methods: An approach for fusing radar and laser , ” in Robotics and Automation (ICRA), 2017 IEEE International Confer ence on . IEEE, 2017, pp. 231–238. [24] M. Michaelis, P . Berthold, D. Meissner , and H.-J. Wuensche, “Heteroge- neous multi-sensor fusion for extended objects in automotiv e scenarios using Gaussian processes and a GMPHD-filter , ” in Sensor Data Fusion: T r ends, Solutions, Applications (SDF), 2017 . IEEE, 2017, pp. 1–6. [25] Y . Xia, K. Granstr ¨ om, L. Sv ensson, and M. F atemi, “Extended target poisson multi-Bernoulli filter , ” arXiv preprint , 2018. [26] J. Williams, “Marginal multi-Bernoulli filters: RFS deriv ation of MHT, JIPDA, and association-based MeMBer, ” IEEE Tr ansactions on Aer ospace and Electr onic Systems , vol. 51, no. 3, pp. 1664–1687, Jul. 2015. [27] Y . Xia, K. Granstr ¨ om, L. Svensson, and A. F . G. Fern ´ andez, “Perfor- mance evaluation of multi-bernoulli conjugate priors for multi-target filtering, ” in Pr oceedings of the International Conference on Information Fusion , Xi’an, China, Jul. 2017. [28] A. F . Garcia-Fernandez, J. W illiams, K. Granstr ¨ om, and L. Svensson, “Poisson multi-Bernoulli mixture filter: direct derivation and imple- mentation, ” IEEE T r ansactions on Aerospace and Electronic Systems , vol. 54, no. 4, Aug. 2018. [29] Y . Xia, K. Granstr ¨ om, L. Svensson, and A. F . G. Fern ´ andez, “ An implementation of the poisson multi-bernoulli mixture trajectory filter via dual decomposition, ” in Proceedings of the International Conference on Information Fusion , Cambridge, UK, Jul. 2018. [30] Y . Xia, K. Granstr ¨ om, L. Svensson, A. F . G. Fern ´ andez, and J. L. W illiams, “Extended target poisson multi-bernoulli mixture trackers based on sets of trajectories, ” in Pr oceedings of the International Confer ence on Information Fusion , Ottaw a, Canada, Jul. 2019. [31] K. Granstr ¨ om, S. Reuter, M. Fatemi, and L. Svensson, “Pedestrian tracking using velodyne data - stochastic optimization for extended object tracking, ” in Pr oceedings of IEEE Intelligent V ehicles Symposium , Redondo Beach, CA, USA, Jun. 2017, pp. 39–46. 13 [32] L. Cament, M. Adams, J. Correa, and C. Perez, “The δ -generalized multi-bernoulli poisson filter in a multi-sensor application, ” in Inter- national Confer ence on Contr ol, Automation and Information Sciences (ICCAIS) , Oct. 2017, pp. 32–37. [33] L. Cament, M. Adams, and J. Correa, “ A multi-sensor , gibbs sampled, implementation of the multi-bernoulli poisson filter , ” in International Confer ence on Information Fusion (FUSION) , Cambridge, UK, Jul. 2018, pp. 2580–2587. [34] K. Granstr ¨ om, L. Svensson, S. Reuter, Y . Xia, and M. Fatemi, “Likelihood-based data association for extended object tracking using sampling methods, ” IEEE T r ansactions on Intelligent V ehicles , vol. 3, no. 1, Mar . 2018. [35] S. Scheidegger , J. Benjaminsson, E. Rosenberg, A. Krishnan, and K. Granstr ¨ om, “Mono-camera 3d multi-object tracking using deep learn- ing detections and pmbm filtering, ” in Pr oceedings of IEEE Intelligent V ehicles Symposium , Changshu, Suzhou, China, Jun. 2018. [36] M. Fatemi, K. Granstr ¨ om, L. Svensson, F . Ruiz, and L. Hammarstrand, “Poisson Multi-Bernoulli Mapping Using Gibbs Sampling, ” IEEE T rans- actions on Signal Processing , vol. 65, no. 11, pp. 2814–2827, Jun. 2017. [37] M. Fr ¨ ohle, C. Lindberg, K. Granstr ¨ om, and H. W ymeersch, “Multisensor poisson multi-bernoulli filter for joint target-sensor state tracking, ” IEEE T ransactions on Intelligent V ehicles , DOI 10.1109/TIV .2019.2938093. [38] D. Simon, Optimal state estimation: Kalman, H infinity , and nonlinear appr oaches . John W iley & Sons, 2006. [39] J. K. Uhlmann, “General data fusion for estimates with unkno wn cross cov ariances, ” in Signal Pr ocessing, Sensor Fusion, and T ar get Recog- nition V , vol. 2755. International Society for Optics and Photonics, 1996, pp. 536–548. [40] R. P . Mahler , “Optimal/robust distributed data fusion: a unified ap- proach, ” in Signal Processing , Sensor Fusion, and T ar get Recognition IX , vol. 4052. International Society for Optics and Photonics, 2000, pp. 128–139. [41] D. Clark, S. Julier, R. Mahler , and B. Ristic, “Robust multi-object sensor fusion with unknown correlations, ” in Sensor Signal Processing for Defence (SSPD) . IET , 2010. [42] M. B. Guldogan, “Consensus Bernoulli filter for distributed detection and tracking using multi-static Doppler shifts. ” IEEE Signal Pr ocess. Lett. , vol. 21, no. 6, pp. 672–676, 2014. [43] M. ¨ Uney , S. Julier , D. Clark, and B. Ristic, “Monte Carlo realisation of a distributed multi-object fusion algorithm, ” in Sensor Signal Processing for Defence (SSPD) . IET , 2010. [44] M. ¨ Uney , D. E. Clark, and S. J. Julier, “Distributed fusion of PHD filters via exponential mixture densities, ” IEEE Journal of Selected T opics in Signal Processing , vol. 7, no. 3, pp. 521–531, 2013. [45] T . Li, J. M. Corchado, and S. Sun, “On generalized cov ariance intersec- tion for distributed PHD filtering and a simple but better alternati ve, ” in 20th International Confer ence on Information Fusion (Fusion) . IEEE, 2017, pp. 1–8. [46] G. Battistelli, L. Chisci, C. Fantacci, A. Farina, and A. Graziano, “Consensus CPHD filter for distributed multitarget tracking, ” IEEE Journal of Selected T opics in Signal Processing , vol. 7, no. 3, pp. 508– 520, 2013. [47] G. Battistelli, L. Chisci, C. Fantacci, A. Farina, and B.-N. V o, “ A verage Kullback-Leibler di vergence for random finite sets, ” in 18th Interna- tional Confer ence on Information Fusion (Fusion) . IEEE, 2015, pp. 1359–1366. [48] C. Fantacci, B.-N. V o, B.-T . V o, G. Battistelli, and L. Chisci, “Consensus labeled random finite set filtering for distributed multi-object tracking, ” arXiv preprint arXiv:1501.01579 , 2015. [49] S. Li, W . Y i, R. Hoseinnezhad, G. Battistelli, B. W ang, and L. K ong, “Robust distributed fusion with labeled random finite sets, ” IEEE T ransactions on Signal Pr ocessing , vol. 66, no. 2, pp. 278–293, 2017. [50] R. P . Mahler, Advances in Statistical Multisource-Multitar get Informa- tion Fusion . Artech House, 2014. [51] J. L. Williams, “ An efficient, variational approximation of the best fitting multi-Bernoulli filter, ” IEEE T ransactions on Signal Processing , vol. 63, no. 1, pp. 258–273, 2015. [52] P . Bromiley , “Products and conv olutions of Gaussian probability density functions, ” T ina-V ision Memo , vol. 3, no. 4, p. 1, 2003. [53] Y . Bar-Shalom, “Multitarget-multisensor tracking: advanced applica- tions, ” Artech House, Norwood, MA , 1990. [54] J. Munkres, “ Algorithms for the assignment and transportation prob- lems, ” J ournal of the society for industrial and applied mathematics , vol. 5, no. 1, pp. 32–38, 1957. [55] F . Bourgeois and J.-C. Lassalle, “ An extension of the Munkres algorithm for the assignment problem to rectangular matrices, ” Communications of the A CM , vol. 14, no. 12, pp. 802–804, 1971. [56] K. Granstr ¨ om and U. Orguner , “Estimation and maintenance of measure- ment rates for multiple extended target tracking, ” in Information Fusion (FUSION), 2012 15th International Conference on . IEEE, 2012, pp. 2170–2176. [57] M. Ester, H.-P . Kriegel, J. Sander, X. Xu et al. , “ A density-based algorithm for discovering clusters in large spatial databases with noise. ” in Kdd , vol. 96, no. 34, 1996, pp. 226–231. [58] S. Blackman and R. Popoli, Design and Analysis of Modern T rac king Systems . Norwood, MA, USA: Artech House, 1999. [59] A. S. Rahmathullah, ´ A. F . Garc ´ ıa-Fern ´ andez, and L. Svensson, “ A metric on the space of finite sets of trajectories for e valuation of multi-target tracking algorithms, ” arXiv pr eprint arXiv:1605.01177 , 2016. [60] ——, “Generalized optimal sub-pattern assignment metric, ” in Informa- tion Fusion (Fusion), 2017 20th International Conference on . IEEE, 2017, pp. 1–8. [61] D. Schuhmacher , B.-T . V o, and B.-N. V o, “ A consistent metric for performance ev aluation of multi-object filters, ” IEEE T ransactions on Signal Processing , vol. 56, no. 8, pp. 3447–3457, Aug. 2008. [62] A. Garcia-Fernandez and L. Svensson, “Spooky effect in optimal ospa estimation and ho w gospa solves it, ” in International Conference on Information Fusion (FUSION) , Ottawa, Canada, Jul. 2019. [63] C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning . MIT Press, 2006. [64] C. M. Bishop, P attern Recognition and Machine Learning . Springer , 2006. Markus Fr ¨ ohle (S’11) received the B.Sc. and M.Sc. degrees in T elematics from Graz Uni versity of T echnology , Graz, Austria, in 2009 and 2012, respectiv ely . He obtained the Ph.D. de gree in Signals and Systems from Chalmers University of T ech- nology , Gothenbur g, Sweden, in 2018. From 2012 to 2013, he was with the Signal Processing and Speech Communication Laboratory , Graz Univ ersity of T echnology . From 2013 to 2018, he was with the Department of Electrical Engineering, Chalmers Univ ersity of T echnology . He joined Zenuity AB in 2019. His current research interests include localization and tracking. Karl Granstr ¨ om (M’08) is a postdoctoral research fellow at the Department of Signals and Systems, Chalmers University of T echnology , Gothenbur g, Sweden. He receiv ed the MSc degree in Applied Physics and Electrical Engineering in May 2008, and the PhD degree in Automatic Control in November 2012, both from Link ¨ oping University , Sweden. He previously held postdoctoral positions at the De- partment of Electrical and Computer Engineering at Univ ersity of Connecticut, USA, from September 2014 to August 2015, and at the Department of Electrical Engineering of Link ¨ oping University from December 2012 to August 2014. His research interests include estimation theory , multiple model estimation, sensor fusion and target tracking, especially for extended tar gets. He receiv ed paper aw ards at the Fusion 2011 and Fusion 2012 conferences. He has organised several workshops and tutorials on the topic Multiple Extended T arget T racking and Sensor Fusion. At Fusion 2018 the International Society for Information Fusion (ISIF) awarded him the ISIF Y oung In vestigator A ward for his contributions to extended target tracking research and his services to the research community . 14 Henk W ymeersch (S’01, M’05, SM’19) obtained the Ph.D. degree in Electrical Engineering/Applied Sciences in 2005 from Ghent Univ ersity , Belgium. He is currently a Professor of Communication Sys- tems with the Department of Electrical Engineer- ing at Chalmers University of T echnology , Sweden. Prior to joining Chalmers, he was a postdoctoral researcher from 2005 until 2009 with the Laboratory for Information and Decision Systems at the Mas- sachusetts Institute of T echnology . Prof. W ymeersch served as Associate Editor for IEEE Communication Letters (2009-2013), IEEE Transactions on W ireless Communications (since 2013), and IEEE Transactions on Communications (2016-2018). His current research interests include cooperative systems and intelligent transportation.