# Ellipse Fitting Based Approach for Extended Object Tracking

Authors: Borui Li; Chundi Mu; Yongqiang Bai; Jianquan Bi; Lei Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

With the increase of sensors’ resolution, traditional object tracking technology, which ignores object’s physical extension, gradually becomes inappropriate. Extended object tracking (EOT) technology is able to obtain more information about the object through jointly estimating both centroid’s dynamic state and physical extension of the object. Random matrix based approach is a promising method for EOT. It uses ellipse/ellipsoid to describe the physical extension of the object. In order to reduce the physical extension estimation error when object maneuvers, the relationship between ellipse/ellipsoid and symmetrical positive definite matrix is analyzed at first. On this basis, ellipse/ellipsoid fitting based approach (EFA) for EOT is proposed based on the measurement model and centroid’s dynamic model of random matrix based EOT approach. Simulation results show that EFA is effective. The physical extension estimation error of EFA is lower than those of random matrix based approaches when object maneuvers. Besides, the estimation error of centroid’s dynamic state of EFA is also lower.


---

## Body

## 1. Introduction

In recent years, object tracking has always been the research focus in information fusion and many other fields [1, 2]. Physical extension of targets can be ignored if sensor’s resolution is relatively low. Traditional object tracking techniques usually ignore target’s physical extension and use centroid’s dynamic state to depict that of the entire target [3]. As the resolution of measurement sensors is gradually improved, this hypothesis becomes no longer suitable [4]. For example, when tracking large surface ships (such as aircraft carriers), each scan may produce several measurement data and the number of these data varies all the time. When some objects compose an intensive formation, such as flight formation in aerial refueling, the measurement sensors may fail to distinguish every individual object. These situations make it quite difficult to establish one-to-one correspondence between an object and a measurement. Traditional OT technology is no longer appropriate in these cases. Thus, it is necessary to measure and estimate object’s physical extension by extended object tracking (EOT) technique to improve the tracking performance [4–8]. Different from traditional OT, EOT estimates object’s physical extension and centroid’s dynamic state jointly. According to different modeling methods of the physical extension, EOT approaches include particle filter [9], probability hypothesis density filter [10, 11], and the approach based on random hypersurface model [12]. Compared with the above methods, random matrix based EOT approach (RMF) is simpler and has a smaller amount of computation, which makes it more practical and promising. RMF is first proposed by Koch [4] and uses ellipse/ellipsoid to describe the physical extension of object. This simplification can be applied in many scenarios, especially in the military fields, as shown in Figure 1.Describing physical extension of object with an ellipse: aircraft carrier (a) and flight formation (b).
(a)
(b)Because ellipse/ellipsoid can be expressed by symmetrical positive definite (SPD) matrix, the extension of object can be further depicted by SPD matrix. RMF assumes that the physical extension evolution is subject to Wishart-related distributions [4, 7]. The basic principle of RMF is to convert the joint posteriori probability distribution to the product of a Gaussian distribution and an inverse Wishart distribution within Bayesian framework. The filtering equations are obtained during the derivation. Feldmann et al. [7] improved measurement noise model and applied RMF in interactive multiple-model algorithm. Lan and Li [5, 6] further improved both of physical extension evolution model and measurement noise model. Li et al. [13] proposed the adaptive RMF based on model parameter adaption. According to the simulation results in literature, physical extension estimation error of RMF significantly increased when object is maneuvering, even utilizing multimodel algorithm [4–7, 13]. The “maneuvering” here refers to object’s actions that influence the physical extension, such as turning [4–7, 13] and significant change of object size.In order to improve EOT performance and reduce the physical extension estimation error during object’s maneuvering, the relationship between ellipse/ellipsoid and SPD matrix is analyzed in this paper. Then ellipse/ellipsoid fitting based approach (EFA) for EOT is proposed. hEFA can be integrated into a switching scheme, as shown in [8]. Utilizing “Givens rotation” [14] for SPD matrix decomposition, EFA divides the ellipse/ellipsoid fitting problem into two subproblems: direction angle estimation and semiaxis lengths estimation. The simulation results show that either the object is maneuvering or not; physical extension estimation error of EFA does not change much. When the object maneuvers, physical extension estimation error of EFA is lower than RMF.The object state is characterized by centroid’s dynamic statex
k and physical extension X
k simultaneously in EOT. Their estimation processes in EFA are mutually influenced and will be introduced successively in the sequel.


## 2. Centroid’s Dynamic State Estimation

Letd denote the spatial dimension, which is usually set as 2 or 3. s is the number of dimensions of centroid’s dynamic state in a spatial dimension. For example, when s
=
2, the centroid’s dynamic state includes position and velocity; when s
=
3, it also includes the acceleration [4–7, 13]. Thus, the centroid’s dynamic state x
k is a random s
d
×
1 vector, while the physical extension X
k   is a random d
×
d SPD matrix representing the ellipse/ellipsoid in d-dimension space.This section will introduce the dynamic model and measurement model ofx
k. Then the centroid’s dynamic state estimation x
k
|
k can be obtained by suitable estimators, for example, standard Kalman Filter.


### 2.1. Dynamic Model

The dynamic model ofx
k used in EOT is quite similar to that of standard Kalman filter [5–7, 13] as follows:
(1)
x
k
=
F
k
x
k
-
1
+
ω
k
,
where F
k
=
F
~
k
⊗
I
d is the state matrix, F
~
k is the state matrix in one spatial dimension, ⊗ is the right Kronecker product [15], and I
d is the d
×
d identity matrix. ω
k is the independent process noise, which is subject to the following Gaussian distribution:
(2)
ω
k
~
𝒩
(
0
,
D
~
k
⊗
I
d
)
,
where D
~
k is the covariance matrix of ω
k in a spatial dimension and 𝒩
(
μ
,
Σ
) is the Gaussian distribution with the expectation of μ and the covariance matrix of Σ.


### 2.2. Measurement Model

At scank, the measurement sensor obtains measurement data z
k
r
(
r
=
1,2
,
…
,
n
k
) with n
k
≥
1, where n
k is independent fromX
k   andx
k and subject to Poisson distribution with mean n
P. The real distribution of these measurement data is set as uniform distribution over the object, which is depicted by ellipse/ellipsoid, with the Gaussian noise of the measurement sensor superimposed [5, 7] as follows:
(3)
z
k
r
~
𝒰
(
X
k
)
+
𝒩
(
0
,
Ψ
k
)
.Measurement model (3) is hard to be directly used in the Bayesian framework or Kalman-related filter. Therefore, the following model of a single measurement is used for approximating [4, 7]:
(4)
z
k
r
=
H
k
x
k
+
υ
k
r
,
where H
k
=
H
~
k
⊗
I
d and H
~
k is the measurement matrix in one spatial dimension. Let 0
s
-
1   be (
s
-
1
)
×
1 zero vector and H
~
k
=
[
1
,
0
s
-
1
]; then, only the position of centroid is measured. υ
k
r is independent from υ
k
j
(
j
=
1,2
,
…
,
n
k
,
j
≠
r
) and is subject to the Gaussian distribution 𝒩
(
0
,
λ
X
k
+
Ψ
k
), whereλ is a scalar constant. λ can be set as 0.25 [7] via equaling the covariance matrices of z
k
r in (3) and (4). Since the probability density of the Gaussian distribution near the expectation is relatively large, only equaling two covariance matrices would lead to excessive measurement data near the ellipse’s center. Therefore, λ can be set as a relatively larger value so that (4) could be closer to the real measurement model as shown in (3).Furthermore, the simplified centroid’s measurement model for estimation ofx
k in EFA is as follows:
(5)
z
-
k
=
H
k
x
k
+
υ
k
,
where H
k is the same as that in (4) and
(6)
z
-
k
=
1
n
k
∑
r
=
1
n
k
z
k
r
.In (5), υ
k
~
𝒩
(
0
,
λ
X
k
-
1
|
k
-
1
+
Ψ
k
) is the simplified measurement noise, and X
k
-
1
|
k
-
1 is the estimation of the physical extension at scan k
-
1
(
k
>
1
), which is used to approximateX
k. Equation (5) can be interpreted as taking the mean value of n
k measurement data as the final measurement result at scan k. It indicates that the distribution of measurement data is influenced by object’s physical extension X
k, which also means that the estimation process of centroid’s dynamic state is influenced by physical extension.


## 2.1. Dynamic Model

The dynamic model ofx
k used in EOT is quite similar to that of standard Kalman filter [5–7, 13] as follows:
(1)
x
k
=
F
k
x
k
-
1
+
ω
k
,
where F
k
=
F
~
k
⊗
I
d is the state matrix, F
~
k is the state matrix in one spatial dimension, ⊗ is the right Kronecker product [15], and I
d is the d
×
d identity matrix. ω
k is the independent process noise, which is subject to the following Gaussian distribution:
(2)
ω
k
~
𝒩
(
0
,
D
~
k
⊗
I
d
)
,
where D
~
k is the covariance matrix of ω
k in a spatial dimension and 𝒩
(
μ
,
Σ
) is the Gaussian distribution with the expectation of μ and the covariance matrix of Σ.


## 2.2. Measurement Model

At scank, the measurement sensor obtains measurement data z
k
r
(
r
=
1,2
,
…
,
n
k
) with n
k
≥
1, where n
k is independent fromX
k   andx
k and subject to Poisson distribution with mean n
P. The real distribution of these measurement data is set as uniform distribution over the object, which is depicted by ellipse/ellipsoid, with the Gaussian noise of the measurement sensor superimposed [5, 7] as follows:
(3)
z
k
r
~
𝒰
(
X
k
)
+
𝒩
(
0
,
Ψ
k
)
.Measurement model (3) is hard to be directly used in the Bayesian framework or Kalman-related filter. Therefore, the following model of a single measurement is used for approximating [4, 7]:
(4)
z
k
r
=
H
k
x
k
+
υ
k
r
,
where H
k
=
H
~
k
⊗
I
d and H
~
k is the measurement matrix in one spatial dimension. Let 0
s
-
1   be (
s
-
1
)
×
1 zero vector and H
~
k
=
[
1
,
0
s
-
1
]; then, only the position of centroid is measured. υ
k
r is independent from υ
k
j
(
j
=
1,2
,
…
,
n
k
,
j
≠
r
) and is subject to the Gaussian distribution 𝒩
(
0
,
λ
X
k
+
Ψ
k
), whereλ is a scalar constant. λ can be set as 0.25 [7] via equaling the covariance matrices of z
k
r in (3) and (4). Since the probability density of the Gaussian distribution near the expectation is relatively large, only equaling two covariance matrices would lead to excessive measurement data near the ellipse’s center. Therefore, λ can be set as a relatively larger value so that (4) could be closer to the real measurement model as shown in (3).Furthermore, the simplified centroid’s measurement model for estimation ofx
k in EFA is as follows:
(5)
z
-
k
=
H
k
x
k
+
υ
k
,
where H
k is the same as that in (4) and
(6)
z
-
k
=
1
n
k
∑
r
=
1
n
k
z
k
r
.In (5), υ
k
~
𝒩
(
0
,
λ
X
k
-
1
|
k
-
1
+
Ψ
k
) is the simplified measurement noise, and X
k
-
1
|
k
-
1 is the estimation of the physical extension at scan k
-
1
(
k
>
1
), which is used to approximateX
k. Equation (5) can be interpreted as taking the mean value of n
k measurement data as the final measurement result at scan k. It indicates that the distribution of measurement data is influenced by object’s physical extension X
k, which also means that the estimation process of centroid’s dynamic state is influenced by physical extension.


## 3. Physical Extension Estimation

Since each scan would produce several measurement data in EOT, the basic idea of physical extension estimation in EFA is to treat the position components of these measurement data as points ind-dimension space. These measurement points could be used for ellipse/ellipsoid fitting, and then the estimated value of object’s physical extension is obtained based on the relationship between ellipse/ellipsoid and SPD matrix.In this section, EFA firstly uses “Givens rotations” [14] to decompose SPD matrix and analyze the relationship between ellipse/ellipsoid and SPD matrix. On this basis, ellipse/ellipsoid fitting problem is divided into the estimation of direction angle and semiaxis lengths, relatively.As most cases in EOT research setd as 2 [4–7, 13], we take d
=
2 as the example for introduction. The principle of EFA when d
=
3 is analogous to that when d
=
2.


### 3.1. Ellipse/Ellipsoid and SPD Matrix

Based on the dynamic model and measurement model ofx
k defined in Section 2, the equation of ellipse/ellipsoid which depicts the object is as follows [5]:
(7)
(
x
→
-
H
k
x
k
)
T
X
k
-
1
(
x
→
-
H
k
x
k
)
=
1
,
where d
×
1 vector x
→ is the coordinates of points on the ellipse/ellipsoid.To facilitate the analysis of the relationship between ellipse/ellipsoid and SPD matrix, the following lemma is introduced at first.Lemma 1.
d
×
d SPD matrix A can be converted into a diagonal matrix by “Givens rotations” with the maximal times ofd
(
d
-
1
)
/
2. The diagonal elements of the diagonal matrix are the eigenvalues of  A [14].Whend
=
2, according to Lemma 1, X
k could be decomposed into
(8)
X
k
=
R
(
φ
k
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
)
R
T
(
φ
k
)
,
where R
(
φ
k
)
=
[
cos
⁡
φ
k
-
sin
φ
k
sin
φ
k
cos
⁡
φ
k
] is rotation matrix, σ
k
,
1
2 and σ
k
,
2
2 are eigenvalues of X
k, and
(9)
φ
k
=
arctan
(
-
ρ
±
1
+
ρ
2
)
,
ρ
=
X
k
(
1,1
)
-
X
k
(
2,2
)
2
X
k
(
1,2
)
.An SPD matrix can be decomposed into the form shown in (8) by singular value decomposition, eigenvalue decomposition, and other means. However, these means could only ensure that R
(
φ
k
) is an orthogonal matrix. By Givens rotations, it could be guaranteed that R
(
φ
k
) is a rotation matrix. If the determinant of an orthogonal matrix is +1, then it is a rotation matrix.With (7) in mind, the following three elements are needed to completely identify an ellipse/ellipsoid in d-dimension space:(a)
coordinates of the center point, which are identified by position component ofx
k,
(b)
direction angle, which is identified by the rotation angle of “Givens rotations” in Lemma1: d
(
d
-
1
)
/
2   direction angles are needed at most to describe an ellipse/ellipsoid in d-dimension space,
(c)
the semiaxis lengths are equal to the square roots ofX
k’s eigenvalues.Therefore,X
k and x
k could jointly depict any ellipse/ellipsoid in d-dimension space. When the direction angles and semiaxis lengths of ellipse/ellipsoid are identified, X
k could be uniquely determined according to Lemma 1.


### 3.2. Direction Angle Estimation

Considering thatX
k could depict all of the ellipses/ellipsoids centering at the origin, according to (8), any ellipse X
k on a 2-dimension (2D) plane could be obtained by rotation of ellipse depicted by X
0
=
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
), where σ
k
,
1
2 and σ
k
,
2
2 are eigenvalues of X
k. X
0 is regarded as the “reference ellipse” in EFA, whose axes are parallel to two coordinate axes, relatively. If we assume that σ
k
,
1
2
≥
σ
k
,
2
2
>
0, then the direction angle φ
k in (8) is the intersection angle between the major axis of ellipseX
k   and the x-axis, while the counter-clockwise direction is set as the positive direction. Therefore, the method to estimate ellipse’s direction angle φ
k when d
=
2 is as follows.(a)
Centralize measurement data:(10)
z
k
r
,
*
=
z
k
r
-
H
k
x
k
|
k
,
r
=
1,2
,
…
,
n
k
.
(b)
Calculate the covariance matrixΣ
z
* of z
k
r
,
*, r
=
1,2
,
…
,
n
k.
(c)
Solve the 2 × 1 feature vectorr
→
z
,
max
⁡ corresponding to the larger eigenvalue of Σ
z
*. Then the intersection angle between r
→
z
,
max
⁡ and x-axis is the estimated value of φ
k as follows:
(11)
φ
^
k
=
arctan
(
r
→
z
,
max
⁡
(
2
)
r
→
z
,
max
⁡
(
1
)
)
.Step (a) of direction angel estimation indicates that physical extension estimation is influenced by centroid’s dynamic state. According to the principle of principal component analysis [16], r
→
z
,
max
⁡   is the direction of the largest variance of z
k
r
,
*. In addition, when utilizing the measurement model shown in (3), the direction of X
k’s major axis could be approximately considered the direction with the greatest changes of z
k
r
,
*. Therefore, EFA uses r
→
z
,
max
⁡   as the estimated value of the major axis direction of X
k, which is also the estimation of direction angle φ
k.


### 3.3. Semiaxis Lengths Estimation

The method of estimating semiaxis lengths in EFA is divided into two steps as follows:(a)
identifying the upper and lower bounds of semiaxis lengths;
(b)
selecting a suitable step size that is selected for traversing, and then the semiaxis lengths of an ellipse are estimated according to corresponding criterion.Whend
=
2, the semiaxis lengths of ellipseX
k are σ
k
,
1
≥
σ
k
,
2
>
0. Their upper bounds can be determined as follows:(a)
z
k
r
,
* is projected to the directions of two feature vectors of Σ
z
*;
(b)
the maximal value of absolute projected coordinates could be taken as the upper boundsσ
k
,
1
max
⁡ and σ
k
,
2
max
⁡, respectively.The lower bounds of semiaxis lengths could be identified by the covariance matrixΨ
k of measurement sensor’s noise. Let σ
Ψ
,
a
2
≥
σ
Ψ
,
b
2
>
0 denote the eigenvalues of Ψ
k. Then the lower bound of semiaxis lengths could be set as σ
k
,
1
min
⁡
=
σ
Ψ
,
a and σ
k
,
2
min
⁡
=
σ
Ψ
,
b.Once the upper and lower bounds of semiaxis lengths are identified, the method to estimate semiaxis lengths is as follows.(a)
Determine traversal step size according to the required accuracy and operating speed.
(b)
Starting fromσ
k
,
1
min
⁡ and σ
k
,
2
min
⁡, each pair of semiaxis lengths candidates σ
k
,
1 and σ
k
,
2 that meets the conditions of (12) is traversed. And the corresponding estimated X
k
|
k is calculated according to (8) as follows:
(12)
σ
k
,
1
≥
σ
k
,
2
σ
k
,
1
∈
[
σ
k
,
1
min
⁡
,
σ
k
,
1
max
⁡
]
,
σ
k
,
2
∈
[
σ
k
,
2
min
⁡
,
σ
k
,
2
max
⁡
]
.
(c)
Calculaten
den based on (7), which is the number of measurement points on and inside the ellipse.
(d)
Ifn
den
/
n
k
>
t
den, the density of measurement points is calculated as ξ
den
∝
n
den
/
σ
k
,
1
σ
k
,
2, which is the number of measurement points in the unit area of ellipse, and 0
<
t
des
⁡
<
1 is a preset threshold; else, ξ
den   is set to an infinitesimal.
(e)
If there is no candidate value of semiaxis lengths that meets (12), advance to Step (f); else, return to Step (b).
(f)
The candidate value of semiaxis lengths which has the largest correspondingξ
den   is the estimated valueσ
^
k
,
1   andσ
^
k
,
2.In Step (d),ξ
den could be replaced by the lower bound of the Wilson interval [17] as follows:
(13)
ξ
wilson
=
(
ξ
~
den
+
κ
2
2
n
den
-
κ
ξ
~
den
(
1
-
ξ
~
den
)
n
den
+
κ
2
4
n
den
2
)
×
(
1
+
κ
2
n
den
)
-
1
,
where κ is set as 1.96 under the confidence level of 95% and ξ
~
den is the normalization result of all valid ξ
den. When using ξ
den, the improper situation that n
des
⁡ and σ
k
,
1
σ
k
,
2 are all quite small but that ξ
den is large is easy to occur. ξ
wilson simultaneously considers n
den and ξ
den; that is, it is required that the density is maximized while the number of measurement points on and inside the ellipse should not be too small. Thus, the aforementioned improper situation could be effectively avoided. When n
k is large, the effect of ξ
wilson and ξ
den is more or less equivalent. Since n
k in EOT is usually not very large, ξ
wilsonwould perform better than ξ
den in EFA, yet resulting in more complex calculations.To take into account higher precision and faster calculation speed simultaneously, we could use a relatively larger step size to obtainσ
^
k
,
1 and σ
^
k
,
2, and then repeat the above algorithm with a smaller step size in the interval near σ
^
k
,
1 and σ
^
k
,
2. Thus, a more efficient estimation result can be achieved through a relatively smaller amount of computation.After calculating the direction angle and semiaxis lengths, the estimated valueX
k
|
k of X
k could be obtained according to (8).


### 3.4. Physical Extension Estimation Whend
=
3

Whend
=
3, according to Lemma 1, we have
(14)
X
k
=
R
3
(
φ
k
,
3
)
R
2
(
φ
k
,
2
)
R
1
(
φ
k
,
1
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
,
σ
k
,
3
2
)
×
[
R
3
(
φ
k
,
3
)
R
2
(
φ
k
,
2
)
R
1
(
φ
k
,
1
)
]
T
,
where R
·
(
φ
k
,
·
)   is the rotation matrix having counterclockwise rotation around a certain coordinate axis. For example, by counterclockwise rotation around x-axis, it should be [
1
0
0
0
cos
⁡
φ
k
,
·
-
sin
φ
k
,
·
0
sin
φ
k
,
·
cos
⁡
φ
k
,
·
].Equation (14) shows that maximally three Givens rotations are needed to convert X
k to a diagonal matrix when d
=
3. Similar to that when d
=
2, we may assume that σ
k
,
1
2
≥
σ
k
,
2
2
≥
σ
k
,
3
2
>
0 and adopt the following sequence of rotation matrix:
(15)
X
k
=
R
z
(
φ
k
,
3
)
R
y
(
φ
k
,
2
)
R
x
(
φ
k
,
1
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
,
σ
k
,
3
2
)
×
[
R
z
(
φ
k
,
3
)
R
y
(
φ
k
,
2
)
R
x
(
φ
k
,
1
)
]
T
,
where R
x
(
φ
k
,
1
)is the rotation matrix with the angleφ
k
,
1   around x-axis at the counterclockwise direction, which is the same with R
y
(
φ
k
,
2
) and R
z
(
φ
k
,
3
). Therefore, the ellipsoid’s direction angles φ
k
,
j
(
j
=
1,2
,
3
) could be estimated as follows.(a)
Centralize measurement data according to (10).
(b)
Projectz
k
r
,
* to the y-z 2D plane. The estimation method of direction angle when d
=
2 is used to obtain φ
^
k
,
1.
(c)
Similarly, projectz
k
r
,
* to 2D plane of x-z and x-y; then, calculate φ
^
k
,
2 and φ
^
k
,
3, respectively.Whend
=
3, the semiaxis lengths estimation method with d
=
2 is completely suitable. The difference is that ξ
den
∝
n
den
/
σ
k
,
1
σ
k
,
2
σ
k
,
3 represents the number of measurement points in the unit volume of the ellipsoid.Once the direction angles and semiaxis lengths are obtained, the estimated valueX
k
|
k could be calculated according to (15).


## 3.1. Ellipse/Ellipsoid and SPD Matrix

Based on the dynamic model and measurement model ofx
k defined in Section 2, the equation of ellipse/ellipsoid which depicts the object is as follows [5]:
(7)
(
x
→
-
H
k
x
k
)
T
X
k
-
1
(
x
→
-
H
k
x
k
)
=
1
,
where d
×
1 vector x
→ is the coordinates of points on the ellipse/ellipsoid.To facilitate the analysis of the relationship between ellipse/ellipsoid and SPD matrix, the following lemma is introduced at first.Lemma 1.
d
×
d SPD matrix A can be converted into a diagonal matrix by “Givens rotations” with the maximal times ofd
(
d
-
1
)
/
2. The diagonal elements of the diagonal matrix are the eigenvalues of  A [14].Whend
=
2, according to Lemma 1, X
k could be decomposed into
(8)
X
k
=
R
(
φ
k
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
)
R
T
(
φ
k
)
,
where R
(
φ
k
)
=
[
cos
⁡
φ
k
-
sin
φ
k
sin
φ
k
cos
⁡
φ
k
] is rotation matrix, σ
k
,
1
2 and σ
k
,
2
2 are eigenvalues of X
k, and
(9)
φ
k
=
arctan
(
-
ρ
±
1
+
ρ
2
)
,
ρ
=
X
k
(
1,1
)
-
X
k
(
2,2
)
2
X
k
(
1,2
)
.An SPD matrix can be decomposed into the form shown in (8) by singular value decomposition, eigenvalue decomposition, and other means. However, these means could only ensure that R
(
φ
k
) is an orthogonal matrix. By Givens rotations, it could be guaranteed that R
(
φ
k
) is a rotation matrix. If the determinant of an orthogonal matrix is +1, then it is a rotation matrix.With (7) in mind, the following three elements are needed to completely identify an ellipse/ellipsoid in d-dimension space:(a)
coordinates of the center point, which are identified by position component ofx
k,
(b)
direction angle, which is identified by the rotation angle of “Givens rotations” in Lemma1: d
(
d
-
1
)
/
2   direction angles are needed at most to describe an ellipse/ellipsoid in d-dimension space,
(c)
the semiaxis lengths are equal to the square roots ofX
k’s eigenvalues.Therefore,X
k and x
k could jointly depict any ellipse/ellipsoid in d-dimension space. When the direction angles and semiaxis lengths of ellipse/ellipsoid are identified, X
k could be uniquely determined according to Lemma 1.


## 3.2. Direction Angle Estimation

Considering thatX
k could depict all of the ellipses/ellipsoids centering at the origin, according to (8), any ellipse X
k on a 2-dimension (2D) plane could be obtained by rotation of ellipse depicted by X
0
=
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
), where σ
k
,
1
2 and σ
k
,
2
2 are eigenvalues of X
k. X
0 is regarded as the “reference ellipse” in EFA, whose axes are parallel to two coordinate axes, relatively. If we assume that σ
k
,
1
2
≥
σ
k
,
2
2
>
0, then the direction angle φ
k in (8) is the intersection angle between the major axis of ellipseX
k   and the x-axis, while the counter-clockwise direction is set as the positive direction. Therefore, the method to estimate ellipse’s direction angle φ
k when d
=
2 is as follows.(a)
Centralize measurement data:(10)
z
k
r
,
*
=
z
k
r
-
H
k
x
k
|
k
,
r
=
1,2
,
…
,
n
k
.
(b)
Calculate the covariance matrixΣ
z
* of z
k
r
,
*, r
=
1,2
,
…
,
n
k.
(c)
Solve the 2 × 1 feature vectorr
→
z
,
max
⁡ corresponding to the larger eigenvalue of Σ
z
*. Then the intersection angle between r
→
z
,
max
⁡ and x-axis is the estimated value of φ
k as follows:
(11)
φ
^
k
=
arctan
(
r
→
z
,
max
⁡
(
2
)
r
→
z
,
max
⁡
(
1
)
)
.Step (a) of direction angel estimation indicates that physical extension estimation is influenced by centroid’s dynamic state. According to the principle of principal component analysis [16], r
→
z
,
max
⁡   is the direction of the largest variance of z
k
r
,
*. In addition, when utilizing the measurement model shown in (3), the direction of X
k’s major axis could be approximately considered the direction with the greatest changes of z
k
r
,
*. Therefore, EFA uses r
→
z
,
max
⁡   as the estimated value of the major axis direction of X
k, which is also the estimation of direction angle φ
k.


## 3.3. Semiaxis Lengths Estimation

The method of estimating semiaxis lengths in EFA is divided into two steps as follows:(a)
identifying the upper and lower bounds of semiaxis lengths;
(b)
selecting a suitable step size that is selected for traversing, and then the semiaxis lengths of an ellipse are estimated according to corresponding criterion.Whend
=
2, the semiaxis lengths of ellipseX
k are σ
k
,
1
≥
σ
k
,
2
>
0. Their upper bounds can be determined as follows:(a)
z
k
r
,
* is projected to the directions of two feature vectors of Σ
z
*;
(b)
the maximal value of absolute projected coordinates could be taken as the upper boundsσ
k
,
1
max
⁡ and σ
k
,
2
max
⁡, respectively.The lower bounds of semiaxis lengths could be identified by the covariance matrixΨ
k of measurement sensor’s noise. Let σ
Ψ
,
a
2
≥
σ
Ψ
,
b
2
>
0 denote the eigenvalues of Ψ
k. Then the lower bound of semiaxis lengths could be set as σ
k
,
1
min
⁡
=
σ
Ψ
,
a and σ
k
,
2
min
⁡
=
σ
Ψ
,
b.Once the upper and lower bounds of semiaxis lengths are identified, the method to estimate semiaxis lengths is as follows.(a)
Determine traversal step size according to the required accuracy and operating speed.
(b)
Starting fromσ
k
,
1
min
⁡ and σ
k
,
2
min
⁡, each pair of semiaxis lengths candidates σ
k
,
1 and σ
k
,
2 that meets the conditions of (12) is traversed. And the corresponding estimated X
k
|
k is calculated according to (8) as follows:
(12)
σ
k
,
1
≥
σ
k
,
2
σ
k
,
1
∈
[
σ
k
,
1
min
⁡
,
σ
k
,
1
max
⁡
]
,
σ
k
,
2
∈
[
σ
k
,
2
min
⁡
,
σ
k
,
2
max
⁡
]
.
(c)
Calculaten
den based on (7), which is the number of measurement points on and inside the ellipse.
(d)
Ifn
den
/
n
k
>
t
den, the density of measurement points is calculated as ξ
den
∝
n
den
/
σ
k
,
1
σ
k
,
2, which is the number of measurement points in the unit area of ellipse, and 0
<
t
des
⁡
<
1 is a preset threshold; else, ξ
den   is set to an infinitesimal.
(e)
If there is no candidate value of semiaxis lengths that meets (12), advance to Step (f); else, return to Step (b).
(f)
The candidate value of semiaxis lengths which has the largest correspondingξ
den   is the estimated valueσ
^
k
,
1   andσ
^
k
,
2.In Step (d),ξ
den could be replaced by the lower bound of the Wilson interval [17] as follows:
(13)
ξ
wilson
=
(
ξ
~
den
+
κ
2
2
n
den
-
κ
ξ
~
den
(
1
-
ξ
~
den
)
n
den
+
κ
2
4
n
den
2
)
×
(
1
+
κ
2
n
den
)
-
1
,
where κ is set as 1.96 under the confidence level of 95% and ξ
~
den is the normalization result of all valid ξ
den. When using ξ
den, the improper situation that n
des
⁡ and σ
k
,
1
σ
k
,
2 are all quite small but that ξ
den is large is easy to occur. ξ
wilson simultaneously considers n
den and ξ
den; that is, it is required that the density is maximized while the number of measurement points on and inside the ellipse should not be too small. Thus, the aforementioned improper situation could be effectively avoided. When n
k is large, the effect of ξ
wilson and ξ
den is more or less equivalent. Since n
k in EOT is usually not very large, ξ
wilsonwould perform better than ξ
den in EFA, yet resulting in more complex calculations.To take into account higher precision and faster calculation speed simultaneously, we could use a relatively larger step size to obtainσ
^
k
,
1 and σ
^
k
,
2, and then repeat the above algorithm with a smaller step size in the interval near σ
^
k
,
1 and σ
^
k
,
2. Thus, a more efficient estimation result can be achieved through a relatively smaller amount of computation.After calculating the direction angle and semiaxis lengths, the estimated valueX
k
|
k of X
k could be obtained according to (8).


## 3.4. Physical Extension Estimation Whend
=
3

Whend
=
3, according to Lemma 1, we have
(14)
X
k
=
R
3
(
φ
k
,
3
)
R
2
(
φ
k
,
2
)
R
1
(
φ
k
,
1
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
,
σ
k
,
3
2
)
×
[
R
3
(
φ
k
,
3
)
R
2
(
φ
k
,
2
)
R
1
(
φ
k
,
1
)
]
T
,
where R
·
(
φ
k
,
·
)   is the rotation matrix having counterclockwise rotation around a certain coordinate axis. For example, by counterclockwise rotation around x-axis, it should be [
1
0
0
0
cos
⁡
φ
k
,
·
-
sin
φ
k
,
·
0
sin
φ
k
,
·
cos
⁡
φ
k
,
·
].Equation (14) shows that maximally three Givens rotations are needed to convert X
k to a diagonal matrix when d
=
3. Similar to that when d
=
2, we may assume that σ
k
,
1
2
≥
σ
k
,
2
2
≥
σ
k
,
3
2
>
0 and adopt the following sequence of rotation matrix:
(15)
X
k
=
R
z
(
φ
k
,
3
)
R
y
(
φ
k
,
2
)
R
x
(
φ
k
,
1
)
diag
⁡
(
σ
k
,
1
2
,
σ
k
,
2
2
,
σ
k
,
3
2
)
×
[
R
z
(
φ
k
,
3
)
R
y
(
φ
k
,
2
)
R
x
(
φ
k
,
1
)
]
T
,
where R
x
(
φ
k
,
1
)is the rotation matrix with the angleφ
k
,
1   around x-axis at the counterclockwise direction, which is the same with R
y
(
φ
k
,
2
) and R
z
(
φ
k
,
3
). Therefore, the ellipsoid’s direction angles φ
k
,
j
(
j
=
1,2
,
3
) could be estimated as follows.(a)
Centralize measurement data according to (10).
(b)
Projectz
k
r
,
* to the y-z 2D plane. The estimation method of direction angle when d
=
2 is used to obtain φ
^
k
,
1.
(c)
Similarly, projectz
k
r
,
* to 2D plane of x-z and x-y; then, calculate φ
^
k
,
2 and φ
^
k
,
3, respectively.Whend
=
3, the semiaxis lengths estimation method with d
=
2 is completely suitable. The difference is that ξ
den
∝
n
den
/
σ
k
,
1
σ
k
,
2
σ
k
,
3 represents the number of measurement points in the unit volume of the ellipsoid.Once the direction angles and semiaxis lengths are obtained, the estimated valueX
k
|
k could be calculated according to (15).


## 4. Simulation Results and Discussion

### 4.1. Scenario 1

Scenario 1, as shown in Figure2, is commonly used in RMF research [5, 7, 13]. The object is an aircraft carrier with the length of about 350 m and the width of about 100 m. An ellipse with the semiaxis lengths of 175 m and 50 m is used to depict its physical extension. The object starts from the origin and performs uniform linear motion along the path. Then it makes three turns as shown in Figure 2. The velocity of object is set as 27 knots (about 50 km/h). The measurement period is set as T
=
10
s and the covariance matrix of measurement sensor’s noise is Ψ
k
=
diag
⁡
(
4
0
2
,
1
5
2
) m2. The number of measurement points of each scan is subject to Poisson distribution with the mean value of n
p = 25.Figure 2
Object trajectory in Scenario 1.The following four EOT approaches are compared throughN
MC
=
600 Monte Carlo simulations:(a)
Koch’s RMF [4]: single-model RMF approach and time decay constant τ is set as 8
T;
(b)
Feldmann’s RMF-MM [7]: multimodel RMF approach and parameters are set according to the literature [7];
(c)
EFA-den:ξ
den   is selected for semiaxis lengths estimation, and threshold value t
den is set as 0.5;
(d)
EFA-wilson:ξ
wilson is selected for semiaxis lengths estimation, and threshold value t
den is set as 0.5.To compare the performance of physical extension estimation, root mean square error (RMSE) of physical extension is defined as follows [7]:
(16)
RMS
E
X
=
∑
h
=
1
N
MC
tr
⁡
[
(
X
k
|
k
h
-
X
k
)
2
]
N
MC
,
where X
k
|
k
h is the estimated value of physical extension at scan k in hth Monte Carlo simulation.Physical extension RMSE curves in Scenario 1 are shown in Figure3. When the object is maneuvering (making three turns), physical extension RMSEs of RMFs increase while those of EFAs do not change significantly. Physical extension RMSEs of EFAs are lower than Koch’s RMF but higher than Feldmann’s RMF-MM when the object is not maneuvering. Physical extension RMSE of EFA-wilson is lower than that of EFA-den.Figure 3
RMSEs of physical extension in Scenario 1.The centroid’s position and velocity RMSE curves in Scenario 1 are shown in Figures4 and 5. The position and velocity RMSEs of two EFAs are nearly the same and lower than two RMFs.Figure 4
RMSEs of centroid’s position in Scenario 1.Figure 5
RMSEs of centroid’s velocity in Scenario 1.


### 4.2. Scenario 2

In Scenario 2, the trajectory of object is shown in Figure6. From the origin, the object takes uniform linear motion for some time and then keeps taking uniform circular motion. Four EOT approaches and other simulation parameters for comparison are the same as those in Scenario 1.Figure 6
Object trajectory in Scenario 2.RMSE curves of physical extension are shown in Figure7. In the beginning, the object is not maneuvering, and physical extension RMSE of Feldmann’s RMF-MM is lower than EFAs. When the object is taking uniform circular motion (maneuvering), the physical extension RMSEs of two EFAs become lower than two RMF approaches.Figure 7
RMSEs of physical extension in Scenario 2.RMSE curves of centroid’s position and velocity are shown in Figures8 and 9. Similar to Scenario 1, the position and velocity RMSEs of two EFAs are quite close to each other and lower than the RMFs.Figure 8
RMSEs of centroid’s position in Scenario 2.Figure 9
RMSEs of centroid’s velocity in Scenario 2.


## 4.1. Scenario 1

Scenario 1, as shown in Figure2, is commonly used in RMF research [5, 7, 13]. The object is an aircraft carrier with the length of about 350 m and the width of about 100 m. An ellipse with the semiaxis lengths of 175 m and 50 m is used to depict its physical extension. The object starts from the origin and performs uniform linear motion along the path. Then it makes three turns as shown in Figure 2. The velocity of object is set as 27 knots (about 50 km/h). The measurement period is set as T
=
10
s and the covariance matrix of measurement sensor’s noise is Ψ
k
=
diag
⁡
(
4
0
2
,
1
5
2
) m2. The number of measurement points of each scan is subject to Poisson distribution with the mean value of n
p = 25.Figure 2
Object trajectory in Scenario 1.The following four EOT approaches are compared throughN
MC
=
600 Monte Carlo simulations:(a)
Koch’s RMF [4]: single-model RMF approach and time decay constant τ is set as 8
T;
(b)
Feldmann’s RMF-MM [7]: multimodel RMF approach and parameters are set according to the literature [7];
(c)
EFA-den:ξ
den   is selected for semiaxis lengths estimation, and threshold value t
den is set as 0.5;
(d)
EFA-wilson:ξ
wilson is selected for semiaxis lengths estimation, and threshold value t
den is set as 0.5.To compare the performance of physical extension estimation, root mean square error (RMSE) of physical extension is defined as follows [7]:
(16)
RMS
E
X
=
∑
h
=
1
N
MC
tr
⁡
[
(
X
k
|
k
h
-
X
k
)
2
]
N
MC
,
where X
k
|
k
h is the estimated value of physical extension at scan k in hth Monte Carlo simulation.Physical extension RMSE curves in Scenario 1 are shown in Figure3. When the object is maneuvering (making three turns), physical extension RMSEs of RMFs increase while those of EFAs do not change significantly. Physical extension RMSEs of EFAs are lower than Koch’s RMF but higher than Feldmann’s RMF-MM when the object is not maneuvering. Physical extension RMSE of EFA-wilson is lower than that of EFA-den.Figure 3
RMSEs of physical extension in Scenario 1.The centroid’s position and velocity RMSE curves in Scenario 1 are shown in Figures4 and 5. The position and velocity RMSEs of two EFAs are nearly the same and lower than two RMFs.Figure 4
RMSEs of centroid’s position in Scenario 1.Figure 5
RMSEs of centroid’s velocity in Scenario 1.


## 4.2. Scenario 2

In Scenario 2, the trajectory of object is shown in Figure6. From the origin, the object takes uniform linear motion for some time and then keeps taking uniform circular motion. Four EOT approaches and other simulation parameters for comparison are the same as those in Scenario 1.Figure 6
Object trajectory in Scenario 2.RMSE curves of physical extension are shown in Figure7. In the beginning, the object is not maneuvering, and physical extension RMSE of Feldmann’s RMF-MM is lower than EFAs. When the object is taking uniform circular motion (maneuvering), the physical extension RMSEs of two EFAs become lower than two RMF approaches.Figure 7
RMSEs of physical extension in Scenario 2.RMSE curves of centroid’s position and velocity are shown in Figures8 and 9. Similar to Scenario 1, the position and velocity RMSEs of two EFAs are quite close to each other and lower than the RMFs.Figure 8
RMSEs of centroid’s position in Scenario 2.Figure 9
RMSEs of centroid’s velocity in Scenario 2.


## 5. Conclusions

Extended object tracking is able to estimate centroid’s dynamic state and physical extension of objects jointly and provide more information about the object. Therefore, it has become a research focus of information fusion and maneuvering object tracking in recent years. RMF uses simple models and could obtain iterative filtering equations within the Bayesian framework. Its main drawback is the significant rise of physical extension estimation error during object’s maneuvering. In order to overcome the shortcoming of RMF, an EOT approach based on ellipse/ellipsoid fitting is proposed on the basis of the analysis of the internal correlation between ellipse/ellipsoid and SPD matrix. The simulation results show that EFA is insensitive to whether the object is maneuvering or not. Physical extension estimation error of EFA is lower than the two RMFs during object’s maneuvering. Centroid’s dynamic state estimation error of EFA is lower than RMF. EFA-wilson performs better than EFA-den in physical extension estimation.



---

*Source: 632815-2014-05-08.xml*