Improving Orbital Uncertainty Realism Through Covariance Determination in GEO

The reliability of the uncertainty characterization, also known as uncertainty realism, is of the uttermost importance for Space Situational Awareness (SSA) services. Among the different sources of uncertainty related to the orbits of Resident Space Objects (RSOs), the uncertainty of dynamic models is one of the most relevant ones, although it is not always included in orbit determination processes. A classical approach to account for these sources of uncertainty is the consider parameters theory, which consists in including parameters in the underlying dynamical models whose variance aims to represent the uncertainty of the system. However, realistic variances of these consider parameters are not known a-priori. This work presents a method to infer the variance of the consider parameters, based on the distribution of the Mahalanobis distance of the orbital differences between predicted and estimated orbits, which theoretically shall follow a χ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2$$\end{document} distribution under Gaussian assumption. This paper presents results in a simulated scenario focusing on Geostationary (GEO) regimes. The effectiveness and traceability of the uncertainty sources is assessed via covariance realism metrics.


Introduction
Orbital uncertainty plays a key role for the provision of Space Situational Awareness (SSA) services, including catalogue maintenance, risk assessment or maneuver detection, among others. Therefore its adequate characterization, also defined This article belongs to the Topical Collection: Advanced Maui Optical and Space Surveillance Technologies (AMOS 2021) Guest Editors: Lauchie Scott, Ryan Coder, Paul Kervin, Bobby Hunt.

