A More Reliable Orbit Initialization Method for LEO Precise Orbit Determination Using GNSS

: Precise orbit determination (POD) using GNSS has been rapidly developed and is the mainstream technology for the navigation of low Earth orbit (LEO) satellites. The initialization of orbit parameters is a key prerequisite for LEO POD processing. For a LEO satellite equipped with a GNSS receiver, su ﬃ cient discrete kinematic positions can be obtained easily by processing space-borne GNSS data, and its orbit parameters can thus be estimated directly in iterative manner. This method of direct iterative estimation is called as the direct approach, which is generally considered highly reliable, but in practical applications it has risk of failure. Stability analyses demonstrate that the direct approach is sensitive to oversized errors in the starting velocity vector at the reference time, which may lead to large errors in design matrix because the reference orbit may be signiﬁcantly distorted, and eventually cause the divergence of the orbit parameter estimation. In view of this, a more reliable method, termed the progressive approach, is presented in this paper. Instead of estimating the orbit parameters directly, it ﬁrst ﬁts the discrete kinematic positions to a reference ephemeris in the form of the GNSS broadcast ephemeris, which construct a reference orbit that is smooth and close to the true orbit. Based on the reference orbit, the starting orbit parameters are computed in su ﬃ cient accuracy, and then the ﬁnal orbit parameters are estimated with a high accuracy by using discrete kinematic positions as measurements. The stability analyses show that the design matrix errors are reduced in the progressive approach, which would assure more robust orbit parameter estimation than the direct estimation approach. Various orbit initialization experiments are performed on the KOMPSAT-5 and FY3C satellites. The results have fully veriﬁed the high reliability of the proposed progressive approach.


Introduction
With the great progress of Global Navigation Satellite System (GNSS), high-precision orbit determination using space-borne GNSS data has been the most mainstream technology for the navigation of low Earth orbit (LEO) satellites, due to its global coverage, abundant observations and low-cost [1,2]. Only with a GNSS receiver onboard LEO satellites, a lot of continuous GNSS measurements are generated without time and geography limitation. By processing the space-borne GNSS data, precise orbits at centimeter level can be achievable, hence the technique of precise orbit determination (POD) using GNSS has been applied to a wide range of LEO satellites nowadays [3][4][5].
In the dynamical orbit determination, the orbit parameters, including the position and velocity at the reference time and several piece-wise dynamical coefficients, such as the atmospheric drag coefficient, the solar radiation pressure coefficient and empirical acceleration coefficients, should be initialized first by fitting the discrete positions, before being accurately estimated [6]. This process, which is called "orbit initialization", is a key prerequisite for the successful operation of LEO Consequently, a more reliable orbit initialization method is important to avoid such fatal error and guarantee stable POD processing.
This paper proposes a more reliable orbit initialization method to improve the stability of the precise LEO orbit estimation. The new method first constructs a reference orbit for a LEO satellite by estimating a set of orbit parameters in the form of the GNSS broadcast ephemeris based on the SPP/PPP discrete kinematic positions. For the convenience, this set of orbit parameters is referred as the LEO reference ephemeris. Then, the LEO reference ephemeris can be used to generate a series of tabular positions and velocities. Although the LEO reference ephemeris can't describe the orbit accurately and the errors of generated tabular positions may be of hundreds or thousands of meters, the resultant reference orbit is nevertheless smooth and close to the true orbit, such that the precise orbit parameters can be estimated stably with the starting position and velocity vectors at the reference epoch of guaranteed accuracy. In a step-by-step approach, the estimated orbit parameters in the previous step are used to compute the starting orbit state for the next direct estimation, and accurate orbit parameters are estimated directly by taking the SPP/PPP discrete kinematic positions as observations for the POD. In a summary, the orbit parameters are obtained by adopting a step-wise manner in the new method, so it is called "the progressive approach".
In the following section, the direct approach will be briefly introduced first. Then, the progressive approach is presented in detail. Next, the corresponding methods and tests will be designed to have an in-depth understanding on the estimation stability of these two approaches. After that, these two approaches will be applied to the orbit initialization experiments of real LEO missions such as KOMPSAT-5 and FY3C. Some cases will be compared and discussed in detail to validate the high stability of the progressive approach. Finally, some conclusions are made.

