Next Article in Journal
Methodology for Addressing Infectious Aerosol Persistence in Real-Time Using Sensor Network
Previous Article in Journal
Unified Chassis Control of Electric Vehicles Considering Wheel Vertical Vibrations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Labeled GM-PHD Filter for Explicitly Tracking Multiple Targets †

1
College of Energy and Electrical Engineering, Hohai University, Nanjing 210098, China
2
Laboratory of Array and Information Processing, College of Computer and Information, Hohai University, Nanjing 210098, China
*
Author to whom correspondence should be addressed.
This paper is an extension of our conference paper Gao, Y.; Fang, N.; Jiang, D.; Fu, W.; Guo, S. An explicit track continuity algorithm for the GM-PHD filter. In Proceedings of the 31th Chinese Control and Decision Conference, Nanchang, China, 3–5 June 2019.
Submission received: 24 April 2021 / Revised: 2 June 2021 / Accepted: 3 June 2021 / Published: 7 June 2021
(This article belongs to the Section Physical Sensors)

Abstract

:
In this study, an explicit track continuity algorithm is proposed for multitarget tracking (MTT) based on the Gaussian mixture (GM) implementation of the probability hypothesis density (PHD) filter. Trajectory maintenance and multitarget state extraction in the GM-PHD filter have not been effectively integrated to date. To address this problem, we propose an improved GM-PHD filter. In this approach, the Gaussian components are classified and labeled, and multitarget state extraction is converted into multiple single-state extractions. This provides the identity label of the individual target and can shield against the negative effects of clutter in the prior density region on the estimates, thus realizing the integration of trajectory maintenance with state extraction in the GM-PHD filter. As no additional associated procedures are required, the overall real-time performance of the proposed filter is similar to or slightly lower than that of the basic GM-PHD filter. The results of numerical experiments demonstrate that the proposed approach can achieve explicit track continuity.

1. Introduction

1.1. Background and Motivation

Multitarget tracking (MTT) jointly estimates the target state and the number of targets simultaneously from a sequence of measurements [1]. In MTT, the time-varying number of targets causes data associations between state and measurement sets of targets, which are uncertain. The key challenges encountered involve detection uncertainty, clutter, and data association uncertainty. Therefore, accurately extracting the states of targets and their trajectories has become a critical issue and a hot topic in the field of MTT. The most popular technologies for MTT are the joint probabilistic data association (JPDA) [2,3], multiple hypothesis tracking (MHT) [4], and random finite set (RFS) [1,5].
Since RFS avoids the traditional data association, it has attracted significant attention. Based on RFS and under the framework of Bayesian filtering, many approximations of Bayes multitarget filters have been proposed, mainly including a probability hypothesis density (PHD) filter [1,6,7], cardinalized PHD (CPHD) filter [8], and multi-Bernoulli (MeMBer) filter [9]. These filters are not multitarget trackers, which only estimate target states at individual time instants as opposed to multitarget trajectories. While these filters were not designed with the aim of estimating the trajectories of targets [10], they have been used in many applications [11,12,13].
To provide the information of each track, the notion of labeled RFS is introduced in [14,15]. The generalized labeled multi-Bernoulli (GLMB) filter [14,15] leads to an analytic solution to the Bayes multitarget tracker, which significantly improves the accuracy of multitarget state extraction and explicitly accommodates the estimation of target trajectories. Recently, labeled RFS-based filters and smoothers [10,16,17,18,19] have been further developed, for generating track estimates, which is also the focus of this paper.
A Gibbs-GLMB filter [10] is a computational efficient implementation of the GLMB filter that combines the prediction and update into a single step. A forward-backward labeled multi-Bernoulli smoother in [17] is proposed for MTT, which can improve both the cardinality estimation and the state estimation, and the major computational complexity is linear with the number of targets. The trajectory Poisson multi-Bernoulli filter and the trajectory Poisson multi-Bernoulli mixture filter [18] have better filtering accuracy and real-time performance than those of the Gibbs-GLMB filter. A onetime step lagged Bayes multitarget smoother using the δ - GLMB density [19], which also inherently outputs targets trajectories, outperforms the Gibbs-GLMB filter on both the estimates and target number and state at the cost of higher computational complexity. However, when targets are in a dense clutter, or misdetected, generating the correct multitarget trajectories is difficult in [16,17,18,19], although the computational complexity problem has been improved in [17,18].
Thus, it is necessary to design a relatively low computing cost filter that can provide explicit trajectories whilst considering the intense clutter environment.

1.2. Brief Survey of Related Work

From the viewpoint of practical applicability, the PHD filter is suitable for applications demanding real-time results [20,21,22]. Unlike the full Bayesian filter, the PHD filter iteratively propagates the first moment of the multitarget posterior density. It is a promising technique in terms of computational complexity to solve the MTT problem. However, it does not explicitly accommodate the trajectories of targets [10,21,23]. Furthermore, the PHD filter is known to produce highly uncertain estimates of the target number [21,23]. In recent years, several modifications have been proposed for the trajectories of targets [22,24,25,26,27,28,29,30].
In [29], an explicit track continuity algorithm for the sequence Monte Carlo PHD (SMC-PHD) filter [31] is proposed, which not only shields against the negative effects of clutter but provides the identities of the individual targets. However, when several newborn targets appear simultaneously in one newborn region, it is unable to extract the state estimates of different observations from these newborn targets. Compared to the SMC-PHD filter, the Gaussian mixture PHD (GM-PHD) filter [32] has the advantages of simple state extraction and low computational cost, which is more suitable for the requirement of real-time scenes. Thus, the focus of this study is the GM-PHD filter.
For the estimation of the target trajectory of the GM-PHD filter, recently proposed solutions involve integrating the track identity into the GM-PHD filter. In each iteration process, each Gaussian component is assigned a label to indicate its correlation with one track. These methods can be mainly classified into two groups. One group [24,25] involves running an additional association method based on the estimates outputted by the filter. The other solution [26,27,28,29] extends a label to the state estimate. However, when targets are in close proximity to each other, in a dense clutter, or misdetected, maintaining the correct multitarget trajectories is difficult, although the accuracy of the estimates of the target number has been increased significantly by improving multitarget state extraction strategies [33,34,35,36,37]. Thus, one can conclude that in the GM-PHD filter, track maintenance and multitarget state extraction have not been effectively integrated. Recently, a Gaussian mixture trajectory PHD filter was proposed in [30], which is able to provide trajectory estimates for live targets, without adding labels or tags. It adapts the estimator for the GMPHD filter for the sets of trajectories; thus, its tracking accuracy is the same as that of the PHD filter. When clutters appear in the high prior density region of the target, it is difficult to shield against the interference of clutters.
An instructive idea is that the Gaussian components approximating the posterior information of one target have the same identity label, which will not be changed throughout the overall filtering process. In the multitarget state extraction, only one state estimate can be extracted from the updated survival Gaussian components with the same label. Then, the explicit trajectories can be obtained by simply connecting the state estimates with the same label, which requires no additional track management. To achieve this concept, in this paper, we propose an explicit track continuity algorithm for the GM-PHD filter.

1.3. Main Contributions

The work in this paper is based on a published conference paper [38]. The major contributions of this paper are as follows:
  • First, each Gaussian component is assigned a label, and the identity labels of Gaussian components are divided into three classes; once it is determined as a confirmed component, its label will be unchanged throughout the whole filtering process. This is convenient for track maintenance; that is, state estimates with the same label at different times belong to the same track.
  • Second, based on the measurements selected by eliminating most clutter, the updated components are obtained, and their parameters are stored in four different query tables. Based only on these query tables, the multitarget state extraction is converted to multiple single-target state extractions providing the identity label, and the negative effects of clutter in the prior density region on the estimates can be prevented. Thus, we realize the integration of trajectory maintenance with state extraction in the GM-PHD filter.
  • Third, since the component updating process is based on the selected measurements and no additional associated process is required, the overall real-time performance of the proposed filter is similar to or slightly lower than that of the basic GM-PHD filter. The performance of the proposed approach is verified by simulations designed in a linear scenario.
The remainder of this paper is organized as follows. The technical background is provided in Section 2, while the specific designing of the improved GM-PHD filter is presented in Section 3. The quantitative experiments are described in Section 4, and concluding remarks are provided in Section 5.

2. Background

2.1. GM-PHD Filter

