Precise Orbit Determination Based on the Unscented Transform for Optical Observations

A satellite moving near the Earth is affected by various forces, and the dynamic model mathematically describes the forces. In this study, the dynamic model was designed considering the perturbations caused by the Earth’s asymmetric gravitational field, the atmospheric drag, the solar radiation pressure, and the lunar/solar gravity (Vallado 2001). The dynamic model developed includes a JGM3 model for the Earth gravitational field, an exponential model for atmospheric density, and an analytical approximation for lunar/solar gravity. The Runge-Kutta 5^th method is also employed for a numerical integrator. Note that the detailed equations are not presented in this section.

The optical observations are based on the topocentric equatorial coordinate system (I_tJ_tK_t), and the geometric relationship between the satellite and the ground station is illustrated in Fig. 2 (Vallado 2001). The I_tJ_tK_t coordinate system, which moves the center of the geocentric equatorial coordinate system (IJK) to the ground station, rotates with the Earth, but the axes are fixed. The $\widehat{I}$ and ${\widehat{I}}_t$ axes are parallel to the vernal equinox.

Satellites orbit around the Earth, and there is a difference between the geocentric and topocentric RA/DEC. When optical observations are performed, the topocentric RA/ DEC values (α_t, δ_t) are extracted from the image. α_t is the angle measured from the ${\widehat{I}}_t$ axis to the projected point after projecting the satellite on the I_tJ_t plane, and δ_t is the angle measured from the point projected on the I_tJ_t plane to the position of the satellite. The position vector of the satellite with respect to the ground station, is expressed as $\overrightarrow{\rho }=\overrightarrow{r}-{\overrightarrow{r}}_{site}$, where $\overrightarrow{r}$ and ${\overrightarrow{r}}_{site}$ are the position vectors in the geocentric equatorial coordinate system of the satellite and the ground station, respectively. If the components of $\overrightarrow{\rho }$ in the It J_tK_t coordinate system are represented by (ρ_I, ρ_J, ρ_K), the observations (α_t, δ_t) are expressed as Eq. (1) by Fig. 2.

There are several errors affecting observation data due to the characteristics of optical observations such as light travel time, annual aberration, and diurnal aberration. These errors must be corrected because they could cause several tens of arcsec error in the case of LEO, hence correction models were developed for each error.

As the light speed is finite, a delay occurs by the light travel time from the observed object to the observatory. That is, when the position of the object is estimated using the observation data at the time t, the position of the object is not at the time t but t–τ where τ is the light travel time.

Annual aberration is a phenomenon caused by the movement of the earth’s observer by the Earth revolving around the sun. Because the Earth revolves in an elliptical orbit and the velocity changes periodically, the aberration also changes periodically. The aberration constant κ ′ , which means the maximum displacement of the observed object due to the annual aberration, is 20.4960”. The variations of RA/DEC due to the annual aberration, ΔRA_A and ΔDEC_A, are given by Eq. (2) (Smart & Green 1977).

where λ is geocentric longitude of the sun in the ecliptic plane and ϵ_T is ecliptic inclination.

Diurnal aberration is a phenomenon caused by the velocity of an observer on the surface of the Earth as the Earth rotates. Therefore, it is influenced not only by the observational time but also by the latitude of the observatory. The effect of the diurnal aberration is 1/100 times that of the annual aberration, and the effect is the largest in the equator due to the fastest rotating speed. The diurnal aberration constant is defined as k = cosϕ where ϕ is the latitude of the observatory, and the variation of RA/DEC due to the diurnal aberration is given by Eq. (3) (Smart & Green 1977).

where h is an hour angle, obtained by subtracting RA of object from the local sidereal time.

The use of estimation algorithms can be classified as a batch processing method or a real time processing method. Because both the observation and the OD are performed on the ground, the batch processing method is primarily used for the optical monitoring system, in which observed values are extracted through data processing after the observation images are acquired. The batch estimation algorithm, or batch filter, is the OD technique that estimates a state vector at the epoch time after performing iterative calculation using a set of observed data acquired over a certain period.

