Validated Uncertainty Propagation for Estimation and Measure
Association, Application to Satellite Tracking
Charlotte Govignon
1 a
, Elliot Brendel
2 b
and Julien Alexandre Dit Sandretto
3 c
1
´
Ecole des Ponts ParisTech, Cit
´
e Descartes, 6 et 8 avenue Blaise Pascal, 77420 Champs-sur-Marne, France
2
Thales Land and Air Systems, Paris, France
3
ENSTA Paris, Institut Polytechnique de Paris, 828 Boulevard des Mar
´
echaux, Palaiseau 91120, France
Keywords:
Estimation, Association, Intervals, Contractors, Satellite Tracking, Radar.
Abstract:
In this paper, we present a new uncertainty propagation algorithm based on interval arithmetic, and its appli-
cations for space surveillance. Using validated simulation, the goal is to explore the benefits of a set-based
approach of estimation and data association for satellite tracking. The presented algorithm capitalises on the
position measures of a satellite to improve its estimation and reduce uncertainties on its trajectory. Our ap-
proach also contributes to data association by computing the precision required for a measure to belong to a
given track with confidence-levels. This paper illustrates the contributions of this new algorithm with several
scenarios of orbit determination and satellite tracking and their numerical simulations.
1 INTRODUCTION
Space surveillance of low-earth objects has become a
major issue, notably in a military context, and can be
addressed by ground radars, providing detection and
tracking of satellites of interest. In order to perform
the tracking of these objects, two crucial steps must
be accomplished: estimation and association.
The estimation problem consists in using the in-
formation provided by a set of measures and their re-
spective uncertainties to compute the state of a dy-
namical system, for instance position and velocity
coordinates, while minimizing the state uncertainty.
Within this context, the simulation of the considered
dynamical system is often needed and can be used to
improve the estimation accuracy. The simulation of
a dynamical system usually requires the computing
of the solution of a differential equation represent-
ing its dynamics. An other issue in the estimation
problem occurs when multiple systems evolve in the
observed environment. In that case, the sensor pro-
viding the measures can produce data coming from
different sources. Hence, it is essential to filter the
measures between those coming from the considered
system and the others, in order to prevent the estimate
a
https://orcid.org/0009-0003-7703-0785
b
https://orcid.org/0000-0002-0458-4993
c
https://orcid.org/0000-0002-6185-2480
from corruption. This issue is called association.
Data association techniques can be divided in two
categories. First, instantaneous decision techniques
which assign a measure to an object at the time when
the measure is generated. In these techniques, Near-
est Neighbor (NN) (Blackman and Popoli, 1999) ap-
proaches consist in assigning the considered measure
to its closest object by computing a chosen quantity
representing a distance between the measure and its
surrounding objects. Such quantities can be deter-
ministic (e.g. the Cartesian distance) or statistical
(e.g. the Mahalanobis distance (Verhagen and Teu-
nissen, 2017)). When multiple measures and multi-
ple objects are considered at the same time, Global
Nearest Neighbor (GNN) approaches (Blackman and
Popoli, 1999) use optimization techniques such as the
Hungarian method (Mills-Tettey et al., 2007) to max-
imize association scores between measures and ob-
jects. Usually, NN and GNN’s efficiency decreases
when the number of measures and objects increases.
Probabilistic Data Association (PDA) (Blackman
and Popoli, 1999; Bar-Shalom and Tse, 1975; Fort-
mann et al., 1983) approaches are instantaneous de-
cision techniques which consist in the computation
of the probability of each association hypothesis and
the selection of a subset of these hypotheses. Then,
a combination of the selected hypotheses is evalu-
ated, assigned to an object and used to improve its
660
Govignon, C., Brendel, E. and Alexandre Dit Sandretto, J.
Validated Uncertainty Propagation for Estimation and Measure Association, Application to Satellite Tracking.
DOI: 10.5220/0012939400003822
Paper published under CC license (CC BY-NC-ND 4.0)
In Proceedings of the 21st International Conference on Informatics in Control, Automation and Robotics (ICINCO 2024) - Volume 1, pages 660-667
ISBN: 978-989-758-717-7; ISSN: 2184-2809
Proceedings Copyright © 2024 by SCITEPRESS – Science and Technology Publications, Lda.

