Modeling and assessment of five-frequency BDS precise point positioning

Since its full operation in 2020, BeiDou Satellite Navigation System (BDS) has provided global services with highly precise Positioning, Navigation, and Timing (PNT) as well as unique short-message communication. More and more academics focus on multi-frequency Precise Point Positioning (PPP) models, but few on BDS five-frequency PPP models. Therefore, this study using the uncombined and Ionospheric-Free (IF) observations develops five BDS five-frequency PPP models and compares them with the traditional dual-frequency model, known as Dual-frequency IF (DF) model. Some biases such as Inter-Frequency Biases (IFB) and Differential Code Bias (DCB) are also addressed. With the data collected from 20 stations, the BDS dual- and five-frequency PPP models are comprehensively evaluated in terms of the static and simulated kinematic positioning performances. Besides, the study also analyzes some by-product estimated parameters in five-frequency PPP models such as Zenith Troposphere Delay (ZTD). The results of experiment show that five-frequency PPP models have different levels of improvement compared with the DF model. In the static mode, the one single Five-Frequency IF combination (FF5) model has the best positioning consequent, especially in the up direction, and in the simulated kinematic mode, the Three Dual-frequency IF combinations (FF3) model has the largest improvement in convergence time.


Introduction
In the early 1990s, Precise Point Positioning (PPP) was first proposed by R.R Anderle (Malys & Jensen, 1990;Zumberge et al., 1997). Compared with the traditional double differenced relative positioning PPP, using different combinations of signals can develop more observation models with more estimable parameters. At the beginning, the Ionospheric-Free (IF) model was mainly used with the corrections of satellite orbits and clocks. However, due to the errors of satellite orbits and clocks as well as other error corrections, such as satellite and receiver antenna phase center corrections, etc., no substantial progress was made in PPP until the emergence of Global Positioning System (GPS) satellite's precise orbit and clock products provided by International GNSS (Global Navigation Satellite System) Service (IGS) analysis center in the 1990s (Kouba & Héroux, 2001). Since then, PPP technology has been extensively studied, and the mathematical model has been continuously expanded. Due to its global coverage, all-weather, and high-precision, GNSS has become an important way for people to obtain position and time information. The early GPS and Global Navigation Satellite System (GLONASS) satellites transmitted two frequencies (Pan et al., 2019). However, with the completion of Galileo navigation satellite system (Galileo), BeiDou Navigation Satellite System (BDS) and Quasi-Zenith Satellite System (QZSS), GNSS has made a tremendous progress and the new generation of satellite navigation systems all use three or more frequencies to broadcast signals, even the new-generation GPS satellites, namely Block IIF satellites, transmit a third civilian signal L5 besides the legacy signals L1 and L2. The broadcast of more signal frequencies has gradually Open Access Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: wqx@cumt.edu.cn 1 School of Environment and Surveying and Mapping, China University of Mining and Technology, 1 Daxue Road, Xuzhou 221100, China Full list of author information is available at the end of the article made dual-frequency PPP model to triple-frequency or even quad-frequency PPP model.
As a global satellite navigation system developed by China, BDS has been highly valued by the government since the 1980s. The development process of BDS is divided into three phases: BeiDou-1 Navigation Satellite System (BDS-1), BeiDou-2 Navigation Satellite System (BDS-2), and BeiDou-3 Navigation Satellite System (BDS-3) (Yang et al., 2011). As the Table 1 shows, BDS-2 satellites are currently transmitting three signals on B1I, B2I and B3I frequencies while BDS-3 satellites can also broadcast several new signals, containing B1C, B2a, B2b and B2(B2a + B2b), besides the legacy B1I and B3I signals . Moreover, BDS-3 can provide unique short-message communication which BDS-1 and BDS-2 do not possess (Yang et al., 2019). Multi-frequency signals of BDS have gradually attracted the attention of scholars, Guo et al. (2016), Li et al. (2018), and Li, et al. (2020b)) developed the BDS-2 triple-frequency PPP models composed of B1I, B2I and B3I signals. Su and Jin (2019), Zhu et al. (2021), and Zhang et al. (2020) analyzed and evaluated the BDS-2/BDS-3 triple-frequency PPP models with B1I, B3I and one of B1C, B2b and B2a signals. Jin and Su (2020) and Li et al. (2020a) further validated BDS-3 quad-frequency PPP models with B1I, B3I, B1C and B2a signals. Furthermore, Su and Jin (2021) derived the Galileo five-frequency PPP model and verified its positioning performance. Although Inter-Frequency Clock Bias (IFCB), caused by the disagreement among satellite clock estimates of different frequency observations, was proved existing in GPS satellites, (Montenbruck et al., 2012), Cai et al. (2016), and Pan et al. (2017a) made a detailed study on the IFCB changes of different combinations with BDS-3 signals, and found that the IFCB of the BDS-3 satellites has no significant changes.
More and more academics focus on multi-frequency PPP models (Duong et al. 2020), but few studies are related to BDS five-frequency PPP model. Yuan et al. (2021) proposed a real-time cycle slip detection and repair method for BDS-3 five-frequency data, and Mi et al. (2021) researched the characteristics of receiver-related biases between BDS-3 and BDS-2 for five frequencies. With the addition of frequencies, the structure of parameters estimated by the model changes to a certain extent. For example, the code biases of the code observation on other frequencies are numerically different from the first two code observations, so it is necessary to introduce some additional parameters named Inter-Frequency Biases (IFB) (Deng et al., 2020) in five-frequency PPP models. For developing and analyzing five-frequency PPP performance with BDS-3 signals and coping with multi-frequency PPP technology in the future, this study establishes five BDS five-frequency PPP models with the uncombined and Ionospheric-Free (IF) observations, namely FF1, FF2, FF3, FF4, and FF5 models, and compares them with the traditional dualfrequency IF model, known as DF model. In particular, some biases including IFB and satellites Differential Code Bias (DCB) correction will be addressed.
This study contributes to the modeling and assessment of five-frequency float PPP with BDS. First, the mathematical and stochastic models of BDS five-frequency PPP are derived after the introduction of general observation models. The characteristics of five-frequency PPP models and DF model are then discussed. Then we will present the dataset and PPP processing strategies in experiment and assess five-frequency PPP models and DF model in terms of the static and simulated kinematic positioning performances. Furthermore, the study also analyzes some by-product estimated parameters in five-frequency PPP models such as ionospheric and tropospheric delays, receiver IFB, etc. Brief conclusions are finally given.