The LS batch filter which is used generally includes a process of linearizing a nonlinear system. Thus, there is a disadvantage that if the initial guess value is not correct or the observed data is insufficient, it is possible to occur large estimation error or even divergence (Lee & Alfriend 2007). To overcome these drawbacks, a batch filter based on the UT has been devised (Park et al. 2010). The UT performs nonlinear transformation by generating 2L +1 points based on the probability distribution, where L is the dimension of the estimated state vector. Unlike the LS batch filter, there is no linearization process and therefore it is not sensitive to initial guesses (Park et al. 2010). In this study, a scaled UT method is employed for the nonlinear transform. Three scaling parameters, α, β, and κ, are required. α is an index indicating the range over which the sample is distributed, and β and κ are parameters for calculating the weight of the sample and the covariance matrix. κ is given by κ = 3 – L (Grewal & Andrews 2008). The sigma point generated based on the state vector X_k using the initial covariance matrix P_initial is defined as Eq. (4) (Haykin 2001),

where λ is a parameter calculated by the scaling parameter and can be expressed by Eq. (5) (Grewal & Andrews 2008).

By introducing these UT techniques into a batch filter, the UT batch filter has been developed (Park et al. 2009;Park et al. 2010). In the batch filter, because the state vector at the epoch time is estimated using all observational data of a certain period, it is not necessary to propagate the state vector and the covariance matrix at any times except for the epoch time. Therefore, the sigma point, the state vector, and the covariance matrix at a particular time t_k(χ_i,k, X_k, and P_k respectively) are the same as the estimated values (${\widehat{\chi }}_{i,k},\hspace{0.17em}{\widehat{X}}_k,\hspace{0.17em}\text{and}\hspace{0.17em}{\widehat{P}}_k$; Eq. (6); Park et al. 2010).

After the sigma points are generated at the epoch time, they are propagated to each observational time by the dynamic model F, and the calculated observations at each time can be obtained using the measurement model function h, expressed by Eq. (7) and (8), respectively.

In these equations, y is a computed observation for sigma points, and subscript k means a time t_k (k = 1, ..., N) where N means the total number of observation data. The nonlinear computed observation Y is a weighted sum of y.

where $w_i^m$ is a weight factor for each sigma point.

Therefore, the estimated state vector at the epoch time is updated by (Park 2009).

where ΔZⁱ is the O-C residual matrix, the subtracting the calculated measurement (Yⁱ) from the actual measurement ( Y ) (Park et al. 2010) and Kⁱ is the gain matrix defined as Eq. (11)

where ${\overline{P}}_i^Y$ is the covariance of observation matrix and ${\overline{P}}_i^{xY}$ is the cross-correlation matrix of the state vector and the calculated observations. Note that the superscription indicates the number of iterations because the algorithm is iterated until the convergence criterion is satisfied. The condition for terminating the iteration is given by

where RMS is root mean square, W is the weight matrix, n is the number of different types of observations, N is the total number of measurements, and ϵ is the value of convergence criterion (Vallado 2001).

In order to verify the dynamic model, a state vector (${\overrightarrow{r}}_{epoch}$ =[–4,981.8931, 842.9225, 5,118.8623] km, ${\overrightarrow{v}}_{epoch}$ =[–4.8963, 2.2328, –5.1458] km/s) was propagated using MATLAB and GMAT for 24 hours from the epoch time (February 18, 2016, 18:00:00) with 10 min intervals. In the simulation environment of MATLAB and GMAT, the Earth gravitational field model (JGM3 10 × 10), the atmospheric drag coefficients (2.2), the solar radiation pressure coefficient (1.8), and the integrator (Runge-Kutta 5^th with an accuracy of 10^–9) were set equally. However, the atmospheric drag model (exponential model in MATLAB, and Jacchia Roberts in GMAT) and the lunar/solar gravity (analytic approximation in MATLAB, and DE421 in GMAT) were set differently because there is no identical model in both software. JGM3 10 × 10 was also chosen because 70 × 70 had a propagation position error of 0.5 m for 24 hours, but the calculation took nearly twice as long. To verify the performance of the developed propagator, the propagation error was checked by considering only the gravitational field, and then all perturbations, including the atmospheric drag, solar radiation pressure, and the lunar/ solar gravity. The left and right plots of Fig. 3 show the error over 24 hours when propagating the state vector considering only the Earth gravitational field model and all perturbations, respectively. Due to the position error being within 3 m and the speed error being within 0.003 m/s for the gravity only case, the performance of the JGM3 model was verified. Considering all perturbations, however, the position error was increased up to 100 m and the speed error was increased to 0.15 m/s. The difference was caused because the employed air drag model and exponential model is a simple and static model with low accuracy, while the Jacchia Roberts model, embedded in the GMAT, reflects various phenomena, such as variation of the Earth magnetic field and hydrogen density at high altitudes (Vallado 2001). Thus, the performance of the propagator developed in MATLAB was verified.

