Optimal Control-Based Track-to-Track Correlation with Optical Measurements

Alessia De Riz
Politecnico di Milano, Department of Aerospace Science and Technology (DAER), Via G.
La Masa 34, 20156, Milan, Italy
Riccardo Cipollone
Politecnico di Milano, Department of Aerospace Science and Technology (DAER), Via G.
La Masa 34, 20156, Milan, Italy

Pierluigi Di Lizia
Politecnico di Milano, Department of Aerospace Science and Technology (DAER), Via G.
La Masa 34, 20156, Milan, Italy

ABSTRACT

Continuous efforts are underway to pioneer innovative techniques and solutions to improve the surveillance
and monitoring of Resident Space Objects (RSOs). However, with the growing interest in cislunar space, the
scientific community is actively mobilizing to extend Space Situational Awareness (SSA) and Space Surveil-
lance and Tracking (SST) activities to this novel domain.
In this context, object catalog maintenance and update stand as essential. In this work, emphasis is placed
on track-to-track correlation, a pivotal component for building and preserving a reliable catalog. While
well-established techniques exist to tackle this issue in the two-body dynamics framework, many rely on its
specific features, making their possible extension to the cislunar domain unfeasible.
A control distance metric-based formulation, explored in the literature as a maneuver detection and cor-
relation metric, emerges as the most suitable method as it is free from strong assumptions depending on
a specific dynamics model. In this regard, this work aims to explore a novel procedure that leverages an
energy-optimal control metric, tailoring it to the application via its transversality conditions.

1. INTRODUCTION

In recent years, space has acquired significant attention from both government agencies and private compa-
nies, driven not only by its vast potential and resources but also by its increasing management complexity.
The near-Earth environment is densely populated with space assets providing a wide array of services,
including telecommunications, navigation, and climate monitoring. These vital functions are increasingly
endangered by the risk of collisions, involving either active satellites occupying congested orbital slots or
space debris fragments. Consequently, there is an ever-increasing need to mitigate the risk of exacerbating
the Kessler syndrome [1] through the development and enhancement of SST tools together with the estab-
lishment of robust space debris mitigation guidelines [2]. Recently, the scientific community has turned its
focus to the cislunar environment, a relatively unexplored region of space between the Earth and the Moon.
This new area of interest is motivated by both the desire for further space exploration and the potential
exploitation of in situ resources, leading to an expected increase in satellite populations in this region. The
near-Earth precedent serves as a cautionary tale, underscoring the importance of proactive measures to
prevent similar issues. In this context, object catalog maintenance and update become crucial. Tracking
both active and inactive objects allows for monitoring population trends and generating vital statistics, as
demonstrated by efforts like ESA’s Annual Space Environment Report [3].
As part of the cataloging pipeline, track-to-track correlation or association plays a critical role when a new
Resident Space Object (RSO) is detected and needs to be integrated into the catalog. At the core of this
process lies the concept of track, which denotes a series of consecutive observations of the same target taken

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
by a single sensor within a restrained time-frame, usually insufficient to estimate a reliable orbit. When it
comes to measurement processing, an acquired track is typically labeled as an Uncorrelated Track (UCT)
until it can be matched with a cataloged object that is confirmed to have generated it. While well-established
techniques exist to tackle this issue in the two-body dynamics framework, many of them rely on its specific
features, making their possible extension to the cislunar domain unfeasible.
Among these, one can find techniques leveraging the conservation of the Keplerian motion characteristic
quantities, as first introduced by Taff and Hall [4] and further developed by Gronchi [5] within the two-body
integral method. Another key concept, first introduced by Milani [6] to address the problem of asteroid
identification based on very short arc observations, and extensively employed in the SSA community, is the
Admissible Region (AR). Since the range and range rate are left undetermined in an optical track-to-track
correlation problem, the idea behind the AR is to derive conditions on those parameters to bind the solution
space. These conditions are typically retrieved making hypotheses on the target’s dynamics. This tool has
been extended and adapted through various approaches, for instance employing different ways to sample the
AR, in Tommei [7] a Delaunay triangulation approach is used, in Pirovano [8] a Differential Algebra inter-
pretation is implemented while in De Mars [9] a uniform distribution approximated as a Gaussian Mixture
Model is employed.
On the other hand, a control distance metric-based formulation, explored in the literature as a maneuver
detection and correlation metric, emerges as the most suitable method as it is free from strong a priori
approximations. An early application of this approach can be found in Holzinger [10], devising a technique
aimed at determining the minimum required control effort to fit a given observation, using Distance Control
Metrics instead of the more common Mahalanobis distance. This work paved the way for further devel-
opments and applications, such as in Pastor [11] where impulsive burns are considered for both maneuver
detection and correlation purposes. To narrow down the observable space in which searching for the min-
imum, this formulation is commonly paired with an admissible region approach, as seen in the works of
Serra [12] and Siminski [13]. However, defining the constraints to bound such a region presents challenges
in scenarios that deviate from the classic two-body dynamics.
In this regard, this work aims to explore a novel procedure that leverages an energy optimal control metric,
tailoring it to the application via its typical transversality conditions, used to replace the definition of a
properly constrained admissible region. This technique offers a dual advantage: on the one hand, it enables
the association of tracks linked by purely ballistic dynamics, on the other hand, it concurrently allows for
robustness to maneuvers. This capability significantly aids catalog maintenance, effectively mitigating the
occurrence of duplicated objects. Furthermore, such a technique implies no assumption on the target orbital
dynamics, enabling the possibility of its extension to a cislunar framework.

