# Unknown Clutter Estimation by FMM Approach in Multitarget Tracking Algorithm

Authors: Ning Lv; Feng Lian; Chongzhao Han
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Finite mixture model (FMM) approach is a research focus in multitarget tracking field. The clutter was treated as uniform distribution previously. Aiming at severe bias caused by unknown and complex clutter, a multitarget tracking algorithm based on clutter model estimation is put forward in this paper. Multitarget likelihood function is established with FMM. In this frame, the algorithms of expectation maximum (EM) and Markov Chain Monte Carlo (MCMC) are both consulted in FMM parameters estimation. Furthermore, target number and multitarget states can be estimated precisely after the clutter model fitted. Association between target and measurement can be avoided. Simulation proved that the proposed algorithm has a good performance in dealing with unknown and complex clutter.


---

## Body

## 1. Introduction

Multiple target tracking (MTT), for its important theoretical significance and widely applied engineering background, is a research focus in tracking field in recent years [1–3]. In MTT problems, measurement set of sensors contains not only the measurements from target, but also massive clutters from interference caused by meteorological phenomena, electromagnetism environment, and false target. In addition, measurements from target and clutter cannot be distinguished usually. How to estimate target number and multitarget states using these mixture measurements is the key.Until now, MTT solutions could be concluded in two classes. The first is data association solutions, such as the nearest neighbor (NN) method [4], the joint probabilistic data association (JPDA) method [5], and the multiple hypothesis tracking (MHT) method [6]. In this class, corresponding relationship between target and measurement should be established before target number and multitarget states estimate. While the second solution dealing with MTT problems is random finite sets (RFS) method based solutions, such as the probability hypothesis density (PHD) filter [7] and the cardinalized probability hypothesis density (CPHD) filter [8]. Target state set can be updated with measurement set directly. MTT in frame of RFS does not need data association consequently.The mentioned MTT solutions are all based on the known clutter model. For randomness of clutter, its distribution model usually includes clutter number and clutter position. If there is not so many interference factors in surveillance area, clutter number could be considered to obey Poisson distribution, while clutter position obeys uniform distribution. But in many actual scenes, especially ground and sea level surveillance, even battlefield surveillance, the clutter model appears to be unknown and more complex as results of complicated landform, jamming station, and unidentified interference source such as electronic countermeasure systems. In this case, assumptions of Poisson and uniform distribution which the clutter model satisfied will lead to severe bias estimated by the filter.A novel MTT algorithm based on clutter model estimation is put forward in the light of the problems mentioned above. In this algorithm, multitarget likelihood function is established with finite mixture model (FMM) [9], whose parameters can be estimated by the algorithms of expectation maximum (EM) and Markov Chain Monte Carlo (MCMC). Furthermore, target number and multitarget states can be estimated as well as the clutter model fitted. Similarly with the RFS based solution, association between target and measurement can be avoided in this algorithm. Compared with the MTT solution without clutter model fitting, it can be proved from a simulation that the algorithm proposed in this paper is more efficient.


## 2. Problem Description

### 2.1. Target Motion Model

Suppose that the sensor was monitoring a fixed region. For situation of spontaneous birth, spawned by existent targets and extinction, the number of target detected over the surveillance region varies with time. Supposing thatt
k is the number of existing targets at time k, we model the motion of the multitarget system as
(1)
x
k
+
1
i
=
f
k
i
(
x
k
i
)
+
ω
k
i
,
i
=
1
,
…
,
t
k
,
where x
k
i
=
[
x
k
i
,
x
˙
k
i
,
x
¨
k
i
,
y
k
i
,
y
˙
k
i
,
y
¨
k
i
]
T represents the state vector of target i, including position, velocity, and acceleration information and X
k
=
{
x
k
1
,
…
,
x
k
T
k
} represents multitarget state set at time k, while ω
k
i represents the process noise vector, obeying Gaussian distribution ω
k
i
~
𝒩
(
0
,
Q
k
i
), where Q
k
i denotes the variance matrix of process noise.


### 2.2. Measurement Model

The measurement set received by sensor at timek could be represented as Z
k
=
{
z
k
1
,
…
,
z
k
n
k
}, where n
k denotes the number of measurement at time k. Sensor measurements are generally regarded as a mixture of target-originated measurement and clutter-originated measurement. Measurements from different targets and clutters are statistically independent in this paper without additional illustration.


#### 2.2.1. Target-Originated Measurement

Supposing that measurementj is originated from target i,
(2)
z
k
j
=
h
k
(
x
k
i
)
+
υ
k
j
,
z
k
j
∈
Z
k
,
where h
k
(
·
) denotes the measurement function of sensor and υ
k
j denotes the measurement noise vector, obeying Gaussian distribution υ
k
j
~
𝒩
(
0
,
R
k
) as well, where R
k denotes the variance matrix of measurement noise.


#### 2.2.2. Clutter Model

Assume that at timek
ℱ
c
,
k denotes the model of clutter position distribution, and measurement j is originated from clutter, then
(3)
z
k
j
~
ℱ
c
,
k
,
z
k
j
∈
Z
k
,
where ℱ
c
,
k is unknown and varies with time.So far as mentioned above, at timek, the information needed to be estimated includes the distribution model of clutter position ℱ
c
,
k, target number t
k, and multitarget state set X
k. In addition, the number model of clutter need not be estimated, because this algorithm is adapted to any variation of the clutter number.


## 2.1. Target Motion Model

Suppose that the sensor was monitoring a fixed region. For situation of spontaneous birth, spawned by existent targets and extinction, the number of target detected over the surveillance region varies with time. Supposing thatt
k is the number of existing targets at time k, we model the motion of the multitarget system as
(1)
x
k
+
1
i
=
f
k
i
(
x
k
i
)
+
ω
k
i
,
i
=
1
,
…
,
t
k
,
where x
k
i
=
[
x
k
i
,
x
˙
k
i
,
x
¨
k
i
,
y
k
i
,
y
˙
k
i
,
y
¨
k
i
]
T represents the state vector of target i, including position, velocity, and acceleration information and X
k
=
{
x
k
1
,
…
,
x
k
T
k
} represents multitarget state set at time k, while ω
k
i represents the process noise vector, obeying Gaussian distribution ω
k
i
~
𝒩
(
0
,
Q
k
i
), where Q
k
i denotes the variance matrix of process noise.


## 2.2. Measurement Model

The measurement set received by sensor at timek could be represented as Z
k
=
{
z
k
1
,
…
,
z
k
n
k
}, where n
k denotes the number of measurement at time k. Sensor measurements are generally regarded as a mixture of target-originated measurement and clutter-originated measurement. Measurements from different targets and clutters are statistically independent in this paper without additional illustration.


### 2.2.1. Target-Originated Measurement

Supposing that measurementj is originated from target i,
(2)
z
k
j
=
h
k
(
x
k
i
)
+
υ
k
j
,
z
k
j
∈
Z
k
,
where h
k
(
·
) denotes the measurement function of sensor and υ
k
j denotes the measurement noise vector, obeying Gaussian distribution υ
k
j
~
𝒩
(
0
,
R
k
) as well, where R
k denotes the variance matrix of measurement noise.


### 2.2.2. Clutter Model

Assume that at timek
ℱ
c
,
k denotes the model of clutter position distribution, and measurement j is originated from clutter, then
(3)
z
k
j
~
ℱ
c
,
k
,
z
k
j
∈
Z
k
,
where ℱ
c
,
k is unknown and varies with time.So far as mentioned above, at timek, the information needed to be estimated includes the distribution model of clutter position ℱ
c
,
k, target number t
k, and multitarget state set X
k. In addition, the number model of clutter need not be estimated, because this algorithm is adapted to any variation of the clutter number.


## 2.2.1. Target-Originated Measurement

Supposing that measurementj is originated from target i,
(2)
z
k
j
=
h
k
(
x
k
i
)
+
υ
k
j
,
z
k
j
∈
Z
k
,
where h
k
(
·
) denotes the measurement function of sensor and υ
k
j denotes the measurement noise vector, obeying Gaussian distribution υ
k
j
~
𝒩
(
0
,
R
k
) as well, where R
k denotes the variance matrix of measurement noise.


## 2.2.2. Clutter Model

Assume that at timek
ℱ
c
,
k denotes the model of clutter position distribution, and measurement j is originated from clutter, then
(3)
z
k
j
~
ℱ
c
,
k
,
z
k
j
∈
Z
k
,
where ℱ
c
,
k is unknown and varies with time.So far as mentioned above, at timek, the information needed to be estimated includes the distribution model of clutter position ℱ
c
,
k, target number t
k, and multitarget state set X
k. In addition, the number model of clutter need not be estimated, because this algorithm is adapted to any variation of the clutter number.


## 3. Multitarget Likelihood Function Based on FMM

Given the measurement setZ
k and condition of independence, the multitarget likelihood function at time k can be described by FMM [9, 10] as
(4)
ℒ
k
(
Z
k
;
ψ
k
)
=
∏
j
=
1
n
k
(
ℱ
c
,
k
(
z
k
i
;
ψ
c
,
k
)
+
ℱ
t
,
k
(
z
k
i
;
ψ
t
,
k
)
)
.Finite mixture model is considered an effective method dealing with multitarget tracking problem, especially under complex unknown clutter environment. Distribution model of target-originated measurement and clutter measurement could be described as a superposition of some normative distribution. Taking complex clutter, for example, it can be considered as a superposition of uniform and finite Gaussian models. By estimating the parameters of these potential models, we can get the multitarget state ultimately.Formula (4) could be explained as follows:ℱ
c
,
k
(
·
;
ψ
c
,
k
) denotes the distribution model of clutter position at time k. Considering the complexity of clutter distribution model, multiple Gaussian models and one uniform model will be used to fit this clutter model:
(5)
ℱ
c
,
k
(
·
;
ψ
c
,
k
)
=
π
c
,
k
1
𝒰
(
·
)
+
∑
i
=
2
g
c
,
k
π
c
,
k
i
𝒩
(
·
;
θ
c
,
k
i
)
,
where ψ
c
,
k
=
{
g
c
,
k
,
π
c
,
k
1
,
…
,
π
c
,
k
g
c
,
k
,
θ
c
,
k
2
,
θ
c
,
k
g
c
,
k
}, g
c
,
k denotes the number of clutter model, π
c
,
k
i denotes the weights of clutter model, and θ
c
,
k
i
=
(
μ
c
,
k
i
,
Σ
c
,
k
i
), composed of the mean vector and the covariance matrix, denotes the parameters of clutter model.ℱ
t
,
k
(
·
;
ψ
t
,
k
) denotes the distribution model of target-originated measurement position at time k. Assuming that the number of existing targets in surveillance region is t
k,
(6)
ℱ
t
,
k
(
·
;
ψ
t
,
k
)
=
∑
i
=
1
g
t
,
k
π
t
,
k
i
f
t
,
k
i
(
·
;
θ
t
,
k
i
)
,
where ψ
t
,
k
=
{
g
t
,
k
,
π
t
,
k
1
,
…
,
π
t
,
k
g
t
,
k
,
θ
t
,
k
t
,
θ
t
,
k
g
t
,
k
} and g
t
,
k
=
t
k, which means that the number of models is equal to the number of targets. π
t
,
k
i denotes the weights of target model, f
t
,
k
i
(
·
;
θ
t
,
k
i
) denotes the distribution model of measurement from a single target, and θ
t
,
k
i denotes the relevant parameters.Lettingg
k
=
g
c
,
k
+
g
t
,
k, all parameters to be estimated in the mixture model can be represented as ψ
k
=
{
g
k
,
π
k
1
,
…
,
π
k
g
k
,
θ
k
1
,
θ
k
g
k
}. consequently, (4) can be rewritten as
(7)
ℒ
k
(
Z
k
;
ψ
k
)
=
∏
j
=
1
n
k
∑
i
=
1
g
k
π
k
i
f
k
i
(
z
k
i
;
θ
k
i
)
,
∑
i
=
1
g
k
π
k
i
=
1
,
where π
k
i
n
k is the number of targets at time k.EM and MCMC are the approaches most widely consulted in FMM parameters estimation. Until now, EM approach is considered a standard algorithm, but it is sensitive to initial value. The iteration will converge to some local extremum caused by worse initialization. Meanwhile the convergence rate will be affected. Comparatively, MCMC belongs to a stochastic algorithm, which is insensitive to initial value. In addition, for some complex situation, MCMC algorithm will lead a better global Convergence. In fact, if the Markov chain is long enough, MCMC approach can obtain massive information of posterior distribution, so local extremum could be avoided. Above all, robustness of MCMC is better than EM, nevertheless, computational complexity of MCMC is larger than EM. Now we will present these two approaches, respectively, in Sections4 and 5.