In MTT, the state and observation sets are time-varying and unknown. The collection of the states of targets at time k is represented as an RFS X k = { x k ( 1 ) , , x k ( i ) , , x k ( | X k | ) } , where | X k | denotes the number of targets, and x k ( i ) denotes the state of an individual target. Similarly, the collection of available measurements is an RFS Z k = { z k ( 1 ) , , z k ( i ) , , z k ( | Z k | ) } , where | Z k | is the number of measurements and z k ( i ) is a measurement from a target or due to a clutter. The PHD [6] is the first-order moment of multitarget posterior density.
Assuming that each target follows a linear Gaussian dynamical model and the sensor has a linear Gaussian measurement model [32],
f k | k 1 ( x k | x k 1 ) = N ( x k ; F k 1 x k 1 , Q k 1 )
g k ( z k | x k ) = N ( z k ; H k x k , R k )
where N ( ; m , P ) denotes the Gaussian density with mean m and covariance P , F k 1 is the state transition matrix, Q k 1 is the process noise covariance, H k is the observation matrix, and R k is the measurement noise covariance.
The birth intensity function of the new targets at time k [32] is γ k ( x ) = i = 1 J γ , k w γ , k ( i ) N ( x ; m γ , k ( i ) , P γ , k ( i ) ) , where J γ , k denotes the number of newborn Gaussian components, and w γ , k ( i ) , m γ , k ( i ) , and P γ , k ( i ) are the weight, mean, and covariance, respectively, of the i th newborn Gaussian component. Suppose that the survival probability p S , k = p S , k ( x k ) , the detection probability p D , k = p D , k ( x k ) . Based on these assumptions [32], the PHD filter can be efficiently approximated by mixed Gaussian components. w is the weight of Gaussian density and the posterior intensity D k 1 | k 1 ( x ) at time k 1 [32] is a Gaussian mixture form with J k 1 survival components as follows:
D k 1 | k 1 ( x ) = i = 1 J k 1 w k 1 ( i ) N ( x ; m k 1 ( i ) , P k 1 ( i ) )
The predicted intensity D k | k 1 ( x ) at time k is also a Gaussian mixture with J k | k 1 components calculated as follows:
D k | k 1 ( x ) = i = 1 J k | k 1 w k | k 1 ( i ) N ( x ; m k | k 1 ( i ) , P k | k 1 ( i ) ) = D S , k | k 1 ( x ) + γ k ( x )
where D S , k | k 1 ( x ) = i = 1 J k 1 w S , k | k 1 ( i ) N ( x ; m S , k | k 1 ( i ) , P S , k | k 1 ( i ) ) is the intensity of survival targets, w S , k | k 1 ( i ) = p S , k w k 1 ( i ) , m S , k | k 1 ( i ) = F k 1 m k 1 ( i ) , P S , k | k 1 ( i ) = F k 1 P k 1 ( i ) ( F k 1 ) T + Q k 1 .
When a new set Z k arrives, the posterior intensity D k | k at time k [32] is updated and can be described as follows:
D k | k ( x ) = ( 1 p D , k ) D k | k 1 ( x ) + z Z k D D , k ( x ; z )
where
D D , k ( x ; z ) = i = 1 J k | k 1 w k | k ( i ) ( z ) N ( x ; m k | k ( i ) ( z ) , P k | k ( i ) )
w k | k ( i ) ( z ) = p D , k w k | k 1 ( i ) g k ( z | m k | k 1 ( i ) ) λ c ( z ) + p D , k l = 1 J k | k 1 w k | k 1 ( l ) g k ( z | m k | k 1 ( l ) )
g k ( z | m k | k 1 ( i ) ) = N ( z ; H k m k | k 1 ( i ) , S k ( i ) )
m k | k ( i ) ( z ) = m k | k 1 ( i ) + K k ( i ) ( z H k m k | k 1 ( i ) )
P k | k ( i ) = [ I K k ( i ) H k ] P k | k 1 ( i )
S k ( i ) = R k + H k P k | k 1 ( i ) ( H k ) T
K k ( i ) = P k | k 1 ( i ) ( H k ) T ( H k P k | k 1 ( i ) ( H k ) T + R k ) 1
For more details on the GM-PHD filter, we refer the reader to the original study [32].
To reduce the computational load, the increasing Gaussian components must be pruned and merged. Then, the posterior intensity D k | k is represented as follows:
D k | k ( x ) = i = 1 J k w k ( i ) N ( x ; m k ( i ) , P k ( i ) )
The multiple-target state estimates can be extracted based on the threshold w t h (this is generally 0.5) from D k | k ( x ) . Figure 1 shows the overall procedure of the basic GM-PHD filter.

2.2. Problem Formulation

The basic GM-PHD filter in Figure 1 does not explicitly accommodate target trajectory estimation. Trajectory continuity and state estimates essentially influence and interact with each other. In MTT, detection uncertainty, false state estimates, and closely spaced targets have detrimental effects on track continuity. Figure 2 graphically depicts these phenomena.
As shown in Figure 2, the clutter rate was relatively high, and the detection probability was 0.9 in the current tracking scenario. In Figure 2a, at times k = 7 and k = 8, the target was misdetected, and its posterior intensity was weakened. This caused its inability to obtain timely state estimates even if the target was redetected at times k = 9 and k = 10. Managing the trajectories of the state estimates at times k = 6 and k = 11 was a challenge. In Figure 2b, at times k = 60, k = 61, k = 62, and k = 65, there were false state estimates caused by clutter points appearing in the high-prior density region of the target. In Figure 2c, two targets were closely spaced, and their state estimates were nearby. It is difficult to identify the state estimate with which they should be associated. The above phenomena make it difficult to provide a correct target trajectory based on the basic GM-PHD filter.
A trajectory is a time-sequence of state estimates with the same label, such as { ( x ^ k 1 , ) , ( x ^ k , ) , ( x ^ k + 2 , ) } , wherein provides the means to identify the trajectory or track of one target. Furthermore, even if there is no state estimate at time k + 1, the correct track continuity can still be obtained. To achieve explicit trajectories based on the GM-PHD filter, we proposed an algorithm that labels all Gaussian components, which we called the labeled GM-PHD (LGM-PHD) filter.

3. The Proposed Multitarget Tracking Algorithm

Assuming the target birth intensity is known a priori, a “one-to-one” principle for the multitarget state extraction of the GM-PHD filter was proposed herein.
The “one-to-one” principle for the GM-PHD filter: In the multitarget state extraction, only one state estimate can be extracted from the updated Gaussian components with the same label, originating from the survival ones.
To abide by the “one-to-one” principle, each Gaussian component is assigned a label, and these labels are divided into three categories in the LGM-PHD filter. In this paper, the identity labels of the Gaussian components are divided into three classes, the labels of newborn Gaussian components, the labels of unconfirmed Gaussian components (approximating the posterior information of targets having their own unconfirmed trajectory), and the labels of confirmed components (approximating the posterior information of targets having their own confirmed trajectory). The attribute of a label is determined based on the value range as follows:
label   of   a   newborn   component , if   V n e w , label   of   an   unconfirmed   component , if   V u n < V n e w , label   of   a   confirmed   component , if   0 < V u n .
where V u n is a value far greater than the total number of targets that could possibly appear in this scenario, V n e w V u n , and V n e w V u n . In general, V u n = 200 and V n e w = 400 .
r max and r max t are used to count the confirmed and unconfirmed trajectories, respectively. At time k = 0, r max is initialized as 0 and r max t is initialized as V u n . When a new confirmed trajectory is initiated, r max is incremented by one; when a new unconfirmed trajectory is initiated, r max t is also incremented by one.
The overall procedure of the LGM-PHD filter contains seven steps—initialization, prediction, observation selection, updating, state extraction, pruning, and merging and capping—as shown in Figure 3. When compared with the basic GM-PHD filter in Figure 1, the submodules marked with “ ” are either modified or proposed in this paper.

3.1. Initialization and Prediction

