Real-Time BDS-3 Clock Estimation with a Multi-Frequency Uncombined Model including New B1C/B2a Signals

: The global system of BDS (BeiDou Navigation Satellite System), i.e., BDS-3, is characterized with a multi-frequency signal broadcasting capability, which was demonstrated as beneﬁcial for GNSS (Global Navigation Satellite System) data processing. However, research on real-time BDS-3 clock estimation with multi-frequency signals is quite limited, especially for the new B1C and B2a signals. In this study, we developed models for BDS-3 multi-frequency real-time data processing, including the uncombined model for clock estimation and the GFIF (Geometry-Free Ionosphere-Free) combined model for IFCB (Inter-Frequency Clock Bias) determination. Based on the models, simulated real-time numerical experiments with about 80 global IGS (International GNSS Service) network stations are conducted for validation and analysis. The results indicate that: (1) the uncombined model with multi-frequency signals can achieve comparable accuracy with the traditional dual-frequency IF model in terms of clock estimation, and the double-differenced clock STDs (Standard Deviations) are generally less than 0.05 ns with post-processed clocks as a reference; (2) unlike the B1C and B1I/B3I signals, the satellite IFCBs generated from multi-frequency clock estimation show apparent temporal variations for B2a and B1I/B3I signals, further investigation with GFIF models conﬁrm the variations mainly result from the errors of receiver antenna corrections. Therefore, we addressed the feasibility of the uncombined model and the importance of accurate antenna information in the multi-frequency data processing.


Introduction
Real-time GNSS (Global Navigation Satellite System) applications with the PPP (Precise Point Positioning) technique call for corresponding high-quality satellite orbit and clock solutions. To serve the real-time end-users at worldwide scales, IGS (International GNSS Service) started providing real-time GNSS orbit and clock products with the RTS (Real-Time Service) project in 2013 [1]. Unlike the satellite orbit commonly generated through orbit integration [2], the satellite clock should be estimated with real-time observation streams due to the short-term variations. To improve the accuracy and efficiency of real-time clock solutions, research has been devoted to various aspects during the past decades, including estimation models [3,4], processing efficiency [5,6] and quality control [7].
While the traditional clock estimation usually incorporates dual-frequency signals with an IF (Ionosphere-Free) combination, the situation becomes different with the evolving of GNSS constellations. With the modernization of GPS (Global Positioning System) and GLONASS (GLObal NAvigation Satellite System) constellations and development of BDS (BeiDou Navigation Satellite System) and Galileo constellations, increasing satellites, including GPS block IIF and IIIA satellites, part of GLONASS-M and all GLONASS-K satellites, and all Galileo and BDS satellites, are broadcasting multi-frequency signals.
Multi-frequency signals were demonstrated as beneficial for data processing, such as cycle slip detection and repair [8,9] and PPP rapid integer ambiguity resolutions [10,11]. However, with more signals involved, the processing complexity increases accordingly since more linear combinations can be formed. Therefore, the UC (Un-Combined) model is preferred over the IF combined model in multi-frequency data processing [12,13].
The BDS, as a constellation transmitting multi-frequency signal, has gone through three phases, including the experimental system BDS-1, the regional system BDS-2, and the global system BDS-3. The BDS-3, completed in July 2020, consists of 24 MEO satellites, three IGSO satellites, and three GEO satellites, among which, the MEO and IGSO satellites are broadcasting OS (Open Service) signals on five distinct frequency bands [14,15]: legacy B1I (1561.098 MHz) and B3I (1268.520 MHz), and new B1C (1575.420 MHz), B2a (1176.450 MHz) and B2b (1207.140 MHz). In [11], the contribution of BDS-3 multi-frequency signals on PPP integer ambiguity resolutions was investigated. In [16], BDS-3 precise orbit determination with B1C and B2a dual-frequency signals was assessed. In [17], the real-time BDS-2/BDS-3 combined clock estimation with overlapping B1I and B3I signals was analyzed. However, the aforementioned research is limited to either IF combination or B1I and B3I signals.
In this study, we will develop the UC model for real-time BDS-3 clock estimation with multi-frequency signals and assess the consistency of BDS-3 multi-frequency signals with the GFIF (Geometry-Free Ionosphere-Free) combination [18]. Since receivers with the B2b signal track capability are quite limited for the current IGS network, the following models and experiments will be restricted to the other four signals, i.e., B1I/B3I/B1C/B2a. The remainder of this paper is organized as follows. After a detailed derivation of mathematical models for real-time BDS-3 clock estimation with UC multi-frequency signals, simulated real-time numerical experiments with IGS network observation data will be conducted and analyzed. Finally, some conclusions and suggestions will be given.