Five-frequency PPP models for BDS
In this section, the functional and stochastic models of one uncombined model (FF1) and three IF models (FF2, FF3, FF4) and one single IF combination model (FF5) will be deduced with five frequencies based on the general observation model. B1C, B2a, B1I, B3I, B2a/B2b are the first to fifth frequencies, respectively. At the same time, the characteristics of developed five-frequency PPP

General observation model
The code and carrier phase observations on a single frequency are as follows (Leick et al., 2015): where the superscript s and subscript r represent satellite and receiver, respectively; p s j and l s j are the code and phase 'Observed Minus Computed' (OMC) values, respectively; j denotes frequency (j = 1, 2, 3, 4, 5); µ s r is the unit vector of direction; x represents the vector of position correction to the a priori position; dt r and dt s indicate the receiver and satellite clock offsets; T r means the tropospheric delay of the signal path; γ j γ j = f 2 1 /f 2 j is the ionospheric factor; I s r,1 denotes the slant ionospheric delay at f 1 frequency; N s r,j represents the integer phase ambiguity; d r,j and d s j are the Uncalibrated Code Biases (UCDs) of receiver and satellite; b r,j and b s j mean the Uncalibrated Phase Delays (UPDs) of receiver and (1) p s r,j = µ s r · x + dt r − dt s + T r + γ j · I s r,1 + d r,j − d s j + ε s r,P j l s r,j = µ s r · x + dt r − dt s + T r − γ j · I s r,1 + N s r,j + b r,j − b s j + ξ s r,L j and d s IF m,n and d r,IF m,n represent the UCD of receiver and satellite after IF combination.
In this study, we use precision ephemeris products generated by the GeoForschungsZentrum (GFZ) German Research Centre for Geosciences based on the B1I/B3I IF combination to correct the satellite orbit and clock errors satellite; and ε s r,P j and ξ s r,L j represent the measurement noises of the code and phase observations. For convenience, the following notations are defined: where f m and f n represent different frequencies (m, n = 1, 2, 3, 4, 5; m � = n) , α mn and β mn represent frequency factors; d s,m,n DCB and d r,m,n DCB represent the Differential Code Bias (DCB) of satellite and receiver, respectively, IF m,n = α m,n · d s m + β m,n · d s n , d r,IF m,n = α m,n · d m,n + β m,n · d r,n (Kouba & Héroux, 2001). As such, the precise satellite clock correction dt s IF 3,4 contains a specific linear combination of B1I and B3I code biases: In addition to the precise satellite clock correction, code biases d s IF 1,2 − d s j should be corrected by satellite DCB products provided by the Multi-GNSS Experiment (MGEX) in code equations while in phase equations d s IF 1,2 will be mapped into ambiguities with UPD treated as float values (Guo et al., 2015). With precise satellite orbit, clock and DCB corrections, the code and carrier phase observations can be simplified as:

FF1: model using five-frequency UC observables
The FF1 model directly uses the corrected code and phase observation equations, which avoids noise amplification caused by a linear combination and retains the slain ionospheric delay (Zhao et al., 2019). The FF1 model is as follows: (5) p s r,j = µ s r · x + dt r + T r + γ j · I s r,1 + d r,j + ε s r,P j l s r,j = µ s r · x + dt r + T r − γ j · I s r,1 + (N s r,j + d s where d r,j j = 1, 2 is absorbed by the receiver clock offset and ionospheric delay, namely a and b, respectively: However d r,j ( j ≥ 3 ) cannot be completely absorbed into the ionospheric estimate. As such, an additional IFB parameter is required to compensate for the effect. Eventually, FF1 model can be written as: with

FF2: model with four dual-frequency IF combined observables
The FF2 model utilizes four dual-frequency IF combinations of B1C/B2a, B1C/B1I, B1C/B3I, and B1C/B2, which can be expressed as: where m = 1, n = 2, 3, 4, 5.The combined model of FF2 can be written as: (9) p s r,IF m,n = α m,n · p s r,m + β m,n · p s r,n l s r,IF m,n = α m,n · l s r,m + β m,n · l s r,n

Stochastic models
The precision of satellite measurements generally can be quantified as a function form related to the satellite elevation angle or signal-to-noise ratio (Wang et al., 2002;Satirapod & Luansang, 2008). Based on the stochastic model of the satellite elevation angle ( E ), observation errors, namely σ , will be quantified as: σ 2 = a 2 + b 2 / sin 2 E , where a = b are constants, generally set to be 0.003 m for carrier phase and 0.3 m for code observations. However, considering the relatively low orbit and clock accuracy of BDS satellites, the phase and code observation precisions are set to 0.006 m and 0.6 m respectively (Cai et al., 2015). With the assumption of uncorrelated observations on all five frequencies, the variance-covariance matrix of code in FF1 model will be written: where σ P represents the variance of code observations. Note that the vaiance-covariance matrices of phase observations are accordingly formed with σ P replaced by σ L . According to the error propagation law, the variancecovariance matrices of code in other models read: p s r,IF 1,2,3,4,5 = µ s r · x + dt r + T r + ε s r,IF 1,2,3,4,5 l s r,IF 1,2,3,4,5 = µ s r · x + dt r + T r + N s r,IF 1,2,3,4,5 + ξ s r,IF 1,2,3,4,5 with

Characteristics of dual-and five-frequency PPP models
In order to facilitate the analysis and assessment of the positioning performance for five-frequency PPP models, this study uses traditional DF PPP model as contradistinction. The characteristics of dual-and five-frequency PPP models are clearly exhibited in Table 2. Different models utilize different combinations of five frequencies, so the noise amplification factors of different models are unequal. The noise amplification factor of DF model with � , A 4 = diag � e 2 1 +e 2 2 +e 2 3 +e 2 4 +e 2 5 � B1I/B3I is 3.527 bigger than those of FF2 model except B1C/B1I which may affect positioning performance. The new combination models FF3 and FF4 proposed in study have nice noise amplification factors. FF5 model with all signals has a smaller noise amplification factor of 1.888 which is closed to the value of 1 in FF1 model.

Positioning performance evaluation
In this section, the data processing strategies used in the experiment will be first introduced. Then, the observation data selected from the MGEX stations are used to verify the static and simulated kinematic PPP performances in terms of convergence time and positioning accuracy. Some by-product estimated parameters, including ionospheric and tropospheric delays, and IFB will be presented. It is worth mentioning that although B1I and B3I signals are from BDS-2 satellites, the experimental data are only collected from BDS-3 satellites considering the equity for all models.

Data processing strategies
Since the BDS-3 system achieved its full operation capacity only recently, the number of stations that can receive five frequencies at the same time is limited worldwide. This study uses 20 stations provided by the MGEX of the IGS organization and the observation data with a 30 s sample interval for a ten-day period of Day of the Year (DOY) 283-292, 2021. Figure 1 shows the geographical distribution of the selected stations, which can receive five frequencies from BDS-3 satellites. In Fig. 1, different  Relativistic effect Corrected (Leick, A. et al., 2015) Phase windup Corrected (Wu J.T., 1992) Earth rotation Corrected (G. Petit and B. Luzum, 2010) Station coordinates Static: estimated as constants; simulated kinematic: estimated as white noise process

Receiver clock Estimated as white noises
Tropospheric delay Zenith Hydrostatic Delays (ZHD) are corrected with Saastamoinen model, and Zenith Wet Delays (ZWD) are estimated as random walk (1 × 10 -9 m 2 /s) (J. Saastamoinen, 1973) Ionospheric delay Estimated as random walk (1 × 10 4 m 2 /s) only in FF1 (Su & Jin, 2019) Receiver inter-frequency bias Absorbed by receiver clock in FF5 models or estimated as constants in other models Ambiguities Estimated as constant shapes and colors indicate different receivers and antenna type of stations, obviously the main antenna type is JAVRINGANT_G5T and JAVRINGANT_DM. Table 3 summarizes the processing strategy of BDS five-frequency PPP models, including dual-frequency PPP models. The precision ephemeris products provided by GFZ analysis center are used to correct satellite orbit and clock offsets (Deng et al., 2017), and the BDS code and phase observation precision are set to 0.6 m and 0.006 m, respectively. For BDS-3 satellites, the satelliteinduced code biases need not be corrected (Zhang et al., 2017. Due to the lack of receiver Phase Center Offset (PCO) and the satellite Phase Center Variation (PCV) values for BDS in igs14.atx file, their corrections are not applied. The position coordinates are modeled as constants and white noise in static and simulated kinematic PPP modes, respectively. The true values of station coordinates are referred to IGS solutions, and the difference  Fig. 2 Observation residuals for FF5 model at stations pots and wind (DOY 283, 2021) between the true value and the estimated coordinate is defined as the positioning error.

Static PPP performance
If the PPP model does not fully consider the biases, the observation residuals will not follow the zero-mean normal distribution (Nie et al., 2020;Pan et al., 2017b). Taking the FF5 model as an example, Fig. 2 shows the code and phase residuals of stations pots and wind for FF5 model on DOY 283, 2021, where the left and right figures represent respectively the observation residuals and the probability distribution of residuals. As shown in the figure, the code residuals are at a level of larger than meter l while the phase residuals reach centimeter level, where different colors represent different satellites. In the right figures the red line stands for a zero-mean normal distribution, which indicates the accuracy of FF5 PPP models.
This section will assess the static PPP performance in terms of convergence time and positioning accuracy. If the three-dimensional positioning errors, namely east (E), north (N) and up (U) direction are less than the 10 cm at the current epoch and the following twenty epochs (Li & Zhang, 2014;Lou et al., 2016) (DOY 283, 2021) coordinate components with an improvement of 3.8%, 2.6%, 31.2% and 11.5%, respectively. Different from the dual-and triple-frequency models, the positioning performance of uncombined model FF1 is poorer than other models, which may be because the increase in frequencies leads to more ambiguity parameters in FF1 model.  Fig. 1. It can be seen that the number of BDS-3 satellites observed worldwide can meet the positioning requirement of at least 4 satellites, even in the Antarctic area 9 BDS-3 satellites can be observed at least, but satellites that can be observed in South America and Africa is fewer. The visible satellite number and PDOP represent current visibility and infer the performance of BDS-3 as a global satellite system. Figures 10 and 11 present Zenith Troposphere Delays (ZTD), including corrected ZHD and estimated ZWD delay series, for different PPP models at stations arht and pots on DOY 283, 2021, where different colors represent different models. The RMS of troposphere errors for different PPP models are also shown in the figure. We use the troposphere products provided by MGEX with a sampling rate of 5 min as reference values, which have a typical formal error of 1.5-5 mm. It can be found that the  estimated ZTD for different PPP models at one station changes slightly in one day and the differences between the different models are small yet the performance of FF1 model is better, which illustrates the estimated ZTDs for the five-frequency PPP models are accurate. The IFB parameters of receiver are estimated as the constant in the BDS five-frequency PPP models. Take the FF1 model as the example, Fig. 12 exhibits the estimated IFB series at stations godn and gods on DOY 283/2021. It can be see that the IFB time series are stable over time and it's reasonable to model the IFB parameters as the constants within one day in the BDS five-frequency PPP models. Besides, the stations in same position have different IFB estimates owing to the type of receivers. Due to the lack of the receiver DCB correction corresponding to stations godn and gods in MGEX product, Fig. 13 shows the receiver IFB estimates difference with respect to MGEX reference IFBs at stations arht and bogt. The difference of IFB1 and reference is not close to zero value due to the combination factor of B1C/B1I valued -54.25. It can be seen that the estimated IFB values are very close to references which can further confirm our derivation of the estimable IFB parameters.

Simulated kinematic PPP performance
After modelling the receiver coordinates as a white noise process, the same data used in static PPP performance are reprocessed in simulated kinematic mode. Figures 14  and 15 display the simulated kinematic positioning errors of stations met3 and urum for DF, FF1, FF2, FF3, FF4 and FF5 PPP models on DOY 283, 2021. Unlike in static PPP, five-frequency PPP models make the convergence much faster than the dual-frequency PPP. Table 4 summarizes the median of RMS and convergence time for all models and Fig. 16 illustrates the boxplot of the convergence time where the RMS statistical method of kinematic PPP and the convergence criterion are the same as static PPP. As we can see, five-frequency PPP models improves both the convergence time and positioning accuracy, but the improvement of convergence time is more apparent. Especially with FF3 model it can reach 52 min, an improvement of 23.3%. As for three-dimensional RMS, FF5 model can achieve the best positioning accuracy in N and U directions valued 1.8 cm and 4.9 cm, improved by 18.2% and 9.3%, respectively; FF3 model can achieve 2.1 cm in E direction, improved greatly by 25.0%.

Conclusions
Using five signals of BDS-3, this contribution develops five BDS five-frequency PPP models, namely FF1, FF2, FF3, FF4, and FF5 models and evaluates the corresponding performances compared with the traditional dualfrequency IF model, namely DF model. Moreover, we describe some biases in detail emerging in five-frequency PPP models such as satellite DCB correction and IFB, etc. Besides, ionospheric and tropospheric delays, IFB parameters estimated are also evaluated.
The results show that the five-frequency PPP models have different positioning performance in static and simulated kinematic modes. In static scenarios, the positioning accuracy of the up component is improved greatly in the five-frequency PPP models which also contribute to convergence time. FF5 model, the best one, can achieve respectively 1.5, 1.0, 2.3 cm in three-dimensional positioning errors and 28.4 min in convergence time, improved by 3.8%, 2.6%, 31.2% and 11.5%, respectively. As to simulated kinematic scenarios, the five-frequency PPP models improve not only three-dimensional positioning errors but also convergence time, especially for FF3 model 52 min, improved by 23.3%.
Overall, compared with dual-frequency IF PPP model, BDS five-frequency PPP models can improve to a certain degree not only three-dimensional positioning errors but also convergence time However, one should be aware that the positioning performance of five-frequency PPP models is limited due to the limitations of the float ambiguity solution. Therefore, we will focus on the integer ambiguity solution with multi-frequency in the future.

Appendix
In this part, we will give the frequency factors in tripleand five-frequency IF combination.