3
The Journal of the Astronautical Sciences (2022) 69:1394-1420 as uncertainty realism, is of paramount importance. Uncertainty realism focuses on correctly representing the Probability Density Function (PDF) of the orbital state. Under Gaussian assumption, uncertainty realism can be reduced to covariance realism, requiring only the two first moments of the PDF for the proper characterization of the system uncertainty, gathered in a state and its associated covariance matrix.
Many Orbit Determination (OD) processes in Space Surveillance and Tracking (SST) are based on weighted batch least-squares algorithms, which rely on available and sufficient measurements to produce an estimate (orbital state and covariance matrix). During this process, the dynamics governing the system is usually assumed deterministic, with measurement noise as the only source of uncertainty. Thus, the obtained covariance matrix is able to account only for the measurements noise, being known as the noise-only covariance [18]. However, one of the main sources of uncertainty during OD and subsequent propagation arises from the lack of knowledge in the underlying dynamical models, which is typically disregarded. The impact of this uncertainty is crucial not only for the state estimation but also for the time evolution of its associated uncertainty due to its inherent correlation with the state variables. This leads to an overly-optimistic noise-only covariance matrix time evolution and, eventually, the loss of covariance realism.
Therefore, it is customary for SSA and particularly for SST purposes to characterize and determine such uncertainty and its effects, which is commonly known as Uncertainty Quantification (UQ). Two fundamental problems can be distinguished for UQ: the propagation of uncertainty and the inverse problem (model and parameter uncertainty) [19]. The former is focused on how to propagate forward an initially given PDF of a state, accurately and efficiently. The inverse problem, on the contrary, consists in assessing the differences between the observed behavior of a system and the underlying models and parameters used to represent it. Regarding the uncertainty in the modelling, a common approach is to drop the deterministic assumption of the equations of motion, recurring to stochastic noise models such as Brownian motion, Ornstein-Uhlenbeck or Gauss-Markov processes [6,19,21]. The other target of the inverse problem is the parameter uncertainty, whose objective is to represent the uncertainty in specific terms of the dynamic or measurement equations. This is the core of the work at hand, specifically focusing on uncertain parameters that can be estimated during the OD process as dynamic parameters, not only for the quantification of its uncertainty but also to represent the relationship between the uncertain parameters and the state variance.
There are other techniques conceived to improve covariance realism without focusing on the sources of uncertainty and their modelling. For instance, Gaussian and linear assumptions are held longer in time when the state is represented in mean orbital elements, allowing the covariance matrix to represent more realistically the uncertainty of Monte Carlo simulations [11,26]. Other typical representations of the state and covariance in non-linear reference frames that are able to slow down the realism degradation upon propagation are being widely studied, such as the QtW frame in [12]. However, despite the advantage of representing the state of the object in orbital elements for covariance realism, we concern ourselves to Cartesian state representation. Even though the proposed methodology is agnostic to the state representation, the operational goal of this work leads us to remain using Cartesian coordinates as is customary in most SST operational scenarios. Other approaches suggest the use of empirical covariance matrices to include all residuals of the estimation process in the covariance computation, regardless of whether the uncertainty has been specifically modelled or not [2,3,8]. This proposal claims to account more accurately for noise time-variations rather than process noise or consider parameter analysis, at the expense of the physical interpretation of the uncertainty.
In an operational environment, simple and robust techniques are required to improve covariance realism since, as previously discussed, the nominal covariance determination methods tend to provide optimistic results. The most applied solutions are: • Kalman filters with process noise matrix to introduce the model uncertainty into the system. These UQ methods are gaining acceptance over stochastic acceleration methods in the current state of the art since they can account for both dynamic model and parameter uncertainty. However, a physically-based derivation of a process noise is rather challenging [22]. Other authors suggest calibration methods to estimate such process noise [4], but typical solutions lack the physical meaning of the applied correction and are not suitable for batch processing, the common framework of SST. • Scaling techniques which inflate the covariance by certain factors. Some authors propose the computation of such scaling based on increasing the initial position uncertainty to match the velocity error [7] whereas other options explore the usage of the Mahalanobis distance of the orbital differences to find the scale factor [13]. However, a common drawback of artificially increasing the covariance is that the physical meaning of the correction is lost, not being able to understand the contributions of each source of uncertainty. These sort of methods are used nowadays in operation centers such as Space Operations Center (CSpOC) [19].
An additional option to be discussed is the consider analysis theory, which is a classical approach for parameter uncertainty analysis in the equations of motion for OD processes [18]. It consists in extending the state space by including parameters in the dynamic models, such as atmospheric force, solar radiation pressure force or measurement models. These parameters are assumed to follow a Gaussian distribution with a null mean and a certain variance. A null mean allows to maintain an unbiased estimation, whereas the uncertainty of the parameter is accounted for in the state covariance matrix. The consider parameter theory is thus designed to tackle the parameter uncertainty of the previously described UQ inverse problem, inflating the estimated covariance to account for such parameter uncertainty. This theory is compatible with both batch and filter applications such as in the Schmidt-Kalman filter [10,28] and provides the clear advantage of assessing the effect of specific uncertain parameters in the process, maintaining the physical meaning of the covariance correction as opposed to scaling factors techniques. However, one of the main drawbacks of the consider parameter theory is that realistic variances of such parameters are not normally known. This can lead either to overly-sized or underestimated state covariance matrices, failing to model the uncertainty of the system and thus, not achieving covariance realism. This work presents a novel methodology to determine the variance of the consider parameters included in the dynamical model of an RSO. It is based on the orbital differences (position, velocity, or both) between estimated and predicted orbits projected into the covariance space (i.e. Mahalanobis distance), which under Gaussian assumptions shall follow a 2 distribution to achieve covariance realism. In other words, if the orbital differences are normally distributed and correctly represented by a covariance matrix, the Mahalanobis distance distribution is equivalent to the squared sum of a normal distributions from each component of the orbital differences (i.e. a 2 distribution), provided that no correlation exists between them. Therefore, a minimization process can be designed to obtain the consider parameter variances that provide the best match between the observed Mahalanobis distance distribution and the expected 2 one.
A similar analysis based on the consider parameter theory to improve the covariance realism is performed in [15], a precursor work for this study. There, it is proposed to correct the noise-only covariance with a least squares fitting to a so-called observed covariance, this latter being obtained from Monte Carlo orbital differences aggregation. This approach has a main drawback, which is that to compute such observed covariance, orbital differences at distinct orbital positions are mixed from orbit estimates based on different observation scenarios. This issue is mitigated by the normalization obtained with the Mahalanobis distance concept, which is the cornerstone of the methodology presented in this paper.
In the work at hand, the covariance determination methodology is applied to GEO RSOs, continuing the efforts of previous studies that applied the proposed methodology to Low Earth Orbit (LEO) RSOs for drag and range bias uncertainty [1]. The realism of the determined covariance matrices remains as the cornerstone of this study, and thus specific covariance realism metrics such as the covariance containment are analyzed. For a geostationary orbit, two of the main sources of uncertainty that come into play are related to the Solar Radiation Pressure (SRP) and the time bias of the sensors, which are analysed in the presented work.
The paper is structured as follows: in Sect. 2, the consider parameter theory in the context of batch least-squares estimation is revisited, as well as the linear covariance propagation. In Sect. 3, the proposed models for the consider parameters and the covariance determination methodology is described, together with the validation approach. In Sect. 4, the results of the simulation campaigns carried out for validation are presented and discussed. Finally Sect. 5 contains the conclusions and future steps of this study.