Materials and Methods
The raw GNSS code and phase observation equations can be expressed as: where Ci s r and Li s r denote the raw code and phase observations in meters, which were corrected for satellite and receiver PCVs (Phase Center Variations), dry part of the tropospheric delay and relativistic effects; ρ Li denotes the geometric distance between receiver and satellite antenna phase centers; cdt r and cdt s denote the receiver and satellite clock offsets; dt r,Ci and dt r,Li denote the code and phase hardware delays on the receiver end; d s Ci and d s Li denote the code and phase hardware delays on the satellite end; f i is the phase frequency; I s r,Lj is the 1st order ionospheric delay on Lj frequency; m s r is the map function of wet tropospheric delay; zwd r is the zenith wet tropospheric delay; λ i is the phase wavelength in meters; w s r denotes the phase wind-up effect in cycles; N s r,Li denotes the phase ambiguity in cycles; s r,Ci and s r,Li denote the code and phase residuals. The geometric distance between receiver and satellite antenna phase centers can be expressed as: where r s and r r denote the coordinates of satellite CoM (Center of Mass) and station MM (Monument Marker) in ECF (Earth-Centered-Fixed) frame; dr s PCO and dr r,PCO denote the satellite and receiver antenna PCO (Phase Center Offset) corrections; dr tides is the station tidal displacements, including solid Earth tides, ocean loading and pole tides; dr ARP is the station ARP (Antenna Reference Point) correction.
Among the coordinate correction terms in Equation (2), only the PCO correction term is frequency-dependent. Additionally, while the station tidal displacements and ARP correction can usually be accurately modeled, the PCO values from either ground calibration or network estimation may be subjected to various errors. As the PCO values are commonly defined in the SCF (Satellite-Centered-Fixed) frame and ENU (East-North-Up) frame for satellite and receiver respectively, Equation (2) can be further linearized as: where e s r is the line-of-sight unit vector from receiver to satellite; R ECF SCF is the rotation matrix from SCF frame to ECF frame; R ECF ENU is the rotation matrix from the ENU frame to the ECF frame.

Dual-Frequency IF Model
In terms of BDS-3 data processing, the IGS ACs (Analysis Centers) commonly incorporate B1I and B3I dual-frequency signals. To be consistent with the IGS ACs' convention, the dual-frequency IF combined clock with B1I/B3I signals, which assimilates code hardware delays on corresponding frequencies, is introduced with definition as: Inserting Equation (4) into Equation (1), the re-parameterized observation equations are obtained for B1I/B3I dual-frequency signals with IF combined model: where C26 s r and L26 s r denote the IF combinations of code and phase observations; ρ L26 denote the IF combination of geometric distance, λ LC is the IF combined phase wavelength in meters,N s r,L26 denotes the IF combined phase ambiguity in cycles with hardware delays assimilated. Once the geometric distance is corrected with known satellite and receiver coordinates, Equation (5) becomes the traditional dual-frequency IF combined model for clock estimation. Since the 1st order ionospheric delay is eliminated, unknown parameters are largely reduced for the IF model.

Multi-Frequency UC Model
Similarly, with the IGS dual-frequency IF combined clock datum, the observation equations can be obtained for B1I/B3I/B1C/B2a quad-frequency signals: where * s r, * is the modified parameter with hardware delays assimilated; IFCB r,Li and IFCB s Li are the receiver and satellite IFCBs (Inter-Frequency Clock Biases). Since the IFCB can in principle be expressed as linear combinations of the hardware delays on corresponding frequencies, the resulting IFCB estimates will be close to zero once the raw observations are corrected with external observation-specific biases. With geometric distances fixed with known satellite and receiver coordinates, Equation (6) becomes the UC model for multi-frequency clock and IFCB estimation.

