Louis Guerlin
email: louis.guerlin@renault.com

Benjamin Pannetier

Maxime Derome

Michèle Rombaut

SPIE Proceedings: Study on group target tracking to counter swarms of drones

Keywords: Bayesian inference, Computer simulations, Detection and tracking algorithms, Electronic filtering, Kinematics, Radar, Sensors, Statistical analysis, Stochastic processes, Target acquisition, Target detection, Time metrology

This article focuses on Group Target Tracking (GTT) to counter swarms of drones using Random Finite Sets (RFSs) and Random Matrix (RM) approaches. Tracking swarms of drones is analog to tracking extended targets that are characterized by their continuously evolving shape and composition. Extended target tracking for groups of targets finds various applications in the literature because detecting and tracking each individual target of a group is computationally demanding and unnecessary if the group itself can be modeled. Elliptic shapes offers a suitable representation for most groups, and their inference is quite inexpensive with the random matrix approach. Indeed, they are efficient when coupled with random finite sets based filters, which represents the current state of the art for Bayesian multi-target tracking. In this work, a practical implementation of a labeled Poisson/Multi-Bernoulli filter using random matrices for group target tracking is proposed. This study compares several random matrix prediction and update algorithms with and without random finite sets based filters on several dataset.

INTRODUCTION

Group Target Tracking (GTT) for swarms of drones is a special case of extended target tracking where the extended target corresponds to a non-rigid group of targets whose kinematic parameters are similar. [START_REF] Granstrom | Extended object tracking: Introduction, overview and applications[END_REF] Such a group sees its composition and shape evolve through time which adds some complexity, but tracking an entire group instead of each individual is much more efficient both in terms of computational complexity and localization errors. [START_REF] Koch | Bayesian approach to extended object and cluster tracking using random matrices[END_REF] Group target tracking is well documented in the scientific community and most of the existing algorithms rely on the combination of a Multi Target Tracking (MTT) algorithm with a non rigid-extended target model. [START_REF] Granstrom | Extended object tracking: Introduction, overview and applications[END_REF][START_REF] Lan | Tracking of extended object or target group using random matrix -Part I: New model and approach[END_REF][START_REF] Feldmann | Tracking of extended objects and group targets using random matrices[END_REF][START_REF] Granström | New prediction for extended targets with random matrices[END_REF] The field of Multi Target Tracking (MTT) consists in the detection and kinematic estimation of objects of interest in an area monitored by at least one sensor. In addition, the track continuity of detected targets must be enforced throughout the whole process despite the presence of false alarms and missed detection, highlighting the need for a robust and efficient solution. Multiple Hypothesis Tracking (MHT), Joint Probability Density Association Filter (JPDAF) [START_REF] Bar-Shalom | Multitarget-multisensor tracking: principles and techniques[END_REF] and Random Finite Sets (RFS) based filters [START_REF] Mahler | Statistical multisource-multitarget information fusion[END_REF] are amongst the most popular MTT approaches.

Random Finite Sets (RFSs) based trackers has drawn a lot of attention from the scientific community during last decades because they offer a mathematical formalism from which optimal MTT algorithms can be derived. With RFSs, the MTT problem is redefined into a set theorical one: the goal is to study sets whose cardinal are finite but subject to random effects. As any random variable, the estimation of the cardinal of a RFS can be done with the help of Bayesian formalism, provided that a proper description of the set probability density function exists. The point process theory offers such concepts and Finite Set STatistics (FISSTs) are a sufficient simplification of this field for MTT: [START_REF] Mahler | Statistics 101 for multisensor, multitarget data fusion[END_REF][START_REF] Mahler | Statistics 102 for multisource-multitarget detection and tracking[END_REF] under certain assumptions, the filters derived with RFSs and FISSTs are optimal, as they minimize Cramér Rao lower bound. [START_REF] Mahler | Multitarget bayes filtering via first-order multitarget moments[END_REF] The first successful implementations of RFS filters were based on Poisson Point Processes (PPP) to model the multi-target system. [START_REF] Mahler | Multitarget bayes filtering via first-order multitarget moments[END_REF] A Poisson Point Process (PPP) is characterized by its first statistical moment: the Probability Hypothesis Density (PHD), also called intensity function. This statistical moment allows the estimation of the cardinal of the multi-target system since its integral over the observed space gives the expected number of targets present in the system. Moreover, since it is a statistical moment with a known model, it can be estimated recursively using Bayes and Chapman-Kolmogoroff equations. [START_REF] Mahler | Statistical multisource-multitarget information fusion[END_REF] What is even more remarkable is the fact that the PHD implicitly embeds the kinematic information about the targets, hence they are simultaneously estimated during the prediction and update procedure. [START_REF] Vo | The gaussian mixture probability hypothesis density filter[END_REF] Latest filters involve Multi Bernoulli (MB) RFSs to model tracked targets. The cardinal of a single Bernoulli RFS is either zero or one, hence it characterizes a single target with a probability of existence and a state probability density function. Using one Bernoulli density per tracked target leads to a Multi Bernoulli (MB) density. Once again, estimating the cardinal of a MB density using Bayes and Chapman-Kolmogoroff equations, [START_REF] Mahler | Statistical multisource-multitarget information fusion[END_REF] leads to the implicit estimation of the state of the targets. Indeed, each individual Bernoulli density embeds a state probability density function. Several implementations of MB filters can be found in the literature including the delta-Generalized Labeled Multi Bernoulli filter, [START_REF] Vo | Labeled random finite sets and multi-object conjugate priors[END_REF] the Labeled Multi Bernoulli filter, [START_REF] Reuter | The labeled multi-bernoulli filter[END_REF] or the Generalized Labeled Multi Bernoulli filter: [START_REF] Vo | An efficient implementation of the generalized labeled multibernoulli filter[END_REF] these are labeled implementations, meaning each target has a track history. However, filters relying solely on MB RFSs are often heavier compared to their Poisson based counterparts: the Poisson Multi Bernoulli Mixture [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF] (PMBM) filter offers a good trade-off between efficiency and versatility, thus it will be considered in this article.

The Poisson Multi Bernoulli Mixture (PMBM) filter brings solutions to most of the issues encountered in MTT for countering drones, but, like most RFS based filters implementations, it only considers point targets. Hence, tracking a group of drones only with a PMBM filter may lead to augmented processing time, coalescence or large localization errors. A simple but effective solution consists in reducing the group of drones to a single extended target. An elegant and efficient answer to group target tracking is the Random Matrix approach 1, 2, 4 (RM). The main idea of this approach is to consider the shape of the group as an ellipse, hence it can be represented by a positive semidefinite matrix subject to modeling and measurement noises. Other methods exist for extended target tracking, but RM versatility, ease of use and Bayesian integration make it a highly regarded solution in RFS based extended target tracking. [START_REF] Lan | Tracking of extended object or target group using random matrix -Part I: New model and approach[END_REF][START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF] This article focuses on Group Target Tracking (GTT) to counter swarms of drones using RFS and RM approaches. The original contributions of this work are the adaptation of existing PMB algorithm to the GTT problem and the comparison between several RM prediction and update algorithms [START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF] for GTT using both simulated radar dataset and real radar dataset. To present these achievements, the first parts of this paper introduce briefly the RFS and RM theory. With this background, a group target tracking filter is presented in the third section: the extended Labeled PMB. Its results are presented in the last section of this work.

FROM RANDOM FINITE SETS TO THE POISSON MULTI BERNOULLI FILTER

Random Finite Sets (RFS) and Finite Set STatistics (FISST) might be challenging to apprehend. These methods consider point targets. The first subsection brings some background about Bayesian Multiple Target Tracking with these tools: more progressive and in depth overview are available in the literature. 7-9

Bayesian inference and Random Finite Sets (RFSs)

