# Robust Multisensor MeMBer Filter for Multiple Extended-Target Tracking

Authors: Xiaoke Lu; Zhiguo Zhang; Qing Li; Jinping Sun
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

This paper develops a robust extended-target multisensor multitarget multi-Bernoulli (ET-MS-MeMBer) filter for enhancing the unsatisfactory quality of measurement partitions arising in the classical ET-MS-MeMBer filter due to increased clutter intensities. Specifically, the proposed method considers the influence of the clutter measurement set by introducing the ratio of the target likelihood to the clutter likelihood. With the constraint of the clutter measurement set, it can obtain better multisensor measurement partitioning results under the original two-step greedy partitioning mechanism. Subsequently, the single-target multisensor likelihood function for the clutter case is derived. Simulation results reveal a favorable comparison to the ET-MS-MeMBer filter in terms of accuracy in estimating the target cardinality and target state under conditions with increased clutter intensities.


---

## Body

## 1. Introduction

Multiple-target tracking (MTT) [1] estimates the number and states of moving targets based on the sensor observations. Traditional MTT algorithms adopt a point target observation model [2], which assumes that one target generates at most one measurement per time step. However, in the high-resolution sensor system, one target may occupy multiple-resolution cells and, thus, can produce multiple measurements. Such targets are known as extended targets [3]. The measurements generated by the extended target provide information not only about the target’s motion state but also the target’s shape. This shape information, e.g., the size and structure of the extended target, can further facilitate the target classification and recognition task [4].The conventional multiple-target tracking algorithm involving a data association step [5] has been widely used in the fields of medical image processing [6], network defense [7], and battlefield environment monitoring [8]. However, for extended-target tracking, the association between measurements and targets becomes more complicated and suffers from a “combinatorial explosion” problem. Therefore, the random finite set (RFS) [9, 10] has received much attention in the MTT domain for the merit of avoiding explicit data association steps. This theory leads to the development of various MTT algorithms, including the probability hypothesis density (PHD) [11], cardinalized PHD (CPHD) [12], arithmetic average multi-Bernoulli (AA-MB) [13], and multisensor multitarget multi-Bernoulli (MS-MeMBer) [14] filters.The abovementioned filters based on the RFS have further been applied to the field of extended-target tracking [15–23]. Based on the PHD filter, an extended-target PHD (ET-PHD) filter was proposed in [15], and a Gaussian implementation of it was given in [16]. As the ET-PHD filter can only estimate the centroid state of the extended target, the PHD filters with shape estimation were proposed in [17, 18]. In [17], a Gaussian inverse Wishart PHD (GIW-PHD) filter was proposed, where the target kinematical state is modeled as a Gaussian distribution and the target extension is assumed to follow an inverse Wishart distribution. As an improvement to the GIW-PHD filter, the GGIW-PHD filter proposed in [18] can additionally obtain the estimation of target measurement rates. Besides the PHD filter, other filters based on the RFS [19–23] have also been modified to track the extended target. In [19], an extended-target CPHD filter that can estimate the centroid state of extended targets was proposed. Compared with the ET-PHD filter, this filter achieves better filtering results at the cost of higher computational complexity. Based on the labeled RFS, the Gamma Gaussian inverse Wishart (GGIW) implementation of the labeled multi-Bernoulli (LMB) filter was introduced in [20]. In [21–23], the Poisson multi-Bernoulli Mixture (PMBM) filter was generalized to the extended-target tracking case.For the multisensor case, the extended-target MS-MeMBer (ET-MS-MeMBer) filter and its Gaussian inverse Wishart (GIW) mixture implementation was proposed in [24]. The ET-MS-MeMBer filter assumes that the observation of extended targets obeys the approximate Poisson-Body (APB) model [25], and the update process of the MS-MeMBer filter is modified accordingly. However, this method ignores the constraint of the clutter measurement set in the multisensor measurement partitioning process, which may degenerate the quality of the multisensor measurement partition and, consequently, the filtering performance of the MS-MeMBer filter in high-clutter-density scenarios. To solve this problem, a robust ET-MS-MeMBer filter is proposed in this paper. By introducing the ratio of the target likelihood to the clutter likelihood, the influence of the clutter measurement set is considered. Therefore, better measurement partitions can be obtained under the original two-step partitioning mechanism [14]. The single-target multisensor likelihood function for the clutter case is also derived. Simulation results show that compared to the ET-MS-MeMBer filter, the proposed filter has higher accuracy in estimating the target cardinality and target state under high-clutter-density conditions.The rest of this paper is organized as follows.Section 2 gives a review of the ET-MS-MeMBer filter. The proposed filter is derived in Section 3. The simulation results are presented in Section 4, and the conclusions are given in Section 5.