3 2 Background
This section briefly describes the consider parameter theory in the context of the common batch least-squares processes in SST, as well as the linear propagation of the covariance.

Consider Parameter Theory in Batch Least-Squares Algorithm for Orbit Determination
The complete description of the consider parameter theory (or consider covariance analysis, as termed by some authors) can be found in [18,24]. For brevity, only the final derivation in the nominal batch least-squares process is described next. Let us define the estimated state vector as where (t) , (t) and n y are the position, velocity and estimated state dimension, respectively. (t) represents the estimated parameters, either applied to the force or the measurement models. Typical examples of these parameters are the drag coefficient ( C D ) in LEO or the solar radiation pressure coefficient ( C R ) in GEO. Following the nominal OD process [18], the noise-only covariance is where y corresponds to the Jacobian of the observations with respect to the estimated state, and is the weighting matrix containing the confidence of each measurement and the possible correlation among the measurements. The consider parameters to be modelled in our system can be gathered in the consider parameter vector where n c is the number of consider parameters. They are defined to follow a normal distribution, i.e, This definition allows the expected value of the orbit estimation to remain unbiased provided that the consider parameters have null mean and their variances are uncorrelated with the measurement noise. On the contrary, the covariance of the estimation is affected. The consider covariance, which is defined as the covariance of a where n c and n y are the number of consider parameters and the state vector dimension, respectively. Therefore, the consider covariance is obtained as the noise-only covariance plus a covariance correction, which depends linearly on the consider parameter variances.

Linear Covariance Propagation
The consider covariance is obtained at estimation epoch and models the effect of the uncertainty of the consider parameters that affect the estimation, regardless whether their uncertainty affects the force or the measurements model. However, our interest in enhancing the covariance realism extends also to the covariance propagation for SST purposes. In this paper, we assume linear propagation of the covariance matrix, as it is generally applied in most operational scenarios in SST. More complex and accurate uncertainty propagation methods are out of the scope of this work since Gaussianity is a cornerstone assumption in the proposed methodology. In this respect, Michael's normality tests can be applied to assess data linearity [17,20]. A complete derivation of linear propagation theory can be found in many well-known references, such as [18]. In the end, to account for the effect of the main dynamic parameters in the propagation of the state, it is required to integrate the variational equations, whose solution is the Extended State Transition Matrix (ESTM) where ∈ ℝ n y ×n y , ∈ ℝ 6×6 , ∈ ℝ 6×n p , ∈ ℝ n p ×n p and • n p is the number of dynamical parameters to consider during propagation, which in this case corresponds to the estimated dynamical parameters, excluding position and velocity, hence, n p = n y − 6.
• t, t 0 corresponds to the state transition matrix, which relates the position and velocity at any time t with respect to the initial state at time t 0 .
• t, t 0 is the so-called sensitivity matrix, which contains the partial derivatives of the state vector with respect to the model dynamical parameters, which are assumed to be constant.
The ESTM can be computed by solving numerically its associated partial differential equations as shown in [18]. Linear covariance propagation is then carried out as By using the ESTM, the effect of the uncertainty of the dynamical model parameters is mapped into the position and velocity covariance not only at estimation, but also along the propagation.

Methodology
In this section, we define first the consider parameters proposed to model the uncertainty in GEO. Then, we describe the concept of the Mahalanobis distance and its connection to the 2 probability distribution and develop the proposed procedure to effectively calculate the consider parameter variances that define the uncertainty levels of the dynamic model and compute the consider covariance. Finally, the validation procedure is described.

Consider Parameter Models
As we mentioned previously, two consider parameters of interest are included in our analysis: in the SRP force and the time bias of the sensor. This subsection describes their proposed models.

Solar Radiation Pressure
Following the classical definition [18], the SRP acceleration equation with the first consider parameter reads where P SRP is the solar radiation pressure, C R is the solar radiation coefficient, A is the cross-sectional area, m is the mass of the object, r is the distance to the Sun, AU expresses the magnitude of an astronomical unit and c SRP is the consider parameter associated to the SRP, which is devised to follow a Gaussian distribution as defined (10) in Eq. (4) with standard deviation SRP . The aim of this parameter is to model the uncertainty that can be associated to most of the components of Eq. (11). The mass may vary for maneuverable satellites, the cross-section area is assumed constant along the estimation and propagation, which is not true in general; the solar radiation pressure and the solar radiation coefficient are affected by the Sun's behavior (solar cycles) and satellite surface variability (light reflection and absorption), which are not modelled for most SST applications.

