# Extended Target Shape Estimation by Fitting B-Spline Curve

Authors: Jin-long Yang; Peng Li; Hong-wei Ge
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Taking into account the difficulty of shape estimation for the extended targets, a novel algorithm is proposed by fitting the B-spline curve. For the single extended target tracking, the multiple frame statistic technique is introduced to construct the pseudomeasurement sets and the control points are selected to form the B-spline curve. Then the shapes of the extended targets are extracted under the Bayes framework. Furthermore, the proposed shape estimation algorithm is modified suitably and combined with the probability hypothesis density (PHD) filter for multiple extended target tracking. Simulations show that the proposed algorithm has a good performance for shape estimate of any extended targets.


---

## Body

## 1. Introduction

In the traditional low resolution sensor system, each target is tracked as a single point source; that is, its extension is assumed to be neglectable in comparison with sensor resolution. With the increase of the resolution of modern radars and other detection equipment, the echo signal of a target may be distributed in a different range resolution cell; thus, the measurement is no longer equivalent to a point; that is, a single target may generate multiple measurements. Such target is referred to as an extended target in [1–4]. Recently, extended target tracking (ETT) is a hot topic in the field including the short-range applications or maritime surveillance, which has drawn a considerable attention [5–9].In the conventional extended target tracking, the measurements are modeled as a spatial distribution model in [1], and two examples—a point target with more measurement sources and an object with infinitely thin stick—are used to prove the effectiveness of the approach. Poisson process with a spatially dependent rate parameter is introduced in [2], assuming that each target produces measurements with Poisson distributed random number. It is considered that in this measurement model, the target is sufficiently far away from the sensor, and the measurements resemble a point cluster rather than a geometric structure [7]. Random matrix (RM) is proposed in [8], which has been used to track elliptical target extension [9]. Another method is random hypersurface model (RHM) [3] which is employed for modeling the target extent. However, these methods can only effectively achieve the shape estimate for the target with similar ellipse shape. They cannot effectively estimate the irregular shape of the extended target. In [4], star-convex target extension estimation method is proposed based on RHM under the condition that a measurement source is assumed to be an element of a randomly scaled version of the shape boundary. Moreover, the one-dimension probability density needs to be specified in advance in star-convex shape estimation and it is assumed to be independent of the shape.To solve the aforementioned problem, a novel shape estimation algorithm based on the B-spline curve fitting is proposed in this paper, and then the proposed shape estimation method is integrated into the framework of extended target probability hypothesis density (ET-PHD) filter for multiple extended target tracking [5, 6]. Simulations show that the proposed algorithm has a good performance for shape estimate of any extended targets.


## 2. Backgrounds

### 2.1. Kalman Filter

Assume the state equation and the measurement equation of a single target in two-dimensional plane are given by(1)
x
k
+
1
=
F
x
k
+
G
w
k
(2)
y
k
=
H
x
k
+
v
k
,
where x
k and y
k are the state vector and the observation vector at time k, respectively. F and H are the transition matrix and the measurement matrix, respectively. w
k and v
k are the process noise and the observation noise and are uncorrelated Gaussian white noise vectors with covariance matrixes Q
k and R
k, respectively.Suppose thatx
k
∣
k and the covariance P
k
∣
k are optimal estimations at time k in the fusion center; then, the recursive steps of the Kalman Filter (KF) at time k
+
1 are as follows [10].(1)
Prediction of state and covariance:(3)
x
k
+
1
∣
k
=
F
x
k
∣
k
,
(4)
P
k
+
1
∣
k
=
F
P
k
∣
k
F
T
+
G
Q
k
G
T
.
(2)
Calculating gain:(5)
K
k
+
1
=
P
k
+
1
∣
k
H
T
[
S
k
+
1
∣
k
]
-
1
,
S
k
+
1
∣
k
=
H
P
k
+
1
∣
k
H
T
+
R
k
+
1
.
(3)
Update of state and covariance:(6)
x
k
+
1
∣
k
+
1
=
x
k
+
1
∣
k
+
K
k
+
1
(
y
k
+
1
-
H
x
k
+
1
∣
k
)
,
(7)
P
k
+
1
∣
k
+
1
=
[
I
-
K
k
+
1
H
]
P
k
+
1
∣
k
.


### 2.2. The B-Spline Curve