## 4. FMM Parameters Estimation by EM Approach

### 4.1. Parameters Initialization

Parameters in FMM should be well initialized as much as possible before estimation in EM algorithm. Initializations of clutter model and target-originated measurement model were discussed, respectively.


#### 4.1.1. Initialization of Clutter Model

Considering that the clutter model could vary with time, the initialization process at timek is listed as follows.Inheritance of the value estimated at timet
k
-
1 is
(8)
μ
c
,
k
i
(
0
)
=
μ
^
c
,
k
-
1
i
,
Σ
c
,
k
i
(
0
)
=
Σ
^
c
,
k
-
1
i
,
i
=
1
,
…
,
g
^
c
,
k
-
1
.Then add some clutter points randomly; the number of the clutter points isg
0. The initialization average μ
c
,
k
i
(
0
) can be got from these chosen points in surveillance region S. The covariance matrix could be initialized as Σ
c
,
k
i
(
0
)
=
σ
2
I, i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0, where
(9)
σ
2
=
1
10
d
trace
(
1
n
k
∑
j
=
1
n
k
(
z
k
j
-
z
¯
k
)
(
z
k
j
-
z
¯
k
)
T
)
,
where z
¯
k
=
∑
j
=
1
n
k
z
¯
k
j
/
n
k denotes the average of measurement data.Consequently, the number of clutter models in initialization processg
c
,
k
(
0
)
=
g
^
c
,
k
-
1
+
g
0. In order to ensure the parameters of clutter models converging to the true value, the value of g
0 should be bigger than the number in reality.


#### 4.1.2. Initialization of Target-Originated Measurement Model

At timek, targets totally are composed of survival targets and spontaneous birth targets and spawned by existent targets. So the number of target models in initialization process is
(10)
g
t
,
k
(
0
)
=
g
s
,
k
(
0
)
+
g
b
,
k
(
0
)
+
g
p
,
k
(
0
)
,
where g
s
,
k
(
0
) denotes the number of survival targets, g
b
,
k
(
0
) denotes the number of spontaneous birth targets, and g
p
,
k
(
0
) denotes the number of spawned by existent targets.The initialization of target-originated measurement model at timek can be represented as
(11)
ℱ
t
,
k
(
·
;
ψ
c
,
k
(
0
)
)
=
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
+
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
+
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
.(a) To Survival Targets. Measurement noise obeys Gaussian distribution, as shown in (2), so f
s
,
k
i
(
·
;
θ
s
,
k
i
)
~
𝒩
(
·
;
μ
s
,
k
i
,
Σ
s
,
k
i
). Furthermore, measurement model initialization of survival targets is
(12)
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
=
∑
i
=
1
g
s
,
k
(
0
)
π
s
,
k
i
(
0
)
𝒩
(
·
;
μ
s
,
k
i
(
0
)
,
Σ
c
,
k
i
(
0
)
)
.
Let g
s
,
k
(
0
)
=
g
^
t
,
k
-
1.The initialization averageμ
s
,
k
i
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
) could be got from prediction of target state. The covariance is known, the same with measurement noise of the sensor, Σ
s
,
k
i
(
0
)
=
R
k
,
i
=
1
,
…
,
g
s
,
k
(
0
).(b) To Spontaneous Birth Targets. As is mentioned in [11–13], assume that the initial position model obeys Gaussian distribution either; that is,
(13)
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
=
∑
i
=
1
g
b
,
k
(
0
)
π
b
,
k
i
(
0
)
𝒩
(
·
;
μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)
)
,
where parameters of model initialization g
b
,
k
(
0
), μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
) are all known, according to prior information.(c) To Spawned by Existent Targets. Supposing that each original target can create p
k new targets at most, model initialization can be represented as [11–13]
(14)
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
=
∑
j
=
1
p
k
∑
i
=
1
g
p
,
k
(
0
)
π
p
,
k
i
j
(
0
)
𝒩
(
·
;
μ
p
,
k
i
j
(
0
)
,
Σ
p
,
k
i
j
(
0
)
)
.The target number and state average can be assumed asg
p
,
k
(
0
)
=
p
k
g
s
,
k
(
0
), μ
p
,
k
i
,
j
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
)
+
d
p
.
k
j, where i
=
1
,
…
,
g
s
,
k
(
0
), j
=
1
,
…
,
p
k.The parametersd
p
.
k
j, Σ
p
,
k
i
j
(
0
) could be valued according to prior information similarly.The weights of models could be briefly treated as the same in initialization process:(15)
π
k
i
(
0
)
=
1
g
k
(
0
)
,
i
=
1
,
…
,
g
c
,
k
(
0
)
+
g
t
,
k
(
0
)
.


### 4.2. Estimation of FMM Parameters

#### 4.2.1. EM Approach

As FMM Parameters Estimation with EM approach, the loss variables can be treated as element labelingE
k
=
{
e
k
1
,
…
,
e
k
n
k
}, where e
k
1 is a vector of n
k dimension. The element e
k
i
j
=
1 or 0 illustrates whether z
k
j is originated from ith element of the FMM or not. So the complete dataset at time k is Y
k
=
{
Z
k
,
E
k
}, and logarithm likelihood function of the complete data is
(16)
log
⁡
𝒞
k
(
Y
k
;
ψ
k
)
=
∑
i
=
1
g
k
∑
j
=
1
n
k
e
k
i
j
{
log
⁡
⁡
π
k
i
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
)
}
.When the number of modelsg
k is known, EM algorithm can be iterated by E-step and M-step. But the number of models is unknown in this paper. Some criteria will be used for number estimation, such as minimum message length (MML) criterion and Bayesian information criterion (BIC) [10, 14]. Here in this paper, for the number of models in initialization being more than reality, the technology of model merging and pruning could be inserted into each step of iteration in EM algorithm. In this way, the number of FMM can be estimated. This method is more intuitive than those based on criteria. Now, this modified EM algorithm is described as follows:(a) E-Step.Conditional expectations of 𝒞
k
(
Y
k
;
ψ
k
) will be iterated. The expectations of loss data e
k
i
j at tth step is
(17)
E
ψ
k
(
t
)
[
e
k
i
j
;
Z
k
]
=
p
r
ψ
k
(
t
)
{
e
k
i
j
=
1
∣
Z
k
}
=
τ
k
i
(
z
k
j
;
ψ
k
(
t
)
)
=
τ
k
i
j
(
t
)
,
where τ
k
i
j
(
t
) represents the posterior probability of z
k
j belonging to model i. The formula for computing τ
k
i
j
(
t
) is
(18)
τ
k
i
j
(
t
)
=
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
∑
i
=
1
g
k
(
t
)
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
.According to the formula above, the conditional expectations oflog
⁡
𝒞
k
(
Y
k
;
ψ
k
) is
(19)
Q
(
ψ
k
;
ψ
k
(
t
)
)
=
E
ψ
k
(
t
)
{
log
⁡
⁡
𝒞
k
(
Y
k
;
ψ
k
)
∣
Z
k
}
=
∑
i
=
1
g
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
{
log
⁡
⁡
π
k
i
(
t
)
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
}
.In fact, the loss variable has reflected the association relationship between the measurement and target, that is, why association process can be avoided in this algorithm.(b) M-Step. The value of ψ
k will be estimated throught global maximum of Q
(
ψ
k
;
ψ
k
(
t
)
). By solving ∂
Q
k
(
ψ
k
,
ψ
k
(
t
)
)
/
∂
ψ
k
=
0, the weights of all models will be reestimated as
(20)
π
k
i
(
t
)
=
1
n
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
.
Estimation for the average of FMM is
(21)
μ
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
z
k
j
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
i
=
1
,
…
,
g
k
(
t
)
.Considering the estimation for the covariance of FMM, because the covariance of target-originated measurement model is known or initialized, only covariance of clutter model needs to be estimated here:(22)
Σ
c
,
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
T
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
where i
=
1
,
…
,
g
c
,
k
(
t
).


### 4.3. Component Management Step