The Kalman filter aims at estimating a vector variate with Bayes and Chapman-Kolmogorov equations. It leads to the optimal unbiased filter to estimate a single target state vector. The goal of Random Finite Sets (RFS) and FInite Sets STatistics (FISSTs) is to bridge the gap between the single vector Kalman filter and a set theoric equivalent of this Kalman filter for the Multi Target Problem (MTT). Indeed, the specific MTT problem is more complex because of the extra association problem, the presence of false alarms, missed detections, and the unknown number of targets. To model these phenomenon, the MTT problem can be redefined into a set theoric one, where the studied random variables are Random Finite Sets (RFSs), and they are estimated with Bayes and Chapman-Kolmogorov equations. For instance, these sets can represent the false alarm set, or the target set, that are both subject to random effects. An RFS is a countable set whose cardinal is a random variable. Hence the goal is to estimate the true cardinal of the RFS with the help of Bayesian inference, using a suitable probability density functions and a set decomposition of the problem. It is remarkable that the association problem and state estimation of the targets are implicitly solved when estimating the cardinal a RFS with a Bayes filter.

As a random variable, the cardinal of a RFS can be modeled by any probability law depending on the application. Hence, rewriting the Bayes recursion with RFSs is straightforward: using the notation of Table 1, the Chapman-Kolmogorov equation becomes

p k|k-1 (X k |Z {1:k-1} ) = f k|k-1 (X k |X k-1 ) • p k-1|k-1 (X k-1 |Z {1:k-1} )δX k-1 , (1)
and for the Bayes equation the result is

p k|k (X k |Z {1:k} ) = l k (Z k |X k ) • p k|k-1 (X k |Z {1:k-1} ) l k (Z k |X k ) • p k|k-1 (X k |Z {1:k-1} )δX k . (2)
These equations involve unordered sets instead of vectors. This raise several questions, such as how to compute the integral of a set probability density function or how to perform any other regular calculus with a set theoric approach. The answers to these questions can be found in the Finite Set STatistics (FISSTs) mathematical toolbox [START_REF] Mahler | Statistics 101 for multisensor, multitarget data fusion[END_REF][START_REF] Mahler | Statistics 102 for multisource-multitarget detection and tracking[END_REF] that allows an easy manipulation of RFSs and equivalent representations of their probability density laws to derive optimal MTT filters. [START_REF] Mahler | Statistical multisource-multitarget information fusion[END_REF] Next subsection, the most documented RFSs probability distribution functions are introduced. They are essential in the construction of the Poisson Multi Bernoulli filter.

Common Random Finite Sets (RFSs) densities

In order to derive the Poisson Multi Bernoulli (PMB) filter, two RFSs density probability laws need to be introduced here: the Poisson and the Bernoulli probability densities. a) Poisson Probability Density Function: A Poisson RFS, is assumed when the elements of a set X = {x 1 , ..., x n } are independently and identically distributed and the cardinal of the set follows a Poisson probability law with expected value N . The associated probability function is given by

p p (X = {x 1 , ..., x n }) = e -N n i=1 I(x i ) , (3)
where I(x) is the intensity or Probability Hypothesis Density (PHD) of this RFS. The PHD is the first statistical moment of a Poisson RFS, it is also called intensity because it can be regarded as the number of target per unit of space. Hence, its integral is the expected number of elements in X:

I(x)dx = N . (4)
Since the intensity function is the first statistical moment of a Poisson RFS and it is linked to the cardinal of the RFS, its estimation can be done with Bayesian recursion and it lead to the estimation of the cardinal of the RFS.

Two common intensity functions are the constant intensity function or the Gaussian Mixture 11 (GM) intensity function; while a constant intensity is useful for homogeneously distributed processes, a Gaussian Mixture intensity is considered when the studied process is not homogeneous anymore. For instance, if a sensor covering 3 units of space is subject to independently and homogeneously distributed false alarms with a rate of I F A = 2 false alarm per unit of space, then I F Adx = 6 false alarms are expected to occur at each time step.

Hence, a Poisson RFS offers good results and versatility. However, using only the intensity of a Poisson RFS to describe the MTT problem doesn't question the existence of the elements of a set. b) Bernoulli Probability Density Function: The Bernoulli RFS is an interesting alternative to the Poisson RFS, in that it represents a singleton set, and emphasizes on the existence probability of its element. Hence the probability density function of a Bernoulli RFS X is given by

p b (X) =      1 -r X = {∅} rg(x) X = {x} 0 |X| ≥ 2 , ( 5
)
where r is the existence probability and g(x) is the state probability density of the corresponding element in case of existence. Since the set contains at most one element, the probability associated to the existence of two or more elements is null.

This parameterization is more refined than the previous one, and estimating a Bernoulli RFS boils down to estimating its cardinal. However, one Bernoulli RFS is necessary for each singleton set, it more computationally demanding than a Poisson RFS: in this paper, this model is applied to detected and tracked targets only.

Since the MTT problem often involves more than one target, if all tracked targets are subject to Bernoulli probability density functions then the set of these targets is subject to the union of their Bernoulli probability densities: a Multi Bernoulli density. Thus, the probability density function for a set X = {x 1 , ...x n } that is the union of n Bernoulli RFSs can be expressed as

p mb ({x 1 , ..., x m }) = m i=1 1 -r i •   1≤i1 =... =im≤n m j=1 r ij 1 -r ij g ij (x j )   , (6)
where m ≤ n. This equation mirror the probability of {x 1 , ..., x m } to exist in the set {x 1 , ...x n }. Since a set is unordered, the sum is taken over all the permutation of m elements of the Multi Bernoulli RFS X. As in the probability density function 5, r ij and g ij stand for the existence and state probability density of each target.

The goal of the PMB filter is to combine the lightness of Poisson RFSs with the extra information offered by Bernoulli RFSs, improving the precision and performance of the filter.

The Poisson Multi Bernoulli (PMB) filter

The Poisson Multi Bernoulli (PMB) filter construction involves several RFSs and assumptions for the prediction and update phases. [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF] In this subsection, only the most important results are presented, the demonstration can be found in other publications. [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF][START_REF] García-Fernández | Poisson multi-bernoulli mixture filter: direct derivation and implementation[END_REF] a) Prediction: Regarding the prior probability derivation, three RFSs are assumed to coexist:

• X β , the RFS of new targets subject to a Poisson process of intensity I β ,

• X U k-1 , the RFS of unobserved and supposedly existing targets also subject to a Poisson process of intensity

I U k-1 ,
• X S k-1 , the Multi Bernoulli RFS of detected and surviving targets at time k -1.

For X U k-1 and X S k-1 the surviving probability of each target when in state x is noted p s (x). Moreover, the evolution process is supposed to be Markovian for both of these RFSs, and the associated transition density is denoted f k|k-1 (x). Also, it should be noted that X β being the RFS of predicted new targets, it is directly integrated in the unobserved targets RFS during the prediction. With the probability density functions described in previous section, the aim of the prediction step is to estimate this probability density function:

p k|k-1 (X k ) = X U X S =X k p p (X U )p mb (X S ) , (7)
where X U includes X β . With equations and 3 and 6, this can be rewritten as

p k|k-1 (X k ) = X U X S =X k e -I U k|k-1 (x)dx x∈X U I U k|k-1 (x)   |X S | i=1 1 -r i •   1≤i1 =... =im≤|X S |   |X S | j=1 r ij k|k-1 1 -r ij k|k-1 g ij k|k-1 (x j )       , (8)
hence only a few parameters need to be estimated in the Bayesian recursion:

I U k|k-1 , (r i k|k-1 , g i k|k-1 ) i∈{1...n k|k-1 } . ( 9
)
Using Chapman-Kolmogorov equation 1 with the evolution model f k|k-1 (x), the survival probability p s (x), the birth intensity I β and the posterior density parameters