2. METHODOLOGY

This section outlines the core of the methodology and the primary tools that constitute it, starting with a
brief description of the minimum-energy Optimal Control Problem (OCP) formulation relevant to the case
at hand.

2.1 Minimum-energy Optimal Control Problem framework

The developed technique aims at solving a correlation problem between two optical UCTs. Thus, the
only available information consists of the target’s angular coordinates contained within the UCT itself and
provided in terms of Right Ascension (α) and Declination (δ) pairs which also coincide with the boundary
conditions (BCs) of the problem. The choice of optical data over radar is driven by the fact that, for future
extensions into the cislunar domain, optical sensors currently offer the most suitable option for space-based
SST applications. To estimate the state of the target within an Earth-Centered Inertial (ECI) reference
frame, the following relationships are generally employed:

r = Rs + ρs,
v = Vs + ρṡ + ρ̇s, (1)
s = (cos α cos δ, sin α cos δ, sin δ),

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
where r is the target’s position and v its velocity, while Rs and Vs correspond to the ground station position
and velocity and s is the line-of-sight unit vector. However, angular information alone is not sufficient to solve
the problem. To address this, a track compression technique is employed to extract additional information
about the angular velocities (α̇, δ̇). The procedure selected for the case at hand involves a second-order
regression, performed at the midpoint of the UCT, which allows all the information contained in the track
to be condensed into a single point. The result is a vector that includes both the angular coordinates and
their respective time derivatives, denoted as a = [α, δ, α̇, δ̇], which is commonly referred to in the literature
as an attributable [14].
To link two condensed UCTs, an energy-optimal control problem is selected and its cost is expressed as a
thrust energy:
∫ tf
1 T
J= u (τ ) u(τ ) dτ ,
2 t0

where u(t) is the acceleration profile while t0 and tf define the initial and final epoch respectively. This
assumption of optimality is considered reliable for two reasons: the inherent tendency to save fuel during
any controlled operation and the primary goal of correlation, which seeks to identify associations through
trajectories and underlying expenses that are usually compatible with ballistic motion.
The resulting OCP formulation is the following:


 ẋ = f (t, x(t), u(t))
∫ tf 

1 T h(x(t0 )) = a0
min J = u (τ ) u(τ ) dτ s.t. (2)
u(τ ) 2 
 h(x(tf )) = af
t0 

t0 , tf given ,

where h(x(t)) is the transformation from a state x at time t to the correspondent attributable vector,
x(t0 ) and x(tf ) stand respectively for the state at initial time and final time, a0 and af are the considered
attributables while ẋ = f (t, x(t), u(t)) represents the problem’s dynamics:
[ ] [ ]
ṙ v
ẋ = f (t, x(t), u(t)) = = . (3)
v̇ − r3 r + u
µ