estimate. This approach is usually suited for envi-
ronments with risk of confusion between multiple ob-
jects. On the other hand, delayed decision techniques
wait until several measures are generated to perform
an association taking into account the information
of the batch of the measures. For example, Multi-
Hypothesis Tracking (MHT) (Blackman and Popoli,
1999; Reid, 1979) approaches consist in building and
updating a set of association hypotheses until their re-
spective estimated probabilities lead to the removal of
one of the previously considered hypotheses. These
techniques are efficient but show high combinatorial
cost. Furthermore, the long time needed to take a de-
cision can be a drawback in operating environment.
Other delayed decision approaches consist in the
combination of gating (Collins and Uhlmann, 1992)
and estimation (Verhagen and Teunissen, 2017; Val-
lado, 2001; Julier and Uhlmann, 2004) methods. Gat-
ing consist in selecting a subset of measures verifying
a condition (e.g. all the measures whose distance with
the estimate is under a chosen threshold). Then, esti-
mation techniques, such as Extended Kalman Filter
(Blackman and Popoli, 1999; Verhagen and Teunis-
sen, 2017; Vallado, 2001; Julier and Uhlmann, 2004)
or Least Squares optimization (Verhagen and Teunis-
sen, 2017; Vallado, 2001; Gill and Murray, 1978) are
used to update the estimate with the subset of selected
measures. At last, the new estimate is used to realize
a final gating step, and the new selected measures are
considered as assigned to the object.
Most of the described association techniques are
probabilistic approaches, and they are often based
on the hypothesis that Gaussian distributions describe
well enough the several uncertainties (from the mea-
sures and the estimates) all along the trajectories of
the considered objects. Yet, a ground radar can only
provide measures for low-earth satellites when these
objects are in its field of view. Therefore, due to earth
rotation and satellite dynamics (see 2.2), measures
of an observed satellite can be spaced out by several
hours, causing a loss of precision in the estimate of
the satellite state, leading to increased uncertainties
and, potentially, confusion with other objects. More-
over, it has been documented that a Gaussian distri-
bution propagated along an orbital trajectory can lose
its Gaussian nature (Horwood and Poore, 2014), in
particular with large uncertainties, for example when
the object generated a limited number of measures.
In that case, the previously described techniques be-
come unsuitable. Conversely, deterministic set-based
propagation and association approaches, using classi-
cal set theory functions such as intersection, inclusion
and contractors, provide a guaranteed enclosure of all
the possible values of the objects states. The obtained
enclosures can be used to decide, with guarantee, if a
measure belongs to an object or not. Such techniques
are therefore relevant in a space surveillance context.
The approaches developed in this paper are applied to
association, estimation and tracking of satellites ob-
served by a ground radar. Performances and bene-
fits of these approaches are discussed in connection
to this space surveillance application.
2 ORBIT DETERMINATION
2.1 Orbit Description
The state of a satellite is usually given with its vectors
of position and speed. These six coordinates can be
expressed in earth-centered Cartesian references such
as Earth-Centered Intertial (ECI) or Earth-Centered
Earth-Fixed (ECEF). Alternatively, it can be con-
venient to describe a satellite state in the Keple-
rian or the Equinoctial frames (Vallado, 2001), es-
pecially when the dynamics are unperturbed, since
these frames provide a more direct handling of the
parameters of the ellipse. Modified equinoctial ele-
ments (Walker et al., 1985) (mEOE) are derived from
Keplerian elements (a, e, i, Ω, ω, θ) as follows:



p = a

1 −e
2

, L = Ω + ω + θ
f = e cos(ω + Ω), g = e sin (ω + Ω)
h = tan

i
2

cosΩ, k = tan

i
2

sinΩ
(1)
where a denotes the semi-major axis, e the eccentric-
ity, i the inclination, Ω the longitude of the ascend-
ing node, ω the argument of the periapsis, θ the true
anomaly. The element p is called the semi-parameter
and the element L the true longitude. This reference is
useful for trajectory analysis because of the stability
of the state components in time.
2.2 Dynamics of a Satellite
The dynamics of a satellite in mEOE is given by
Equation (2) (Walker et al., 1985):





















˙p =
2p
w
q
p
µ
∆
t
˙
f =
q
p
µ

((w+1)cos L+ f )∆
t
−g j∆
n
w
+ ∆
r
sinL

˙g =
q
p
µ

((w+1)sin L+g)∆
t
+ f j∆
n
w
−∆
r
cosL

˙
h =
q
p
µ
s
2
∆
n
2w
cosL,
˙
k =
q
p
µ
s
2
∆
n
2w
sinL
˙
L =
√
µp
r
2
+
1
w
q
p
µ
(hsinL −k cos L)∆
n
(2)
where w = 1 + f cos L + g sin L, s
2
= 1 + h
2
+ k
2
, and
j = h sin L −k cosL. ∆
t
, ∆
r
and ∆
n
are the non-two-
Validated Uncertainty Propagation for Estimation and Measure Association, Application to Satellite Tracking
661