Sensor Clock Time Bias
The second consider parameter aims to represent the variability in the time bias present in the sensors when time-stamping the measurements. It is characterized as where t denotes time and c TB is the consider parameter associated to the time bias, also following the definition of Eq. (4) with standard deviation TB . This consider parameter is affecting the measurements model and is associated to each observation.

Covariance Determination Method
This subsection describes the combination of the consider parameter theory with the concept of Mahalanobis distance and defines the proposed methodology to determine the consider covariance.

Mahalanobis Distance and Consider Covariance Determination
The Mahalanobis distance is a well-known statistical metric that describes how far a state (t) is from a certain reference ref (t) , projected into the covariance space [16]. This is: where and ref are the covariance of both the state and the reference, respectively. However, the covariance of the reference is generally several orders of magnitude smaller, which enables us to neglect it many applications. For the sake of clarity, the reference covariance is omitted in the following equations, although it is considered in the computations when applicable. In order to introduce the consider parameter effect, we recall the definition of the state vector of Eq. (1) and combine Eq. (10) with Eq. (13), which yields Thus, Eq. (14) allows the computation of the squared Mahalanobis distance at any propagation epoch as a function of the consider parameter variances included inside (12)

3
matrix . It is worth mentioning that the projection of the orbital differences in the covariance space allows to combine samples derived from ODs at different epochs and conditions, which is of paramount importance for an operational scenario.

Consider Parameter Variance Determination
It remains to describe the method to compute the variance of the proposed consider parameters that is applied in this work following the previous line of research of Cano et al. [1]. Under linear and Gaussian assumptions, this is, when the differences between the state and the reference are normally distributed and the covariance is representative of such distribution (i.e. realistic), a population of squared Mahalanobis distances should follow a 2 distribution, whose detailed characteristics may be found in D'Agostino and Stephens [5]. Additionally, the 2 distribution requires the aggregated normal distributions to be uncorrelated. Consequently, to reduce the correlation between the state variables, the Mahalanobis distance is computed in the Intrack, normal , and cross-track (TNW) local frame [25]. The in-track axis is defined as tangential to the orbit and parallel to the velocity. The normal axis is placed in the orbit plane, perpendicular to the in-track axis, and the cross-track axis is perpendicular to the orbit plane. Using local frames aligned with the satellite motion allows to reduce the correlation of orbital differences and helps to trace the impact of the uncertainty sources modeled in the analysis. Equation (14) allows us to compute such Mahalanobis distance at any desired epoch during the propagation arc by comparing an estimated and later propagated (predicted) orbit against a reference. Therefore, if a sufficient number of predicted orbits is available, it is possible to look for the variance of the consider parameters represented in matrix so that the squared Mahalanobis distance population resembles the theoretical 2 distribution. In the work presented here, this optimization problem is based on minimizing the differences between two Cumulative Distribution Functions (CDFs), the observed one and the theoretical 2 one. In this case, the theoretical 2 distribution should have as many Degrees Of Freedom (DOF) as components of the orbital differences are included in the Mahalanobis distance. This minimisation process can be adapted to any desired number of consider parameters, being able to combine the two proposed for our analysis (SRP and time bias). The cost function to minimize in order to determine the consider parameter variances is defined next for completion. It consists in grouping the observed squared Mahalanobis distances in N b bins, comparing their joined probability against the theoretical 2 one as where CDF E corresponds to the empirical CDF of the ith bin, and CDF 2 the CDF of the 2 distribution. Considering more statistically robust CDF comparison metrics was out of the scope of the present work. However, the application of other metrics such as Cramer-von-Mises or Kolmogorov-Smirnov [5] is an open line of research.