## 2. The ET-MS-MeMBer Filter

### 2.1. The Prediction Step

Let x andX=x1,…,xMk−1|k−1 denote a single-target state and a multitarget state, where Mk−1|k−1 is the cardinality of targets at time k−1. Assume that, at time k−1, the multitarget posterior density is a multi-Bernoulli(1)πk−1|k−1X=rk−1|k−1i,pk−1|k−1ixi=1Mk−1|k−1,where each rk−1|k−1i,pk−1|k−1ix represents an independent Bernoulli density with probability rk−1|k−1i and distribution pk−1|k−1ix.Then, the predicted multitarget density is also a multi-Bernoulli(2)πk|k−1X=rs,k|k−1i,ps,k|k−1ixi=1Mk−1|k−1∪rb,kj,pb,kjxj=1Mb,k,where rb,kj,pb,kjx and Mb,k are the Bernoulli density and the cardinality of newborn targets at time k. The survival Bernoulli density rs,k|k−1i,ps,k|k−1ix follows(3)rs,k|k−1i=rk−1|k−1ipk−1|k−1i,ρsv,ps,k|k−1ix=fk|k−1x|⋅,ρsvpk−1|k−1ixpk−1|k−1ix,ρsv,where ρsv is the probability that a target can survive to the next time, ⋅,⋅ is the inner product function, and fk|k−1x|⋅ is the state transition function.


### 2.2. The Distance Partition and Multisensor Measurement Partition

The update step of the ET-MS-MeMBer filter requires the multisensor measurement partitions. Since each extended target may generate multiple measurements at each scan, our first task is to divide each sensor’s measurements into groups that contain measurements from the same target. Here, we utilize a distance partition method [17] which is based on the general belief that measurements from the same extended target are spatially close to each other.LetZi=zi1,…,zin be the collection of all measurements generated on sensor i at time k, where zil1≤l≤n is the l-th element in Zi and n is the number of the measurements. The distance partition of Zi can be denoted as Di=wi0,wi1,…,wiDi−1, where(4)wil∩wim=∅0≤l≠m≤Di−1,wi0∪wi1∪⋯∪wiDi−1=Zi.⋅ is the cardinality function. Each disjoint subset wij (j=0 for the clutter measurement subset and 0<j≤Di−1 for the target measurement subset) in partition Di contains several measurements in Zi and satisfies wij≥0. Di is the collection of all possible distance partitions Di.Multisensor measurement partition aims to group all sensors’ measurements from the same target into the same subset. We consider a multisensor measurement setZ1:S=Z1,…,ZS, where S is the number of sensors. Suppose that, at time k, the predicted multitarget density contains Mk|k−1 Bernoulli components. For a given set of distance partitions D1,…,DS, the multisensor measurement partition can be defined as P=W1:S0,W1:S1,…,W1:SMk|k−1, where(5)W1:Sl=w1l1,…,wSlS1≤l≤Mk|k−1,andw1l1∈D1,…,w1lS∈DS,W1:Sl∩W1:Sj0≤l≠j≤Mk|k−1,W1:S0∪W1:S1∪⋯∪W1:SMk|k−1=Z1:S.Here,W1:S0 is the clutter measurement set generated by all S sensors, and W1:Sll>0 is the target measurement set corresponding to the l-th Bernoulli component. P is defined as the collection of all possible multisensor measurement partitions P. Finally, we define the mapping function TW1:Sl=i,li|wili∈W1:Sl, where i,li specifies the sensor label i and the index li of wili in Di.


### 2.3. The Update Step

