Development, Demonstration and Validation of the Deep Space Orbit Determination Software Using Lunar Prospector Tracking Data

Orbit determination is procedure to obtain parameters that are related to the motion and the tracking data of the spacecraft. The procedure requires a system model and an estimation algorithm to predict the tracking data and estimate parameters using actual tracking data. The system model can be divided into a dynamic model and a measurement model. The dynamic model includes the force models acting on the spacecraft and a numerical integrator for the orbit propagation. The measurement model includes the tracking algorithm and hardware characteristics to simulate the tracking data. In the DSODS, each model is implemented in two submodules: OP module and DS module. Furthermore, the OD module consists of an estimation algorithm and some useful subroutines. Kim et al. (2017) presented an overview of the DSODS, including its construction, data flow and user interfaces of the submodules. Fig. 1 shows the brief data flow of DSODS and subroutines of the module, whose details are explained subsequently.

In case of the OP module, a general mission analysis tool (GMAT) is employed as a third-party software component using a MATLAB-executable (MEX) interface (Kim et al. 2017). The GMAT was developed by NASA GSFC and has been used in missions such as the Lunar Crater Observation and Sensing Satellite (LCROSS) and the Acceleration, Reconnection, Turbulence and Electrodynamics of the Moon’s Interaction with the Sun (ARTEMIS). The force model and numerical integrator (propagator) are determined and validated using multiple test scripts (Hughes et al. 2014). The OP module receives the settings for the orbital prediction, such as a gravity model and accuracy of the propagator, from the OD module, and transfers the propagated states set, i.e., ephemeris, in the form of a satellite tool kit (STK) format.

Alternatively, the DS module is developed by Yonsei to simulate the tracking data. NASA’s deep space network (DSN), near earth network (NEN) and universal space network (USN) are usually employed for deep space missions (Policastri et al. 2015), and these stations provide various types of tracking data such as tracking angles, un-ramped two-way range and Doppler, and the difference between the uplink and the downlink frequency divided by the transponding ratio of the spacecraft (Chang 2015). In the DSODS, un-ramped two-way range and Doppler are modeled by the DS module, where the corrections for tropospheric delay and antenna offset delay are involved. I should be noted that the ramped range and Doppler, and the corrections for ionospheric delay and charged particle delay are being developed and will be added to the DS module. The DS module simulates tracking data by utilizing the configuration for the DS, the ephemeris and the mesh points are transferred from OD module, and the simulated tracking data is returned to the OD module.

In the OD module, a batch least squares algorithm is implemented. The batch type estimation is a post-processing method that treats all the observed data of an arc to determine the states at a specific epoch, while the sequential estimation algorithm estimates the states sequentially at each observation time. Batch estimation is a simple, more accurate and robust technique compared to sequential estimation algorithms, like the Kalman filter (Crassidis & Junkins 2011). The functionality requirements and subroutine algorithms are addressed in the following sections.

The OD module is required to provide some necessary range of operations. Table 1 presents the functionality requirements of the OD module (Song et al. 2016). The functionality code categorizes the functionality requirement, and can consist of some sub-codes. In case of the ODF-0 (to read and write information from/to a text file), for instance, there are two sub-codes to print out and load information, and each function is developed within the OD module. Note that the functionality requirements and the correlated codes are subject to change should the details of the mission requirement and operating concepts change.

The OD module consists of multiple subroutines to estimate parameters, improve performance and support additional analysis. In this section, some representative subroutines, the weighted batch least squares, the specific/iterative data editor and the data handler are briefly explained.

The batch least squares algorithm estimates parameters using a system model (F):(1)

In case of the OD, Y is the tracking data, є is the measurement error vector, and X is a ‘solve-for’ vector, which is to be estimated such as the position or velocity vector (Schutz et al. 2004). The system model is usually composed of the dynamic model and the measurement model.

The estimation can be referred to as an optimization, which minimizes the performance index (I), weighted sum of the squares of the measurement residual and the difference from a priori states:(2)

where ${\text{P}}_{{\text{ΔX}}_{\text{0}}}$ and X₀ are the a priori covariance matrix and ‘solve-for’ vector, respectively, and W is the weight matrix. The a priori covariance matrix and weight matrix correspond to the a priori uncertainty and measurement noise:(3)

where σ_j is the standard deviation of the j^th element of the ‘solve-for’ vector, n is the dimension of the ‘solve-for’ vector, ϱ_k is the noise level of the k^th tracking data measurement and m is total number of tracking data measurements.