The component management step of model is completed by the technology of model merging and pruning [12].Notification. Models between clutter and target-originated measurement cannot merge each other.To clutter model, the process of merging and pruning is listed as follows:(a) Merging Step.Given the merging threshold U, let the set I
k
=
{
1
,
…
,
g
c
,
k
(
t
)
}, l
=
0, circulatel
=
l
+
1, j
=
argmax
i
∈
I
k
⁡
π
c
,
k
i; take a model i
∈
I
k, for example,ifd
k
i
j
(
t
)
≤
U, where
(23)
d
k
i
j
(
t
)
=
(
μ
c
,
k
i
(
t
)
-
μ
c
,
k
j
(
t
)
)
T
(
Σ
c
,
k
j
(
t
)
)
-
1
(
μ
c
,
k
i
(
t
)
-
μ
c
,
k
j
(
t
)
)
.Then leti
∈
L
k, and merge the models in set L
k, with the following merging formula:
(24)
π
~
c
,
k
l
(
t
+
1
)
=
∑
i
∈
L
k
π
c
,
k
i
(
t
)
,
μ
~
c
,
k
l
(
t
+
1
)
=
1
π
~
c
,
k
l
(
t
+
1
)
∑
i
∈
L
k
π
c
,
k
i
(
t
)
μ
c
,
k
i
(
t
+
1
)
,
Σ
~
c
,
k
l
(
t
+
1
)
=
1
π
~
c
,
k
l
(
t
+
1
)
×
∑
i
∈
L
k
π
c
,
k
i
(
t
)
(
Σ
c
,
k
i
(
t
)
+
μ
~
c
,
k
l
(
t
+
1
)
-
μ
c
,
k
i
(
t
)
)
×
(
μ
~
c
,
k
l
(
t
+
1
)
-
μ
c
,
k
i
(
t
)
)
T
.LetI
k
,
t
+
1
=
I
k
,
t
-
L
k
,
t; repeat this process until the set I
k
=
Φ; the merging ends.(b) Pruning Step.If the weight n
k
π
~
c
,
k
<
D
c
(
i
=
1
,
…
,
l
), where D
cdenotes the pruning threshold of clutter model, the corresponding model should be pruned.Then letl
=
l
-
1. Finally we can get the clutter model after merging and pruning, with the target number g
c
,
k
(
t
+
1
)
=
l.To target-originated measurement model, the process of merging and pruning is similar to the clutter model above. But the difference is that a target can not generate more measurements except one, the pruning gateD
t should be set far less than D
c in clutter model.Iteration among the three steps above circulates untilℒ
k
(
Y
k
;
ψ
k
(
t
+
1
)
)
-
ℒ
k
(
Y
k
;
ψ
k
(
t
)
) is under the convergence threshold ε of EM or MCMC algorithm. Then we can get ψ
^
k, the estimation of parameter value in FMM at the time k. Furthermore, the estimation of clutter model ℱ
^
c
,
k
(
·
;
ψ
^
c
,
k
), target number g
^
t
,
k and average value of target-originated measurement model μ
^
t
,
k
i
(
i
=
1
,
…
,
g
^
t
,
k
) are all obtained.FMM Parameters Estimation by EM approach over Section4 could be described by Algorithm 1.Algorithm 1:FMM Parameters Estimation by EM approach.
Input: the convergence threshold ε, the estimation of clutter model parameters at time t
k
-
1,
μ
c
,
k
i
(
0
)
=
μ
^
i
c
,
k
-
1
,
Σ
c
,
k
i
(
0
)
=
Σ
^
c
,
k
-
1
i
,
i
=
1
,
…
,
g
^
c
,
k
-
1.
Initialization step:
(A) Initialization of Clutter Model:
sample some clutter points randomly𝒫
k
=
{
p
k
i
}
i
=
1
g
0, the number of the clutter points is g
0,
For i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0, do,
μ
c
,
k
i
(
0
)
=
p
k
¯
,
Σ
c
,
k
i
(
0
)
=
σ
2
I
,
i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0
where σ
2
=
(
1
/
10
d
)
trace
(
(
1
/
n
k
)
∑
j
=
1
n
k
(
z
k
j
-
z
¯
k
)
(
z
k
j
-
z
¯
k
)
T
)
end for  i. g
c
,
k
(
0
)
=
g
^
c
,
k
-
1
+
g
0.
(B) Initialization of Target-Originated Measurement Model:
To survival targets,g
s
,
k
(
0
)
=
g
^
t
,
k
-
1, for i
=
1
,
…
,
g
s
,
k
(
0
), do
μ
s
,
k
i
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
),  Σ
s
,
k
i
(
0
)
=
R
k
,
(
i
=
1
,
…
,
g
s
,
k
(
0
)
)
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
=
∑
i
=
1
g
s
,
k
(
0
)
π
s
,
k
i
(
0
)
𝒩
(
·
;
μ
s
,
k
i
(
0
)
,
Σ
c
,
k
i
(
0
)
). end for i.
To Spontaneous Birth Targets,g
b
,
k
(
0
), μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)  will be set according to prior
information, fori
=
1
,
…
,
g
b
,
k
(
0
), do
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
=
∑
i
=
1
g
b
,
k
(
0
)
π
b
,
k
i
(
0
)
𝒩
(
·
;
μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)
). end for i.
To Spawned by Existent Targets,g
p
,
k
(
0
)
=
p
k
g
s
,
k
(
0
), for  i
=
1
,
…
,
g
s
,
k
(
0
), j
=
1
,
…
,
p
k, do
μ
p
,
k
i
,
j
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
)
+
d
p
.
k
j, ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
=
∑
j
=
1
p
k
∑
i
=
1
g
p
,
k
(
0
)
π
p
,
k
i
j
(
0
)
𝒩
(
·
;
μ
p
,
k
i
j
(
0
)
,
Σ
p
,
k
i
j
(
0
)
)
.
end for j; end for i.
g
t
,
k
(
0
)
=
g
s
,
k
(
0
)
+
g
b
,
k
(
0
)
+
g
p
,
k
(
0
).g
k
(
0
)
=
g
c
,
k
(
0
)
+
g
t
,
k
(
0
)
π
k
i
(
0
)
=
1
/
g
k
(
0
). Set t
:
=
0.
Repeat:
Expectation-step: calculate the conditional expectation of missing-data  E
k
=
{
e
k
1
,
…
,
e
k
n
k
}
For i
=
1
,
…
,
g
k, j
=
1
,
…
,
n
k, do
E
ψ
k
(
t
)
[
e
k
i
j
;
Z
k
]
=
τ
k
i
j
(
t
)
=
(
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
)
/
(
∑
i
=
1
g
k
(
t
)
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
). end for j; end for i.
calculate the conditional expectation of complete-data log likelihood given𝒫
k and ψ
k
(
t
)
Q
(
ψ
k
;
ψ
k
(
t
)
)
=
E
ψ
k
(
t
)
{
log
⁡
𝒞
k
(
Y
k
;
ψ
k
)
|
Z
k
}
=
∑
i
=
1
g
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
{
log
⁡
π
k
i
(
t
)
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
}
Maximization-step: require the global maximization of Q
(
ψ
k
;
ψ
k
(
t
)
) with respect to ψ
k over the
parameter space to give the updated estimate ψ
^
k
(
t
)
=
arg
⁡
ma
x
ψ
k
Q
(
ψ
k
;
ψ
k
(
t
)
).
for i
=
1
,
…
,
g
k, j
=
1
,
…
,
n
k, do
π
k
i
(
t
)
=
1
n
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
μ
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
z
k
j
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
i
=
1
,
…
,
g
k
(
t
),∑
c
,
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
T
∑
j
=
1
n
k
τ
k
i
j
(
t
)
end for j; end for i.
Component management step: manage the components of the set ψ
^
k
(
t
)
=
{
π
^
k
i
(
t
)
,
μ
^
k
i
(
t
)
,
Σ
^
k
i
(
t
)
}
i
=
1
g
k
according to the merging and pruning strategy described in IV.C.
Update g
k
(
t
+
1
) and let ψ
k
(
t
+
1
)
=
{
π
k
i
(
t
+
1
)
,
μ
k
i
(
t
+
1
)
,
Σ
k
i
(
t
+
1
)
}
i
=
1
g
k
(
t
+
1
). denote the  managed
component number and parameter set. Set t
:
=
t
+
1.
Until
|
ℒ
k
(
Y
k
;
ψ
k
(
t
+
1
)
)
-
ℒ
k
(
Y
k
;
ψ
k
(
t
)
)
|
<
ε
ℒ
k
(
Y
k
;
ψ
k
(
t
)
).
Output: the set of estimated parameters ψ
^
k
=
ψ
k
(
t
+
1
)


## 4.1. Parameters Initialization

Parameters in FMM should be well initialized as much as possible before estimation in EM algorithm. Initializations of clutter model and target-originated measurement model were discussed, respectively.


### 4.1.1. Initialization of Clutter Model

Considering that the clutter model could vary with time, the initialization process at timek is listed as follows.Inheritance of the value estimated at timet
k
-
1 is
(8)
μ
c
,
k
i
(
0
)
=
μ
^
c
,
k
-
1
i
,
Σ
c
,
k
i
(
0
)
=
Σ
^
c
,
k
-
1
i
,
i
=
1
,
…
,
g
^
c
,
k
-
1
.Then add some clutter points randomly; the number of the clutter points isg
0. The initialization average μ
c
,
k
i
(
0
) can be got from these chosen points in surveillance region S. The covariance matrix could be initialized as Σ
c
,
k
i
(
0
)
=
σ
2
I, i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0, where
(9)
σ
2
=
1
10
d
trace
(
1
n
k
∑
j
=
1
n
k
(
z
k
j
-
z
¯
k
)
(
z
k
j
-
z
¯
k
)
T
)
,
where z
¯
k
=
∑
j
=
1
n
k
z
¯
k
j
/
n
k denotes the average of measurement data.Consequently, the number of clutter models in initialization processg
c
,
k
(
0
)
=
g
^
c
,
k
-
1
+
g
0. In order to ensure the parameters of clutter models converging to the true value, the value of g
0 should be bigger than the number in reality.


### 4.1.2. Initialization of Target-Originated Measurement Model

At timek, targets totally are composed of survival targets and spontaneous birth targets and spawned by existent targets. So the number of target models in initialization process is
(10)
g
t
,
k
(
0
)
=
g
s
,
k
(
0
)
+
g
b
,
k
(
0
)
+
g
p
,
k
(
0
)
,
where g
s
,
k
(
0
) denotes the number of survival targets, g
b
,
k
(
0
) denotes the number of spontaneous birth targets, and g
p
,
k
(
0
) denotes the number of spawned by existent targets.The initialization of target-originated measurement model at timek can be represented as
(11)
ℱ
t
,
k
(
·
;
ψ
c
,
k
(
0
)
)
=
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
+
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
+
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
.(a) To Survival Targets. Measurement noise obeys Gaussian distribution, as shown in (2), so f
s
,
k
i
(
·
;
θ
s
,
k
i
)
~
𝒩
(
·
;
μ
s
,
k
i
,
Σ
s
,
k
i
). Furthermore, measurement model initialization of survival targets is
(12)
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
=
∑
i
=
1
g
s
,
k
(
0
)
π
s
,
k
i
(
0
)
𝒩
(
·
;
μ
s
,
k
i
(
0
)
,
Σ
c
,
k
i
(
0
)
)
.
Let g
s
,
k
(
0
)
=
g
^
t
,
k
-
1.The initialization averageμ
s
,
k
i
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
) could be got from prediction of target state. The covariance is known, the same with measurement noise of the sensor, Σ
s
,
k
i
(
0
)
=
R
k
,
i
=
1
,
…
,
g
s
,
k
(
0
).(b) To Spontaneous Birth Targets. As is mentioned in [11–13], assume that the initial position model obeys Gaussian distribution either; that is,
(13)
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
=
∑
i
=
1
g
b
,
k
(
0
)
π
b
,
k
i
(
0
)
𝒩
(
·
;
μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)
)
,
where parameters of model initialization g
b
,
k
(
0
), μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
) are all known, according to prior information.(c) To Spawned by Existent Targets. Supposing that each original target can create p
k new targets at most, model initialization can be represented as [11–13]
(14)
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
=
∑
j
=
1
p
k
∑
i
=
1
g
p
,
k
(
0
)
π
p
,
k
i
j
(
0
)
𝒩
(
·
;
μ
p
,
k
i
j
(
0
)
,
Σ
p
,
k
i
j
(
0
)
)
.The target number and state average can be assumed asg
p
,
k
(
0
)
=
p
k
g
s
,
k
(
0
), μ
p
,
k
i
,
j
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
)
+
d
p
.
k
j, where i
=
1
,
…
,
g
s
,
k
(
0
), j
=
1
,
…
,
p
k.The parametersd
p
.
k
j, Σ
p
,
k
i
j
(
0
) could be valued according to prior information similarly.The weights of models could be briefly treated as the same in initialization process:(15)
π
k
i
(
0
)
=
1
g
k
(
0
)
,
i
=
1
,
…
,
g
c
,
k
(
0
)
+
g
t
,
k
(
0
)
.