At time k , newborn targets appear spontaneously according to a Poisson point process with the intensity function
γ k ( x ) = i = 1 J γ , k w γ , k ( i ) N ( x ; m γ , k ( i ) , P γ , k ( i ) )
where each newborn Gaussian component is assigned a label γ , k ( i ) = V n e w + i , i = 1 , , J γ , k . According to (15), the initial intensity D 0 | 0 is described as follows: { w γ , 0 ( i ) , m γ , 0 ( i ) , P γ , 0 ( i ) , γ , 0 ( i ) } i = 1 J γ , 0 .
With this assumption, at time k − 1, the posterior probability density D k 1 | k 1 ( x ) is approximated by { w k 1 ( j ) , m k 1 ( j ) , P k 1 ( j ) , k 1 ( j ) } j = 1 J k 1 , where k 1 ( j ) is the identity label of the j th survival Gaussian component. Then, according to (4), the priori density at time k is approximated by { w k | k 1 ( i ) , m k | k 1 ( i ) , P k | k 1 ( i ) , k | k 1 ( i ) } i = 1 J k | k 1 , where J k | k 1 = J k 1 + J γ , k . The i th predicted Gaussian component is labeled as follows:
k | k 1 ( i ) = { k 1 ( i ) , i = 1 , , J k 1 γ , k ( i J k 1 ) , i = J k 1 + 1 , , J k | k 1
To obtain the time-sequence of states with the same label, the labels of conformed components remain unchanged throughout the overall filtering process. In addition, the labels of all types of components can be inherited, and only components with the same label can be merged.

3.2. Observation Selection

To achieve robust track maintenance, clutter should be eliminated. Considering the overall real-time performance of the filter, it is necessary to eliminate most clutter before the updating procedure, to select effective observations. Given the zero-mean-measurement Gaussian white noise with the covariance matrix R k = d i a g ( [ σ x 2 σ y 2 ] T ) , the threshold d ( a ) is calculated by d ( a ) = a [ σ x σ y ] T [ 0 0 ] T , in which is the Euclidean distance and a is the confidence coefficient. In this gating algorithm, a = 6 . Initializing the effective observation set Z k , e f as , based on all the predicted Gaussian components at time k , Z k , e f = { z k , e f ( j ) } j = 1 | Z k , e f | is obtained by the operations in Algorithm 1.
Algorithm 1. The observation selection for the GM-PHD filter algorithm
1: Given { w k | k 1 ( i ) , m k | k 1 ( i ) , P k | k 1 ( i ) } i = 1 J k | k 1 and Z k = { z k ( j ) } j = 1 | Z k |
2: Initialization Z k , e f =
3: for j = 1 , , | Z k | do
4:  for i = 1 , , J k | k 1 do
5:       d ( j , i ) = z k ( j ) H k m k | k 1 ( i )
6:      if d ( j , i ) d ( a )
7:       Z k , e f = [ Z k , e f , z k ( j ) ]
8:      break
9:     end
10:   end
11:  end

3.3. Updating

First, four tables U w , U m , U P , and U , all with size J k | k 1 × ( 1 + | Z k , e f | ) , are created to store parameters w , m , P , and of each updated Gaussian component. They are uniformly called query table U .
Second, based on Z k , e f , the posterior intensity D k | k at time k can be described as follows: D k | k ( x ) = ( 1 p D , k ) D k | k 1 ( x ) + z k , e f ( j ) Z k , e f D D , k ( x ; z k , e f ( j ) ) , where D D , k ( x ; z k , e f ( j ) ) = i = 1 J k | k 1 w k | k ( i ) ( z k , e f ( j ) ) N ( x ; m k | k ( i ) ( z k , e f ( j ) ) , P k | k ( i ) ) . Based on the i th predicted Gaussian component and z k , e f ( j ) , we obtain the parameter set { w k | k ( i ) ( z k , e f ( j ) ) , m k | k ( i ) ( z k , e f ( j ) ) , P k | k ( i ) , k | k ( i ) } of the updated Gaussian component via (5)–(12), where k | k ( i ) = k | k 1 ( i ) . Based on this parameter set, the elements in the i th row and the ( j + 1 ) th column of U w , U m , U P , and U are obtained as follows:
U w ( i , j + 1 ) = w k | k ( i ) ( z k , e f ( j ) ) U m ( i , j + 1 ) = m k | k ( i ) ( z k , e f ( j ) ) U P ( i , j + 1 ) = P k | k ( i ) U ( i , j + 1 ) = k | k ( i )
( 1 p D , k ) D k | k 1 ( x ) is approximated by the parameter set { ( 1 p D , k ) w k | k 1 ( i ) , m k | k 1 ( i ) , P k | k 1 ( i ) , k | k 1 ( i ) } i = 1 J k | k 1 . Based on this parameter set, the elements in the i th row and the 1 th column of U w , U m , U P , and U are obtained as follows:
U w ( i , 1 ) = ( 1 p D , k ) w k | k 1 ( i ) U m ( i , 1 ) = m k | k 1 ( i ) U P ( i , 1 ) = P k | k 1 ( i ) U ( i , 1 ) = k | k 1 ( i )
Third, four tables, U w , U m , U P , and U , are created, which are collectively called the U
U w = U w , U m = U m , U P = U P , U = U
When a new confirmed trajectory or a new unconfirmed trajectory is established, its information is documented and stored in U .

3.4. State Extraction

In this section, by means of U , certain Gaussian components with the same label are associated with one effective observation, and if the maximum weight in these components satisfies some conditions, one state estimate is extracted.
Complete target status information contains state estimate x ^ , track label , and occurrence time k . Thus, at time k , the state estimate set X ^ k , the track label set k , s t , and the occurrence time set T k , s t are all initialized as . The specific steps in state extraction process are as follows.
Step 1: First, obtain the maximum value w ¯ in U w ( : , 2 : e n d ) , where U w ( : , 2 : e n d ) represents the elements from the second column to the last column of each row in U w . U w ( i , j ) denotes the element in the j th column of the i th row in U w . Here, w ¯ = U w ( i , j ) is the maximum, ( i , j ) = arg max ( U w ( i , j ) ) , i = 1 , , J k | k 1 , and j = 2 , , 1 + | Z k , e f | . Very rarely, w ¯ may be in two or more different locations in U w .
Next, if w ¯ < w b , jump to Step 4. w b is the threshold of the basic weight. In general, w b is set to 0.02. If w ¯ w b , based on the row and column numbers of w ¯ , obtain the unique label set w ¯ = { w ¯ ( i ) = U ( r n ( i ) , c n ( i ) ) } i = 1 | w ¯ | , and initialize C N = . C N is used to store the column numbers of observations based on which some tracks will be established or maintained. In most cases, the label set only has one element, i.e., | w ¯ | = 1 .
Step 2: For i = 1 , , | w ¯ | , determine the attribute of w ¯ ( i ) , and execute the following actions for the Gaussian components with w ¯ ( i ) .
  • First, obtain x ^ by x ^ = H k U m ( r n ( i ) , c n ( i ) ) . The row numbers of elements with w ¯ ( i ) in U are denoted as R N w ¯ ( i ) .
  • Next, the attribute of w ¯ ( i ) will be used to determine which of the following three cases will be executed.
Case A: w ¯ ( i ) V n e w , the label of newborn components.
As several newborn targets may appear simultaneously in one newborn region, newborn Gaussian components with the same label may contain information on several newborn targets. Therefore, several components with large weights in U w ( R N w ¯ ( i ) , 2 : ( 1 + | Z k , e f | ) ) can be extracted as the state estimates of different observations, and new tracks can be established accordingly. The specific steps are as follows:
① The sequence number set C N w l of the different observations with large weights in U w ( R N w ¯ ( i ) , 2 : ( 1 + | Z k , e f | ) ) is obtained as follows:
C N w l = { j | U w ( R N w ¯ ( i ) , j ) w l , j = 2 , , 1 + | Z k , e f | } C N w l = { c n w l ( j 1 ) } j 1 = 1 | C N w l | = u n i ( C N w l )
where w l is the threshold of large weights and u n i ( A ) is an operator that obtains the unique elements in set A . In general, w l is set to 0.4.
Additionally, C N = [ C N , C N w l ] .
② If | C N w l | 1 , for j 1 = 1 , , | C N w l | , | C N w l | state estimates will be extracted for the large weight observations and | C N w l | new confirmed tracks will be established. The j 1 th target status information is obtained as follows:
r max = r max + 1 , r n w l = arg max ( i 1 R N w ¯ ( i ) U w ( i 1 , c n w l ( j 1 ) ) ) , X ^ k = [ X ^ k , H k U m ( r n w l , c n w l ( j 1 ) ) ] , k , s t = [ k , s t , r max ] , T k , s t = [ T k , s t , k ] .
In (21), the label of the new confirmed track is given by r max = r max + 1 , r n w l is the row number of the maximum weight in U w ( R N w ¯ ( i ) , c n w l ( j 1 ) ) and the single-target state estimate is H k U m ( r n w l , c n w l ( j 1 ) ) , which is stored in state estimate set X ^ k and labeled as r max . The information on the j 1 th new confirmed track is established by modifying U w , U m , U P , and U as follows:
U w ( R N w ¯ ( i ) , : ) = 0 , U w ( r n w l , c n w l ( j 1 ) ) = U w ( r n w l , c n w l ( j 1 ) ) , U ( R N w ¯ ( i ) , : ) = r max . } if j 1 = 1 . U w ( end + 1 , : ) = 0 , U w ( end , c n w l ( j 1 ) ) = U w ( r n w l , c n w l ( j 1 ) ) , U ( end + 1 , : ) = r max , U m ( end + 1 , : ) = U m ( r n w l , : ) , U P ( end + 1 , : ) = U P ( r n w l , : ) . } if j 1 > 1 .
where U ( end , j ) denotes the elements in the j th column of the last row in U , and U ( end + 1 , : ) denotes adding one row at the bottom of U . If j 1 > 1 , it means at least one state estimate has been extracted at time k . In other words, two or more confirmed trajectories will be simultaneously established from one prior newborn region. Thus, one additional row should be added in U to establish the posterior information of the individual track. To shield against interference from closely spaced components with another label, only the elements in U ( r n w l , c n w l ( j 1 ) ) can be reserved.
③ Additionally, some components whose weights are between w s and w l in U w ( R N w ¯ ( i ) , 2 : ( 1 + | Z k , e f | ) ) should also be associated with different observations. The sequence number set C N w b l of these observations with small weights is given as follows:
C N w b l = { j | w s U w ( R N w ¯ ( i ) , j ) < w l , j = 2 , , 1 + | Z k , e f | } C N w b l = u n i ( C N w b l ) C N w b l = { c n w b l ( j 2 ) } j 2 = 1 | C N w b l | = C N w b l C N w l
where w s is the threshold of small weights. In general, w s is set to 0.1.
Additionally, C N = [ C N , C N w b l ] .
④ If | C N w b l | 1 , for j 2 = 1 , , | C N w b l | , | C N w b l | unconfirmed tracks will be established for the small weight observations. The information on the j 2 th unconfirmed track is obtained by specific steps as follows:
r max t = r max t + 1
r n w b l = arg max ( i 1 R N w ¯ ( i ) U w ( i 1 , c n w b l ( j 2 ) ) )  
U w ( R N w ¯ ( i ) , : ) = 0 , U w ( r n w b l , c n w b l ( j 2 ) ) = U w ( r n w b l , c n w b l ( j 2 ) ) , U ( R N w ¯ ( i ) , : ) = r max t . } if | C N w l | = 0 & j 2 = 1 . U w ( end + 1 , : ) = 0 , U w ( end , c n w b l ( j 2 ) ) = U w ( r n w b l , c n w b l ( j 2 ) ) , U ( end + 1 , : ) = r max t , U m ( end + 1 , : ) = U m ( r n w b l , : ) , U P ( end + 1 , : ) = U w ( r n w b l , : ) . } if ( | C N w l | = 0 & j 2 > 1 ) | ( | C N w l | 1 ) .
where U ( end , j ) and U ( end + 1 , : ) are the same as those in (22).
Note that no state estimate is extracted here, since state estimates are only extracted from newborn components with large weight or from confirmed components. Storing the parameter set in U makes it easy to determine the observation associated with one updated component with the maximum weight, such as C N w l and C N w b l .
Case B: V u n w ¯ ( i ) < V n e w , the label of unconfirmed components.
① If w ¯ w l or ( w m w ¯ < w l ) & ( w ¯ ( i ) k 1 , u n ) & ( w ¯ ( i ) k 2 , u n ) , new target status information is obtained as follows:
r max = r max + 1 , X ^ k = [ X ^ k , H k U m ( r n ( i ) , c n ( i ) ) ] , k , s t = [ k , s t , r max ] , T k , s t = [ T k , s t , k ] .
where k 1 , u n is the label set of unconfirmed components with small weights at time k 1 , w b < w m < w l . In general, w m is set to 0.2.
The information on this new confirmed track is established by modifying U w and U as follows:
U w ( R N w ¯ ( i ) , : ) = 0 , U w ( R N w ¯ ( i ) , c n ( i ) ) = U w ( R N w ¯ ( i ) , c n ( i ) ) , U ( R N w ¯ ( i ) , : ) = r max .
Additionally, C N = [ C N , c n ( i ) ] .
② Otherwise, w ¯ ( i ) is stored in k , u n , i.e., k , u n = [ k , u n , w ¯ ( i ) ] .
Case C: 0 w ¯ ( i ) < V u n , the label of confirmed components. The single-target state estimate is extracted.
① Preprocessing state estimate. If the observation of one target drifts from its real location, or one target is misdetected at time k 1 , the distribution of this target’s posterior information will diverge, and the weights of components associated with its observation at time k will be weakened. Thus, if w ¯ < w m , r n w m = arg max ( i 1 R N w ¯ ( i ) U w ( i 1 , 1 ) ) . Furthermore, if H k U m ( r n ( i ) , c n ( i ) ) H k U m ( r n w m , 1 ) > d ( 1 ) , x ^ = H k U m ( r n w m , 1 ) . The definition of d ( ) is provided in Section 3.2.
Remark 1.
As shown inFigure A1, the 1 st column elements in U represent the updated information for misdetection at time k . When there is some deviation, i.e., w ¯ < w m , H k U m ( r n w m , 1 ) is the optimal state estimate since U m ( r n w m , 1 ) is the update of the posterior component with w ¯ ( i ) at time k 1 , and U w ( r n w m , 1 ) is the maximum weight of Gaussian components with w ¯ ( i ) .
② Survival target status information is obtained as follows:
X ^ k = [ X ^ k , x ^ ] , k , s t = [ k , s t , w ¯ ( i ) ] , T k , s t = [ T k , s t , k ] .
Additionally, C N = [ C N , c n ( i ) ] .
③ Modifying the weights of the Gaussian components. The larger w ¯ is, the more concentrated the posterior distribution is, and the smaller the radius of the effective region is. The smaller w ¯ is, the more dispersed the posterior distribution is, and the larger the radius of the effective region is. By calculating the distances between the current x ^ and the components with the same label via (30), the weight is modified via (31).
d ( i 1 , j ) = x ^ H k U m ( R N w ¯ ( i ) ( i 1 ) , j ) , i 1 = 1 , , | R N w ¯ ( i ) | , j = 2 , , 1 + | Z k , e f |
U w ( R N w ¯ ( i ) ( i 1 ) , j ) = a 1 w ¯ , ( w ¯ 0.5 & d ( i 1 , j ) > d ( 3 ) ) | ( 0.4 w ¯ < 0.5 & d ( i 1 , j ) > d ( 4 ) ) | ( 0.2 w ¯ < 0.4 & d ( i 1 , j ) > d ( 5 ) ) | ( w ¯ < 0.2 & d ( i 1 , j ) > d ( 6 ) )
In the effective region, the weights of Gaussian components remain unchanged; however, they are reduced to a 1 w ¯ outside the region. Here, a 1 = 0.3 . Dividing the effective region of the posterior distribution based on the weights can further shield against interference, such as detection uncertainty and closely spaced measurements originating from clutter or the other target.
  • According to the “one-to-one” principle, the elements in U ( R N w ¯ ( i ) , : ) are employed in the state extraction; then, they will be discarded as follows:
    U w ( R N w ¯ ( i ) , : ) = 0