where µ is the Earth’s gravitational parameter. The OCP is solved through an indirect approach, augmenting
the state with an adjoint one λ, also called co-state vector, resulting in a Two-Point Boundary Value Problem
(TPBVP). By leveraging the problem’s Hamiltonian, defined as:
( )
1 T µ
H = l(t, x(t), u(t)) + λ(t)T f (t, x(t), u(t)) = u u + λr · v + λv · − 3 r + u , (4)
2 r
where λr and λv respectively refer to the first and second 3 components of the adjoint state, one can retrieve
the relation u(t) = −λv (t), from the null Hamiltonian derivative with respect to the control (necessary
condition for optimality). As for its partial derivatives with respect to the co-state Hλ and state Hx , they
are employed to define the TPBVP equations:



 ẋ = Hλ

λ̇ = −H
x
(5)

 h(x(t0 )) = a0


h(x(tf )) = af .

The initial value for λ(t) has to be retrieved to solve the problem. Then by integrating over time, the
target continuous control profile is achieved along with the corresponding energy expenditure. When dealing

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
with the BCs of an optical track correlation problem, two approaches can be used to solve the TPBVP
and determine the optimal expense. The first enforces transversality conditions, leveraging the available
observables to obtain an optimal solution. The second approach instead involves a nested optimization
process, searching for the minimum optimal expense within a space defined by the missing observables. In
this instance, the first approach is used to fully describe the state at the initial time.

2.2 Transversality conditions approach

Referring to Eq. (1), it can be seen that two variables, more precisely range and range rate, remain unknown,
making it impossible to define the target’s state. Many of the methodologies in the literature use other
techniques, alongside an effort metric formulation, such as the AR to search for these missing unknowns.
For the case at hand, transversality conditions, an inherent connection between co-state and BCs stemming
from the OCP indirect formulation, are instead used to compensate for this problem. Since neither x(t0 ) nor
x(tf ) are fixed but only partially constrained, the Pontragyn Maximum Principle (PMP) provides additional
information, helping in the definition of additional equations.
Transversality conditions are in fact necessary conditions of optimality given constraints and costs on the
state at both initial and final time [10]. Thus, λ0 and λf values must satisfy both the hard constraints
represented by BCs and soft ones, denoted by initial or terminal costs values, if present. Their overall
expression can thus be formulated as:

∂V0 T ∂g0 T
λ0 = − − ν0
∂x0 ∂x0
(6)
∂Vf T ∂gf T
λf = + νf ,
∂xf ∂xf

where V0 represents the costs at initial time V (x0 , t0 ) while Vf corresponds to the one at final time V (xf , tf ).
The variables ν are Lagrange multipliers introduced in the formulation to enforce the BCs, expressed as initial
and final constraints g0 and gf , respectively standing for g (x0 , t0 ), g (xf , tf ).
Let us now analyze the form of Eq. (6) for the case under analysis. The term related to the cost is not
present while the BCs are formulated using the information on the initial and final attributables combined
with the relationship linking Cartesian to spherical coordinates.
The constraints g(x0 , t0 ) and g(xf , tf ) are therefore described as vectors, with the following structure:

g0 = ∆ a0 (x0 ) = a0 (x0 ) − a0
(7)
gf = ∆ af (xf ) = af (xf ) − af ,
where the first term embeds the coordinates change while the second, overlined, the measurements set of the
two attributables. Each component of the BCs general expression, denoted as g (x), is detailed below:
(y)
g1 = α − α = atan −α
( x )
z
g2 = δ − δ = asin √ −δ
x2 + y 2 + z 2
(8)
xvy − yvx
g3 = α̇ − α̇ = 2 − α̇
y + x2
√ (( ) )
(y2 + x2 ) y2 + x2 vz + (−yvy − xvx ) vz
g4 = δ̇ − δ̇ = − δ̇ .
(y2 + x2 ) z2 + y4 + 2x2 y2 + x4

In this way, exploiting Eq. (6), an analytical expression for λ0 and λf can be retrieved, which leaves as
unknowns the values of the related Lagrange multipliers ν and the state at initial time x0 .
The problem is then iteratively solved via single shooting, once a first guess for both initial state x0 and
Lagrange multipliers [ν0 ,νf ] are defined. To retrieve a suitable initial guess of the state, an Initial Orbit

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
Determination (IOD) process between the two selected attributables is performed.
The non-linear system built to retrieve the OCP solution is composed of 14 equations in 14 unknowns:

 λ (tf ) − λf = 0