Suppose that, at timek, the predict multitarget density is a multi-Bernoulli consisting of Mk|k−1 independent Bernoulli component:(6)πk|k−1X=rk|k−1j,pk|k−1jxj=1Mk|k−1.Then, the posterior density can be approximated by a multi-Bernoulli as follows:(7)πk|kX=∪P∈P∪j=1Mk|k−1rk|kP,j,pk|kP,jx.The updated probabilityrk|kP,j and the distribution pk|kP,jx follow(8)rk|kP,j=αPrk|k−1jpk|k−1j,η1−rk|k−1j+rk|k−1jpk|k−1j,η,W1:Sj=∅,αP,W1:Sj≠∅,pk|kP,jx=pk|k−1jxηxpk|k−1j,η,W1:Sj=∅,pk|k−1jxfW1:Sj|x∫pk|k−1jxfW1:Sj|xdx,W1:Sj≠∅,where ηx is the probability that a target is not detected by any sensor and fW1:Sj|x is the single-target multisensor likelihood function. αP is the weight of P in P, which can be calculated as the score of P divided by the sum of all partitioning scores in P:(9)αP=KP∏j=1Mk|k−1φW1:Sjj1∑M∈PKM∏j=1Mk|k−1φW1:Sjj1,where(10)φW1:Sjju=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxfW1:Sj|xdx,W1:Sj≠∅,KP=∏i=1SCiwi00×∏j=1Mk|k−1∏i,li∈TW1:Sjwili!.KP and φW1:Sjj1 quantify how likely the set W1:S0 is a clutter set, and the set W1:Sjj>0 is associated with the j-th Bernoulli component, respectively. Cinx is the n-th derivative of the clutter probability generating function Cix=eλix−λi, and λi is the clutter intensity of sensor i.


## 2.1. The Prediction Step

Let x andX=x1,…,xMk−1|k−1 denote a single-target state and a multitarget state, where Mk−1|k−1 is the cardinality of targets at time k−1. Assume that, at time k−1, the multitarget posterior density is a multi-Bernoulli(1)πk−1|k−1X=rk−1|k−1i,pk−1|k−1ixi=1Mk−1|k−1,where each rk−1|k−1i,pk−1|k−1ix represents an independent Bernoulli density with probability rk−1|k−1i and distribution pk−1|k−1ix.Then, the predicted multitarget density is also a multi-Bernoulli(2)πk|k−1X=rs,k|k−1i,ps,k|k−1ixi=1Mk−1|k−1∪rb,kj,pb,kjxj=1Mb,k,where rb,kj,pb,kjx and Mb,k are the Bernoulli density and the cardinality of newborn targets at time k. The survival Bernoulli density rs,k|k−1i,ps,k|k−1ix follows(3)rs,k|k−1i=rk−1|k−1ipk−1|k−1i,ρsv,ps,k|k−1ix=fk|k−1x|⋅,ρsvpk−1|k−1ixpk−1|k−1ix,ρsv,where ρsv is the probability that a target can survive to the next time, ⋅,⋅ is the inner product function, and fk|k−1x|⋅ is the state transition function.


## 2.2. The Distance Partition and Multisensor Measurement Partition

The update step of the ET-MS-MeMBer filter requires the multisensor measurement partitions. Since each extended target may generate multiple measurements at each scan, our first task is to divide each sensor’s measurements into groups that contain measurements from the same target. Here, we utilize a distance partition method [17] which is based on the general belief that measurements from the same extended target are spatially close to each other.LetZi=zi1,…,zin be the collection of all measurements generated on sensor i at time k, where zil1≤l≤n is the l-th element in Zi and n is the number of the measurements. The distance partition of Zi can be denoted as Di=wi0,wi1,…,wiDi−1, where(4)wil∩wim=∅0≤l≠m≤Di−1,wi0∪wi1∪⋯∪wiDi−1=Zi.⋅ is the cardinality function. Each disjoint subset wij (j=0 for the clutter measurement subset and 0<j≤Di−1 for the target measurement subset) in partition Di contains several measurements in Zi and satisfies wij≥0. Di is the collection of all possible distance partitions Di.Multisensor measurement partition aims to group all sensors’ measurements from the same target into the same subset. We consider a multisensor measurement setZ1:S=Z1,…,ZS, where S is the number of sensors. Suppose that, at time k, the predicted multitarget density contains Mk|k−1 Bernoulli components. For a given set of distance partitions D1,…,DS, the multisensor measurement partition can be defined as P=W1:S0,W1:S1,…,W1:SMk|k−1, where(5)W1:Sl=w1l1,…,wSlS1≤l≤Mk|k−1,andw1l1∈D1,…,w1lS∈DS,W1:Sl∩W1:Sj0≤l≠j≤Mk|k−1,W1:S0∪W1:S1∪⋯∪W1:SMk|k−1=Z1:S.Here,W1:S0 is the clutter measurement set generated by all S sensors, and W1:Sll>0 is the target measurement set corresponding to the l-th Bernoulli component. P is defined as the collection of all possible multisensor measurement partitions P. Finally, we define the mapping function TW1:Sl=i,li|wili∈W1:Sl, where i,li specifies the sensor label i and the index li of wili in Di.