## 4.1.1. Initialization of Clutter Model

Considering that the clutter model could vary with time, the initialization process at timek is listed as follows.Inheritance of the value estimated at timet
k
-
1 is
(8)
μ
c
,
k
i
(
0
)
=
μ
^
c
,
k
-
1
i
,
Σ
c
,
k
i
(
0
)
=
Σ
^
c
,
k
-
1
i
,
i
=
1
,
…
,
g
^
c
,
k
-
1
.Then add some clutter points randomly; the number of the clutter points isg
0. The initialization average μ
c
,
k
i
(
0
) can be got from these chosen points in surveillance region S. The covariance matrix could be initialized as Σ
c
,
k
i
(
0
)
=
σ
2
I, i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0, where
(9)
σ
2
=
1
10
d
trace
(
1
n
k
∑
j
=
1
n
k
(
z
k
j
-
z
¯
k
)
(
z
k
j
-
z
¯
k
)
T
)
,
where z
¯
k
=
∑
j
=
1
n
k
z
¯
k
j
/
n
k denotes the average of measurement data.Consequently, the number of clutter models in initialization processg
c
,
k
(
0
)
=
g
^
c
,
k
-
1
+
g
0. In order to ensure the parameters of clutter models converging to the true value, the value of g
0 should be bigger than the number in reality.


## 4.1.2. Initialization of Target-Originated Measurement Model

At timek, targets totally are composed of survival targets and spontaneous birth targets and spawned by existent targets. So the number of target models in initialization process is
(10)
g
t
,
k
(
0
)
=
g
s
,
k
(
0
)
+
g
b
,
k
(
0
)
+
g
p
,
k
(
0
)
,
where g
s
,
k
(
0
) denotes the number of survival targets, g
b
,
k
(
0
) denotes the number of spontaneous birth targets, and g
p
,
k
(
0
) denotes the number of spawned by existent targets.The initialization of target-originated measurement model at timek can be represented as
(11)
ℱ
t
,
k
(
·
;
ψ
c
,
k
(
0
)
)
=
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
+
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
+
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
.(a) To Survival Targets. Measurement noise obeys Gaussian distribution, as shown in (2), so f
s
,
k
i
(
·
;
θ
s
,
k
i
)
~
𝒩
(
·
;
μ
s
,
k
i
,
Σ
s
,
k
i
). Furthermore, measurement model initialization of survival targets is
(12)
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
=
∑
i
=
1
g
s
,
k
(
0
)
π
s
,
k
i
(
0
)
𝒩
(
·
;
μ
s
,
k
i
(
0
)
,
Σ
c
,
k
i
(
0
)
)
.
Let g
s
,
k
(
0
)
=
g
^
t
,
k
-
1.The initialization averageμ
s
,
k
i
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
) could be got from prediction of target state. The covariance is known, the same with measurement noise of the sensor, Σ
s
,
k
i
(
0
)
=
R
k
,
i
=
1
,
…
,
g
s
,
k
(
0
).(b) To Spontaneous Birth Targets. As is mentioned in [11–13], assume that the initial position model obeys Gaussian distribution either; that is,
(13)
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
=
∑
i
=
1
g
b
,
k
(
0
)
π
b
,
k
i
(
0
)
𝒩
(
·
;
μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)
)
,
where parameters of model initialization g
b
,
k
(
0
), μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
) are all known, according to prior information.(c) To Spawned by Existent Targets. Supposing that each original target can create p
k new targets at most, model initialization can be represented as [11–13]
(14)
ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
=
∑
j
=
1
p
k
∑
i
=
1
g
p
,
k
(
0
)
π
p
,
k
i
j
(
0
)
𝒩
(
·
;
μ
p
,
k
i
j
(
0
)
,
Σ
p
,
k
i
j
(
0
)
)
.The target number and state average can be assumed asg
p
,
k
(
0
)
=
p
k
g
s
,
k
(
0
), μ
p
,
k
i
,
j
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
)
+
d
p
.
k
j, where i
=
1
,
…
,
g
s
,
k
(
0
), j
=
1
,
…
,
p
k.The parametersd
p
.
k
j, Σ
p
,
k
i
j
(
0
) could be valued according to prior information similarly.The weights of models could be briefly treated as the same in initialization process:(15)
π
k
i
(
0
)
=
1
g
k
(
0
)
,
i
=
1
,
…
,
g
c
,
k
(
0
)
+
g
t
,
k
(
0
)
.


## 4.2. Estimation of FMM Parameters

### 4.2.1. EM Approach

As FMM Parameters Estimation with EM approach, the loss variables can be treated as element labelingE
k
=
{
e
k
1
,
…
,
e
k
n
k
}, where e
k
1 is a vector of n
k dimension. The element e
k
i
j
=
1 or 0 illustrates whether z
k
j is originated from ith element of the FMM or not. So the complete dataset at time k is Y
k
=
{
Z
k
,
E
k
}, and logarithm likelihood function of the complete data is
(16)
log
⁡
𝒞
k
(
Y
k
;
ψ
k
)
=
∑
i
=
1
g
k
∑
j
=
1
n
k
e
k
i
j
{
log
⁡
⁡
π
k
i
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
)
}
.When the number of modelsg
k is known, EM algorithm can be iterated by E-step and M-step. But the number of models is unknown in this paper. Some criteria will be used for number estimation, such as minimum message length (MML) criterion and Bayesian information criterion (BIC) [10, 14]. Here in this paper, for the number of models in initialization being more than reality, the technology of model merging and pruning could be inserted into each step of iteration in EM algorithm. In this way, the number of FMM can be estimated. This method is more intuitive than those based on criteria. Now, this modified EM algorithm is described as follows:(a) E-Step.Conditional expectations of 𝒞
k
(
Y
k
;
ψ
k
) will be iterated. The expectations of loss data e
k
i
j at tth step is
(17)
E
ψ
k
(
t
)
[
e
k
i
j
;
Z
k
]
=
p
r
ψ
k
(
t
)
{
e
k
i
j
=
1
∣
Z
k
}
=
τ
k
i
(
z
k
j
;
ψ
k
(
t
)
)
=
τ
k
i
j
(
t
)
,
where τ
k
i
j
(
t
) represents the posterior probability of z
k
j belonging to model i. The formula for computing τ
k
i
j
(
t
) is
(18)
τ
k
i
j
(
t
)
=
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
∑
i
=
1
g
k
(
t
)
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
.According to the formula above, the conditional expectations oflog
⁡
𝒞
k
(
Y
k
;
ψ
k
) is
(19)
Q
(
ψ
k
;
ψ
k
(
t
)
)
=
E
ψ
k
(
t
)
{
log
⁡
⁡
𝒞
k
(
Y
k
;
ψ
k
)
∣
Z
k
}
=
∑
i
=
1
g
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
{
log
⁡
⁡
π
k
i
(
t
)
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
}
.In fact, the loss variable has reflected the association relationship between the measurement and target, that is, why association process can be avoided in this algorithm.(b) M-Step. The value of ψ
k will be estimated throught global maximum of Q
(
ψ
k
;
ψ
k
(
t
)
). By solving ∂
Q
k
(
ψ
k
,
ψ
k
(
t
)
)
/
∂
ψ
k
=
0, the weights of all models will be reestimated as
(20)
π
k
i
(
t
)
=
1
n
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
.
Estimation for the average of FMM is
(21)
μ
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
z
k
j
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
i
=
1
,
…
,
g
k
(
t
)
.Considering the estimation for the covariance of FMM, because the covariance of target-originated measurement model is known or initialized, only covariance of clutter model needs to be estimated here:(22)
Σ
c
,
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
T
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
where i
=
1
,
…
,
g
c
,
k
(
t
).


## 4.2.1. EM Approach

As FMM Parameters Estimation with EM approach, the loss variables can be treated as element labelingE
k
=
{
e
k
1
,
…
,
e
k
n
k
}, where e
k
1 is a vector of n
k dimension. The element e
k
i
j
=
1 or 0 illustrates whether z
k
j is originated from ith element of the FMM or not. So the complete dataset at time k is Y
k
=
{
Z
k
,
E
k
}, and logarithm likelihood function of the complete data is
(16)
log
⁡
𝒞
k
(
Y
k
;
ψ
k
)
=
∑
i
=
1
g
k
∑
j
=
1
n
k
e
k
i
j
{
log
⁡
⁡
π
k
i
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
)
}
.When the number of modelsg
k is known, EM algorithm can be iterated by E-step and M-step. But the number of models is unknown in this paper. Some criteria will be used for number estimation, such as minimum message length (MML) criterion and Bayesian information criterion (BIC) [10, 14]. Here in this paper, for the number of models in initialization being more than reality, the technology of model merging and pruning could be inserted into each step of iteration in EM algorithm. In this way, the number of FMM can be estimated. This method is more intuitive than those based on criteria. Now, this modified EM algorithm is described as follows:(a) E-Step.Conditional expectations of 𝒞
k
(
Y
k
;
ψ
k
) will be iterated. The expectations of loss data e
k
i
j at tth step is
(17)
E
ψ
k
(
t
)
[
e
k
i
j
;
Z
k
]
=
p
r
ψ
k
(
t
)
{
e
k
i
j
=
1
∣
Z
k
}
=
τ
k
i
(
z
k
j
;
ψ
k
(
t
)
)
=
τ
k
i
j
(
t
)
,
where τ
k
i
j
(
t
) represents the posterior probability of z
k
j belonging to model i. The formula for computing τ
k
i
j
(
t
) is
(18)
τ
k
i
j
(
t
)
=
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
∑
i
=
1
g
k
(
t
)
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
.According to the formula above, the conditional expectations oflog
⁡
𝒞
k
(
Y
k
;
ψ
k
) is
(19)
Q
(
ψ
k
;
ψ
k
(
t
)
)
=
E
ψ
k
(
t
)
{
log
⁡
⁡
𝒞
k
(
Y
k
;
ψ
k
)
∣
Z
k
}
=
∑
i
=
1
g
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
{
log
⁡
⁡
π
k
i
(
t
)
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
}
.In fact, the loss variable has reflected the association relationship between the measurement and target, that is, why association process can be avoided in this algorithm.(b) M-Step. The value of ψ
k will be estimated throught global maximum of Q
(
ψ
k
;
ψ
k
(
t
)
). By solving ∂
Q
k
(
ψ
k
,
ψ
k
(
t
)
)
/
∂
ψ
k
=
0, the weights of all models will be reestimated as
(20)
π
k
i
(
t
)
=
1
n
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
.
Estimation for the average of FMM is
(21)
μ
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
z
k
j
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
i
=
1
,
…
,
g
k
(
t
)
.Considering the estimation for the covariance of FMM, because the covariance of target-originated measurement model is known or initialized, only covariance of clutter model needs to be estimated here:(22)
Σ
c
,
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
T
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
where i
=
1
,
…
,
g
c
,
k
(
t
).