The optimal solution is evaluated by the Newton-Rapson method, that uses 1^st order approximations under the assumption of differentiability of the system model and the suitability of the ca priori states. Through several iterations, the ‘solve-for’ vector is updated using the variate of the i^th iteration (ΔX_i) (Long et al. 1989):(4)

where H is the Jacobian matrix, a partial derivative of the system model, that is numerically evaluated by that central differencing method in the OD module.

The convergence criterion is the cost reduction ratio (γ):(5)

The estimation sequence stops if the ratio becomes less than a small tolerance (ε), set as 10^-3 by default, and also terminates when the number of iterations reach a designated maximum number.

The covariance matrix (P), which corresponds to the precision of the estimation, is evaluated after convergence by the following expression:(6)

Data editing is done to reject any abnormal observation data that diminishes the accuracy of the solution from the orbit determination process. There are two types of data editors in the developed software. The first one is the specific editor, where the user designates a specific type, site and date of the observation. It is useful when the user has the empirical skills and when there are known problems with some observatories or dates. Another editor is the iterative editor, which is automatically executes every iteration. The iterative editor rejects the observation data if the data does not meet the conditions that are calculated by the measurement residuals and the states of the iteration (Long et al. 1989).(7)

In the above expression, w_j is a weight component corresponding to the j-th observation data, y_j is a j-th measurement residual, k is a constant multiplier and K is an additive constant. The predicted root mean square, PRMS, is evaluated by following form.(8)

Fig. 2 is an example of iterative editing using the weighted residual method, where two abnormal observation data points, which are circumscribed, are edited.

All information about the orbit determination including the dynamic model, tolerance, maximum iteration number and the estimated result are saved in the orbit determination data handler. The handler can generate a residual plot, compare the estimation results to the ephemeris and print the information to an ASCII file. The handler also can read the output file, load the configuration and the result from the file and utilize this information for the next orbit determination. Fig. 3 shows the residual plot and the output ASCII file.

The orbit determination environment including the force model, the propagator options, the a priori covariance matrix and the weights for the measurement, is compiled in Table 2. The force model and the propagator are implemented in the GMAT, and the measurement weight is assigned based on the nominal noise level of the S-band tracking data (Chang 2015). The ‘solve-for’ parameters are threedimensional position and velocity vectors, and the range and Doppler biases of each tracking station. Fig. 4 shows the 8 hrs' tracking data, un-ramped two-way range and Doppler data provided by the 2 DSN stations, the DSS-24 (Goldstone, USA) and the DSS-42 (Canberra, Australia), also the DSS- 27 (Goldstone, USA) provides only the un-ramped two-way Doppler data. The NASA PDS¹ ephemeris is considered as the true state, with an accuracy of approximately 1 km in position and 1 m/s in velocity².

A set of a priori state vectors are generated by adding pseudo Gaussian errors to the true states, where the errors correspond to an a priori covariance matrix (Table 2). The a priori error is set as 1 km for position and 0.001 km/s for velocity, in accordance to the ephemeris accuracy of the Lunar Prospector. Note that the a priori error of biases are ignored, i.e., a priori biases are fixed as zero since the robustness is checked only for the position and velocity vector of the spacecraft. The orbit solutions are evaluated using the generated a priori states set of 20 state vectors, and are statistically analyzed.

Fig. 5 shows the statistical covariance ellipsoid and the distribution of the solutions, where the major axis of the ellipsoid is normal to the direction towards the Earth. Since the tracking data, range and Doppler, is measured in the radial direction, uncertainty in the radial direction is much smaller than the along-track and cross-track directions. The solutions are distributed with a three-dimensional standard deviation of 18.66 m, which satisfies the mission requirement of the KPLO. The precision, however, is predicted as ~ 4 m by the estimated covariance matrix (as mentioned in Section 2.3), which includes the orbit solutions within a 5σ ellipsoid. Since the statistical consistence of the developed software has already been confirmed by that the estimated covariance matrix and the distribution of the actual solutions are consistent under the simplified environment (Lee et al. 2016), the causes that under-estimate the covariance matrix in comparison to the solution distribution in this study are analyzed considering the characteristics of the actual system. The first cause is the process noise, which denotes the unmodeled phenomena in the dynamic and measurement models. The covariance matrix evaluated by a batch type estimator does not contain process noise, but instead a a correlation matrix, which presents uncertainty of system-related parameters, ‘consider parameter’. The system model parameters, such as the geopotential coefficients and the solar radiation coefficients, can be used as the consider parameter, and their uncertainties can be modeled in the correlation matrix. In this analysis, however, no consider parameter is employed, resulting in an additional estimation error. The other cause is the inaccuracy of the a priori covariance and weight matrices. Since the uncertainty of a priori states and the noise level of the tracking data are not exactly determined, the covariance matrix has a corresponding error. Furthermore, the truncation error and round-off error also contribute to the estimation error.

