# Physics and Human-Based Information Fusion for Improved Resident Space Object Tracking

Maintaining a catalog of Resident Space Objects (RSOs) can be cast in a typical Bayesian multi-object estimation problem, where the various sources of uncertainty in the problem - the orbital mechanics, the kinematic states of the identified objects, the data sources, etc. - are modeled as random variables with associated probability distributions. In the context of Space Situational Awareness, however, the information available to a space analyst on many uncertain components is scarce, preventing their appropriate modeling with a random variable and thus their exploitation in a RSO tracking algorithm. A typical example are human-based data sources such as Two-Line Elements (TLEs), which are publicly available but lack any statistical description of their accuracy. In this paper, we propose the first exploitation of uncertain variables in a RSO tracking problem, allowing for a representation of the uncertain components reflecting the information available to the space analyst, however scarce, and nothing more. In particular, we show that a human-based data source and a physics-based data source can be embedded in a unified and rigorous Bayesian estimator in order to track a RSO. We illustrate this concept on a scenario where real TLEs queried from the U.S. Strategic Command are fused with realistically simulated radar observations in order to track a Low-Earth Orbit satellite.

There are no comments yet.

## Authors

• 3 publications
• 13 publications
• 1 publication
• ### Multi-modal Tracking for Object based SLAM

We present an on-line 3D visual object tracking framework for monocular ...
03/14/2016 ∙ by Prateek Singhal, et al. ∙ 0

• ### Robust video object tracking via Bayesian model averaging based feature fusion

In this article, we are concerned with tracking an object of interest in...
04/02/2016 ∙ by Yi Dai, et al. ∙ 0

• ### On the loss of Fisher information in some multi-object tracking observation models

The concept of Fisher information can be useful even in cases where the ...
03/26/2018 ∙ by Jeremie Houssineau, et al. ∙ 0

• ### Learnable Online Graph Representations for 3D Multi-Object Tracking

Tracking of objects in 3D is a fundamental task in computer vision that ...
04/23/2021 ∙ by Jan-Nico Zaech, et al. ∙ 0

• ### A Common Operating Picture Framework Leveraging Data Fusion and Deep Learning

Organizations are starting to realize of the combined power of data and ...
01/16/2020 ∙ by Benjamin Ortiz, et al. ∙ 10

• ### Exploiting Source-Object Network to Resolve Object Conflicts in Linked Data

Considerable effort has been made to increase the scale of Linked Data. ...
04/28/2016 ∙ by Wenqiang Liu, et al. ∙ 0

• ### Toward Biochemical Probabilistic Computation

Living organisms survive and multiply even though they have uncertain an...
11/09/2015 ∙ by Jacques Droulez, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Space Situational Awareness (SSA) concerns itself with having the actionable knowledge to predict, deter, avoid, operate through, recover from, or attribute cause to the loss, disruption, degradation, or denial of space services, capabilities, or activities. Maintaining a catalog of so-called RSOs can be cast as a typical multi-object detection and tracking problem. The individual state of each object consists of kinematic and relevant model parameters to the construction of the catalog and/or its exploitation for various activities related to space, such as mission planning, conjunction analysis to prevent collisions, etc. Typically, the individual state of an orbiting object describes a kinematic relationship (Kepler elements, or position/velocity coordinates), but may include additional characteristics such as its attitude, its ballistic coefficient, etc.

The main challenge to the resolution of most detection and tracking problems is to appropriately address the various sources of uncertainty involved in its different components, due to the limited knowledge in the signal processing chain affecting the acquisition of information on the population of objects. The SSA problem is no different, and the relevant sources of uncertainty can be sorted in two broad categories:

1. To first order, the perturbing accelerations on the RSO population are known, gravity being the dominant perturbing source. However, no orbital propagator is able to produce the exact dynamical state of an RSO, since the various physical perturbations affecting orbiting trajectories are either approximated, discarded on purpose for the sake of algorithmic efficiency, or simply unknown to the space analyst. In addition, satellites owned by others may transition to new orbital paths unknown to the analyst. New RSOs entering some specific orbital region, either through a collision or a new launching event, are not all accounted for and thus the number of RSOs remain uncertain.

2. The mechanisms through which observations (or opinions) are produced by the data sources are partially known, at best. Traditional, physics-based sensors like radars, telescopes, or cameras are subject to missed detections, false positives, observation noise, from which some statistical description is usually available. On the other hand, the TLEs produced by the U.S. Strategic Command (USStratCom) are publicly available, but no information regarding the accuracy of the orbital elements or the confidence in the object labeling are available.

The Bayesian estimation framework is a very popular approach to solve multi-object detection and tracking problems, and the origin of the overwhelming majority of modern tracking algorithms. Its key features are that

1. the uncertainty on each component of the system is modeled with a random variable and characterized with a probability distribution, and

2. a probability distribution can be updated sequentially, or filtered, with the availability of new information regarding the corresponding random variable.

Bayesian filters then maintain a probabilistic description of the population of objects and propagate it through time, updating it whenever new observations are collected from a data source, also described through a probabilistic model. The characteristics of a specific Bayesian filter primarily depends on the nature of the probabilistic description representing the population of objects; most notably, through individual tracks (Blackman, 2004), labeled random finite sets (Vo and Vo, 2013), or stochastic populations (Houssineau and Clark, 2016; Delande et al., 2016). Applications of these various multi-object filtering techniques to the context of SSA can be found in DeMars et al. (2012), in Jones and Vo (2015); Jones et al. (2016), and in Delande et al. (2018, 2017), respectively.