Step 3: Based on the components associated with effective observations, whose sequence numbers are stored in C N , some tracks are established or maintained. Thus, after the i = 1 , , | w ¯ | loop is completed, the elements in U ( : , C N ) will be discarded as follows:
U w ( : , C N ) = 0
Next, return to Step 1.
Step 4: Here, n r U denotes the number of rows in U , U w , U m , U P , and U , which all have size n r U × ( 1 + | Z k , e f | ) . Converting these matrices into 1 × ( n r U ( 1 + | Z k , e f | ) ) column by column, the parameter set of modified Gaussian components is obtained as follows: { w ˜ k | k ( j ) , m ˜ k | k ( j ) , P ˜ k | k ( j ) , ˜ k | k ( j ) } j = 1 n r U ( 1 + | Z k , e f | ) .
This completes the state extraction operation at time k .
The overall procedure of the state extraction of the LGM-PHD filter is given in Appendix A.
Scheme 1 illustrates the process for extracting multiple single-target state estimates based on U . As the elements in each row in U have the same label, they are recorded at the left of each row in U w , as { 401 } , to simplify the space. In Scheme 1a, at time k , there are two classes of confirmed components: { 1 , 2 } ; one unconfirmed component: { 203 } ; four classes of newborn components: { 401 , 402 , 403 , 404 } ; seven effective measurements selected via Section 3.2. At time k 1 , r max and r max t are incremented to r max = 2 and r max t = 203 . Initialize X ^ k = , k , s t = , and T k , s t = .
First, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.8256 , which is U w ( 5 , 2 ) , as shown in Scheme 1a. As w ¯ > w b , obtain U ( 5 , 2 ) = 1 and initialize C N = . Here, | w ¯ | = 1 , x ^ = H k U m ( 5 , 2 ) and R N w ¯ ( i ) = [ 5 , 6 ] . As U ( 5 , 2 ) < V u n , it is the label of the confirmed component. Then, based on (29), survival target status information is obtained. C N = [ C N , 2 ] . Via (32) and (33), U w ( [ 5 , 6 ] , : ) = 0 and U w ( : , 2 ) = 0 , Scheme 1b is obtained. It is obvious that the interference of U w ( 5 , 8 ) caused by z k , e f ( 7 ) with state estimation is prevented, which represents the clutter close to z k , e f ( 1 ) .
Second, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.5262 , which is U w ( 1 , 5 ) , as shown in Scheme 1b. As w ¯ > w b , obtain U ( 1 , 5 ) = 401 and initialize C N = . Here, | w ¯ | = 1 , x ^ = H k U m ( 1 , 5 ) and R N w ¯ ( i ) = 1 . As U ( 1 , 5 ) > V n e w , it is the label of the newborn component. Based on (20), C N w l = 5 , and C N = [ C N , 5 ] . Then, information on one new target status is established via (21), and r max is incremented to 3. The information on this new confirmed track is established by (22).
Additionally, based on (23), C N w b l = 6 , and C N = [ C N , 6 ] . Then, one new unconfirmed track is established via (24–26). Here, r max t is incremented to 204. Via (32) and (33), U w ( 1 , : ) = 0 and U w ( : , [ 5 , 6 ] ) = 0 , Scheme 1c is obtained.
Third, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.5086 , which is U w ( 8 , 4 ) , as shown in Scheme 1c. As w ¯ > w b , obtain U ( 8 , 4 ) = 203 and initialize C N = . Here, | w ¯ | = 1 , x ^ = H k U m ( 8 , 4 ) and R N w ¯ ( i ) = 8 . As V u n < U ( 8 , 4 ) V n e w , it is the label of the unconfirmed component. Since w ¯ > w l , a new target status information is obtained by (27) and C N = [ C N , 4 ] . Via (32) and (33), U w ( 8 , : ) = 0 and U w ( : , 4 ) = 0 , Scheme 1d is obtained.
Fourth, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.4087 , which is U w ( 3 , 7 ) , as shown in Scheme 1d. As w ¯ > w b , obtain U ( 3 , 7 ) = 403 and initialize C N = . Here, | w ¯ | = 1 , x ^ = H k U m ( 3 , 7 ) and R N w ¯ ( i ) = 3 . As U ( 3 , 7 ) > V n e w , it is the label of the newborn component. Based on (20), C N w l = 7 , and C N = [ C N , 7 ] . Then, information on one new target status is established via (21), and r max is incremented to 4. The information on this new confirmed track is established by (22). Via (32) and (33), U w ( 3 , : ) = 0 and U w ( : , 7 ) = 0 , Scheme 1e is obtained.
Fifth, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.2777 , which is U w ( 7 , 3 ) , as shown in Scheme 1e. As w ¯ > w b , obtain U ( 7 , 3 ) = 2 and initialize C N = . Here, | w ¯ | = 1 , x ^ = H k U m ( 7 , 3 ) and R N w ¯ ( i ) = 7 . As U ( 7 , 3 ) < V u n , it is the label of the confirmed component. Then, based on (29), survival target status information is obtained. C N = [ C N , 3 ] . Via (32) and (33), U w ( 7 , : ) = 0 and U w ( : , 3 ) = 0 , Scheme 1f is obtained.
Sixth, the maximum weight in U w ( : , 2 : e n d ) is w ¯ = 0.0003 in Scheme 1f. As w ¯ < w b , step 1 to step 3 of the state extraction at time k is completed. X ^ k , k , s t , and T k , s t are all updated.
Remark 2.
In the proposed filter, since the components approximating one survival target can be identified, the weight threshold for state extraction is reduced compared with the threshold in the basic GM-PHD filter. This ensures that the survival target can be tracked when the proper measurement error occurs, such as U w ( 7 , 3 ) = 0.2777 in Scheme 1. Additionally, only one state estimate can be extracted from the survival Gaussian components with the same label. This ensures that the interference of clutter on the number of state estimates, such as U w ( 5 , 8 ) = 0.4703 in Scheme 1, can be prevented.
Remark 3. 
The behavior of labeling state estimates in the proposed filter is determined by employing the weights and labels of the updated components, which requires no additional computational processing. Furthermore, the updating process and query tables in the state extraction algorithm are all based on the selected observations. Thus, the real-time performance of the LGM-PHD filter is guaranteed.