I U k-1|k-1 , (r i k-1|k-1 , g i k-1|k-1 ) i∈{1...n k-1|k-1 }
leads to the predicted parameters

I U k|k-1 (x k ) = I β (x k ) + p s f k|k-1 (x k |x k-1 ) • I U k-1|k-1 (x k-1 ) dx k-1 g i k|k-1 (x k ) = f k|k-1 (x k |x k-1 ) • g i k-1|k-1 (x k-1 ) dx k-1 ∀i ∈ 1, n k-1|k-1 r i k|k-1 = p s r i k-1|k-1 ∀i ∈ 1, n k-1|k-1 n k|k-1 = n k-1|k-1 , (10)
where n k-1|k-1 denotes the cardinal of the MB mixture. In this paper, the probability of survival is assumed to be a constant, p s (x) = p s . b) Correction: When new measurements are available the prediction can be corrected with new information using Bayes equation. The full proof for the following equations is not presented here due to a space limitation: other papers cover this subject thoroughly with FInite Sets STatistics (FISSTs), [START_REF] Mahler | Statistical multisource-multitarget information fusion[END_REF][START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF] and some assumptions regarding the measurement model:

• measurements might belong to the RFS of false alarms. It is supposed to be an homogeneous Poisson Point Process with intensity I F A : integrated over the observation area it gives the expected number of false alarms.

• a new or already tracked target can generate at most one measurement: in particular, it can be missed.

• the likelihood of the measurements is noted

l k (z j k |x k ), where Z k = {z j k } j∈ 1,m k is the set of the m k measurements received at time k.
These assumptions are considered valid in this article too.

Since a tracked or yet unobserved target can generate at most one measurement, when no measurement is available the correction step must be carried on accordingly. Hence, for the RFS X U of unobserved and supposedly existing targets, the update is straightforward, only the missed detection hypothesis is considered here, otherwise they wouldn't be unobserved:

I U k|k (x k ) = (1 -p d )I U k|k-1 (x k ) , (11)
that is lowering the predicted intensity of the undetected targets by a factor (1 -p d ) where p d is the constant probability density.

Tracked targets might be missed too, leading to the validation of the prediction for the state detection probability such as

g i,m k +1 k|k (x k ) = g i k|k-1 (x k ) r i,m k +1 k|k = (1 -p d ) • r i k|k-1 w i,m k +1 w i,m k +1 = 1 -r i k|k-1 + (1 -p d ) • r i k|k-1 , (12)
where the superscript m k + 1 is bound to the missed detection hypothesis and i ∈ 1, n k|k-1 . w i,m k +1 is the likelihood of associating one measurement with a known tracked target distribution.

Second, in case of detection the correction step consists in updating the prior knowledge with the information brought by a new measurement z k : it can be seen as finding its best possible association with a tracked target or a new target. In this context, the scenario that leads to the the highest marginal probability for the global association is considered the best. [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF] As stated earlier, a measurement can emerge from either a tracked target, or a new target, or it can be a false alarm. When a tracked target is detected, its existence probability is set to one and its state probability density can be derived from Bayes equation 2 using FISSTs 15 such as:

g i,j k|k (x k ) = p d • l k (z j k |x k ) • g i,j k|k-1 (x k ) w i,j r i,j k|k = 1 w i,j = p d l k (z j k |x k ) • g i,j k|k-1 (x k )dx k , (13)
for i ∈ 1, n k|k-1 and j an index over the m k measurements at time k. w i,j is the probability of association between measurement j and target i without considering any other association.

In addition to the tracked target, for each measurement a new Bernoulli hypothesis is created, based on the set of unobserved targets:

g i ,j k|k (x k ) = p d • l k (z j k |x k ) • I U k|k-1 (x k ) w i ,j r i ,j k|k = w i ,j I F A + w i ,j w i ,j = p d l k (z j k |x k ) • I U k|k-1 (x k ) dx k , (14)
with i ∈ n k|k-1 , n k|k-1 + m k , which gives m k hypothetical new targets whose existence probability takes into account the false alarm intensity of the sensor I F A and the likelihood of measurement j being generated by I U k|k-1 , the unobserved targets intensity.

After the update of the Multi Bernoulli density with equations 12, 13 and, 14, the resulting density is not a MB anymore, it is a MB Mixture. Indeed, all Bernoulli components of the corrected density can't exist at the same time, but they may be combined from several manners to form new MB densities as shown in figure 1.

In order to obtain a MB density after the update, only the best assignment between targets and measurements should be retained here. To find this best assignment, the marginal probability of association must be computed, Multi Bernoulli RFS with Loopy Belief Propagation (LBP) for instance. [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF][START_REF] Williams | Approximate evaluation of marginal association probabilities with belief propagation[END_REF] Computing the marginal association probabilities can be seen as normalization.

x 1 k|k-1 x 2 k|k-1 Poisson RFS z 2 k z 1 k z 2 k z 1 k n.d. n.d. x∈X U k|k-1 I U (x) z 2 k z 1 k w 1,2 w 2,1 w 2,2 w 1,1 w 1,3 w 2,3 w 3,2 w 3,1 Association Matrix w 1,2 w 2,1 w 2,2 w 1,1 w 1,3 w 2,3 w 3,2 w 3,1 0 0 Meas. z 1 k Target x 1 k Target x 2 k New target x 3 k New target x 4 k Meas. z 2 k n.
The resulting normalized association weights ŵi,j can be compared to keep only the best assignment between the tracked targets, the new targets, and the measurements: this can be done using Munkres algorithm. [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] In the end of the process, only one Multi Bernoulli (MB) density is retrieved from the correction step. Hence the Poisson Multi Bernoulli (PMB) filter is closed under Bayes equation: the prior density is the same as the posterior density. This is a desirable characteristic for an iterative filter.

The Poisson Multi Bernoulli filter is very effective for point targets with a regular Gaussian state probability density. This article is interested in tracking groups of point targets, thus these groups are not a point anymore. Hence, next section discusses the state probability density, measurement likelihood, and evolution model that allow to track groups of targets.

GROUP TARGETS AND THE RANDOM MATRIX THEORY

The main problem of RFS filters for MTT is the point targets assumption regarding the extent of the tracked objects. To avoid wrong associations, coalescence or phantom movement of closely spaced targets evolving together, extended target tracking can be considered for group targets. Indeed, a group of targets can be seen as an extended target whose shape is non rigid.

Reminders about the Random Matrix theory

Amongst extended target tracking methods, the Random Matrix theory is a relatively inexpensive stochastic approach, 1, 2 this motivates its application to model groups of targets. a) Correction: In the Random Matrix approach, the correction step relies on Bayes equation, hence the measurement likelihood, and the prior density probability of the extended target must be modeled.

When the size of a target is not negligible compared with the resolution of the sensor, it generates multiple measurements. It is called an extended target and its extent X should be estimated, with some assumptions regarding its probability distribution. The shape of an extended target can be estimated jointly with its kinematic vector x using Bayesian inference. The resulting joint density is

p(x k , X k |Z {1:k} ) = p(x k |X k , Z {1:k} )p(X k |Z {1:k} ) , (15)
where Z {1:k} is composed of k consecutive sensor swipes, and the swipe Z t at time t ∈ 1, k consists in m measurements:

Z t = {z j t } j∈ 1,m .
It is important to emphasize on the decomposition offered by equation 15: the first density represent the regular state vector while the second is the random matrix representation of the extent, a d × d positive semidefinite covariance matrix. So as to model properly these probability densities, they can be deduced from the measurement likelihood, assuming a Gaussian measurement noise. [START_REF] Koch | Bayesian approach to extended object and cluster tracking using random matrices[END_REF] In a first approximation, only one extended target is tracked and it returns m measurement at each time step: it can be a swarm of drones for instance. The mean and the covariance of the associated point cloud can be computed as