One of the key challenges of Bayesian estimation is to maintain an appropriate representation of the uncertainty on the components of the system, especially when information about their state is scarce. A typical example in the context of SSA is the initial orbit determination procedure through which the possible values for the dynamical state (position, velocity coordinates) of a newly-detected RSO are restricted to a specific subspace known as the admissible region in DeMars and Jah (2013). The physical considerations behind the initial orbit determination procedure do not provide information on whether any dynamical state within the admissible region is likelier than the other, but it does not indicate that all states within the admissible region are equally likely either111It turns out that they are not, since a more restrictive initial orbit determination procedure can lead to a constrained admissible region (DeMars et al., 2012).. However, a probabilistic interpretation of the RSO’s initial state can lead to a uniform probability distribution on the admissible region, thus producing a description of the RSO’s state that is not inferred from the information available to the analyst. It thus appears that probability density functions are inappropriate tools to describe admissible regions, as argued in a recent study (Worthy III and Holzinger, 2017b). Another example is the exploitation of human-based or semantic data sources such as TLEs or natural languages statements: they could provide a wealth of information regarding RSOs, but the lack of statistical information on their accuracy/truthfulness makes their probabilistic representation, and thus their integration to a Bayesian tracking filter, difficult and largely unexplored to this day.

Alternatives to the standard probabilistic representation of uncertainty exist, such as fuzzy logic, imprecise probabilities, possibility theory, fuzzy random sets, and Dempster-Shafer theory (Zadeh, 1965; Walley, 1991; Dempster, 1967; Shafer, 1976; Dubois and Prade, 1983; Yen, 1990; Friedman and Halpern, 2001); recently, the Dempster-Shafer theory has been exploited to approach the initial orbit determination problem with admissible regions in Worthy III and Holzinger (2017a). While these methods have been widely used to describe the uncertainty referring to a fixed unknown state, their exploitation to the estimation of dynamical systems is less straightforward, limiting their applicability to sequential estimation problems – that is, to the design of detection and tracking filters.

A recent alternative is given by the outer probability measures as introduced in Houssineau and Bishop (2017). They aim at proposing a “prejudice-free” representation of the uncertainty on a system through a natural construction that is derived from the available information, and nothing more. Built from fundamental tools of measure theory, o.p.m.s are compatible with the Bayesian filtering framework and can be integrated in a complex system where some uncertain components are described with pdfs, while others are described with o.p.m.s.

In this paper, we will focus on the representation of uncertainty for SSA data sources supporting a RSO tracking algorithm. More specifically, we will show that radar observations (physics-based information) and TLEs (human-based information) can both be represented by o.p.m.s and integrated into a Bayesian filtering algorithm, where the information on the RSO’s state is maintained with a probability distribution. Sec. 2 presents the concepts of o.p.m. and possibility function, and Sec. 3 describes the filtering equations of the single-object Bayesian tracking algorithm exploiting possibility functions. Sec. 4 focuses on the modeling of the data sources relevant to this paper, i.e., a radar with Doppler effect and a “TLE-generator”. Then, Sec. 5 describes the construction of a target tracking scenario with realistically-simulated radar observations and real TLE data provided by the USStratCom, and gives a detailed implementation of the single-target tracking algorithm. Sec. 6 presents the filtering results. A discussion on future works follows in Sec. 7, and Sec. 8 concludes.

## 2 Uncertain variables and o.p.m.s

### 2.1 Random and non-random uncertainty

Assume some system, whose state is described by some state space , on which some operator – say, a space analyst – possesses some knowledge. We shall decompose the uncertainty on the system’s state in two different parts:

• The random uncertainty is the inherent component, that exists as a natural disposition of the system and independently of the analyst, sometimes known as the aleatoric uncertainty,

• The non-random uncertainty is the extraneous component, that exists only with respect to (w.r.t.) the information possessed by the analyst, sometimes known as the epistemic uncertainty.

The random component can be seen as the remaining uncertainty on a system’s state, once the analyst possesses all the information on the system that could be acquired. It is unclear if physical systems described on the macroscopic level, such as an RSO’s dynamical state, can be anything but deterministic, or if there are random effects that cannot be explained assuming perfect knowledge of the physical laws governing the orbital motion model. It is, in the end, irrelevant to the estimation problem that focuses on reducing the non-random – or epistemic – uncertainty.

### 2.2 Possibility functions

The random uncertainty, inherent to the system, shall be characterized by a probability distribution on the state space .222For the sake of simplicity, we shall assume that admits a density w.r.t. some reference measure on , also called , and we shall handle pdfs rather than probability distributions, whenever applicable, in the rest of the paper. The non-random uncertainty reflects the information on the system’s state possessed by the analyst. In the simplest case, it is maintained by a single function satisfying

 ∫Bp(x)dx≤maxx∈Bf(x), (1)

for any appropriate333The subset must be measurable, for the integral to be well-defined. For the sake of simplicity this requirement on subsets, whenever applicable, will be omitted in the rest of the paper. subset of . The function is assumed to have maximum for the sake of convenience. In other words, while the pdf describing the system’s state is inaccessible to the analyst, they maintain an upper bound reflecting the information they possess about the system. A simple example of upper bound is given in Fig. 1.

The non-negative, measurable functions on with maximum are also called possibility functions (Houssineau, 2015), or possibility distribution as in Dubois and Prade (1983). Box-shaped possibility functions such as illustrated in Fig. 1 are convenient to represent negative evidence of the system’s state, i.e., to exclude regions of the state space as possible values without inferring anything else – most notably, without assuming that all the remaining values are equally probable.

In many situations, such as the modeling of a physics-based sensor, more refined possibility functions are necessary. Similarly as for probability distributions, we define Gaussian possibility functions as functions on , , of the form

 ¯N(x;μ,S)=exp(−12(x−μ)tS−1(x−μ)),

for some and some positive-definite matrix , where denotes the matrix transposition. For the sake of convenience, we will refer to and

as the mean and variance of

, although these concepts do not bear the statistical meaning of their namesakes for probability distributions. A family of Gaussian possibility functions is illustrated in Fig. 2.