body perturbations in the radial, tangential and nor-
mal directions, respectively. These perturbations ag-
gregate all other orbital perturbations acting on the
satellite, like the J
k
effects (k ∈ {2, 3, ...}) which
come from the non-sphericity of the Earth, or the at-
mospheric drag for instance. With only J
2
effects,
which induce the strongest perturbation from the non-
sphericity of the Earth, the perturbations are:









∆
J
2
r
= −
3µJ
2
R
2
e
2r
4

1 −
12 j
2
s
4

∆
J
2
t
= −
12µJ
2
R
2
e
r
4

j(h cosL−k sin L)
s
4

∆
J
2
n
= −
6µJ
2
R
2
e
r
4

(1−h
2
−k
2
) j
s
4

(3)
where J
2
is the second order zonal coefficient, R
e
the
earth radius, µ the standard gravitational parameter
and r =
p
w
is the radius. Without perturbations, all
equinoctial elements are constant except the true lon-
gitude L which is close to linear (see Equation (2)).
2.3 Satellite Tracking by Ground Radar
In order to detect new objects, a ground radar scans
periodically the space. Then, measure association is
performed: every produced measure is compared to
the estimates of known objects in order to decide if
the measure was generated by a known object or by
a previously unknown one. If the measure belongs to
a known object, its estimate is updated and improved
by the information of the measure. If the measure
belongs to a new object, an estimate of this object is
computed and the object will be considered as known
for the future measures performed by the radar.
In this paper, the considered radar produces mea-
sures of the position of detected objects in the ECI
frame. This simplification is justified by the fact
that measures performed in a Radar frame, such as
Azimuth-Elevation-Range, can easily be converted in
a Cartesian frame using classical formulas (Vallado,
2001). The measures provided by the radar come with
a covariance matrix, also in ECI, representing the un-
certainties of the associated measures. Mathemati-
cally, each of these measures can be modelled by a
random variable X
mes
∈ R
3
normally distributed with
mean x
true
∈ R
3
(the true position of the satellite in
the ECI frame) and covariance P ∈ R
3×3
(depending
on the radar performances and the satellite position).
3 SET-BASED PROPAGATION
One of the most used tool to handle uncertainty en-
closed by sets is the interval arithmetic.
3.1 Interval Arithmetic
By using interval arithmetic (Jaulin et al., 2001), un-
certainties can be represented with a deterministic set-
based approach. An interval [a] = [a, a] defines the
set of a ∈ R such that a ≤ a ≤ a. IR denotes the
set of all intervals over R. The extension of mathe-
matical operators for interval arithmetic is defined as
[a]∗[b] = {a ∗b|a ∈[a], b ∈[b]} with ∗a binary oper-
ator on real numbers. Such operations can sometimes
be expressed via the bounds of the intervals, but not
always. For instance, the sum of two intervals [a] and
[b] gives

a + b, a + b

. Elementary functions can also
be extended to interval arithmetic. A tube denotes the
image of a function f : t ∈ R 7→ [ f (t)] ∈ IR. Further-
more, the Cartesian product of n intervals is called an
interval vector or a box [a] ∈ IR
n
.
An interval contractor associated to the set A ⊂
R
n
is an operator C : [a] ∈ IR
n
7→C([a]) ∈ IR
n
satis-
fying the contractance property (C([a]) ⊆ [a]) and the
completeness property (C([a]) ∩A = [a] ∩A). Associ-
ated to a constraint, a contractor allows to reduce the
size of an interval while preserving the subset of real
numbers verifying the constraint.
3.2 Validated Simulation
Numerical integration consists in approximating a so-
lution of an ordinary differential equation (ODE),
˙y = f (y), using an integration scheme (such as Euler
or Runge-Kutta). Interval arithmetic can be applied
to numerical integration in order to provide validated
numerical integration methods, also named reachabil-
ity analysis or guaranteed simulation.
Such validated integration methods use a guaran-
teed integration scheme, providing a discretization of
time, t
0
⩽ ··· ⩽ t
end
, and a computation of enclosures
of the set of states of the system y
0
, . . . , y
end
.
A guaranteed integration scheme consists in
an integration method Φ( f , y
j
, t
j
, h), approximat-
ing the exact solution y(t
j+1
;y
j
), i.e., y(t
j+1
;y
j
) ≈
Φ( f , y
j
, t
j
, h), where y
j
is an initial value, t
j
the ini-
tial time, and h the step-size (t
j+1
= t
j
+ h), and a
truncation error function LTE
Φ
( f , y
j
, t
j
, h), such that
y(t
j+1
;y
j
) = Φ( f , y
j
, t
j
, h) + LTE
Φ
( f , y
j
, t
j
, h). Tak-
ing into account the inaccuracy of the used integration
scheme, in contrast with usual integration methods, a
validated numerical integration method is a two-step
method which, firstly, computes an enclosure [ ˜y
j
] of
the solution over the time interval [t
j
, t
j+1
] to bound
LTE
Φ
( f , y
j
, t
j
, h), then, computes a tight enclosure of
the solution for the final time instant t
j+1
. Multiple
methods exist to perform these steps, such as Taylor
series and Runge-Kutta, see (Nedialkov et al., 1999;
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
662

Alexandre dit Sandretto and Chapoutot, 2016) and the
references therein for more details. The step-size can
be fixed (h is a constant) or adaptive (h is computed
at each step). Finally, functions R and
˜
R are provided
by validated numerical integration methods
R :