z k = 1 m m j=1 z j k Z k = 1 m m j=1 (z j k -z k )(z j k -z k ) T (16)
Since each measurement is subject to a Gaussian noise, the measurements likelihood can be noted 2

l k (Z k |m, x k , X k ) = m j=1 N (z j k ; (H k ⊗ I d )x k , X k ) ∝ N (z k ; (H k ⊗ I d )x k , X k m )W(Z k ; m -1, X k ) (17)
where (H k ⊗ I d ) is the observation matrix, and W(•) denotes the Wishart density associated to the extent. A Wishart density is a d dimensionnal χ 2 law, hence it models the uncertainty of the empirical covariance computed from Gaussian measurements. Regarding (H k ⊗ I d ), ⊗ denotes the Kronecker product, I d is the identity matrix and H k corresponds to a 1 × 3 linear observation vector when the acceleration, the speed, and the position are studied. In this paper, the sensor returns only the position:

H k ⊗ I d = 1 0 0 ⊗ I d (18)
The prior probability density of the extended target, p k|k-1 (x k , X k |Z k-1 ), is assumed to be separable in an extent matrix density and a state vector density as in equation 15. Also, in order for the filter to be closed under Bayes equation, it is necessary to use the densities that are conjugate to the Gaussian Wishart measurement likelihood. This leads to the Gaussian Inverse Wishart (GIW) prior:

p k|k-1 (x k , X k |Z k-1 ) = N (x k ; µ k|k-1 , P k|k-1 ⊗ V k|k )IW(X k ; ν k|k-1 , V k|k-1 ) (19)
where {µ k|k-1 , P k|k-1 } are the predicted parameter for the state vector probability density, and {ν k|k-1 , V k|k-1 } are the predicted parameters for the extent probability density. It is important to note that the Gaussian prior of x k is based on X k : this simplifying assumption helps with the computation and it holds since, for now, the shape is assumed to evolve slowly between two sensor acquisitions.

Using the likelihood 17 and the prior probability density 19 leads to the posterior probability density 2

p k|k (x k , X k |Z k ) = N (x k ; µ k|k , P k|k ⊗ V k|k )IW(X k ; ν k|k , V k|k ) ( 20
)
with the updated parameters from Table 2.

The Gaussian Inverse Wishart extended target model is closed under Bayes equation. As discussed in section 2.3, it is a desirable characteristic and it also needs to be the case after the prediction step. b) Prediction: For the prediction step, two simplifying assumptions are necessary: 2 the evolution process is markovian and the extent is independent from the state vector. Thus, the markovian evolution density can be written as

f k+1|k (x k+1 , X k+1 |x k , X k , Z {1:k} ) = f k+1|k (x k+1 |x k , X k+1 )f k+1|k (X k+1 |X k ) (21)
where only the predicted state vector x k+1 depends on the predicted extent X k+1 . This last observation raises an issue when using this density within the Chapman-Kolmogorov equation. Indeed in

p k+1|k (x k+1 , X k+1 |Z {1:k} ) = f k+1|k (x k+1 |x k , X k+1 )f k+1|k (X k+1 |X k ) p k|k (x k |X k , Z k )p k|k (X k |Z k )dx k dX k (22)
it is desirable to have two distinct integrals for simplification purposes. Hence, it can be assumed that the object extension evolution has no impact on the state vector during the prediction:

2 f k+1|k (x k+1 |X k+1 , x k ) ≈ f k+1|k (x k+1 |X k , x k ).
With this approximation, the single integral of the prediction step becomes separable

p k+1|k (x k+1 , X k+1 |Z {1:k} ) = f k+1|k (x k+1 |X k , x k )p k|k (x k |X k , Z k )dx k f k+1|k (X k+1 |X k )p k|k (X k |Z k )dX k (23)
The first integral corresponds to the predicted state vector and the second to the predicted extent.

Regarding the Markov transition density for the state vector integral, Van Keuk's model can be assumed. [START_REF] Koch | Multiple hypothesis track maintenance with possibly unresolved measurements[END_REF] It is a linear model with a Gaussian noise on the acceleration, its evolution matrix is given by

F k+1|k ⊗ I d =   1 ∆ t 1 2 ∆ 2 t 0 1 ∆ t 0 0 e -∆ t θ   ⊗ I d ( 24
)
where δ t is the acquisition time interval, and θ is a maneuver correlation time. With this model the dynamic covariance noise matrix is expressed as

Q k+1|k ⊗ X k =   0 0 0 0 0 0 0 0 a 2 (1 -e -2 ∆ t θ )   ⊗ X k (25)
a is the acceleration noise root mean square. This dynamic noise covariance matrix is assumed to be dependent on the current extent, which is justified because the bigger the group, the harder its dynamic evolution is to describe. For instance, a large group of targets will likely break into smaller groups because of the individual behaviors of its targets. [START_REF] Koch | Bayesian approach to extended object and cluster tracking using random matrices[END_REF] Thereby, with the help of equations 23 and 20, the state vector prediction process is similar to a standard Kalman prediction

p k+1|k (x k+1 |X k |Z {1:k} ) = N (x k+1 ; (F k+1|k ⊗ I d )x k , Q k+1|k ⊗ X k )N (x k ; µ k|k , P k|k ⊗ X k )dx k = N (x k+1 ; µ k+1|k , P k+1|k ⊗ X k ) , (26)
where Table 2 present the prediction equations.

As for the extent density prediction, an heuristic approach can be consider to induce a formulae closed under Chapman-Kolmogorov. Even if a more formal proof exists, the heuristic approach is simpler and based on a few hypothesis: 2

• the extent's density should be an Inverse Wishart p k+1|k (X k+1 |X k ) = IW(X k+1 ; ν k+1|k , V k+1|k ), in order to get a closed prediction equation.

• the extent's shape doesn't change much between two acquisitions, its expectation remains the same:

E[X k+1 ] = E[X k ], where E[X k ] = V k|k /(ν k|k -2d -2).
• Since ν k|k can be assimilated to the precision of the estimated extent, and because the precision of the estimation should decrease with the prediction step, a temporal decay τ is introduced. [START_REF] Koch | Bayesian approach to extended object and cluster tracking using random matrices[END_REF] With these assumptions, the parameters of the predicted extent probability density can be approximated as in Table 2.

Finally, the full random matrix model is closed under Bayes and Chapman-Kolmogorov equation with the Gaussian Inverse Wishart probability density as posterior and prior probability densities. However, numerous assumptions were made in order to achieve this goal, some of which are questionable for group target tracking with a radar: possible improvements are discussed in next section.

Modelling Groups as Ellipses

A group of targets is analog to an extended target in some aspects: the group size is superior to the resolution cell of the sensor so it returns more than one measurement per acquisition, and the kinematic states of its targets are similar. However a group has its own specificities with respect to a rigid extended target: its shape might evolve through time, as well as its composition and its density. That is why, some constraints must be defined in order to avoid the tracking of large and incoherent groups. For instance, using a limit size for a group avoid large uncertainty on the state vector covariance. Taking the specificities of groups into account leads to the refinement of the RM model. a) Time evolving number of measurements: The first hypothesis that can be questioned is the constant amount of measurement m returned by the extended target at each time step:

l k (Z k |m, x k , X k ) can be rewritten l k (Z k |γ k , x k , X k )
, where γ k is the expected number of measurements for a group. [START_REF] Granström | New prediction for extended targets with random matrices[END_REF] Counting the average number of occurrences γ k of an event happening during a time interval ∆ t is the typical application of a Poisson probability density. This Poisson density is supposed independent of the shape of the group target. [START_REF] Gilholm | Poisson models for extended target and group tracking[END_REF][START_REF] Gilholm | Spatial distribution model for tracking extended objects[END_REF] Hence the measurement likelihood 17 can be augmented