g =0 (9)
 0
gf = 0 ,

where a further condition on the final co-state is added to the already mentioned BCs. Both the final state
x(tf ) and co-state λ(tf ) are the outcome of the OCP dynamics propagation flows.

x(tf ) = ϕx (tf ; x0 , λ0 , t0 )
λ(tf ) = ϕλ (tf ; x0 , λ0 , t0 )

2.3 Uncertainties and correlation

When values coming from experimental measurements are considered as variables inside a function, their
uncertainty, mainly due to measurement limitations (i.e. instrument precision), propagates across the func-
tion, affecting its output. In order to take this phenomenon into account, the input quantities have to
be considered as distributions that can be approximately described by a finite number of their statistical
moments. This section deals with how the uncertainty characterizing the angular measurements involved in
the problem translates into the final results, generating distributions of solutions.
The first step is to derive the covariance relative to the whole attributable a since α and δ accuracy can be
directly retrieved from the instrumentation. They are in fact known and taken as standard deviation σ of
the measurements, thus describing the corresponding distribution, assumed Gaussian, by mean value and
covariance. To propagate them to the angular velocities, the process leverages a Monte Carlo simulation of
the previously introduced track compression. Firstly, the available normal distribution statistics of the angles
composing the tracks under analysis are used to sample a population. Then, by applying the compression
on each sample, a set of (α̇, δ̇) is obtained and employed to extract covariance. Having all the data required,
the attributable covariance matrix Σa is built as a 4x4 diagonal.
The underlying idea related to the uncertainty integration is that of solving the optimal control problem
introduced in the previous sections to find what is considered as the nominal expense ∆Vnom . Then, the
obtained solution in terms of state x0,ref is propagated using a ballistic dynamical model, free from additional
accelerations, from the initial t0 to the final tf time. Once the two states x0,ref and xf,ref are retrieved, the
corresponding attributables are computed by projecting them onto the measurement space. This information
is used to solve again an OCP on the ballistic trajectory, taking uncertainty into account, to compute an
approximate background expense related to the only problem uncertain BCs, N0 (∆Vb , Σ∆Vb ).
If the nominal solution ∆Vnom can be explained by the background distribution N0 , then the hypothesis of
a ballistic motion linking the two BCs is confirmed.
To understand this, one possible approach is considering the nominal expense ∆Vnom as a deterministic
value, and test it against N0 through a Z-score metric. The Z-score is a statistical measure that quantifies
the distance between a data point x and the mean of a dataset µ in terms of standard deviations:

x−µ
z= . (10)
σ
To determine the probability that x lies within a certain range of the distribution, the cumulative distribution
function (CDF) of the standard normal distribution Z, is used. Therefore the probability that a sample X
of the distribution is going to be higher than the value of x, meaning that x is contained in the distribution
itself is computed as:

P (X ≥ x) = 1 − P (X ≤ x) . (11)

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
This metric may be coupled with an additional figure of merit, widely employed in the field of probabilistic
data association, the Squared Mahalanobis Distance (SMD). The SMD is a statistical distance between two
Gaussian distributions, often used to test the association between them. It is defined as:

SM D(x) = (x − x)T (Σx + Σx )−1 (x − x) ∼ χ2 (x) , (12)

where x is a reference distribution mean value with associated covariance Σx , x represents a sample distri-
bution mean value with covariance Σx , and it is assumed that (x − x) ∼ N (0, Σx + Σx ).
The SMD is then distributed as a χ2 with nx degrees of freedom, that is the size of the multivariate function
considered, being the sum of nx normalized Gaussian distributions. Association is performed when the SMD
value is smaller than the quantile of the χ2 (nx ) having a confidence level of αL [15].
A correlation index based on the SMD and defined as PSM D = SM D(x)/χ2 can then be introduced. Differ-
ently from the first one, this second score is based on the comparison between two expense distributions:
the background one N0 and the one representing the OCP solution N (∆V ∗ , Σ∆V ∗ ), no more described as
the deterministic value ∆Vnom . For the case at hand, due to the low orders of magnitude of the involved
variable, the logarithm of the expense is employed to avoid incurring in numerical errors, as suggested in
[12].
For the computation of all the ∆V distributions involved in the aforementioned processes, both a classical
Monte Carlo simulation and an Unscented Transform (UT) are employed and compared. This is because the
former, with a sufficient number of samples, is considered as benchmark to represent any kind of distribution
but it is at the same time computationally expensive. As for the latter, it comes into play as a faster
alternative when the distribution can be described as Gaussian. In this case, once properly set, it should
generate outputs that are comparable to a Monte Carlo.

