1. Introduction
The Earth Rotation Parameters (ERP), encompassing polar motion (PM) and its rates, along with UT1-UTC and its rate (length of day, LOD), describe the spatiotemporal variations in the Earth’s rotation vector [
1]. These parameters play a crucial role in defining a global reference frame and are determined through space geodetic techniques such as Very-Long-Baseline Interferometry (VLBI), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), Satellite Laser Ranging (SLR), and Global Navigation Satellite System (GNSS) [
2]. Due to the characteristics of the different techniques, these techniques are sensitive to different aspects of the Earth system [
3]. The GNSS technology is predominant in PM estimates because of the dense and uniform global distribution of the observing network, the continuous and high-precision observations, as well as the wide and stable satellite constellation [
4,
5]. The latest IGS Repro3 ERP products, contributing to the most recent ITRF2020, were generated without incorporating data from the Beidou Navigation Satellite System (BDS). [
6,
7]. The global BDS, namely BDS-3, has completed its construction which consists of hybrid constellations including Geostationary Earth Orbit (GEO), Inclined Geosynchronous Orbit (IGSO), and Medium Earth Orbit (MEO) satellites [
8]. It is worth noting that the BDS-3 MEO satellites were made by two entities including the China Academy of Space Technology (CAST) and the Shanghai Engineering Center for Microsatellites (SECM) [
9,
10]. The network of ground observing stations for BDS-3 are evenly distributed globally now. The five open service signals of BDS-3 satellites include B1I, B3I, B1C, B2a, and B2b [
11]. Notably, B1I and B3I signals are also transmitted by the regional BDS (BDS-2). The B1C, B2a, and B2b signals, unique to BDS-3, are designed to be compatible and interoperable with signals from other GNSSs [
12]. For the usage of an ionosphere-free combination, the official interface control document recommends B1C/B2a and B1I/B3I frequency combinations [
13,
14]. This research aims to assess the ERP derived from observations based on both BDS-3 B1I/B3I and B1C/B2a frequency combinations and seeks to identify the optimal dual-frequency combination for BDS-3 ERP estimates for the future ITRF definition with a consideration of BDS-3.
Many studies related to ERP derived by BDS-3 are based on B1I/B3I dual-frequency combination because most of BDS station receivers only support this frequency couple at the initial stage [
9,
15,
16]. Duan et al. investigated the influence of the solar radiation pressure (SRP) model on BDS orbit, geocenter motion, and ERP estimation as well as estimated the box-wing parameters based on B1I/B3I combination observations, concentrating on BDS satellites with a maximum PRNs of C37 due to receiver limitations from Multi-GNSS Experiment (MGEX) stations at the outset [
9]. Afterwards, Peng et al. enhanced ERP estimates from BDS satellites, incorporating box-wing models estimated based on B1I/B3I observations with PRNs extending up to C46 [
15]. Fang et al. showcased the advantages of hybrid constellation construction for ERP determination, utilizing IGSO and MEO satellites from BDS with B1I/B3I observations [
16]. Jia et al. further investigated the benefits of hybrid BDS constellation construction, considering GEO, IGSO, and MEO satellites, along with a hybrid SRP model [
17]. Moreover, He et al. compared the ERP estimated by BDS, GPS, and Galileo and analyzed the features of their systematic bias in ERP [
18]. In addition, Wang et al. studied the ERP-obtained discontinuous VLBI observations and BDS-3 based on B1I/B3I measurements [
19]. However, these studies for ERP derived from BDS-3 only utilized B1I/B3I observations without a consideration of using BDS-3 B1C/B2a observations.
Recently, many researchers studied the superiority of the BDS-3 B1C/B2a combination on the usage of Precise Point Positioning (PPP), clock offset, and orbit estimation. Yan et al. investigated the superiority of B1C/B2a observation quality for early BDS-3 [
20]. Li et al. performed a comparative analysis of BDS-3 orbit determination using B1I/B3I and B1C/B2a observations across nine MGEX stations. Their findings revealed that the accuracy of the BDS-3 orbit generated from B1C/B2a observations was improved by approximately 9% compared with the B1I/B3I-based orbits [
21]. Ye et al. evaluated BDS-3 orbits estimated using B1I/B3I and B1C/B2a observations, comparing them with orbits from other MGEX analysis centers. They also validated these orbits using Satellite Laser Ranging (SLR) data, demonstrating that the B1C/B2a-based BDS-3 orbits exhibit superior performance compared with the B1I/B3I-based orbits [
22]. He et al. further enhanced BDS-3 MEO orbits with an ambiguity resolution using B1C/B2a observations. Their evaluation results indicated that a higher percentage of ambiguity fixing, lower SLR residuals, and enhanced performance for PPP applications were achieved by using B1C/B2a combination [
23]. Geng et al. conducted a comprehensive comparison of data quality, precise orbit determination accuracy, PPP performance, and clock offset product accuracy based on B1C/B2a and B1I/B3I observations. They concluded that B1C/B2a-based results consistently exhibited better performance [
24]. However, these studies based on B1C/B2a dual-frequency combination observations did not cover the analyses of the ERP estimates, which need be further investigated. As of November 2021, with over 140 globally distributed MGEX stations supporting B1C/B2a signals [
24], there is a significant advantage for ERP estimates based on BDS-3 B1C/B2a observations.
In this research, we present our results as follows. Initially, the observation quality is comprehensively analyzed. Subsequently, the orbit determination results are assessed through orbit overlapping, a comparison with external products, and validation using SLR data. Following this, the ERP derived from BDS-3 B1C/B2a-based observations are thoroughly evaluated. Finally, we provide conclusions based on the analysis results.
2. Methodology
We selected 120 global evenly distributed MGEX stations that support the B1I/B3I and B1C/B2a observations of BDS-3, during the period from 1 October to 31 October 2023. The distribution of the chosen BDS-3 stations is illustrated in
Figure 1. Our ERP data processing schemes are outlined in
Table 1. Our data processing was performed on BSW54 [
25]. The dataset spans one month, covering days of the year (DOY) 274 to 304 in 2023, and we utilized double-difference observations with ionosphere-free (IF) combinations. To achieve high-precision orbit and ERP results, we implemented a two-step process: initially generating 1-day arc solutions, followed by stacking consecutive normal equations of 3 days to obtain the final 3-day arc solutions [
26]. ERP estimation was performed with a 24 h resolution, represented by epochs at noon and midnight, utilizing piece-wise linear functions for parameterization. For reference frame selection, we opted for IGS20, which is a realization of ITRF2020 for GNSS [
2]. The datum definition was achieved by using IGS20 core stations with minimal constraints including NNR and NNT to align with IGS20. To maintain consistency, we employed antenna models from IGS20.atx for both receivers and satellites, and the phase center variations (PCV) of BDS satellites were not taken into consideration. In this study, the SRP model employed ECOM2 with seven parameters [
27]. Sub-daily variations were corrected using the Desai–Sibois model, which is recommended by the 3rd IGS reprocessing campaign [
28]. Because of the correlation between orbital elements and UT1-UTC, the GNSS technology could not directly estimate the UT1-UTC parameter [
29]. Consequently, we tightly constrained the UT1-UTC of the first epoch to the a priori value from IERS EOP 20C04, which originates from VLBI technology. For the troposphere estimation, we use the Vienna Mapping Function (VMF3) grid product as the a priori model and zenith wet delay was parametrized as piecewise constants with a temporal resolution of 2 h. We use Chen–Herring model as an a priori model and 24 h resolution for troposphere gradients estimates.
3. Observation Quality
In this section, we conducted an analysis of observation quality, focusing on the Signal-to-Noise Ratio (SNR) and multipath errors for four representative stations. The selected example for this analysis is the C46 satellite.
Table 2 provides a summary of receiver types for the 120 selected stations. The receivers at these stations can be categorized into eight types, representing four different brands. To illustrate, we selected one representative station for each brand from the most commonly used receiver types.
The SNR is a critical indicator for characterizing the signal strength and the receiver’s ability to track signals. We analyzed the SNR for the four frequencies of BDS-3 from four representative stations with an elevation mask angle of 10°.
Figure 2 illustrates that SNR increases with the elevation angle, with the MET3 station exhibiting the best signal tracking performance. MET3 particularly excels in B2a, where the SNR can reach up to 60 dB-Hz. Moreover,
Table 3 provides the mean SNR for the four BDS-3 frequencies at the selected stations. In
Table 3, the mean SNR is approximately 3 dB-Hz stronger for a B2a signal than that for a B3I signal. We can conclude that the B2a signal demonstrates the highest SNR, followed by B3I, while B1C and B1I show comparable SNR values.
The multipath is one of obvious errors in GNSS observations. The multipath can be calculated as follows [
40]:
where
and
(
) denote the double frequencies.
is the multipath of frequency i and
stands for the pseudorange measurement of frequency i.
and
represent the wavelengths.
and
stand for the phase observation of frequency i and j.
denotes the phase ambiguity bias.
Figure 3 illustrates the multipath analysis of four different BDS-3 frequencies from four stations. The selected stations are from the MGEX network, without obstruction, and are situated far from water areas. The figure depicts that the multipath errors for most epochs are within 0.5 m, and these errors decrease as the elevation angle rises. The mean values and STD of multipath for different BDS-3 frequencies from the four stations are presented in
Table 4. The table shows that the mean multipath for each station is close to 0 m. In terms of STD, the B3I signal exhibits the least multipath, approximately 0.17 m, and the value of B2a signal is close to B3I. The B1I signal has the most significant multipath, which is about 0.26 m.
4. Ambiguity Resolution
Ambiguity resolution is a crucial step in achieving GNSS data processing with phase observations. Successfully fixing phase ambiguities to integers in double-difference data processing can be challenging due to observation noise [
23]. Double-difference ambiguities are generally categorized into wide-lane (WL) ambiguity and narrow-lane (NL) ambiguity. These two types of ambiguity need to be resolved sequentially. WL ambiguities are typically resolved using the Hatch–Melbourne–Wübbena (HMW) combination, followed by the calculation of NL ambiguities using ionosphere-free and WL ambiguities. The formulas for calculating WL and NL ambiguities are presented below [
23,
41]:
where
and
represent the two frequencies we selected for ionosphere-free combination.
and
are the phase observations while
and
are the code observations of the two frequencies.
and
represent the wavelength of WL and NL, which can be calculated by the wavelength
and
.
denote the WL ambiguity and
denote the NL ambiguity.
is the ionosphere-free combination ambiguity.
When
and
are fixed to integers, the double-difference ambiguity is resolved successfully.
Figure 4 presents the percentage of successfully resolved WL ambiguity and NL ambiguity using BDS-3 B1I/B3I and B1C/B2a combination. The figure indicates that the B1C/B2a achieves a higher percentage of ambiguity fixing compared with the B1I/B3I combination. In detail, the average percentage of fixed WL ambiguities is 89.1% for the B1I/B3I combination and increases to 93.8% for the B1C/B2a solutions. In the case of NL ambiguities, the mean percentage of fixed ambiguities is 68.8% for B1I/B3I and is notably higher at 79.1% for the B1C/B2a combination. Considering all the fixed ambiguities, including both WL and NL, the mean percentage is 78.9% for the B1I/B3I combination, while the B1C/B2a combination demonstrates a 7.6% improvement with a mean percentage of fixed ambiguities.
In addition, we conducted a comparison of the residual distributions for WL ambiguity and NL ambiguity between B1I/B3I and B1C/B2a combinations, as illustrated in
Figure 5. The black lines in
Figure 5 denote the probability density function (PDF) of the residuals, which follows the residual distribution. The plots indicate that the PDFs of WL and NL for both B1I/B3I and B1C/B2a combinations closely adhere to normal distribution. The residual histograms show symmetric distributions centered around zero cycles, particularly for NL ambiguities. The residuals of WL and NL ambiguities represent the disparities between the real values and the fixed values. The mean value, STD, and RMS of the ambiguity residuals are exhibited in
Table 5. It presents the average residuals for both WL ambiguity and NL ambiguity, indicating that the mean residuals are near to zero. The RMS values of the residuals for WL ambiguity and NL ambiguity are 0.069 and 0.057 cycles for B1I/B3I, whereas those for B1C/B2a are slightly smaller at 0.066 and 0.055 cycles, respectively. Consequently, the B1C/B2a combination demonstrates superior performance in ambiguity fixing compared with B1I/B3I.
7. Conclusions
The ERP constitute a crucial component in defining the global reference frame, and GNSS serves as a vital technology for obtaining ERP, including PM, its rates, and LOD. The latest IGS Repro3 ERP products provided the IGS contribution to the latest ITRF2020, but without consideration of the BDS observations. In this study, the ERP and the orbits generated from the BDS-3 B1I/B3I observations and B1C/B2a observations were calculated in the ITRF2020. The evaluation of these parameters was conducted, and the validation of 22 BDS-3 MEO orbits using SLR was also performed.
Firstly, we evaluated the signal quality from several indicators. SNR is an important index of characterizing the signal strength. The B2a signal has the largest SNR, then the second one is B3I and the SNR of B1C is comparable with B1I. The other indicator of observation quality is multipath. The ranking of multipath errors is B3I < B2a < B1C < B1I. Ambiguity resolution is one of crucial steps in GNSS data processing with phase observations, and the B1C/B2a combination has a higher percentage of ambiguity fixing than the B1I/B3I combination. The mean percentage of all fixed ambiguities including WL and NL is 78.9% for the B1I/B3I combination and the mean percentage for the B1C/B2a combination is higher, which is 86.5%.
Subsequently, the precision of the satellite orbits was compared by using orbit overlapping. The RMS for the B1C/B2a-based BDS-3 orbit overlapping was calculated for radial component, along component, cross component, and 1D component, yielding values of 1.58 cm, 2.09 cm, 1.86 cm, and 1.86 cm, respectively. It is noteworthy that the RMS values for the BDS-3 orbit overlapping, employing the B1I/B3I combination, are slightly higher compared with the B1C/B2a solution. This discrepancy may be thanks to the higher percentage of ambiguity fixing to integers in the B1C/B2a combination.
Satellite Laser Ranging, as one of key space technologies, can realize the orbit accuracy validation of GNSS satellites independently. The mean RMS of SLR residuals is 5.72 cm for 3-day arc orbits of 22 BDS-3 satellites from the B1C/B2a solution and 4.40 cm for BDS-3 satellite orbits without consideration of satellites C22, C43, and C44, which are both slightly better than the B1I/B3I combination and CODE MGEX products.
Finally, we assessed the accuracy of ERP estimates from the aspect of internal and external accuracy, particularly focusing on the formal errors in ERP estimation and the disparities between the estimated ERP values and the 20C04 products. The mean formal errors of the 3-day arc ERP estimates derived from B1C/B2a combination observations are approximately 10 μas for PM, 12 μas/d for PM rates, and 0.7 μs/d for LOD. Comparatively, the utilization of B1C/B2a-based solutions results in a reduction of mean formal errors by about 17% for pole coordinates, 22% for LOD, and 21% for pole coordinate rates when compared with B1I/B3I-based solutions.
Regarding the external accuracy, the mean differences of the 3-day arc ERP estimates from the BDS-3 B1C/B2a combination, with reference to IERS 20C04, are 48 μas for X PM, 19 μas for Y PM, 17 μas/d for X PM rate, 5 μas/d for Y PM rate, and 4.0 μs/d for LOD. The utilization of the BDS-3 B1C/B2a combination leads to a noteworthy improvement, with mean differences reduced by 56% for X PM, 54% for Y PM, 39% for the X PM rate, 64% for the Y PM rate, and 23% for LOD, in comparison with B1I/B3I. In the aspect of STD, the outcomes for the 3-day arc ERP estimates from the BDS-3 B1C/B2a combination are as follows: 63 μas for X PM and 47 μas for Y PM, 123 μas/d for the X PM rate, 101 μas/d for the Y PM rate, and 22.0 μs/d for LOD. The B1C/B2a solutions result in STDs of residuals that are approximately 3% for X PM, 2% for Y PM, 1% for X PM rate, 6% for Y PM rate, and 3% for LOD, better than the B1I/B3I solutions.
Therefore, the results demonstrate the excellent signal quality for the BDS-3 B1C/B2a and its superiority in the ERP estimation when comparing with the B1IB3I combination. Since the number and distribution of MGEX stations supporting B1C/B2a signals tracking are inferior to the B1I/B3I tracking stations, the performance of ERP derived from BDS-3 B1C/B2a signals maybe better with more global MGEX stations supporting BDS-3 B1C/B2a signals tracking in the future or adding Low Earth Orbit satellite observations.