R → IR
n
t 7→ [y]
,
˜
R :

IR → IR
n
[t, t] 7→ [˜y]
(4)
with for a given t
i
, R(t
i
) = {y(t
i
;x
0
) : ∀y
0
∈[y
0
]}⊆[y],
and
˜
R([t, t]) = {y(t; y
0
) : ∀y
0
∈ [y
0
] ∧∀t ∈ [t, t]} ⊆ [ ˜y].
3.3 Confidence Contractor
Intervals and probabilities are strongly connected. In-
deed, the universe of a distribution is an interval, and
a given confidence level is associated to a confidence
interval. A confidence interval is a set S for which the
probability of the given random variable to be in this
set is equal to the given probability P. In other words,
let ˆx be a single observed sample of the quantity X. In
statistics, an observed data allows to compute a confi-
dence interval (Neyman, 1937), that is to say an inter-
val which may contain the actual value, with respect
to a given confidence level. For example, considering
a confidence level CL = 95%, one can define the con-
fidence interval C
95%
. This interval can be obtained
by observation (statistical approach) or with the help
of a known distribution (probabilistic approach). A
new measure ˆx coming from the (same) experiment
will be in the associated confidence interval such that
ˆx ∈C
95%
95% of the time.
We use in our algorithm the confidence contractor
proposed in (Alexandre Dit Sandretto, 2019):
Cbc([x]|f
X
, cc) : IR → IR
[x] 7→ [x] ∩[y]
(5)
with [y] defined such that Pr(x ∈ [y]) =
R
[y]
f
X
(x)dx =
cc ([y] is the confidence interval), Pr(x ∈ [y]) stands
for “probability that x is in [y]”, with x following the
distribution f
X
, and cc being the confidence coeffi-
cient (0 ≤ cc ≤ 1). For example, cc = 0.68 for a 68%
confidence level. In the domain of ODEs, confidence
contractor has been used to produce potential clouds
(Fuchs and Neumaier, 2009) in (Alexandre Dit San-
dretto, 2019).
4 ASSOCIATION, ESTIMATION
Based on interval arithmetic and validated simulation,
our approach uses measures of the satellite’s position
to improve the estimation of its trajectory and rule out
incompatible measures and tracks.
4.1 Definition of Measures
To measure a quantity x
true
∈R
n
, instruments provide
a measure with its associated uncertainty, in either a
deterministic (Equation (6)) or a probabilistic (Equa-
tion (7)) approach. In the deterministic case, the mea-
sure x
mes
∈ R
n
and uncertainty ∆x
mes
∈ R
n
+
generate
an interval vector that necessarily includes x
true
:
x
true
∈ [x
mes
±∆x
mes
]. (6)
In the probabilistic approach, the measure is repre-
sented as a random variable X
mes
normally distributed
with mean x
true
and covariance P ∈ R
n×n
:
X
mes
∼ N (x
true
, P). (7)
In the probabilistic approach, the uncertainty is car-
ried by the covariance matrix P. Usually, the prob-
abilistic approach is preferred to represent uncertain
measures. However, the downside of this method is
its lack of accuracy to describe rare events and sin-
gularities that have a low probability to be drawn. In
contrast, the deterministic, or set-based approach al-
lows to take the entire measure interval into consid-
eration, which is significantly more robust with worst
case scenarios. For this reason, we choose to consider
measures as real valued intervals.
For our application, the ground radar provides
probabilistic uncertainties. Thus, a conversion from
the uncertainty covariance P to a deterministic uncer-
tainty ∆x
m
is required and explained in section 5.
4.2 Representation with Contractors
Two cases are considered: an ideal case where com-
plete observation of the system is provided and a more
realistic case where this observation is only partial.
Complete Observation of the System
Consider a dynamical system of state x
true
(t) ∈ R
n
at
time t ∈ R, following the dynamics ˙x
true
= f (x
true
)
with f : R
n
→R
n
. The state is measured with a deter-
ministic approach, as in Equation (6).
On one hand, measures of the dynamical system
provide an enclosure of the state. At time t, the mea-
sure x
m
t
∈R
n
and its associated uncertainty ∆x
m
t
∈R
n
+
are obtained, and give the following inclusion:
x
true
(t) ∈ [x
m
t
±∆x
m
t
]. (8)
On the other hand, validated numerical integration of
the dynamical system provides an interval estimate
encompassing the state of the system. At time t, the
estimate [ ˆx
t
] of the system is given by:
[ ˆx
t
] := R(t) ∈ IR
n
, (9)
Validated Uncertainty Propagation for Estimation and Measure Association, Application to Satellite Tracking
663