Validation Chain
In this section, we describe a process to generate the required data to apply the proposed covariance determination method in a simulated scenario for validation purposes.
The presented methodology can be validated by checking if the obtained consider parameter variances are representative of the uncertainty present in the system. To confirm this, we run a Monte Carlo iterative scheme in which we simulate a population of orbits. Under simulation conditions, we are able to impose the uncertainty levels in the time bias and SRP. Thus, we condition the outputs of the optimization process to known, controlled values for validation purposes. The steps of the data generation are described below. For further clarification of this scheme, the different kinds of orbits in a relative propagation timeline is shown in Fig. 1 exemplifying a single Monte Carlo iteration.
• From a GEO reference state associated to an existing RSO and a certain epoch ( t 0 , estimation epoch), we perform an orbit propagation process to obtain the reference orbit, which is used to obtain the orbital differences required for the Mahalanobis distance in Eq. (14). The details of the RSO can be found in Table 1. The dynamic model of this propagation is deterministic, without including any additional perturbation in the SRP or time bias models. See Table 2 for model details on the dynamic model characteristics. • In a Monte Carlo process, we generate a population of orbits sampling the Gaussian distribution of the consider parameters models, and introducing each sample as a constant perturbation for each orbit. We then perform the following steps per iteration, also depicted in Fig. 1:   Fig. 1 Monte Carlo scheme showing the process to generate the predicted and reference orbits, the typical propagation arc and analysis epochs, and the operational reference orbit 1. We propagate the reference state backwards 28 days considering a SRP consider parameter perturbation sample in the dynamics. We refer to this propagated orbit as the simulated orbit. 2. We generate tracks of the simulated orbits using a model of a ground-based telescope, whose details can be found in Table 3. We introduce the time bias uncertainty on each simulated observation, again selecting a perturbation sample from the consider parameter model defined in Sect. 3.1.2. Along the generation of these simulated observations, the measurement noise is also included. 3. We perform an OD process based on the simulated tracks. We have inputs from a period spanning between t 0 -28 days and t 0 , where t 0 corresponds to the time of estimation. The estimation output (state and noise-only covariance)   Table 2. This orbit is known as the predicted orbit.
• From one iteration to the next, the simulation time frame is shifted one day forward, so that the generation of observations in the OD process is aligned in time as in an operational scenario. 400 samples, which correspond to more than 1 year of daily tracking data for an object, are expected to suffice for the statistical computations of the proposed methodology, as was found in the simulations performed in previous studies [1]. • Operational reference orbit: We propose an alternative orbit to use for the computation of the orbital differences of Eq. (14) other than the aforementioned reference orbit. The operational reference orbit consists in the output of an OD whose determination arc includes the propagation epochs of interest for each Monte Carlo iteration. In other words, if we want to analyse the current orbit between t 0 +7 days and t 0 +21 days, we can use as reference the output of another OD whose arc ranges from t 0 days up to t 0 +28 (see Fig. 1). The reason behind this choice is that, for most SST operational environments, precise orbits whose ephemeris can be assumed perfect are not commonly accessible. For instance, external sources of precise GNSS orbit data is not available when the targets are non-collaborative objects such as space debris. It is important to remark that, the estimation corresponding to this operational reference orbit has certain noiseonly covariance due to the measurements uncertainty. Thus, it is necessary to include its noise-only covariance inside the Mahalanobis distance computation of Eq. (14).
The Mahalanobis distances can be calculated once all the necessary data for all the orbits is available at the desired analysis epochs. Then, the proposed covariance determination methodology can be applied. This process is summarised in Algorithm 1.

Results and Discussion
This section is divided in three different subsections. The first one describes the characteristics of the simulated environment, providing details on the dynamical model used and the expected uncertainty levels. The second section describes the results of the validation chain of the proposed covariance determination methodology in several test scenarios. Finally, the effectiveness that the proposed methodology achieves in terms of covariance realism is further analysed.

Simulation Environment
For the sake of completeness, this section characterizes the conditions under which the validation test sequence has been performed. Table 1 contains the information corresponding to the reference RSO. Table 2 defines the applied dynamic model. Finally, Table 3 describes the parameters of the simulated telescope used for the generation of orbit tracks (set of measurements).
The tracking strategy, using the telescope described in Table 3, consists on series of three measuring intervals per night. Every night, the telescope takes measurements at intervals of 15 min located at the beginning, middle and end of the night, with a time-step of 20 s. The exact position of the intervals within the night may vary slightly due to shadow conditions. This approach is taken to resemble telescope availability in an operational scenario.
It remains to describe the chosen uncertainty variances to introduce during the validation campaign. To this end, it is customary to resort to the literature for expected uncertainty levels in operational environments for the target uncertainty sources that have been modelled as consider parameters. Regarding the SRP uncertainty, according to Montenbruck and Gill [18] an annual variation of 3.3% is estimated. Other studies performed for Ajisai and LAGEOS geodetic satellites [9] identified variations ranging between 1% and 3%. Thus, for the purposes of the work at hand, a 5% standard deviation in SRP has been considered a realistic uncertainty for validation. Nonetheless, higher uncertainty levels are also analysed in Sect. 4 for validation purposes. The characterization of the time bias behavior is more complex, mostly due to its strong dependence on the type of instrument. From Lee and Pogge [14] and teindorfer, Kirchner et al. [23], typical uncertainty levels ranging from 1 to more than 100 ms can be extracted. Similarly to the SRP case, 100 ms of standard deviation has been selected as a realistic uncertainty for validation, despite reaching higher values in the presented results.