Multi-Frequency GFIF Model
Based on Equation (1), the GFIF combination, which is commonly used for multifrequency signal consistency analysis, can be easily accessed for B1I/B3I/B1C and B1I/B3I/B2a triple-frequency signals: where C * s r and L * s r are the GFIF combinations of code and phase observations; ρ L * is the GFIF combination of geometric distance; B L * is the GFIF combination of phase ambiguity in meters; s r,C * and s r,L * denote the GFIF combinations of code and phase residuals. The GFIF combinations of geometric distances can be further written as: Combining Equations (3) and (8), the GFIF combination of geometric distances can be simplified as: As indicated by Equation (9), only frequency-dependent terms, i.e. PCO corrections, are remaining in GFIF combinations. Within IGS ACs' routine multi-GNSS data processing, antenna-related corrections are commonly extracted from the igs14.atx file. However, the receiver antenna calibrations are absent in igs14.atx for BDS-3 signals. As a result, the corresponding GPS L1/L2 values are used for adjacent BDS-3 frequencies. Unlike the situation at the receiver end, the satellite antenna information for BDS-3 B1I/B3I/B1C signals is released in the igs14.atx file. The missing values for B2a signals can be filled with the corresponding B2b values considering the frequency proximity. Obviously, the PCO errors at either satellite or receiver end will be projected into line-of-sight observations. The estimation of either satellite clock or IFCB will then be contaminated by the PCO errors.
Once the geometric distance for each frequency is known precisely, the GFIF combination of geometric distances can be safely reduced from Equation (7). Then, the GFIF combinations for B1I/B3I/B1C and B1I/B3I/B2a signals can be derived as: With IFCBs and ambiguity terms remaining, Equation (10) becomes the GFIF model for IFCB estimation. While the IFCB can be generated from either the UC model as Equation (6) or the GFIF model as Equation (10), the results will be identical only if the applied corrections, mainly PCO corrections, are consistent. With the GFIF combination commonly formed without correcting PCOs, the resulting IFCB will derivate from that with the UC model. The disagreement may further increase once incorrect PCOs are applied in the UC model.

Clock and IFCB Estimation
Finally, the overall flow chart for real-time BDS-3 clock and IFCB estimation is illustrated in Figure 1. The whole processing consists of two main procedures, i.e., data preprocessing and parameter estimation. The raw observations are firstly corrected for biases and errors. In the data preprocessing, the phase cycle slips are detected with the MW (Melbourne-Wübbena) and GF (Geometry-Free) combination. Then, with cleaned raw observations, the parameter estimation with various models developed above, including the IF model, the UC model and the GFIF model, are conducted to obtain the real-time clock and IFCB solutions. To reduce the side effect of remaining bad observations after data preprocessing, the post-fit residuals are edited with predefined empirical thresholds. The corresponding observation noises for linear combinations can be mapped from that of raw observations following the law of error propagation. Obviously, linear correlations exist between the satellite and receiver clock/IFCB parameters. To remove the intrinsic rank deficiencies, clock and IFCB datums, such as the constellation zero-mean constraint, should be introduced.