3.5. Pruning

To reduce the computational cost, reducing the number of Gaussian components by pruning is necessary.
Step 1: Prune confirmed and unconfirmed components that have no state estimates for successive n n o time steps. In our filter, n n o = 6 .
① The label set k , n o of confirmed and unconfirmed components with no state estimates at time k is:
k , n o = { k , n o ( i ) } i = 1 | k , n o | = u n i ( { ˜ k | k ( j ) } j = 1 n r U ( 1 + | Z k , e f | ) ) { γ , k ( i ) } i = 1 J γ , k k , s t
② For i = 1 , , | k , n o | , check the occurrence number of k , n o ( i ) in { k 1 , n o } k 1 = k n n o + 1 k 1 . If k , n o ( i ) appears in k 1 , n o at each time k 1 = k n n o + 1 , , k 1 , { w ˜ k | k ( j ) = 0 | ˜ k | k ( j ) = k , n o ( i ) } j = 1 n r U ( 1 + | Z k , e f | ) . This procedure improves accuracy.
Step 2: Based on the truncation threshold w e t = 10 5 , a certain number of the components with the strongest weights are retained in I ,
I = { j = 1 , , n r U ( 1 + | Z k , e f | ) | w ˜ k | k ( j ) > w e t }
Step 3: Obtain { w k | k ( j ) , m k | k ( j ) , P k | k ( j ) , k | k ( j ) } j = 1 J k | k , where J k | k = | I | , ( w k | k ( j ) , m k | k ( j ) , P k | k ( j ) , k | k ( j ) | j = 1 ) = ( w ˜ k | k ( j ) , m ˜ k | k ( j ) , P ˜ k | k ( j ) , ˜ k | k ( j ) | j = I ( 1 ) ) , among others.

3.6. Merging and Capping

In the basic GM-PHD filter, some of the Gaussian components that are very close together will be approximated by a single Gaussian [32]. This means that, in the merging procedure, the posterior information of targets in close proximity will interfere. This will result in the loss of posterior information, approximated by components with small weights, of some targets. To avoid such loss, in the proposed filter, only the components with the same label can be merged.
First, let d T h be the merging threshold. In general, d T h = 4 . Set J k = 0 , and I 1 = { i 1 } i 1 = 1 J k | k . The specific steps in merging process are as follows.
Step 1: Let J k = J k + 1 , j = arg max i 1 I 1 ( w k | k ( i 1 ) ) . Then, k | k ( j ) is obtained.
Step 2: The sequence number set I 2 is given as follows:
I 2 = { i 1 I 1 | k | k ( i 1 ) = k | k ( j ) , ( m k | k ( i 1 ) m k | k ( j ) ) T ( P k | k ( i 1 ) ) 1 ( m k | k ( i 1 ) m k | k ( j ) ) d T h }
Step 3: Then, a merged Gaussian component is given as follows:
w k ( J k ) = i 1 I 2 w k | k ( i 1 )
m k ( J k ) = 1 w k ( J k ) i 1 I 2 w k | k ( i 1 ) m k | k ( i 1 )
P k ( J k ) = 1 w k ( J k ) i 1 I 2 w k | k ( i 1 ) ( P k | k ( i 1 ) + ( m k ( J k ) m k | k ( i 1 ) ) ( m k ( J k ) m k | k ( i 1 ) ) T )
k ( J k ) = k | k ( j )
Step 4: Execute I 1 = I 1 I 2 . If I 1 is not empty, skip to Step 1.
Step 5: The output of merging procedure is obtained as follows: { w k ( j ) , m k ( j ) , P k ( j ) , k ( j ) } j = 1 J k .
Next, if J k > J max 1 , replace { w k ( j ) , m k ( j ) , P k ( j ) , k ( j ) } j = 1 J k by those of the J max 1 Gaussian components with the largest weights, and the posterior intensity for propagating is obtained as follows: { w k ( j ) , m k ( j ) , P k ( j ) , k ( j ) } j = 1 J k . In this paper, J max 1 = 150 .

3.7. Steps in the Proposed Algorithm

