Dynamic GPS-based LEO orbit determination with 1 cm precision using the Bernese GNSS Software

The Astronomical Institute of the University of Bern (AIUB) has been performing GPS-based Precise Orbit Determination (POD) for a large variety of Low Earth Orbit (LEO) satellites since two decades. Traditionally, LEO orbits have been generated by a reduced-dynamic POD strategy using the Bernese GNSS Software, replacing an explicit modeling of non-gravitational forces by dedicated empirical orbit parametrizations. This LEO POD strategy can be advanced by two main developments: on the one hand, use is made of the GNSS Observation-Speciﬁc Bias (OSB) and clock products provided by the Center for Orbit Determination in Europe (CODE), allowing for the resolution of single-receiver GNSS carrier-phase ambiguities. On the other hand, the main focus of this article, a reﬁned satellite non-gravitational force modeling strategy is con-structed to reduce the amount of empirical parameters used to compensate for force modeling deﬁciencies. LEO POD is ﬁrst performed for Sentinel-3, a satellite formation currently consists of two identical satellites -3A and -3B, which experience a similar in-ﬂight environment and allow for direct POD performance comparisons. A third satellite Swarm-C, which ﬂies at a lower altitude and has a more sophisticated surface geometry, is selected to validate the robustness of the new POD strategy. As a result, both the internal consistency checks and external orbit validations suggest superior orbit quality obtained for the three satellites for a time span of 1.5 years (7 June, 2018 to 31 December, 2019). The ambiguity resolution adds strong constraints to the orbits and the satellite non-gravitational force modeling leads to more tightly constrained (towards zero) pseudo-stochastic empirical parameters. The ﬁnal orbit solutions agree with external orbit solutions (cid:58)(cid:58)(cid:58)(cid:58) and independent satellite laser ranging measurements at levels of sub-cm, indicating approximately 20% improvement w.r.t. the nominal reduced-dynamic orbit solutions. This suggests potential beneﬁts to the space geodesy community that always pursues best-possible satellite orbits.


Precise Orbit Determinations in BSW
This section introduces the 6 different POD :::: orbit solutions (Tab.1) that can be generated by BSW. The new features of the proposed non-gravitational force modeling POD strategy will be elaborated subsequently.
Firstly, a kinematic POD strategy is fully independent of LEO satellite force models. A kinematic orbit is an ephemeris at discrete measurement epochs since all positions are determined solely from a high-low satellite-to-satellite geometric positioning. Therefore it requires a minimum number (normally ≥ 5 to guarantee redundancy) of tracked GPS satellites for solving four unknown parameters (3 coordinates and 1 clock offset) of a receiver. The kinematic orbit quality is heavily dependent on the performance of GPS receivers and no solutions are available for epochs experiencing large data outliers or gaps (Yunck, 1996). In BSW, a typical kinematic orbit in the Earth-Centered Inertial (ECI) reference system is related to an epoch-wise trajectory of the antenna phase center position r leo , which is modeled as r leo (t leo ) = R(t leo ) · ( r leo,e,0 (t leo ) + δ r leo,e,ant (t leo )) Phase Center Variation (PCV) that can be created through ground experiments, or currently through a Residual Approach ::::::: residual :::::::: approach using the in-flight GPS data (Jäggi et al., 2009). A typical scientific application of a kinematic trajectory is gravity field recovery from data of non-dedicated gravity missions (??) :::::::::::::::::::::::::::::::::::::::::::::::: (Jäggi et al., 2016;Teixeira da Encarnação et al., 2020).
Secondly, contrary to a pure kinematic orbit, a dynamic orbit is a particular solution fully dependent on the equation of motion and the underlying force models, e.g. the Earth gravity field. A typical representation of a dynamic orbit described in the ECI reference system and its initial conditions can be given by r leo (t leo ) = r leo,0 (t leo ; a, e, i, Ω, ω, u 0 ; Q 1 , ..., Q d ) + δ r leo,ant (t leo ) r = −GM r r 3 + f (t, r,˙ r, Q 1 , ..., Q d ) r(t 0 ) = r(a, e, i, Ω, ω, u 0 ; t 0 ) r(t 0 ) =˙ r(a, e, i, Ω, ω, u 0 ; t 0 ) (2) note that here r leo,0 denotes the LEO CoM coordinate in the ECI frame, a, e, i, Ω, ω, u 0 are the six osculating Keplarian elements of the orbit at t 0 , GM is the gravitational constant times mass of the Earth, Q 1 , ..., Q d indicate d empirical parameters used to compensate for force modeling deficiencies. This satellite trajectory can be described by a particular solution of the equation of motions w.r.t. satellite force models and empirical parameters Q, e.g. co-estimated parameters that are scaling dynamic force models. It is difficult to determine an ideal dynamic orbit for a LEO satellite which orbits the Earth in such a heavily perturbed environment that perturbations might vary significantly even in a short orbit arc.
number of Piece-wise Constant Accelerations (PCAs) to ensure that a satellite trajectory is continuous and differentiable at any epoch. PCAs can be first characterized by a priori known statistical properties, e.g. a priori variances σ 2 p and spacing time ∆ t . The equation of motion of a reduced-dynamic orbit in the ECI reference frame can be represented bÿ r = −GM r r 3 + f (t, r,˙ r, Q 1 , ..., Q d , P 1 , ..., P s ) where, compared with Eq. 2 and given the same initial conditions, Q 1 , ..., Q d are often set as periodic and/or once-per-arc constant accelerations in three directions defined by the local orbital :::::::: reference : frame (i.e. radial, along-track and cross-track directions). P 1 , ..., P s are the s pseudo-stochastic parameters to compensate for force modeling deficiencies. BSW was traditionally used without explicit modeling of non-gravitational forces, a reduced-dynamic POD approach which was very successful to generate LEO orbit solutions of high quality, e.g. (Jäggi et al., 2007;Bock et al., 2011). It required, however, relatively loose constraints to fully compensate for the not explicitly modeled non-gravitational forces with PCAs, which are usually set up over intervals ranging from about 5 to 15 minutes (Jäggi et al., 2006). Thanks to the new non-gravitational force modeling capabilities, the uncertainties of satellite dynamics are significantly reduced, which may lead to more tightly constrained (towards zero) pseudostochastic parameters. Especially when it comes to a reliable radial leveling e.g.
for altimetry satellites, the use of empirical and pseudo-stochastic parameters should be carefully revised/limited, since they will allow to degrade the orbit if offset problems of any kinds (e.g. PCO and CoM) exist.
Therefore in our research, advances are made in the explicit modeling of non-gravitational forces and reducing the heavy dependence on empirical parameters, suggesting a more dynamic orbit solution. It will be shown that onceper-arc constant accelerations can be removed and PCAs can be more more tightly constrained towards zero, i.e. the a priori standard deviation (STD) σ p is reduced by a factor of 10 specifically for the associated satellites in this article.
8 An overview of the new satellite dynamic modeling and POD processing strategy is summarized in Tab.2. Details for each non-gravitational force modeling will be elaborated in Sect. 3. The a priori STD set up of PCAs for the Swarm-C satellite will be doubled as compared to the Sentinel-3 satellites due to stronger perturbations. A kine-