Simulation Results
The main objective of the validation phase is to verify that the uncertainty introduced in the Monte Carlo simulation is correctly retrieved as consider parameter variances in the optimization process. To verify this, we propose a sequence of validation tests, summarized in Table 4. Only the most relevant test cases are included.
The tests become gradually more complex. The first three series (I to III) represent the core validation pipeline, in which we subject the dynamic model to perturbations significantly higher than the ones suggested by the literature (Sect. 4.1) to verify the robustness of the presented covariance determination method in highlyperturbed environments where data Gaussianity is stressed to the limits of the methodology assumptions. In the first series, we restrain our analysis to the estimation epoch t 0 , starting from null perturbation. We then proceed to include perturbations as uncertainty in the system. In the second series, we extend the analysis to a future epoch, located at 15 days of propagation into the future. The third series is representative of the final complexity of the methodology, combining Mahalanobis distance calculations from epochs located all over an interval of 14 days, with a time-step of one day. Including different epochs of the propagation arc in the same optimization process has the objective of computing a single consider parameter variance that is able to improve the realism of the covariance in the complete interval of interest, which is a desired operational feature. The fourth series is conceptually equal to the third, but show a case with reduced perturbation levels in line with the discussion Sect. 4.1, looking for values of uncertainty closer to the ones expected in the literature.
The analyses were performed discarding velocity components (in the minimization phase), due to the accumulation of small errors in the in-track velocity after propagation. Apart from the velocity components being several orders of magnitude smaller than the position ones, the expected precision of the in-track velocity was found to be the smallest covariance term among the velocity components. This lead to ill-conditioning of the covariance and an exponential increase in the Mahalanobis distance for high perturbations after more than 3 days of propagation. Therefore, only 4 DOF remain for the 2 distribution, namely the position differences (TNW frame) and the SRP coefficient. Nonetheless, since a series of different position ephemeris are considered, the dynamics of the system are fully regarded in the computations.
The overall results of the validation test sequence, paying special attention to the last ones due to their increased level of complexity, has shown to provide satisfactory results when computing the consider parameter variances. The deviations between the input perturbations and output consider parameter variances did not exceed 11% for any of the perturbations in any of the cases.
The starting point of the validation (Case I-A) explores the ideal situation in which no perturbation is present in the dynamic model. Then, if the measurement noise is properly characterized, the noise-only covariance is expected to be representative of the uncertainty of the system, and the squared Mahalanobis distance distribution should directly resemble a 2 one. Figure 2 shows the results of this case, where this fact can be appreciated. The retrieved distribution is very close to the 2 distribution without applying any consider parameter correction. This is a clear indication that in the sole presence of measurement noise, the noise-only covariance is able to represent the uncertainty of the system. Moreover, It can be noted that the number of DOF in this case is 7. As opposed to the rest of the tests, series I cases are analysed at t 0 , not suffering from in-track velocity error accumulation during propagation.
On the contrary, the noise-only covariance no longer represents the uncertainty of the system when dynamic model uncertainty is introduced. Case I-B shows the results of including both consider parameter perturbations simultaneously, performing two-variable optimization at the estimation epoch. As seen in Table 4 the methodology is capable of computing consider parameter variances within a margin lower than a 10%.
However, we discussed that from an operational perspective it is desirable to extend this analysis into the prediction region. In the second series (II-A, II-B), we establish the analysis epoch at t 0 +15, proving that the quality of the results obtained in the previous series can be maintained when propagating up to this time horizon. Additionally, we introduce the usage of operational reference 1 3 The Journal of the Astronautical Sciences (2022) 69:1394-1420 orbits in the case II-B. Thus, we include a new layer of realism by not relying on an absolutely true representation of the reference orbit, but rather on the result of an OD process with an associated covariance matrix. As was observed in [1], obtaining appropriate results using the operational reference orbit is of high relevance in operational scenarios, allowing the methodology to function using only the observations, not requiring other external information sources.
The third series of cases extend the analysis epoch from a single one to a full interval spanning for 14 days between t 0 +7 and t 0 +21. Table 4 shows the most relevant subcases, combining both consider parameters and comparing against reference and operational orbits. Mahalanobis distances are calculated at specific points within this period, at time steps of one day. The main leap forward in this series is that a singular optimized consider parameter variance is able to recover the 2 behavior of the squared Mahalanobis distance population for a wide propagation interval, obtaining consider parameter variances very similar to the introduced perturbations as can be seen in Table 4 and recovering the 2 behaviour as can be observed in Fig. 3. Again, accurate results are obtained independently on the orbit used to obtain the orbital differences. Besides the fact that adding  Figure 4 corresponds to the same case scenario as in Fig. 3, but in this case maintaining the noise-only covariance for the Mahalanobis distances computations. We find that the values of the distribution are far from resembling the 2 behavior. This figure has been included to show the inability of the noise-only covariance to represent the uncertainty of the system when model uncertainty is present. The high values of the distribution indicates that, without the consider parameter variances correction, the covariance is overly-optimistic and the observed orbital differences are much larger than the standard deviation present in the covariance.
The fourth series of tests does not introduce new elements into the methodology. It serves to prove that the model is sensitive to lower perturbation levels which are aligned with those suggested by the literature. Figure 5 contains the fitting results of the last test case, IV-A, which is representative of the full complexity of the methodology proposed in this work. Again, the differences between the input perturbation and the results of the consider parameter variances is lower than a 10%, showing that the methodology is able to recover the theoretical 2 behavior and thus improve the covariance realism.