3. NUMERICAL SIMULATION

Regarding the numerical simulations both synthetic and real test cases are performed. This section is
dedicated to the description of the scenarios as well as the presentation of the related results. The orbital
regimes considered throughout the simulations are GEO and MEO.
3.1 Synthetic scenarios
For the synthetic cases, the setup remains consistent across the different orbital regimes. The process begins
by selecting a reference ground station and retrieving the objects’ Two-Line Elements (TLEs) from Space-
Track website. The TLE epoch is then set as the starting point of the observation window, which spans 48
hours. The object’s dynamics are subsequently propagated throughout this observation window, with a time
interval ∆t of 40 seconds between consecutive observations. Visibility constraints relative to the selected
ground station are also implemented, taking into account the geometric visibility and the lighting conditions
of both the target and the ground station. Each orbital state that meets these conditions is then projected
onto the measurement space to obtain observations, represented as Right Ascension and Declination pairs.
To maintain realism in the simulation, only five measurements per track are considered, given the typically
long duration of the retrieved tracks. Once the dataset is generated, track compression is performed, and a
list of attributables is computed for subsequent testing.
In these synthetic scenarios, two targets from the GEO and MEO orbital regimes are analyzed. The first
object, with NORAD ID 44475, is ESA’s European Data Relay System (EDRS), situated in the GEO regime.
The second target, with NORAD ID 43566, is a satellite from the Galileo constellation, located in the MEO
regime. The orbital elements (OE) of these targets are detailed below (with the exception of true anomaly
θ) to give an idea of their shape and orientation:

OEGEO = [42165 km, 3.7e − 5, 0.03 deg, 140.8 deg, 101 deg]
OEM EO = [29602 km, 1.4e − 4, 57.3 deg, 0.395 deg, 47.6 deg]

As an initial step, the outcomes of all pairs of measurements in the dataset are tested, with the first angular
measurement in the set fixed as a0 . Meanwhile, af varies with each iteration, spanning the remaining values.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
The objective is to evaluate the quality of the results as a function of the time interval between measurements,
different dynamical models, and the presence or absence of measurement noise.
Fig. 1 and Fig. 2 display the nominal expenses ∆Vnom for the GEO and MEO cases, respectively. The crosses
represent the scenario without uncertainties using a Keplerian model for integration, the diamonds represent
the same dynamical model with added measurement noise, the circles correspond to the scenario using the
SGP4 model, and the triangles indicate the case with both measurement noise and SGP4 integration.
As the graphs illustrate, the expenses for the Keplerian case without noise are the lowest, with the exception
of a couple of outliers. These outliers can be attributed to the quality of the initial state x0 obtained from
the IOD process, which is used as the initial guess for the iterations. Similarly, in scenarios where the SGP4
model is applied, the expenses remain low, with magnitudes comparable to the Keplerian case, meaning
perturbations do not affect the trajectory significantly, given the considered time spans and orbital data
accuracy. Here, another outlier can be noticed in Fig. 1, which can also be traced back to the x0 guess.
When measurement noise is introduced, the results begin to change. In these scenarios, noise corresponding
to the sensor’s accuracy is added to the measurement dataset. For both Keplerian and SGP4 models,
the resulting expenses increase, though not uniformly across all attributable pairs. Two main observations
can be drawn from this: first, in some instances, the Keplerian model shows higher expenses compared
to the SGP4 model. This may occur because the measurement noise sometimes counteracts the effects of
dynamical perturbations, especially in orbital regimes that are relatively free from high-impact perturbations
like atmospheric drag. Second, the impact of measurement noise is closely tied to the random values selected
from the normal distribution generated for the noise. When these values are close to the mean, the final
expenses tend to resemble those in scenarios where noise is not included.