Fig. 2 also illustrates two particular possibility functions. The indicator function denotes full knowledge about the system, known to have state almost surely, since the pdf with mass concentrated in is the only pdf on satisfying (1) when . Conversely, the indicator denotes the total absence of information about the system, since every pdf on satisfy (1) when . The latter possibility represents a truly uninformative uncertainty, for one does not infer anything regarding the system’s state, and in particular does not infer that all the states in are equally probable: indeed, the uniform probability distribution on does satisfy (1) when , but so do any other probability distribution on .

### 2.3 Uncertain variables and o.p.m.s

In the general case, the information on some system’s state maintained by an analyst can assume a complex form built from a combination of possibility functions. More formally, it can be represented through a probability distribution on the space of possibility functions such that

 ∫Bp(x)dx≤∫maxx∈Bf(x)P(df), (2)

for any subset of . In practical problems, a finite combination of possibility functions is often sufficient to describe the uncertainty regarding many systems; in the scope of this paper, we will see in Sec. 4 that a single possibility function is enough to describe various data sources relevant to the context of SSA.

The information given by induces an outer probability measure on given by

 ¯P(B)=∫maxx∈Bf(x)P(df), (3)

for any subset of . The scalar reads as the possibility that the system’s state is in . This is to be contrasted with the probability measure characterizing the inherent randomness of the object’s state, independently of the information possessed by the analyst, and given by

 P(X∈B)=∫Bp(x)dx. (4)

From (2) it follows that

 1−¯P(X∖B)≤P(X∈B)≤¯P(B), (5)

for any subset of , where denotes the set difference. In other words, an analyst with information represented by cannot (in the general case) evaluate the probability for the object to lie within , but they may bound its value.

Note that an o.p.m. represents the uncertainty of an analyst in the state of some system that might not be random. In that case, the system has a fully deterministic behavior and the mass of the pdf in (4) is concentrated to some state – unknown to the analyst, as they are unaware of the nature of . For this reason, the notion of uncertain variables is introduced as a replacement for random variables to describe uncertain systems (Houssineau and Bishop, 2017).

### 2.4 A simple example of o.p.m.s

We shall describe two distinct levels of information, supported by the same subsets , of the state space illustrated in Fig. 3.

Suppose first that the information maintained by the analyst is based solely on the possibility function – that is, the probability measure in (3) gives weight to the singleton . Following (3), the information is represented by the o.p.m.

 ¯P1(A)=[A∩(B∪B′)≠∅], (6)

where is the Iverson bracket444 is the event is true, and zero otherwise.. The analyst wishes to gauge the probability that the object’s state lies within : from Eq. (6) we have

 ¯P1(X∖(B∪B′))=0,¯P1(B∪B′)=1,

and from Eq. (5) the analyst concludes that

 1≤P(X∈B∪B′)≤1,

that is, the object lies in almost surely. The analyst then wishes to inquire about the presence of the object in some smaller region : from Eq. (6) we have

 ¯P1(X∖C)=1,¯P1(C)=1,

and from Eq. (5) the analyst concludes that

 0≤P(X∈C)≤1,

i.e., the analyst is clueless about the distribution of the object within . In conclusion, the analyst infers that the object is in with certainty, but nothing more.

Suppose now that the information maintained by the analyst is based on the two possibility functions , , with respective weights and , where . Following (3), the information is now represented by the o.p.m.

 ¯P2(A)=α[A∩B≠∅]+(1−α)[A∩B′≠∅]. (7)

The analyst wishes to gauge the probability that the object’s state lies within : from Eq. (7) we have

 ¯P2(X∖B)=1−α,¯P2(B)=α,

and from Eq. (5) the analyst concludes that

 α≤P(X∈B)≤α,

that is, the object lies in with probability . The analyst then wishes to inquire about the presence of the object in some smaller region : from Eq. (7) we have

 ¯P2(X∖C)=1,¯P2(C)=α,

and from Eq. (5) the analyst concludes that

 0≤P(X∈C)≤α.

i.e., the analyst is clueless about the distribution of the object within . A similar reasoning applies in region ; in conclusion, the analyst infers that the object is in with probability , in with probability , but nothing more.

### 2.5 From o.p.m.s to pdfs

We see in the example illustrated in the previous section that provides a more refined level of information on the system’s state than , because it imposes tighter constraints on the system’s underlying pdf . In the extreme case where an o.p.m. imposes constraints on all the singletons of the state space, i.e., when the associated probability measure in (3) has support in the possibility functions of the form , where , then it can be shown that induces a pdf on such that

 ∫Bp∞(x)dx=∫maxx∈Bf(x)P∞(df)=¯P∞(B), (8)

for any subset of . That is, the o.p.m. is equivalent to the probability measure in (4), and thus the uncertainty of the analyst in the system’s state reduces to the system’s inherent random component characterized by . In other words, the analyst has nothing more to learn about the system’s state. In particular, if has its mass concentrated in some point , the analyst concludes that the system is deterministic and has state almost surely.

This connection between o.p.m.s and pdfs is of practical importance. It highlights the fact that the distinction between random and non-random uncertainty in some component of a complex system can be purposely ignored for the sake of simplicity; this component is then represented with a usual random variable and described with a level of information that amounts to a probability distribution. As we shall see in Sec. 3, it allows for the representation of each uncertain component of a detection and tracking problem – object’s state, prediction model, observation model specific to each sensor – with an uncertain variable, that may or may not reduce to a random variable, in a single coherent Bayesian filtering framework.

## 3 Bayesian filtering with uncertain variables

### 3.1 Conditional possibility functions and o.p.m.s

The formulation of Bayesian inference algorithms rely on the definition of appropriate conditional

pdfs describing the motion and the observation of the object of interest. It is then natural to introduce conditional possibility functions of the form describing an uncertain variable given the realization in a state space of another uncertain variable and verifying

 maxx∈Xf(x|x′)=1

for any . Conditional possibility functions verify the same sort of property as conditional pdfs: letting be a possibility function describing , it holds that

 f(x)=maxx′∈X′f(x|x′)f′(x′) (9)

is a possibility function describing the state of in , and that

 f′(x′|x)=f(x|x′)f′(x′)maxx′∈X′f(x|x′)f′(x′) (10)

is a possibility function describing the state of given a realization of . Note that Eq. (9) is the analogue of the prediction formula, Eq. (10

) is the analogue of Bayes’ theorem, for possibility functions. Following the approach considered in

Houssineau (2018), conditional o.p.m.s can also be defined as

 ¯P(B|X′=x′)=∫maxx∈Bf(x|x′)P(df|X′)

where is a probability measure on the space of conditional possibility functions555The conditioning on indicates the nature of the possibility functions in the support of but does not make it a function of ..

### 3.2 Bayesian filtering equations for o.p.m.s

Suppose that an analyst studies some object with state in some space , described at time ( respectively (resp.) ) by the uncertain variable (resp. ). Besides, some sensor with observations666Since TLEs are not generated by a physics-based sensor, a more accurate term to describe one TLE point would perhaps be “opinion”. For the sake of simplicity, the terms “sensor” and “observation” are usually employed throughout this paper regardless of the nature of the data source. in some space observes the scene, and is described by the uncertain variable .

Suppose that the information that the analyst possesses on the object’s state at time is represented by an o.p.m. , and the information they possesses on the object’s evolution between time and is represented by a conditional o.p.m. , depending on a realization of . The predicted information the analyst possesses on the object’s state at time is then described by the o.p.m.

 ¯Pk|k−1(B)=∫maxx∈B(maxx′∈X′f(x|x′)f′(x′))Mk(df|Xk−1)Pk−1(df′), (11)

Suppose that the information that the analyst possesses on the sensor’s observation process is represented by a conditional o.p.m. , depending on a realization of , and that an observation is collected from the sensor. The updated information the analyst possesses on the object’s state at time is then described by the o.p.m.

 ¯Pk(B|Yk=y)=∫maxx∈B(f(y|x)f′(x))Ok(df|Xk)Pk|k−1(df′)∫maxx∈X(f(y|x)f′(x))Ok(df|Xk)Pk|k−1(df′). (12)

The equations (11), (12) are the analogues of the time prediction and data update equations for pdfs, but for o.p.m.s; that is, they are the Bayesian filtering equations for an estimation problem in which the three sources of uncertainty – the object’s state, the object’s evolution, the sensor’s observation process – are described with arbitrary uncertain variables, that may or may not reduce to random variables.

Detection and tracking problems in the context of SSA present many challenges, whose resolution with Bayesian filters may benefit from approaching these three types of uncertainty with uncertain variables (this is briefly discussed in Sec. 7). In the context of this paper, however, we shall simply illustrate the benefit from modeling the sensors’ observation process with a single possibility function; in particular, both the information on the object’s state and the object’s evolution will be represented with usual pdfs.

### 3.3 A simpler case of Bayesian filtering equations

We shall suppose, then, that the o.p.m. on the object’s state at time reduces to a pdf on , and the o.p.m. on the object’s evolution reduces to a transition kernel on , depending on a realization of . The prediction formula (11) then reduces to

 pk|k−1(x)=∫mk(x|x′)pk−1(x′)dx′, (13)

where is the pdf describing the predicted information on the object’s state. As expected, we obtain the usual prediction formula for random variables.

We shall also suppose that the o.p.m. , representing the information the analyst possesses on the sensor’s observation process, is based on a single possibility function on . The update formula (12) then reduces to

 pk(x|y)=hk(y|x)pk|k−1(x)∫hk(y|x′)pk|k−1(x′)dx′, (14)

where is the pdf describing the updated information on the object’s state777The explicitly dependency of over past observations is omitted here for the sake of simplicity. As expected, we obtain an update equation with a very similar structure than the usual Bayes’ update rule, except that the observation process is described by a possibility function on rather than a likelihood function , i.e., a pdf on .

## 4 Modeling of Ssa data sources with o.p.m.s

In this section, we shall focus on the modeling of two data sources relevant to the context of SSA, namely, a radar ground station and USStratCom’s catalog. Both data sources provide point-like observations in a well-defined observation space. They differ, however, on the available information about their observation process: a statistical description of the accuracy of a radar-like sensor is usually available, while no measure of uncertainty is provided (at least publicly) alongside TLEs. We will see in that section that both data sources can be represented in a natural way with a single possibility function.

### 4.1 Radar observations

We shall simulate, in the scope of this paper, a radar ground station with Doppler effect. Suppose that a radar is exploited at some time relevant to the scenario. The observation space describes a radar-like return , providing range, azimuth, elevation, and range rate coordinates in the topocentric frame centered on the radar ground station on the surface of Earth.

As usual for radar- or telescope-like sensors, the analyst is assumed to have access to some statistical description of the accuracy of the observation process, given by some positive-definite covariance matrix . The possibility function on in Eq. (14) is simply given by the Gaussian possibility

 hrad.k(y|x)=¯N(y;ok(x),Sk) (15a) =exp(−12(y−ok(x))tS−1k(y−ok(x))), (15b)

where is the mapping transforming coordinates from the target state space to the observation space . Note that this possibility function does not characterize the residual error of the radar observation process, a complex phenomenon that involves physical properties of the radar that may be known from the manufacturer, but also unknown perturbations due to the atmosphere, the object’s shape, etc. Perhaps the radar observation is fully deterministic, or perhaps it possesses some aleatoric uncertainty that even an observer with perfect knowledge of the radar, the object, and the middle between them would still be unable to eliminate. In any case, a characterization of the radar’s noise profile is unavailable to the analyst, who is only able to bound the error in the observation process with (15).

Note that the possibility (15) is proportional to the likelihood