where R(·) is defined in Equation (4), and R(t
0
) = x
0
with x
0
the initial value of the system, at time t
0
. The
true state at time t is also included in [ ˆx
t
]:
x
true
(t) ∈ [ ˆx
t
]. (10)
Using Equations (8) and (10), one concludes that the
state of the dynamical system is included in the inter-
section of the measure and the estimate:
x
true
(t) ∈ [x
m
t
±∆x
m
t
] ∩[ ˆx
t
]. (11)
Further, if the exact time t ∈ R of the measure is
not known precisely, it can be approximated by a time
interval. Assume that the exact time t is bounded by
t
≤t ≤ t, then instead of using R(·) to provide the es-
timate as in Equation (9), the function
˜
R(·) generates
an estimate [ ˆx
t
] :=
˜
R(t) ∈ IR
n
, that takes account of
the uncertainty on the time of the measure.
Partial Observation of the System
In some cases, the intervals measured from the dy-
namical system are only partial, or in a different ref-
erence frame. Let us consider that the sensor can only
make observations y
m
t
∈ R
p
of the system x
true
(t) ∈
R
n
, with p ≤ n. A function g : R
n
→ R
p
gives the
constraint y
m
t
= g(x
true
(t)), or equivalently
g(x
true
(t)) −y
m
t
= 0, (12)
which represents the compatibility of the system x
true
with the measure y
m
t
at time t. If a measure y
m
t
does
not verify this constraint, then it cannot come from the
system x
true
. Further, if the measure y
m
t
comes from
the system x
true
, the constraint in Equation (12) pro-
vides more information on the state of system x
true
(t).
Given a partial measure interval [y
m
t
] ∈IR
p
and an
estimated interval [ˆx
t
] ∈ IR
n
of a system at time t, the
interval counterpart of Equation (12) is as follows:

g([ ˆx
t
]) ⊂ [y
m
t
]
g
−1
([y
m
t
]) ⊂ [ ˆx
t
]
(13)
The smallest subset of [ ˆx
t
] that respects Equation (13)
is computed with a forward-backward contractor
(Benhamou F. and J-F., 1999).
As a result, interval contractors allow to use the
partial information provided by the measures when
the measures are only partial and cannot be used to di-
rectly intersect the estimated interval (Equation (11)).
4.3 Association of Measures
Comparing the estimate of a dynamical system with
new measures helps reducing the uncertainty of the
estimate and contributes to the association issue.
Our algorithm (see Algorithm 1) ensures that a
measure is compatible with the simulated dynami-
cal system to associate them. Since the simulation
is validated, the estimated interval vector necessarily
includes every possible states of the satellite. So, if
a measure interval does not intersect with the simu-
lation interval associated with that instant, this nec-
essarily means that the measurement is incompatible
with that track. Nonetheless, it is possible that certain
measures give a non-empty intersection with the sim-
ulation estimate, but do not really correspond to the
same system. Hence the significance of introducing
confidence levels in the measurements.
Data: Estimate [ ˆx
t
], Measure [x
m
t
±∆x
m
t
]
Result: Association decision at time t
if [ ˆx
t
] ∩[x
m
t
±∆x
m
t
] ̸=
/
0 then
Association because of compatibility;
Improvement of the estimate:
[ ˆx
t
] ← [ ˆx
t
] ∩[x
m
t
±∆x
m
t
];
else
Incompatibility, no association ;
end
Algorithm 1: Measure Association.
Once the algorithm ruled out incompatible mea-
sures, it uses the additional information of the mea-
sure to increase the precision on the state of the dy-
namical system. The system is included in both the
measure interval and the estimated interval. Conse-
quently, it is included in their intersection. This gives
an upgrade of the estimate by reducing the uncertainty
on the state of the system after the measure.
5 RESULTS
The presented method has been implemented in the
DynIbex framework (Alexandre dit Sandretto and
Chapoutot, 2016). It provides interval arithmetic
tools such as contractors and offers a validated nu-
merical integration procedure based on Runge-Kutta
methods with adaptive step-size. The presented sim-
ulations were computed using the implicit midpoint
method. The local truncation error was set to 10
−7
.
This choice of parameters granted a good trade-off
between computation time and precision.
The objects were simulated in mEOE. The semi-
parameter p was normalised by the Earth radius and
the time unit for the integration was hours instead of
seconds to facilitate the computations. The measures
provided by the ground radar were on the position of
the satellite, in the Cartesian reference ECI.
The initial estimate of each simulated satellite is
supposed to be normally distributed with a mean vec-
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
664