For completeness, we summarized the key steps of the LGM-PHD filter as follows.
Step 1: Prediction. The predicted components { w k | k 1 ( i ) , m k | k 1 ( i ) , P k | k 1 ( i ) } i = 1 J k | k 1 are obtained via (4) and their labels are obtained via (16). Then, obtain the parameter set { w k | k 1 ( i ) , m k | k 1 ( i ) , P k | k 1 ( i ) , k | k 1 ( i ) } i = 1 J k | k 1 of the predicted components.
Step 2: Observation selection. The effective observation set Z k , e f is obtained by eliminating the clutter in the low prior region via Algorithm 1.
Step 3: Updating. Based on Z k , e f , the predicted components are updated via (5) to (12), and the parameters w , m , P , and associated with each component are stored in U and U , respectively.
Step 4: State extraction. As per the steps in 3.4, based on U and U , the multitarget state extraction is converted into multiple single-target state extractions that can provide the identity label, and the occurrence time for each estimate, i.e., X ^ k , k , s t , and T k , s t are obtained. Then, converting the four U values into 1 × ( n r U ( 1 + | Z k , e f | ) ) column by column, the parameter set { w ˜ k | k ( j ) , m ˜ k | k ( j ) , P ˜ k | k ( j ) , ˜ k | k ( j ) } j = 1 n r U ( 1 + | Z k , e f | ) of the modified Gaussian components is obtained.
Step 5: Pruning. As per the steps in Section 3.5, { w k | k ( j ) , m k | k ( j ) , P k | k ( j ) , k | k ( j ) } j = 1 J k | k are obtained.
Step 6: Merging and capping. The closely spaced components with the same label are merged as per the steps in Section 3.6 to obtain { w k ( j ) , m k ( j ) , P k ( j ) , k ( j ) } j = 1 J k for propagation.
After Step 4, based on X ^ k , k , s t , and T k , s t , explicit track maintenance is achieved by connecting the state estimates with the same label in sequence.
Remark 4. 
The proposed filter does not specifically define track disappearance, and each track is left as it is, whether it has disappeared or been misdetected, except that no state estimate appears in seven successive time steps. When a target has been redetected, as long as confirmed Gaussian components for this target are still available for propagation, this target can still have clues to be tracked, and the state estimate will be extracted from these components only if its posterior information satisfies the necessary conditions.

4. Simulation and Results Analysis

The experiment presented in this section aimed to test the performance of the proposed filter. We compared the proposed filter with the basic GM-PHD filter [32] and the Gibbs-GLMB filter [10]. Without loss of generality, we employed the MTT simulation model given in [39]. In a two-dimensional scenario over the region [ 1000 1000 ]   m × [ 1000 1000 ]   m , each target moved according to the following linear Gaussian dynamics:
x k = [ 1 Δ 0 0 0 1 0 0 0 0 1 Δ 0 0 0 1 ] x k 1 + σ ω 2 [ Δ 2 / 4 Δ 3 / 2 0 0 Δ 3 / 2 Δ 2 0 0 0 0 Δ 2 / 4 Δ 3 / 2 0 0 Δ 3 / 2 Δ 2 ]
where x k = [ x 1 , k x 2 , k x 3 , k x 4 , k ] T , [ x 1 , k x 3 , k ] T , and [ x 2 , k x 4 , k ] T are, respectively, the position and velocity of the target at time k . Δ = 1   s   is the sampling time, and σ ω = 5   m / s 2 . Targets can appear and disappear in the scene at any time with the survival probability p S , k = 0.95 . For the Gibbs-GLMB filter, the birth model is a labeled multi-Bernoulli with parameters { r B , k ( i ) , p B , k ( i ) } i = 1 4 , where i = ( k , i ) , r B , k ( i ) = 0.03 and p B , k ( i ) = N ( x ; m B ( i ) , P B ) with m B ( 1 ) = [ 0 0 0 0 ] T , m B ( 2 ) = [ 400 0 600 0 ] T , m B ( 3 ) = [ 800 0 200 0 ] T , m B ( 4 ) = [ 200 0 800 0 ] T , and P B = d i a g ( [ 10 10 10 10 ] T ) 2 . For the GM-PHD filter, the newborn targets appear spontaneously according to a Poisson point process with intensity function γ k ( x k ) = i = 1 4 w γ , k ( i ) N ( x ; m γ , k ( i ) , P γ , k ( i ) ) , which are given by m γ , k ( 1 ) = m B ( 1 ) , m γ , k ( 2 ) = m B ( 2 ) , m γ , k ( 3 ) = m B ( 3 ) , m γ , k ( 4 ) = m B ( 4 ) , P γ , k = P B , and w γ , k ( i ) = r B , k ( i ) . Additionally, γ , k ( 1 ) = 401 , γ , k ( 2 ) = 402 , γ , k ( 3 ) = 403 , and γ , k ( 4 ) = 404 .
The linear target-originated measurement is given as follows:
y k = [ 1 0 0 0 0 0 1 0 ] x k + [ υ x , k υ y , k ]
where υ x , k and υ y , k represent mutually independent zero-mean Gaussian white noise with standard deviations of σ x = 10   m   and σ y = 10   m   , respectively. The clutter is uniformly distributed over the region [ 1000 1000 ]   m × [ 1000 1000 ]   m with an average rate of λ points per scan, i.e., κ = λ / 2000 2 , for the basic GM-PHD filter. To ensure a fair comparison, the basic GM-PHD filter employs the gating algorithm that is used in Gibbs-GLMB filter, and G t h = 9 .
The maximum number of Gaussian components is J max = 100 in the proposed filter and the basic GM-PHD filter. The optimal subpattern assignment (OSPA) metric [40] is used to evaluate the tracking performance of different filters, and the order and cut-off are set as p = 1 and c = 100 , respectively. The presented result is the average of 250 independent repeated experiments but with independently generated observations for each trial. All experiments were run in MATLAB R2013b on a computer with a 3.3 GHz CPU and 8 GB RAM.
Example 1: Detailed performance comparison under λ = 100 and p D , k = 0.9 . As shown in Figure 4, Figure 5 and Figure 6, the number of clutter points is λ = 100 , and the detection probability is p D , k = 0.9 . Figure 4 shows the track estimates obtained from the proposed filter for a single trial [38]. The average OSPA errors and computing time for 250 sample runs of each filter are shown in Figure 5 and Figure 6.
Figure 4a shows that the proposed filter can achieve correct trajectory maintenance in MTT in the presence of dense clutter. As the proposed filter does not supplement the state estimate for the misdetected target, there are discontinuities in the tracks. However, this does not affect the trajectory continuity, and the redetected targets can be tracked in time. Additionally, as shown in Figure 4b, one false track is around H k m γ , k ( 2 ) , one false track is around H k m γ , k ( 3 ) and two false tracks are around H k m γ , k ( 4 ) . This is because the proposed filter has no special procedure to account for the false estimates around the newborn locations. Similarly, they do not affect the continuity of other established tracks, and will disappear over time.
In Figure 5, the average OSPA distances of the proposed filter, the Gibbs-GLMB filter and the basic GM-PHD filter were 22.522 m, 20.497 m, and 31.990 m, respectively. Compared with the basic GM-PHD filter, the average tracking performance of the proposed filter and the Gibbs-GLMB filter increased by 29.60% and 35.92%, respectively. Clearly, the tracking performance of the LGM-PHD filter was slightly lower than that of the Gibbs-GLMB filter. This is mainly caused by the phenomena shown in Figure 4b. However, this did not affect the survival target tracking.
As shown in Figure 6, the average computing times in seconds of a single run (100 time steps) of our MATLAB implementations were 4.60 s (the proposed filter), 30.21 s (the Gibbs-GLMB filter), and 4.96 s (the basic GM-PHD filter), respectively. Compared with the basic GM-PHD filter, the real-time performance of the proposed filter and the Gibbs-GLMB filter increased by 7.26% and −509%, respectively. As expected, the real-time performance of the proposed filter was better than that of the Gibbs-GLMB filter.
Example 2: Performance comparison under various clutter rates for p D , k = 0.95 . In this example, the performance of the proposed filter was evaluated with various clutter rates. The clutter rate λ changed from 1 to 100, and the probability of detection was set to p D , k = 0.95 . Figure 7 shows the comparison results in terms of the average OSPA errors for 250 sample runs of the different algorithms under different clutter rates, and their computing times are shown in Figure 8.
From Figure 7 and Figure 8, we can see that the proposed filter achieved a better OSPA distance than the basic GM-PHD filter and kept a similar real-time performance as the clutter rate increased. Compared with the basic GM-PHD filter, with λ = 20 , the average tracking performances of the proposed filter and the Gibbs-GLMB filter increased by 15.75% and 33.14%, respectively; with λ = 100 , they increased by 31.96% and 37.14%, respectively. The advantage of the proposed filter can be fully demonstrated in the dense clutter environment.
Example 3: Performance comparison under various probabilities of detection for λ = 100 . In this example, various probabilities of detection were employed to evaluate the effectiveness of the proposed filter. The probability of detection p D , k changed from 0.8 to 1, and the clutter rate was kept unchanged at λ = 100 . The comparison results of the average OSPA errors for 250 sample runs of the different algorithms under different probabilities of detection are shown in Figure 9. Figure 10 shows the corresponding computing times of each filter.
Clearly, the proposed filter achieved a better OSPA distance than the basic GM-PHD filter and maintained a similar real-time performance at each probability of detection. As the probability of detection increased, the OSPA distance of the proposed filter decreased. Compared with the basic GM-PHD filter, with p D , k = 0.84 , the average tracking performances of the proposed filter and the Gibbs-GLMB filter increased by 23.03% and 29.13%, respectively; with p D , k = 0.96 , they increased by 33.21% and 35.79%, respectively.