i.e., the probabilistic description usually adopted for a radar described by , characterizing the radar’s noise profile. The difference is important: while the possibility is a dimensionless function on , the likelihood is a pdf on and has dimension as the inverse of the unit volume of the observation space . The value of the likelihood in Eq. (16) thus scales with the changes on the reference measure on – say, if the azimuth and elevation angles are measured in arc-seconds rather than in degrees – while the possibility in Eq. (15) does not (Delande et al., 2017). This has important consequences in the context of multi-target detection and tracking problems, as briefly discussed in Sec. 7.

### 4.2 TLEs

We shall query, in the scope of this paper, real TLE data produced by the USStratCom888Available at https://www.space-track.org.. For the sake of simplicity, we shall convert the original orbital elements of the TLEs to Cartesian coordinates in the reference Earth-Centered Inertial (ECI) frame, so that a TLE is considered as a point-like observation in the observation space . Unlike radar-like observations, little is known about the generation process behind the TLEs, let alone a statistical description of their residual error, i.e., the discrepancy between a RSO’s dynamical state and the generated TLE .

In order to produce some quantitative description of the TLEs’s accuracy, we performed a brief analysis on the satellite 0E0E operated by Planet Labs, on the -day-long period starting from Sep 2, 2017 at midnight. Using a batch least-square method101010Details on the orbital propagator employed are given in the description of the tracking algorithm in Sec. 5., an orbital trajectory for that satellite was produced from Global Positioning System (GPS) points provided by Planet Labs, considered as ground truth for the analysis to come. We then queried the USStratCom for the TLE points available relating to that satellite over that period, and studied the discrepancies of the TLEs as explained in the following paragraphs.

The natural observation space of a TLE is six-dimensional, since it provides information on the full dynamical state of the observed RSO. In this paper, we shall propose a simple possibility function of the observation process that aggregates the information extracted from a TLE into two meaningful dimensions, capturing salient features of the RSO’s orbital regime. Both the specific angular momentum and the specific orbital energy are conservative quantities of an orbital trajectory, in the two-body problem and assuming Keplerian motion. Denoting by (resp.

) the vector of position (

resp. velocity) coordinates of the RSO with state , and are given by

 h(x) =p×v, (17a) ϵ(x) =−μe||p||+||v||22, (17b)

where denotes Earth’s gravitational constant.

Because TLEs are expected to provide accurate information on the RSO’s orbital plane, the first component on which we shall build the possibility function is the normalized dot product

 δang.(y|x)=h(y)||h(y)||⋅h(x)||h(x)||, (18)

quantifying the angle offset between the orbital planes described by some TLE and some RSO’s state . The study of the angle offset (18) on the training data is depicted in Fig. 4.

The analysis of the angle offset (18) on the training period (Fig. 4) suggests that the TLEs do indeed provide accurate information on the RSO’s orbital plane, with an extremely small deviation of magnitude below the nominal (and maximal) value of . We can then represent this information with a possibility function on the interval as shown in Fig. 5.

The possibility function in Fig. 5 reflects our acquired information on the observation process; since it is limited in scope, we follow a cautious approach and model a possibility that has the maximum value of around the values consistent with the training data, but decreases linearly towards zero over a larger region that covers five times the deviation observed on the training data.

Similarly to the radar model, this possibility function does not characterize the residual error in the RSO’s orbital plane as given by a TLE point, like a typical likelihood function would. Such a descriptive model remains largely inaccessible through the modest analysis in Fig. 4, limited to a small period of time only and to a single satellite only. Whether a statistical profile of a TLE’s residual error could be inferred from a large scale analysis is an interesting question beyond the scope of this paper.

A similar analysis can be performed on the specific orbital energy in Eq. (17b). The second component on which we shall build the possibility function is the energy offset

 δen.(y|x)=−μe(||p(y)||−1−||p(x)||−1)+12(||v(y)||2−||v(x)||2), (19)

which represents the difference in specific orbital energies between the orbits described by some TLE and some RSO’s state . The study of the energy offset (19) on the training data yields the results in Fig. 6.

The energy offset on the training period suggests an underestimation from the collected TLEs w.r.t. to the (supposed) ground truth, that amounts to some , with a deviation of about . Similarly to the study of the angle offset, we can represent this information with a possibility function on as shown in Fig. 7.

Interpreting the information produced by a TLE solely through the orbital plane and the specific orbital energy of the observed RSO, we can then design a possibility function for a “TLE-generator” as

 htle(y|x)=htleang.(δang.(y|x))htleen.(δen.(y|x)). (20)

It is easy to see that in Eq. (20) is indeed a possibility function; in particular, it reaches the maximum value of on the TLE space .

## 5 Scenario implementation

This section details the implementation of the Bayesian estimation scenario on which the fusion of TLEs and radar-like observations will be exploited for the detection and tracking of a single orbiting RSO.

### 5.1 Ground truth and observations

As for the training data on which the analysis of the TLEs was performed in Sec. 4, the reference RSO is Planet Labs’ 0E0E satellite (the corresponding CAT ID assigned by the USStratCom is 41609). The filtering period covers approximately one day, from Sep 6, 2017 00:11:47 (UTC) to Sep 7, 2017 00:00:00 (UTC). A ground truth trajectory was produced for that period from a batch least-square method, exploiting the GPS points provided by Planet Labs.

The TLE points relevant to the filtering period were queried from the USStratCom; their state are given in Tab. 1. A radar with Doppler effect is simulated in order to generate radar observations, following the characteristics of LeoLabs’ PFISR given in Nicolls et al. (2017). Its Field of View (FoV) is set as the ball centered on Fairbanks (Alaska, U.S.A.) with a radius of , and the RSO is observed during eight brief windows over the filtering period (see Fig. 8). The noise covariance matrix of the radar is set as the diagonal matrix

, where the standard deviations are set to the values

, , , .

### 5.2 Filter implementation

The object state space describes the RSO’s position and velocity coordinates in the reference ECI frame. The time flow of the filtering scenario is built as follows: the filtering period is split in even time lapses of , to which the five collection dates of TLEs are added. The resulting time flow is indexed with with