Figure 1: Simulation of a trajectory of a satellite for 6 hours,
represented by the boxes of uncertainties on its position.
tor x
0
and a covariance matrix P. x
0
depends on the
scenario and P is chosen diagonal such that P
11
=
594.817, P
66
= 1.258 ·10
−11
and P
ii
= 5.948 ·10
−12
for i = {2, . . . , 5}. This covariance matrix is repre-
sentative of the uncertainties on the state of a satellite
after a few measures. It is satisfactory to demonstrate
the contribution of our algorithm but in operating en-
vironment, distinct covariance matrices for the satel-
lites of interest would be considered, depending on
the measures used to compute their estimate, and a
given time at which the matrices are evaluated.
Then, a reference box is computed using the
confidence contractor of Equation (5) with a 99%
confidence level: ∆x
ref
= [−x
ref
, +x
ref
] where x
ref
=

100, 10
−6
, 10
−6
, 10
−6
, 10
−6
, 1.45408 ·10
−5

. This
represents an uncertainty of ±100 meters on the semi-
parameter p, and approximately ±100 meters on the
arc of circle represented by the true longitude L when
the altitude of the satellite is 500 kilometers.
The potential clouds associated to the 5% and 95%
confidence levels are given by ∆x
5%
= 0.21∆x
ref
and
∆x
95%
= 0.645∆x
ref
. Then, with x
0
an initial state of a
satellite, the initial box associated to the n% potential
cloud (n = 5, 95) is given by x
0
+ ∆x
n%
.
5.1 Sc. 1: Propagation of Uncertainties
The first scenario consists in the propagation of the
5% and 95% potential clouds of a satellite for 6 hours,
using the function R defined in Equation (4).
Figure 2 shows the results on the components p
and L. For an easier comparison between the two po-
tential clouds, the intervals are translated using the
center of the 95% interval. After 6 hours, the 95%
potential cloud associated to the semi-parameter p is
about 600 meters wide, three times greater than the
5% potential cloud. Along the propagation, the or-
der of magnitude stays the same for the two poten-
tial clouds associated to p. In contrast, the widths of
Figure 2: Propagation of the 5% (blue) and 95% (red) po-
tential clouds on the p and L components.
the potential clouds associated to L increase faster and
after 6 hours, the circular arcs corresponding to such
intervals are about 7 kilometers for the 5% potential
cloud and 27 km for the 95% potential cloud. The
specificity of the component L is discussed in 5.2. Fi-
nally, it can be observed that, as expected, the 5% in-
tervals are always included in the 95% intervals.
5.2 Sc. 2: Effect of Measures on
Estimation Uncertainties
In this scenario, a ground radar provided an estimate
of the state of the satellite at an initial time, and prop-
agated its estimate for 12 hours with the function R
(see Equation (4)) before being able to measure its
position again. When the satellite gets in the scope of
the radar again, a series of 10 measures of its position
over 5 minutes are taken. Using our algorithm, these
measures are added to the estimation process.
0 5 10 15
Time (h)
0
0.5
1
1.5
2
Width of p (m)
10
4
p
0 5 10 15
Time (h)
0
1
2
3
4
Width of f (adim)
10
-3
f
0 5 10 15
Time (h)
0
1
2
3
4
Width of g (adim)
10
-3
g
0 5 10 15
Time (h)
0
2
4
6
8
Width of h (adim)
10
-4
h
0 5 10 15
Time (h)
0
2
4
6
Width of k (adim)
10
-4
k
0 5 10 15
Time (h)
0
0.05
0.1
0.15
Width of L (rad)
L
Propagation of a satellite during 15 hours
with estimation improvement using 10 measures at 12 hours
Figure 3: Propagation of the uncertainties, with (blue) and
without (black) addition of 10 measures after 12 hours (red).
Validated Uncertainty Propagation for Estimation and Measure Association, Application to Satellite Tracking
665