Containment Analysis
To obtain a physical interpretation and a proper visual representation of the effectiveness of the proposed covariance determination method for improving the covariance realism, covariance containment tests, such as the one proposed in [27], are provided in this section. To evaluate whether the covariance is representative of the orbital differences (i.e. realistic), the Mahalanobis distance can be used as a metric to see the number of points that lay inside a k − ellipsoid ( k = 1, 2, 3, 4 ) and compare it against the theoretical expected fraction for a multivariate Gaussian distribution of as many Degrees Of Freedom (DOF) as components of the state vector used for the Mahalanobis distance. In all tests presented previously, only position and the SRP coefficient of the state vector were included in the Mahalanobis distance computation. Gaussianity tests of the orbital  Table 5 summarises the evolution of the containment metric at different epochs, comparing both noise-only and consider covariance matrices versus the 4 DOF theoretical containment results. The results found in Table 5 were obtained using the optimum consider parameter variances found in the results of case III-B. The results of this test case are chosen for this analysis due to a high resemblance to an operational scenario, as it maintains a realistic amount of data, it uses the operational reference orbit for comparison and achieves sufficient accuracy.
It is clear from Table 5 that the proposed methodology leads to an improvement in covariance containment when model uncertainty is present. The results at all epochs are close to the theoretically expected ones, particularly for 3 and 4 . There are several implications of these results. Firstly, the average containment along the complete propagation interval resembles closely its theoretical expectation. Recalling the operational goal of the methodology, this shows that a unique consider parameter variance is able to improve substantially the covariance realism in the interval of interest for GEO scenarios. Secondly, the degree of similarity of the containment in all sigma ellipsoids indicates that the consider covariance is not over-sized, showing that the proposed methodology is able to tackle the model uncertainty properly and maintain the traceability of its sources. On the contrary, it is directly observed the lack of covariance realism in the absence of consider parameter correction when model errors are present. Noiseonly covariance matrices are too optimistic and fail to represent the PDF of the state. Another way to visualize this phenomenon can be found in Figs. 6 and 7, where we compare the amount of orbital differences included in the 3 ellipsoid (in green) in both cases: with consider parameter and without consider parameter, respectively. The orbital differences are shown in the relative frame TNW (intrack, normal and cross-track directions) Figure 6 shows that the orbital differences are larger in the in-track direction. The consider covariance is elongated in such direction, representing the actual distribution. The ellipsoidal shape is not easily discerned in Fig. 6 due to the significantly larger standard deviation of the in-track component. It is clearly observed that the consider parameter methodology correction with our determined variance elongates the ellipsoid in the in-track direction, which is the one showing a higher dispersion.
For further analysis, Figs. 9, 8 and10 are included to show the effect on the position covariance that enables the consider parameter correction, focusing on a single orbit and separating the contribution of each of the consider parameters included in the present analysis. The figures compare the standard deviation of the noise-only covariance and the consider covariance in the three main directions of the local TNW reference frame. Again, the optimum consider parameter variances obtained in case III-B (29.4% and 982.75 ms for SRP and clock-bias consider parameter standard deviations, respectively) have been used for consistency.
SRP model correction causes a fast covariance growth in in-track and normal directions as compared to the noise-only, mildly affecting the out-of-plane    Fig. 8. However, at t 0 its effects are barely visible in the state covariance. This is caused by the consider parameter correction, which at t 0 for the SRP is mostly applied to the SRP coefficient variance, and is transmitted to the state via the linear propagation of the ESTM (which includes the complete augmented state with the SRP coefficient), due to the correlation between the SRP coefficient term and the state vector. This causes the observed growth in the in-track and normal directions. Clock bias correction, on the contrary, provokes a remarkable growth of the in-track covariance at t 0 without any further time evolution, as obsered in Fig. 9. This is expected due to the nature of the clock bias perturbation, which is only affecting the measurements model, being observed mostly at the output of the estimation. It is concentrated in the in-track component, since any error in the observation time is directly translated into an uncertainty in the satellite position in the direction of the motion. For this reason the in-track covariance appears constant in the CB-only case, since the orbital oscillations from the noise-only covariance are several orders of magnitude smaller.
When both corrections are combined (Fig. 10), the clock bias provides not only a higher initial value of the in-track covariance, but also a lower uncertainty limit for this direction along the orbit oscillations. Again, the growth of both intrack and normal directions with time is dominated by the SRP correction. Fig. 8 Covariance evolution in TNW frame from estimation epoch up to 21 days of propagation for a sample orbit comparing 2 cases: without applying any correction ("n" prefix), with a SRP consider parameter variance of 29.4% ("se" prefix) 1 3