Assume that readers are familiar with the concepts of B-spline curves. A smooth subsection curve can be obtained by fitting the control point set. The B-spline curve of orderl can be described as [11]
(8)
S
(
t
)
=
∑
i
=
1
n
μ
i
·
N
i
,
l
(
t
)
,
where μ
i in U
k
=
[
μ
1
,
μ
2
,
…
,
μ
n
]
T is a control point and N
i
,
l
(
t
) is the B-spline basis function, which is defined over a knot vector T
=
{
t
0
,
…
,
t
n
+
l
}. The basis function can be recursively defined as [11, 12]
(9)
N
i
,
l
(
t
)
=
t
-
t
i
t
i
+
l
-
1
-
t
i
N
i
,
l
-
1
(
t
)
+
t
i
+
l
-
t
t
i
+
l
-
t
i
+
1
N
i
+
1
,
l
-
1
(
t
)
,
N
i
,
1
(
t
)
=
{
1
,
if
t
i
≤
t
≤
t
i
+
1
0
,
otherwise
.
When l
=
3, we can obtain that
(10)
N
i
,
3
(
t
)
=
{
(
t
-
t
i
)
2
(
t
i
+
1
-
t
i
)
(
t
i
+
2
-
t
i
)
,
t
∈
[
t
i
,
t
i
+
1
)
(
t
-
t
i
)
(
t
i
+
2
-
t
)
(
t
i
+
2
-
t
i
)
(
t
i
+
2
-
t
i
+
1
)
+
(
t
-
t
i
+
1
)
(
t
i
+
3
-
t
)
(
t
i
+
2
-
t
i
+
1
)
(
t
i
+
3
-
t
i
+
1
)
,
t
∈
[
t
i
+
1
,
t
i
+
2
)
(
t
i
+
3
-
t
)
2
(
t
i
+
3
-
t
i
+
1
)
(
t
i
+
3
-
t
i
+
2
)
,
t
∈
[
t
i
+
2
,
t
i
+
3
)
0
,
otherwise
.


## 2.1. Kalman Filter

Assume the state equation and the measurement equation of a single target in two-dimensional plane are given by(1)
x
k
+
1
=
F
x
k
+
G
w
k
(2)
y
k
=
H
x
k
+
v
k
,
where x
k and y
k are the state vector and the observation vector at time k, respectively. F and H are the transition matrix and the measurement matrix, respectively. w
k and v
k are the process noise and the observation noise and are uncorrelated Gaussian white noise vectors with covariance matrixes Q
k and R
k, respectively.Suppose thatx
k
∣
k and the covariance P
k
∣
k are optimal estimations at time k in the fusion center; then, the recursive steps of the Kalman Filter (KF) at time k
+
1 are as follows [10].(1)
Prediction of state and covariance:(3)
x
k
+
1
∣
k
=
F
x
k
∣
k
,
(4)
P
k
+
1
∣
k
=
F
P
k
∣
k
F
T
+
G
Q
k
G
T
.
(2)
Calculating gain:(5)
K
k
+
1
=
P
k
+
1
∣
k
H
T
[
S
k
+
1
∣
k
]
-
1
,
S
k
+
1
∣
k
=
H
P
k
+
1
∣
k
H
T
+
R
k
+
1
.
(3)
Update of state and covariance:(6)
x
k
+
1
∣
k
+
1
=
x
k
+
1
∣
k
+
K
k
+
1
(
y
k
+
1
-
H
x
k
+
1
∣
k
)
,
(7)
P
k
+
1
∣
k
+
1
=
[
I
-
K
k
+
1
H
]
P
k
+
1
∣
k
.


## 2.2. The B-Spline Curve

Assume that readers are familiar with the concepts of B-spline curves. A smooth subsection curve can be obtained by fitting the control point set. The B-spline curve of orderl can be described as [11]
(8)
S
(
t
)
=
∑
i
=
1
n
μ
i
·
N
i
,
l
(
t
)
,
where μ
i in U
k
=
[
μ
1
,
μ
2
,
…
,
μ
n
]
T is a control point and N
i
,
l
(
t
) is the B-spline basis function, which is defined over a knot vector T
=
{
t
0
,
…
,
t
n
+
l
}. The basis function can be recursively defined as [11, 12]
(9)
N
i
,
l
(
t
)
=
t
-
t
i
t
i
+
l
-
1
-
t
i
N
i
,
l
-
1
(
t
)
+
t
i
+
l
-
t
t
i
+
l
-
t
i
+
1
N
i
+
1
,
l
-
1
(
t
)
,
N
i
,
1
(
t
)
=
{
1
,
if
t
i
≤
t
≤
t
i
+
1
0
,
otherwise
.
When l
=
3, we can obtain that
(10)
N
i
,
3
(
t
)
=
{
(
t
-
t
i
)
2
(
t
i
+
1
-
t
i
)
(
t
i
+
2
-
t
i
)
,
t
∈
[
t
i
,
t
i
+
1
)
(
t
-
t
i
)
(
t
i
+
2
-
t
)
(
t
i
+
2
-
t
i
)
(
t
i
+
2
-
t
i
+
1
)
+
(
t
-
t
i
+
1
)
(
t
i
+
3
-
t
)
(
t
i
+
2
-
t
i
+
1
)
(
t
i
+
3
-
t
i
+
1
)
,
t
∈
[
t
i
+
1
,
t
i
+
2
)
(
t
i
+
3
-
t
)
2
(
t
i
+
3
-
t
i
+
1
)
(
t
i
+
3
-
t
i
+
2
)
,
t
∈
[
t
i
+
2
,
t
i
+
3
)
0
,
otherwise
.


## 3. Single Extended Target Shape Estimation

In this section, the Bayesian filter framework is introduced for single extended target state and shape estimates. Assume the state equation and the measurement equation are the same as (1) and (2). We define x
k as the center of the extended target and X
k as the control matrix including the shape information. They are not related to each other. Thus we can estimate them with two parallel KFs, which can modify them by recursion update. And the detailed steps of the proposed algorithm are described as follows.Step 1.
At timek
=
0, initialize parameters x
0
∣
0, P
0
∣
0, X
0
∣
0, and Δ
0
∣
0, where P
0
∣
0 and Δ
0
∣
0 denote the state covariance and the shape covariance, respectively.Step 2.
Whenk
≥
1, implement the KF for state estimate:(2.1)
prediction of the state and covariance according to (3) and (4), respectively;
(2.2)
update of the state by the latest measurement setY
k
+
1:
(11)
x
k
+
1
∣
k
+
k
=
x
k
+
1
∣
k
+
K
k
+
1
(
1
|
Y
k
+
1
|
∑
i
=
1
|
Y
k
+
1
|
y
k
+
1
,
i
-
H
x
k
+
1
∣
k
)
,where Y
k
+
1
=
{
y
k
+
1
,
i
}
i
=
1
,
…
,
n and |
Y
k
+
1
| denotes the number of the measurements at time k
+
1. y
k
+
1
,
i denotes the ith measurement in Y
k
+
1. The covariance P
k
+
1
∣
k
+
1 can be updated by (7).Step 3.
Construct a new pseudomeasurement setZ
k
+
1.(3.1)
LetY
~
k
+
1
=
{
y
k
+
1
,
i
-
H
x
k
+
1
∣
k
∣
y
k
+
1
,
i
∈
Y
k
+
1
}
i
=
1
|
Y
k
+
1
|, and then add Y
~
k
+
1 to the pseudomeasurement set Z
k, that is, Z
^
k
+
1
=
Z
k
∪
Y
~
k
+
1.
(3.2)
Setm as the largest number of pseudomeasurements. If |
Z
^
k
+
1
|
≤
m, then let the new pseudomeasurement set Z
k
+
1
=
Z
^
k
+
1. Otherwise, delete the |
Z
^
k
+
1
|
-
m foregoing measurements from Z
^
k
+
1, and then let Z
k
+
1
=
Z
^
k
+
1.Step 4.
Update the shape of the extended target by the pseudomeasurement setZ
k
+
1.
(4.1) Divide the interval[
0,2
π
] into n equal angles, and generate a fixed angle set ℵ
=
{
θ
i
}
i
=
1
n. Then partition the pseudomeasurement set by the angles, and define D
k
+
1
,
i as the measurement set in the area of the ith partition. That is,
(12)
D
k
+
1
,
i
=
{
z
k
+
1
,
j
∣
d
L
i
(
z
k
+
1
,
j
)
<
d
,
C
(
z
k
+
1
,
j
)
=
ture
}
,
d
L
i
(
z
k
+
1
,
j
)
=
∥
B
1
z
k
+
1
,
j
(
1
)
+
B
2
z
k
+
1
,
j
(
2
)
∥
B
1
2
+
B
2
2
,
C
(
z
k
+
1
,
j
)
=
{
true
,
L
-
i
:
A
1
z
k
+
1
,
j
(
1
)
+
A
2
z
k
+
1
,
j
(
2
)
>
0
,
mlmmmmmmm
θ
i
∈
(
0
,
π
]
or
L
-
i
:
A
1
z
k
+
1
,
j
(
1
)
+
A
2
z
k
+
1
,
j
(
2
)
<
0
,
mmmmmmmmll
θ
i
∈
[
π
,
2
π
)
false
,
otherwise
,
where d is the width of partition. z
k
+
1
,
j
∈
Z
k
+
1, z
k
+
1
,
j
=
[
z
k
+
1
,
j
(
1
)
,
z
k
+
1
,
j
(
2
)
], z
k
+
1
,
j
(
1
), and z
k
+
1
,
j
(
2
) denote the position of x- and y-coordinate at time k
+
1. L
i denotes a line through the origin of coordinate and along the ith partition angle direction. d
L
i
(
z
k
+
1
,
j
) denotes the perpendicular distance between z
k
+
1
,
j and line L
i. The parameters B
1 and B
2 of line L
i satisfy B
1
/
B
2
=
-
tan
(
θ
i
). L
-
i denotes a line through the origin of coordinate, and it is perpendicular to line L
i. The parameters A
1 and A
2 of line L
-
i satisfy A
2
/
A
1
=
tan
(
θ
i
). Figure 1 shows the ith angle partition and the measurements in the area enclosed by points O1, O2, O3, and O4 which constitute the measurement set D
k
+
1
,
i.
(4.2) CalculateE
(
∥
D
k
,
i
∥
L
i
) which is the expectation of the perpendicular distance between the elements of D
k
+
1
,
i and line L
-
i. Then the radial extension distance ρ
k
+
1
,
i can be obtained by
(13)
ρ
k
+
1
,
i
=
2
1
|
D
k
+
1
,
i
|
∑
j
=
1
|
D
k
+
1
,
i
|
∥
A
1
z
k
+
1
,
j
(
1
)
+
A
2
z
k
+
1
,
j
(
2
)
∥
A
1
2
+
A
2
2
.
In this paper, we set the ρ
k
+
1
,
i as the control point of the extended target shape parameter and define X
k
∣
k
=
[
{
ρ
k
,
i
}
i
=
1
n
]
T as the control matrix of the shape, which includes the shape information.
(4.3) Shape estimation by implementing the one-dimension (1-D) KF: assume the shape control matrixX
k
∣
k and the shape covariance Δ
k
∣
k have been obtained. Then the recursion steps of shape estimates are as follows.
Prediction:(14)
X
k
+
1
∣
k
=
X
k
|
k
,
Δ
k
+
1
∣
k
=
Δ
k
∣
k
+
H
-
∥
R
k
+
1
∥
,
where H
-
=
[
1,1
,
…
,
1
]
︸
n
T and R
k
+
1 denotes the measurement noise covariance.
Update:(15)
X
k
+
1
∣
k
+
1
=
X
k
+
1
∣
k
+
(
X
k
∣
k
-
K
^
k
+
1
·
2
E
(
∥
D
k
+
1
,
i
∥
)
)
,
Δ
k
+
1
∣
k
+
1
=
(
I
-
K
^
k
+
1
H
-
T
)
Δ
k
+
1
∣
k
,
K
^
k
+
1
=
diag
⁡
(
{
δ
k
+
1
∣
k
,
i
δ
k
+
1
∣
k
,
i
+
v
k
+
1
}
i
=
1
n
)
,
where the shape covariance Δ
k
+
1
∣
k
+
1
=
[
δ
k
+
1
,
i
]
n
×
1. v
k
+
1 denotes the measurement noise.Figure 1
Angle interval partition.Step 5.
Shape estimation according toX
k
+
1
|
k
+
1 and the B-spline curve fitting technique.(5.1)
Map the control points to the Cartesian coordinates by(16)
U
k
+
1
=
{
μ
i
}
i
=
1
n
=
ℵ
^
X
k
+
1
|
k
+
1
,
ℵ
^
=
[
cos
⁡
(
θ
i
)
,
…
,
cos
⁡
(
θ
n
)
sin
(
θ
i
)
,
…
,
sin
(
θ
n
)
]
2
×
n
.
(5.2)
Produce a closed control point set by adding the elementμ
1
,
μ
2
,
μ
3 to the end of the U
k
+
1, and describe it as U
_
k
+
1
=
{
μ
_
i
}
i
=
1
n
+
3. Then the closed cubic B-spline curve can be obtained by
(17)
S
(
t
)
=
∑
i
=
1
n
+
3
μ
_
k
+
1
·
N
i
,
3
(
t
)
.


## 4. Multiple Extended Target Shape Estimation

### 4.1. Multiple Extended Target PHD Filter

The standard PHD filter for single measurement target tracking has been described in [13–15]; however, it is not suitable for METT. Recently, Mahler has derived a correct equation for extended target PHD (ET-PHD) filter based on the Poisson multitarget measurement model [5]. The Gaussian mixture implement of the multiple extended target PHD filter is presented in [6] and referred to as ET-GM-PHD filter. The detailed filter process of ET-GM-PHD can refer to [6, 16, 17].


### 4.2. Multiple Extended Target Tracking Algorithm

In this section, we combine the proposed shape estimation algorithm into the framework of ET-GM-PHD filter, which can effectively achieve the multiple extended target tracking with different shape estimation. We refer to this algorithm as Shape-ET-GM-PHD, and its steps are as follows.(1) Prediction. Assume that the state vector x and its shape vector X are independent, and the survival and detection probabilities are independent of them; that is, p
S
,
k
(
x
,
X
)
=
p
S and p
D
,
k
(
x
,
X
)
=
p
D. Assume that v
k
-
1
(
x
,
X
) denotes the joint posterior PHD function at time k
-
1 and is approximated by the Gaussian mixture distributions. Then the predicted PHD function v
k
∣
k
-
1
(
x
,
X
) can be described as
(18)
v
k
∣
k
-
1
(
x
,
X
)
=
p
S
∬
p
k
∣
k
-
1
(
x
,
X
∣
x
′
,
X
′
)
v
k
-
1
∣
k
-
1
(
x
,
X
)
d
x
′
d
X
′
=
p
S
∬
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
p
k
∣
k
-
1
(
X
∣
X
′
)
×
v
k
-
1
∣
k
-
1
(
x
,
X
)
d
x
′
d
X
′
≈
p
S
∑
j
=
1
J
k
∣
k
-
1
w
k
-
1
∣
k
-
1
∫
N
(
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
+
σ
-
k
-
1
∣
k
-
1
I
~
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
)
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
d
x
′
·
∫
Ψ
(
X
k
-
1
∣
k
-
1
∣
X
~
k
-
1
∣
k
-
1
,
Δ
k
-
1
∣
k
-
1
)
p
k
∣
k
-
1

(
X
∣
X
′
)

d
X
′

,
where
(19)
∫
N
(
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
+
σ
-
k
-
1
∣
k
-
1
I
~
)
×
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
d
x
′
=
N
(
x
k
∣
k
-
1
∣
m
k
∣
k
-
1
(
j
)
,
P
k
∣
k
-
1
(
j
)
+
σ
-
k
∣
k
-
1
I
~
)
,
∫
Ψ
(
X
k
-
1
∣
k
-
1
∣
X
~
k
-
1
∣
k
-
1
,
Δ
k
-
1
∣
k
-
1
)
p
k
∣
k
-
1
(
X
∣
X
′
)
d
X
′
=
∏
i
=
1
Num
N
(
X
k
∣
k
-
1
(
i
)
∣
X
~
k
∣
k
-
1
(
j
,
i
)
,
Δ
k
∣
k
-
1
(
j
,
i
)
)
,
I
~
=
diag
⁡
[
1,1
,
0,0
], and σ
-
k
-
1
∣
k
-
1 denotes the mean of the shape standard variance, One has m
k
∣
k
-
1
(
j
)
=
F
k
∣
k
-
1
m
k
-
1
∣
k
-
1
(
j
), P
k
∣
k
-
1
(
j
)
=
F
k
∣
k
-
1
P
k
-
1
∣
k
-
1
(
j
)
F
k
∣
k
-
1
T
+
Q
k
∣
k
-
1. Num denotes the divided number of the shape area according to the degree of angle. X
~
k
∣
k
-
1
(
j
,
i
)
=
X
~
k
-
1
∣
k
-
1
(
j
,
i
) and Δ
k
∣
k
-
1
(
j
,
i
)
=
Δ
k
-
1
∣
k
-
1
(
j
,
i
)
+
Q
k denote the control matrix of shape and its predicted covariance, and Q
k denotes the shape process noise covariance. X
~
k
∣
k
-
1
(
j
,
i
) and Δ
k
∣
k
-
1
(
j
,
i
) can also be described as X
~
k
∣
k
-
1
(
j
,
i
)
=
[
{
ρ
~
k
∣
k
-
1
(
j
,
i
)
}
i
=
1
Num
]
T and Δ
k
∣
k
-
1
(
j
,
i
)
=
[
{
σ
k
∣
k
-
1
(
j
,
i
)
}
i
=
1
Num
]
T, where the definition of ρ
~
k
∣
k
-
1
(
j
,
i
) is the same as that of (13). σ
k
∣
k
-
1
(
j
,
i
) denotes the predicted variance of the shape.(2) Update. The updated formula of the intensity function can be described as
(20)
v
k
∣
k
(
x
,
X
)
=
L
Z
k
(
x
,
X
)
v
k
∣
k
-
1
(
x
,
X
)
,
where v
k
∣
k
(
x
,
X
) denotes an intensity function of the state x with shape information X and L
Z
k
(
x
,
X
) denotes the pseudolikelihood function. Z
k denotes the measurement set at time k. When Z
k
=
∅,
(21)
L
Z
k
(
x
,
X
)
≜
1
-
p
D
(
x
,
X
)
+
e
-
r
(
x
,
X
)
p
D
(
x
,
X
)
;
otherwise,
(22)
L
Z
k
(
x
,
X
)
≜
1
-
(
1
-
e
-
r
(
x
,
X
)
)
p
D
(
x
,
X
)
+
e
-
r
(
x
,
X
)
p
D
(
x
,
X
)
×
∑
P
∠
Z
k
ω
P
∑
W
∈
P
r
(
x
,
X
)
|
W
|
d
W
∏
z
∈
W
ϕ
z
(
x
,
X
)
λ
k
c
k
,
where P
∠
Z
k denotes a partition subset P of the measurement set Z
k, W denotes a subset of a partition P, and ∪
W
∈
P
W
=
Z
k. ϕ
z
(
x
,
X
) denotes the measurement likelihood function of one measurement originating from an extended target x, r
(
x
,
X
) is the measurement expectation, and p
D
(
x
,
X
)
=
p
D denotes the detection probability of the sensor. Clutter has a Poisson distribution, and its density can be described as β
F
A
,
k
=
λ
k
c
k
(
z
k
), λ
k denotes the mean number of clutter measurements, and c
k
(
z
k
) is the space distribution of the clutter. Consider
(23)
ω
P
=
∏
W
∈
P
d
W
∑
P
′
∠
Z
k
∏
W
∈
P
′
d
W
′
,
d
W
=
δ
|
W
|
,
1
+
v
k
∣
k
-
1
[
e
-
r
r
|
W
|
p
d
∏
z
∈
W
ϕ
z
β
F
A
,
k
]

,
where δ
i
,
j denotes the Kronecker delta function and r
=
r
(
x
,
X
). β
F
A
,
k denotes the clutter distribution of the measurement space.In the following subsections, assume that the current estimated PHDv
k
∣
k
(
x
,
X
) can be approximated as a Gaussian mixture distribution. The corrected PHD can be described as
(24)
v
k
∣
k
(
x
,
X
)
=
v
k
∣
k
N
D
(
x
,
X
)
+
∑
p
∈
Z
∑
W
∈
p
v
k
∣
k
D
(
x
,
X
,
W
)
,
where v
k
∣
k
N
D
(
x
,
X
) denotes the PHD of the targets without detecting cases, and it can be described as
(25)
v
k
∣
k
N
D
(
x
,
X
)
=
∑
j
=
1
J
k
∣
k
-
1
w
N
D
,
k
∣
k
(
j
)
N
(
x
k
∣
k
∣
m
k
∣
k
(
j
)
,
P
k
∣
k
(
j
)
+
σ
-
k
∣
k
(
j
)
I
~
)
×
∏
i
=
1
Num
N
(
X
k
∣
k
(
i
)
∣
X
~
k
∣
k
(
j
,
i
)
,
Δ
k
∣
k
(
j
,
i
)
)
,where   m
k
∣
k
(
j
)
=
m
k
∣
k
-
1
(
j
), P
k
∣
k
(
j
)
=
P
k
∣
k
-
1
(
j
), Δ
k
∣
k
(
j
,
i
)
=
Δ
k
∣
k
-
1
(
j
,
i
), and X
~
k
∣
k
(
j
,
i
)
=
X
~
k
∣
k
-
1
(
j
,
i
).v
k
|
k
D
(
x
,
X
,
W
) denotes the PHD of the detected target cases and can be described as
(26)
v
k
∣
k
D
(
x
,
X
,
W
)
=
∑
j
=
1
J
k
∣
k
-
1
w
k
∣
k
(
j
)
N
(
x
k
∣
m
k
∣
k
(
j
)
,
P
k
∣
k
(
j
)
+
σ
-
k
∣
k
(
j
)
I
~
)
×
∏
i
=
1
Num
N
(
X
k
(
i
)
∣
X
~
k
∣
k
(
j
,
i
)
,
Δ
k
∣
k
(
j
,
i
)
)
,
where
(27)
m
k
∣
k
(
j
)
=
m
k
∣
k
-
1
(
j
)
+
K
k
(
z
-
k
W
-
H
m
k
∣
k
-
1
(
j
)
)
,
z
-
k
W
=
1
|
W
|
∑
z
k
(
l
)
∈
W
z
k
(
l
)
,
P
k
∣
k
(
j
)
=
(
I
-
K
k
H
)
P
k
∣
k
-
1
(
j
)
,
S
k
=
P
k
∣
k
-
1
(
j
)
+
σ
-
k
∣
k
(
j
)
[
1
0
0
0
]
⊗
I
2
,
K
k
=
P
k
∣
k
-
1
(
j
)
H
T
(
S
k
)
-
1
,
σ
-
k
∣
k
(
j
)
=
1
Num
∑
i
=
1
Num
σ
k
∣
k
(
j
,
i
)
,
w
k
(
j
)
=
w
p
d
W
e
-
γ
(
j
)
(
γ
(
j
)
β
F
A
,
k
)
|
W
|
p
D
L
k
(
j
,
W
)
w
k
∣
k
-
1
(
j
)
,
w
p
=
∏
W
∈
p
d
W
∑
p
′
∠
Z
′
∏
W
′
∈
p
′
d
W
′
,
d
W
=
δ
|
W
|
,
1
+
∑
l
=
1
J
k
|
k
-
1
e
-
γ
(
l
)
(
γ
(
j
)
β
F
A
,
k
)
|
W
|
p
D
L
k
(
l
,
W
)
w
k
∣
k
-
1
(
l
)
,
L
k
(
j
,
W
)
=
∏
i
=
1
Num
L
^
k
(
j
,
W
,
i
)
,
L
^
k
(
j
,
W
,
i
)
=
{
(
L
~
k
(
j
,
W
,
i
)
)
|
W
(
i
)
|
/
Num
,
|
W
(
i
)
|
≠
0
α
·
(
L
~
k
(
j
,
W
,
i
)
)
1
/
Num
,
|
W
(
i
)
|
=
0
,
L
~
k
(
j
,
W
,
i
)
=
1
2
π
σ
k
∣
k
-
1
(
j
,
i
)
·
exp
⁡
(
-
1
2
σ
k
∣
k
-
1
2
(
j
,
i
)
(
ρ
~
k
+
1
∣
k
(
j
,
i
)
-
ρ
~
(
W
,
i
)
)
2
)
,
where α denotes the penalty coefficient. The shape parameters can be obtained by
(28)
Δ
k
∣
k
(
j
,
i
)
=
(
I
Num
-
K
~
k
H
-
T
)
Δ
k
∣
k
-
1
(
j
,
i
)
,
X
~
k
∣
k
(
j
,
i
)
=
X
~
k
∣
k
-
1
(
j
,
i
)
+
K
~
k
(
X
~
k
(
W
)
-
X
~
k
∣
k
-
1
(
j
,
i
)
)
,
K
~
k
=
I
Num
Δ
k
∣
k
-
1
(
j
)
(
I
Num
S
~
k
)
-
1
,
S
~
k
=
Υ
Num
(
H
P
k
∣
k
-
1
H
T
+
R
k
|
W
|
)
+
Δ
k
∣
k
-
1
(
j
)
,
where I
Num denotes the identity matrix with Num order. R
k denotes the measurement noise of a measurement source. Υ
Num
(
P
) denotes the decomposed function which can decompose the matrix P as a set with Num variances. X
~
k
(
W
)
=
[
{
ρ
~
(
W
,
i
)
}
i
=
1
Num
]
T denotes the pseudocontrol point matrix which is used to update the control point matrix X
~
k
∣
k
(
j
,
i
), where ρ
~
(
W
,
i
) can be obtained by
(29)
ρ
~
(
W
,
i
)
=
{
1
|
W
(
i
)
|
∑
z
∈
W
(
i
)
(
z
-
z
-
k
W
)
T
(
z
-
z
-
k
W
)
,
|
W
(
i
)
|
≠
0
0
,
|
W
(
i
)
|
=
0
,
where |
W
(
i
)
| denotes the measurement number of ith angle direction in cell W.Notice that the D-distance partition method [17] is implemented in the part of measurement partition. Set the maximum distance as the mean size of the shape, and it can be obtained by
(30)
ρ
~
k
,
max
⁡
=
∑
j
=
1
J
k
-
1
∣
k
-
1
w
k
-
1
∣
k
-
1
(
j
)
∑
i
=
1
Num
2
ρ
~
k
-
1
∣
k
-
1
(
j
,
i
)
Num
·
∑
j
=
1
J
k
-
1
∣
k
-
1
w
k
-
1
∣
k
-
1
(
j
)
.Finally, shapes of multiple extended targets are extracted according toX
k
+
1 and the B-spline curve fitting technique described in Section 2.2.


## 4.1. Multiple Extended Target PHD Filter

The standard PHD filter for single measurement target tracking has been described in [13–15]; however, it is not suitable for METT. Recently, Mahler has derived a correct equation for extended target PHD (ET-PHD) filter based on the Poisson multitarget measurement model [5]. The Gaussian mixture implement of the multiple extended target PHD filter is presented in [6] and referred to as ET-GM-PHD filter. The detailed filter process of ET-GM-PHD can refer to [6, 16, 17].


## 4.2. Multiple Extended Target Tracking Algorithm

In this section, we combine the proposed shape estimation algorithm into the framework of ET-GM-PHD filter, which can effectively achieve the multiple extended target tracking with different shape estimation. We refer to this algorithm as Shape-ET-GM-PHD, and its steps are as follows.(1) Prediction. Assume that the state vector x and its shape vector X are independent, and the survival and detection probabilities are independent of them; that is, p
S
,
k
(
x
,
X
)
=
p
S and p
D
,
k
(
x
,
X
)
=
p
D. Assume that v
k
-
1
(
x
,
X
) denotes the joint posterior PHD function at time k
-
1 and is approximated by the Gaussian mixture distributions. Then the predicted PHD function v
k
∣
k
-
1
(
x
,
X
) can be described as
(18)
v
k
∣
k
-
1
(
x
,
X
)
=
p
S
∬
p
k
∣
k
-
1
(
x
,
X
∣
x
′
,
X
′
)
v
k
-
1
∣
k
-
1
(
x
,
X
)
d
x
′
d
X
′
=
p
S
∬
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
p
k
∣
k
-
1
(
X
∣
X
′
)
×
v
k
-
1
∣
k
-
1
(
x
,
X
)
d
x
′
d
X
′
≈
p
S
∑
j
=
1
J
k
∣
k
-
1
w
k
-
1
∣
k
-
1
∫
N
(
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
+
σ
-
k
-
1
∣
k
-
1
I
~
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
)
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
d
x
′
·
∫
Ψ
(
X
k
-
1
∣
k
-
1
∣
X
~
k
-
1
∣
k
-
1
,
Δ
k
-
1
∣
k
-
1
)
p
k
∣
k
-
1

(
X
∣
X
′
)

d
X
′

,
where
(19)
∫
N
(
x
k
-
1
∣
k
-
1
∣
m
k
-
1
∣
k
-
1
(
j
)
,
P
k
-
1
∣
k
-
1
(
j
)
+
σ
-
k
-
1
∣
k
-
1
I
~
)
×
p
k
∣
k
-
1
(
x
∣
x
′
,
X
′
)
d
x
′
=
N
(
x
k
∣
k
-
1
∣
m
k
∣
k
-
1
(
j
)
,
P
k
∣
k
-
1
(
j
)
+
σ
-
k
∣
k
-
1
I
~
)
,
∫
Ψ
(
X
k
-
1
∣
k
-
1
∣
X
~
k
-
1
∣
k
-
1
,
Δ
k
-
1
∣
k
-
1
)
p
k
∣
k
-
1
(
X
∣
X
′
)
d
X
′
=
∏
i
=
1
Num
N
(
X
k
∣
k
-
1
(
i
)
∣
X
~
k
∣
k
-
1
(
j
,
i
)
,
Δ
k
∣
k
-
1
(
j
,
i
)
)
,
I
~
=
diag
⁡
[
1,1
,
0,0
], and σ
-
k
-
1
∣
k
-
1 denotes the mean of the shape standard variance, One has m
k
∣
k
-
1
(
j
)
=
F
k
∣
k
-
1
m
k
-
1
∣
k
-
1
(
j
), P
k
∣
k
-
1
(
j
)
=
F
k
∣
k
-
1
P
k
-
1
∣
k
-
1
(
j
)
F
k
∣
k
-
1
T
+
Q
k
∣
k
-
1. Num denotes the divided number of the shape area according to the degree of angle. X
~
k
∣
k
-
1
(
j
,
i
)
=
X
~
k
-
1
∣
k
-
1
(
j
,
i
) and Δ
k
∣
k
-
1
(
j
,
i
)
=
Δ
k
-
1
∣
k
-
1
(
j
,
i
)
+
Q
k denote the control matrix of shape and its predicted covariance, and Q
k denotes the shape process noise covariance. X
~
k
∣
k
-
1
(
j
,
i
) and Δ
k
∣
k
-
1
(
j
,
i
) can also be described as X
~
k
∣
k
-
1
(
j
,
i
)
=
[
{
ρ
~
k
∣
k
-
1
(
j
,
i
)
}
i
=
1
Num
]
T and Δ
k
∣
k
-
1
(
j
,
i
)
=
[
{
σ
k
∣
k
-
1
(
j
,
i
)
}
i
=
1
Num
]
T, where the definition of ρ
~
k
∣
k
-
1
(
j
,
i
) is the same as that of (13). σ
k
∣
k
-
1
(
j
,
i
) denotes the predicted variance of the shape.(2) Update. The updated formula of the intensity function can be described as
(20)
v
k
∣
k
(
x
,
X
)
=
L
Z
k
(
x
,
X
)
v
k
∣
k
-
1
(
x
,
X
)
,
where v
k
∣
k
(
x
,
X
) denotes an intensity function of the state x with shape information X and L
Z
k
(
x
,
X
) denotes the pseudolikelihood function. Z
k denotes the measurement set at time k. When Z
k
=
∅,
(21)
L
Z
k
(
x
,
X
)
≜
1
-
p
D
(
x
,
X
)
+
e
-
r
(
x
,
X
)
p
D
(
x
,
X
)
;
otherwise,
(22)
L
Z
k
(
x
,
X
)
≜
1
-
(
1
-
e
-
r
(
x
,
X
)
)
p
D
(
x
,
X
)
+
e
-
r
(
x
,
X
)
p
D
(
x
,
X
)
×
∑
P
∠
Z
k
ω
P
∑
W
∈
P
r
(
x
,
X
)
|
W
|
d
W
∏
z
∈
W
ϕ
z
(
x
,
X
)
λ
k
c
k
,
where P
∠
Z
k denotes a partition subset P of the measurement set Z
k, W denotes a subset of a partition P, and ∪
W
∈
P
W
=
Z
k. ϕ
z
(
x
,
X
) denotes the measurement likelihood function of one measurement originating from an extended target x, r
(
x
,
X
) is the measurement expectation, and p
D
(
x
,
X
)
=
p
D denotes the detection probability of the sensor. Clutter has a Poisson distribution, and its density can be described as β
F
A
,
k
=
λ
k
c
k
(
z
k
), λ
k denotes the mean number of clutter measurements, and c
k
(
z
k
) is the space distribution of the clutter. Consider
(23)
ω
P
=
∏
W
∈
P
d
W
∑
P
′
∠
Z
k
∏
W
∈
P
′
d
W
′
,
d
W
=
δ
|
W
|
,
1
+
v
k
∣
k
-
1
[
e
-
r
r
|
W
|
p
d
∏
z
∈
W
ϕ
z
β
F
A
,
k
]

,
where δ
i
,
j denotes the Kronecker delta function and r
=
r
(
x
,
X
). β
F
A
,
k denotes the clutter distribution of the measurement space.In the following subsections, assume that the current estimated PHDv
k
∣
k
(
x
,
X
) can be approximated as a Gaussian mixture distribution. The corrected PHD can be described as
(24)
v
k
∣
k
(
x
,
X
)
=
v
k
∣
k
N
D
(
x
,
X
)
+
∑
p
∈
Z
∑
W
∈
p
v
k
∣
k
D
(
x
,
X
,
W
)
,
where v
k
∣
k
N
D
(
x
,
X
) denotes the PHD of the targets without detecting cases, and it can be described as
(25)
v
k
∣
k
N
D
(
x
,
X
)
=
∑
j
=
1
J
k
∣
k
-
1
w
N
D
,
k
∣
k
(
j
)
N
(
x
k
∣
k
∣
m
k
∣
k
(
j
)
,
P
k
∣
k
(
j
)
+
σ
-
k
∣
k
(
j
)
I
~
)
×
∏
i
=
1
Num
N
(
X
k
∣
k
(
i
)
∣
X
~
k
∣
k
(
j
,
i
)
,
Δ
k
∣
k
(
j
,
i
)
)
,where   m
k
∣
k
(
j
)
=
m
k
∣
k
-
1
(
j
), P
k
∣
k
(
j
)
=
P
k
∣
k
-
1
(
j
), Δ
k
∣
k
(
j
,
i
)
=
Δ
k
∣
k
-
1
(
j
,
i
), and X
~
k
∣
k
(
j
,
i
)
=
X
~
k
∣
k
-
1
(
j
,
i
).v
k
|
k
D
(
x
,
X
,
W
) denotes the PHD of the detected target cases and can be described as
(26)
v
k
∣
k
D
(
x
,
X
,
W
)
=
∑
j
=
1
J
k
∣
k
-
1
w
k
∣
k
(
j
)
N
(
x
k
∣
m
k
∣
k
(
j
)
,
P
k
∣
k
(
j
)
+
σ
-
k
∣
k
(
j
)
I
~
)
×
∏
i
=
1
Num
N
(
X
k
(
i
)
∣
X
~
k
∣
k
(
j
,
i
)
,
Δ
k
∣
k
(
j
,
i
)
)
,
where
(27)
m
k
∣
k
(
j
)
=
m
k
∣
k
-
1
(
j
)
+
K
k
(
z
-
k
W
-
H
m
k
∣
k
-
1
(
j
)
)
,
z
-
k
W
=
1
|
W
|
∑
z
k
(
l
)
∈
W
z
k
(
l
)
,
P
k
∣
k
(
j
)
=
(
I
-
K
k
H
)
P
k
∣
k
-
1
(
j
)
,
S
k
=
P
k
∣
k
-
1
(
j
)
+
σ
-
k
∣
k
(
j
)
[
1
0
0
0
]
⊗
I
2
,
K
k
=
P
k
∣
k
-
1
(
j
)
H
T
(
S
k
)
-
1
,
σ
-
k
∣
k
(
j
)
=
1
Num
∑
i
=
1
Num
σ
k
∣
k
(
j
,
i
)
,
w
k
(
j
)
=
w
p
d
W
e
-
γ
(
j
)
(
γ
(
j
)
β
F
A
,
k
)
|
W
|
p
D
L
k
(
j
,
W
)
w
k
∣
k
-
1
(
j
)
,
w
p
=
∏
W
∈
p
d
W
∑
p
′
∠
Z
′
∏
W
′
∈
p
′
d
W
′
,
d
W
=
δ
|
W
|
,
1
+
∑
l
=
1
J
k
|
k
-
1
e
-
γ
(
l
)
(
γ
(
j
)
β
F
A
,
k
)
|
W
|
p
D
L
k
(
l
,
W
)
w
k
∣
k
-
1
(
l
)
,
L
k
(
j
,
W
)
=
∏
i
=
1
Num
L
^
k
(
j
,
W
,
i
)
,
L
^
k
(
j
,
W
,
i
)
=
{
(
L
~
k
(
j
,
W
,
i
)
)
|
W
(
i
)
|
/
Num
,
|
W
(
i
)
|
≠
0
α
·
(
L
~
k
(
j
,
W
,
i
)
)
1
/
Num
,
|
W
(
i
)
|
=
0
,
L
~
k
(
j
,
W
,
i
)
=
1
2
π
σ
k
∣
k
-
1
(
j
,
i
)
·
exp
⁡
(
-
1
2
σ
k
∣
k
-
1
2
(
j
,
i
)
(
ρ
~
k
+
1
∣
k
(
j
,
i
)
-
ρ
~
(
W
,
i
)
)
2
)
,
where α denotes the penalty coefficient. The shape parameters can be obtained by
(28)
Δ
k
∣
k
(
j
,
i
)
=
(
I
Num
-
K
~
k
H
-
T
)
Δ
k
∣
k
-
1
(
j
,
i
)
,
X
~
k
∣
k
(
j
,
i
)
=
X
~
k
∣
k
-
1
(
j
,
i
)
+
K
~
k
(
X
~
k
(
W
)
-
X
~
k
∣
k
-
1
(
j
,
i
)
)
,
K
~
k
=
I
Num
Δ
k
∣
k
-
1
(
j
)
(
I
Num
S
~
k
)
-
1
,
S
~
k
=
Υ
Num
(
H
P
k
∣
k
-
1
H
T
+
R
k
|
W
|
)
+
Δ
k
∣
k
-
1
(
j
)
,
where I
Num denotes the identity matrix with Num order. R
k denotes the measurement noise of a measurement source. Υ
Num
(
P
) denotes the decomposed function which can decompose the matrix P as a set with Num variances. X
~
k
(
W
)
=
[
{
ρ
~
(
W
,
i
)
}
i
=
1
Num
]
T denotes the pseudocontrol point matrix which is used to update the control point matrix X
~
k
∣
k
(
j
,
i
), where ρ
~
(
W
,
i
) can be obtained by
(29)
ρ
~
(
W
,
i
)
=
{
1
|
W
(
i
)
|
∑
z
∈
W
(
i
)
(
z
-
z
-
k
W
)
T
(
z
-
z
-
k
W
)
,
|
W
(
i
)
|
≠
0
0
,
|
W
(
i
)
|
=
0
,
where |
W
(
i
)
| denotes the measurement number of ith angle direction in cell W.Notice that the D-distance partition method [17] is implemented in the part of measurement partition. Set the maximum distance as the mean size of the shape, and it can be obtained by
(30)
ρ
~
k
,
max
⁡
=
∑
j
=
1
J
k
-
1
∣
k
-
1
w
k
-
1
∣
k
-
1
(
j
)
∑
i
=
1
Num
2
ρ
~
k
-
1
∣
k
-
1
(
j
,
i
)
Num
·
∑
j
=
1
J
k
-
1
∣
k
-
1
w
k
-
1
∣
k
-
1
(
j
)
.Finally, shapes of multiple extended targets are extracted according toX
k
+
1 and the B-spline curve fitting technique described in Section 2.2.


## 5. Simulations

Assume that there is an extended target making a uniform motion in a two-dimensional simulation scenario, and the state equation and the measurement equation are the same as (1) and (2), respectively, where x
k
=
[
x
,
v
x
,
y
,
v
y
]
Tdenotes the target state, which contains position and velocity information. The state transition matrix F
=
[
1
1
0
0
0
1
0
0
0
0
1
1
0
0
0
1
] and measurement matrix H
=
[
1
0
0
0
0
1
0
0
]
T. The process noise covariance Q
k
=
diag
⁡
[
σ
w
1
2
,
σ
w
2
2
] with σ
w
1
=
σ
w
2
=
1, and the measurement noise covariance R
k
=
diag
⁡
[
σ
v
1
2
,
σ
v
2
2
] with σ
v
1
=
σ
v
2
=
1. The initial covariance P
0
∣
0
=
diag
⁡
[
5
,
1
,
5
,
1
] and the shape noise covariance is Δ
0
∣
0
=
[
1,1
,
…
,
1
]
︸
n
T. The other parameters ℵ
=
{
(
2
π
/
n
)
·
i
}
i
=
0
n, n
=
20, d
=
1.5, and m
=
45. The way to generate measurements is the same as that of the RM method. The real measurement is assumed as the scattering center, that is, a measurement source. The observation measurements are generated from the scattering center with measurement noises. Notice that the measurement noise is assumed to be small compared to the target extent.Example 1 (star-shaped extended target).
Figure2 shows the shape estimation by the proposed algorithm and the RM method [9], respectively. It is clear that the proposed algorithm has a higher accuracy than that of the RM method. The reason is that the RM method can only estimate the elliptical shape of the target.Figure 2
Shape estimation by the proposed algorithm and the RM method.Figure3 shows the measurements of the extended target from the 1st frame to the 20th frame. As can be seen, the measurement noise is big which makes the shape hard to be estimated from the single frame. However, the proposed algorithm can extract the accurate shape by multiple frame statistic technique and the B-spline curve fitting. It is shown that the proposed algorithm has a good capacity of resisting disturbance of the noise.Figure 3
The real measurements of the extended target (1~20 frames).Figure4 shows the average shape estimate. It is clear that the proposed algorithm can obtain the shape features effectively.Figure 4
Average shape estimate.Example 2 (Y-shaped extended target).
Figure5 shows the shape estimation by the proposed algorithm and the RM method, respectively. Figure 6 shows the measurements of the extended target from the 1st to the 20th frames. Figure 7 shows the average shape estimate. As can be seen, the proposed algorithm has a better performance for shape estimate of irregular extended targets.Figure 5
Shape estimation by the proposed algorithm and the RM method.Figure 6
The real measurements of the extended target (1~20 frames).Figure 7
Average shape estimate.Example 3 (multiple extended target tracking).
The scenario of multiple extended target tracking is the same as that of [17]. There are four targets and two of them cross at time k
=
56, and one target is spawned at time k
=
66. Assume that the measurement noise covariance and the process noise covariance are R
k
=
diag
⁡
[
0.7,0.7
] and Q
k
=
diag
⁡
[
1,1
], respectively. Shape noise covariance is Q
k
=
[
0.05,0.05
,
…
,
0.05
]
︸
Num, Num
=
12, and ω
=
0.01. The probabilities of target survival and detection are p
S
=
0.99 and p
D
=
0.99, respectively. The birth process is modeled on a Poisson RFS with Gaussian mixture intensity:
(31)
D
b
(
x
,
X
)
=
∑
j
=
1
2
ω
b
(
j
)
N
(
x
∣
m
b
(
j
)
,
P
b
(
j
)
+
σ
-
b
I
~
)
×
∏
i
=
1
Num
N
(
X
b
(
j
)
∣
X
~
b
(
j
,
i
)
,
Δ
b
(
j
,
i
)
)
,
where m
b
(
1
)
=
[
250
m
,
250
m
,
0
m
/
s
,
0
m
/
s
]
T, m
b
(
2
)
=
[
-
250
m
,
-
250
m
,
0
m
/
s
,
0
m
/
s
]
T, P
b
(
1
)
=
P
b
(
2
)
=
diag
⁡
[
100,100,15,15
], and ω
b
(
1
)
=
ω
b
(
2
)
=
0.1, X
~
b
(
1,1
)
=
⋯
=
X
~
b
(
1
,
Num
)
=
X
~
b
(
2,1
)
=
⋯
=
X
~
b
(
2
,
Num
)
=
5, Δ
b
(
1,1
)
=
⋯
=
Δ
b
(
1
,
Num
)
=
Δ
b
(
2,1
)
=
⋯
=
Δ
b
(
2
,
Num
)
=
1.The intensity of the spawned targets isv
β
(
x
)
=
0.05
N
(
x
∣
ξ
,
Q
β
)
×
∏
i
=
1
Num
N
(
X
∣
ξ
,
Δ
β
), where Q
β
=
diag
⁡
[
100,100,400,400
] and Δ
β
=
[
1,1
,
…
,
1
]
︸
Num.Figure8(a) shows the tracking results of the targets in the whole tracking area, and Figures 8(b), 8(c), and 8(d) show the shape estimates of different targets in special tracking area. As can be seen from Figure 8(b), the shapes of Targets 1 and 2 are assumed as ellipses; although they are not accurate, they can be updated approximately to the real shapes of the targets at several time points. The reason is that the shape parameters in the proposed shape estimate method are updated at each time when the latest measurements arrive. Moreover, the Gaussian mixture technique is employed, which can approximately fit the real shape distribution of the extended target. In Figure 8(c), we can see that the shape estimates are not accurate when the targets cross each other, but they can also be updated approximately to the real target shapes. Figure 8(d) shows the shape estimate of a spawned target; it is showed that the proposed algorithm also has a good performance for shape estimate of the spawned target.Shape estimate results.
(a)
Measurements and tracking results
(b)
Shape estimates of Targets 1 and 2
(c)
Shape estimates of cross targets
(d)
Shape estimates of the spawned targetFigure9 shows the number estimates and Figure 10 shows the accuracy statistic by the OSPA distance [18]. Notice that OSPA distance sharply increases at 56th and 66th seconds. The reason is that the targets make a cross with each other at 56th second and a spawned target appears. Generally, we can see that the proposed algorithm has a good performance for multiple extended target tracking with shape estimates.Figure 9
Number estimates.Figure 10
OSPA distance.


## 6. Conclusions

In this paper, a novel shape estimation algorithm is proposed based on the B-spline curve fitting. The multiple frame statistic technique is introduced to construct the pseudomeasurement sets. The selected control points are used to form the B-spline curve, and then the curve yields the shapes of the extended targets. Moreover, the proposed shape estimation algorithm is modified suitably and combined with the probability hypothesis density (PHD) filter for multiple extended target tracking. Simulations show that the proposed algorithm has a good performance for shape estimates of extended targets.



---

*Source: 741892-2014-06-25.xml*