corresponding to the scenario’s initial date. In addition, the duration since epoch J2000 (in

) is given at time by .

#### 5.2.1 Time prediction step

Since the law on describing the RSO’s state is not parameterizable in a straightforward manner, the practical implementation of the Bayesian filtering equations follows a Sequential Monte Carlo (SMC) approach. At time , the posterior pdf is approximated by a set of weighted particles such that

 pk−1(⋅)≃N∑i=1w(i)k−1δx(i)k−1(⋅),

where is the Dirac function at , and with .

The prediction kernel in Eq. (13) aims at describing a Low-Earth Orbit (LEO) trajectory between epochs and , and is constructed as follows. Assuming the object has state at some epoch , the acceleration vector is given by

 ˙v(t)=amod.(p(t),v(t),t),

where the mapping computes the orbital acceleration terms modeled in the scope of this paper. In addition to the central term of Earth’s gravitational pull, it includes the following perturbations: the zonal/tesseral effects up to order and degree , the gravitational pull of the Sun and the Moon, the solar radiation pressure (assuming an Area-to-Mass Ratio (AMR) of , and a spherical shape with a radiation pressure coefficient of ), and a drag term based on the MSIS86 atmospheric model (assuming a ballistic coefficient of ). An additional acceleration term, accounting for the unmodeled perturbations, is built as

 aϵ(p(t),v(t),t,ω(tk−1))=feciric(p,v)((t−tk−1)ω(tk−1)),

where is a zero-mean Gaussian noise with standard deviation on each component in the object’s Radial-Intrack-Crosstrack (RIC) frame, and is the mapping that transforms a vector in the object’s RIC frame to the reference ECI frame. Denoting by

 aorb.(p(t),v(t),t,ω(tk−1))=amod.(p(t),v(t),t)+aϵ(p(t),v(t),t,ω(tk−1))

the total acceleration term describing the orbital dynamics, the time-derivative of the target state is then given by

 [˙p(t),˙v(t)]=[v(t),aorb.(p(t),v(t),t,ω(tk−1))].

Given a particle , the predicted particle is then computed through

 x(i)k|k−1=[p(i)(tk−1),v(i)(tk−1)]+∫tktk−1[v(i)(t),aorb.(p(i)(t),v(i)(t),t,ω(i)(tk−1))]dt. (21)

Since the stochastic component is constant throughout the time interval , Eq. (21) can be solved with a usual numerical integrator (we used Matlab®’s ode45).

#### 5.2.2 Data update step

At time , the predicted pdf is thus approximated by the set . Three cases are now to be considered, depending on the availability of corrective data:

Case 1. No observation is available
This is by far the most frequent case, and also the most straightforward to process. Since no additional information is available on the RSO the prior is set as the posterior , i.e.

 x(i)k=x(i)k|k−1,w(i)k=w(i)k−1,

for any .

Case 2. A TLE is available
This is also straightforward to process. The update equation (14) becomes

 x(i)k=x(i)k|k−1,w(i)k=htle(y|x(i)k|k−1)w(i)k−1∑Nj=1htle(y|x(i)k|k−1)w(i)k−1,

for any .

Case 3. A radar observation is available
The data update equation in this case is similar to the TLE’s above, except that the radar possibility is substituted to the TLE possibility . However, because of the high accuracy of radar observations, implementing the data update mechanism leads to quick degeneracy among the particles, especially when the particle cloud has spread significantly after a long period without observations (Delande et al., 2017). Our approach is to find a parametrization of the predicted pdf leading to a more robust data update mechanism unaffected by the scarcity of particles.

Since pdfs representing orbital states can hardly be parametrized in a simple manner in the Cartesian coordinates , we exploit an alternative space – namely, the spherical coordinates in the sensor’s local frame – in which a Gaussian approximation is more valid. The data update procedure can be summarized as follows (Delande et al., 2017):

1. Transform in the spherical frame : ,

2. Approximate as a Gaussian distribution:

,

3. Update Gaussian distribution, using observation : ,

4. Sample resulting distribution: ,

5. Transform back in the object space : .

Note that the data update step in the spherical frame is a simple linear Kalman update, since the radar is modeled with a Gaussian possibility in Eq. (15).

In any case, the particles are resampled, following the data correction step, if the efficiency ratio falls below a threshold set at .

#### 5.2.3 Initial orbit determination

No prior information is assumed on the RSO’s state prior to its first observation, and an initial orbit determination procedure is implemented in order to determine the posterior distribution following the first radar observation. We followed the admissible region approach developed in DeMars and Jah (2013), where admissible values for the unobserved angular rates are established from the initial radar observation , exploiting internal energy constraints (Tommei et al., 2007; Farnocchia et al., 2010). We employed a similar SMC implementation of the admissible region approach as in our previous works; more details can be found in Delande et al. (2017).

## 6 Filtering results

In order to illustrate the exploitation of TLEs as a complementary data source, we ran the Bayesian tracking filter described in Sec. 5 on two parallel scenarios based on the same ground truth trajectory: one exploits the simulated radar observations only, the other exploits the radar observations and the TLE points collected from the USStratCom’s catalog. The filters’ outputs are averaged over 20 runs, where the observations generated for each run follow the radar model presented in Sec. 5.

We first compare the output of the filter through a Maximum A Posteriori (MAP) estimate of the RSO

’s state. Because of the frequent resampling and the scarcity of data corrective steps, the posterior particle weights are often uniformly distributed; added to the fact that the statistical moments of a

pdf on have little sense in the context of SSA111111For example, even if the particle cloud spreads alongside a ray of close orbital trajectories, the weighted mean of the particle distribution may belong to an entirely different orbit. it is not straightforward to derive a MAP estimate directly from the posterior pdf . Similarly to the processing of radar observations presented in Sec. 5, we use the expression of the posterior pdf in spherical coordinates (still in the reference ECI frame), assumed Gaussian distributed, to derive a MAP estimate of the RSO’s state. The procedure can be summarized as follows:

1. Transform to spherical coordinates: ,

2. Compute the MAP in spherical coordinates as the weighted mean of the resulting distribution: ,

3. Transform back the MAP estimate to Cartesian coordinates: .

Note that the resulting point is not, strictly speaking, the MAP estimate of the RSO’s state, but will be considered as such for the purpose as this analysis. The MAP estimate of the filter’s output is depicted in Fig. 8.

We also compare the output of the filter in the RIC frame of the ground truth trajectory; the result is depicted in Fig. 9 for position coordinates, and in Fig. 10 for velocity coordinates.

By construction, the output of the two filters are nearly identical, from the first pass of the RSO in the radar FoV to the first data correction with a TLE point. Since the radar does not provide information on the RSO’s angular rate, the initial estimation on the RSO’s state is poor and the associated uncertainty grows significantly as soon as the RSO leaves the radar FoV and observations are unavailable. As expected, the following data correction steps with a TLE improve the quality of the MAP estimate (see Fig. 8), and sharpen the distribution significantly (see Figs 9, 10).

Due to the cautious approach in the modeling of the TLEs (see Sec. 4), the information gain of the following TLE points is much less significant, though not inexistent, once the RSO has entered the radar FoV for the third time (from around time step , onwards). This suggests that a more informative TLE model than presented in Sec. 4, built on a larger set of training data and/or on different physical quantities, may benefit the estimation to a larger extent.

## 7 Future works

The introduction of o.p.m.s to the context of SSA opens a wealth of possibilities for the development of Bayesian estimation algorithms, as the modeling of uncertain components in a typical SSA target tracking scenario can be revisited through o.p.m.s (and the associated possibility functions) rather than pdfs.

Many differences between the two approaches are apparent in multi-target detection and tracking scenarios, where marginalized quantities play an important role in the computation of the relative weight (i.e., probability of existence) of individual tracks. For example, because possibility functions are dimensionless, a marginalization term such as the denominator in Eq. (14), assessing the matching between the estimated object and the observation , does not scale whether the distances in the radar observation space are measured in meters or kilometers, the angles in degrees or radians, etc. The usual denominator does, however, since the pdf in Eq. (16) has units and thus scales with the change of reference measure on the observation space. This feature is often overlooked, but raises issues when such a data association event must be compared to another uncertain event that either scales differently (e.g., the association of the same estimated object to a telescope-like observation ), or does not scale at all (e.g., a missed detection event) with arbitrary changes in physical units (Delande et al., 2017).

A natural follow-up on our current work would be to refine the modeling of the TLEs as an uncertain data source. The analysis presented in this paper could be extended to larger period of time and to a wider range of RSOs covering different orbital regimes; studies on the accuracy of TLEs such as in Kelso (2007); Frueh and Schildknecht (2012) would also provide ground for a further refinement of the TLE model. Another natural follow-up would be to propose an o.p.m.-based representation of other data sources already exploited in the context of SSA, such as telescopes, cameras, etc. Less conventional data sources can be found in the numerous natural language statements related to specific launch events, the ill-formatted data collected from unknown/untrusted amateur telescopes, etc. Natural language statements have already been approached with o.p.m.s in the design of target tracking algorithms in Bishop et al. (2018+), and this study could spearhead the integration of a range of data sources relevant to the context of SSA.

Another promising lead to follow is to explore an o.p.m.-based representation of the propagated information, currently represented with a usual pdf in Eqs (13), (14). The initial orbit determination procedure (DeMars and Jah, 2013) leaves no real knowledge of the distribution of the RSO’s state within the admissible region: this is a straightforward example where the information possessed by the analyst is represented in a more natural and less prejudiced manner with a possibility function equal to one within the admissible region and zero outside, rather than a uniform pdf supported by the admissible region.

Orbital propagators designed for practical SSA tracking algorithms always operate on limited information. They approximate physical effects affecting the orbital trajectory (e.g. limited order in the zonal/tesseral components of the Earth’s gravitational pull), discard them on purpose (e.g. the Earth’s shadow is ignored) or, less conspicuously, ignore perturbation effects whose existence is simply unknown to the space analyst. Given the complexity of the orbital dynamics and the many sources of uncertainty shaping the design of an orbital propagator, a probabilistic representation of the prediction model appears over-descriptive and could benefit from an o.p.m.-representation as well.

## 8 Conclusion

In this paper, we introduce a new representation of uncertainty for Bayesian detection and tracking algorithms in the context of SSA, based on outer probability measures rather than probability density functions. Less descriptive than pdfs, o.p.m.s do not characterize the distribution of an uncertain system’s state; rather, they aim at proposing a prejudice-free representation of the system’s state, matching the information possessed by some operator – say, a space analyst – about that system and nothing more. We show that each uncertain component of a single object Bayesian estimation problem – the estimated state, the prediction model, the data source(s) generating observations/opinions – can be represented either with an o.p.m. or a pdf in a single coherent Bayesian estimation framework.

More versatile than pdfs, o.p.m.s allow for the modeling of uncertain components for which statistical information is scarce or inexistent. In particular, we show that TLEs can be treated as a data source in a Bayesian estimation problem through a simple o.p.m.-based model, and fused with radar observations in a RSO tracking algorithm. This concept is then illustrated on a scenario where TLE points queried from the USStratCom catalog and simulated radar observations are exploited in order to track a LEO satellite.

## Acknowledgements

The authors wish to thank Vivek Vittaldev (Planet Labs) for his help in accessing and intepreting GPS data collected by Planet Labs, and Marcus Bever (University of Texas at Austin) for the extraction of USStratCom’s TLEs and Planet Labs’ GPS points, and their formatting into a common reference frame.

## References

• Bishop et al. (2018+)