Conclusions and Future Work
The results shown in this work indicate that the presented covariance determination methodology is capable of accurately capturing the model error present within the dynamics of an RSO on a simulated SST scenario. It has been successfully applied to GEO, considering SRP and sensor time bias uncertainty. We have tested the robustness of the model by enforcing large perturbation levels, while also ensuring that sensitivity to lower values is maintained. The deviation between the uncertainty introduced as perturbation inputs and the determined consider parameters has remained lower than a 11% throughout the development of this work. Additionally, successful results have been obtained when estimated orbits are used as reference to compute the Mahalanobis distance, indicating the operational suitability of the methodology for operational scenarios. Relevant metrics for covariance realism assessment such as the covariance containment tests have shown that the proposed methodology is able to determine a realistic covariance, applicable to the complete propagation region of interest in GEO and without over-sizing.
The presented methodology allows to maintain the traceability of the different sources of uncertainty, being capable of determining the uncertainty accordingly Fig. 9 Covariance evolution in TNW frame from estimation epoch up to 21 days of propagation for a sample orbit comparing 2 cases: without applying any correction ("n" prefix), with a clock bias consider parameter variance of 982.75 ms ("cb" prefix) in different consider parameters. Regarding the impact of the corrections induced by the retrieved consider parameter variances on the covariance matrices, it was observed that the time bias uncertainty accumulates in the in-track component at estimation epoch, but does not introduce a significant degree of growth rate in the long term. On the contrary the SRP uncertainty correction, though small at t 0 , increases the growth rate in the in-track and normal direction with time.
The next natural step is to use real data from a GEO RSO to further study the capabilities of the methodology. However, the efforts carried out towards realistic simulated scenarios, and the real data results of previous lines of research suggest that the applicability of the proposed methodology for GEO in a real scenario is possible [1]. Further research is under development regarding the optimization process. Robust statistics widely used for the comparison of CDFs such as Cramervon-Mises or Kolmogorov-Smirnov are being tested as cost functions to retrieve the optimum consider parameter variance, allowing also to determine bounds for variances that are statistically consistent to a 2 behavior. Additional analysis needs to be performed to increase the number of consider parameters for different regimes, intending to quantify as much as possible all sources of uncertainty. Although the work presented here focused on covariance realism, analyzing the two first moments of the state PDF, the impact of the proposed methodology in terms of uncertainty realism by retaining higher order moments of the distribution must be considered. Fig. 10 Covariance evolution in TNW frame from estimation epoch up to 21 days of propagation for a sample orbit comparing 2 cases: without applying any correction ("n" prefix), with both consider consider parameter variances applied simultaneously ("2c" prefix)