The Direct Approach
The equation of the motion of a LEO satellite can be expressed as: ..
where (r, . r, .. r) are the position, velocity and acceleration vectors of the satellite in a geocentric inertial coordinate frame, respectively, and p represents the dynamical coefficients in the force models, including the atmospheric drag coefficient C d and the solar radiation pressure coefficient C r . To account for those un-modeled forces that act on the LEO satellite or for the inaccurate force models, the empirical accelerations w = (a R ,a A ,a C ) in the radial (R), along-track (A) and cross (C) directions are customarily introduced in the POD processing, where some empirical parameters, including the one-cycle-per-revolution coefficients in the R/A/C directions, are estimated [17]. These to-be-estimated empirical parameters are also treated as a part of the vector p. f (·) is a function expression that represents the acceleration .. r, and is related to the position, velocity and dynamical coefficients of the satellite. According to Equation (1), the first-order differential equation of the motion of the satellite is formed: A dynamical orbit can be expressed by a set of orbit parameters in the form of the position r 0 and velocity . r 0 at a reference time t 0 and the dynamical coefficients p. In order to express the orbit more precisely, the dynamical coefficients p are always set to be piece-wise constant (PWC), thus they are called PWC parameters. All to-be-estimated orbit parameters are abbreviated as a vector: X = r 0 , . r 0 , p . According to differential Equation (2), the position r at any epoch t can be expressed by an implicit function: where, it is difficult to compute the explicit formulation of g(·) because of the high non-linearity of differential Equation (2). However, it is obvious that, using a series of discrete kinematic positions (r(t 1 ), r(t 2 ), . . . , r(t n )) at the subsequent epochs (t 1 ,t 2 , . . . ,t n ) as measurements, the orbit parameters X can be estimated.
In the estimation process, setting the starting state of the orbit parameters X = X* first, observation Equation (3) can be linearized as: (4) where, dX* is the estimated correction vector to the starting state X*. r*(t) = g(X*,t) is the computed position based on the starting state X*. Φ(X * , t, t 0 ) = ∂g(X, t)/∂X X=X * is the state transition matrix, which represents the first-order partial derivative of the position with respect to the estimated orbit parameters. o(X*,dX*,t) is the remainder of linearization, which is relative to the starting state X* and the correction vector dX*. Although g(·) cannot be expressed by an explicit formulation, r*(t) and Φ(X * , t, t 0 ) can be computed by the numerical integration of differential Equation (2) [18]. dX* is the unknown parameter to be estimated, so the remainder o(X*,dX*,t) is unknown and cannot be computed. It should be neglected in order to estimate dX* by the least squares method.
According to Equation (4), the orbit parameters are estimated iteratively, and the iterative estimation process of the direct approach is designed as Figure 1. The processing scheme mainly consists of three steps. First, the starting state of the orbit parameters is initialized as X*, where the position r 0 and velocity . r 0 at the reference time t 0 are valued directly by the interpolation of the SPP/PPP discrete kinematic positions, and other dynamical coefficients are always set as 0.0 or their nominal values such as 2.2 for the drag coefficient. Then, based on the starting state X*, the predicted positions r*(t) and state transition matrix Φ(X * , t, t 0 ) within the whole estimated arc are computed by the numerical integration of the motion equations. The positions r*(t i ) and state transition matrix Φ(X * , t, t 0 ) at each measurement epoch (t 1 ,t 2 , . . . ,t n ) can be computed by an interpolation method. Thus, error Equation (4) can be constructed. Finally, the normal equation is obtained, and the correction dX* is estimated by the least squares method. Generally, the second and third steps should be performed iteratively, until the convergence condition is satisfied. Obviously, the oversized errors in the starting state X* may cause the divergence of the iterative estimation of orbit parameters.