The measurement error models, the light travel time effect and the annual/diurnal aberration, were verified using the OWL-Net data that shall be introduced in section 4. We analyzed the OWL-Net measurement residuals computed by the CPF with and without the corrections. The ephemeris contains approximately 145 daily-state vectors (X, Y, and Z positions) of the target, Cryosat-2, and the position error of the ephemeris is several tens of meters (NERC Space Geodesy Facility 2019). Because the velocity is not provided by the CPF, the position vector at the desired time was obtained by interpolation using the Lagrangian method. Table 1 presents the standard deviation (1σ) values of the measurement residuals before and after the correction. After correction, the RA 1σ was reduced from 8–16 to approximately 5–6 arcsec and the DEC 1σ from 4–10 to approximately 2–7 arcsec, which corresponds to the expected quality level of the OWL-Net observation data (Choi et al. 2018). Fig. 4 shows the measurement residuals of the first arc of the observation. The black dots indicate the corrected residuals, whereas the gray dots represent the uncorrected residuals. The uncorrected residuals tend to be biased owing to the aberration, and the residual became zero-mean by correcting the aberration effect. Note that the light travel time effect was 3.4–5.2 ms, the annual aberration was up to 50 arcsec, and the effect of the diurnal aberration was approximately 1/1,000 of the annual aberration.

In order to verify the LS and UT batch filter, artificial optical observation data were created using GMAT, and numerical simulations were performed using cases with various initial orbit errors. In order to realize a virtual LEO satellite, the initial true orbital elements were set to 7,195.0983 km for the semi-major axis (SMA), 1.2780 × 10^–5 for the eccentricity (ECC), 97.5067° for the inclination (INC), 342.7249° for the RA of ascending node (RAAN), 207.9929° for the argument of perigee (AOP), and 286.1304° for the true anomaly (TA). The reference state vector (${\overrightarrow{r}}_{epoch}$ = [–4,981.8931, 842.9225, 5,118.8623] km, ${\overrightarrow{v}}_{epoch}$ = [–4.8963, 2.2328, -5.1458] km/s) converted from the orbital elements at the epoch time (18 February 2016, 18:00:00) was propagated for approximately 8 min with printing the position vector of the ground station in Daejeon, designated as the geodetic longitude of 36.3976°, the geodetic latitude of 127.3757°, and the altitude of 139 m. For the propagation, the Runge-Kutta 9^th method was used for the numerical integration and four perturbations were considered: (1) Earth gravitational field model (EGM96 70 × 70), (2) atmospheric drag model (Jacchia-Roberts), (3) solar radiation pressure, and (4) lunar/solar gravity (DE421). The atmospheric drag coefficient was set to 2.0, and the solar radiation coefficient was set to 1.0.

Using the propagated trajectory, the artificial optical observation data composed of topocentric RA and DEC were generated. The optical observation environment was set at 10 Hz for the observation data frequency, 4 seconds for the camera exposure time, 16 seconds for the observational time interval, and 10° for the cutoff angle. Because optical observations can only be performed at night without sunlight, the altitude of the sun was also considered by calculating the solar ephemeris and the observation data were printed only if the solar altitude was less than –10°. Fig. 5 shows the generated virtual optical observational data. RA distributes from approximately 285° to 150° and DEC is spread from approximately 80° to –40°. The total arc length of the observational data is approximately 8 min, and it consists of 570 points of data, extracted from 27 pictures.

POD simulations were conducted using the artificial optical observation data, setting the weight for the observation data as 3 arcsec in both the RA/DEC, and the convergence criterion as 10^–3. To include the systematic uncertainty in the simulation, a dynamic model for POD simulation was set as one different to that used for generating the observation data: the JGM3 10 × 10 model for Earth gravitational field; the exponential model for atmospheric density, the analytical approximation for lunar/solar gravity, and the Runge-Kutta 5^th method for the numerical integrator.

Various initial errors were employed for the simulations, and the effect of the initial error on the accuracy of POD was analyzed. Table 2 shows the initial position/velocity error and the initial covariance matrix of the six cases, and Fig. 6 shows the result of POD according to the initial orbit error. It shows three-dimensional position error at the epoch time with respect to the artificial ephemeris for each iteration, which changes as the iterative computation. Table 3 and Table 4 show the estimated position error between the POD results of LS/UT batch filter and the true ephemeris in the topocentric RTN (radial, transverse, and normal) coordinate system. R is the line of sight of the satellite as viewed from the ground station, T is the direction perpendicular to the R direction in the orbital plane, and N is the direction perpendicular to the orbital plane.

In all the cases, both filters converged within five iterations, and the estimation error from the true ephemeris was approximately 10 m. The artificial LEO satellite orbits a nearly circular orbit with an altitude of 817 km, and therefore, a Gaussian observational noise of 3 arcsec corresponds to a distance error of approximately 12 m on the sky. Therefore, both estimation algorithms converge to the error value corresponding to the observation data quality, and as the optical observation cannot acquire the three-dimensional distance information of the satellite, the position error of the POD result is concentrated in the radial direction. Further, the position error of the UT batch filter is 1–4 m (8%–33%) lower than the error of the LS batch filter in all the cases, thereby indicating improved performance because the UT batch filter is more robust to nonlinear systems than the LS batch filter.

To analyze the effect of the initial orbit error on the POD, various initial error levels were designated. In this study, several initial orbit determination (IOD) results and propagated two-line element (TLE) were employed as initial guessed state vectors to verify that the POD can be performed independently. Gauss’s technique was employed for IOD, which requires three pairs of observational data, $\left[\left({\alpha }_1,\hspace{0.17em}{\delta }_1\right),\left({\alpha }_2,\hspace{0.17em}{\delta }_2\right),\hspace{0.17em}\left({\alpha }_3,\hspace{0.17em}{\delta }_3\right)\right]$, and determines the state vector, $\left({\overrightarrow{r}}_2,{\overrightarrow{v}}_2\right)$, at the intermediate point. The first arc of the observational data consists of a total of 447 data points that were used for IOD, and there were various combinations of observational data: the 232^nd data point was fixed as the intermediate point, $\left({\alpha }_2,\hspace{0.17em}{\delta }_2\right)$, and the first point, $\left({\alpha }_1,\hspace{0.17em}{\delta }_1\right)$, and the last point, $\left({\alpha }_3,\hspace{0.17em}{\delta }_3\right)$, were selected from the remaining data. Each IOD result was propagated to the epoch time through the GMAT, and the three-dimensional position error was computed by comparing to the CPF ephemeris, where Fig. 8 shows an error of up to 700 km with a 1σ value of 40.2464 km. To analyze the effect of the initial orbit error on the POD, three combinations were selected from approximately 50,000 combinations and used as the initial guess state vector of the POD. Table 7 lists the position error of the selected initial guesses with respect to CPF in the topocentric RTN coordinates. The position errors of cases 1–3 are concentrated in the radial direction of the ground station because the optical observations cannot acquire the 3D distance information of the satellite. In case 4, however, the concentration is not apparent because the propagated TLE was used as the initial state vector. The initial covariance matrix was set as P_initial = diag([1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/s]²).

Table 8 shows the RMS of the post-fit residuals of POD, and Figs. 9 and 10 show the post-fit residuals of case 1 and case 4, respectively, where the gray dots represent the post-fit residuals of the LS batch filter, and the black dots represent those of the UT batch filter. In cases 2–4, the post-fit residuals converged to several arcsec at both the LS and UT batch filters. However, in case 1, where the initial position error is more than 100 km, the LS batch filter converged to the level of several tens of arcsec. In other words, the UT batch filter converged to the residual level corresponding to the data quality, irrespective of the initial orbit error level, but the LS batch filter exhibited degraded OD performance when the initial orbit error was very large. However, in case 4 which used TLE information, the results of both LS and UT batch filters were similar and the standard deviation value was small.

Fig. 11 shows the estimated position error relative to the CPF ephemeris for each iteration, and Tables 9 and 10 show the topocentric RTN position error of the LS and the UT batch filter, respectively. The UT batch filter converged to approximately 300 m in all cases, regardless of the initial error. However, the LS batch filter only converged to the error of several tens of meters in case 3 which had small position error and in case 4 which used TLE information, while it converged to errors of several hundreds of meters in cases 1–2.