l k (Z k |γ k , x k , X k ) ∝ l k (m k |γ k )l k (z k |x k , X k )l k (Z k |X k ) = P(m k ; γ k )N (z k ; (H k ⊗ I d )x k , X k m )W(Z k ; m k -1, X k ) (27)
Since the conjugate prior of this Poisson distribution is a Gamma distribution, the prior probability density function can be rewritten:

p k|k (x k , X k |Z k ) = G(γ k ; α k|k , κ k|k )N (x k ; µ k|k , P k|k ⊗ V k|k )IW(X k ; ν k|k , V k|k ) = GGIW(γ k , x k , X k ; α k|k , κ k|k , µ k|k , P k|k , ν k|k , V k|k ) ( 28
)
where α k|k is the expected number of measurement arriving at an expected rate κ k|k , such as E(γ k ) = α k|k /κ k|k . This proposed model fits the group target structure better, thus it will be retained here. This model is a Gamma Gaussian Inverse Wishart (GGIW) density. [START_REF] Lundquist | An extended target CPHD filter and a gamma gaussian inverse wishart implementation[END_REF] Its prediction and update procedures for the Gaussian and Inverse Wishart parameters don't differ from the equations in Table 2, but two more parameters needs to be estimated throughout the Bayesian recursion: α and κ.

α and κ update equations are available in Table 2. For the prediction step, exponential forgetting is retained here for α k+1|k and κ k+1|k . This exponential forgetting is tuned such as γ k+1 evolve over a time window w e = η η-1 where η can be seen as a forgetting factor. [START_REF] Lundquist | An extended target CPHD filter and a gamma gaussian inverse wishart implementation[END_REF] Hence the prediction step for the gamma density maintains the expected number of measurements, E(γ k ) = E(γ k+1 ). b) Measurement model: Here, the measurement likelihood is further refined. Indeed, equation 17 implicitly assumes that the extent of a group target is predominant with respect to the measurement error. When performing group target tracking with a radar, this hypothesis may be questioned. Several methods exist to solve this problem. [START_REF] Granstrom | Extended object tracking: Introduction, overview and applications[END_REF] In this paper the approximation holds in a matrix B k that describes the distortion of the measured extension

X k : 3 l k (Z k |γ k , x k , X k ) = m j=1 N (z j k ; (H k ⊗ I d )x k , B k X k B T k ) ( 29
)
This formulation is flexible enough to model an approximation of the additive Gaussian noise, indeed when

B k = (cE[X k ] + R k ) 1 2 E[X k ] -1 2 (30)
then the resulting estimated error covariance is

B k X k B T k = (cE[X k ] + R k ) 1/2 E[X k ] -1/2 X k E[X k ] -T /2 (cE[X k ] + R k ) T /2 ≈ cX k + R k (31)
where c is a scalar weight, and R k is the covariance matrix of the sensor measurement error. This new model resembles the additive Gaussian noise model, [START_REF] Feldmann | Tracking of extended objects and group targets using random matrices[END_REF] without the inconvenience of an assumed independence between the state vector and the extent densities. The value of c has already been discussed, [START_REF] Vivone | Converted measurements random matrix approach to extended target tracking using x-band marine radar data[END_REF] it will be set to 1 4 for group targets in this paper. However, when tracking drones, some targets might be under-resolved or point targets. This leads to a poor extent estimation. We propose to set c to the value c = 1 4 (1 -e 1-α k|k β k|k ), so that it varies between 0 and 1 4 as a function of the expected number of measurements α k|k β k|k , ignoring the extent when necessary.

Table 2. Prediction and update cycles tested in this paper. Time indexes for F , Q, A, υ, K, N , and S are omitted.

Original algorithm [START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF] Algorithm we propose for swarms of drones

Prediction α k|k-1 = α k-1|k-1 /η α k|k-1 = α k-1|k-1 /η β k|k-1 = β k-1|k-1 /η β k|k-1 = β k-1|k-1 /η µ k|k-1 = (F ⊗ I d )µ k-1|k-1 µ k|k-1 = (F ⊗ I d )µ k-1|k-1 P k|k-1 = F P k-1|k-1 F T + Q P k|k-1 = F P k-1|k-1 F T + Q ν k|k-1 = 2d + 2 + e -∆t/θ Υ ν k|k-1 = 2υ(Υ-1)(Υ-2) Υ(Υ+υ) + 2d + 4 V k|k-1 = e -∆ t /θ (ν k|k-1 -2d-2) Υ V k-1|k-1 V k|k-1 = υ(ν k|k-1 -2d-2) Υ AV k-1|k-1 A T Υ = (ν k-1|k-1 -2d -2) Υ = (ν k-1|k-1 -2d -2) Update α k|k = α k|k-1 + m k α k|k = α k|k-1 + m k β k|k = β k|k-1 + 1 β k|k = β k|k-1 + 1 µ k|k = µ k|k-1 + Kε µ k|k = µ k|k-1 + Kε P k|k = P k|k-1 -KHP k|k-1 P k|k = P k|k-1 -KSK T ν k|k = ν k|k-1 + m k ν k|k = ν k|k-1 + m k V k|k = V k|k-1 + N + Z k V k|k-1 = V k|k-1 + N + B -1 Z k B -T ε = z k -(H ⊗ I d )x k|k-1 ε = z k -(H ⊗ I d )x k|k-1 X = V k|k-1 ν k|k-1 -2d-2 X = V k|k-1 ν k|k-1 -2d-2 S = HP k|k-1 H T + X m k S = HP k|k-1 H T + |B| d/2 m k K = P k|k-1 HS -1 K = P k|k-1 HS -1 N = X 1/2 S -1/2 εε T S -T /2 X T /2 N = S -1 εε T B = (cX + R) 1/2 X -1/2
c) Evolution model for the shape: The last improvement presented in this paper is the evolution process of a group target. A simple transition might be that the group should expand at each time step, which can be interpreted as a loss of precision due to the prediction step. The model associated with this transition probability density is a Wishart density, since it is conjugate to the Inverse Wishart density from the update. Thus, the usage of Chapman-Kolmogorov equation gives:

3 p k|k-1 (X k |Z {1:k-1} ) = f k|k-1 (X k |X k-1 )p k-1|k-1 (X k-1 |Z {1:k-1} )dX k-1 = W(X k ; υ k , A k X k-1 V T k-1|k-1 )IW(X k-1 ; ν k-1|k-1 , V k-1|k-1 )dX k-1 (32)
where υ k is the degree of freedom, which can be related to the precision of the prediction, and A k is the evolution matrix. In this equation, the integration step results in a Generalized Beta Type II density, 2, 3 which can be approximated via moment matching to an Inverse Wishart density IW(X k ; ν k|k-1 , V k|k-1 ).The equation for these predicted parameters are presented in Table 2. The demonstration of these results is out of the scope of this paper, more details are available in the aforementioned papers. [START_REF] Granstrom | Extended object tracking: Introduction, overview and applications[END_REF][START_REF] Koch | Bayesian approach to extended object and cluster tracking using random matrices[END_REF][START_REF] Lan | Tracking of extended object or target group using random matrix -Part I: New model and approach[END_REF] The evolution model adopted here is an expansion model, hence A k can be expressed as

A k = ρ 0 0 ρ ( 33
)
where ρ is an expansion coefficient. It is set to 0.7 in this paper.

In this paper, we propose to gather these improvements in order to track swarms of drones. Next section discusses the integration of the RM model inside the Poisson Multi Bernoulli filter.

GROUP TARGET POISSON MULTI BERNOULLI FILTER (GT-PMB)

In section 2.3, the Poisson Multi Bernoulli (PMB) filter was introduced. This section aims at bridging the gap between the point target PMB filter and the Random Matrix (RM) approach for group target tracking proposed in previous section. More details about the rigid extended target PMBM filter are available in the literature. [START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF]

Parameterization of the Group Poisson Multi Bernoulli filter

Section 2.3 introduced I U , {(r i , g i )} i∈ 1,n as sufficient parameters for the description of the multi target system with a Poisson Multi Bernoulli (PMB) filter. In order to track groups of targets, section 3.2 introduced an appropriate state probability density: the Gamma Gaussian Inverse Wishart 1 (GGIW) probability density.

The GGIW density is closed under Chapman-Kolmogorov and Bayes equations. Moreover, it is well suited for the Multi Target Tracking (MTT) problem with non rigid groups of targets as presented in previous section. Thus, interfacing this density with the PMB filter is straightforward at least for the tracked targets:

g i (x k , X k , γ k ) = GGIW(x k , X k , γ k ; Ξ i ) ∀i ∈ 1, n (34)
where Ξ = {α, κ, µ, P, ν, V} is a meta-parameter that gather the GGIW density parameters.

For unobserved but supposedly existing targets, the RFS X U of section 2.3, the parameter to estimate is an intensity function I U . The integral of I U over the observed space should be the expected number of unobserved targets. Hence, I U is defined as a mixture of GGIW densities:

I U (x k , X k , γ k ) = n U i=1 λ U,i GGIW(x k , X k , γ k ; Ξ U,i ) ( 35
)
where λ U,i represents the weight associated to the component i of the mixture. If all the weights λ U,i are equal to one, then the expected number of unobserved target is n U . It should be noted that I β , the intensity of the new targets, is also described by a GGIW mixture, more on this subject can be read in next section.

Using GGIW densities to describe undetected and tracked group targets means that the prediction and update equations presented in Table 2 are still relevant here.

Derivation of the Group Target Poisson Mutli Bernoulli Filter

In this section, equation 10 to 14 are interfaced with the GGIW state probability density to obtain the Group Target Multi Bernoulli prediction and correction equations. a) Prediction: Regarding the prediction step, Multi Bernoulli (MB) and Poisson densities can be dealt separately, leading to a direct derivation for each RFS. For the MB density, the prediction translates to the application of the Markov transition density presented in sections 3.1 and 3.2:

g i k|k-1 (x k , X k , γ k ) = GGIW(x k , X k , γ k ; Ξ i k|k-1 ) ∀i ∈ 1, n k|k-1 r i k|k-1 = p s r i k-1|k-1 ∀i ∈ 1, n k|k-1 (36)
for each tracked target i, with the predicted parameters Ξ i k|k-1 from Table 2. The existence probability is computed according to equation 10 with the constant survival probability p s . For the two Poisson RFSs X U and X β , the prediction step corresponds to their union. For instance, as shown in the prediction step of figure 2 at t = 1, the union of X U and X β corresponds to the concatenation of I β with I U t=0|0 . In general case, at time k, the birth intensity and the corrected unknown intensity are handled as explained below.

t = 0 t = 0 t = 1 t =
The birth intensity I β of new targets, represents the places where new targets are the most likely to appear. As shown in figure 2, the birth intensity I β remains the same at time t = 0 and t = 1. To put it into equation, it is a mixture of n β GGIW components with birth weights λ β,i and parameters Ξ β,i , as shown in equation 35:

I β (x k , X k , γ k ) = n β i=1 λ β,i • GGIW(x k , X k , γ k ; Ξ β,i ) (37)
I U k|k is the updated intensity of unobserved targets: it is subject to the Markov transition model. Since I U k|k is a mixture of GGIW densities, applying the prediction equation 10 leads to:

I U pred (x k , X k , γ k ) = p s f k|k-1 (x k , X k , γ k |x k-1 , X k-1 , γ k-1 ) • I U k-1|k-1 (x k-1 , X k-1 , γ k-1 )dx k-1 dX k-1 dγ k-1 = n U k-1|k-1 i=1 (p s • λ U,i k-1|k-1 ) f k|k-1 (x k , X k , γ k |x k-1 , X k-1 , γ k-1 ) • GGIW(x k-1 , X k-1 , γ k-1 ; Ξ U,i k-1|k-1 )dx k-1 dX k-1 dγ k-1 (38)
hence it is decomposed into the prediction of the weight of each mode, (p s • λ U,i k-1|k-1 ), and the state prediction of each mode of the GGIW density which is carried on as shown in Table 2.

Finally, the full predicted intensity for the RFSs subject to a Poisson density is:

I U k|k-1 (x k , X k , γ k ) = I β (x k , X k , γ k ) + p s f k|k-1 (x k , X k , γ k |x k-1 , X k-1 , γ k-1 ) • I U k-1|k-1 (x k-1 , X k-1 , γ k-1 )dx k-1 dX k-1 dγ k-1 = I β (x k , X k , γ k ) + I U pred (x k , X k , γ k ) (39)
Therefore, after the prediction steps,

I U k|k-1 (•) is composed of n U k|k-1 = n U k-1|k-1 + n β components.
As depicted in figure 35, when initializing the algorithm at t = 0, n U 0|-1 = n β . To summarize , an unobserved mixture of group targets is created from a mixture of new targets. Its weights λ U,i are initialized as λ β,i then they decrease iteration after iteration until they drop bellow a threshold and are destroyed.

During its lifespan, the mixture of unknown group targets can generates new tracked targets if measurements exist, as explained in next section. b) Correction: The correction step is straightforward for the Poisson RFS X U in case of non detection. The predicted mixture is validated and the weights w U,i k|k are adjusted according to the non detection probability (1 -p d )

I U k|k (x k , X k , γ k ) = (1 -p d )I U k|k-1 (x k , X k , γ k ) = n U k|k-1 i=1 (1 -p d )w U,i k|k-1 GGIW(x k , X k , γ k ; Ξ U,i k|k-1 ) = n U k|k i=1 w U,i k|k GGIW(x k , X k , γ k ; Ξ U,i k|k ) (40)
Regarding the MB RFS of detected and tracked targets, several hypotheses need to be considered here: the missed, the detected, and the new target hypotheses. Using the RM update to create these hypotheses implies that the measurements should be organized into clusters.

As discussed in section 3.2, a size limit should be imposed on the groups of targets in order to avoid large error covariance matrices. Hence the clustering of measurements can rely on a parametric algorithm such as the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm. [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] The resulting clusters have a mean z j k and a covariance Z j k computed with equations 16, based on m j k measurement where j ∈ 1, M k is an index on the set of M k clusters at time k. Other approaches exist to cluster the measurement space, a good example is the multiple partitioning method. [START_REF] Granstrom | Extended target tracking using a gaussian-mixture phd filter[END_REF] Once the measurements are clustered, the update of existing and new targets can be carried on. The update of already tracked targets is straightforward for the missed detection hypothesis. It corresponds to the validation of the predicted state density, and the update of the existence probability,

g i,M k +1 k|k (x k , X k , γ k ) = GGIW(x k , X k , γ k ; Ξ i k|k-1 ) i.e. Ξ i,M k +1 k|k = Ξ i k|k-1 r i,M k +1 k|k = (1 -p d )r i k|k-1 w i,M k +1 = 1 -r i k|k-1 + r i k|k-1 (1 -p d ) (41)
where w i,M k +1 is the weight of the hypothesis, ∀i ∈ 1, n k|k-1 , with the indexing strategy from section 2.3.

Then for each measurement j ∈ 1, M k and each tracked target i ∈ 1, n k|k-1 , the updated parameters Ξ i,j k|k are computed according to the update strategy chosen in Table 2. This lead to the correction equations:

g i,j k|k (x k , X k , γ k ) = GGIW(x k , X k , γ k ; Ξ i,j k|k ) r i,j k|k = 1 w i,j = r i,j k|k-1 p d P(m j k ; α i k|k-1 /κ i k|k-1 ) × N (z j k ; (H k ⊗ I d )x i k|k-1 , S k ) × W(Z j k ; ν k|k-1 , V k|k-1 ν k|k-1 -2d -2 ) ( 42
)
where P(•) denotes the Poisson probability density function.