The Progressive Approach
The progressive approach doesn't estimate the orbit parameters directly, but it first fits the discrete kinematic positions to a reference ephemeris Y, which has the same parameters in the GNSS broadcast ephemeris [19]: where, A is the semi-major axis, i 0 , M 0 and Ω 0 are the orbital inclination, mean anomaly and right ascension of the ascending node at the reference time t 0 , respectively, e is the eccentricity, and w is the argument of perigee. ∆n represents the difference between the average and computed angular velocity. i and .
Ω are the rates of orbital inclination and the right ascension of the ascending node, respectively. C uc and C us are the amplitudes of the cosine and sine correction terms to the latitude angle, respectively. C rc and C rs are those to the orbital radius, and C ic and C is to the orbital inclination. The first six parameters √ A, i 0 , e, Ω 0 , w, M 0 ) are the Keplerian orbit elements, which can be computed directly by from the position and velocity vectors. The following three parameters (i, Ω, . ∆n) are the supplements to the main Keplerian orbit elements. The last six parameters (C uc , C us , C rc , C rs , C ic , C is ) are to account for the perturbation effects on the latitude angle, radius and inclination. Using the SPP/PPP discrete kinematic positions as measurements, the parameters in the reference ephemeris can be estimated. The starting state of the six main Keplerian orbit elements √ A, i 0 , e, Ω 0 , w, M 0 ) are computed from the position and velocity r 0 , . r 0 at the reference time t 0 , and other parameters are set as 0.0. Regarding the starting state of the reference ephemeris Y* as the nominal state, observation Equation (3) can be linearized as: ....
where, dY* is the estimated correction to the starting state Y*. r*(t) is the position computed from the starting state Y*. ϕ(Y*,t,t 0 ) is the state transition matrix, which represents the first-order partial derivative of the position with respect to the reference ephemeris. Both r*(t) and ϕ(Y*,t,t 0 ) can be computed by using an explicit analytical formula [20]. υ(Y*,dY*,t) is the remainder of linearization. Similarly, the unknown remainder υ(Y*,dY*,t) should be neglected in order to estimate dY* by the least squares method.
The estimation processing scheme of the progressive approach is shown in Figure 2. The processing scheme mainly consists of five steps. First, the starting state of Keplerian orbit elements is computed from the position and velocity vectors at an arbitrary epoch. Next, the parameters in the reference ephemeris are estimated iteratively by using only the SPP/PPP discrete kinematic positions as measurements. This generates a reference orbit. Although the position errors of the reference orbit may be at the level of hundreds or even thousands of meters, the determined reference orbit is smooth and continuous, and close to the true orbit. Then, the orbit parametersX 1 are estimated by using the tabular positions of reference orbit as measurements. Unlike the SPP/PPP discrete kinematic positions, the reference orbit is smooth and follows the orbit dynamics. So it can provide a starting position and velocity with enough accuracy to estimate a set of orbit parameters reliably that coincide with the tabular positions of reference orbit. Finally, the estimated orbit parametersX 1 are used as the starting state, and the raw SPP/PPP discrete kinematic positions are used as measurements, to determine the final orbit parametersX iteratively. Because the starting state of the orbit parametersX 1 , which represents the reference orbit, is close to the true orbit, the errors in the starting state would have negligible impact on the iterative estimation of final orbit parameters.

Effect of the Starting State Error
According to Figure 1, the iterative estimation of orbit parameters in the direct approach is at the risk of divergence if the errors in the starting state X* are too large. In the linear Equation (4), X* appears in three terms: (1) L(X*), (2) B(X*) and (3) o(X*,dX*). L(X*) is the difference between the observed and computed positions, which is used as measurements to estimate the correction d X*. B(X*) is the design matrix, and o(X*,dX*) is the linearization remainder which should be small enough to neglect. If X* has large errors, the neglect of remainder o(X,dX*) and errors in B(X*) will severely distort error Equation (5), which may lead to a divergence of the iterative estimation. The neglected remainder o(X*,dX*) can be treated as the linearization error, and the errors in the design matrix B(X*) is called as the design error, which are both originated from the errors of the starting value X*. If the start value X* is not equal to the true value, the linearization error cannot be avoided, but the iterative estimation can reduce the linearization error, assuming that there is no error of the design matrix B(X*) or the design error is controlled well. Therefore, the design error is essential to the success of iterative estimation. Setting the true value of the orbit parameters as T x and the true design matrix as M, the error of the starting value X* is dX* = T x − X*, and the error of the design matrix B(X*) is M − B(X*). It is difficult to express the effect of the design error if using M − B(X*). The expression is more meaningful to describe the design error in the form of: Here, ERR B is a 3-dimensional vector that expresses the influence of the design matrix B(X*) on the position. The larger ERR B is, the more probably the divergence of iterative estimation occurs. According to Equation (7), the error of starting value dX* decides the design error. If dX* is small enough, ERR B is also not enough to affect the iterative estimation. If dX* is large, ERR B will probably be too large to impact the success of the orbit parameter estimation.
According to Figure 2, the stability of the progressive approach depends on the estimation of parameters in the reference ephemeris. Only when the reference orbit in the form of the reference ephemeris is generated accurately, can the orbit parameters converging to true values be estimated. The generation of the reference orbit is also a process of iterative estimation. Similarly, setting the true parameters in the reference ephemeris as K x , and the true design matrix as D, the error of the starting value Y* is dY* = K x − Y*, and the design error can be formed as:

Stability Analysis Tests
The starting state error dX* consists of three parts: the starting position error dr*, the starting velocity error d /s will lead to the maximum design errors of ±200, ±20,000, ±500,000, and ±2,000,000, respectively. The figures show that the design error increases by approximately proportion to the squares of the starting position and velocity error, and it accumulates rapidly as the orbit integration moves forward. In addition, the design error caused by the starting dynamical coefficient error is relatively small, which is only within the range of ±50.0 m. According to the patterns in Figure 3, the design error could reach tens or hundreds of kilometers level if the start position error is up to several or tens of kilometers, while it can soar up to the equivalent level just with the starting velocity error at sub-or several m/s level. Obviously, the design error is more sensitive to the starting velocity error. The main reason for the higher sensitivity to velocity error is that, not only could the velocity error itself produce design errors, but also it can cause the position error accumulated over time away from the reference epoch. Although a rigorous mathematical equation is difficult to derive, an approximate correlation between the starting velocity error at the reference epoch and the position error at an epoch ∆t away from the reference epoch may be given as ∆p = ∆ν · ∆t to express the influence of the starting velocity on the position error, where ∆ν is the starting velocity error, ∆p represents the position error. Assuming that the starting velocity error is 0.5, the position error will increase up to 43.2 km approximately after 24 h of orbit integration (= 86,000.0 s). For the LEO satellites at the altitude of hundreds of kilometers, the design error of tens or hundreds of kilometers is likely to cause the divergence of orbit parameter estimation, as a matter of experimental experiences, because it may make the estimated orbit deviate severely from the true orbit. Nevertheless, the staring position error of kilometer level at the reference time can be avoided easily as long as the GNSS pseudo-range and carrier-phase measurements are processed with a simple outlier detection scheme in SPP/PPP, but unfortunately, it is very difficult to prevent the starting velocity error of subor several m/s at the reference time absolutely. On the one hand, the position errors of SPP discrete kinematic results are usually at meter level, so the velocity error is easily to reach sub-or several m/s if the velocity is computed by the interpolation of kinematic positions. On the other hand, although the discrete position accuracy of PPP results is always expected at centimeter or decimeter level [6], and the Doppler measurements also can be used for high-precision velocity determination, it is still difficult to preclude the velocity error of sub-or several m/s occurring at certain epochs completely. It is because that, if only a few GNSS satellites are tracked with poor geometric strength, or only less accurate Doppler measurements are available for some LEO satellites at some arcs, the velocity error at certain epochs may still be up to sub-or several m/s [21]. Because it is critical to have an accurately enough starting velocity at the reference epoch for the subsequent POD, one should avoid choosing a reference epoch at which the velocity error happens to be large. However, there would be no external precise orbits as reference to make such choice, the more reliable progressive approach is proposed in this paper to make sure the determined velocity at the reference epoch have sufficient accuracy. Therefore, the orbit parameter estimation is at the risk of having oversized starting velocity error, no matter SPP or PPP discrete kinematic results are used as measurements. It also can be concluded that the oversized starting velocity error is the main possible factor causing the divergence of orbit parameter estimation.
Similarly, the design errors in the X/Y/Z direction and the 3D total error |ERR A | for the nine cases (P1/P2/P3/P4/V1/V2/V3/V4/H1) in reference ephemeris estimation are shown in In the least squares estimator, if the design error of the error equation is too large, the normal equation may become singular, which would result in the divergence of parameter estimation. All in all, the smaller the design error, the more stable the parameter estimation. Comparing Figure 4 with Figure 3, it is obvious that the design errors in the reference ephemeris estimation caused by the same starting position and velocity errors P1/P2/P3/P4/V1/V2/V3/V4 are far less than those in the orbit parameter estimation. Even if the starting velocity error is 10 m/s, the corresponding design error in the reference ephemeris estimation is still less than ±30 km, which is far less than that of ±2000 km in the orbit parameter estimation. It demonstrates that the starting state error has less negative effect on reference ephemeris estimation than on orbit parameter estimation. Although the design error in the reference ephemeris estimation in the case H1 is more than that in orbit parameter estimation, it is still less than ±10 km, which is small enough to perform stable reference ephemeris estimation. Furthermore, the design error in the reference ephemeris estimation would not grow gradually with the orbital integration time. Therefore, it can be concluded that the progressive approach is more stable and reliable than the direct approach to the orbit initialization when there is an oversized starting velocity error.

