Novel Bayesian Track-Before-Detection for Drones Based VB-Multi-Bernoulli Filter and a GIGM Implementation

RADIOENGINEERING, VOL. 29, NO. 2, JUNE 2020 397 DOI: 10.13164/re.2020.0397 SIGNALS Novel Bayesian Track-Before-Detection for Drones Based VB-Multi-Bernoulli Filter and a GIGM Implementation Ibrahim M. SALIM1, Mohamed BARBARY2, Mohamed H. ABD ELAZEEM1 1 Dept. of Electronics and Communication, Arab Academy for Science Technology and Maritime Transport, Egypt 2 Dept. of Electrical Engineering, Alexandria University, Egypt [email protected], mbarbary3[email protected], m[email protected] Submitted August 26, 2019 / Accepted March 10, 2020 Abstract. Joint detection and tracking of drones is a challenging radar technology; especially estimating their states with unknown measurement variances. The Bayesian track-before-detect (TBD) approach is an efficient way to detect low observable targets. In this paper, we proposed a new variational Bayesian (VB)-TBD technique for drones based on Multi-Bernoulli filter, which implemented with unknown measurement variances. Current implementation includes an analytical Gaussian inverse Gamma mixtures solution, which applied to estimate augmented kinematic drones state under the same circumstance. The results demonstrate that the proposed filter is more accurate than other Multi-Bernoulli filters in cardinality estimation. The proposed algorithm estimates the fluctuated parameters for each drone and it has no difficulty in handling the crossing of multiple drones. The Optimal Subpattern Assignment (OSPA) distances of proposed algorithm under different SNR are less than the other filters. It can be seen that at SNR (–5dB), the proposed algorithm and the other filters settle to errors 51 m, 125 m and 200 m, respectively. Keywords Drones tracking, Track-Before-Detect (TBD), MultiBernoulli filter, Variational Bayesian (VB) approximation 1. Introduction In recent years, cardinality-balanced multi-target multi-Bernoulli (CBMeMBer) filter has been widely used for multi-target tracking with excellent performance [1–5]. The CBMeMBer recursion is a tractable approximation of Bayesian multi-target recursion under low clutter density. It directly propagates approximate posterior density of the targets. The key advantage of multi-Bernoulli lies in extracting the reliable and inexpensive state estimates. This algorithm only performs well under given detection parameters; otherwise, it may decrease the tracking performance greatly. In this case, the common approach applies a threshold and treats those cells of exceeding the threshold. It is acceptable if signal-to-noise ratio (SNR) is high. The suitability of conventional radar systems to detect and identify micro-drones is a matter of significant importance to national security. It is being investigated while the use of such platforms is becoming more and more widespread. This is expected to be challenging, as microdrones have low radar cross-section (RCS), fly at low altitude, and speed in comparison with more conventional targets for military applications [6–10]. For low observable targets such as drones, there is a considerable advantage of using unthresholded data in simultaneous detection and tracking, which is known as track-before-detect (TBD). Many methods have been proposed to deal with the multitarget filtering under given parameters of low SNR [2], [3], [9]. However, they all have the same problem of larger computation and only adapted slowly to measurement variances. Normally the measurement variances are unknown and time-varying in real small targets tracking scenarios such as drones. In [2], the authors proposed a MultiBernoulli based on TBD filtering solution that can accommodate a linear Gaussian function and unknown detection profile, but with a given measurement variance. Recently, tracking filters based on variational Bayesian (VB) approximation method [11–17] have been applied to estimate the state of linear Gaussian process whose measurement variances are unknown. However, these algorithms depend on constant detection probability, which is unsuitable for drones tracking. In this paper, a new VB-CBMeMBer filter with Bayesian TBD is proposed to cope with jointly unknown detection probability and measurement variances. The detection profile is prior unknown due to variation of drone parameters such as small RCS [8–10]. Therefore, we introduce a VB approximation into the framework of CBMeMBer-TBD filter instead of the recent filters [2], [11–17]. Further, we derive a closed-form solution for time-varying drones tracking in the circumstance of jointly unknown detection parameters. Moreover, a Gaussian Inverse Gamma mixtures implementation is applied to estimate the augmented kinematic state. The rest of the paper is organized as follows. In Section 2, we give a formulation of drones’ TBD, and in Section 3 we present the updated formula of new recursions. The analytical implementation of GIGM is given in Section 4 and the numerical results and conclusions are given in Section 5, 6. 398 I. SALIM, M. BARBARY, M. ABD ELAZEEM, NOVEL BAYESIAN TRACK-BEFORE-DETECTION FOR DRONES … 2. Problem Formulation of Drones’ TBD Suppose at time k, there are N(k) drones states xk,1,⋯, xk,N(k), taking values in a state space χ , the dynamics of drone target ℓ are generally represented as   ,1, ,, 1,2,, kkkk F wNk    xx (1) where xk,ℓ=󰇟xk,ℓẋk,ℓ yk,ℓẏk,ℓIk,ℓ󰇠 is the state of drone target ℓ and comprises the position, velocity and intensity (RCS) of the drone target at instant k, Fk denotes the state transition function and wk is the process noise. Then, drones state (Xk) at time k, are finite sets, Xk = {xk,1,⋯, xk,N(k)} ⊂ χ. We now describe a drones’ observation model for TBD. Drone target returns measured by the radar are assumed to fluctuate according to the return amplitude fluctuation models [13]. The 2D images of the surveillance region (with nx   ny resolution cells) provided by the sensor can be denoted as [4]:     1 0 () ,, ,, , 1 , : Target is present, : target prese tNno ( )() , ,. Nk ij ij kk kk ij k ij k A            H H z  xh x (2) where A(xk,ℓ󰇜 is the complex echo of drone target x k,ℓ, hk (i,j)(xk,ℓ󰇜 is a measurement matrix and represents the contribution intensity of drone target ℓ to the cell (i,j), νk (i,j)  (0,Rk) is a Gaussian measurement noise cell, and the covariance matrix Rk is diagonal. If the point drone is assumed, the drone intensity can be approximated as (3) according to the point spread function [2]:     2 2 ,, , ,22 exp 2Σ2Σ xk yk ij xy kk ix jy             hx (3) where Σ is the amount of blurring introduced by the sensor. The complete set of measurements at instant k and up to instant k is denoted as Zk = {zk (i,j):i = 1…nx, j = 1…ny}. Assuming of an independent measurement noise from cell to cell. Here the fluctuation models are incorporated into the likelihood function, to account for the target return fluctuations. Under the further assumption of Gaussian background noise, the PDFs can be expressed as [2]:         ,, , 1,, , 1 ;, , Nk ij ij ij kkkkk ij A         z Hz xhxR       ,, 0 ,;0, ij ij kk ij   z HzR. (4) The drones state at time k is naturally modelled by a Random Finite Set (RFS) Xk = {xk,1,⋯, xk,N(k)} ⊂χ, χ  n. The surveillance region is divided into D resolution cells denoted by V1,V2,…,VD  n/2, a single drone target with state x occupies a set of resolution cells denoted by U(x). Here, we only consider the case that drones are rigid bodies. It means that the regions affected by different drones without overlap, i.e., U(x)  U(x’) = 0 if x x’. Assume that the measurements in different cells are mutually independent condition of the drones state Xk. The proposed probability density of measurement data Zk conditioned on the drones state Xk at time k by the two likelihoods is given by:          ,, ,, ,, (,) , |,, |,, | , kk k kkkk ij ij kk ij ij ijU ij U kz fgA g A                    xX xX x x x X ZXRA ZR x zx z R RR (5) where          , , , ,, |,, , , , | ij k ij zij ijU k ij A gA     x zxR xR zR    , , ,1 (,) | Dij kk ij ij f   zZR R . (6) The drone target state prior and process noise are assumed to be known, while the drone target measurement variance Rk and probability of drone target detection pD,k(x) are assumed unknown due to fluctuation of drone RCS. In addition, the process noise, Rk, as a function of initial target state value are independent of each other. The objective of this paper is to propose a state estimation algorithm for the nonlinear system with jointly unknown Rk. Unlike the conventional state estimation, what we considered here is to infer sequentially unobserved state xk together with jointly unknown measurement variances statistics Rk according to a set of measurements. The variational Bayesian (VB) approximation method has been applied to obtaining the state estimation of non-linear systems with unknown Rk. The main idea is to approximate the jointposterior distribution of the states and Rk by factorising the fixed form distribution and constructing the recursive expressions. 3. Variational Bayesian MultiBernoulli-TBD Recursions The basic idea of this improvement approach is to joint unknown drones parameters (i.e. measurement variance (R)) with unknown detection probability (a) of TBD scheme into a single drone target state. Estimating the augmented state of drones will yield information of particular drone target, such as number of drones, as well as individual kinematic states, jointly unknown R. In the following, we derive VB Multi-Bernoulli-TBD recursion for augmented state model, which propagates cardinality distribution of drones and intensity function of augmented state. We propose a modified filter with unknown parameters based on VB and TBD algorithm. It can be derived by the original filter in a specially chosen state space. The drones’ profiles are accommodated by incorporating unknown R with the state variable (x). It is achieved by putting a kinematic state into augmented variable (corre- RADIOENGINEERING, VOL. 29, NO. 2, JUNE 2020 399 sponding to unknown R). The filter should be able to estimate unknown drones’ parameters according to how well the measurements fit to an underlying drone target model. An explicit specification of single drone target model with new augmented state space is given with  ⁄,, 󰆸 ,  ⁄| 󰆸,  ⁄󰇛|󰇜, where  ⁄,, 󰆸 is the transition density at time k, for drone target augmented state with unknown drone target parameters , given previous value  󰆸,. 3.1 VB-CBMeMBer-TBD Prediction A multi-Bernoulli RFS [1–3] is united by a fixed number of independent Bernoulli RFS with existence probability r(i)  (1,0) and probability density p(i), i = 1,…,M, where M is the number of Bernoulli RFS. Thus a multiBernoulli RFS is completely described by the multi-Bernoulli parameter set  = {r(i),P(i)(x,R)}M i=1. Assuming that the state vector x and measurement covariance R are independent, and the survival probability are independent of them, so ps,k(x,R) = ps,k, we also assume that posterior multi-target density at time step k – 1 is represented by multi-Bernoulli parameter set as [16]:  󰇥󰇡  󰇛󰇜 , 󰇛󰇜 , 󰆸󰇢󰇦  (7) where r(i) k–1 is existence probability and p(i) k–1(.) is state distribution of the ith Bernoulli component. Mk–1 denotes the number of posterior hypothesized tracks at time k – 1, then the predicted multi-target density is multi-Bernoulli given by a union of birth and persisting or predicted components [16]:             1Γ, /1 ,/1 ,/1 Γ,Γ, 11 ,, kk MM ii ii kk pkk pkk k k ii rp rp      (8) where       1 ,/ 1 ,/ 1 1 ,k M ii pk k pk k i rp    and       Γ, Γ,Γ,1 ,k M ii kk i rp  denote the parameter set of multi-Bernoulli RFS for the surviving targets and spontaneous births, respectively. Mk–1 and MΓ,k denote predicted hypothesized track number of surviving targets and spontaneous births, respectively, then    ,/ 1 1 1 , , iii pk k k k sk rrpp  , (9)     1, 1 ,1 , 1 ,|., (, ) , i kk sk k i pk k i sk k fpp p pp    xR xR (10) where , represents inner product operation. The total number of predicted hypothesized tracks is Mk/k–1 = Mk–1 + MΓ,k. 3.2 VB-CBMeMBer-TBD Update If at time k, the predicted multi-target density is multiBernoulli in form of k/k–1=       /1 /1 /1 1 ,kk M ii kk kk i rp    , given the measurement vector of an image observation at time k, (Zk), the proposed updated Multi-Bernoulli parameters are        /1 1 , , ; kk k M ii kk k zi rp     Z zxRz (11) where               |1 |1 |1 |1 |1 ,,, , 1,, , ii ikk kk kiii kk kk kk rp rrrp    xR xR zxR xR     , 1 , (, ) ,(,) i iDk ki z kk | p; pg  Pz zxR xR xR (12) where P(i) D,k(x,RZk) = p(i) k/k–1(x,R)gz(x,R). Note that, R and ,  ⁄| 󰆸 are unknown in detail. Moreover, when the latest measurement Zk is available, the joint posterior distribution P(i) D,k(x,RZk) are intractable analytically and we cannot arrive at the solution. In order to calculate the posterior distribution with unknown Rk, the VB approximation method [16] is proposed to achieve the approximation solution. Now, assume that P(i) D,k(x,RZk) can be approximated by two intensity functions ,󰇛󰇜 and ,󰇛󰇜, while ,󰇛󰇜 and ,󰇛󰇜 are approximate posterior for x and R respectively. Therefore, P(i) D,k(x,Rz) ,󰇛󰇜,󰇛󰇜. In order to obtain the best approximation of intensity function P(i) D,k(x,Rz), the variational calculus is employed to minimize following Kullback–Leibler (KL)- divergence as [12–14, 16]:    ,, ,, arg m ˆin , ˆ kk kk  xR xR qxqR q q R (13) where           ,, ,, , , , , log d , d. , | kk k i kk ki D Dk k KL z |           P Pz R xR x xR qxqR qxqR qxqR x R R xR x R (14) The solution is obtained iteratively by optimization with respect to only one of the multiplicative factors in 󰇛∙󰇜 and fixing all the others to their last estimated values. The analytical solutions for  ,󰇛󰇜 and  ,󰇛󰇜 with such a procedure are given as follows [13]:     , ,,1 ˆ:1 log log , , ˆ| k kDkk zc     P R x q x qx xR     , ,,1: ˆ1 log log , | ˆk kDkk zc     P x Rq R qR xR (15) where  and  are constants with respect to variables  and , respectively. The expected values on the right-hand sides of (17) are taken by using the last estimated versions of ,󰇛󰇜 and ,󰇛󰇜 to obtain the new values which yield a convergent recursion [13]. Owe to coupling of ,󰇛󰇜 and ,󰇛󰇜, they cannot be solved directly. However, similar to the derivation of variational calculus in [12], [13], we can obtain the results that intensity function ,󰇛󰇜 and ,󰇛󰇜 are integral forms in Gaussian and inverse Gamma distributions, respectively, as, 400 I. SALIM, M. BARBARY, M. ABD ELAZEEM, NOVEL BAYESIAN TRACK-BEFORE-DETECTION FOR DRONES … ,󰇛󰇜󰇛; ,󰇜d,   2 ,,,,, 1 Inv - Gamma ; , d d kklklklkl l           R qR (16) where superscript d denotes the dimension of . Inv-Gamma(;α,β) denotes the inverse Gamma distribution of variable, whose degree of freedom is  and scalar parameter is , so   1 Inv Gamma ; , exp Γ           R R R ,    1 0 exp d ,ttt          ,1 ,1 , , diag / , , / ˆkkkkdkd   R . (17) 4. GIGM Implementation Based TBD In the following, the Gaussian Inverse-Gamma mixtures (GIGM) implementation of the VB-Multi-Bernoulli TBD (VB-CBMeMBer-TBD) recursion is described. The solution is obtained on GIGM-TBD distributions, the model of a proposed pdf is represented by,                 1 ,,,,, ,, 11 ~, , ,,,,, k k i k M ii kkk i M J iijijijijij k k k k kl kl ji rp rwm p         XxR (18) where Jk (i) is number of GM components of the ith drone; wk (i,j) is weight of the jth component; αk,l (i,j) and βk,l (i,j) are Inverse-Gamma distribution’ parameters used for estimating unknown Rk estimation; mk (i,j) and pk (i,j) are Gaussian mean and covariance of the jth component, respectively; Furthermore, suppose in an intermediate stage at time k – 1 and estimation of CBMeMBer can be approximated by an unnormalized mixture of GIGM-TBD distributions, we have:             1 , 2 ,,, 1, 1, 1, 1 ,,, 1111 1 Inv - Gamma ; , ; , , ij i k dij ij ij kl kl iij k ij ij kk l j l kk pwNmp               J xR x (19) where , 󰇛,󰇜 and , 󰇛,󰇜 are Inverse-Gamma distribution’ parameters which using for the unknown Rk–1 estimation, and      ,, ,1,1 1, 1,, 1,1 1, diag , , ii ikkd kii kkd          R    . Furthermore, suppose in an intermediate stage at time k – 1 and estimated MultiBernoulli for drones with unknown parameters can be approximated by an unnormalized VB-CBMeMBer-TBD posterior density, we have:                 1 1 1 111 1 ,, , , 111 1,1, 111 ,~ , , ,, , . , k k i k M ii kkk i M J d iijijijij kkk klkl lji rp rw                 XR xR x (20) 4.1 GIGM VB CBMeMBer-TBD Prediction Suppose that the intensity of birth RFS is an unnormalized mixture of GIGM-TBD distributions. Thus, the posterior birth model of the ith drone target is a multi-Bernoulli distribution with parameter set of      Γ, ,, 1 ,,k M ii kk i rp   xR  where p(i) Γ,k, i = 1,…,MΓ,k is a GIGM-TBD distribution, Then, we have             , ,2 , ,,, , ,, ,, ,, ,, 1 ,,, 1 Inv - Gamma ; , ,; , i k ij J di iiji jijij kl jij kl kl l kkkk j pwNmp               xR x (21) where wΓ,k (i,j), αΓ,k,l (i,j) and βΓ,k,l (i,j) denote weights and parameters of the jth component of birth model corresponding to the ith drone target. The survival probability is state independent, i.e. pS,k(x) = pS,k. Then, predicted (multi-Bernoulli) multi-target posterior density   /1 , kk    XR               1Γ, ,/ 1 ,/ 1 Γ,Γ, 11 ,, ,, kk MM ii ii pk k pk k k k ii rp rp    xR xR can be computed by a birth model as     ,/ 1 1 , ii P kk k Sk rrp  , (22)             1 , ,,, ,/ 1 1 ,/ 1 ,/ 1 1 2 ,,, , / 1, , / 1, , / 1, 1 , ; , Inv Gamma ; , , i k ij J iijijij pkk k pkk pkk j d ij ij ij pk k l pk k l pk k l l pwNmp              xR x     ,, ,/ 1 1 1 , ij ij pkk k k mFm      T ,, ,/ 1 1 1 1 1 , ij ij pk k k k k k pFpF     Q        ,,,,,, , / 1, 1, , / 1, 1, , ij ij ij ij ij ij pkkl l kl pkkl l kl      (23) where l = 1,…,d and   , ,/ 1 diag ˆij Pkk R       ,, , , ,/ 1,1 ,/ 1,1 ,/ 1, ,/ 1, /,, / ij ij ij ij P kk Pkk Pkk d Pkk d      , ρl (i,j) is a fading factor parameter, that actually affects estimation of drone target variance, ρl (0,1]. When ρl is large, it represents a gentle drone variance fluctuation from time k – 1 to k, in contrast ρl is small, the variances will be nonstationary and have strong fluctuations[11–15]. 4.2 GIGM VB-CBMeMBer-TBD Updating Suppose that the predicted VB CBMeMBer-TBD posterior density        /1 /1 /1 /1 1 ,,, kk M ii kk kk kk i rp     XR xR is RADIOENGINEERING, VOL. 29, NO. 2, JUNE 2020 401 given at time k and each drone target probability density p(i) k/k–1, i = 1,…,Mk/k–1 is comprised of a GIGM-TBD distributions, we have.             / , 1 2 ,,, 1, ,,, /1 /1 /1 /1 1, 1, 1 1 Inv - Gam ,; ma . ;, , i i j kk J dij ij ij kk l kk l kk l l iijijij kk kk kk kk j pwNmp               xR x (24) Then, the approximation of VB CBMeMBer-TBD updating         /1 1 ,, , ; kk k M ii kkk i rp     zZ XR z xRz can be computed as follows:           |1 |1 |1 , , 1, ii ikk k kiii kk kk k rz rz rr z    R R e e                / , 1 /1 ,,, /1 1 , / 2 ,, 1 , 1 , ,, 1 ; , ,; Inv - Gamma ; . , i k j i kk i k Jij ij ij kk k k ij kJij dij ij ij kk j kl kl kl l wNm p p w                x xRz (25) The updating parameters can be obtained by following iteration calculations. First set      ,0 , /1 ij ij kkk mm  ,   ,0 , /1 ij ij kkk pp  ,     ,, ,/1, 0.5 ij ij kl kk l    and       ,0 , ,/1, ij ij kl kk l   , for l = 1,…,d. Then iteration is following        ,, , ,,1 ,2 , ,, , ,1 ,2 , diag , , , ij n ij n ij n ij n kk kd kij ij ij kk kd         R, (26)             /1 T ,,, , /1 /1 /1 1 , ; , , i kk in k J ij ij ij ij n kk k kk k kk k k j z WNzHm Hp H       R R e (27)       1 TT ,, , , /1 /1 , ij n ij ij ij n kkkkkkkkk K pHHpH      R (28)      ,,, , /1 /1 ij n ij ij n ij kkkkkkk mzmK zHm      , (29)       ,,, /1 1 ij n ij n ij kkkkk pKHp     , (30)           , T ,,, , /1 /1 /1 ,; ; , , ij n k ij ij ij ij n kk k kk k kk k k wz wNzHm Hp H    xR R (31)          2 ,1 , , ,, , 1 2 1 2 ij n ij n ij n kl kl k k i T ij n kk k ii zHm Hp H    (32) where n  [1,N], N denotes the maximum number of iterations. If the estimation state satisfies    ,,1 ij n ij n kk mm    , then iteration ends, and the parameters are obtained as   ,,ij ij n kk WW,       ,, ij ij n kk mm,       ,,ij ij n kk pp and       ,, ,, ij ij n kl kl  . Otherwise, it continues iterations from (26) – (32). 5. Numerical Results This section presents numerical results of a new VBCBMeMBer-GIGM-TBD filter. Linear TBD scenarios are used to illustrate and examine the performance against the recent VBCBMeMBer [15] and CBMeMBer-TBD [2] filters. In [6], the authors proposed a practical scenario to simulate radar back-scattering of moving drone targets. The goal was to collect and process RCS data of drones and birds at high frequencies (K-band and W-band). The experimental setup has been intentionally chosen such that it corresponds to a realistic scenario for a drone/bird detection radar system. A complex RCS is normally computed by coherently combining cross sections of the simple target shapes. In [8], [9], the authors proposed a facet model to simulate radar back-scattering of a moving complex target based on POFACET program. Three drones of different sizes (DJI Phantom 3 Standard, DJI Inspire 1 and DJI S900 Hexacopter) are shown in Fig. 1. The statistical distribution of RCS for each target falls within a certain range which is useful for predicting the performance of drone detection radar. The real drone RCS values can be used to reconstruct the drones’ signatures at any aspect in angular extent summarized in Tab. 1. The RCS data of different drone target models are used to predict true probability density function (pdf). Figure 1 shows the Cumulative Density Function (CDF) of statistical histogram method for different statistical models. Note that  for actual drones are unknown to the filter. Thus, the measurement data is simulated according to a state dependent  and a of each drone target model, whose measurement variance peaks at the edge of surveillance region and tapers off at the origin. We demonstrate the performance of VB-CBMeMBer-TBD filter and compare it with CBMeMBer-TBD, and VBCBMeMBer filters via following examples. Consider a two-dimensional scenario with an unknown and time varying number of drones in the noise. A time varying number of drones is observed in noise region with surveillance dimensions of [–250,250] m  [–250,250] m. Assume that there are three drones emerging and disappearing successively. The real drones’ trajectories are shown in Fig. 2. The sample image of a simulation scenario is shown in Fig. 3 at times k = 20, 60 s, it shows the simulated noisy data at various SNR levels based on RCS with the middle point drones. Three targets appear at the position with circle marked, but they are buried by the background noise. In following dynamical and measurement models, state T ,,,, ,, ,,,, k k xk xk yk yk k a ppppw      R comprises unknown detection probability ak, measurement variances Rk, vector of planar position and velocity T ,,,, ,,, xk xk yk yk pppp    , and turn to rate wk. 402 I. SALIM, M. BARBARY, M. ABD ELAZEEM, NOVEL BAYESIAN TRACK-BEFORE-DETECTION FOR DRONES … Fig. 1. Drones geometry model, Overlaid RCS histograms and the corresponding CDF plots: (a) DJI Phantom 3 Standard, (b) DJI Inspire 1, (c) DJI S900. No. Drones type RCS (dBsm*)  (m) 1 DJI Phantom 3 −20 10 m–6 m 2 DJI Inspire 1 −15 8 m–5 m 3 DJI S900 −8 5 m–2 m Tab. 1. Drones parameters and average modal RCS. * dBsm: Decibel relative to one square meter. No. Simulation parameters Value 1 , 0.99 2  5 m/s 2 3 ∆ 1 s 4 Observation region 250 m × 250 m 5 Equivalent image 250 × 250 pix 6 ∆∆ 1 m 7 Source intensity (󰇜 30 8 Blurring factor ( 󰇜 1 9 Weight threshold (󰆒) 100 10 merging threshold (󰆒󰇜 4 m 11 J max [15] 100 12   1 13 T max 100 14 weight threshold 󰇛󰇜 10 –3 15  [2] 5 m 16  0.95 17  0.7 18 OSPA distance parameters p = 1, c = 300 Tab. 2. Simulation parameters of the proposed filter. Fig. 2. True drones tracks in xy plane used in simulations. Fig. 3. Measurement frames of drones at different time steps. The single-drone target transition model is linear Gaussian and specified by 43 22 22 2 3 22 2 22 42 , 0 2 kkv II II III                FQ (33) where In and 0n denote n  n identity and zero matrices,  is sampling period, and  is standard deviation of the noise. The probability of target survival is pS,k. The birth CBMeMBerGIGM-TBD is set to make Gaussian mean vectors corresponding to actual starting points. The birth process of  , 3, 1, 2, 3 i k Ji   is a multi-Bernoulli RFS with density        31,2,3 Γ,Γ,Γ,, 1 , , 0.02 ii kkk k i rp r     ,       ,,,, ,; , iiii kkkk Px wNxmp  R ,     32 ,, ,, ,, 1 IG ; , d iii kl kl kl l        ,   1 , 150m, 0, 150m, 0 k m   ,   2 , 200m, 0m/s, 50m, 0m/s , k m     3 , 50 m, 0m/s, 200m, 0m/s k m   ,      1,2,3 ,,,, diag 100,1 ,1 00,1 , 0.2, 1, iii kkkk pw           ,1 ,1 ,2 ,2 ,3 ,3 diag / , / , / . iiiiiii kkkkkkk    R We use a common observation model for all the experiments. At each time step k, it produces an image consisting an array of cells with a scalar intensity. In these demonstrations, the observation region is a 250 m  250 m, and equivalent image is a 250  250 pix array. The array index will be treated as an order pair of integers (a,b), where 1  a, b  250. The measurement variance Rk is unknown for all filters. The point spread function, Hi,k, is given by (3), and (px,k, py,k) being the position of state xk. A typical observation is shown in Fig. 3. At each time step of VB-CBMeMBer-TBD and other filters, pruning and merging of Gaussian components are performed in each hypothesized track. We summarized the all parameters of simulation in Tab. 2. Additionally, pruning of hypothesized tracks is performed by weight threshold and maximum of tracks are kept. Taking VB-CBMeMBer [11], [12] and CBMeMBer- RADIOENGINEERING, VOL. 29, NO. 2, JUNE 2020 403 Fig. 4. True drones track in x and y directions versus time (s). TBD [2] filters for comparison, the same pruning and merging of Gaussian components are performed on the posterior intensity. The selection of fading factor ρl is difficult, because the stability of noise estimation is controlled by ρl. We find that the performance of CBMeMBer-TBD is better when ρl = 0.95 and R = 5 m [2], VB-CBMeMBer using a = 0.7 for all drones [16]. Simulation results are obtained by 200 Monte Carlo runs and shown in Figs. 4–7. The proposed filter complexity was measured by recording the computational time of each time step of the filtering algorithm. The mean and standard deviation of the run time for filter iteration (a single prediction and update step) in seconds are 4.156 s and 0.0916 s, respectively. Figure 4 shows the three drones trajectories (x and y coordinates vs. time), they are estimated in a single run by different filters with unknown detection probability and measurement variance. The results indicate that proposed filter provides accurate tracking performance. However, the other filters encounter occasional tracking false and drops because the detection probability is not known for VB-CBMeMBerGIGM filter, and the measurement variances is not known for CBMeMBer-TBD filter. These results demonstrate that the proposed filter is able to correctly track the motions of individual drone target and identify the various births and deaths throughout. The filter has no difficulty in handling crossings of two drones (DJI Inspire 1, DJI Phantom 3) at time t = 47, so as two drones (DJI S900, DJI Inspire 1) at time t = 56. The proposed algorithm estimates the fluctuated parameters for each drone target, it has better adaptive ability for jointly unknown drones parameters such as estimated  than other algorithms as shown in Fig. 5. Figure 6 shows the number of drone target estimated by VBCBMeMBer, CBMeMBer-TBD and the proposed algorithms. It is clear that the proposed algorithm has better performance in estimating drone target number than CBMeMBer-TBD and VB-CBMeMBer under unknown a and R. This is due to the fact that both a and R are correctly estimated by the proposed algorithm. The performances of other algorithms decrease greatly with erroneous σ and a values. The reason is that assumed CBMeMBer-TBD is inaccurate for value σ = 5 m. While, VB-CBMeMBer is neither accurate at a = 0.7, due to the same reason. Figure 7 Fig. 5. Estimated measurement variances. Fig. 6. The number of drones estimated.  Fig. 7. OSPA of cardinality against time (s). 404 I. SALIM, M. BARBARY, M. ABD ELAZEEM, NOVEL BAYESIAN TRACK-BEFORE-DETECTION FOR DRONES … Algorithms Cardinality/ Time(s) OSPA distance/Time 2 /25 s 3 /40 s 3 /55 s m /25 s m /40 s m /55 s VB-CBMeMBer-TBD 2 2.9 2.8 5 10 25 CBMeMBer-TBD 1.9 2.7 2.4 10 50 100 VB-CBMeMBer 1.8 2.5 2 25 100 180 Tab. 3. Comparison of algorithms at different tracking times. shows the comparison of OSPA distances. In former methods, the individual state estimates are means of corresponding posterior densities. We summarized the results of Figs. 4–7 in Tab. 3. 6. Conclusion This work considered the method of multi-BernoulliTBD filter for unknown drones’ detection parameters. To tackle the problems of drones tracking, an improved VBCBMeMBer-TBD algorithm is proposed in this paper. The VB-CBMeMBer-TBD is approximated by GIGM implementation. Therefore, it is capable of estimating the drone kinematic state augmented by unknown parameters. The simulation results show that the proposed algorithm has better performance than recently tracking filters. Acknowledgments Research was supported by AASTMT in Cairo. References [1] OUYANG, C., JI, H., LI, C. Improved multi-target multi-Bernoulli filter. IET Radar Sonar Navigation, 2012, vol. 6, no. 6, p. 458–464. DOI: 10.1049/iet-rsn.2011.0377 [2] VO, B.-N, VO, B.-T, PHAM, N., et al. Joint detection and estimation of multiple objects from image observations. IEEE Transactions on Signal Processing, 2010, vol. 58, no. 10, p. 5129–5141. DOI: 10.1109/TSP.2010.2050482 [3] RISTIC, B., VO, B.-T., VO, B.-N., et al. A tutorial on Bernoulli filters: Theory, implementation and applications. IEEE Transactions on Signal Processing, 2013, vol. 61, no. 13, p. 3406–3430. DOI: 10.1109/TSP.2013.2257765 [4] VO, B.-T., SEE, C. M., MA, N., et al. Multi-sensor joint detection and tracking with Bernoulli filter. IEEE Transactions on Aerospace and Electronic System, 2012, vol. 48, no. 2, p. 1385–1402. DOI: 10.1109/TAES.2012.6178069 [5] VO, B.-N, VO, B.-T, HOSEINNEZHAD, R., et al. Robust multiBernoulli filtering. IEEE Journal of Selected Topics in Signal Processing, 2013, vol. 7, no. 3, p. 399–409. DOI: 10.1109/JSTSP.2013.2252325 [6] RAHMAN, S., ROBERTSON, D. A. In-flight RCS measurements of drones and birds at K-band and W-band. IET Radar Sonar Navigation, 2019, vol. 13, no. 2, p. 300–309. DOI: 10.1049/ietrsn.2018.5122 [7] PATEL, J. S., FIORANELLI, F., ANDERSON, D. Review of radar classification and RCS characterisation techniques for small UAVs or drones. IET Radar, Sonar, and Navigation, 2018, vol. 12, no. 9, p. 911–919. DOI: 10.1049/iet-rsn.2018.0020 [8] ZONG, P., BARBARY, M. Improved multi-Bernoulli filter for extended stealth targets tracking based on sub-random matrices. IEEE Sensors Journal, 2016, vol. 16, no. 5, p. 1428–1447. DOI: 10.1109/JSEN.2015.2499268 [9] BARBARY, M., ZONG, P. A novel stealthy target detection based on stratospheric balloon-borne positional instability due to random wind. Radioengineering, 2014, vol. 23, no. 4, p. 1192–1202. [10] DUQUE DE QUEVEDO, A., IBANEZ URZAIZ, F., GISMERO MENOYO, J., et al. Drone detection and radar-cross-section measurements by RAD-DAR. IET Radar, Sonar, and Navigation, 2019, vol. 13, no. 9, p. 1437–1447. DOI: 10.1049/ietrsn.2018.5646 [11] YANG, J., GE, H. Adaptive probability hypothesis density filter based on variational Bayesian approximation for multi-target tracking. IET Radar, Sonar, and Navigation, 2013, vol. 7, no. 9, p. 959–967. DOI: 10.1049/iet-rsn.2012.0357 [12] WU, X., HUANG, G., GAO, J. Adaptive noise variance identification for probability hypothesis density-based multi-target filter by variational Bayesian approximations. IET Radar, Sonar, and Navigation, 2013, vol. 7, no. 8, p. 895–903. DOI: 10.1049/ietrsn.2012.0291 [13] SARKKA, S., NUMMENMAA, A. Recursive noise adaptive Kalman filtering by variational Bayesian approximations. IEEE Transactions on Automatic Control, 2009, vol. 54, no. 3, p. 596–600. DOI: 10.1109/TAC.2008.2008348 [14] XINBO GAO, DACHENG TAO, XUELONG LI, et al. Multisensor centralized fusion without measurement noise covariance by variational Bayesian approximation. IEEE Transactions on Aerospace and Electronic System, 2011, vol. 47, no. 1, p. 718–727. DOI: 10.1109/TAES.2011.5705702 [15] SCHUHMACHER, D., VO, B.-N, VO, B.-T. A consistent metric for performance evaluation of multi-object filters. IEEE Transactions on Signal Processing, 2008, vol. 56, no. 8, p. 3447–3457. DOI: 10.1109/TSP.2008.920469 [16] YANG, J., GE, H. An improved multi-target tracking algorithm based on CBMeMBer filter and variational Bayesian approximation. Signal Processing, 2013, vol. 93, p. 2510–2515. DOI: 10.1016/j.sigpro.2013.03.027 [17] QIU, H., HUANG, G., GAO, J. Variational Bayesian labeled multi-Bernoulli filter with unknown sensor noise statistics. Chinese Journal of Aeronautics, 2016, vol. 29, no. 5, p. 1378–1384. DOI: 10.1016/j.cja.2016.05.002 About the Authors ... Ibrahim SALIM was born in Cairo 1986. Salim received his B.Sc. in July 2008 from the Egyptian Military Technical College, Radar and Communication Dept. He is a Master’s student at AASTMT in Cairo.He is working to develop new trends in detection and tracking multi-target. Mohamed BARBARY received the B.Sc degree in Electronics and Communications from the Faculty of Engineering, Alexandria University, Egypt, in 2003, and the Ph.D. degree in College of Electronic and Information Engineering, NUAA, Nanjing-P.R. China, in 2016. His work is focusing on stealth targets detection and tracking. Mohammed ABD EL-AZEEM received his B.Sc. in 1985 from the Egyptian Military Technical College Communications Dept. and Ph.D. from University of Kent in 1996. He is currently a Vice Dean for engineering faculty and Professor in Communication and Electronic Dept. in AASTMT in Cairo. His work is focusing on security of communications systems and networks. 