5. Conclusions

The paper proposed a computationally efficient LGM-PHD filter that can explicitly accommodate the estimation of target trajectories in MTT. Based on the proposed query matrices, state estimates with identity labels and occurrence times were extracted, and no additional associated procedures were required. Three numerical and simulations analyses demonstrated that the LGM-PHD filter could effectively improve the tracking performance as compared with the basic GM-PHD filter, which was slightly lower than that of the Gibbs-GLMB filter. As the proposed filter could maintain a similar to or slightly lower real-time performance than that of the basic GM-PHD filter, it could be an alternative filter explicitly accommodating target trajectories in real-time processing applications. We should point out that the LGM-PHD filter is not suitable for low detection scenarios. In our future work, we will investigate multisensor technology [12,13,39] for improving tracking performance.

Author Contributions

Y.G. performed the simulations and wrote the paper; D.J. and C.Z. revised this paper and offered some useful suggestions in methodologies and simulations; S.G. offered revision during the whole process. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Fundamental Research Funds for the Central Universities under grant B200202165 and by National Natural Science Foundation of China under grant 61971179.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1. The overall procedure of the state extraction. The operations in the red boxes are step 1, the operations in the black boxes are step 2, the operation in the blue box is step 3, and the operations in the green box are step 4.
Figure A1. The overall procedure of the state extraction. The operations in the red boxes are step 1, the operations in the black boxes are step 2, the operation in the blue box is step 3, and the operations in the green box are step 4.
Sensors 21 03932 g0a1

References

  1. Mahler, R.P.S. Statistical Multisource-Multitarget Information Fusion; Artech House: Norwood, MA, USA, 2007. [Google Scholar]
  2. Blackman, S.; Popoli, R. Design and Analysis of Modern Tracking Systems; Artech House: Norwood, MA, USA, 1999. [Google Scholar]
  3. Bar-Shalom, Y.; Willett, P.K.; Tian, X. Tracking and Data Fusion: A Handbook of Algorithms; YBS Publishing: Storrs, CT, USA, 2011. [Google Scholar]
  4. Blackman, S. Multiple hypothesis tracking for multiple target tracking. IEEE Aerosp. Electron. Syst. Mag. 2004, 19, 5–18. [Google Scholar] [CrossRef]
  5. Jiang, D.; Li, M.; Gao, Y.; Gao, Y.; Fu, W.; Han, Y. Time-Matching Random Finite Set-Based Filter for Radar Multi-Target Tracking. Sensors 2018, 18, 4416. [Google Scholar] [CrossRef] [Green Version]
  6. Mahler, R.P.S. Multi-target Bayes filtering via first-order multi-target moments. IEEE Trans. Aerosp. Electron. Syst. 2003, 39, 1152–1178. [Google Scholar] [CrossRef]
  7. Gao, Y.; Jiang, D.; Liu, M. Particle-gating SMC-PHD filter. Signal Process. 2017, 130, 64–73. [Google Scholar] [CrossRef]
  8. Mahler, R.P.S. PHD filters of higher order in target number. IEEE Trans. Aerosp. Electron. Syst. 2007, 43, 1523–1543. [Google Scholar] [CrossRef]
  9. Vo, B.T.; Vo, B.N.; Cantoni, A. The cardinality balanced multi-target multi-Bernoulli filter and its implementations. IEEE Trans. Signal Process. 2009, 57, 409–423. [Google Scholar]
  10. Vo, B.N.; Vo, B.T.; Hoang, H. An Efficient Implementation of the Generalized Labeled Multi-Bernoulli Filter. IEEE Trans. Signal Process. 2017, 65, 1975–1987. [Google Scholar] [CrossRef] [Green Version]
  11. Wang, S.; Bao, Q.; Chen, Z. Refined PHD Filter for Multi-target Tracking under Low Detection Probability. Sensors 2019, 19, 2842. [Google Scholar] [CrossRef] [Green Version]
  12. Chai, L.; Yi, W.; Jiang, X.; Kong, L. A Distributed PHD Filter for On-line Joint Sensor Registration and Multi-target Tracking. In Proceedings of the International Conference on Information Fusion, Ottawa, ON, Canada, 2–5 July 2019; pp. 1–6. [Google Scholar]
  13. Li, T.; Mallick, M.; Pan, Q. A parallel Filtering-Communication-Based Cardinality Consensus Approach for Real-Time Distributed PHD Filtering. IEEE Sens. J. 2020, 20, 13824–13832. [Google Scholar] [CrossRef]
  14. Vo, B.T.; Vo, B.N. Labeled Random Finite Sets and Multi-Object Conjugate Priors. IEEE Trans. Signal Process. 2013, 61, 3460–3475. [Google Scholar] [CrossRef]
  15. Reuter, S.; Vo, B.T.; Vo, B.N.; Dietmayer, K. The labelled multi-Bernoulli filter. IEEE Trans. Signal Process. 2014, 62, 3246–3260. [Google Scholar]
  16. Tian, Z.; Liu, W.; Ru, X. Multi-Target Localization and Tracking Using TDOA and AOA Measurement Based on Gibbs-GLMB Filtering. Sensors 2019, 19, 5437. [Google Scholar] [CrossRef] [Green Version]
  17. Liu, R.; Fan, H.; Li, T.; Xiao, H. A Computationally Efficient Labeled Multi-Bernoulli Smoother for the Multi-Target Tracking. Sensors 2019, 19, 4226. [Google Scholar] [CrossRef] [Green Version]
  18. García-Fernández, Á.F.; Svensson, L.; Williams, J.L.; Xia, Y.; Granström, K. Trajectory Poisson Multi-Bernoulli Filters. IEEE Trans. Signal Process. 2020, 68, 4933–4945. [Google Scholar] [CrossRef]
  19. Liang, G.; Li, Q.; Qi, B.; Qiu, L. Multitarget Tracking Using One Time Step Lagged Delta-Generalized Labeled Multi-Bernoulli Smoothing. IEEE Access. 2020, 8, 28242–28256. [Google Scholar] [CrossRef]
  20. Reuter, S.; Danzer, A.; Stübler, M.; Scheel, A.; Granström, K. A fast implementation of the Labeled Multi-Bernoulli filter using gibbs sampling. In Proceedings of the IEEE Intelligent Vehicles Symposium (IV), Los Angeles, CA, USA, 11–14 June 2017; pp. 765–772. [Google Scholar]
  21. Beard, M.; Reuter, S.; Granström, K.; Vo, B.T.; Vo, B.N.; Scheel, A. Multiple Extended Target Tracking With Labeled Random Finite Sets. IEEE Trans. Signal Process. 2016, 64, 1638–1653. [Google Scholar] [CrossRef] [Green Version]
  22. García-Fernández, Á.F.; Svensson, L. Trajectory probability hypothesis density filter. In Proceedings of the International Conference on Information Fusion, Cambridge, UK, 10–13 July 2018; pp. 1430–1437. [Google Scholar]
  23. Granström, K.; Willett, P.; Bar-Shalom, Y. Approximate Multi-Hypothesis Multi-Bernoulli Multi-Object Filtering Made Multi-Easy. IEEE Trans. Signal Process. 2016, 64, 1784–1797. [Google Scholar] [CrossRef] [Green Version]
  24. Lin, L.; Bar-shalom, Y.; Kirubarajan, T. Tracking labeling and PHD filter for multitarget tracking. IEEE Trans. Aerosp. Electron. Syst. 2006, 42, 778–795. [Google Scholar] [CrossRef]
  25. Yang, J.; Ji, H. A novel track maintenance algorithm for PHD/CPHD filter. Signal Process. 2012, 92, 2371–2380. [Google Scholar] [CrossRef]
  26. Li, Y.; Xiao, H.; Wu, H.; Fu, Q.; Hu, R. Modified Labeled Particle Probability Hypothesis Density Filter for Joint Multi-target Tracking and Classification. In Proceedings of the International Conference on Signal Processing and Communication Systems, Cairns, QLD, Australia, 14–16 December 2015; pp. 1–5. [Google Scholar]
  27. Li, T.; Sun, S.; Corchado, J.M.; Siyan, M.F. A particle dyeing approach for track continuity for the SMC-PHD filter. In Proceedings of the International Conference on Information Fusion, Slamanca, Spain, 7–10 July 2014; pp. 1–8. [Google Scholar]
  28. Cao, C.; Zhao, Y.; Pang, X.; Suo, Z.; Chen, S. An efficient implementation of multiple weak targets tracking filter with labeled random finite sets for marine radar. Digital Signal Process. 2020, 101, 102710. [Google Scholar] [CrossRef]
  29. Gao, Y.; Jiang, D.; Liu, M.; Fu, W. An Explicit Track Continuity Algorithm for the SMC-PHD filter. In Proceedings of the International Conference on Communication Technology, Chengdu, China, 27–30 October 2017; pp. 1702–1708. [Google Scholar]
  30. García-Fernández, Á.F.; Svensson, L. Trajectory PHD and CPHD Filters. IEEE Trans. Signal Process. 2019, 67, 5702–5714. [Google Scholar] [CrossRef] [Green Version]
  31. Vo, B.N.; Singh, S.; Doucet, A. Sequential Monte Carlo methods for multi-target filtering with random finite sets. IEEE Trans. Aerosp. Electron. Syst. 2005, 41, 1224–1245. [Google Scholar]
  32. Vo, B.N.; Ma, W.K. The Gaussian mixture probability hypothesis density filter. IEEE Trans. Signal Process. 2006, 54, 4091–4104. [Google Scholar] [CrossRef]
  33. Zhang, H.; Ge, H.; Yang, J.; Yuan, Y. A GM-PHD algorithm for multiple target tracking based on false detection irregular window. Signal Process. 2016, 120, 537–552. [Google Scholar] [CrossRef]
  34. Choi, M.E.; Seo, S.W. Robust Multitarget Tracking Scheme Based on Gaussian Mixture Probability Hypothesis Density Filter. IEEE Trans. Veh. Technol. 2016, 65, 4217–4229. [Google Scholar] [CrossRef]
  35. Dehkordi, M.; Azimifar, Z. Refined GM-PHD tracker for tracking targets in possible subsequent missed detections. Signal Process. 2015, 116, 112–126. [Google Scholar] [CrossRef]
  36. Gao, L.; Liu, H.; Liu, H. Probability hypothesis density filter with imperfect detection probability for mutli-target tracking. Optik 2016, 127, 10428–10436. [Google Scholar] [CrossRef]
  37. Zhang, H.; Gao, L.; Xu, M.; Wang, Y. An improved probability hypothesis density filter for multi-target tracking. Optik 2019, 182, 23–31. [Google Scholar] [CrossRef]
  38. Gao, Y.; Fang, N.; Jiang, D.; Fu, W.; Guo, S. An explicit track continuity algorithm for the GM-PHD filter. In Proceedings of the 31st Chinese Control and Decision Conference, Nanchang, China, 3–5 June 2019; pp. 1072–1079. [Google Scholar]
  39. Vo, B.N.; Vo, B.T. An implementation of the Multi-sensor Generalized Labeled Multi-Bernoulli Filter via Gibbs Sampling. In Proceedings of the International Conference on Information Fusion, Xi’an, China, 10–13 July 2017; pp. 1–8. [Google Scholar]
  40. Schuhmacher, D.; Vo, B.T.; Vo, B.N. A consistent metric for performance evaluation in multi-object filters. IEEE Trans. Signal Process. 2008, 56, 3447–3457. [Google Scholar] [CrossRef] [Green Version]