Figure 3 shows that L is significantly impacted by
the series of measures, in comparison to other com-
ponents of the state of the satellite. Indeed, the un-
certainty over L drops with the addition of measures.
The precision of the other components also improves
after the addition of measures, compared to the sce-
nario where no measure was added. This results from
the precision gain of L spreading to the other state’s
components through the dynamics equations.
The impact of the measures on the component
L is due to L encoding the position of the satellite
along the orbit described by the five other compo-
nents. Therefore, this component is strongly linked to
the position of the satellite, hence the steep decrease
of its uncertainties when a position measure is added.
Further, the true longitude L is the component which
is the most subject to variations over time, since it is
the only non-constant one when the dynamics is un-
perturbed (see Equation (2)). Therefore, it makes the
propagated uncertainties grow faster for the true lon-
gitude L. By focusing on the impact of each measure
of the series on L, it gives Figure 4.
11.9 11.95 12 12.05 12.1
Time (h)
-0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
Width of L (rad)
Influence of measures on the true longitude L
Simulation with intermediate measures
Simulation without intermediate measures
Instants of measures
Figure 4: Focus on the effect of each measure on L.
Figure 4 helps raising two points. First, not every
measure helps improving the estimate of the satellite’s
position. Here, the first measure of the series does not
reduce the estimate’s uncertainty. Mathematically, it
means that the estimation interval was included in the
measure interval. The radar is not precise enough to
bring new information. Secondly, most of the preci-
sion given by the series of measures comes from the
first measures. This could allow computing the num-
ber of measures that would be necessary to achieve a
certain precision on the position of the tracked satel-
lite, depending on the accuracy of the ground radar.
5.3 Sc. 3: Potential Clouds for the
Measure Association Issue
The 5% and 95% potential clouds for two satellites
are propagated using validated simulation and the
˜
R
function (see Equation (4)) with the initial states:













p
1
= 6890939, p
2
= 6961989,
f
1
= 0.0014125, f
2
= 0.0014182,
g
1
= −0.0014160, g
2
= −0.0014102,
h
1
= 0.57763, h
2
= −0.57763,
k
1
= 0.0095098, k
2
= −0.0095068,
L
1
= −0.32425, L
2
= −0.31903.
(14)
At last, the trajectories of the two satellites can be
compared.
6.85 6.9 6.95
x (m)
10
6
-8
-6
-4
-2
0
2
4
6
y (m)
10
5
6.86 6.88 6.9 6.92 6.94
x (m)
10
6
2.9
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
z (m)
10
5
-5 0 5
y (m)
10
5
2.8
2.9
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
z (m)
10
5
Intersection with the 95% confidence level
Figure 5: Intersection of position boxes in the ECI frame for
the propagation of the 95% potential cloud of two satellites.
6.85 6.9 6.95
x (m)
10
6
-8
-6
-4
-2
0
2
4
6
y (m)
10
5
6.86 6.88 6.9 6.92 6.94
x (m)
10
6
2.9
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
z (m)
10
5
-6 -4 -2 0 2 4
y (m)
10
5
2.8
2.9
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
z (m)
10
5
No intersection with the 5% confidence level
Figure 6: No intersection of position boxes in the ECI frame
for the propagation of the 5% potential cloud of two satel-
lites at the same time that Figure 5.
Figures 5 and 6 show a sample of the trajectories
where the propagated boxes intersect for a 95% confi-
dence level but not for a 5% confidence level. In that
respect, the sensor performances can be evaluated in
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
666