Orbit Initialization Experiements
The orbit initialization experiments for comparison of the two approaches will be performed on processing the real-data of KOMPSAT-5 and FY3C satellites. KOMPSAT-5 is the first satellite in Korea that provides 1 m resolution synthetic aperture radar (SAR) images [22]. FY3C is a Chinese meteorological satellite equipped with space-borne GPS/BDS receiver [23]. For the KOMPSAT-5 satellite, the space-borne GPS data is used to compute the discrete kinematic positions by SPP, and the SPP results are used for orbit initialization. For the FY3C satellite, the SPP results based on the space-borne GPS/BDS data are used to estimate initial orbit parameters.
The daily status of the orbit initialization experiments of KOMPSAT-5 on 2016/42~51 and FY3C on 2013/285~294 are summarized in Table 1, where "DA" and "PA" represent the direct approach and progressive approach, respectively. The sign of "×" indicates that the orbit initialization fails, and the sign of "

Analysis and Discussion
Detailed analyses and discussions are carried out for these dates where the failure of the orbit initialization occurs. According to Figure 1, if using the direct approach, the orbit parameters are estimated iteratively only at the third step. At each iteration step, the resulting orbit parameters are used to compute the position at any epoch and the corresponding 3D position errors are computed by using the precise orbits. According to Figure 2, if using the progressive approach, the reference ephemeris is estimated iteratively at the second step, and the orbit parameters are estimated iteratively at steps 4 and 5. Likewise, the resulting reference ephemeris or orbit parameters of each iteration process at steps 2, 4 and 5 are all used to generate the position sequence and the 3D position errors are computed, too. The 3D position error statistics (RMS) for KOMPSAT-5 and FY3C on their respective failure dates are all shown in Table 2 for each iteration of these two approaches, as well as the mean starting position and velocity errors in the XYZ directions, where "2.X", "4.X" and "5.X" represent the X-th iteration at steps 2, 4 and 5, respectively. As can be seen clearly, the starting position errors on these days are all less than 3.0 m, but the velocity errors are up to 1~3 m/s, so the oversized velocity errors have caused the estimation divergence in the direct approach. During the iterative process, the position error in each iteration of the direct approach increases until the orbit parameter estimation diverges, while that in the progressive approach decreases gradually until the final orbit parameters are determined. The table shows the step-by-step fitting of orbit parameters in the progressive approach, and also illustrates its reliability for the orbit initialization. Table 2.
The position error statistics during the iterative process with the direct and progressive approaches. In order to analyze the divergent process in the direct approach, taking the orbit initialization experiment of KOMPSAT-5 on 2016/50 as an example, the position errors in the radial (R), along-track(T) and cross(N) directions and the 3D position error after each iteration step of the orbit parameter estimation are listed in Table 3. As can be observed clearly, the position error is increasing until the estimation diverges in the 6-th iteration, when the estimation iteration process proceeds. After the 5-th iteration, the 3D position error has been more than 50 km, which is a sign that the estimated orbit of KOMPSAT-5 at the altitude of~590 km is severely deviated from the true orbit, so in the 6-th iteration the orbit parameters cannot be estimated. Correspondingly, the estimated altitude of KOMPSAT-5 after each iteration step is shown in Figure 5. As can be seen, at the end of the 5-th iteration, an altitude deviation of more than 50 km is resulted, which causes divergence in the 6-th iteration with a negative altitude estimated. The orbit initialization experiment of KOMPSAT-5 on 2016/50 demonstrates that the direct approach is unreliable when the starting velocity has large errors. Diverge... As for the progressive approach, the iterative estimation of parameters in reference ephemeris at step 2 is the key for the success of orbit initialization. The experiments above have demonstrated the high stability of the reference ephemeris estimation when facing with the starting velocity error of 1~10 m/s. To analyze its estimation process further, considering FY3C on 2013/287 as an example, the position differences in R/T/N directions between the reference orbit obtained at step 2 and the precise orbit are shown in Figure 6. As can be observed clearly, the position error of the reference orbit is within the range of ±1000 m and varies smoothly and continuously. The corresponding altitude also varies smoothly and periodically between 825 km and 845 km, which indicates the estimated reference orbit is sufficiently close to the true orbit. It is based on the smooth, continuous reference orbit that the estimation of orbit parameters X 1 at step 4 could be performed well, although the position error is still at 500 m level. The estimation at step 5 uses the resulting orbit parametersX 1 at step 4 as the starting state and the raw SPP/PPP discrete kinematic positions as measurements. The accuracy of the starting stateX 1 is accurate enough for convergence of the estimation at step 5. All in all, the progressive approach is more reliable then the direct approach for orbit initialization.