In case 2, the estimated position error was approximately 1,860 m, which is more than three times larger than cases 3–4 results even though the initial position error is 20 km, smaller than that of case 4, and the post-fit residual converged to several arcsec. It might be caused by initial velocity error which is not evaluated by comparing to the CPF data, and the initial guesses were compared to the TLE for analyzing the velocity error indirectly. Tables 11 and 12 exhibit the initial error with respect to TLE. In case 3, the three-dimensional position and velocity errors were 74 km and 78 m/s, which is an order of 1,000 times. In case 1 and case 2, however, the position errors were approximately 120 km and 71 km, while the velocity errors were 605 m/s and 169 m/s, where the velocity error is 2 to 5 times larger than 1/1,000 of the position error. In summary, the initial guesses of cases 1–2 have large velocity error relative to the position error. It can also be concluded that UT batch filter robustly converged regardless of initial velocity error, whereas the LS batch filter is sensitive to the initial velocity error, converging to the local minimum value.

In order to analyze the influence of the initial covariance matrix on the POD, the initial covariance was variously set while the propagated mean of all IOD results was used as an initial value. The initial position error in the topocentric RTN coordinate system is [14.9560, 3.8612, 0.1869] km, where the error is concentrated in the radial direction of the ground station. The size of the initial covariance matrix is diag (0.3² × [1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/ s]2) for case 1, diag (0.5² × [1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/s]2) for case 2, diag (12 × [1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/s]2) for case 3, diag (22 × [1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/s]2) for case 4, and diag (3² × [1 km, 1 km, 1 km, 10^–3 km/s, 10^–3 km/s, 10^–3 km/s]2) for case 5.

Table 13 shows the RMS of the post-fit residual. In cases 4–5, i.e., with large covariance, both LS and UT batch filters converged similarly to 5 arcsec level in RA and 3 arcsec level in DEC. In cases 1–3, however, the smaller the initial covariance matrix size, the greater the post-fit residual, up to several tens of arcsec in the LS batch filter. In other words, the UT batch filter converged to the residual level corresponding to the data quality in all cases, irrespective of the initial covariance matrix size, but the LS batch filter exhibited degraded OD performance when the initial covariance matrix was small. The initial covariance matrix implies the uncertainty of the initial state vector. When the initial covariance matrix is smaller than the error level of the initial state vector, the LS batch filter converged to the wrong states.

Fig. 12 shows the position errors of POD relative to CPF for each iteration, where the upper and bottom are the LS and UT results, respectively. In the case of LS batch filter, the estimation error was several hundred meters in cases 4–5, whereas the position error increased from 1 km to 6 km in cases 1–3. Conversely, the UT batch filter converged to several hundreds of meters, although the estimation error slightly differs for each case. It can be concluded that reliable performance is only achievable for the LS only if the initial covariance is properly designated, but the UT converges robustly regardless of the initial covariance matrix.

In this section, arcs 1–4 of the OWL-Net data were used to analyze the effects of observation arc length and the number of data points on POD. The arc length is the time between the first data point and the last data point of each arc, which means a single observation image.

The initial guess state vectors were determined in the same way as above: fixing the middle point as the 232^nd point for arc 1, the 141st point for arc 2, the 171^st point for arc 3, and the 103^rd point for arc 4, and choosing the combinations for the Gauss method. Table 14 shows the position error relative to the CPF of the initial guess state vector in the topocentric RTN coordinate system. The initial guesses were selected to contain initial position error of 14.3 km, and the initial covariance matrix is set as P_initial = diag (32 [1 km,1 km,1 km,10^–3 km/s, 10^–3 km/s, 10^–3 km/s]²).

Table 15 shows the RMS of the post-fit residual of the POD. After estimating the precise orbits of all the arcs, the residual converged to the level of 4–6 arcsec in RA, and 2–4 arcsec in DEC, corresponding to the data quality. Fig. 13 shows the POD result as a position error relative to the CPF for each iteration, where the upper shows the LS results and the bottom shows the UT results. Although the measurement residuals converged in all cases, the estimated position vectors did not. For UT batch filter, the estimated position error converged to several hundreds of meters in the case of arcs 1–3. For LS batch filter, the estimated position error was several hundreds of meters in the case of arcs 1 and 3 where the number of observational data points is more than 440. However, in arc 2, where the number of data points is less than 400, the error converged to approximately 1 km level. As the length of arc 2 is approximately 20 seconds longer than arc 3, but the number of data are approximately 50 fewer, the LS batch filter can be analyzed to be more sensitive to the number of data points in the arc than the arc length. The length of the arc 4, meanwhile, is very short, about 60 seconds, and the number of data points is 157. In that arc, both LS and UT batch filter did not converge, and the estimated position errors were 10 km and 6 km, respectively. In conclusion, the POD result could not converge with LS or UT batch filter if more than 400 points of data cannot be obtained for more than 3 min.