Bishop, A. N., Houssineau, J., Angley, D., Ristić, B., 2018+. Spatio-temporal tracking from natural language statements using outer probability theory.

• Blackman (2004) Blackman, S. S., Jan. 2004. Multiple Hypothesis Tracking for Multiple Target Tracking. Aerospace and Electronic Systems Magazine, IEEE, 5–18.
• Delande et al. (2018) Delande, E. D., Frueh, C., Franco, J., Houssineau, J., Clark, D. E., 2018. Novel Multi-Object Filtering Approach For Space Situational Awareness. Journal of Guidance, Control, and Dynamics 41 (1), 59–73.
• Delande et al. (2016) Delande, E. D., Houssineau, J., Clark, D. E., 2016. Multi-object filtering with stochastic populations. ArXiv preprint, arXiv:1501.04671v2.
• Delande et al. (2017) Delande, E. D., Houssineau, J., Franco, J., Frueh, C., Clark, D. E., Feb. 2017. A new multi-target tracking algorithm for a large number of orbiting objects. In: 2017 AAS/AIAA Spaceflight Mechanics Meeting.
• DeMars and Jah (2013)

DeMars, K. J., Jah, M. K., 2013. Probabilistic Initial Orbit Determination Using Gaussian Mixture Models. Journal of Guidance, Control, and Dynamics 36 (5), 1324–1335.

• DeMars et al. (2012) DeMars, K. J., Jah, M. K., Schuhmacher, D., 2012. Initial orbit determination using short-arc angle and angle rate data. Aerospace and Electronic Systems, IEEE Transactions on 48 (3), 2628–2637.
• Dempster (1967) Dempster, A. P., 1967. Upper and lower probability inferences based on a sample from a finite univariate population. Biometrika 54 (3-4), 515–528.
• Dubois and Prade (1983) Dubois, D., Prade, H., 1983. Ranking fuzzy numbers in the setting of possibility theory. Information sciences 30 (3), 183–224.
• Farnocchia et al. (2010) Farnocchia, D., Tommei, G., Milani, A., Rossi, A., 2010. Innovative methods of correlation and orbit determination for space debris. Celestial Mechanics and Dynamical Astronomy 107 (1-2), 169–185.
• Friedman and Halpern (2001) Friedman, N., Halpern, J. Y., 2001. Plausibility measures and default reasoning. Journal of the ACM 48 (4), 648–685.
• Frueh and Schildknecht (2012) Frueh, C., Schildknecht, T., 2012. Accuracy for Two-Line-Element Data for Geostationary and High-Eccentricity Orbits. Journal of Guidance, Control, and Dynamics 35 (5).
• Houssineau (2015) Houssineau, J., 2015. Representation and estimation of stochastic populations. Ph.D. thesis, Heriot Watt University.
• Houssineau (2018) Houssineau, J., 2018. Detection and estimation of partially-observed dynamical systems: an outer-measure approach. arXiv:1801:00571.
• Houssineau and Bishop (2017) Houssineau, J., Bishop, A., 2017. Smoothing and filtering with a class of outer measures. ArXiv preprint arXiv:1704.01233.
• Houssineau and Clark (2016) Houssineau, J., Clark, D. E., 2016. Multi-target filtering with linearised complexity. arXiv:1404.7408v2.
• Jones and Vo (2015) Jones, B. A., Vo, B.-N., Jan. 2015. A Labelled Multi-Bernoulli Filter for Space Object Tracking. In: 2015 AAS/AIAA Spaceflight Mechanics Meeting. pp. AAS 15–413.
• Jones et al. (2016) Jones, B. A., Vo, B.-T., Vo, B.-N., Sep. 2016. Generalized Labeled Multi-Bernoulli Space-Object Tracking with Joint Prediction and Update. In: AAS/AIAA Astrodynamics Specialist Conference.
• Kelso (2007) Kelso, T. S., Jan. 2007. Validation of SGP4 and IS-GPS-200D Against GPS Precision Ephemerides. In: 2007 AAS/AIAA Space Flight Mechanics Conference.
• Nicolls et al. (2017) Nicolls, M., Vittaldev, V., Ceperley, D., Creus-Costa, J., Foster, C., Griffith, N., Lu, E., Mason, J., Park, I.-K., Rosner, C., Stepan, L., Sep. 2017. Conjunction Assessment For Commercial Satelitte Constellations Using Commercial Radar Data Sources. In: Advanced Maui Optical and Space Surveillance Technologies Conference.
• Shafer (1976) Shafer, G., 1976. A mathematical theory of evidence. Vol. 1. Princeton university press Princeton.
• Tommei et al. (2007) Tommei, G., Milani, A., Rossi, A., 2007. Orbit determination of space debris: admissible regions. Celestial Mechanics and Dynamical Astronomy 97 (4), 289–304.
• Vo and Vo (2013) Vo, B.-T., Vo, B.-N., 2013. Labeled Random Finite Sets and Multi-Object Conjuguate Priors. Signal Processing, IEEE Transactions on 61 (13), 3460–3475.
• Walley (1991) Walley, P., 1991. Statistical reasoning with imprecise probabilities.
• Worthy III and Holzinger (2017a) Worthy III, J. L., Holzinger, M. J., Sep. 2017a. Uninformative Prior Multiple Target Tracking Using Evidential Particle Filters. In: Advanced Maui Optical and Space Surveillance Technologies Conference.
• Worthy III and Holzinger (2017b) Worthy III, J. L., Holzinger, M. J., 2017b. Use of Uninformative Priors to Initialize State Estimation for Dynamical Systems. Advances in Space Research 60 (7), 1373–1388.
• Yen (1990) Yen, J., 1990. Generalizing the Dempster-Schafer theory to fuzzy sets. IEEE Transactions on Systems, Man and Cybernetics 20 (3), 559–570.
• Zadeh (1965) Zadeh, L. A., 1965. Fuzzy Sets. Information and Control 8, 338–353.