Summary and Conclusions
The sufficiently accurate and stable orbit initialization is an important prerequisite for LEO precise orbit determination using GNSS. This paper presents two orbit initialization methods, including the traditional direct approach and the proposed progressive approach. In the direct approach, the orbit parameters, including the position and velocity vector at the reference time and all piecewise dynamic coefficients, are estimated directly and iteratively from the SPP/PPP discrete kinematic positions. The analysis tests on the stability have demonstrated that the orbit parameters estimation is at the risk of divergence, which is mainly caused by the oversized starting velocity error, which may be difficult to prevent completely, because many factors could cause the inaccurate GNSS pseudo-range or Doppler measurements at the reference epoch, including the weak signal, tracking interruption, ionospheric disturbances, abnormal antenna, etc. Hence, the progressive approach is proposed to eliminate or weaken the risk of POD divergence, where the reference ephemeris in form of the parameters in GNSS broadcast ephemeris, instead of the conventional orbit parameters, are estimated from the SPP/PPP discrete positions first. The estimated reference ephemeris is then used to generate a series of tabular positions and velocities with the position error of hundreds or thousands of meters, indicating that the reference orbit is sufficiently close to the true orbit. Then, these tabular positions and velocities can be used as measurements to estimate the orbit parameters by using the direct approach. Finally, the generated orbit parameters are used as the starting state for the next direct estimation, and then more accurate orbit parameters can be determined. The stability analysis tests have validated that, the starting state error, especially the oversized starting velocity error, has less negative effect on the estimated reference ephemeris than on the orbit parameters, so the progressive approach is more stable and reliable than the direct approach.
The orbit initialization experiments are carried out for the KOMPSAT-5 and FY3C missions. The experiment results show that the initialization process fails on some days if using the direct approach, but the initial orbit parameters can be generated reliably on these days using the progressive approach. It indicates that the proposed progressive approach is more reliable and adaptable for LEO precise orbit determination, especially for the automated orbit determination of LEO satellites in a large-scale constellation.
Author Contributions: In this study, X.G. and J.S. designed the algorithm, performed experiments, analyzed data and wrote this paper. F.W. designed the software, performed experiments and edited the manuscript. X.L. supervised its analysis and edited the manuscript. All authors have read and agreed to the published version of the manuscript.