Non-gravitational Force Modeling
A LEO satellite normally experiences more complex perturbations than satellites at higher altitudes. The main perturbations are gravitational forces that can be numerically computed based on various supporting models (Tab.2).
The modeling of non-gravitational forces is more challenging since it often relies on various external products which unfortunately can not perfectly represent the real in-flight perturbation environment (Doornbos, 2012). Therefore more empirical parameters are necessitated to address the model imperfections, and often the co-estimation of dimensionless scale factors during a POD process is conducted to absorb deficiencies of modeled forces and satellite macro-model.
In this article, only SRP and AF are scaled, the ERP is not scaled since it will impact the orbit leveling particularly in radial direction (Montenbruck et al., 2018a;Hackel, 2019). The overall non-gravitational forces can be given by where SRP, the Earth REFlectivity radiation pressure (REF), the Earth EMis-siviTy radiation pressure (EMT) and AF are the surface forces considered. As stated above, a description of LEO satellites in terms of flat-plate macro-models is widely accepted for non-gravitational force modeling. This article uses the Sentinel-3 8-plate macro-model introduced in (Fernández, 2019a), which has been widely used by the CPOD and QWG community. The Swarm-C satellite has a more lengthy complex geometry, which can be described by a 15-plate macro-model (Montenbruck et al., 2018b;Hackel, 2019 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 3.1. Solar Radiation Pressure SRP originates from the interaction between photons and satellite surface materials. In general it can be divided into three categories -absorption, specular reflection and diffuse reflection -which are determined by the characteristics of surface materials. SRP is causing an acceleration given by where the index i enumerates all plates of the macro-model, the satellite mass m can be extracted from an official mass table, P 1AU is the solar radiation pressure at the distance of 1 AU (astronomical unit), r Sun,s is the distance between the Sun and the satellite, and f s denotes the so-called geometric shadowing factor which takes into account a few impacting factors e.g. the proportions of Sun radiation absorbed by the atmosphere, eclipse, satellite and Earth shadowing, etc. (Doornbos, 2012;Hackel et al., 2017). The focus of computing SRP is modeling C S,i , a vectorial radiation pressure coefficient for a certain plate i. It is given by Doornbos (2012); Montenbruck and Gill (2012) (note Eq. 6 removes i for the sake of readability) where, the fractions of diffuse reflection (p d ), specular reflection (p s ), and absorption (p a ) of photons for the short-wavelength visible radiation are described by the macro-model and they sum to 1. n denotes the unit vector of a plate's surface normal that usually points outwards, A is the surface area of the plate, arrays generate power, the accumulation of energy will increase the temperatures of particularly non-solar plates. Normally, spacecraft engineers are aiming for a thermal balance and protect these plates with special materials such as polyimide. Therefore, we assume that all the absorbed photons will be spontaneously re-emitted according to Lambert's cosine law, as marked by the blue term of Eq. 6, according to the formula it is then exactly the same with C S,d   orbit arc the Sentinel-3A satellite's beta angle is 23.2 • and that for the Swarm-C satellite is −28.5 • . In fact, the Sentinel-3 formation has quite stable beta angles only ranging between 23.2 • and 34.5 • during the 1.5 years due to its Sun-synchronous orbit (inclination 98.6 • ), whereas the Swarm-C satellite's beta angles vary from −79.4 • to 81.2 • due to its more polar orbit (inclination 87.4 • ).
Variations of the beta angles are visible in Fig.9.

Earth Radiation Pressure
The Earth's energy budget accounts for a balanced situation where the Earth reflects and emits nearly all incoming solar radiation back into the outer space.
ERP is caused by 1) scattered short-wavelength visible solar radiation and 2) emitted long-wavelength thermal infrared radiation of the Earth. A few scientific satellites such as NASA's Aqua and Terra, have been continuously measuring the radiosity of the Earth. In BSW the monthly Clouds and the Earth's Radiant Energy System (CERES) S4 grid products, obtained from the Aqua and Terra satellites with a spatial resolution of 2.5 • × 2.5 • , are used to compute the corresponding ratios of radiative flux to the incoming solar irradiance, which is set to 1372 W/m 2 at 1 AU (Wielicki et al., 1996). In the work done by Hackel et al. (2017); Montenbruck et al. (2018b), the CERES-S4 grids were represented by zonal coefficients of Legendre polynomials to compromise between computation efficiency and modeling accuracy. Please note that there are different types of CERES associated products, e.g. Vielberg and Kusche (2020) made use of the CERES hourly SYN1deg data for a more refined modeling ERP.
All possible monthly CERES-S4 data from July 2002 to September 2019 are retrieved to generate the reflectivity and emissivity grids. The variations of the mean of monthly grids are displayed in Fig.2. It reveals that both reflectivity and emissivity grids change significantly from month to month, nevertheless they are rather stable for the same month from year to year, within a difference of merely 0.01. It is also interesting to see that the reflectivity has been decreasing during the past years. During the selected period, there are four months (August 2002, January 2011, December 2011, and August 2014 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62 63 64 data gaps. The CERES-S4 products also have a latency of a few months, which might obstruct a near-real-time POD processing, e.g., for the Sentinel RSR (Fernández, 2019a) ::::::::::::::::::::: (Fernández et al., 2019). To overcome these deficiencies, we average the monthly products from all available years. The averaged grid for June, and the difference between the averaged and the specific June 2018 grid, are depicted in Fig.3. The seasonal reflectivity and emissivity changes are visible particularly for the high-latitude regions during the polar nights and days. In addition, an averaged monthly grid can not fully describe the dynamic variation within a month, therefore a linear interpolation is performed between the current monthly grid and its neighboring monthly grid, which is selected as either the previous month or the next month depending on the day of month :: for ::::: orbit ::::::::::: computation. When comparing with the ERP modeling for the Sentinel-3A satellite using the specific monthly grid e.g. for June 2018, the modeled forces using the averaged June grid and performing the additional linear interpolation 16   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 between two months (May and June) show discrepancies at levels within 1% (selected period: 7 June, 2018, 00:00-03:20, consistent with Fig.4), which can be easily handled by PCAs and do not impact the reduced-dynamic LEO POD solutions at a visible manner. Therefore we use the averaged products for orbit computations. The total ERP can be described as where j is the index of a grid with N bins, N = 72 × 144 for a resolution of 2.5 • × 2.5 • . The computations of P REF and P EM T need to modify a few aspects based on Eq. 5: firstly, all radiations originate from the top of the Earth's atmosphere (ToA, 30 km), rather than from the Sun. Secondly, P REF depends on the illumination status of the Earth, whereas P EM T does not. Thirdly, for the computation of C R,i , Eq. 6 can be used when the Sun-satellite vector is replaced by the ToA element-satellite vector. In addition, the computation of C E,i has to specifically use the material characteristics (again absorption, specular reflection and diffuse reflection) for the long-wavelength infrared radiation. Fig. 4 shows that ERP has the largest component in radial direction. The scale factors for ERP are not estimated (fixed to 1) otherwise potential erroneous PCO or CoM offsets will turn into radial orbit shifts, which might be problematic for particularly altimetry missions e.g. Sentinel-3.

Aerodynamic Force
The thermosphere consists of neutral atoms and charged particles that are interacting with the satellite surfaces. Two component forces can be distinguished by definition, drag is the projection of AF onto the velocity direction with respect to the atmosphere and lift is the portion of AF perpendicular to the velocity direction. Nevertheless lift normally only accounts for a small proportion, e.g. for the Swarm-C satellite it is at a level of 1 nm/s 2 , comparing to its drag at a level of 50 nm/s 2 . AF is heavily dependent on the orbit altitude and the dominating force for LEO satellites flying at very low altitude, e.g. 17   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  the CHAMP, GOCE, GRACE and Swarm satellites (Visser et al., 2009;Mao, 2019a). AF can be described as where C D,L denotes the coefficients for drag and lift, which can be modeled by algorithms such as Goodman, Sentman and SESAM (Doornbos, 2012;Girardin, 2016;Pilinski et al., 2013); v i represents the relative velocity between a satellite plate and the atmosphere. The upper thermosphere can be very dynamic such that the Horizontal Wind Models (HWM), which reflect the time-varying atmospheric circulation dynamics, are needed to compute a most realistic relative velocity (Drob et al., 2015). As for the Swarm-C satellite, this contributes roughly a few nm/s 2 to AF particularly in cross-track direction for polar orbits.
A precise modeling of AF necessitates high-precision atmospheric models.   spheric density models is less than 10% during the selected period (Bruinsma, 2015), and no significant impacts on the POD performances are witnessed in BSW when estimating scale factors and PCAs.

Results and Discussion
This section first includes a quality assessment of the associated GPS data, followed by internal consistency checks and external validations to different orbit solutions.

Summary and Outlook
This article investigates the latest development of the Bernese GNSS Software for low Earth orbit satellite precise orbit determination. The focus of the new strategy is based on a refined dynamic modeling of satellites. Three main non-gravitational forces -solar radiation pressure, Earth radiation pressure and aerodynamic force -are modeled as the sum of independent plate-wise surface force acting on a satellite macro-model. The non-solar plates are assumed to spontaneously re-emit all absorbed energy after receiving radiation pressures from both the Sun and the Earth. The modeling of Earth radiation pressure further takes into account of the Earth's Radiant Energy System S4 grid products, and a set of monthly products are created by averaging data from the past 19 years. A linear interpolation has to be done between two neighboring months.
Besides, the GNSS Observation-Specific Bias (OSB) and ambiguity fixed clock products provided by the Center for Orbit Determination in Europe (CODE) allow for single-receiver integer ambiguity resolution which further constrains the associated orbit solutions.
GPS data of three Earth observation satellites -Sentinel-3A, Sentinel-3B and Swarm-C -are processed for a time span of 1.5 years between June 2018 and December 2019. Firstly, the empirical accelerations, that are supposed to compensate for deficiencies in the satellite dynamic modeling, can be partially replaced by the non-gravitational force modeling. For instance, statistics of the mean of three dimensional (radial, along-track and cross-track ::::::: referring :: to :::: the :::: local ::::::: orbital :::::::: reference :::::: frame) empirical accelerations for the Sentinel-3A satellite can be significantly reduced from −38.6/ − 1.4/ − 31.0 nm/s 2 to below 1 nm/s 2 in the three directions. The standard-deviations are also reduced to much lower levels due to tighter constraints (by a factor of 10, Tab.1)and : , ::: and ::::: more :::::::::::: importantly, the non-gravitational force modeling. For the situations where the supporting models for non-gravitational forces can not fully represent the real in-flight environment, the co-estimated corresponding scale factors will adjust the associated forces, e.g. the over-performed aerodynamic forces mod- 37   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 eling can be adjusted by scale factors smaller than 1. This is consistent with LEO POD performance using the GHOST software package ( Van den IJssel et al., 2020). The modeling of solar radiation pressure has been very stable and no extra scale factors are estimated for the Earth radiation pressure to maintain its neutral contributions to shift orbits particularly in radial direction. In addition, the integer ambiguity fixing brings more tightly constrained geometry to the satellite positions.
Finally, the best possible orbits can be obtained by combining all these benefits. The satellite laser ranging validations show orbit precisions of 9.2, 9.2 and 9.0 ::: 8.8 mm for the Sentinel-3A, Sentinel-3B and Swarm-C satellites, respectively. The different comparisons again prove a superior performance of the new POD strategy, however also indicate a few small issues to the accuracy of the coordinates of the satellite Center of Mass (CoM) ::::: center ::: of ::::: mass : and instrument reference points.

Acknowledgment
This research was funded by the Swiss National Science Foundation (Grant