Scheme 1. Steps of extracting multiple single-target state estimate based on U . (af): the successive steps of extracting multiple single-target state estimates.
Scheme 1. Steps of extracting multiple single-target state estimate based on U . (af): the successive steps of extracting multiple single-target state estimates.
Sensors 21 03932 sch001aSensors 21 03932 sch001b
Figure 1. Overall procedure of the basic GM-PHD filter.
Figure 1. Overall procedure of the basic GM-PHD filter.
Sensors 21 03932 g001
Figure 2. Phenomena affecting track continuity of the basic GM-PHD filter. (a) detection uncertainty; (b) false state estimates; (c) closely spaced targets.
Figure 2. Phenomena affecting track continuity of the basic GM-PHD filter. (a) detection uncertainty; (b) false state estimates; (c) closely spaced targets.
Sensors 21 03932 g002
Figure 3. Overall procedure of the LGM-PHD filter.
Figure 3. Overall procedure of the LGM-PHD filter.
Sensors 21 03932 g003
Figure 4. Track estimates of the proposed filter and ground truths for a sample run ( λ = 100 , p D , k = 0.9 ). (a) The true and estimated trajectories of targets. Different trajectories estimated by the proposed filter are denoted by different colored dots. (b) The true and estimated trajectories of targets in x and y coordinates, respectively. © 2021 IEEE.
Figure 4. Track estimates of the proposed filter and ground truths for a sample run ( λ = 100 , p D , k = 0.9 ). (a) The true and estimated trajectories of targets. Different trajectories estimated by the proposed filter are denoted by different colored dots. (b) The true and estimated trajectories of targets in x and y coordinates, respectively. © 2021 IEEE.
Sensors 21 03932 g004
Figure 5. Comparison of different algorithms for 250 sample runs ( λ = 100 , p D , k = 0.9 ). (a) Average Optimal subpattern assignment (OSPA) distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Figure 5. Comparison of different algorithms for 250 sample runs ( λ = 100 , p D , k = 0.9 ). (a) Average Optimal subpattern assignment (OSPA) distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Sensors 21 03932 g005
Figure 6. Computing times of different algorithms for 250 sample runs ( λ = 100 , p D , k = 0.9 ).
Figure 6. Computing times of different algorithms for 250 sample runs ( λ = 100 , p D , k = 0.9 ).
Sensors 21 03932 g006
Figure 7. Comparison of different algorithms under various clutter rates for 250 sample runs ( p D , k = 0.95 ). (a) Average OSPA distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Figure 7. Comparison of different algorithms under various clutter rates for 250 sample runs ( p D , k = 0.95 ). (a) Average OSPA distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Sensors 21 03932 g007
Figure 8. Computing times of different algorithms under various clutter rates for 250 sample runs ( p D , k = 0.95 ).
Figure 8. Computing times of different algorithms under various clutter rates for 250 sample runs ( p D , k = 0.95 ).
Sensors 21 03932 g008
Figure 9. Comparison of different algorithms under various probabilities of detection for 250 sample runs ( λ = 100 ). (a) Average OSPA distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Figure 9. Comparison of different algorithms under various probabilities of detection for 250 sample runs ( λ = 100 ). (a) Average OSPA distances of different filters; (b) OSPA locations of different filters; (c) OSPA cardinalities of different filters.
Sensors 21 03932 g009
Figure 10. Computing times of different algorithms under various probabilities of detection for 250 sample runs ( λ = 100 ).
Figure 10. Computing times of different algorithms under various probabilities of detection for 250 sample runs ( λ = 100 ).
Sensors 21 03932 g010
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Gao, Y.; Jiang, D.; Zhang, C.; Guo, S. A Labeled GM-PHD Filter for Explicitly Tracking Multiple Targets. Sensors 2021, 21, 3932. https://0-doi-org.brum.beds.ac.uk/10.3390/s21113932

AMA Style

Gao Y, Jiang D, Zhang C, Guo S. A Labeled GM-PHD Filter for Explicitly Tracking Multiple Targets. Sensors. 2021; 21(11):3932. https://0-doi-org.brum.beds.ac.uk/10.3390/s21113932

Chicago/Turabian Style

Gao, Yiyue, Defu Jiang, Chao Zhang, and Su Guo. 2021. "A Labeled GM-PHD Filter for Explicitly Tracking Multiple Targets" Sensors 21, no. 11: 3932. https://0-doi-org.brum.beds.ac.uk/10.3390/s21113932

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop