Abstract
Real-time precise global navigation satellite system (GNSS) orbit and clock products play a key role for real-time GNSS-based applications, both in the scientific and industrial communities. Different from the typical two-step procedure to generate orbit and clock solutions separately, we estimate the real-time orbit and clock products simultaneously using a Kalman filter. For this purpose, we developed a GNSS data processing software that can run in pseudo-real-time mode with RINEX files and is ready to run in real-time mode once given the real-time observation stream. Meanwhile, a quasi-orbit-fixed solar radiation pressure (SRP) model is developed. In order to verify the performance of the software and the new SRP model, several experiments with a global network of 60 tracking stations over a time span of three months were conducted to generate real-time Global Positioning System (GPS) orbit and clock products. Then, the results were assessed in terms of accuracy and efficiency, both critical for real-time precise GNSS applications. Compared to the International GNSS Service (IGS) final orbits, the real-time GPS orbit accuracy was 2.82 cm, 5.45 cm, and 5.47 cm in the radial, along-track, and cross-track components, respectively. The precision of the clock product in terms of standard deviation (STD) value was about 0.1 ns. Moreover, the average execution time per epoch was usually less than 1.0 s, which ensures the high efficiency of the processing.
1 INTRODUCTION
Real-time precise point positioning (PPP) and ambiguity resolution like the PPP real-time kinematic (RTK) technique (Wuebbena et al., 2005; Zumberge et al., 1997) have gained a wide range of usage in the field of real-time GNSS (global navigation satellite system) applications, especially in the natural and anthropogenic hazard monitoring and warning systems that are deployed to mitigate the risk of various disasters, such as the earthquakes, tsunamis, landslides, and dam-breaks, to name a few (Bock & Melgar, 2016). To fulfill these life- and treasure-critical engineering projects, the real-time generation of high-precision GNSS orbit and clock products is indispensable, as they provide the space-time reference frame for these applications. Thus, orbit determination and clock estimation for GNSS satellites in real time is currently a crucial research focus in the GNSS community.
There are mainly two types of institutions that are interested in computing real-time precise GNSS orbit and clock products. The first type includes commercial companies who are willing to broadcast the product stream via the communication satellites or the Internet to serve their consumers for profit (Dixon, 2006; Leandro et al., 2011). The other type consists of GNSS-related institutes and universities or their voluntary associations, among which the International GNSS Service (IGS) is the most famous. They provide the product mainly for the purpose of scientific research with public benefit as a by-product (Caissy et al., 2012; Kuang et al., 2019; Laurichesse et al., 2013).
As for the strategy to generate real-time orbit and clock products, there mainly exist two approaches. The first approach is to get the real-time orbit information directly from the predicted orbit, and then to estimate the real-time clock product using the filtering technique (Caissy et al., 2012; Dixon, 2006; Leandro et al., 2011). In contrast, the second approach estimates the real-time orbit and clock products simultaneously within a Kalman filter (Kuang et al., 2019; Laurichesse et al., 2013). Obviously, both approaches have their advantages and disadvantages.
Because of the availability of high-precision predicted orbit, the first approach has the advantage of low complexity and high efficiency when estimating the clock product in the filter step. However, we dismiss this approach due to the following drawbacks. First, if the predicted orbit is retrieved from an external agency, this situation will seriously weaken the independence of the service provider. Second, if the predicted orbit is self-generated, the whole processing system can be broadly divided into two parts: one responsible for the orbit prediction using a batch estimator and the other can, thus, only estimate the clock corrections with a sequential estimator. This strategy has both fragility and redundancy. On the issue of fragility, the service provider must ensure that both the batch estimator and the sequential estimator are operational without disruptions. In terms of redundancy, the service provider must process the data from two station networks or process the data of one network twice.
In principle, the second approach is the most natural and rigorous method to generate real-time GNSS orbit and clock products. However, the increased computational complexity when including the orbit parameters in the Kalman filter becomes a major concern for service providers, prohibiting widespread use of the optimal approach. This situation is the motivation of our research.
With the availability of high-precision geopotential models and planetary ephemerides, the conservative forces acting on GNSS satellites are no longer problematic for precise orbit determination (POD). Due to the orbit height being about 20,000 kilometers above the Earth’s surface, the effect of atmospheric drag on GNSS satellites can be totally ignored. The perturbing effect from the solar radiation pressure (SRP) becomes the largest error source in GNSS orbit modeling. With a magnitude of about 100 nm/s2, this perturbing acceleration is a crucial factor impacting orbit quality. Due to the complexity and variation of the solar wind (King-Hele, 1983), the SRP perturbation is very difficult to model. The large solar panels of GNSS satellites and the imperfect attitude control system onboard make the situation even worse.
From the very beginning, SRP modeling for GPS satellites was conducted by the satellite manufacturer using the so-called analytical method (Fliegel & Gallini, 1996; Fliegel et al., 1992). For this method to work, the optical and thermal parameters for each satellite surface component must be measured before the satellite launch. The SRP acceleration for an orbiting satellite can then be computed according to the illumination condition of each satellite component. The physical meaning of the analytical method is very clear; a complete modeling of SRP perturbation using this method is possible but very difficult (Ziebart et al., 2002). Another method to deal with the SRP perturbation is to set up several dynamical parameters empirically according to the characteristics of the SRP effect. The Empirical CODE Orbit Model (ECOM) developed at the Center for Orbit Determination in Europe (CODE) IGS analysis center are representative work of the empirical method (Arnold et al., 2015; Beutler et al., 1994; Springer et al., 1999).
ECOM models are widely used in IGS analysis centers and have demonstrated excellent performance, so we have adopted the reduced ECOM model (Springer et al., 1999) as the default in our real-time data analysis. Meanwhile, a new empirical SRP model was developed for real-time precise orbit determination. The frame of this new SRP model is quasi-orbit-fixed (QOF) rather than spacecraft-fixed like the ECOM model series. After a long data analysis, we assess the validity of this new empirical SRP model.
In this article, we will briefly describe the methodology of our research work, such as the error sources that must be considered in GNSS measurements and their treatments in measurement modeling, the various perturbation forces acting on GNSS satellites and their modelings, and parameter configurations in preprocessing and filtering stages. All of these are of great significance for high-precision GNSS orbit and clock estimation. The modeling work of the new SRP model is also a focus. Starting from some reasonable assumptions, we derive a formula-based empirical SRP model.
This article is mainly structured into five sections. After the general introduction in Section 1, we start with a review of the methodology underlying the research work in Section 2, where the new quasi-orbit-fixed SRP model for GNSS satellites is presented. The processing strategies used for the generation of real-time GNSS orbit and clock products are described and tabulated in Section 3. Then, in Section 4, some experiments taking the GPS constellation, for instance, are illustrated in detail, followed by the resulting assessment of orbit and clock results in terms of accuracy and efficiency. Finally, the research work is concluded with Section 5.
2 METHODOLOGY
Both the observation equation and equation of motion are fundamental to the estimation procedure involving the POD problem. The observation equation describes the relationship between the system state parameters and can constrain the orbit information geometrically. On the other side, the equation of motion represents the force field surrounding a satellite and dynamically constrains the satellite’s trajectory. Thus, the fundamental principles of both parts deserve a brief description.
In the application of real-time data analysis, Kalman filters are extensively used. In our software, a variant of the Kalman filter using the matrix factorization technique, the U-D filter, is employed due to its good numerical stability and high efficiency. See Bierman (1977) and Gibbs (2011) for more information about the U-D filter. The Estimation Subroutine Library (ESL) in close connection with Bierman (1977) can now be accessed at https://netlib.org/a/esl.tgz.
2.1 Observation Equation
To eliminate the first-order ionospheric path delay, an ionosphere-free (IF) linear combination is used in our processing. Thus, we need dual-frequency pseudorange and carrier-phase measurements for each receiver-satellite link. The corresponding observation equations for these raw measurements can be formulated as:
1
where 
is the composition of various frequencyindependent terms. The various symbols used in Equation (1) are explained as follows:
s: satellite index
r: receiver index
ρ: geometrical distance from the center of mass of the satellite to the antenna reference point of the receiver with relativistic effects and station tidal displacements considered
c: speed of light in vacuum
dt: receiver/satellite clock error
T: tropospheric path delay
A: receiver/satellite antenna phase center correction
I: ionospheric path delay
λ: wavelength of the carrier radio wave
N: initial carrier-phase ambiguity
W: wind-up effect
ϵ: measurement noise and unmodeled errors
The IF linear combination is also extensively used for the generation of various IGS products. So IGS conventions define that the satellite clock must be referred to the IF linear combination of basic P1 and P2 observations, as is the case with the satellite clock definition in broadcast ephemerides. If a receiver cannot provide P1 and P2 measurements concurrently, the consistency problem due to signal hardware delay should be dealt with and differential code bias (DCB) correction(s) must be applied. To minimize the dependence on external resources as far as possible, however, we don’t use the P1C1 DCB product from the IGS. Instead, we use pseudorange measurement C1 rather than P1, since all receivers can provide this observation type. Although the P1 observable has a lower measurement noise than C1 observables, considering the small weight of pseudorange measurements compared with the carrier-phase observation, this adoption will mainly change the definition of clock indication with little influence on the precision of the clock product.
2.2 Equation of Motion
The equation of motion for a GNSS satellite can be written as (Dach et al., 2015):
2
with:
G: gravitational constant
M⊕: the Earth’s mass
: satellite’s position and its time derivativesa: perturbing acceleration acting on the satellite
p: dynamical parameters to be estimated.
The variables r, 
, and 
in Equation (2) are all time dependent. For brevity and readability, however, the time argument is omitted. Besides, the composite perturbation term a can be further detailed as:
3
with:
ageo: the Earth’s non-spherical gravity field
atid: tide-related potentials
a3rd: gravitational attractions from Moon/Sun/planets
arel: relativistic correction
asrp: solar radiation pressure
Given the initial state of the satellite (position r0 and velocity 
) at a specific epoch t0, a typical initial value problem (IVP) arises:
4
With the help of a numerical integrator, the IVP of Equation (4) can be easily solved. The resulting approximate reference orbit can then be used for the linearization of observation equations.
For the problem of dynamic POD, solving another set of differential equations (i.e., variational equations) is vitally important. Closely associated with each equation of motion, there exists one, and only one, set of variational equations (Montenbruck & Gill, 2000). The solutions of variational equations are the so-called state transition matrix and sensitivity matrix, which provide the essential information for the correction of the not-so-correct initial condition r0 and 
and specific dynamical parameters p. This process of orbit improvement based on the differential correction method is the very meaning of precise orbit determination. Because of the close connection between these two sets of differential equations, it is necessary to solve the variational equations simultaneously with the equation of motion (Beutler et al., 2005).
For more information about the theory of dynamic POD, see Montenbruck and Gill (2000), Beutler et al. (2005), and Tapley et al. (2004).
2.3 Quasi-Orbit-Fixed SRP Model
First of all, let us make several reasonable assumptions to simplify the problem of SRP modeling:
The solar flux at one point is proportional to the inverse square of the distance from the point to the Sun considering that the Sun radiates energy spherically.
The solar radiation pressure is a parallel force field near the Earth due to the fact that the orbit size of GNSS satellites is several orders of magnitude smaller than one astronomical unit.
The solar panels of GNSS satellites are always perpendicular to the assumed SRP force field.
The mass and illuminated cross-section area of GNSS satellites remain constant.
Each GNSS satellite has a circular orbit.
Based on these assumptions, the direction of the SRP acceleration is always from the Sun to the Earth and, thus, only magnitude is considered here. The magnitude of the SRP acceleration of a satellite can be expressed as:
5
with the explanation of symbols as follows:
A: satellite illuminated cross-section area
m: satellite mass
r⊙,⊕: distance from the Sun to the Earth
r⊙,s: projection of the distance from the Sun to the satellite in the direction from the Sun to the Earth
P⊙,⊕: intensity of SRP at a distance of r⊙,⊕
CR: an adjustable coefficient
Csrp: another adjustable coefficient
From Equation (5), we can conclude that the SRP acceleration of a satellite is a periodic function whose period coincides with the satellite’s orbital period. With the help of the relation:
6
the periodic component in the SRP acceleration of the satellite can be revealed by:
7
where rs is the radius of the satellite circular orbit, β is the elevation of the Sun with regard to the satellite orbital plane, and Δu is the difference between the argument of satellite’s latitude and the argument of the Sun’s latitude in the satellite orbital plane. Considering that rs is orders of magnitude smaller than r⊙,⊕, only the main terms in both the numerator and denominator can be reserved. Thus, the periodic component in Equation (7) can be simplified as:
8
From Equation (8), we know that the periodic component of the SRP acceleration is a trigonometric function with Δu as argument, considering that the change of other terms is so slow compared to the change rate of Δu. This is the underlying reason that trigonometric functions are used in the modeling work of SRP perturbations (Arnold et al., 2015; Beutler et al., 1994; Springer et al., 1999). Thus, in the D direction from the Earth to the Sun, the SRP acceleration can be modeled as:
9
where φD is used to compensate for the potential lag of the satellite solar panels when tracking the Sun. In addition to the component in the D direction, we further assume that, in reality, there are some projections or leakages of the SRP acceleration in other directions, and these components are as relatively stable as the D component in space. By taking into account the symmetry of the geometric relationship between the Sun, the Earth, and the satellite orbit, we define the other two directions Y and X (see Figure 1). Like the D component, the components of SRP acceleration in the Y and X directions are modeled as:
10
Geometric relationship between Sun-Earth-satellite triplet and orientation of SRP reference frame
After conducting POD experiments based on the above prototype SRP model, we found that the best orbit accuracy could be achieved with the following model:
11
Due to the fact that the distance variation from the Sun to the satellite is also affected by the Sun’s elevation with regard to the satellite orbital plane, we can consider this effect by incorporating β into the D component of the SRP model. According to Equation (8), the revised model becomes:
12
Until now, we have assumed that each GNSS satellite has a circular orbit. If a GNSS satellite is injected in orbit correctly, the circular orbit is a pretty good approximation. However, in the case of launch accidents like Galileo satellites E201 and E202, the satellite orbit might have a large eccentricity. We can also consider the orbit eccentricity into our SRP model. Because the cos β cos Δu and cos β sin Δu can be substituted by (rs · eD)/ rs and cos β · (rs · eX)/ rs in the case of circular orbit, the generalized model can be written as:
13
where a is the semi-major axis of the satellite orbit, rs is the satellite’s geocentric position vector, and eD and eX are unit vectors in the D and X directions, respectively. The inclusion of both the Sun’s elevation β and the satellite orbit eccentricity only leads to a revision of SRP model in the D direction, since this is the primary direction that is formula based.
The above content is how the new SRP model is conceived, and in the following is a formal definition of the reference frame in which the new SRP model is represented. The three orthogonal directions are defined as:
14
where rs and vs are geocentric position and velocity vectors of the satellite and r⊙ is the geocentric position vector of the Sun, all in one inertial reference system. Thus, eh is the normal unit vector of the satellite orbital plane, eD is the unit vector pointing from the Earth to the Sun, and eX is in the intersecting line of the Earth’s terminator plane and the satellite orbital plane, whose direction is consistent with the requirement that eY be on the same side of the orbital plane as the vector eh. Figure 1 displays the geometry of the Sun-Earth-satellite system and the orientation of the newly defined SRP reference frame. This orthogonal reference frame is rather different from that of ECOM models. While the reference frame for ECOM models rotates constantly together with the satellite attitude (Montenbruck et al., 2015), our reference frame is nearly fixed to the satellite orbit when the satellite revolves around the Earth.
There is an issue with the definition of the quasi-orbit-fixed reference frame for the GLONASS constellation. Because the GLONASS satellite orbit has an inclination of about 65° and the obliquity of the ecliptic is about 23.5°, the maximum value of the Sun’s elevation with regard to the GLONASS satellite orbital plane can be up to 90°. When the Sun is perpendicular to the satellite orbital plane, the definition of the Y and X directions is singular. We will address this problem of singularity later when the GLONASS orbit and clock products are estimated. Due to the smaller orbit inclinations of the BeiDou, Galileo, and GPS satellite systems, this is not a problem for them.
3 PROCESSING STRATEGY
Whether from measurement modeling to orbit modeling or from preprocessing to filtering, there exist many options for the configuration of processing strategies. To make the generation of precise GNSS orbit and clock products possible, it is necessary to illustrate these key pieces of information in some detail.
3.1 Measurement Modeling
The carrier-phase observable can be measured by a geodetic receiver with a precision at the centimeter, or even millimeter, level. In order to fully exploit the potential of these high-precision measurements to obtain an accurate network solution, the error sources larger than 1 cm residing in the measurements must be modeled carefully, either for compensation in data reduction or for estimation in data filtering. An overview of the various error sources and their treatments is listed in Table 1.
Error Sources in Measurement Modeling
Through extensive studies by many researchers over the past several decades, the error sources in GNSS measurements and their modelings have become standardized to a large extent. Take the following three points, for example. First, the correction formulae for tidal and loading effects on station displacements are given in the International Earth Rotation and Reference Systems Service (IERS) conventions (Petit & Luzum, 2010). Second, for the GPS constellation, the antenna phase center corrections on both the satellite end and receiver end rely on the official IGS antenna phase center model contained in the file igsyy_wwww.atx (Rothacher & Schmid, 2010). Third, the modeling work of troposphere hydrostatic a-priori delay and mapping functions is also summarized in the IERS conventions. The models in our research work are aligned with the latest IERS conventions and IGS standards as much as possible.
In addition to the error sources in Table 1, multipath is also a major error source for GNSS measurements. It is closely related to the receiver operating environment and, thus, the impact is very hard to model for compensation or estimation. So, the multipath is merged into receiver noises and is directly neglected in data processing.
3.2 Orbit Modeling
For the high-level accuracy of a few centimeters required in GNSS precise orbit determination, the forces acting on the satellites should be modeled carefully. Thanks to the high orbit of GNSS satellites, the force models for them are not so complicated compared to the satellite orbits in an atmosphere (King-Hele, 1964). The most important force models used for orbit integration in our software are listed in Table 2. Some other perturbing accelerations are not considered at all because of their much smaller magnitude and modeling difficulty, like Earth’s albedo, antenna thrust, and thermal re-emission.
Processing Strategies in Orbit Modeling
With the work of force modeling completed, the numerical integrator also plays a significant role in orbit integration. For the real-time orbit estimation based on the filtering technique, the work of orbit integration in the phase of state propagation can be conveniently done using a single-step integrator. The classical Runge-Kutta fourth-order numerical integrator with a step size of 30 s is deployed in our software.
3.3 Quality Control
Quality control (QC) for GNSS measurements is indispensable for high-precision orbit and clock estimation in which cycle-slip detection is the major task. Because real-time preprocessing can only rely on the data of zero age, the QC in such a case, thus, becomes much more difficult. Nevertheless, we have taken three measures to deal with the QC problem in real-time mode.
First, the elevation cut-off angle is set to 15 degrees to exclude low-elevation data that is more likely contaminated by cycle slips due to the low level of signal strength and the high probability of ionospheric scintillation. Then, the classical method for cycle-slip detection based on Melbourne-Wuebbena (MW) and geometry-free (GF) linear combinations is performed (Blewitt, 1990). Because of the large noise of pseudorange measurements and the unpredictable variation of ionospheric delays, however, the small cycle slips cannot be detected reliably in this step. Thus, detecting large cycle slips (> 5 cycles) is the primary task at this moment. In the third and last step, the screening of post-fit residuals is conducted to remove the inconsistencies after filtering all the network data of the current epoch. The threshold values used in the last two steps are given in Table 4.
Filter Parameterization and Tuning
Processing Configurations in Experiments
For cycle slips detected in the second step and anomaly carrier-phase measurements screened in the third step, the corresponding ambiguities are re-initialized thereafter. Through the layer-by-layer QC procedures, it is hoped that network data processing becomes immune to cycle slips.
3.4 Filter Settings
As stated earlier in Section 2, the observation equation describes the relationship between the system state parameters. However, it is not possible to estimate all parameters at the same time. To remove the singularity, some parameters must be considered as known so that the resulting solution can be tied to a specific reference frame and become definite and meaningful. In our case, the reference stations are considered as fiducial points and their coordinates are fixed to the SINEX solution. Meanwhile, the Earth rotation parameters are fixed to the C04 product.
With the reference frame specified, other state parameters can be set up and estimated in the filter. The initial noise and the process noise are key features of the Kalman filter and play a significant role in the filter solution’s accuracy and stability, so the filter should be tuned carefully. The value choice for the initial noise has a larger degree of freedom and, thus, it is not a tough problem. However, this is not the case with process noise. To account for the statistical information distortion caused by process model uncertainty with the passage of time, process noise is applied when performing the Kalman filter time update procedure (Francisco, 1996). For this reason, the selection of process noise values must be consistent with the reality and can be derived through a great deal of experiments.
Take SRP parameters, for example. Due to the differences between the new SRP model and the reduced ECOM model, such as the reference frame for SRP modeling and the number of SRP parameters, the process noise for them is not the same. To determine the optimal process noise for each SRP model, different values are tested using 14-day data analysis. The daily accuracies of these real-time orbits are displayed in Figure 2. Based on the criterion of continuous generation of stable high-precision orbit products, the process noise for the new SRP model is set to 10−13 m/s2 and the process noise is set to 10−12 m/s2 for the parameters of the reduced ECOM model. The filter parameterization and tuning in our data processing are summarized in Table 3.
Determination of optimal process noises for SRP parameters; all process noise values are in units of m/s2.
4 EXPERIMENTS
Based on the fundamental principles described thus far, a software prototype for real-time GNSS orbit and clock estimation was developed. To demonstrate the performance of the software in terms of accuracy and efficiency and to assess the contribution of the newly developed SRP model to orbit accuracy, some experiments generating real-time orbit and clock products for the GPS constellation were conducted. A global network consisting of 60 IGS reference stations was selected and the time span of the experiment data was 3 months (the fourth quarter of the year 2019). Figure 3 shows the global distribution of the tracking stations.
Distribution of the 60-station network for experiments
The GNSS measurements stored in the RINEX files are processed epoch by epoch to simulate real-time data processing. For comparison with IGS final products, the generated orbit and clock products were stored in SP3-c fromat and the final orbit and clock products in SP3 format from the IGS analysis center, GeoForschungsZentrum (GFZ), were considered to be the benchmark. It should be noted that the so-called product accuracy in the following sections is just a reasonable approximation because there was no orbit and clock truth available. For the reason of orbit convergence, data processing began the day before the experiment time span began.
To evaluate the performance of the new SRP model, three experiments were conducted. The only difference between these experiments was the use of different SRP models. According to the SRP models used, the three experiments were labeled ecom (the reduced ECOM model in Springer et al. [1999]), circ (the model in Equation [11]), and beta (the model in Equation [12]), respectively. Because of the small orbit eccentricity for all GPS satellites, we didn’t experiment with the generalized SRP model in Equation (13). Based on the tests in Section 3.4, the process noise for the reduced ECOM model was set to 10−12 m/s2 and the process noise was set to 10−13 m/s2 for the SRP parameters in the other two experiments. Table 4 lists the processing configurations for the experiments.
4.1 Accuracy
With the GFZ final products as the benchmark, the real-time generated GPS orbit and clock products were assessed day by day. The satellites suffering from maneuvers or clock problems were excluded from the comparison. Considering the GPS constellation as a whole, the daily orbit accuracies given by root-mean-square (RMS) values in the radial, along-track, and cross-track (RTN) directions of the orbital frame for these three experiments are illustrated in Figure 4. For the experiment with the reduced ECOM model, the mean orbit accuracy for the whole GPS constellation during the experiment time period was 4.30 cm, 6.34 cm, and 5.57 cm in radial, along-track, and cross-track directions, respectively. For the experiment using the new SRP model from Equation (11), the mean RMS values of the orbit product were 2.82 cm, 5.45 cm, and 5.47 cm in the three directions, respectively. When further considering the effect of the change of β angle, the orbit mean RMS values were 2.80 cm, 5.47 cm, and 5.59 cm in the three directions, respectively. According to these results, the orbit products with the new SRP model have better accuracy than those with the reduced ECOM model, especially in the radial direction. The satellite-specific as well as constellation-mean accuracy comparison between the three experiments are given in Figure 6.
Accuracy comparison for experiments with different SRP models considering the GPS constellation as a whole
Satellite-specific accuracy comparison for experiments with different SRP models
The orbit products from experiments circ and beta are at about the same level of accuracy. This proves that the effect of the change of the beta angle is very small and the simpler SRP model of Equation (11) would suffice for the GPS constellation. So for experiment circ, more detailed statistics of product accuracy for individual satellites are presented in Figure 7.
Daily accuracy time series for individual satellites (color dots) and constellation mean (black lines) in experiment circ
Because both the new SRP model and the reduced ECOM model were designed for GNSS satellites in yaw-steering mode, the SRP models probably wouldn’t work for eclipsing satellites for some time, at least until the satellites returned to the nominal attitude after exiting the Earth’s shadow. This phenomenon may heavily degrade the orbit quality of the satellites in eclipse seasons. So, another accuracy assessment was conducted considering all non-eclipsing satellites as a whole this time. The comparison of orbit accuracies for the three experiments is plotted in Figure 5.
Accuracy comparison for experiments with different SRP models considering non-eclipsing satellites as a whole
In the case of non-eclipsing satellites as a whole, the mean orbit accuracy for experiment ecom was 4.17 cm, 6.39 cm, and 5.86 cm in radial, along-track, and cross-track directions, respectively. For the experiment circ, the three numbers were, respectively, 2.66 cm, 5.39 cm, and 5.50 cm. Then, the three values for the experiment beta were 2.65 cm, 5.43 cm, and 5.66 cm, respectively. From these specific numbers, we know that there exists an orbit accuracy decrease in the radial direction for eclipsing satellites. However, no apparent accuracy degradations are observed in the other two directions. The fact that Figure 5 has a very similar pattern as Figure 4 also indicates that the orbit accuracy of eclipsing satellites is not affected seriously by passing the Earth’s shadow.
The precision of the GPS clock products from these three experiments is also plotted in Figure 4 and Figure 5. In the case of all GPS satellites as a whole, the mean standard deviation (STD) values for the three experiments ecom, circ, and beta during the experiment time span were 3.70 cm, 3.08 cm, and 3.13 cm, respectively. When regarding the non-eclipsing satellites as a whole, the three numbers were 3.55 cm, 2.95 cm, and 2.98 cm. For each experiment, the improvement of clock accuracy after excluding eclipsing satellites from the product comparison was consistent with the change of orbit accuracy in the radial direction, as expected.
4.2 Efficiency
Processing efficiency is a critical indicator of real-time GNSS product generation because the timeliness of real-time products has a direct impact on the GNSS applications relying on these products. Before an efficiency assessment, the performance of the computer platform must be specified. The server computer on which the experiments were conducted was equipped with an Intel® Core™ i9-10850K CPU (3.6 GHz) and 62 GB RAM. Note that the software generating the real-time GNSS orbit and clock products is single-threaded.
To verify processing efficiency, the time series of epoch execution time was recorded. Figure 8 shows this time series to be about 264,960 epochs, and the average execution time per epoch was usually less than 1.0 s. There exist steep increases in execution time at day boundaries due to the decompression of compressed RINEX files, so they were removed from the series.
Statistics of epoch execution time
5 CONCLUSION
Real-time GNSS orbit and clock products are the prerequisite for the realization of real-time precise GNSS applications. For this purpose, a rigorous and optimal approach to generate these real-time products was deployed in our research with a prototype software developed. In addition, a new quasi-orbit-fixed SRP model and its generalized models were developed. To assess the performance of the new software and the new SRP model, three experiments with a tracking network of 60 stations and a time span of three months were conducted to simulate the generation of real-time GPS orbit and clock products.
After comparing with the final products from the IGS GFZ Analysis Center, the orbit accuracies for experiments using different SRP models were illustrated. With the new SRP model for circular satellite orbits, the mean RMS values for the real-time GPS orbit product were 2.82 cm, 5.45 cm, and 5.47 cm in radial, along-track, and cross-track directions, respectively. Compared to the orbit accuracy of 4.30 cm, 6.34 cm, and 5.57 cm in the three directions using the reduced ECOM model, the new SRP model makes a substantial contribution to the generation of real-time precise orbit products, especially in the radial and along-track directions.
The precision of the GPS clock product was also obviously improved from 3.70 cm to 3.08 cm. When considering the effect of the change of the β angle, the GPS orbit and clock products had no further apparent accuracy improvement. In terms of efficiency, the average epoch execution time of our software was less than 1.0 s, which is a good indicator for real-time GNSS applications. Thus, from the viewpoint of accuracy and efficiency, this research demonstrates the potential of this rigorous method to generate real-time precise GNSS orbit and clock products effectively, if the data collection and product dissemination are well operated.
HOW TO CITE THIS ARTICLE
Liu, P., & Chen, J. (2022). Real-time precise GPS orbit and clock estimation with a quasi-orbit-fixed solar radiation pressure model. NAVIGATION, 69(4). https://doi.org/10.33012/navi.549
CONFLICT OF INTEREST
The authors declare no potential conflict of interests.
ACKNOWLEDGMENTS
This work is mainly funded by Program of Shanghai Academic/Technology Research Leader; The Key R&D Program of Guangdong province (Grant No. 2018B030325001); National Natural Science Foundation of China (Grant No. 11673050); and the National Key R&D Program of China (Grant No. 2018YFB0504300). The authors would like to thank the IGS for providing the RINEX data. The authors also thank the anonymous reviewers for their comments to the manuscript. All of the figures in this paper were prepared using generic mapping tools (Wessel et al., 2013).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.