Fig. 1: GEO object ∆Vnom distribution in time for Fig. 2: MEO object ∆Vnom distribution in time for
different combinations of dynamical model and different combinations of dynamical model and
measurements noise. measurements noise.

Following the previous considerations, the assessment of the technique’s behavior is presented through an
analysis of the best-case and worst-case scenarios.
Starting with the Keplerian case without associated noise, both Fig. 3a and Fig. 4a indicate that the distance
metric-based correlation probability for each dataset pair is high. A comparison between a Monte Carlo
approach and an UT method is also presented. This choice is primarily motivated by the need to evaluate the
quality of the results using an UT-based uncertainty propagation, which is less computationally demanding
than a traditional Monte Carlo approach. As shown in the graphs, the results are quite comparable: in
some cases, the probabilities derived from both methods coincide, while in others, the UT method tends
to slightly overestimate or underestimate the probability. However, even in these cases, the correlation
probability remains well above 50%, considered as a reasonable threshold for the correlation test. It is also
worth noting that the presented results do not appear to be significantly influenced by the time interval
between the two considered attributables. To further validate the correctness of the correlation hypothesis,
an additional check using the alternative SMD metric is conducted. As shown in Fig. 3b and Fig. 4b, the

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
SMD index for each case remains well below the prescribed threshold, corresponding to 1 as in a typical χ2
test. The lower the value of this kind of index is, the more likely the association between BCs. This means
that the proposed test cases confirm the correlation process is reliable. For this evaluation, only the UT
approach was employed, as it is proven comparable to the Monte Carlo in terms of accuracy, while being
less computationally intensive.
In the more realistic scenario, which includes dynamic propagation using SGP4 and embedded measurement
noise, Fig. 5a and Fig. 6a demonstrate, as expected, a lower Z-score-based correlation probability. Despite
this, in every case, the probability calculated using both the Monte Carlo and UT methods remains above
50%. The results show a pattern similar to previous cases, where the UT outcomes do not exactly align with
the Monte Carlo results. Specifically, in the GEO case of Fig. 5a, the UT method generally underestimates
the correlation probability, while in the MEO case of Fig. 6a, it tends to overestimate it overall. For this
scenario, an SMD-based correlation test is also conducted, with the results shown in Fig. 5b and Fig. 6b.
The threshold is complied in these cases as well, further increasing confidence in the correctness of the
measurement association procedure.

Fig. 3a: GEO Z-score-based correlation probability Fig. 3b: GEO SMD-based correlation index
with a Monte Carlo and UT approach (keplerian leveraging UT approach (keplerian dynamics, no
dynamics, no noise). noise).

Fig. 4a: MEO Z-score-based correlation probability Fig. 4b: MEO SMD-based correlation index
with a Monte Carlo and UT approach (keplerian leveraging UT approach (keplerian dynamics, no
dynamics, no noise). noise).

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Fig. 5a: GEO Z-score-based correlation probability
Fig. 5b: GEO SMD-based correlation index
with a Monte Carlo and UT approach (sgp4
leveraging UT approach (sgp4 dynamics, noise).
dynamics, noise).

Fig. 6a: MEO Z-score-based correlation probability

Fig. 6b: MEO SMD-based correlation index
with a Monte Carlo and UT approach (sgp4
leveraging UT approach (sgp4 dynamics, noise).
dynamics, noise).

3.2 Real scenarios

A fundamental step in validating the technique is testing it with real data. It is important to note that
using real data introduces limitations in flexibility, as the track frequency available may not match that of
the simulated cases. In this analysis, TDMs from different ground stations, employing various sensors with
differing accuracies, are employed. The datasets are characterized as follows:

a) GEO case
• NORAD ID 44034: Single TDM containing 8 short tracks spanning a 12-hour time window.
• NORAD ID 19397: Two TDMs separated by a 1-day period, each covering a 2-hour observation
window.
b) MEO case
• NORAD ID 41861: Single TDM containing observations spanning a 3-hour window.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 • NORAD ID 43567: Two TDMs separated by an almost 3-day gap. The tracks in both TDMs
span less than 3 hours.