Through the stress test, the robustness of the DSODS can be demonstrated within an a priori error margin of about 1 km. Moreover, the estimated covariance matrix can be regarded as an indicator of the precision of the orbit solution.

A comparison between solutions achieved through different systems is one of the methods to assess the accuracy of an orbit solution (Schutz et al. 2004). To confirm whether the developed software can prove suitably reliable for a mission, the orbit solution is evaluated using Lunar Prospector tracking data propagated over 24 hr, and is compared to the NASA PDS ephemeris.

In this study, an orbit determination is executed once in a day over the lunar orbit using one-day long tracking data, range and Doppler data, since the mission operation concept of the KPLO is not fixed yet. The arc length is 24 hr on November 25, 1998, using the tracking data provided by several DSN stations, as configured in Fig. 6. Un-ramped two-way range and Doppler data was provided by 5 DSN stations, DSS-16 in Goldstone, USA, and DSS-24, DSS-42, DSS-61 and DSS-66 in Madrid, Spain, with additional unramped two-way Doppler data supplied by a single DSN station, DSS-27 in Goldstone, USA.

After the orbit determination, the standard deviations of the residuals are 4.28 m for range and 1.13 mm/s for Doppler, which correspond to the noise level of the S-band tracking data (Fig. 7). The system model however, should be improved since the residuals are not completely random distributions, i.e., include the existence of a systematic error. The position difference of the orbit solution with respect to the NASA PDS ephemeris is presented in Fig. 8. The variation by axis is conspicuous in the orbital frame, i.e., radial-transverse-normal frame, although not in the original inertial frame. The radial difference is the smallest, while the transverse-normal frame differences show an increase from 2 to 4 times after 24 hr. This is also due to the characteristics of the tracking system which measure the range and Doppler in the radial direction. Note that although the difference accumulates with time due to the nonlinearity of the dynamic system, the RMS of the difference over 24 hr is 39.51 m which is significantly lower than the orbit accuracy of the Lunar Prospector (1 km). Through solution comparison, the performance of the DSODS is deemed appropriate to support the KPLO.

Overlapping is another technique to measure the accuracy of the orbit solution. The differences between separated arcs that share some part of tracking data are evaluated and considered indicators of orbit accuracy (Schutz et al. 2004). In this analysis, an 8 hrs' arc is designated and neighboring arcs are separated with 6 hrs' interval, that share 2 hr of tracking data. The tracking data used for the overlap analysis is the data of November 25, 1998, same as that used in the previous section.

Fig. 9 shows the residuals of each arc, where two types of boxes, solid and dash, represent overlapped tracking data. The standard deviations of the range residuals are 1.40 m, 3.64 m and 1.93 m for each arc, and that of Doppler are 0.33 mm/s, 0.43 mm/s and 0.92 mm/s (Table 3). The residuals correspond to the nominal noise level of the S-band tracking data, signifying that the orbit solution is well converged by iterative process. However, the difference between the residuals of arcs is at least double due to a systematic error which needs further improvement. Fig. 10 is the position variation between the overlapped arcs with RMS position differences of 36.68 m and 50.14 m, with most differences concentrated in the transverse and normal directions. These differences correspond to the results from the previous section, confirming the consistency of the software and validating the accuracy of the DSODS to several tens of meters.

In addition, Fig. 11 presents the position difference of each arc compared to the NASA PDS ephemeris, with the RMS differences of 54.41 m, 20.23 and 63.55 m, respectively. These differences are also of the same order as the previous section. Here too the differences in the radial direction are minimum and almost same for all arcs, while those in the transversenormal directions have considerable variations between arcs. The discordance between the residual level and the position difference is interpreted as a systematic property caused by the shortage of information in the transverse and normal directions due to the characteristics of the tracking system.