terms of required confidence level and decisions can
be made depending on this analysis. For example, in
that case, if a sensor can only ensure that a collision
will never happen with a 5% confidence level, one can
decide to use complementary sources of measures in
order to improve the estimation of the satellites, or to
maneuver one of the satellites in order to increase the
confidence level associated to the collision.
6 CONCLUSION AND FUTURE
WORKS
In this paper, we proposed a validated algorithm to
estimate the trajectory of dynamical systems taking
account of uncertainties on the initial state. Our ap-
proach uses measures and their uncertainties to quan-
tify likelihood to belong to a certain track. Unlike
Monte Carlo methods, our approach provides math-
ematically guaranteed results that do not ignore low-
probability cases, while taking advantage of probabil-
ity distributions to supplement the information pro-
vided by measures. Further, Monte Carlo methods
require a substantial number of simulations to be reli-
able, compared to set-based validated methods which
only need one simulation to have as much information
with more guarantees.
As future works, implementing the retro-
propagation feature of DynIbex to our algorithm
would allow using the information given by a
measure to reduce the uncertainty of the track be-
forehand. Moreover, studying the impact of a series
of measures depending on the number of measures,
the state of the system, the precision of the estimate
and the precision of the sensor could help computing
the optimal times for taking measures, to gain a
maximum of information on the system.
REFERENCES
Alexandre Dit Sandretto, J. (2019). Confidence-based con-
tractor, propagation and potential cloud for differen-
tial equations. Summer Workshop on Interval Methods
2019.
Alexandre dit Sandretto, J. and Chapoutot, A. (2016). Vali-
dated Explicit and Implicit Runge-Kutta Methods. Re-
liable Computing, 22.
Bar-Shalom, Y. and Tse, E. (1975). Tracking in a cluttered
environment with probabilistic data association. Au-
tomatica, 11(5):451–460.
Benhamou F., Goualard F., G. L. and J-F., P. (1999). Re-
vising hull and box consistency. In Proceedings of
the 16th International Conference on Logic program-
ming, page 230–244. The MIT Press.
Blackman, S. and Popoli, R. (1999). Design and Analysis of
Modern Tracking Systems. Artech House radar library.
Artech House.
Collins, J. and Uhlmann, J. (1992). Efficient gating in
data association with multivariate gaussian distributed
states. IEEE Transactions on Aerospace and Elec-
tronic Systems, 28(3):909–916.
Fortmann, T., Bar-Shalom, Y., and Scheffe, M. (1983).
Sonar tracking of multiple targets using joint proba-
bilistic data association. IEEE journal of Oceanic En-
gineering, 8(3):173–184.
Fuchs, M. and Neumaier, A. (2009). Autonomous robust
design optimisation with potential clouds. Int. J. Reli-
ability and Safety, 3:23–34.
Gill, P. E. and Murray, W. (1978). Algorithms for the so-
lution of the nonlinear least-squares problem. SIAM
Journal on Numerical Analysis, 15(5):977–992.
Horwood, J. T. and Poore, A. B. (2014). Gauss von mises
distribution for improved uncertainty realism in space
situational awareness. SIAM/ASA Journal on uncer-
tainty Quantification, 2(1):276–304.
Jaulin, L., Kieffer, M., Didrit, O., and Walter, E. (2001).
Applied Interval Analysis. Springer London.
Julier, S. and Uhlmann, J. (2004). Unscented filtering
and nonlinear estimation. Proceedings of the IEEE,
92(3):401–422.
Mills-Tettey, G. A., Stentz, A., and Dias, M. B. (2007).
The dynamic hungarian algorithm for the assignment
problem with changing costs. Robotics Institute, Pitts-
burgh, PA, Tech. Rep. CMU-RI-TR-07-27.
Nedialkov, N. S., Jackson, K., and Corliss, G. (1999). Val-
idated solutions of initial value problems for ordi-
nary differential equations. Appl. Math. and Comp.,
105(1):21 – 68.
Neyman, J. (1937). Outline of a theory of statistical estima-
tion based on the classical theory of probability. Phil.
Trans. R. Soc. Lond. A, 236(767):333–380.
Reid, D. (1979). An algorithm for tracking multiple targets.
IEEE Transactions on Automatic Control, 24(6):843–
854.
Vallado, D. A. (2001). Fundamentals of astrodynamics and
applications, volume 12. Springer Science & Busi-
ness Media.
Verhagen, S. and Teunissen, P. J. (2017). Least-squares
estimation and kalman filtering. In Teunissen, P. J.
and Montenbruck, O., editors, Springer Handbook of
Global Navigation Satellite Systems. Springer Inter-
national Publishing.
Walker, M. J., Ireland, B., and Owens, J. (1985). A set
modified equinoctial orbit elements (errata 1986). Ce-
lestial mechanics, 36(4):409–419.
Validated Uncertainty Propagation for Estimation and Measure Association, Application to Satellite Tracking
667