The analysis outcomes are presented in the same order as the datasets listed above. Starting with the GEO
case of NORAD ID 44034, the technique’s performance assessment on the single TDM is shown in Fig. 7a
and Fig. 7b. For all the real cases, the attributable a0 retrieved from the first track remains fixed, while af
iterates through the remaining tracks within the TDM being tested.
Fig. 7a compares the Z-score correlation probabilities derived from both the Monte Carlo-based and UT-
based methods. In this instance, the UT formulation tends to slightly underestimate the probability, but
overall, all tested cases still yield good results, with probabilities exceeding 50%. An SMD-based check,
reported in Fig. 7b, further confirms the correctness of the measurement association hypothesis, as all index
values satisfy the imposed threshold of 1.
To investigate the method’s behavior with increasing time interval between observations, the information
from the two TDMs associated with NORAD ID 19397 is used. Since the individual TDMs cover observation
periods that are similar to the previous simulation in magnitude, it is deemed more insightful to present
results corresponding to the one-day interval. This analysis thus includes only one case, summarized in
Tab. 1 in terms of results. Correlation probability is lower with respect to the previous cases, but still within
the acceptability limits and comparable within the two methodologies applied. Similarly, the PSM D satisfies
the threshold thus confirming correlation.

100

90 1

80
0.8
70

60
0.6
50

40
0.4
30

20 0.2

10

0 0
0 1 2 3 4 5 6 7 8 9 0 2 4 6 8 10 12

Fig. 7a: NORAD ID 44034 Z-score-based correlation Fig. 7b: NORAD ID 44034 SMD-based correlation
probability with a Monte Carlo and UT approach. index leveraging UT approach.

∆t [hours] PMC [%] PUT [%] idxSMD [-]

NORAD ID 19397 25.2 63.44 67.04 0.2485

Table 1: NORAD ID 19397 Z-score-based correlation probability and SMD-based correlation index.

Moving onto the MEO dataset, the first presented results are referred to NORAD ID 41861. Since the
observations contained in the TDM are divided by a gap of one hour, the cases adding the most information
about the method’s sensitivity to time span are two, featuring a ∆t of 1.8 hours and 2.5 hours respectively.
The technique’s performance assessment on the single TDM is shown in Fig. 8a and Fig. 8b. Fig. 8a highlights
an almost complete overlap of Monte Carlo-based and UT-based Z-score correlation probabilities for both
case studies. As for the order of magnitude of the correlation probabilities involved, it is comparable to
the previous UT-based ones, as reported in Fig. 7a. Nonetheless, this time the probability associated with
the Monte Carlo simulation is reduced too. Overall though, all tested cases still yield good results, with
probabilities exceeding the target 50%. The SMD-based check, reported in Fig. 8b, confirms once again that
the measurement association hypothesis is correct, as all index values lie below 1.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 100

90 1

80
0.8
70

60
0.6
50

40
0.4
30

20 0.2

10

0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3

Fig. 8a: NORAD ID 41861 Z-score-based correlation Fig. 8b: NORAD ID 41861 SMD-based correlation
probability with a Monte Carlo and UT approach. index leveraging UT approach.

As the last case, the analysis involving the two TDMs belonging to NORAD ID 43567 is presented. Also
for this scenario, the observations composing the two TDMs span a 2-hour-long window, thus the testing is
carried out between the two TDMs rather than analyzing the tracks composing them. This time, the time
interval between the considered attributable is further increased. Tab. 2 shows that, even if the ∆t is of
almost 3 days, results are still more than satisfactory both from the probability and the index perspectives.

∆t [hours] PMC [%] PUT [%] idxSMD [-]

NORAD ID 43567 71.5 82.97 74.09 0.0161

Table 2: NORAD ID 43567 Z-score-based correlation probability and SMD-based correlation index.

4. CONCLUSIONS