## 4.3. Component Management Step

The component management step of model is completed by the technology of model merging and pruning [12].Notification. Models between clutter and target-originated measurement cannot merge each other.To clutter model, the process of merging and pruning is listed as follows:(a) Merging Step.Given the merging threshold U, let the set I
k
=
{
1
,
…
,
g
c
,
k
(
t
)
}, l
=
0, circulatel
=
l
+
1, j
=
argmax
i
∈
I
k
⁡
π
c
,
k
i; take a model i
∈
I
k, for example,ifd
k
i
j
(
t
)
≤
U, where
(23)
d
k
i
j
(
t
)
=
(
μ
c
,
k
i
(
t
)
-
μ
c
,
k
j
(
t
)
)
T
(
Σ
c
,
k
j
(
t
)
)
-
1
(
μ
c
,
k
i
(
t
)
-
μ
c
,
k
j
(
t
)
)
.Then leti
∈
L
k, and merge the models in set L
k, with the following merging formula:
(24)
π
~
c
,
k
l
(
t
+
1
)
=
∑
i
∈
L
k
π
c
,
k
i
(
t
)
,
μ
~
c
,
k
l
(
t
+
1
)
=
1
π
~
c
,
k
l
(
t
+
1
)
∑
i
∈
L
k
π
c
,
k
i
(
t
)
μ
c
,
k
i
(
t
+
1
)
,
Σ
~
c
,
k
l
(
t
+
1
)
=
1
π
~
c
,
k
l
(
t
+
1
)
×
∑
i
∈
L
k
π
c
,
k
i
(
t
)
(
Σ
c
,
k
i
(
t
)
+
μ
~
c
,
k
l
(
t
+
1
)
-
μ
c
,
k
i
(
t
)
)
×
(
μ
~
c
,
k
l
(
t
+
1
)
-
μ
c
,
k
i
(
t
)
)
T
.LetI
k
,
t
+
1
=
I
k
,
t
-
L
k
,
t; repeat this process until the set I
k
=
Φ; the merging ends.(b) Pruning Step.If the weight n
k
π
~
c
,
k
<
D
c
(
i
=
1
,
…
,
l
), where D
cdenotes the pruning threshold of clutter model, the corresponding model should be pruned.Then letl
=
l
-
1. Finally we can get the clutter model after merging and pruning, with the target number g
c
,
k
(
t
+
1
)
=
l.To target-originated measurement model, the process of merging and pruning is similar to the clutter model above. But the difference is that a target can not generate more measurements except one, the pruning gateD
t should be set far less than D
c in clutter model.Iteration among the three steps above circulates untilℒ
k
(
Y
k
;
ψ
k
(
t
+
1
)
)
-
ℒ
k
(
Y
k
;
ψ
k
(
t
)
) is under the convergence threshold ε of EM or MCMC algorithm. Then we can get ψ
^
k, the estimation of parameter value in FMM at the time k. Furthermore, the estimation of clutter model ℱ
^
c
,
k
(
·
;
ψ
^
c
,
k
), target number g
^
t
,
k and average value of target-originated measurement model μ
^
t
,
k
i
(
i
=
1
,
…
,
g
^
t
,
k
) are all obtained.FMM Parameters Estimation by EM approach over Section4 could be described by Algorithm 1.Algorithm 1:FMM Parameters Estimation by EM approach.
Input: the convergence threshold ε, the estimation of clutter model parameters at time t
k
-
1,
μ
c
,
k
i
(
0
)
=
μ
^
i
c
,
k
-
1
,
Σ
c
,
k
i
(
0
)
=
Σ
^
c
,
k
-
1
i
,
i
=
1
,
…
,
g
^
c
,
k
-
1.
Initialization step:
(A) Initialization of Clutter Model:
sample some clutter points randomly𝒫
k
=
{
p
k
i
}
i
=
1
g
0, the number of the clutter points is g
0,
For i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0, do,
μ
c
,
k
i
(
0
)
=
p
k
¯
,
Σ
c
,
k
i
(
0
)
=
σ
2
I
,
i
=
g
^
c
,
k
-
1
+
1
,
…
,
g
^
c
,
k
-
1
+
g
0
where σ
2
=
(
1
/
10
d
)
trace
(
(
1
/
n
k
)
∑
j
=
1
n
k
(
z
k
j
-
z
¯
k
)
(
z
k
j
-
z
¯
k
)
T
)
end for  i. g
c
,
k
(
0
)
=
g
^
c
,
k
-
1
+
g
0.
(B) Initialization of Target-Originated Measurement Model:
To survival targets,g
s
,
k
(
0
)
=
g
^
t
,
k
-
1, for i
=
1
,
…
,
g
s
,
k
(
0
), do
μ
s
,
k
i
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
),  Σ
s
,
k
i
(
0
)
=
R
k
,
(
i
=
1
,
…
,
g
s
,
k
(
0
)
)
ℱ
s
,
k
(
·
;
ψ
s
,
k
(
0
)
)
=
∑
i
=
1
g
s
,
k
(
0
)
π
s
,
k
i
(
0
)
𝒩
(
·
;
μ
s
,
k
i
(
0
)
,
Σ
c
,
k
i
(
0
)
). end for i.
To Spontaneous Birth Targets,g
b
,
k
(
0
), μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)  will be set according to prior
information, fori
=
1
,
…
,
g
b
,
k
(
0
), do
ℱ
b
,
k
(
·
;
ψ
b
,
k
(
0
)
)
=
∑
i
=
1
g
b
,
k
(
0
)
π
b
,
k
i
(
0
)
𝒩
(
·
;
μ
b
,
k
i
(
0
)
,
Σ
b
,
k
i
(
0
)
). end for i.
To Spawned by Existent Targets,g
p
,
k
(
0
)
=
p
k
g
s
,
k
(
0
), for  i
=
1
,
…
,
g
s
,
k
(
0
), j
=
1
,
…
,
p
k, do
μ
p
,
k
i
,
j
(
0
)
=
h
k
(
x
^
k
|
k
-
1
i
)
+
d
p
.
k
j, ℱ
p
,
k
(
·
;
ψ
p
,
k
(
0
)
)
=
∑
j
=
1
p
k
∑
i
=
1
g
p
,
k
(
0
)
π
p
,
k
i
j
(
0
)
𝒩
(
·
;
μ
p
,
k
i
j
(
0
)
,
Σ
p
,
k
i
j
(
0
)
)
.
end for j; end for i.
g
t
,
k
(
0
)
=
g
s
,
k
(
0
)
+
g
b
,
k
(
0
)
+
g
p
,
k
(
0
).g
k
(
0
)
=
g
c
,
k
(
0
)
+
g
t
,
k
(
0
)
π
k
i
(
0
)
=
1
/
g
k
(
0
). Set t
:
=
0.
Repeat:
Expectation-step: calculate the conditional expectation of missing-data  E
k
=
{
e
k
1
,
…
,
e
k
n
k
}
For i
=
1
,
…
,
g
k, j
=
1
,
…
,
n
k, do
E
ψ
k
(
t
)
[
e
k
i
j
;
Z
k
]
=
τ
k
i
j
(
t
)
=
(
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
)
/
(
∑
i
=
1
g
k
(
t
)
π
k
i
(
t
)
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
). end for j; end for i.
calculate the conditional expectation of complete-data log likelihood given𝒫
k and ψ
k
(
t
)
Q
(
ψ
k
;
ψ
k
(
t
)
)
=
E
ψ
k
(
t
)
{
log
⁡
𝒞
k
(
Y
k
;
ψ
k
)
|
Z
k
}
=
∑
i
=
1
g
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
{
log
⁡
π
k
i
(
t
)
+
log
⁡
f
k
i
(
z
k
j
;
θ
k
i
(
t
)
)
}
Maximization-step: require the global maximization of Q
(
ψ
k
;
ψ
k
(
t
)
) with respect to ψ
k over the
parameter space to give the updated estimate ψ
^
k
(
t
)
=
arg
⁡
ma
x
ψ
k
Q
(
ψ
k
;
ψ
k
(
t
)
).
for i
=
1
,
…
,
g
k, j
=
1
,
…
,
n
k, do
π
k
i
(
t
)
=
1
n
k
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
μ
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
z
k
j
∑
j
=
1
n
k
τ
k
i
j
(
t
)
,
i
=
1
,
…
,
g
k
(
t
),∑
c
,
k
i
(
t
)
=
∑
j
=
1
n
k
τ
k
i
j
(
t
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
(
z
k
j
-
μ
c
,
k
i
(
t
)
)
T
∑
j
=
1
n
k
τ
k
i
j
(
t
)
end for j; end for i.
Component management step: manage the components of the set ψ
^
k
(
t
)
=
{
π
^
k
i
(
t
)
,
μ
^
k
i
(
t
)
,
Σ
^
k
i
(
t
)
}
i
=
1
g
k
according to the merging and pruning strategy described in IV.C.
Update g
k
(
t
+
1
) and let ψ
k
(
t
+
1
)
=
{
π
k
i
(
t
+
1
)
,
μ
k
i
(
t
+
1
)
,
Σ
k
i
(
t
+
1
)
}
i
=
1
g
k
(
t
+
1
). denote the  managed
component number and parameter set. Set t
:
=
t
+
1.
Until
|
ℒ
k
(
Y
k
;
ψ
k
(
t
+
1
)
)
-
ℒ
k
(
Y
k
;
ψ
k
(
t
)
)
|
<
ε
ℒ
k
(
Y
k
;
ψ
k
(
t
)
).
Output: the set of estimated parameters ψ
^
k
=
ψ
k
(
t
+
1
)


## 5. FMM Parameters Estimation by MCMC Approach

MCMC is an approach which obtains a Markov chain sampling from posterior distribution. We can extract information of posterior distribution through this Markov chain. Gibbs sampling method is one of MCMC approach. In order to solve the uncertainty of element number in mixture distribution, merging and pruning technology was consulted after each sampling step instead of RJMCMC method [15].


### 5.1. Posterior Distribution of FMM

Parameters to be estimated consist of the weight of modelsπ
j, the element labeling e
i
,
j, variance σ
j
2, and average value μ
j. In order to get posterior distribution in Bayesian frame, prior distribution of these parameters is necessary [10].Let the weight of modelsπ
j obey Dirichelet distribution:
(25)
(
π
1
,
⋯
π
m
)
~
𝒟
(
a
1
+
l
1
,
⋯
a
m
+
l
m
)
,
where a
j
>
0 is constant and l
j is the number of measurements which belong to model j.The element labelinge
i
,
j could be estimated by Bayes formula:
(26)
e
^
i
,
j
=
π
j
f
(
y
i
∣
θ
j
)
∑
j
=
1
k
π
j
f
(
y
i
∣
θ
j
)
,
l
j
=
∑
i
=
1
n
e
^
i
,
j
.Let varianceσ
j
2 obey Wishart distribution:
(27)
σ
j
2
~
Wishart
(
α
0
+
l
j
M
0
,
β
0
+
κ
j
2
N
0
)
κ
j
2
=
∑
i
=
1
n
(
y
i
-
μ
i
)
(
y
i
-
μ
i
)
T
·
e
^
i
,
j
∑
i
=
1
n
e
^
i
,
j
,
where α
0 and β
0 are positive constant, M
0 and N
0 are positive number.The average valueμ
j obey Gaussian distribution:
(28)
μ
j
~
𝒩
(
ξ
j
,
σ
j
2
)
,
where ξ
j
=
∑
i
=
1
n
y
i
·
e
^
i
,
j
/
∑
i
=
1
n
e
^
i
,
j.


### 5.2. Parameters Initialization

Easier than those of EM approach, in MCMC approach, initial values can be randomly selected in parameter space. While the number of models should be larger than expected value.


#### 5.2.1. Gibbs Sampling Method

Circulate(29)
e
i
,
j
(
t
)
=
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
∑
j
=
1
k
max
⁡
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
,
l
j
(
t
)
=
∑
i
=
1
n
e
^
i
,
j
(
t
)
,
μ
j
(
t
)
=
∑
i
=
1
n
y
i
·
e
i
,
j
(
t
)
l
j
(
t
)
,
ξ
j
(
t
)
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
·
e
i
,
j
(
t
)
,
κ
j
(
t
)
2
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
(
y
i
-
μ
j
(
t
)
)
T
·
e
i
,
j
(
t
)
∑
i
=
1
n
e
i
,
j
(
t
)
,
σ
j
(
t
)
2
~
Wishart
(
α
0
+
l
j
(
t
)
M
0
,
β
0
+
κ
j
(
t
)
2
N
0
)
,
μ
j
(
t
)
~
𝒩
(
ξ
j
(
t
)
,
σ
j
(
t
)
2
)
,
(
π
1
(
t
)
,
…
,
π
m
(
t
)
)
~
𝒟
(
a
1
+
l
1
(
t
)
,
…
,
a
m
+
l
m
(
t
)
)
.


### 5.3. Management of Model Number

Choosing proper thresholdsλ
π, λ
μ, and λ
σ, use merging and pruning technology for reference from EM approach. In this way, algorithm becomes more concise ignoring jumping in parameter space.


## 5.1. Posterior Distribution of FMM

Parameters to be estimated consist of the weight of modelsπ
j, the element labeling e
i
,
j, variance σ
j
2, and average value μ
j. In order to get posterior distribution in Bayesian frame, prior distribution of these parameters is necessary [10].Let the weight of modelsπ
j obey Dirichelet distribution:
(25)
(
π
1
,
⋯
π
m
)
~
𝒟
(
a
1
+
l
1
,
⋯
a
m
+
l
m
)
,
where a
j
>
0 is constant and l
j is the number of measurements which belong to model j.The element labelinge
i
,
j could be estimated by Bayes formula:
(26)
e
^
i
,
j
=
π
j
f
(
y
i
∣
θ
j
)
∑
j
=
1
k
π
j
f
(
y
i
∣
θ
j
)
,
l
j
=
∑
i
=
1
n
e
^
i
,
j
.Let varianceσ
j
2 obey Wishart distribution:
(27)
σ
j
2
~
Wishart
(
α
0
+
l
j
M
0
,
β
0
+
κ
j
2
N
0
)
κ
j
2
=
∑
i
=
1
n
(
y
i
-
μ
i
)
(
y
i
-
μ
i
)
T
·
e
^
i
,
j
∑
i
=
1
n
e
^
i
,
j
,
where α
0 and β
0 are positive constant, M
0 and N
0 are positive number.The average valueμ
j obey Gaussian distribution:
(28)
μ
j
~
𝒩
(
ξ
j
,
σ
j
2
)
,
where ξ
j
=
∑
i
=
1
n
y
i
·
e
^
i
,
j
/
∑
i
=
1
n
e
^
i
,
j.


## 5.2. Parameters Initialization

Easier than those of EM approach, in MCMC approach, initial values can be randomly selected in parameter space. While the number of models should be larger than expected value.


### 5.2.1. Gibbs Sampling Method

Circulate(29)
e
i
,
j
(
t
)
=
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
∑
j
=
1
k
max
⁡
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
,
l
j
(
t
)
=
∑
i
=
1
n
e
^
i
,
j
(
t
)
,
μ
j
(
t
)
=
∑
i
=
1
n
y
i
·
e
i
,
j
(
t
)
l
j
(
t
)
,
ξ
j
(
t
)
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
·
e
i
,
j
(
t
)
,
κ
j
(
t
)
2
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
(
y
i
-
μ
j
(
t
)
)
T
·
e
i
,
j
(
t
)
∑
i
=
1
n
e
i
,
j
(
t
)
,
σ
j
(
t
)
2
~
Wishart
(
α
0
+
l
j
(
t
)
M
0
,
β
0
+
κ
j
(
t
)
2
N
0
)
,
μ
j
(
t
)
~
𝒩
(
ξ
j
(
t
)
,
σ
j
(
t
)
2
)
,
(
π
1
(
t
)
,
…
,
π
m
(
t
)
)
~
𝒟
(
a
1
+
l
1
(
t
)
,
…
,
a
m
+
l
m
(
t
)
)
.


## 5.2.1. Gibbs Sampling Method

Circulate(29)
e
i
,
j
(
t
)
=
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
∑
j
=
1
k
max
⁡
π
i
j
p
(
y
i
∣
θ
(
t
-
1
)
)
,
l
j
(
t
)
=
∑
i
=
1
n
e
^
i
,
j
(
t
)
,
μ
j
(
t
)
=
∑
i
=
1
n
y
i
·
e
i
,
j
(
t
)
l
j
(
t
)
,
ξ
j
(
t
)
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
·
e
i
,
j
(
t
)
,
κ
j
(
t
)
2
=
∑
i
=
1
n
(
y
i
-
μ
j
(
t
)
)
(
y
i
-
μ
j
(
t
)
)
T
·
e
i
,
j
(
t
)
∑
i
=
1
n
e
i
,
j
(
t
)
,
σ
j
(
t
)
2
~
Wishart
(
α
0
+
l
j
(
t
)
M
0
,
β
0
+
κ
j
(
t
)
2
N
0
)
,
μ
j
(
t
)
~
𝒩
(
ξ
j
(
t
)
,
σ
j
(
t
)
2
)
,
(
π
1
(
t
)
,
…
,
π
m
(
t
)
)
~
𝒟
(
a
1
+
l
1
(
t
)
,
…
,
a
m
+
l
m
(
t
)
)
.


## 5.3. Management of Model Number

Choosing proper thresholdsλ
π, λ
μ, and λ
σ, use merging and pruning technology for reference from EM approach. In this way, algorithm becomes more concise ignoring jumping in parameter space.


## 6. Multitarget State Estimation

### 6.1. Equivalent Measurement of Target

The real measurement originated from target is hardly obtained because of the effect from clutter in MTT. In many circumstances, the real measurement is always replaced by equivalent measurement of target. Taking JPDA algorithm for example, the equivalent measurement could be obtained by probability weighted moments from actual measurement in the gate. To the algorithm proposed in this paper, association probability of target and measurement is represented by so-called loss variable, so average value of target-originated measurement model, that is,μ
^
t
,
k
i
(
i
=
1
,
…
,
g
^
t
,
k
), represents the equivalent measurement of target i, with variance R
k still. And the so-called gate is extended to the global surveillance region.


### 6.2. Multitarget State Estimation

Suppose that the state and measurement functions are both of the linear systems. That is,(30)
f
k
i
(
x
k
i
)
=
F
k
i
,
h
k
i
(
x
k
i
)
=
H
k
i
(
x
k
i
)
,
i
=
1
,
…
,
t
k
.By substitution of the equivalent measurementμ
^
t
,
k
iinto Kalman filter, state estimation of target i is acquired. The algorithm of Kalman filter is(1)
Prediction Step(31)
x
^
k
∣
k
-
1
i
=
F
k
-
1
i
x
^
k
-
1
i
,
P
k
∣
k
-
1
i
=
Q
k
-
1
i
+
F
k
-
1
i
P
k
-
1
i
(
F
k
-
1
i
)
T
.
(2)
Update Step(32)
K
k
i
=
P
k
∣
k
-
1
i
H
k
T
(
H
k
P
k
∣
k
-
1
i
H
k
T
)
-
1
,
x
^
k
i
=
x
^
k
∣
k
-
1
i
+
K
k
i
(
μ
^
c
,
k
i
-
H
k
x
^
k
∣
k
-
1
i
)
,
P
k
i
=
(
I
-
K
k
i
H
k
)
P
k
∣
k
-
1
i
.If the state function or measurement function cannot meet linear condition, nonlinear filter methods such as extended Kalman filter (EKF) and unscented filter (UF) [16] will make contribution.


## 6.1. Equivalent Measurement of Target

The real measurement originated from target is hardly obtained because of the effect from clutter in MTT. In many circumstances, the real measurement is always replaced by equivalent measurement of target. Taking JPDA algorithm for example, the equivalent measurement could be obtained by probability weighted moments from actual measurement in the gate. To the algorithm proposed in this paper, association probability of target and measurement is represented by so-called loss variable, so average value of target-originated measurement model, that is,μ
^
t
,
k
i
(
i
=
1
,
…
,
g
^
t
,
k
), represents the equivalent measurement of target i, with variance R
k still. And the so-called gate is extended to the global surveillance region.


## 6.2. Multitarget State Estimation

Suppose that the state and measurement functions are both of the linear systems. That is,(30)
f
k
i
(
x
k
i
)
=
F
k
i
,
h
k
i
(
x
k
i
)
=
H
k
i
(
x
k
i
)
,
i
=
1
,
…
,
t
k
.By substitution of the equivalent measurementμ
^
t
,
k
iinto Kalman filter, state estimation of target i is acquired. The algorithm of Kalman filter is(1)
Prediction Step(31)
x
^
k
∣
k
-
1
i
=
F
k
-
1
i
x
^
k
-
1
i
,
P
k
∣
k
-
1
i
=
Q
k
-
1
i
+
F
k
-
1
i
P
k
-
1
i
(
F
k
-
1
i
)
T
.
(2)
Update Step(32)
K
k
i
=
P
k
∣
k
-
1
i
H
k
T
(
H
k
P
k
∣
k
-
1
i
H
k
T
)
-
1
,
x
^
k
i
=
x
^
k
∣
k
-
1
i
+
K
k
i
(
μ
^
c
,
k
i
-
H
k
x
^
k
∣
k
-
1
i
)
,
P
k
i
=
(
I
-
K
k
i
H
k
)
P
k
∣
k
-
1
i
.If the state function or measurement function cannot meet linear condition, nonlinear filter methods such as extended Kalman filter (EKF) and unscented filter (UF) [16] will make contribution.


## 7. Simulation

### 7.1. Scene Generation

The motion model of the targets can be described as(33)
x
^
k
i
=
F
k
-
1
i
x
^
k
-
1
i
+
ω
k
-
1
i
,
i
=
1
,
…
,
t
k
,
where t
k denotes the target number at time k, x
k
i
=
[
x
k
i
,
x
˙
k
i
,
x
¨
k
i
,
y
k
i
,
y
˙
k
i
,
y
¨
k
i
] denotes the state vector of target i, F
k
i denotes the state transition matrix, and ω
k
i
~
𝒩
(
0
,
Q
k
i
) denotes the process noise. The tracking process lasted for 100 s with a sampling interval 1 s. In two-dimensional coordinate system, the surveillance region is [
-
1000
m
,
1000
m
]
×
[
-
1000
m
,
1000
m
]. The flight track is shown in Figure 1.Figure 1
Tracks of targets.In Figure1, the circle denotes starting point, while the rectangular denotes the end points of a track. The solid line is the track of target 1, which was born at time 1 and disappeared at time 70. The dotted line is the track of target 2, which was born at time 20 and disappeared at time 90. The dashed line is the track of target 3, which was spawned by target 2 and disappeared at time 100.Target 1 and target 2 moved with constant acceleration (CA). The state transition matrixF
k
i and variance matrix of process noise Q
k
i [17] are
(34)
F
k
i
=
[
F
CA
F
CA
]
,
Q
k
i
=
[
Q
CA
Q
CA
]
,
where
(35)
F
CA
=
[
1
T
T
2
2
1
T
1
]
,
Q
CA
=
σ
ω
2
[
T
4
4
T
3
2
T
2
2
T
3
2
T
2
2
T
T
2
2
T
1
]
.Target 3 moved with constant velocity (CV). The state transition matrixF
k
3 and variance matrix of process noiseQ
k
3 are
(36)
F
k
3
=
[
F
CV
F
CV
]
,
Q
k
3
=
[
Q
CV
Q
CV
]
,
where
(37)
F
CV
=
[
1
T
0
1
0
0
]
,
Q
CV
=
σ
ω
2
[
T
4
4
T
3
2
0
T
3
2
T
2
0
0
0
0
]
,
where  σ
ω denotes the standard deviation of process noise, σ
ω
=
0.01
m
/
s
2.Supposing that the measurement function is linear,(38)
z
k
=
[
1
0
0
0
0
0
0
0
0
1
0
0
]
x
k
+
υ
k
.The sensor was in the origin of coordinate, with detection probabilityp
D
=
0.98 and measurement noise υ
k
~
𝒩
(
0
,
R
k
), where variance matrix of measurement noise R
k
=
σ
υ
2
I,  σ
υ
=
12.5
m.Assume that the clutter model obeys stable Poisson distribution in this tracking process, with the number of cluttersN
c and parameter λ
c
=
50. 								(39)
P
(
N
c
=
c
k
)
=
e
-
λ
c
λ
c
c
k
c
k
!
,
where c
k is the number of clutters and λ
c
=
50 denotes the average level that the sensor could receive 50 clutters each frame.Assume that the position distribution of clutter model is(40)
ℱ
c
,
k
(
·
;
θ
c
)
=
π
c
1
𝒰
(
·
)
+
∑
i
=
2
4
π
c
i
𝒟
(
·
)
,
∑
i
=
1
4
π
c
i
=
1
,
which is composed of one uniform distribution dispersed over the whole surveillance region and three groups of complex distribution concentrated in elliptic areas. Each group is a superposition of different types clutter model, not limited to uniform and Gaussian distribution, see in Figure 2.Figure 2
Fitted clutter model.Elliptic range of this complex clutter is listed as in Table1.Table 1
Elliptic range of the complex clutter.
Model
Weight
Center(×
1
0
2
m)
Long axis(×
1
0
2
m)
Short axis (×
1
0
2
m)
1
0.3
/
/
2
0.3
[4
,
2]
3
2
3
0.2
[−2, 4]
2
2
4
0.2
[2, −2]
2
2


### 7.2. Parameter Estimation in FMM

The effect of algorithm proposed in this paper is compared with Gaussian mixture PHD (GM-PHD) filter, which directly estimates the number and state of multitarget without clutter model fitting.To spontaneous birth targets, the measurement model is represented by (13), with g
b
,
k
(
0
)
=
2, μ
b
,
k
1
(
0
)
=
[
200,800
]
T, μ
b
,
k
2
(
0
)
=
[
800
,
-
800
]
T, and Σ
b
,
k
1
(
0
)
=
Σ
b
,
k
2
(
0
)
=
diag
⁡
(
[
100,100
]
). To the spawned targets generated from existent targets, the measurement model is represented by (14), with p
k
=
1, d
p
.
k
1
=
[
0,0
]
T, and Σ
p
,
k
i
1
=
diag
⁡
(
[
100,100
]
), i
=
1
,
…
,
t
k.Let the merging threshold of the modelsU
=
4, pruning threshold of clutter model D
c
=
3, and pruning threshold of target-originated measurement model D
t
=
0.5. With the algorithm proposed, estimation of clutter model can approximately converge to real model, as is shown in Figure 2.In Figure2, the solid line illustrates 95% ellipse of the real clutter model, while the dotted line illustrates 95% ellipse of the estimation. Times of Monte Carlo simulations have shown that some estimation results could fit the region which the clutter concentrated according to Table 1, shown as Figure 2. But it does not exclude some results shown as Figure 3, caused by the complexity of clutter model. No matter how complex the clutter is, FMM algorithm would fit it as a linear superposition of one uniform and many Gaussian distributions, but the parameters may differ from each other, as shown by different dotted ellipse in Figure 3. Nevertheless, these situations do not affect the validity of FMM algorithm, with estimation accuracy satisfied.Other estimation result of complex clutter model.
(a)
(b)
(c)
(d)


#### 7.2.1. Effect Comparison of Target Number Estimation

The target number estimated by GM-PHD filter is represented in Figure4, while the algorithm proposed in Figure 5.Figure 4
Target number estimation of GM-PHD filter.Figure 5
Target number estimation of the proposed method.In Figures4 and 5, the solid line illustrates the real target number varying with time, while it circles estimation. Apparently, effect of GM-PHD filter is worse. In circumstances of unknown complex clutter distribution, the assumption of uniform style will lead to serious error in target number estimation. Generally speaking, the more dense clutters the target passing and more time spent in this area, the more fake targets could emerge from estimation. As the algorithm proposed can fit clutter model synchronous with tracking, fake targets in high-density clutter area will decrease. Consequently, estimation accuracy of target number improved apparently. One attention: When target is passing high-density clutter area, the state estimation will get lost provisionally.


#### 7.2.2. Effect Comparison of Target Position Estimation

Position estimation of multitarget by GM-PHD filter is represented in Figure6, while the method is proposed in Figure 7.Figure 6
Target position estimation of GM-PHD filter.Figure 7
Target position estimation of the proposed method.In Figures6 and 7, the solid line illustrates the real target track, while circles illustrate the position estimated. Similarly, the assumption of uniform style will lead to a crowd of fake targets emerging in high-density clutter area.


#### 7.2.3. Evaluation of Effects

Different from single target tracking, root mean square error (RMSE) cannot measure error in multiple targets tracking quantitatively, [18] suggesting that Wasserstein distance can describe error in case of the number of target varying with time in MTT.Suppose thatX
k
=
{
x
k
1
,
…
,
x
k
T
k
} represents actual multitarget state set at time k; estimated state set X
^
k
=
{
x
^
k
1
,
…
,
x
^
k
T
^
k
}, where T
k and T
^
k denote the actual and estimated number of multitarget, respectively. The Wasserstein distance can be defined as
(41)
d
p
,
k
(
X
^
k
,
X
k
)
=
min
⁡
C
k
⁡
∑
i
=
1
|
X
^
k
|
∑
j
=
1
|
X
k
|
C
k
i
j
∥
x
^
k
i
-
x
k
j
∥
p
p
,
where C
k denotes transfer matrix with every element C
k
i
j
≥
0, Σ
j
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
^
|, and Σ
i
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
|. Also, |
·
| represents cardinality and ∥
·
∥
p represents norm with p
=
2. When the set of X
k or X
^
k is empty, Wasserstein distance can be assigned to 0. When the number of elements in X
k and X
^
k is the same, Wasserstein distance is the best association distance.Wasserstein distance varying with time of these two methods is shown in Figures8 and 9. To punish the mistake in target number estimation, Wasserstein distance reaches peak value when the number is wrongly estimated. From Figures 8 and 9, Wasserstein distance of the latter is better than that of former, because the latter has more precise estimation. When the target number is correctly estimated, the value of Wasserstein distance is about 10 m, approximately the square root of trace
(
Σ
b
,
k
i
(
0
)
) in Section 7.2.Figure 8
Wasserstein distance of GM-PHD filter.Figure 9
Wasserstein distance of the proposed method.


## 7.1. Scene Generation

The motion model of the targets can be described as(33)
x
^
k
i
=
F
k
-
1
i
x
^
k
-
1
i
+
ω
k
-
1
i
,
i
=
1
,
…
,
t
k
,
where t
k denotes the target number at time k, x
k
i
=
[
x
k
i
,
x
˙
k
i
,
x
¨
k
i
,
y
k
i
,
y
˙
k
i
,
y
¨
k
i
] denotes the state vector of target i, F
k
i denotes the state transition matrix, and ω
k
i
~
𝒩
(
0
,
Q
k
i
) denotes the process noise. The tracking process lasted for 100 s with a sampling interval 1 s. In two-dimensional coordinate system, the surveillance region is [
-
1000
m
,
1000
m
]
×
[
-
1000
m
,
1000
m
]. The flight track is shown in Figure 1.Figure 1
Tracks of targets.In Figure1, the circle denotes starting point, while the rectangular denotes the end points of a track. The solid line is the track of target 1, which was born at time 1 and disappeared at time 70. The dotted line is the track of target 2, which was born at time 20 and disappeared at time 90. The dashed line is the track of target 3, which was spawned by target 2 and disappeared at time 100.Target 1 and target 2 moved with constant acceleration (CA). The state transition matrixF
k
i and variance matrix of process noise Q
k
i [17] are
(34)
F
k
i
=
[
F
CA
F
CA
]
,
Q
k
i
=
[
Q
CA
Q
CA
]
,
where
(35)
F
CA
=
[
1
T
T
2
2
1
T
1
]
,
Q
CA
=
σ
ω
2
[
T
4
4
T
3
2
T
2
2
T
3
2
T
2
2
T
T
2
2
T
1
]
.Target 3 moved with constant velocity (CV). The state transition matrixF
k
3 and variance matrix of process noiseQ
k
3 are
(36)
F
k
3
=
[
F
CV
F
CV
]
,
Q
k
3
=
[
Q
CV
Q
CV
]
,
where
(37)
F
CV
=
[
1
T
0
1
0
0
]
,
Q
CV
=
σ
ω
2
[
T
4
4
T
3
2
0
T
3
2
T
2
0
0
0
0
]
,
where  σ
ω denotes the standard deviation of process noise, σ
ω
=
0.01
m
/
s
2.Supposing that the measurement function is linear,(38)
z
k
=
[
1
0
0
0
0
0
0
0
0
1
0
0
]
x
k
+
υ
k
.The sensor was in the origin of coordinate, with detection probabilityp
D
=
0.98 and measurement noise υ
k
~
𝒩
(
0
,
R
k
), where variance matrix of measurement noise R
k
=
σ
υ
2
I,  σ
υ
=
12.5
m.Assume that the clutter model obeys stable Poisson distribution in this tracking process, with the number of cluttersN
c and parameter λ
c
=
50. 								(39)
P
(
N
c
=
c
k
)
=
e
-
λ
c
λ
c
c
k
c
k
!
,
where c
k is the number of clutters and λ
c
=
50 denotes the average level that the sensor could receive 50 clutters each frame.Assume that the position distribution of clutter model is(40)
ℱ
c
,
k
(
·
;
θ
c
)
=
π
c
1
𝒰
(
·
)
+
∑
i
=
2
4
π
c
i
𝒟
(
·
)
,
∑
i
=
1
4
π
c
i
=
1
,
which is composed of one uniform distribution dispersed over the whole surveillance region and three groups of complex distribution concentrated in elliptic areas. Each group is a superposition of different types clutter model, not limited to uniform and Gaussian distribution, see in Figure 2.Figure 2
Fitted clutter model.Elliptic range of this complex clutter is listed as in Table1.Table 1
Elliptic range of the complex clutter.
Model
Weight
Center(×
1
0
2
m)
Long axis(×
1
0
2
m)
Short axis (×
1
0
2
m)
1
0.3
/
/
2
0.3
[4
,
2]
3
2
3
0.2
[−2, 4]
2
2
4
0.2
[2, −2]
2
2


## 7.2. Parameter Estimation in FMM

The effect of algorithm proposed in this paper is compared with Gaussian mixture PHD (GM-PHD) filter, which directly estimates the number and state of multitarget without clutter model fitting.To spontaneous birth targets, the measurement model is represented by (13), with g
b
,
k
(
0
)
=
2, μ
b
,
k
1
(
0
)
=
[
200,800
]
T, μ
b
,
k
2
(
0
)
=
[
800
,
-
800
]
T, and Σ
b
,
k
1
(
0
)
=
Σ
b
,
k
2
(
0
)
=
diag
⁡
(
[
100,100
]
). To the spawned targets generated from existent targets, the measurement model is represented by (14), with p
k
=
1, d
p
.
k
1
=
[
0,0
]
T, and Σ
p
,
k
i
1
=
diag
⁡
(
[
100,100
]
), i
=
1
,
…
,
t
k.Let the merging threshold of the modelsU
=
4, pruning threshold of clutter model D
c
=
3, and pruning threshold of target-originated measurement model D
t
=
0.5. With the algorithm proposed, estimation of clutter model can approximately converge to real model, as is shown in Figure 2.In Figure2, the solid line illustrates 95% ellipse of the real clutter model, while the dotted line illustrates 95% ellipse of the estimation. Times of Monte Carlo simulations have shown that some estimation results could fit the region which the clutter concentrated according to Table 1, shown as Figure 2. But it does not exclude some results shown as Figure 3, caused by the complexity of clutter model. No matter how complex the clutter is, FMM algorithm would fit it as a linear superposition of one uniform and many Gaussian distributions, but the parameters may differ from each other, as shown by different dotted ellipse in Figure 3. Nevertheless, these situations do not affect the validity of FMM algorithm, with estimation accuracy satisfied.Other estimation result of complex clutter model.
(a)
(b)
(c)
(d)


### 7.2.1. Effect Comparison of Target Number Estimation

The target number estimated by GM-PHD filter is represented in Figure4, while the algorithm proposed in Figure 5.Figure 4
Target number estimation of GM-PHD filter.Figure 5
Target number estimation of the proposed method.In Figures4 and 5, the solid line illustrates the real target number varying with time, while it circles estimation. Apparently, effect of GM-PHD filter is worse. In circumstances of unknown complex clutter distribution, the assumption of uniform style will lead to serious error in target number estimation. Generally speaking, the more dense clutters the target passing and more time spent in this area, the more fake targets could emerge from estimation. As the algorithm proposed can fit clutter model synchronous with tracking, fake targets in high-density clutter area will decrease. Consequently, estimation accuracy of target number improved apparently. One attention: When target is passing high-density clutter area, the state estimation will get lost provisionally.


### 7.2.2. Effect Comparison of Target Position Estimation

Position estimation of multitarget by GM-PHD filter is represented in Figure6, while the method is proposed in Figure 7.Figure 6
Target position estimation of GM-PHD filter.Figure 7
Target position estimation of the proposed method.In Figures6 and 7, the solid line illustrates the real target track, while circles illustrate the position estimated. Similarly, the assumption of uniform style will lead to a crowd of fake targets emerging in high-density clutter area.


### 7.2.3. Evaluation of Effects

Different from single target tracking, root mean square error (RMSE) cannot measure error in multiple targets tracking quantitatively, [18] suggesting that Wasserstein distance can describe error in case of the number of target varying with time in MTT.Suppose thatX
k
=
{
x
k
1
,
…
,
x
k
T
k
} represents actual multitarget state set at time k; estimated state set X
^
k
=
{
x
^
k
1
,
…
,
x
^
k
T
^
k
}, where T
k and T
^
k denote the actual and estimated number of multitarget, respectively. The Wasserstein distance can be defined as
(41)
d
p
,
k
(
X
^
k
,
X
k
)
=
min
⁡
C
k
⁡
∑
i
=
1
|
X
^
k
|
∑
j
=
1
|
X
k
|
C
k
i
j
∥
x
^
k
i
-
x
k
j
∥
p
p
,
where C
k denotes transfer matrix with every element C
k
i
j
≥
0, Σ
j
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
^
|, and Σ
i
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
|. Also, |
·
| represents cardinality and ∥
·
∥
p represents norm with p
=
2. When the set of X
k or X
^
k is empty, Wasserstein distance can be assigned to 0. When the number of elements in X
k and X
^
k is the same, Wasserstein distance is the best association distance.Wasserstein distance varying with time of these two methods is shown in Figures8 and 9. To punish the mistake in target number estimation, Wasserstein distance reaches peak value when the number is wrongly estimated. From Figures 8 and 9, Wasserstein distance of the latter is better than that of former, because the latter has more precise estimation. When the target number is correctly estimated, the value of Wasserstein distance is about 10 m, approximately the square root of trace
(
Σ
b
,
k
i
(
0
)
) in Section 7.2.Figure 8
Wasserstein distance of GM-PHD filter.Figure 9
Wasserstein distance of the proposed method.


## 7.2.1. Effect Comparison of Target Number Estimation

The target number estimated by GM-PHD filter is represented in Figure4, while the algorithm proposed in Figure 5.Figure 4
Target number estimation of GM-PHD filter.Figure 5
Target number estimation of the proposed method.In Figures4 and 5, the solid line illustrates the real target number varying with time, while it circles estimation. Apparently, effect of GM-PHD filter is worse. In circumstances of unknown complex clutter distribution, the assumption of uniform style will lead to serious error in target number estimation. Generally speaking, the more dense clutters the target passing and more time spent in this area, the more fake targets could emerge from estimation. As the algorithm proposed can fit clutter model synchronous with tracking, fake targets in high-density clutter area will decrease. Consequently, estimation accuracy of target number improved apparently. One attention: When target is passing high-density clutter area, the state estimation will get lost provisionally.


## 7.2.2. Effect Comparison of Target Position Estimation

Position estimation of multitarget by GM-PHD filter is represented in Figure6, while the method is proposed in Figure 7.Figure 6
Target position estimation of GM-PHD filter.Figure 7
Target position estimation of the proposed method.In Figures6 and 7, the solid line illustrates the real target track, while circles illustrate the position estimated. Similarly, the assumption of uniform style will lead to a crowd of fake targets emerging in high-density clutter area.


## 7.2.3. Evaluation of Effects

Different from single target tracking, root mean square error (RMSE) cannot measure error in multiple targets tracking quantitatively, [18] suggesting that Wasserstein distance can describe error in case of the number of target varying with time in MTT.Suppose thatX
k
=
{
x
k
1
,
…
,
x
k
T
k
} represents actual multitarget state set at time k; estimated state set X
^
k
=
{
x
^
k
1
,
…
,
x
^
k
T
^
k
}, where T
k and T
^
k denote the actual and estimated number of multitarget, respectively. The Wasserstein distance can be defined as
(41)
d
p
,
k
(
X
^
k
,
X
k
)
=
min
⁡
C
k
⁡
∑
i
=
1
|
X
^
k
|
∑
j
=
1
|
X
k
|
C
k
i
j
∥
x
^
k
i
-
x
k
j
∥
p
p
,
where C
k denotes transfer matrix with every element C
k
i
j
≥
0, Σ
j
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
^
|, and Σ
i
=
1
|
X
k
|
C
k
i
j
=
1
/
|
X
|. Also, |
·
| represents cardinality and ∥
·
∥
p represents norm with p
=
2. When the set of X
k or X
^
k is empty, Wasserstein distance can be assigned to 0. When the number of elements in X
k and X
^
k is the same, Wasserstein distance is the best association distance.Wasserstein distance varying with time of these two methods is shown in Figures8 and 9. To punish the mistake in target number estimation, Wasserstein distance reaches peak value when the number is wrongly estimated. From Figures 8 and 9, Wasserstein distance of the latter is better than that of former, because the latter has more precise estimation. When the target number is correctly estimated, the value of Wasserstein distance is about 10 m, approximately the square root of trace
(
Σ
b
,
k
i
(
0
)
) in Section 7.2.Figure 8
Wasserstein distance of GM-PHD filter.Figure 9
Wasserstein distance of the proposed method.


## 8. Conclusion

FMM approach can solve unknown clutter problem in MTT. A novel MTT algorithm based on clutter model preestimation is put forward in this paper. In this algorithm, multitarget likelihood function is established with the finite mixture model (FMM), the parameters of which can be estimated by the algorithm of EM and MCMC. These two algorithms were put forward in the paper. Furthermore, target number and multitarget states can be estimated as well as the clutter model fitted. No matter how complex the clutter is, FMM algorithm would fit it as a linear superposition of one uniform and many Gaussian distributions, but the parameters may differ from each other, but it will not affect the validity of FMM algorithm.



---

*Source: 938242-2014-03-16.xml*