## 2.3. The Update Step

Suppose that, at timek, the predict multitarget density is a multi-Bernoulli consisting of Mk|k−1 independent Bernoulli component:(6)πk|k−1X=rk|k−1j,pk|k−1jxj=1Mk|k−1.Then, the posterior density can be approximated by a multi-Bernoulli as follows:(7)πk|kX=∪P∈P∪j=1Mk|k−1rk|kP,j,pk|kP,jx.The updated probabilityrk|kP,j and the distribution pk|kP,jx follow(8)rk|kP,j=αPrk|k−1jpk|k−1j,η1−rk|k−1j+rk|k−1jpk|k−1j,η,W1:Sj=∅,αP,W1:Sj≠∅,pk|kP,jx=pk|k−1jxηxpk|k−1j,η,W1:Sj=∅,pk|k−1jxfW1:Sj|x∫pk|k−1jxfW1:Sj|xdx,W1:Sj≠∅,where ηx is the probability that a target is not detected by any sensor and fW1:Sj|x is the single-target multisensor likelihood function. αP is the weight of P in P, which can be calculated as the score of P divided by the sum of all partitioning scores in P:(9)αP=KP∏j=1Mk|k−1φW1:Sjj1∑M∈PKM∏j=1Mk|k−1φW1:Sjj1,where(10)φW1:Sjju=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxfW1:Sj|xdx,W1:Sj≠∅,KP=∏i=1SCiwi00×∏j=1Mk|k−1∏i,li∈TW1:Sjwili!.KP and φW1:Sjj1 quantify how likely the set W1:S0 is a clutter set, and the set W1:Sjj>0 is associated with the j-th Bernoulli component, respectively. Cinx is the n-th derivative of the clutter probability generating function Cix=eλix−λi, and λi is the clutter intensity of sensor i.


## 3. The Robust ET-MS-MeMBer Filter

### 3.1. The Constraint of the Clutter Measurement Set

In [14], a two-step greedy partitioning mechanism is proposed, which can obtain several high-scoring measurement partitions. In this mechanism, the target measurement set corresponding to each Bernoulli component is first obtained separately, and then, the target measurement sets corresponding to different Bernoulli components are combined to obtain all possible measurement partitions. However, it can be seen from (9) that the two-step greedy partitioning mechanism ignores the influence of the clutter measurement set W1:S0. Considering the constraint of the clutter measurement set, αP can be rewritten as(11)αP=KP/∏i=1SλiZie−λi∏j=1Mk|k−1φW1:Sjj1∑M∈PKM/∏i=1SλiZie−λi∏j=1Mk|k−1φW1:Sjj1,=∏j=1Mk|k−1φW1:Sjj1⋅∏i,li∈TW1:Sjwili!/λiwili∑M∈P∏j=1Mk|k−1φW1:Sjj1⋅∏i,li∈TW1:Sjwili!/λiwili.We introduce the function(12)φ˜W1:Sjj1=φW1:Sjj1⋅∏i,li∈TW1:Sjwili!λiwili,=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxfW1:Sj|x∏i,li∈TW1:Sjwili!λiwilidx,W1:Sj≠∅,which denotes the ratio of the likelihood that W1:Sj is generated by a target to the likelihood that W1:Sj is from the clutter.According to [24], the single-target multisensor likelihood function follows(13)fW1:Sj|x=∏i,li∈TW1:Sjpi,dγiwilie−γiwili!∏z∈wilihiz|xciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi,where γi is the mean of the extended measurements, pi,d is the target detection probability by sensor i, hiz|x is the measurement likelihood of sensor i, and ciz is the clutter intensity.Hence, (12) becomes(14)φ˜W1:Sjj1=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxf˜W1:Sj|xdx,W1:Sj≠∅.In (14), f˜W1:Sj|x denotes the single-target multisensor likelihood function for the clutter case and follows(15)f˜W1:Sj|x=∏i,li∈TW1:Sjpi,dγiwilie−γi∏z∈wilihiz|xλiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi.Finally, (11) becomes(16)αP=∏j=1Mk|k−1φ˜W1:Sjj1∑M∈P∏j=1Mk|k−1φ˜W1:Sjj1.It can be seen from (16) that the influence of the clutter measurement set has been introduced into the two-step greedy partitioning mechanism.


### 3.2. The Numerical Implementation

We define the target state asx=x˜,X, where x˜=px,py,vx,vy,ax,ayT contains the target position, target velocity, and target acceleration and X represents the target shape. Utilizing the GIW mixture from [17, 24], we can derive a numerical implementation of the proposed filter.Suppose that the distributions of the Bernoulli component at timek−1 and the newborn target at time kk are both GIW mixtures as follows:(17)pk−1|k−1ix=∑n=1Jk−1|k−1iwn,k−1|k−1iNx˜;mn,k−1|k−1i,Pn,k−1|k−1i⊗XℐX;vn,k−1|k−1i,Vn,k−1|k−1i,pb,kjx=∑n=1Jb,kjwb,n,kjNx˜;mb,n,kj,Pb,n,kj⊗XℐX;vb,n,kj,Vb,n,kj.Jk−1|k−1i and Jb,kj are the number of the GIW component. ⊗ denotes the Kronecker product. wn,k−1|k−1i and wb,n,kj are the weight of the corresponding GIW component. Nx˜;m,P is a Gaussian function with mean m and variance P. ℐX;v,V is an inverse Wishart distribution with degrees of freedom v and inverse scale matrix V.Then, the distribution of the predicted Bernoulli component is also a GIW mixture, and the legacy Bernoulli component at timek follows(18)rs,k|k−1i=rk|k−1iρsv,(19)ps,k|k−1ix=∑n=1Jk−1|k−1iwn,k−1|k−1iNx˜;mn,k|k−1i,Pn,k|k−1i⊗XℐX;vn,k|k−1i,Vn,k|k−1i.In (19), the Gaussian and the inverse Wishart parameters satisfy(20)mn,k|k−1i=Fk|k−1⊗Idmn,k−1|k−1i,(21)Pn,k|k−1i=Fk|k−1Pn,k−1|k−1iFk|k−1T+Qk|k−1,(22)vn,k|k−1i=vn,k−1|k−1ie−Ts/τ,(23)Vn,k|k−1i=Vn,k−1|k−1iVn,k|k−1i−d−1Vn,k−1|k−1i−d−1.Fk|k−1=1Ts1/2Ts00Ts00e−Ts/θ is the state transition matrix. Qk|k−1 is the process noise covariance matrix. τ is a temporal decay constant. Ts is the sampling period. θ is the maneuver correlation time. Id is the unit matrix of size d.If the distribution of the predicted Bernoulli component satisfies(24)pk|k−1ix=∑n=1Jk|k−1iwn,k|k−1iNx˜;mn,k|k−1i,Pn,k|k−1i⊗XℐX;vn,k|k−1i,Vn,k|k−1i,then the updated Bernoulli density follows(25)rk|kP,j=αPrk|k−1j∏i=1S1−pi,d+pi,de−γi1−rk|k−1j+rk|k−1j∏i=1S1−pi,d+pi,de−γi,W1:Sj=∅,αP,W1:Sj≠∅,(26)pk|kP,jx=∑n=1Jk|k−1jwn,k|k−1jNx˜;mn,k|k−1j,Pn,k|k−1j⊗XℐX;vn,k|k−1j,Vn,k|k−1j,W1:Sj=∅,∑n=1Jk|k−1jwn,k|kP,jNx˜;mn,k|kP,j,Pn,k|kP,j⊗XℐX;vn,k|kP,j,Vn,k|kP,j,W1:Sj≠∅.In (26), the weight of the GIW component follows(27)wn,k|kP,j=wn,k|k−1j∏i,li∈TW1:Sjpi,dγiwilie−γiqzwili∏z∈wili1/λiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi∑m=1Jk|k−1jwn,k|k−1j∏i,li∈TW1:Sjpi,dγiwilie−γiqzwili∏z∈wili1/λiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi,where qzwili is the likelihood ratio.The update of the GIW component is a recursive procedure with(28)Nx˜;mn,k|kP,j,Pn,k|kP,j⊗X⋅ℐX;vn,k|kP,j,Vn,k|kP,j∝Nx˜;mn,k|k−1j,Pn,k|k−1j⊗X⋅ℐX;vn,k|k−1j,Vn,k|k−1j⋅∏i,li∈TW1:Sj∏z∈wilihiz|x,where(29)Nx˜;m,P⋅ℐX;v,V⋅∏z∈wilihiz|x=qzwili⋅Nx˜;m′,P′⋅ℐX;v′,V′.In (28), the updated GIW parameters follow(30)m′=m+K⊗Id⋅Z˜kW−Hk⊗Idm,v′=v+wili,V′=V+N+ZkW,where(31)Z˜kW=1wili∑z∈wiliz,ZkW=∑z∈wiliz−Z˜kWz−Z˜kWT,K=PHkTHkPHkT+1wili−1,N=HkPHkT+1wili−1Z˜kW−Hk⊗IdmZ˜kW−Hk⊗IdmT.Here,Hk is the observation matrix at time k.


## 3.1. The Constraint of the Clutter Measurement Set

In [14], a two-step greedy partitioning mechanism is proposed, which can obtain several high-scoring measurement partitions. In this mechanism, the target measurement set corresponding to each Bernoulli component is first obtained separately, and then, the target measurement sets corresponding to different Bernoulli components are combined to obtain all possible measurement partitions. However, it can be seen from (9) that the two-step greedy partitioning mechanism ignores the influence of the clutter measurement set W1:S0. Considering the constraint of the clutter measurement set, αP can be rewritten as(11)αP=KP/∏i=1SλiZie−λi∏j=1Mk|k−1φW1:Sjj1∑M∈PKM/∏i=1SλiZie−λi∏j=1Mk|k−1φW1:Sjj1,=∏j=1Mk|k−1φW1:Sjj1⋅∏i,li∈TW1:Sjwili!/λiwili∑M∈P∏j=1Mk|k−1φW1:Sjj1⋅∏i,li∈TW1:Sjwili!/λiwili.We introduce the function(12)φ˜W1:Sjj1=φW1:Sjj1⋅∏i,li∈TW1:Sjwili!λiwili,=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxfW1:Sj|x∏i,li∈TW1:Sjwili!λiwilidx,W1:Sj≠∅,which denotes the ratio of the likelihood that W1:Sj is generated by a target to the likelihood that W1:Sj is from the clutter.According to [24], the single-target multisensor likelihood function follows(13)fW1:Sj|x=∏i,li∈TW1:Sjpi,dγiwilie−γiwili!∏z∈wilihiz|xciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi,where γi is the mean of the extended measurements, pi,d is the target detection probability by sensor i, hiz|x is the measurement likelihood of sensor i, and ciz is the clutter intensity.Hence, (12) becomes(14)φ˜W1:Sjj1=1−rk|k−1j+rk|k−1jpk|k−1j,uη,W1:Sj=∅,rk|k−1j∫uxpk|k−1jxf˜W1:Sj|xdx,W1:Sj≠∅.In (14), f˜W1:Sj|x denotes the single-target multisensor likelihood function for the clutter case and follows(15)f˜W1:Sj|x=∏i,li∈TW1:Sjpi,dγiwilie−γi∏z∈wilihiz|xλiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi.Finally, (11) becomes(16)αP=∏j=1Mk|k−1φ˜W1:Sjj1∑M∈P∏j=1Mk|k−1φ˜W1:Sjj1.It can be seen from (16) that the influence of the clutter measurement set has been introduced into the two-step greedy partitioning mechanism.


## 3.2. The Numerical Implementation

We define the target state asx=x˜,X, where x˜=px,py,vx,vy,ax,ayT contains the target position, target velocity, and target acceleration and X represents the target shape. Utilizing the GIW mixture from [17, 24], we can derive a numerical implementation of the proposed filter.Suppose that the distributions of the Bernoulli component at timek−1 and the newborn target at time kk are both GIW mixtures as follows:(17)pk−1|k−1ix=∑n=1Jk−1|k−1iwn,k−1|k−1iNx˜;mn,k−1|k−1i,Pn,k−1|k−1i⊗XℐX;vn,k−1|k−1i,Vn,k−1|k−1i,pb,kjx=∑n=1Jb,kjwb,n,kjNx˜;mb,n,kj,Pb,n,kj⊗XℐX;vb,n,kj,Vb,n,kj.Jk−1|k−1i and Jb,kj are the number of the GIW component. ⊗ denotes the Kronecker product. wn,k−1|k−1i and wb,n,kj are the weight of the corresponding GIW component. Nx˜;m,P is a Gaussian function with mean m and variance P. ℐX;v,V is an inverse Wishart distribution with degrees of freedom v and inverse scale matrix V.Then, the distribution of the predicted Bernoulli component is also a GIW mixture, and the legacy Bernoulli component at timek follows(18)rs,k|k−1i=rk|k−1iρsv,(19)ps,k|k−1ix=∑n=1Jk−1|k−1iwn,k−1|k−1iNx˜;mn,k|k−1i,Pn,k|k−1i⊗XℐX;vn,k|k−1i,Vn,k|k−1i.In (19), the Gaussian and the inverse Wishart parameters satisfy(20)mn,k|k−1i=Fk|k−1⊗Idmn,k−1|k−1i,(21)Pn,k|k−1i=Fk|k−1Pn,k−1|k−1iFk|k−1T+Qk|k−1,(22)vn,k|k−1i=vn,k−1|k−1ie−Ts/τ,(23)Vn,k|k−1i=Vn,k−1|k−1iVn,k|k−1i−d−1Vn,k−1|k−1i−d−1.Fk|k−1=1Ts1/2Ts00Ts00e−Ts/θ is the state transition matrix. Qk|k−1 is the process noise covariance matrix. τ is a temporal decay constant. Ts is the sampling period. θ is the maneuver correlation time. Id is the unit matrix of size d.If the distribution of the predicted Bernoulli component satisfies(24)pk|k−1ix=∑n=1Jk|k−1iwn,k|k−1iNx˜;mn,k|k−1i,Pn,k|k−1i⊗XℐX;vn,k|k−1i,Vn,k|k−1i,then the updated Bernoulli density follows(25)rk|kP,j=αPrk|k−1j∏i=1S1−pi,d+pi,de−γi1−rk|k−1j+rk|k−1j∏i=1S1−pi,d+pi,de−γi,W1:Sj=∅,αP,W1:Sj≠∅,(26)pk|kP,jx=∑n=1Jk|k−1jwn,k|k−1jNx˜;mn,k|k−1j,Pn,k|k−1j⊗XℐX;vn,k|k−1j,Vn,k|k−1j,W1:Sj=∅,∑n=1Jk|k−1jwn,k|kP,jNx˜;mn,k|kP,j,Pn,k|kP,j⊗XℐX;vn,k|kP,j,Vn,k|kP,j,W1:Sj≠∅.In (26), the weight of the GIW component follows(27)wn,k|kP,j=wn,k|k−1j∏i,li∈TW1:Sjpi,dγiwilie−γiqzwili∏z∈wili1/λiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi∑m=1Jk|k−1jwn,k|k−1j∏i,li∈TW1:Sjpi,dγiwilie−γiqzwili∏z∈wili1/λiciz×∏i,∗∉TW1:Sj1−pi,d+pi,de−γi,where qzwili is the likelihood ratio.The update of the GIW component is a recursive procedure with(28)Nx˜;mn,k|kP,j,Pn,k|kP,j⊗X⋅ℐX;vn,k|kP,j,Vn,k|kP,j∝Nx˜;mn,k|k−1j,Pn,k|k−1j⊗X⋅ℐX;vn,k|k−1j,Vn,k|k−1j⋅∏i,li∈TW1:Sj∏z∈wilihiz|x,where(29)Nx˜;m,P⋅ℐX;v,V⋅∏z∈wilihiz|x=qzwili⋅Nx˜;m′,P′⋅ℐX;v′,V′.In (28), the updated GIW parameters follow(30)m′=m+K⊗Id⋅Z˜kW−Hk⊗Idm,v′=v+wili,V′=V+N+ZkW,where(31)Z˜kW=1wili∑z∈wiliz,ZkW=∑z∈wiliz−Z˜kWz−Z˜kWT,K=PHkTHkPHkT+1wili−1,N=HkPHkT+1wili−1Z˜kW−Hk⊗IdmZ˜kW−Hk⊗IdmT.Here,Hk is the observation matrix at time k.


## 4. Experimental Results

In this section, we consider tracking 4 extended targets moving within an observation area of size [−1000 m, 1000 m]× [−1000 m, 1000 m]. The target shape is modeled as an ellipse [17] with a major axis A=20 and a minor axis a=5. The duration time is 80 s. In the simulation, the constant acceleration (CA) model is considered. The target trajectories are shown in Figure 1. The initial position and the survival time of all the targets are shown in Table 1.Figure 1
Target trajectories.Table 1
The true movement of targets.
The initial position of targetsSurvival time[200 m, −800 m, 0 m/s, 35 m/s, 1 m/s2, −1 m/s2](5 s∼40 s)[−200m, −800 m, 0 m/s, 35 m/s, −1 m/s2, −1 m/s2](45 s∼80 s)[−800 m, 200 m, 40 m/s, 20 m/s, 0 m/s2, 0 m/s2](5 s∼40 s)[800 m, 800 m, −25 m/s, −25 m/s, 0 m/s2, 0 m/s2](45 s∼80 s)In (22), the sampling period Ts=1s and the temporal decay constant τ=5s . The maneuver correlation time θ=0.5. In the two-step measurement partitioning method [14], the number of subsets in P is no more than Wmax=6, and the number of partitions in P is no more than Pmax=6. The pruning threshold is set as Tmax=10−5, and the number of the GIW components of each Bernoulli density is up to Jmax=100. The merge threshold of Gaussian components is UG=50, and the merge threshold of inverse Wishart components is UIW=50. We use the optimal subpattern assignment (OSPA) distance [26] as the metric to evaluate the filtering performance. The OSPA distance accounts for estimation errors in both position and cardinality, and in this section, we set the order p=1 and the penalty factor c=100.Now, we consider observations from 4 sensors. The sensors are located at [−900 m, 0 m], [0 m, −900 m], [900 m, 0 m], and [0 m, 900 m], respectively. For every sensor, the detection probabilitypi,d=0.7 and the clutter density parameter λi=80. The measurement noise standard deviations along x-y coordinates of each sensor are εx1=εy1=5m, εx2=εy2=5m, εx3=εy3=10m, and εx4=εy4=10m. The mean of the extended measurements is set as γi=20. Figure 2 shows the measurements collected from sensor 1. The estimated states by using the proposed method and the ET-MS-MeMBer filter are shown in Figures 3 and 4, respectively. We can observe that compared with the ET-MS-MeMBer filter, the proposed method can obtain more accurate state estimation results. We also compare the filtering performance of the ET-MS-MeMBer filter with that of the proposed method via 50 Monte Carlo runs. Figure 5 shows the cardinality estimations of the two methods. The black line is the true cardinality of the targets, and the magenta dotted line and the blue dotted line represent estimated standard deviations of the ET-MS-MeMBer filter and the proposed filter, respectively. We can observe that compared with the ET-MS-MeMBer filter, the proposed method obtains a smaller standard deviation of the estimated cardinality. Figure 6 shows the OSPA errors of the two methods. Note that, for both two methods, the delay in response to target birth and target death leads to the high peaks at the corresponding times. Nevertheless, the OSPA error of the proposed method is still smaller than that of the ET-MS-MeMBer filter.Figure 2
The measurements of sensor 1. The black dot represents the target measurement, and the gray dot represents the clutter.Figure 3
The estimated states of the proposed method.Figure 4
The estimated states of the ET-MS-MeMBer filter.Figure 5
Comparison of the estimated cardinality.Figure 6
Comparison of the OSPA error.We also consider the case of a varying number of sensors withpi,d=0.7 and λi=80 for each sensor i. As before, 50 Monte Carlo runs are performed. Figure 7 shows the time-averaged OSPA error of the two methods, and Figure 8 shows the running time. We observe that as the number of the sensors increases, the time-average OSPA error of both two methods decreases, however, at the cost of the growing computational time. Comparatively, the proposed method performs better.Figure 7
Time-average OSPA errors with different sensor numbers.Figure 8
Running time with different sensor numbers.Figure9 shows the time-averaged OSPA error of the two methods in scenarios with different clutter intensities. The corresponding running time of the two methods is shown in Figure 10. The clutter density λi is varied in a range of {5, 20, 40, 60, 80, 100}. Besides, the number of sensors is 3, and the detection probability is pi,d=0.7. It is observed that the OSPA error and running time of the two methods are similar in the low clutter density case, i.e., λi=5. With the increase of the clutter intensity, the filtering performance of the two methods decreases. However, the proposed method requires less time and has smaller OSPA errors when the clutter intensity is high.Figure 9
Time-average OSPA errors with different clutter intensities.Figure 10
Running time with different clutter intensities.


## 5. Conclusions

In this paper, we propose a robust MS-MeMBer filter for extended-target tracking in clutters. We modify the multisensor measurement partitioning process by considering the constraint of the clutter measurement set to make it robust under conditions with different clutter intensities. Accordingly, we propose the single-target multisensor likelihood function for the clutter case. Simulation results show that compared with the ET-MS-MeMBer filter, the proposed method has smaller OSPA errors with less time consumption under high-clutter-density conditions.



---

*Source: 9942365-2021-05-20.xml*