This paper introduces a novel technique based on a minimum energy optimal control problem to address
track-to-track association, with the goal of enhancing catalog maintenance methodologies. The main focus
is testing the ballistic dynamics hypothesis to correlate two UCTs by comparing the OCP solution expense,
both as a deterministic value and a distribution, with an approximate background distribution, connected
to boundary condition uncertainty. The application of transversality conditions simplifies the problem,
eliminating the need for the formulation of an admissible region and relaxing possible assumptions on the
underlying dynamics. The reported tests highlight the impact of different dynamic models, in terms of
considered perturbations, and the presence of measurement noise on method’s correlation performance.
Uncertainty is embedded in the process via two approaches: Monte Carlo simulation and UT. While the UT-
based approach shows occasional deviations from Monte Carlo results, both procedures proved as reasonable
methods for uncertainty propagation in this specific application. To ensure robust correlation, two metrics,
the Z-score and the more stringent SMD, are formulated and analyzed side by side. Both of them yield
promising results across all scenarios. Future developments aim to extend this technique to lower orbital
regimes, enhance the accuracy of UT-derived distributions, and test the methodology with a larger dataset
of real measurements to identify potential weaknesses. Additionally, the technique will be adapted to include
maneuver detection, extending its application to scenarios involving actively controlled targets. Additional
work is also underway to tune the technique for cislunar space, in an attempt to assess its performance when
a different orbital dynamics is involved.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 ACKNOWLEDGEMENTS

This research has received funding as part of the work developed for the agreement n. 2023-37-HH.0 for the
project “Attività tecnico-scientifiche di supporto a C-SSA/ISOC e simulazione di architetture di sensori per
SST”, established between the Italian Space Agency (ASI) and Politecnico di Milano (POLIMI).

REFERENCES

[1] D. J. Kessler and B. G. Cour-Palais. Collision frequency of artificial satellites: The creation of a debris
belt. Journal of Geophysical Research: Space Physics, 83(A6):2637–2646, 1978.
[2] IADC. Iadc space debris mitigation guidelines. Technical report, Inter-Agency Space Debris Coordina-
tion Committee, 2020.
[3] ESA Space Debris Office. Esa’s annual space environment report. Technical report, ESA, 2024.
[4] L.G. Taff and D.L. Hall. The use of angles and angular rates: I: Initial orbit determination. Celestial
mechanics, 16(4):481–488, 1977.
[5] G. F. Gronchi, L. Dimare, and A. Milani. Orbit determination with the two-body integrals. Celestial
Mechanics and Dynamical Astronomy, 107:299–318, 2010.
[6] A. Milani, G. F. Gronchi, M. De’ Michieli Vitturi, and Z. Knežević. Orbit determination with very
short arcs. i admissible regions. Celestial Mechanics and Dynamical Astronomy, 90:57–85, 2004.
[7] G. Tommei, A. Milani, D. Farnocchia, and A Rossi. Correlation of space debris observations by the
virtual debris algorithm. In Proc. of the Fifth European Conference on Space Debris, volume 30, 2009.
[8] L. Pirovano, R. Armellin, J. Siminski, and T. Flohrer. Data association for too-short arc scenarios with
initial and boundary value formulations. In 20th AMOS Conference, September, pages 17–20, 2019.
[9] K. J. DeMars and M. K. Jah. Probabilistic initial orbit determination using gaussian mixture models.
Journal of Guidance, Control, and Dynamics, 36(5):1324–1335, 2013.
[10] M. J. Holzinger, D. J. Scheeres, and K. T. Alfriend. Object correlation, maneuver detection, and charac-
terization using control distance metrics. Journal of Guidance, Control, and Dynamics, 35(4):1312–1325,
2012.
[11] A. Pastor, G. Escribano, and D. Escobar. Satellite maneuver detection with optical survey observations.
In Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS), 2020.
[12] R. Serra, C. Yanez, and C. Frueh. Tracklet-to-orbit association for maneuvering space objects using
optimal control theory. Acta Astronautica, 181:271–281, 2021.
[13] J. Siminski, H. Fiedler, and T. Flohrer. Correlation of observations and orbit recovery considering
maneuvers. AAS/AIAA Space Flight Mechanics, 2017.
[14] A. Milani, M. E. Sansaturio, and S. R. Chesley. The asteroid identification problem iv: Attributions.
Icarus, 151(2):150–159, 2001.
[15] L. Pirovano, R. Armellin, J. Siminski, and T. Flohrer. Differential algebra enabled multi-target tracking
for too-short arcs. Acta Astronautica, 182:310–324, 2021.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com