Results and Discussion
To validate the developed models for real-time BDS-3 clock and IFCB estimation, RINEX3 (Receiver INdependent EXchange format 3) observation data from about 80 globally distributed IGS network stations ( Figure 2) spanning 1 April to 8 April 2021 are used for the simulated real-time numerical experiments. All of the stations are able to track BDS-3 quad-frequency signals. As mentioned in the previous section, the satellite antenna information for BDS-3 B1I/B3I/B1C signals is available in the igs14.atx file, the corresponding values for the receiver antenna are still missing. However, within the igsR3.atx file for the IGS 3rd reprocessing campaign (http://acc.igs.org/repro3/repro3.html accessed on 29 January 2021), the antenna calibrations of about 40 stations, which are equipped with SEPT POLARX5 or SEPT POLARX5TR receivers, are made available for BDS-3 B1I/B3I/B1C/B2a signals. Taking account of the limited receivers with BDS-3 quad-frequency antenna information, the igs14.atx file is used when conducting clock estimation with the UC model, and the igsR3.atx file is used for IFCB determination with the GFIF model due to the smaller number of unknowns. Furthermore, to be consistent with the IGS ACs' clock definition, the code observations are corrected with the code bias product in OSB (Observation Specific Bias) format from CAS (Chinese Academy of Sciences). The detailed processing strategy and filter settings are summarized in Tables 1 and 2, respectively. The post-processed solutions from IGS ACs are used for clock accuracy assessment with the commonly used DD (Double-Difference) method to remove the clock datum differences. Given the DBD (Day Boundary Discontinuities) of post-processed clocks, the clock STD will be calculated on a daily basis. In all, 24 BDS-3 MEO satellites, well-distributed on three orbital planes, will be the focus of this study. The UC models with B1I/B3I dual-frequency signals, B1C/B2a dualfrequency signals, B1I/B3I/B1C triple-frequency signals, B1I/B3I/B2a triple-frequency signals and B1I/B3I/B1C/B2a quad-frequency signals are conducted for real-time BDS-3 clock and IFCB estimation. For comparison, the IF model with B1I/B3I dual-frequency signals is performed in addition to the five UC models.

Satellite Clock Comparison
The real-time DD clock time series with dual-frequency and multi-frequency signals are depicted in Figures 3 and 4, respectively. Apart from the convergence part on the first day, the DD clocks show an apparent DBD phenomenon. Two potential reasons may contribute to the clock DBD. First, the reference clocks, which are generated on a daily basis, will experience jumps at day boundaries and then impact the resulting DD clocks. Second, the real-time clocks are estimated with orbit fixed to GFZ post-processed product. The GFZ daily orbit is determined with independent 24-h observation data, and therefore will inevitably be affected by the DBD phenomena. The orbit DBD will affect the orbit interpolation around day boundaries and then have an impact on the corresponding real-time clocks.   The daily STD values averaged over the last 6 days are depicted in Figure 5 for all of the six models. The STD values are generally less than 0.05 ns, except for satellite C45 with the B1I/B3I/B1C and B1I/B3I/B1C/B2a UC models. The variation of STD values between each model is less than 0.01 ns with the exception of satellite C45. The largest STD value difference between the six models appears for satellite C45, for which the STD value with the B1C/B2a UC model is about 0.015 ns smaller than that with the B1I/B3I/B2a model. In all, the high accuracy and consistency of the real-time clock are clearly visible for all six The real-time clock solutions are similar for all of the six models with the exception of the B1C/B2a UC model. As for the other five models, while the DD clocks for most satellites fall in +/− 2.0 ns after convergence, satellite C45 shows a derivation of around 3.0 ns, which may have resulted from differences between bias products used in real-time and reference clock generation. With the convergence part neglected, the DD clocks exhibit rather high intra-day stability, which indicates a small STD value compared to the corresponding RMS value. Unlike diversities among satellites with the other five models, the DD clocks with the B1C/B2a UC model fall into two distinct groups, which vary between +1.2 ns and −0.8 ns. When further examining the satellite types it is found that the two groups correspond to satellites manufactured by the CAST (China Academy of Space Technology) and the SECM (Shanghai Engineering Center for Microsatellites), respectively. Since the GFZ post-processed clock is based on B1I/B3I signals, the intra-day stability of DD clocks with the B1C/B2a UC model suggests good consistency between the B1I/B3I and B1C/B2a signals. Moreover, we must emphasize that the satellite PCO corrections for the B2a signal are absent in the igs14.atx file and replaced with corresponding B2b values in the study. Being aware of this kind of model imperfection, the consistency between the B1I/B3I UC model and the B1C/B2a UC model is really inspiring.
The daily STD values averaged over the last 6 days are depicted in Figure 5 for all of the six models. The STD values are generally less than 0.05 ns, except for satellite C45 with the B1I/B3I/B1C and B1I/B3I/B1C/B2a UC models. The variation of STD values between each model is less than 0.01 ns with the exception of satellite C45. The largest STD value difference between the six models appears for satellite C45, for which the STD value with the B1C/B2a UC model is about 0.015 ns smaller than that with the B1I/B3I/B2a model. In all, the high accuracy and consistency of the real-time clock are clearly visible for all six models, which confirms the feasibility of the developed multi-frequency UC models. The constellation-mean STD values over the experimental period are summarized in Table 3. The statistics generally coincide with those inferred from

Satellite IFCB Characteristic
As by-products of clock estimation with multi-frequency signals, the satellite and receiver IFCBs are principally linear combinations of corresponding hardware delays. With the code hardware delays corrected at the satellite end, the resulting satellite IFCBs between the third signal B1C or B2a and the datum signals B1I and B3I are thought to be with a small deviation from zero. The satellite IFCBs between the B1C/B2a and B1I/B3I signals are presented in Figures 6 and 7 for the UC models with multi-frequency signals. The satellite IFCBs are generally less than 0.2 m, which confirms the inference above. While not explicitly shown here, the satellite IFCB differences are well below 1.0 cm between both models. Leaving out the convergence part, no obvious periodical signals are observed for the resulting satellite IFCBs between the B1C and B1I/B3I signals. However, apparent nearly diurnal variations with peak-valley values reaching 3.0 cm are found for the satellite IFCBs between the B2a and B1I/B3I signals. Recalling Equation (6) for multi-frequency data processing, we speculate that one of the main contributions is the inaccuracy of satellite and receiver PCO corrections applied for the B2a signal.  To further analyze satellite IFCB between the B2a and B1I/B3I signals, the IFCB estimation is conducted with the GFIF combination. Considering the handling of satellite and receiver PCO corrections, three different strategies, including both satellite and receiver PCO-corrected (as a benchmark), only receiver PCO-corrected and both satellite and receiver PCO-uncorrected, are adopted to obtain the satellite IFCB. As shown in Figure 8, the generated satellite IFCB with satellite and receiver PCO-corrected are very similar to those in Figure 7. This kind of similarity can be fully explained with the underlying equivalence of the UC and GFIF models. When it comes to the satellite IFCB with only receiver PCO-corrected in Figure 9, no clear differences can be observed compared to the benchmark ones. However, as shown in Figure 10, further ignoring the receiver PCO in the GFIF model leads to obviously more stable satellite IFCB estimations, which suggests the aforementioned temporal variations are mainly originated from incorrect PCO corrections at the receiver end. Moreover, the reduction of nearly diurnal signals in satellite IFCB with PCOs uncorrected suggests that the GFIF combinations of satellite and receiver PCO corrections should be close to, if not exactly equal to zero in Equation (9). In general, the temporal stability of satellite IFCB based on the GFIF model confirms the consistency between the B2a and B1I/B3I signals.
While the satellite IFCB shown in Figure 10 is already much more stable than the original ones in Figure 8, the underlying basis is that the GFIF combinations of PCO corrections can be fully ignored. However, the assumption is hardly valid in theory and the residual small variations for some satellites in Figure 10 confirm its failure. To further investigate the residual signals in the satellite IFCB, we here further use the PCO information in the igsR3.atx file to conduct IFCB estimations with the GFIF model. The resulting satellite IFCB is shown in Figure 11. Because of the relatively smaller number of stations used (about half of that for Figures 8-10), there exists a re-convergence on 6 April for one single satellite and the other satellites experience tiny jumps due to the zeromean constellation constraint. Except for the re-convergence, the temporal variations are generally disappeared for all satellites compared to Figure 10, which clearly demonstrates the consistency between the B2a and B1I/B3I signals. From the above analysis, we can conclude that the deficiency of antenna information will destroy model completeness and then impede GNSS applications with multi-frequency signals.

Conclusions
In this study, we thoroughly discussed the real-time BDS-3 clock estimation with multi-frequency uncombined signals. The multi-frequency UC model is first developed with B1I/B3I IF combined clock datum and then the GFIF model is derived for IFCB estimations. With data from about 80 global IGS network stations over one week, the simulated real-time numerical experiments are conducted for validation and analysis. The results indicate that the satellite clock estimation from multi-frequency UC models achieves high agreement with that from the traditional dual-frequency IF combined model. The STD of DD clocks reaches about 0.03 ns with respect to post-processed solutions. Different from the steady satellite IFCBs between the B1C and B1I/B3I signals, the satellite IFCBs between the B2a and B1I/B3I signals showed apparent periodical variations with multi-frequency UC models. Further investigation with the GFIF model using different PCO correction strategies imply the variations mainly result from imperfect receiver PCO corrections. Therefore, we address the urgency for receiver antenna information to achieve utmost accuracy in the multi-frequency data processing.

Data Availability Statement:
The observation data and precise products used in the research are available on the FTP of Wuhan University (ftp://igs.gnsswhu.cn/pub/ accessed on 29 January 2021).