To get a full update, the creation of new Bernoulli hypotheses must be included. Based on each cluster j ∈ 1, M k , a new Bernoulli component is created with the help of the intensity of unknown target I U :

g i+M k ,j k|k (x k , X k , γ k ) = 1 w i+M k ,j n U k|k-1 i =1 w U,i k|k-1 GGIW(x k , X k , γ k ; Ξ U,i ,j k|k ) r i+M k ,j k|k = w i+M k ,j /(w i+M k ,j + I F A ) w i+M k ,j = r i,j k|k-1 p d n U k|k-1 i =1 P(m j k ; α U,i k|k-1 /κ U,i k|k-1 ) × N (z j k ; (H k ⊗ I d )x U,i k|k-1 , S k ) × W(Z j k ; ν k|k-1 , V k|k-1 ν k|k-1 -2d -2 ) (43)
Ξ U,i ,j k|k denotes the updated GGIW component resulting from the association of the mode i of I U with measurement z j k , and I F A stands for the false alarm intensity of the sensor. As mentioned in 2.3, at this stage, the output of the PMB algorithm is a mixture of Multi Bernoulli densities with unnormalized weights w i,j k|k . Before selecting the best assignment, the marginal probability for each association is computed with a Loopy Belief Propagation algorithm 15 (LBP). These marginalized probabilities can then be compared using Munkres algorithm, [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] hence the resulting posterior is a PMB.

Update Prediction

Input:

{Ξ i k|k-1 , r i k|k-1 } n k|k-1 i=1 ; Input: {Ξ i k|k , r i k|k } n k|k i=1 ; {Ξ U,i k|k , w U,i k|k } n U k|k i=1 {Ξ U,i k|k-1 , w U,i k|k-1 } n U k|k-1 i=1 ; {z j k , Z j k } M k j=1
Output:

{Ξ i k+1|k , r i k+1|k } n k+1|k i=1 ; Output: {Ξ i k|k , r i k|k } n k|k i=1 ; {Ξ U,i k|k , w U,i k|k } n U k|k-1 i=1 {Ξ U,i k+1|k , w U,i k+1|k } n U k+1|k i=1 For i ∈ 1, n k|k-1 : For i ∈ 1, n k|k : Compute w i,M k +1 k|k , r i,M k +1 k|k and Ξ i,M k +1 k|k
as in eq. 41

Compute r i k+1|k and Ξ i k+1|k as in eq. 36 For j ∈ 1, M k :

For i ∈ 1, n U k|k : Compute w i,j k|k , r i,j k|k and Ξ i,j k|k as in eq. 42 Compute r U,i k+1|k and Ξ U,i k+1|k as in eq. 38 For j ∈ 1, M k|k :

For i ∈ n U k|k , n U k|k + n β : i = n k|k + j Compute r U,i k+1|k and Ξ U,i k+1|k as in eq. 37 Compute w i ,j k|k , r i ,j k|k and Ξ i ,j k|k as in eq. 43 n U k+1|k = n U k|k + n β Compute the marginal association from w i,j k|k 15 Find the n k|k best assignments with Munkres algorithm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] For i ∈ 1, n U k|k-1 : Compute w U,i k|k and Ξ U,i k|k as in eq. 40

n U k|k = n U k|k-1
Table 3 summarizes the prediction and update procedures presented in this section. Since only one hypothesis is propagated after the update step, the filter is named Group Target Poisson Multi Bernoulli.

Other algorithms involved in Group Targets Tracking

The Group Target Poisson Multi Bernoulli filter (GT-PMB) relies on several other algorithms in order to simplify the computations. This section introduces them. a) Pruning procedure: Regarding the pruning procedure there is a threshold to set. Any component of I U that has a weight λ U,i bellow this threshold is discarded. We propose to set its value such as any component of the GGIW mixture I U shall survive during n surv iteration. Since each component is initialized with a constant weight λ β , the threshold is computed as

T hresh = λ β p nsurv-1 s (1 -p d ) nsurv . ( 44
)
This threshold remains valid for tracked targets. In this paper the n surv parameter is set depending on the false alarm intensity of the sensor: for the simulated dataset presented below, n surv = 2 provides good results.

b) Recycling unassigned measurements: Moreover, in addition to the traditional way of initializing new targets, successive unassigned measurements are used to create new components for the GGIW mixture I β . Recycling unassigned measurements allow to retrieve approximate speeds that are closer to the reality: the filter converges faster. However this comes with a shortcoming when high false alarm rates are reached: the filter might initializes many false targets. c) Delaying tracked targets: Since the recycling of unassigned measurements might lead to the validation of numerous false alarms, each new target is subject to a delay of n del iterations before its final validation. This leads to a slight lag between the first detection and the display of the target. In practice, with the simulated dataset 4.a, n del = 3 iterations has shown to be an appropriate delay before the validation of a tracked target. d) Data pre-processing: Finally, the measurements are subject to two pre-processing steps: first a polar to Cartesian conversion, then a clustering with DBSCAN algorithm. [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] Moreover, the resulting clusters are gated with tracked targets, using Kirubarajan's gating procedure. [START_REF] Kirubarajan | Ground target tracking with variable structure imm estimator[END_REF] This leads to a computationally less demanding algorithm. Figure 3 is a block diagram that includes all the algorithms involved in the GT-PMB algorithm. The results given by this framework are presented and discussed in next section on a simulated and a real dataset.

Update

RESULTS

Baseline and parameters

The performances of the Group Target Poisson Multi Bernoulli for swarms of drones is evaluated on two datasets: one simulated dataset and one real dataset. In order to compare the results, the Poisson Multi Bernoulli (PMB) filter [START_REF] Williams | Marginal multi-bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer[END_REF][START_REF] García-Fernández | Poisson multi-bernoulli mixture filter: direct derivation and implementation[END_REF] without Random Matrices (RM), and the Extended Target Poisson Multi Bernoulli (ET-PMB) filter with a Random Matrix (RM) model for extended rigid targets [START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF] are the baseline here. They will be compared to the Random Matrix implementation we proposed in sections 3 and 4 of this paper that takes into account the non rigid aspect of the targets: the Group Target Poisson Multi Bernoulli (GT-PMB).

An important remark should be made regarding the ET-PMB filter: [START_REF] Granstrom | Poisson multi-bernoulli mixture conjugate prior for multiple extended target filtering[END_REF] because of the fast evolution of the shape of a group target, the weight proposed in the original paper for the association hypothesis is dropped here: while testing, it lead to very high track discontinuities. It is replaced by the association weight w found in equations 42 and 43, which is less dependent on the shape and more on the cardinal.

The advantages of the Random Matrix approach in comparison with the traditional PMB filter for extended or group targets will be studied in this section. Moreover, the improvements that arise from taking the nonrigidity of groups into account will be discussed. To compare the results, three metrics are retained here: the Generalized Optimal Sub-Pattern Assignement (GOSPA) metric [START_REF] Rahmathullah | Generalized optimal sub-pattern assignment metric[END_REF] (with parameters c = 10, p = 1 and α = 1), the cardinal of the estimated set of targets, and the average Root Mean Squared Error (RMSE).

For the Random Matrix approach, it is important to define the critical size of a group, since a clustering of the sensor output is done before filtering the data. As discussed in previous section, this clustering is done with a DBSCAN algorithm. Hence two parameters are required: the minimum number of elements in a group n min , and the minimum distance between them d min . Here, a group is defined by at least two targets with similar kinematic parameters, separated by a distance of 15 meters for the first dataset and 3 meters for the second dataset.

The first dataset is a simulation where one group and one single target are merging then separating, as shown in figure 4.a. As for the second one, it is the acquisition of a group of person that separate and merge forming a T-shaped path, as can be seen in 4.b.

