A review of real-time multi-GNSS precise orbit determination based on the filter method

Stable and reliable high-precision satellite orbit products are the prerequisites for the positioning services with high performance. In general, the positioning accuracy depends strongly on the quality of satellite orbit and clock products, especially for absolute positioning modes, such as Precise Point Positioning (PPP). With the development of real-time services, real-time Precise Orbit Determination (POD) is indispensable and mainly includes two methods: the ultra-rapid orbit prediction and the real-time filtering orbit determination. The real-time filtering method has a great potential to obtain more stable and reliable products than the ultra-rapid orbit prediction method and thus has attracted increasing attention in commercial companies and research institutes. However, several key issues should be resolved, including the refinement of satellite dynamic stochastic models, adaptive filtering for irregular satellite motions, rapid convergence, and real-time Ambiguity Resolution (AR). This paper reviews and summarizes the current research progress in real-time filtering POD with a focus on the aforementioned issues. In addition, the real-time filtering orbit determination software developed by our group is introduced, and some of the latest results are evaluated. The Three-Dimensional (3D) real-time orbit accuracy of GPS and Galileo satellites is better than 5 cm with AR. In terms of the convergence time and accuracy of kinematic PPP AR, the better performance of the filter orbit products is validated compared to the ultra-rapid orbit products.


Introduction
Providing real-time precise positioning services with Global Navigation Satellite Systems (GNSSs) is a significant development trend. Stable, reliable, and highprecision real-time satellite orbits and clocks are the prerequisites for real-time Precise Point Positioning (PPP) services. To achieve PPP Ambiguity Resolution (AR), the impact of satellite orbit error on GNSS measurements should be less than one-quarter of the Narrow Lane (NL) wavelength, e.g., 2.7 cm for Global Positioning System (GPS) satellites (Laurichesse et al., 2013). Currently, in support of multi-GNSS real-time applications, several Multi-GNSS Experiment (MGEX) and International GNSS Monitoring and Assessment System (iGMAS) analysis centers and research institutions provide multi-GNSS ultra-rapid orbit products. Typically, the Three-Dimensional Root Mean Squared Error (3D RMS) of the predicted orbit part with the ultra-rapid orbit product is approximately 5, 10, 5-10 and 10-20 cm for GPS, GLObal NAvigation Satellite System (GLO-NASS), Galileo navigation satellite system (Galileo), and BeiDou Navigation Satellite System (BDS) Inclined Geo-Synchronous Orbit (IGSO)/Medium Earth Orbit (MEO) satellites, respectively (Hadas et al. 2015;Kazmierski et al., 2018Kazmierski et al., , 2020. However, during eclipse seasons or satellite maneuvering periods, the accuracy of the predicted orbits is obviously degraded or even unavailable Dai et al., 2019;Duan et al., 2019;Laurichesse et al., 2013). For the emerging BDS satellites, some special issues, such as relatively imperfect orbit force models, different attitude modes, and frequent orbit maneuvers, will additionally affect the accuracy and reliability of orbit prediction. For instance, the BDS Geostationary Earth Orbit (GEO) satellite orbital maneuver occurs frequently, usually every 2 to 3 weeks (Song et al., 2009;Prange et al., 2017). During orbital maneuver periods, the ultra-rapid precise orbit products are not generally available. Therefore, it is more challenging to improve the accuracy and availability of real-time BDS satellite orbits compared to other GNSS satellites.
Satellite Precise Orbit Determination (POD) needs two kinds of information, i.e., the satellite geometric observations and satellite orbit dynamics. In essence, a satellite POD method is to establish an optimal trade-off between geometric observations and orbit dynamic information. Dynamic, kinematic, and reduced dynamic orbit determination methods have been proposed by assigning different weights to these two kinds of information. Geometric observations are mainly affected by observation equipment and the environment, and the information on satellite orbit dynamics is affected by the precision of the perturbation force model. Currently, real-time data streams are collected in global tracking networks with gradual data quality improvements (https:// igs. org/ rts/). Additionally, satellite orbit force models have been refined step by step. Considering the above situation, the filtering orbit determination method has the advantage of providing accurate and reliable real-time orbit products using a real-time data stream compared to the ultra-rapid orbit prediction method.
When the GPS was constructed, the filtering orbit determination method was employed to generate the broadcast ephemeris (Parkinson et al., 1996). Subsequently, this method was adopted by most commercial companies to provide real-time, high-precision, and high-reliability global services. In 2000, Real-Time GIPSY-OASIS (RTG) software, namely GIPSY-OASIS stands for GNSS-Inferred Positioning System and Orbit Analysis Simulation Software, was developed to support the Global Differential GPS (GDGPS) system by the Jet Propulsion Laboratory (JPL), which relies on a unit Upper triangular and Diagonal factorization (UD) decomposition filter and Square Root Information filter (SRIF). RTG can produce real-time satellite orbits and clock products, and the 3D accuracy of GPS real-time orbits is better than 10 cm (Bertiger et al., 1997(Bertiger et al., , 2012. By 2020, the User Range Error (URE) under Root Mean Square Error (RMS) of the GPS real-time orbit and clock products can reach 5 cm by RTG (Bertiger et al., 2020). RTG has been employed in several real-time services, including the Wide Area Augmentation System (WAAS) of the Federal Aviation Administration (FAA) Loh et al., 1995), the Japanese Multifunctional Satellite Augmentation System (MASAS), the GDGPS of the National Aeronautics and Space Administration (NASA) (Muellerschoen et al., 2001(Muellerschoen et al., , 2004, and Starfire, which was launched by Navcom (Dixon, 2006). In addition, Newcastle University developed the Auto-BAHN software, namely a software for near real-time GPS orbit and clock computation, by using an Extended Kalman Filter (EKF) based on the BAHN (Dow et al., 1993) software. Auto-BAHN calculates GPS real-time orbits with the observation data from 52 global tracking stations, and the 3D RMS is approximately 15 cm (Zhang et al., 2007). Additionally, the National Centre for Space Studies (Centre National d'Etudes Spatiales or CNES) has employed the Kalman filtering algorithm in real-time satellite orbit and clock determination, and the 3D accuracy of the GPS, GLONASS, Galileo, BDS MEO, and BDS IGSO satellite orbits is 5, 10, 18, 18 and 36 cm, respectively (Kazmierski et al., 2018). Tokyo Marine University and the Japan Aerospace Exploration Agency (JAXA) jointly developed a Multi-GNSS Advanced Demonstration Tool for Orbit and Clock Analysis (MADOCA) software, which is for post-and real-time orbit and clock determination. Real-time calculation with the EKF method can update the orbital state at a rate of 30 s and the satellite clock error at a rate of 1 s with the latency of 5 s (Takasu, 2013). MADOCA software is employed to provide a Centimeter-Level Augmentation Service (CLAS) by the Japanese Quasi-Zenith Satellite System (QZSS), and the 3D RMS of the real-time GPS orbit error is smaller than 4 cm (Allahvirdi-zadeh et al., 2021;Zhang et al., 2019). The CenterPoint RTX, a global precision real-time service system developed by Trimble, also performs UD decomposition filtering to determine multi-GNSS including GPS, GLONASS, Galileo, QZSS and BDS satellite orbit and clock products in real-time (Leandro et al., 2011;Glocker et al., 2012;Chen et al. 2011). GPS and GLONASS real-time precise orbits exhibit One-Dimensional (1D) accuracy of 3 and 6 cm, respectively, and the orbit accuracy during eclipse seasons is significantly higher than that of the International GNSS Service (IGS) Ultra-rapid (IGU) product. Laurichesse et al. (2013) substantiated this conclusion. Dai et al. (2016;2019a, b) and Duan et al. (2019) verified the advantages of the filtering method for real-time POD of GNSS satellites, especially for BDS satellites, because of their imperfect force models.
Currently, the filtering POD for GNSS satellites is adopted in the real-time positioning services mostly by commercial companies, and thus, the detailed description of the relevant processing methods and strategies is not fully available in the literature. In this study, we summarize the development status, key issues, and trends of multi-GNSS real-time filtering POD. First, we compare and analyze the advantages and disadvantages of the two main real-time POD methods, i.e., the ultra-rapid orbit prediction and the real-time filtering methods. Subsequently, the main issues, challenges, and solutions to real-time filtering orbit determination are examined. The fourth section outlines the research progress in the multi-GNSS real-time filtering orbit determination system, the accuracy of the real-time filtering orbit products, and the prospects for further work.

Real-time satellite orbit determination method
Real-time satellite POD approaches can generally be divided into two categories: the short orbit arc prediction based on batch post processing which is called the ultrarapid orbit prediction methods, and the real-time filtering orbit determination methods.

Ultra-rapid orbit prediction method
In the orbit prediction method, satellite orbit state parameters at the reference epoch of long orbit arcs, usually including satellite positions, velocities, and dynamic model parameters, are first calculated by the Least-Squares (LSQ) batch post processing method. Then, with these estimated state parameters, orbit integration is applied to generate predicted orbit products. The data processing flow is shown in Fig. 1. Typically, the dynamic model of GNSS satellites maintains a sufficient precision, and the real-time orbit generated via rapid extrapolation algorithm exhibits a high accuracy. Additionally, the IGS ultra-rapid products are generated with this approach, which are updated every six hours and have a 3D RMS of approximately 3 cm for observed GPS orbits and approximately 5 cm for predicted orbits (realtime orbits).
In the above method, the performance of orbit prediction depends on the accuracy of the orbit state and dynamic force model parameters. Hence, the relevant studies of this method have focused on precise orbit state parameter estimation, force model refinement, and orbit prediction strategies, such as those involving the observed and predicted arc lengths. Lou (2008) proposed  a short-arc orbit parameter rapid update method based on a normal equation combination of recent multiple arcs, which continuously updates the state parameters through a sliding short arc. The 3D RMS of predicted GPS orbits with an update rate of one hour is approximately 5 cm (Lou, 2008). Choi et al. (2013) analyzed the strategy of predicted GPS orbit arcs considering different fitting arc lengths for IGU orbits in 2010. The results demonstrated that employing fitting arc lengths of 40-45 h can yield the optimal orbit prediction accuracy. Additionally, Li et al. (2014) further analyzed the GPS orbit prediction accuracy using different fitting arc lengths and studied its influence on clock estimation and PPP AR. The results indicated that the use of the 42-h fitting arc length is optimal for clock estimation and PPP AR. Meanwhile, the real-time satellite orbit accuracy obtained with this method also depends on the accuracy of the dynamic model and the length of the predicted arc. Thus, the update interval of the IGU orbit was reduced from the original 12 h to 6 h, and several research institutions further reduced to 3 h or even 1 h to shorten the predicted arc (Liu, 2016;Deng et al., 2016;Li et al., 2018;Liu et al., 2019).

Real-time filtering orbit determination method
In the real-time filtering POD method, the orbit state parameters can be updated every epoch as soon as the observations are available and used as input. As a result, the estimated parameters can quickly respond to the reality of the satellite orbit state and are beneficial for handling the abnormal behaviors of satellite dynamics, such as orbit maneuvers. Generally, the Kalman filter or its evolution algorithm is used in most software. The Square Root Information Filter (SRIF) method proposed by the researchers from JPL (Bierman, 1977) was adopted in our process because of its high numerical accuracy and exceptional stability. In filtering recursive calculations, the SRIF algorithm uses the household orthogonal transform to update the measurement and time. The real-time POD process based on SRIF is shown in Fig. 2. The major issues in real-time filtering POD include the refinement of dynamic stochastic models for the different types of satellites, the handling of irregular satellite behavior patterns to ensure the stability and reliability of orbit determination results, the AR, and rapid convergence. These issues are discussed in more detail in Sect. "Key issues of real-time filtering POD".

Comparison between two methods
The two aforementioned methods are compared and analyzed in this section. In normal conditions, the satellite motion tracks are usually smooth, and the motion state can be expressed precisely with a dynamic model. Therefore, the ultra-rapid orbit prediction method can be used to generate high-precision orbit products. However, once the force acting on a satellite becomes abnormal, for example, when the satellite entering the eclipse seasons or in orbital or attitude maneuvers, the orbit state variation is hard to precisely represent using the original dynamic model, resulting in a significant accuracy degradation for the predicted orbits. In addition, the predicted orbit obtained with this method exhibits splicing during a certain arc update period. Jumps occur between adjacent update arcs, particularly for emerging navigation satellites with relatively large orbit modeling errors, and the jumps can reach up to several decimeters Dai Z et al., 2019). In contrast, the filtering orbit determination method is sensitive to the variations in the orbit state. When the dynamic model imposes an evident sudden deviation, the accuracy and reliability of real-time orbits can be improved by adaptively adjusting process noise and the weights between geometric observation information and dynamic information. However, compared to the ultra-rapid orbit prediction method, the real-time filtering POD technique has higher requirements in terms of the efficiency, reliability, and stability of data processing. A more specific comparison of these two methods is listed in Table 1. In addition, two issues for real-time processing should be noticed. Firstly, the collection of global datasets is more difficult for the real-time processing than that for the post-processing, because not all stations provide a real-time data stream. Also, the stability of data streams depend on the quality of the internet and data loss is unavoidable. Fortunately, there are more than 260 globally distributed stations of the MGEX network providing real-time data (Andrea et al., 2020). More attentions should be paid to the effect of data interruption on orbit state parameter reconvergence, which will be described in detail in Sect. "Rapid convergence and reconvergence". Additionally, in the real-time processing the data quality control is more complicated than in the post-processing which allows the data editing recursively. Teunissen (1990Teunissen ( , 1998Teunissen ( , 2018 proposed a DIA method for recursive detection, identification, and adaptation of model misspecifications by combining estimation with testing in dynamic systems that can be used for real-time GNSS processing. An adaptively robust Kalman filter that can resist the influence of both dynamic model errors and the measurement outliers has been well developed by Yang et al., (2001) and Yang & Gao (2006). Fu et al. (2019) and Zuo et al. (2021) developed the efficient quality control methods in the multi-GNSS real-time clock estimation to meet the requirement of high-rate update. Since the data quality control in the real-time POD is nearly the same as that for real-time clock estimation, it is not described in this paper.

Satellite dynamic stochastic model
The satellite dynamic stochastic model is a key issue in the filtering POD, which affects the quality of the filtering ig. 2 Filtering orbit determination flow orbit determination in long-term arcs. In principle, the process noise of the dynamic model should reflect the reality of satellite motion, but precisely determining this motion is challenging. If the process noise is too large, the dynamic information is limited to the estimations of the orbit parameters, and the orbit determination accuracy can be reduced. In contrast, if the process noise is excessively small, the orbital state variance matrix can become saturated, and the state parameters can be insensitive to the current observation information, leading to filtering divergence. At present, empirical values are mostly employed in stochastic models, but this approach does not analyze the dynamic noise characteristics, and the superiority of the filtering orbit determination method cannot be fully achieved (Dai et al., 2018;Duan et al., 2019). Figure 3 shows the average 3D RMS time series of the GPS satellite orbits determined with the filtering orbit determination method when different process noise levels are considered. The detailed noise selection strategies are listed in Table 2. As shown in Fig. 3, the dynamic model information cannot be fully utilized with weak dynamic constraints, and the orbit determination results fluctuate greatly and exhibit a low accuracy. When the dynamic constraints are very strong, even though the orbit accuracy can normally converge in the initial course of filtering, the orbit accuracy can deteriorate or even diverge in the later stage. Thus, determining a proper process noise level is essential to balance the weights between the dynamic model and the observation information. In further research, given the differences in the accuracy of the dynamic model according to the type of GNSS satellites, an optimized dynamic stochastic model should be refined using long-term historical observation data.

Handling of satellite maneuvers
Satellite maneuvers, including attitude and orbital maneuvers, can cause irregular satellite motions, such as entering the eclipse seasons and undergoing orbital adjustment. In this case, the satellite motion state will experience abnormal changes, leading to the accuracy degradation of the dynamic model and even unavailability. If adaptive adjustment is not introduced in the dynamic model, the filtering orbit determination method can have very bad performance or even collapsed.   The satellite maneuver time is important for the realtime orbit determination. The approximate satellite maneuvering epochs are published several months in advance in the Notice Advisory to Navstar Users (NANU) message (https:// www. navcen. uscg. gov). However, this early warning information is published only for GPS at present. Additionally, the maneuvering satellite is also marked as unhealthy in the broadcast ephemeris. However, the repositioned epochs from NANU are so rough that many effective observations may be missed by the users. Furthermore, the unhealthy marks in the broadcast ephemeris are sometimes misidentified or missing, so that the information received by general users is unreliable (Huang et al. 2008). Some efforts have been made to detect satellite maneuvers. In the post-processing POD, the maneuver epoch can be defined iteratively as the closest approach of the arc before and after the maneuver, and this definition has been adopted in the Center for Orbit Determination in Europe (CODE). In addition, the characteristics of orbit mutual differences (Ye et al. 2017) or station OMC (Observation Minus Computed) (Huang et al. 2008) calculated from the broadcast ephemeris can also be used for identifying orbit maneuvers. However, only the start time of orbit maneuvers can be detected by this method since the end time depends on the recovery of the broadcast ephemeris. Qiao & Chen (2018) proposed a carrier phase triple-differenced method to detect the start and end repositioning epochs for satellite orbit maneuvers that can also be used in real-time (Song et al., 2022). In the filtering POD, the DIA method based on the predicted residuals has been proven to be effective for satellite maneuver detection .
With the exact maneuver epochs, a simple method for maintaining a robust filtering POD is to reset all the parameters associated with the maneuvering satellite until the end of a maneuver. In this case, the satellite orbit cannot be determined during the maneuver period. Furthermore, a long reconvergence time that is greater than ten hours is needed after the end of maneuvers. A continuous filtering POD is supposed to improve the recovery of maneuvering satellite orbit, while the effect of the maneuver on the satellite dynamic model must be handled carefully. To compensate for the dynamic force model error, the main methods are divided into a functional model and a stochastic model. In the case of known maneuvering satellite telemetry data, the maneuver force can be modeled as an acceleration of a known magnitude, direction, start time, and end time, which is introduced as an additional force to the dynamic model of POD (Guo et al., 2013). However, the maneuver force model is usually not accurate enough because of the inevitable error in telemetry data. Moreover, the telemetry data are difficult to be obtained in many cases, so the introduction of additional empirical force parameters to eliminate the influence of orbital maneuvers is an alternative solution (Zhang et al., 2015). Some feasible methods have been proposed, including the constant empirical force method (Yoon et al., 2006(Yoon et al., , 2009, the piecewise linear empirical force method, and the pulse empirical force method (Lichten et al. 1987;Moon et al., 2012;Jaggi et al., 2012). Compared to the functional method of modeling the maneuver force, the stochastic process noise method is relatively easier to implement. Several studies have verified that this method can prevent filter divergence by reducing the weight of dynamic information during orbital maneuvers, and the achievable accuracy is similar to that of the kinematic POD method Duan et al., 2019). However, this method is not conducive to rapid reconvergence after orbital maneuvers because the dynamic force model contributes little to POD during maneuvers (Du, 2006). Hence, an improved approach, that adapts both the functional and stochastic models for the filtering orbit determination should be adopted. Yang et al. (2003) proposed an adaptive robust filtering method for satellite POD, which can effectively reduce the influence of observation and dynamic information anomalies on orbit state estimation. Xu et al. (2012) applied this method to GEO satellite orbit determination using simulated observation data, and the orbit accuracy reached approximately 1 m during maneuvers. Dai X et al. (2019) developed an adaptive SRIF POD strategy to handle maneuver events. Once a satellite maneuver is detected, the dynamic process noise of the maneuvering satellite is adaptively adjusted to realize continuous filtering orbit determination. Accordingly, Qing (2018) further adopted additional empirical force parameters to absorb the errors due to the maneuver force in the filtering orbit determination, resulting in continuous orbit determination during the maneuver period. With this approach, a meter-level accuracy of the BDS orbit can be obtained within six hours after the maneuver. In addition, RTGx software relies on a "decoupled" partition algorithm to handle unhealthy satellites in real-time data processing (Bertiger et al., 2020). To ensure high robustness of this method, the state parameters associated with the abnormal satellite are separated from those of additional satellites. This algorithm is generally equivalent to resetting the abnormal satellite-associated parameters and infinitely de-weighting the associated observations.

Rapid convergence and reconvergence
Typically, the cold start of the real-time filtering POD has a long convergence time, such as more than 12 h for GNSS satellites. Meanwhile, a long reconvergence time is required when a satellite is anomalous, such as orbit maneuvers occur. Due to the poor geometry of tracking satellites, the related geometric parameters, including tropospheric correction parameters and station coordinates, contribute little to rapid convergence in orbit determination. When tropospheric parameters or station coordinates are fixed to known values, there are insignificant effect on the fast convergence of the filtering POD (Qing, 2018). The convergence of the GNSS filtering orbit determination mainly depends on the accuracy of the initial satellite orbit states and dynamic parameters (primarily Solar Radiation Pressure (SRP) parameters), and the prior constraints. Typically, the initial satellite position and velocity can be directly computed using the broadcast ephemeris or the ultra-rapid orbit products, and the initial SRP parameters are estimated through dynamic fitting from these two types of orbit products. Fan et al. (2018) analyzed the influence of the initial orbit state with the broadcast ephemeris and ultra-rapid products on the convergence of the real-time BDS filtering POD. The results show that approximately 15 h are required to reach a stable accuracy level when the initial orbit state is from the broadcast ephemeris, while instant convergence can be achieved when the initial orbit state is tightly constrained to the ultra-rapid products. Qing et al. (2018) further examined the influence of the a priori constraints of the initial satellite position, velocity, and SRP parameters on the filtering orbit convergence. The results reveal that the constraint of the initial satellite position and velocity largely affects the convergence in the alongand cross-track directions, and the constraint of SRP parameters mainly affects the convergence in the radial direction. In the implementation of software, to realize rapid convergence and re-convergence, RTGx designed a "snapshot" function to store all the orbit state information in real-time filtering processing. Once the filter restarts, the precise initial orbit state can be obtained by reading the stored snapshot information, and instant convergence can be achieved (Bertiger et al., 2020).

Ambiguity resolution
Ambiguity resolution in real-time orbit determination is challenging. Because ionosphere-free observations are usually applied in the GNSS POD, the ambiguity parameter for ionosphere-free combinations does not exhibit integer properties. Therefore, to fix this integer ambiguity, the ionosphere-free ambiguity should be divided into Wide Lane (WL) and NL parts (Blewitt, 1989). In addition, the original ambiguity parameters include the initial phase delay and the signal hardware delay, which exist in both satellites and receivers. To utilize the integer value of the ambiguity, the unknown delay biases of the receiver and satellite can be eliminated using the double-difference between the satellite and the receiver. Once the double-difference ambiguities are resolved, the mapping relation between the double-difference and the undifferenced ambiguities can be employed to introduce ambiguity fixed conditions into the undifferenced observation equation (Ge et al., 2006;Teunissen, 1995). To improve the success rate and reliability of AR, Ge et al. (2005) proposed a method based on a whole network-independent ambiguity search. This AR method avoids calculating the Uncalibrated Phase Delay (UPD) or the Fractional Cycle Biases (FCB) of satellites and receivers, so it has been widely applied in GNSS data network post processing (Duan et al., 2019;. Most studies have focused on undifferenced AR, such as PPP AR. The key point of this approach is to accurately separate the integer parts from the fractional parts of the undifferenced ambiguity. The proposed methods can be divided into three categories: the UPD/ FCB method (Ge et al., 2008), the integer clock method (Laurichesse & Mercier, 2007), and the decoupled clock method (Collins et al., 2008). The implementation processes of these three methods differ slightly, but in practice, they are essentially equivalent (Geng et al., 2010;Teunissen & Khodabandeh, 2014). Because the undifferenced AR methods do not require independent baseline selection or whole-network independent ambiguity search, the processing efficiency is higher than that of the traditional double-difference ambiguity fixing method. Therefore, this approach has been gradually applied to GNSS network applications (Gong et al., 2018a;Dai Z et al., 2019;Kuang et al., 2021;Geng & Mao, 2021).
The following section demonstrates the recent progress of the UPD/FCB method in the real-time orbit determination for undifferenced AR (Dai et al., 2021). The RMS values of the WL and NL UPD residuals are small for GPS and Galileo, and more than 92% of the WL and NL UPD residuals remain within ± 0.1 cycles, which can meet the requirements of undifferenced AR in the real-time orbit determination. The average success rates of WL AR for GPS and Galileo are 91.84 and 95.64%, respectively. The success rate of WL AR for GPS is lower than that for Galileo due to the larger signal distortion error in the GPS pseudoranges (Gong et al., 2021). In terms of NL AR, the average success rates for GPS and Galileo are 92.45 and 92.94%, respectively. It should be noted that the NL ambiguity can only be fixed after the WL ambiguity is fixed. Thus, the number of WL and NL ambiguity candidates varies, and as a result, the success rate of GPS NL ambiguity fixing is larger than that of WL ambiguity fixing.

The performance of real-time GNSS orbit products
Based on the SRIF method, we developed real-time POD software that includes three main modules, i.e., the data receiving and managing module, the filter processing module, and the product dissemination module, as shown in Fig. 4. Typically, the input data of real-time POD processing include real-time observation data streams, predicted Earth Rotation Parameter (ERP) files provided by the International Earth Rotation and Reference Systems Service (IERS), and broadcast ephemeris or ultra-rapid products. The initial orbit state can be obtained using the broadcast ephemeris or ultra-rapid products, which can contribute to rapid filter convergence. The predicted ERPs from IERS are taken as the initial values in the filtering POD. The Universal Time 1 (UT1) that cannot be determined by GNSS is strongly constrained or fixed to the initial value, while the x-pole, y-pole, and LOD (Length of Day) are adjusted as random walk parameters with the initial constraints of 5 × 10 -4 arcsec, 10 -8 s/s, and process noises of 2 × 10 -5 arcsec and 5 × 10 -12 s/s in the correlation time of 300 s, respectively . The most complicated part of the software is the filter processing module, which comprises several important functions, such as real-time data preprocessing, orbit integration and force modeling, observation signal modeling, SRIF processing, filter quality control, and AR. The TurboEdit algorithm that combines Melbourne-Wubbena (MW) and geometry-free (Blewitt, 1990) observations is used in our real-time data preprocessing module to eliminate outliers and detect cycle slips. The clean observations are then fed into the filter module for WL UPD estimation. The filter process includes time updating and measurement updating, which require the satellite reference orbit and state transition matrix which are integrated from the initial orbit state based on the precise dynamic model. Because the computing time grows significantly as the number of satellites increases, orbit integration is carried out as a parallel process for each satellite to support real-time processing. After filter convergence, the undifferenced ambiguities are resolved, resulting in a significant improvement in the filter solutions. Finally, the real-time precise satellite orbit and clock products are converted into state-space representation correction expressions (SSR), then combined with UPD products and broadcast to support PPP and PPP Real-Time Kinematic (RTK) applications.
To test the real-time POD performance, observations collected at 110 globally distributed MGEX stations on DOY 059-072 of year 2021 were processed with our software to generate the orbit, clock, and UPD products for GPS, Galileo, and GLONASS satellites. The dynamic and observation models adopted in orbit determination are provided in Table 3. The orbit accuracy after filter convergence was evaluated by comparing with the CODE (Center for Orbit Determination in Europe) final products (Table 4). Figure 5 shows the average RMS values of the orbit differences in the along-track, cross-track, and radial directions for GPS, Galileo, and GLONASS satellites. The accuracy of the 1-2 h predicted orbit of the ultra-rapid products from GNSS Research Center of Wuhan University (WHU) (ftp:// 123. 57. 234.5/ pub/ gnss/ produ cts/ mgex/) is also shown in this figure for comparison. Compared to the float solution, the accuracy of the fixed solutions is significantly improved, from (5.7, 4.2, 2.4) to (2.9, 2.0, 2.1) cm and from (7.2, 3.9, 3.2) to (3.2, 2.3, 2.6) cm for the GPS and Galileo satellites in the alongtrack, cross-track and radial directions, respectively. The accuracy of the filter solution with AR is better than that of the ultra-rapid products for both GPS and Galileo, especially in the along-track direction. Because the AR algorithms implemented in our software are dedicated to the Code Division Multiple Access (CDMA) signals, we only generated the float orbit solutions for the GLONASS satellites. The orbit accuracy of the filter solution for the GLONASS satellites is 9.0, 8.6, and 4.7 cm in the alongtrack, cross-track and radial directions, respectively, which is much worse than that for the GPS and Galileo satellites. In this case, the filter orbit of GLONASS  (Petit & Luzum, 2010); ocean tide: FES2004 (Lyard et al., 2006) Solar radiation pressure GPS BLOCK IIA, IIR, III: Empirical CODE Orbit Model 2 (ECOM2) 7-parameter (Arnold et al., 2015) GPS BLOCK IIF, GLONASS: ECOM1 5-parameter (Springer et al., 1999) Galileo: Bow-wing (Rodriguez-Solano et al., 2012) + ECOM1 5-parameter, satellite properties from GSA (2017) Wu et al. (1993) Tropospheric delay Saastamoinen model+random walk process (Saastamoinen, 1972) Ambiguity Average orbit accuracy of filtering solutions and ultra-rapid products for GPS, Galileo, and GLONASS satellites compared with the final products of CODE from DOY 061 to 072, 2021. The 1-2 h predicted part of ultra-rapid products from WHU (wum) are used for comparison. The symbols of G, E and R refer to GPS, Galileo and GLONASS, respectively. The symbols of 'Filter' , 'Filter-AR' and 'Ultra' mean the orbit products of filter solutions without AR, the orbit products of filter solutions with AR and the ultra-rapid orbit products, respectively (   satellite is also worse than that of the ultra-rapid products in the cross-and along-track directions. It means that the observation and dynamic models of GLONASS satellites used in our filter processing need improving. In addition, the Satellite Laser Ranging (SLR) residuals of the filtering orbits are computed for external validation. Since SLR observations for GPS satellites are unavailable, only the SLR residual results for Galileo and GLONASS satellites are shown in Table 4. The outliers larger than 0.5 m were excluded in this analysis. With AR, the Standard Deviation (STD) and RMS of SLR residuals for Galileo satellites are improved from 4.4 and 4.4 cm to 3.9 and 4.2 cm, respectively. For the SLR residuals of GLONASS satellites that are considered to have smaller mean offset, the STD and RMS are much larger than those of the float solutions of Galileo satellites. This result may be attributed to the inaccurate force models of the GLONASS satellites or error correction models that need further investigation. The filtered orbit, clock, and UPD products were further validated with the kinematic PPP-AR at another 40 test stations. The satellite clock and UPD products were also estimated from the same 110 MGEX stations based on the WHU ultra-rapid orbit products, which were used in the PPP-AR experiment for comparison. Both the convergence time (Fig. 6) and accuracy (Table 4) of PPP-AR with the filtering products are better than those obtained with the ultra-rapid products. Since the NL ambiguity time-to-first-fix of the filtering orbit is shorter than that of the predicted orbit, there is a rapid dip convergence for filtering at approximately the 9th minute as shown in Fig. 6.
Currently, as more than 130 GNSS satellites are available, the processing efficiency is increasingly restricting the real-time POD of multi-GNSS. Considering the development of high-performance numerical computation technology, Gong et al., (2018a) proposed an efficient method of real-time network processing by combining the dense linear algebra algorithms and GNSS algorithm optimization. This method can improve the processing efficiency by approximately two orders compared with the traditional method, which is implemented in GNSS network data processing (Fu et al., 2018(Fu et al., , 2019. Based on this method, our experiment was conducted on Linux servers (main frequency of the Central Processing Unit (CPU) is 3.40 GHz; 24 cores; memory is 128 GB). In the experiment, the average number of parameters is 3636, and the average processing time is 6.36 s per epoch. According to the results, the time delay of filtering orbit is within 10 s per epoch. In real time processing, the loss of orbit prediction accuracy with a few minutes delay can usually be negligible. Thus, even considering the subsequent processing of the BDS, this method can fully satisfy the update requirements of satellite orbit parameters in the interval of 30 s.

Summary and outlook
To support real-time precise positioning services, real-time precise satellite orbit and clock products are essential. Compared to the ultra-rapid orbit prediction method, the real-time filtering orbit determination method has some advantages. It can quickly adjust the satellite orbit state information when the satellite motion is abnormal, which can obviously improve the stability and reliability of real-time orbit products. In this paper, the key issues and challenges that should be resolved in the filtering orbit determination method, including the refinement of dynamic stochastic models for different type of GNSS satellites, adaptive adjustment by filter models under satellite orbital maneuvers, rapid convergence, and real-time AR, are discussed. Subsequently, the self-developed software architecture, processing strategy, and the performance of real-time filtering orbit, clock and UPD products are presented. With undifferenced AR, the 3D RMSs of the GPS and Galileo satellite orbits can be improved by approximately 45% over the float solution, and the accuracy is better than 5 cm.
The satellite dynamic and measurement models used in the real-time POD are the same as those in the post-processing and summarized in Zhao et al. (2022). In addition to the improvement of these models, the following works will further improve the quality of GNSS real-time filtering POD.
1. For observations, the interruption of the real-time data stream is unavoidable, which significantly affects the stability, accuracy, and reliability of real-time orbit products. Therefore, the ultra-rapid orbit products can be introduced as virtual observations into filter processing. 2. Better satellite maneuver handling methods, particularly for the eclipse maneuvers and small abnormal dynamic changes that may contribute to orbit errors up to decimeter or meter level, need to be developed. 3. The real-time undifferenced AR algorithm for GLO-NASS satellites should be implemented in the filtering POD. Meanwhile, since the success rate of realtime AR for BDS satellites is much lower than those for the GPS and Galileo satellites, more efforts should be made to improve the performance of BDS satellites. 4. Satellite clock parameters are usually considered as white noise in the satellite orbit determination, which affects the accuracy of the satellite orbits due to a strong correlation between satellite clocks and orbit radial direction errors. Meanwhile, the current atomic clocks onboard GNSS satellites are extremely stable, and this correlation can be decoupled or weakened through precise satellite atomic clock modeling. The addition of a satellite clock model to the real-time filtering orbit determination and its impact on improving the accuracy and initialization time should be further studied. 5. The Earth rotation parameters used in real-time applications are commonly the predicted values for several hours or days ahead, which inevitably degrades the accuracy. In the real-time filtering orbit determination, as soon as new data are obtained, the ERPs, which mainly include the length of day (LOD) and x/y pole motions, can be estimated simultaneously with the orbit parameters in short time intervals, and the results can provide the information for real-time and high-frequency applications. However, the optimization of the stochastic noise and the strategy for ERP estimation in the real-time filtering POD still requires further investigations. 6. While the precision information of real-time products is quite important for real-time positioning applications, it is not fully provided by most realtime services. In the filtering process, the variance information of real-time satellite orbits and clocks can be provided for positioning as the a priori information. Additionally, integrity monitoring for realtime products should be implemented, which is more important for a reliable positioning service.