Discussion with a simulated dataset

The objective of this first simulated dataset is to see how the algorithms react to a nonlinear maneuver with merging, splitting, and an uneven number of targets per group. An added difficulty comes from the fact that the top target is alone at the beginning and the end of the scenario, as shown in figure 4.a. Indeed, the Random Matrix model assumes that several measurements are available for each extended or group target at each time step, which is not always the case. Regarding the group target of this scenario, it is slightly unresolved, since it is the union of only 5 targets.

As for the simulated sensor, the characteristics of the radar used for the real acquisition are reported here. The measurements are 2D polar points converted to Cartesian coordinates, the detection probability is set to p d = 0.9, the range standard deviation varies between σ r = {1, 4} meter, and the angular standard deviation ranges between σ θ = {0.1, 0.6} degrees. The false alarm rate is set to N F A = 40 false alarms per radar swipe. The sampling period of the sensor is ∆ t = 0.5 seconds.

The graphs from figure 5 are the averaged results over 100 Monte Carlo runs. The results of the PMB filter are only given in the case where σ r = 1 meter and σ θ = 0.1 degree, for better readability regarding other results. The simple Gaussian state probability density of the PMB filter leads to an highly underestimated number of targets, as shown in graph 5.b. It should be noted that the PMB filter is tracking the true cardinal of the problem: it should be around 6 targets. Indeed, this approach struggles at properly initializing the targets composing the group and when it does, it fails to maintain them because of wrong association, and phantom movements. Another remark can be made about the RMSE that is lower for the PMB filter. Allowing a higher The impact of the Random Matrix model reduces phantom movements and track discontinuities which can be clearly seen in the GOSPA metric of both the ET-PMB and the GT-PMB filters as shown in graph 5.a. With σ r = 1 and σ θ = 0.1, the GT-PMB and the ET-PMB filter seem to perform similarly. However, the higher the noise, the more the ET-PMB filter struggles to estimate properly the number of targets, as shown in graphs 5.e. Indeed, for the GT-PMB, considering that the sensor noise is non negligible with respect to the shape of the group allows to track single targets more accurately, this lead to better results at the beginning and the end of the scenario, when one group and one single target are tracked.

Moreover, the additive Gaussian noise model of the GT-PMB proposed in section 3.2 shows better tracking performances with the group target, as seen in the middle of the scenario when only one group is tracked on graph 5.e. This improved tracking performance comes with a higher RMSE with respect to the ET-PMB filter: in fact, the ET-PMB filter converges quickly, but since the noise is high, with σ r = 4 and σ θ = 0.6, it fails to maintain the track continuity. Hence, what happen is a lower RMSE for the ET-PMB filter, at the cost of a worst tracking continuity, meaning higher cardinality error, as seen in graphs 5.e and 5.f. In conclusion, the Group Target Poisson Multi Bernoulli filter works well with under-resolved single and group targets, which is very important when tracking swarms of drones: when a drone separates from the swarms, its extent becomes negligible with respect to the sensor noise. Moreover, compared to the ET-PMB filter, it proved to be more robust to higher sensor noise.

Discussion with a real dataset

For this dataset, an initial group of 9 peoples walks away from the radar sensor. Then they separate at a ninety degree angle, group again, and walk back toward the radar sensor. The resulting tracks form a T-shape, as can be seen in the figure 4.b. Most of the characteristics of the radar remain the same as for the simulated dataset. Indeed, the measurements are 2D polar points converted to Cartesian coordinates, and the detection probability is equal to p d = 0.9. Moreover, the range standard deviation is σ r = 1 meter and the angular standard deviation is equal to σ θ = 0.1 degrees. The false alarm rate is low at N F A = 5 false alarms per radar swipe, and the sampling period of the sensor is ∆ t = 1 second. Finally, it should be noted that the ground truth is empirical, which explains the drift in RMSE at the end of the scenario. Once again, the PMB filter struggles with the track continuity. Regarding the ET-PMB and the GT-PMB filters, they are on par with each other. Indeed, this scenario is composed of slowly evolving groups, and the false alarm rate is low. Moreover, the sensor noise error covariance is also low since the targets are moving close to the sensor: it can be approximated as negligible in regard to the extent in this specific scenario. Thus the derivation and update procedures of the original Random Matrix model 2 are sufficient here. Hence, the cardinal and errors are very similar on the graphs presented in figure 6.

CONCLUSION

RFSs have proven to be a very efficient and powerful way of solving the general Multi Target Tracking (MTT) with a Bayesian approach. However, they assume point targets, an hypothesis which can be questioned for many applications. In this work, when point targets are evolving as a groups, high cardinal estimation errors occur, in addition to track discontinuities. The Random Matrix (RM) approach allows to tackle this problem, hence it is implemented in the final solution of this work.

Both of the Random Matrix recursion presented in this paper are performing well in high clutter scenarios. The Group Target Poisson Multi Bernoulli implementation we proposed has proven to be very helpful when the sensor noise is non-negligible with respect to the extent, as discussed in the result part. Unfortunately, the sensor noise from the real dataset is low, hence the two approaches retained here perform similarly in this case.

Since these two approaches lead to labelled approaches, a comparative study using the Generalized Labelled Multi Bernoulli (GLMB) filter with Random Matrix based state probability densities could be interesting. Another extent to this work could be to test several evolution models for the extent and kinematic vector: using interacting multiple models to describe the evolution of the extent should lead to better results.

Figure 3 .

 3 Figure 3. Block diagram of the full Group Target Poisson Multi Bernoulli filter.

Figure 4 .

 4 Figure 4. Scenarios presented in this paper. Scenario a corresponds to the simulated dataset ground truth and scenario b corresponds to the real dataset ground truth.

Figure 5 .

 5 Figure 5. Metrics of the simulated dataset. One line corresponds to one scenario: a-c for σr = 1 meter and σ θ = 0.1 degrees, and d-f for σr = 4 meters and σ θ = 0.6 degrees

Figure 6 .

 6 Figure 6. Metrics for the real dataset.

Table 1 .

 1 Table of the notations used throughout this paper.

	x k	: State vector placeholder	z k	: Single measurement at time k
	X k	: Extent matrix placeholder	Z k	: Empirical measured extent at time k
	X k	: Set of targets at time k	z {1:k} /Z {1:k} : Single/Set of measurements up to time k
	p k|k-1 (•)	: Prior probability density	p k|k (•)	: Posterior probability density
	f k|k-1 (•)	: Markov transition density	l k (•)	: Measurement likelihood
	A k|k-1	: Any predicted parameter	A k|k	: Any updated parameter
	p d	: Detection probability	p s	: Survival probability
	N (x; µ, P)	: Gaussian probability density G(γ; α, κ)	: Gamma probability density
	W(Z; m, X) : Wishart probability density	IW(X; ν, V ) : Inverse Wishart probability density

  d.Figure 1. Association problem example with 2 targets and 2 measurements. n.d. stands for not detected and n.r. for not relevant. In General, using a PMB density as a prior leads to a PMBM density as posterior: LBP, and Munkres algorithms alleviate this problem.

	0 (n.r.)
	0 (n.r.)
	Possible posterior MB densities

  Birth intensity localization example. At t = 0, I β is predicted, then during the update, a measurement arrives and a new Bernoulli component is created. At t = 1 the cycle loops and pruning is applied to I U during the update. For the sake of clarity, speed vectors are omitted.

	Initialization	Update	Prediction	Update	Component of I β ,
					the birth intensity
					Component of I U ,
					the undected intensity
					Detected target
					as a Bernoulli RFS
					Measurement
					Sensor field of view
			1
	Figure 2.

Table 3 .

 3 Summary of the prediction and correction procedures.

ACKNOWLEDGMENTS

The authors want to thank the ONERA, Renault, and the Gipsa-Lab for their financial support.