GPS and BeiDou Differential Code Bias Estimation Using Fengyun-3C Satellite Onboard GNSS Observations

Abstract

Differential code biases (DCBs) are important parameters in GNSS (Global Navigation Satellite System) applications such as positioning as well as ionosphere remote sensing. In comparison to the conventional approach, which utilizes ground-based observations and parameterizes global ionosphere maps together with DCBs, a method is presented for GPS and BeiDou system (BDS) satellite DCB estimation using onboard observations from the Chinese Fengyun-3C (FY3C) satellite. One month worth of GPS and BDS data during March 2015 was exploited and the GPS C1C-C2W and BDS C2I-C7I DCBs were explored. To improve DCB estimation precision, the dual frequency carrier phase measurements leveled by code measurements were used to form basic observation equation. Code multipath errors of the FY3C onboard GPS/BDS observations were assessed and modeled as grid maps, and their impact on DCB estimation was analyzed. By correcting code multipath errors, the stability of DCB estimates was improved by 5.0%, 3.1%, 16.2% and 13.6% for GPS, and BDS geosynchronous orbit satellites (GEOs), inclined geosynchronous satellite orbit satellites (IGSOs) and medium Earth orbit satellites (MEOs), respectively. The monthly stability of FY3C-based DCBs was at the order of 0.1 ns for GPS satellites, 0.2 ns for BDS GEOs and 0.1 ns for BDS IGSOs and MEOs. By comparison to the ground-based DCB products issued by other institutions, FY3C-based DCBs showed stability degradation for BDS C02 and C05 satellites, while, for other satellites, the stability reached a similar or even superior level. The estimated FY3C receiver DCB stability was at the order of 0.2 ns for both GPS and BDS. In addition to the DCB estimates, the obtained vertical total electron content above the FY3C satellite orbit was also investigated and its realism was examined in physical and numerical aspects.

1. Introduction

DCBs (differential code biases) are usually identified as the signal timing biases existing between two different frequencies or channels within one GNSS (Global Navigation Satellite System) system, which originates from different hardware path delays and associates with both GNSS satellites and receivers. As a main error source of GNSS observations, the GNSS satellite and receiver DCBs have been long under investigation since the systematic differences between GPS P1 and P2 as well as C/A and P1 code observations were noticed [1,2,3,4]. GNSS satellite DCBs are essential for basic PNT (positioning, navigation and timing) services, especially for single-frequency users when higher accuracy is required [5,6]. Specifically, they are broadcast in navigation messages in different forms by different GNSS systems, e.g., TGD (timing group delay) in GPS and BDS system [7,8], etc. In addition, DCBs are also important factors in precise applications such as GNSS-based ionosphere sounding [9,10,11], precise point positioning (PPP) and PPP ambiguity resolution as well as GNSS satellite clock estimation [6,12,13]. Due to recent development in GNSS systems, issues on the satellite and receiver DCB calibration and monitoring become more and more important for explorations on multi frequency observations from single or multi GNSS systems.

Although GNSS satellite DCBs can be well determined within an anechoic chamber before launch, they suffer from variations due to causes such as temperature changes and satellite aging [14,15]. Thus, they must be re-calibrated during the satellite service period. Since 1998, the Ionosphere Working Group of the International GNSS Service (IGS) has been generating the GNSS satellite and receiver DCBs by using ground-based GPS data. Generally, in implementation, dual frequency GPS carrier-phase and code observations are employed and formed as geometry-free combinations for DCB estimation together with global ionosphere map (GIM) representation [9,11]. Since DCBs are rather stable, they are usually estimated as constants per each day for each satellite and each receiver. However, the ionosphere vertical total electronic content (VTEC) can be considerably variable in both spatial and temporal domains, thus they are often parameterized as spherical harmonics in solar-geomagnetic reference frame with piece-wise linear connections for representation of temporal variations.

This conventional method allows us to determine GPS satellite DCBs with a fine precision at the level of few tenths of ns when observations from a world-wide tracking network are processed [14,15]. However, for the new emerging GNSS systems, especially for regional systems such as BDS-2, they can hardly provide global coverage in their current constellation status, thus making GIM estimation difficult. To solve this problem, ionosphere delay has to be firstly calculated, e.g., by using GIM generated from GPS data, and then subtracted from the geometry-free combinations, thus leaving only satellite and receiver DCBs for adjustment [16,17]. However, GIM precision is typically 2–8 TECU (1 TECU equals to ${10}^{16}$ electrons/m² and approximates 0.16 m in GPS L1 frequency) [14], and can be even worse during strong solar or geomagnetic activities or in areas with very sparse data, which could bring in extra uncertainties into the observation equations.

The approach for onboard receiver DCB estimation proposed in [18] provides another idea for GNSS satellite DCBs estimation by using onboard measurements. In their implementation, the GPS satellite DCBs are directly taken from the IGS Analysis Center (AC) publishing values, whilst the slant TECs (STEC) are expressed as a product of VTEC and a geometric mapping function proposed in [19] on the basis of spherical symmetry assumption; by differencing between GPS satellites in one epoch, the VTEC parameters can be eliminated and finally only receiver DCB is left for estimation [18]. Practically, this approach can also be adjusted to estimate the GPS satellite DCBs if they are also treated as unknown parameters in constructing the observation equations and have a zero-mean constraint imposed on them to avoid rank deficiency. Following this idea, Lin et al. (2016) estimated GPS satellite DCBs using COSMIC and CHAMP satellite data. Their results indicate that the precision of GPS DCBs estimated from COSMIC satellite is about 0.791 ns, while that from CHAMP is 0.279 ns [20]. Wautele et al. (2017) also adopted this approach by incorporating the Jason-2 satellite onboard GPS data. In their study, the VTEC parameters are retained in normal matrix for estimation; however, only GPS P1 and P2 code observations are used to form geometry-free combinations. The estimated GPS DCBs are on average at 0.3 ns precision level when compared to the DCB solutions from other IGS ACs, while the monthly stability is about 0.1 ns [21].

The Chinese FY3C satellite carries a dual-system, multi-frequency receiver named GNOS (GNSS Occultation Sounder) onboard, which can generate GPS L1/L2 and BDS B1/B2 observations, thus giving a good opportunity for studying DCB determination. In this paper, we investigated DCB estimation by using the onboard data from FY3C satellite. More importantly, dual-system data, i.e., GPS and BDS, were processed and analyzed, and the impact of code multipath errors was considered to further improve the DCB precision. The remainder of the manuscript is organized as followed. Section 2 describes the DCB estimation method using onboard measurements and some practical considerations on DCB estimation with the FY3C measurements. In Section 3, the FY3C pseudorange multipath model as well as its impact on DCB estimation is analyzed. Section 4 demonstrates the DCB estimation results and analysis. Section 5 is the conclusion.

2. Data and Methodology

In this section, the DCB estimation method using onboard measurements is described. The data employed for this research are briefly introduced and the processing strategy on DCB estimation is given. The experiment design is also specified.

2.1. DCB Estimation Method

The onboard dual-frequency GNSS carrier-phase and code observations, which are used for DCB estimation, can be generally expressed as [22]:

$$P_{r_i}^s=R_r^s+c\left(d{t}_{r,k}-d{t}^s\right)+\frac{40.28}{f_i^2}·STEC+c\left(d_{r,i}+d_i^s\right)+{\mathsf{\Delta }}_r^s+M{P}_i$$

(1)

$$L_{r_i}^s=R_r^s+c\left(d{t}_{r,k}-d{t}^s\right)-\frac{40.28}{f_i^2}·STEC+B_{r_i}^s+{\mathsf{\Delta }}_r^s+d{\varphi }_{r_i}^s$$

(2)

where $P_{r_i}^s$ and $L_{r_i}^s$ are the pseudorange and phase observations in meters between GNSS satellite s and receiver r on frequency i (i = 1, 2), which represents for GPS L1/L2 frequency or BDS B1/B2 frequency, respectively; $R_r^s$ is geometry distance; k (k = 1, 2) is the GNSS system of satellite s, i.e., GPS and BDS in this case; $d{t}_{r,k}$ and $d{t}^s$ are the receiver clock offset for the kth system and satellite clock offset for satellite s, respectively; f is frequency, and c denotes light speed in vacuum; $STEC$ is the slant TEC; $d_{r,i}$ and $d_i^s$ are the instrument delays for receiver r and for GNSS satellite s on the ith frequency, respectively; $d{\varphi }_i^s$ is the phase wind-up errors in meter, while ${\mathsf{\Delta }}_r^s$ denotes all the other frequency-irrelevant observations errors such as relativistic errors and antenna phase center offsets and variations; $B_{r_i}^s$ is the phase ambiguity in meter on frequency i; and $M{P}_i$ is code multipath error.

By forming geometry-free combination for both pseudorange and phase observations, the frequency-irrelevant term can be canceled out, the instrument delays differences result in as the satellite and receiver DCB parameters while the ionosphere delay term STEC remains. Thus, the geometry-free combination, which is also called the ionospheric observation, can be denoted as:

$$P_{G_r}^s=P_{r_1}^s-P_{r_2}^s=\alpha ·STEC+DC{B}_{r,k}+DC{B}^s+M{P}_{1,2}$$

(3)

$$L_{G_r}^s=L_{r_1}^s-L_{r_1}^s=-\alpha ·STEC+B_{r_{1,2}}^s+d{\varphi }_{r_{1,2}}^s$$

(4)

where $\alpha =40.28·\left(\frac{1}{f_1^2}-\frac{1}{f_2^2}\right)$, $DC{B}_{r,k}$ is the DCB of receiver r for the kth system while $DC{B}^s$ is the DCB for satellite s; and $M{P}_{1,2}$, $B_{r_{1,2}}^s$ and $d{\varphi }_{r_{1,2}}^s$ are the multipath, ambiguity, and phase wind-up difference between pseudoranges or phase observations on different frequency, respectively. The phase wind-up term $d{\varphi }_{r_{1,2}}^s$ can be accurately calculated from model proposed by [23] when GNSS satellite orbit and receiver position are provided. The code multipath errors can be calculated from constructed grid models, which will be described in Section 3.1. If only pseudorange observables are used for DCB calibration, as seen in Equation (3), only satellite and receiver DCBs and STEC are left to be considered while the ambiguities can be ignored, which is the case in [21]. However, if the high-precision carrier-phase observables are preferred, the ambiguities should be carefully handled and subtracted from Equation (4).

In this study, the carrier-phase observables are used for DCB estimation with the code multipath errors considered. Thus, the carrier-phase ambiguities have to be removed. By adding Equations (3) and (4) together, the ionospheric delay part can be eliminated. Since $B_{r_{1,2}}^s$, $DC{B}_{r,k}$ and $DC{B}^s$ can all be considered as constants for a continuous arc (i.e., a certain time span during which a specific GNSS satellite is tracked) when no cycle slips occur, they can be calculated as one term by averaging $P_{G_r}^s+L_{G_r}^s$ combinations for the whole arc. Thus, we have:

$$B_{r_{1,2}}^s=⟨P_{G_r}^s+L_{G_r}^s⟩+⟨M{P}_{1,2}⟩-DC{B}_{r,k}-DC{B}^s$$

(5)

where $⟨·⟩$ means the average value for an arc. It should be noted that, during the averaging process, the phase wind-up term $d{\varphi }_{r_{1,2}}^s$ should be first calculated and subtracted from $L_{G_r}^s$ observables. Taking Equation (5) back to Equation (4) when $⟨P_{G_r}^s+L_{G_r}^s⟩$ is obtained, the ambiguity term can be eliminated from the carrier-phase ionospheric observable:

$${\tilde{L}}_{G_r}^s=⟨P_{G_r}^s+L_{G_r}^s⟩-L_{G_r}^s=\alpha ·STEC+DC{B}_{r,k}+DC{B}^s-⟨M{P}_{1,2}⟩$$

(6)

where ${\tilde{L}}_{G_r}^s$ is the carrier-phase ionospheric observable leveled to the pseudorange ionospheric observable. The above procedure is usually called the carrier-to-code leveling process. It can be seen that the pseudorange multipath errors contribute to the ionospheric observations through $⟨M{P}_{1,2}⟩$. The above equation shows that the DCB parameters are actually resulted from the low-frequency part of the ionospheric observations, while the high-frequency part represents for STEC variations. As STEC corresponds to the integral of the electronic density along the receiver-to-satellite path, the STEC can be expressed as a function of the VTEC and their relationship is denoted as:

$$STEC=mf·VTEC$$

(7)

where $mf$ is mapping function. Since FY3C satellite orbits above the Earth surface at about 836 km, which is actually at the altitude of plasmasphere, the mapping functions used for ground stations such as in [9,24] are not applicable. As indicated in some studies [18,20,25], the geometric mapping function (F&K) proposed in [19] can be used if spherical symmetry is assumed. The F&K geometric mapping function can be expressed as:

$$mf=\frac{1+\left(R_E+h_{shell}\right)/\left(R_E+h_{leo}\right)}{\cos z+\sqrt{{\left(R_E+h_{shell}\right)}^2/{\left(R_E+h_{leo}\right)}^2-{\left(\sin z\right)}^2}}$$

(8)

where $R_E$ is the Earth radius, z is the GNSS satellite elevation angle observed from low Earth orbit (LEO) satellite, and $h_{shell}$ and $h_{leo}$ are the ionosphere effective height (IEH) and LEO satellite orbit height, respectively. The only unknown parameter in Equation (8) is $h_{shell}$. Different mapping functions and IEH calculation methods are evaluated for space-borne data processing in [25], and it is found that the F&K mapping functions along with IEH from centroid method is more suitable for LEO-based DCB estimation and TEC conversion than other models mentioned in their study, which shows a relative error generally below 10% at 800 km altitude during different solar activities. In our study, the F&K mapping function with centroid-derived IEH is adopted. The IEH is computed using below equation proposed in [25]:

$$h_{shell}=\left(0.0027F_{107}+1.79\right)·h_{leo}-5.52F_{107}+1350$$

(9)

where $F_{107}$ is solar radio flux at 10.7 cm, which is often used for representation of solar activity.

Based on above description, the leveled carrier-to-code ionospheric observables can be expressed in a closed form where only VTEC, GNSS satellite DCBs and receiver DCBs remain unknowns. Normal matrix can be constructed by linearizing and stacking Equation (6), and then these unknowns can be estimated using a least-squares adjustment. Since the GNSS satellite DCBs and receiver DCBs are closely coupled resulting in normal matrix rank deficiency, an external constraint condition is needed for de-correlation. In convention, a zero-mean reference is imposed on all the GPS and BDS satellites involved in the estimation. The zero-mean reference is implemented by constraining the sum of GPS satellite DCBs and BDS satellite DCBs to zero with a very large weight, respectively. As VTEC parameters are estimated for each epoch as process parameters, a pre-removing and recovery procedure is implemented to improve efficiency [26].

It is also noted that, since, in our study, only dual-frequency GPS/BDS data are used for DCB estimation, the above method makes use of the classical definition of the DCBs and investigates the DCBs between two frequencies; however, if multi-frequency data are provided from the space-borne receiver, this method can also be adjusted to estimate the code biases on the third and higher frequencies as in that case they become estimable [27].

2.2. Data Description and Processing Strategy

Using the method introduced in Section 2.1, the GPS and BDS DCBs can be estimated with the LEO onboard measurements. In our study, the FY3C onboard GPS/BDS pseudorange and carrier-phase observations during March 2015 are exploited. These data are issued by National Satellite Meteorological Center, China. For the GPS, the dual-frequency pseudorange C1C/C2W and carrier-phase L1C/L2W are provided, while, for BDS, C2I/C7I and L2I/L7I are provided. This enables estimation of the GPS C1C-C2W DCBs and BDS C2I-C7I DCBs. Both GPS and BDS data are provided in 1 Hz sampling rate, which can be very useful in revealing high-frequency VTEC variations.

For DCB estimation, the FY3C data are processed for each day in 1 s interval without any loss of information. To obtain higher precision DCB estimates, the leveled carrier-to-code observations are selected rather than the pseudorange observables used in [21]. To reduce mapping function error, an elevation mask is set as 40°, which is actually larger than that used in previous studies. This is because that FY3C can provide more data with dual-GNSS system availability than the other LEO satellites, thus the drawback of dramatic data volume decrement due to a high elevation cutoff can be avoided. The weighting strategy is empirically considered as proportional to satellite elevations. In our study, it is implemented as $P=\frac{{\sigma }_0^2}{{\sigma }_i^2}, {\sigma }_i=\frac{{\sigma }_o}{\sin e}$, where P is the observation weight and e is the satellite elevation angle, while ${\sigma }_o$ and ${\sigma }_i$ are the a priori variance and ionospheric observation variance, respectively. For leveled carrier-to-code ionospheric observables, ${\sigma }_o$ is empirically set as 0.1 m since the posterior residuals are generally at this level. The GNSS satellite and receiver DCBs are estimated daily as constants while the VTEC parameters are estimated per each epoch [21]. To acquire better VTEC continuity, a very loose constraint is applied between consecutive VTEC parameters, i.e., $1m/\sqrt{h}$.

For analysis of FY3C multipath impact, DCB estimation with and without the multipath correction model applied are also investigated. Thus, there are in total two different test cases designed in our study, of which the first is denoted as L4-raw using raw leveled phases and the second L4-cor using corrected leveled phases.

3. FY3C Pseudorange Multipath and Its Impact on DCB Estimation

In this section, the multipath modeling for FY3C GPS and BDS code observations is briefly introduced. Then, its impact on the DCB estimation is investigated through a daily DCB estimation case.

3.1. Multipath Modeling and Validation

As mentioned in Section 2, the code multipath errors can contribute to the leveled carrier-to-code ionospheric observables, thus affecting the DCB estimation. Multipath behaviors can be quite different and more complicated for onboard data compared to ground stations. The causes are various, and are mainly related to cross-talk between different antennas (e.g., navigation antenna and occultation antennas), different configurations of the solar array panel and different antenna designs. Usually, the multipath pattern can be quite stable for a specific onboard antenna, as indicated by the multipath maps of CHAMP, GRACE, SAC-C and TerraSAR-X satellites shown in [18]. In such case, it is possible to model the multipath errors as grid maps in the antenna-fixed reference frame (ARF).

For BDS-2 signals, different multipath patterns are found for different constellations, i.e., BDS GEOs, IGSOs and MEOs. As indicated by ground-based BDS data, strong low-frequency elevation-dependent variations are observed for BDS IGSOs and MEOs, and can be well modeled as piece-wise linear functions or polynomials with respect to elevation angles [28,29], while the GEOs multipath generally exhibits similar signatures to IGSOs but shows poorer fitness [29]. These are seen as the multipath errors induced by the BDS satellites. Thus, for FY3C BDS code multipath modeling, different multipath models should be estimated for different BDS constellations.

For FY3C GPS and BDS code multipath model construction, the multipath errors are first calculated using the multipath combinations (e.g., [30]). Please note that the pseudorange observation noise also contributes to the calculated multipath errors, as these two terms are un-isolatable from the multipath combinations. Then, the multipath maps are estimated as piece-wise linear functions in a 2° × 2° resolution in the ARF for each frequency of GPS, BDS GEOs, BDS IGSOs and BDS MEOs, respectively. Through this parameterization, the estimated multipath model actually contains combined contributions from both the near-field multipath and the BDS-satellite induced multipath errors. The detailed patterns of the model can be found in our former study [31]. The FY3C multipath patterns are quite complicated, which could be attributed to the combined effects from BDS-2 satellite-induced errors, FY3C solar panel rotations as well as cross-talk effects from the occultation antennas, etc.

Figure 1 shows the FY3C GPS/BDS multipath RMS (root mean squares) values against the satellite elevation angles. The red dots in left panels represent the RMS values of raw GPS L1 or BDS B1 frequency code multipath errors (MP1), while the blue dots represent the MP1 RMS values with multipath grid map applied; the red and blue dots in right panels are the original GPS L2 or BDS B2 frequency multipath (MP2) RMS and the grid model corrected MP2 RMS values, respectively. As only observations with elevations larger than 40° are used for carrier-to-code leveling, the statistics are based on observations with elevation between 40° and 90°. In Figure 1, it can be found that: (1) The code multipath error generally reduces while elevation grows; however, abnormal large multipath errors are observed for elevations close to 90° in all the panels. This is because, firstly, there are fewer data for observations near 90° elevations; and, secondly there are frequent sudden pseudorange jumps when elevation reaches near 90°. These jumps can be as large as 2–5 m, and the reason is still under investigation. (2) Improvement can be observed for all cases. For GPS, the improvement is more prominent in lower elevations, while, for BDS, in higher elevations. Especially for the BDS MEOs, the MP1 RMS can be reduced from 0.6–0.9 m to about 0.3 m when elevations are larger than 70°, while MP2 RMS reduced from 0.4–0.6 m to under 0.3 m. This can be attributed to the BDS satellite-induced biases as shown in [31]. (3) For GPS, BDS GEOs and BDS IGSOs, their average improvement is generally about 5–10%, while, for BDS MEOs, the MP1 and MP2 improvement is 56.4% and 34.2%, respectively.

3.2. Multipath Impacts on DCB Estimation

To demonstrate the impact of multipath errors on DCB estimation, a daily solution on DOY (day of year) 062/2015 is taken as an example. It can be obviously seen from Equations (5) and (6) that the multipath error affects the leveled phase observation through its arc-averaging value, which indicates that it can bias the leveled phase measurements for the entire arc. If no code multipath correction is applied, the calculated $P_G+L_G$ series would show larger scatters, resulting in larger standard deviation (STD) of its arc-average $⟨P_G+L_G⟩$.

As the multipath errors are most prominent for BDS MEO satellites, the C14 $P_G+L_G$ series for a continuous arc is shown in Figure 2. The red lines show the original $P_G+L_G$ series while the blue lines are with multipath correction applied; the grey line is the elevation. For illustration, their average values, which are mainly the DCBs and ambiguities, are subtracted. The $P_G+L_G$ combination is used for leveling the phase observations; thus, the smaller the variations are, the better precision the leveled carrier-to-code observables are. Significant improvement can be found for corrected series, especially for high elevations, as they suffer from larger biases. After multipath correction, the $P_G+L_G$ series is more consistent, indicating smaller STD. For this arc, the value of $⟨P_G+L_G⟩$ and its STD is 101,425,887.203 ± 0.620 m; however, with multipath model correction, the average and STD is 101,425,887.328 ± 0.438 m. This means that, in this arc, the multipath errors result in 0.12 m (≈0.4 ns) bias for the leveled phase observables. The improvement in STD indicates that $⟨P_G+L_G⟩$ is calculated with higher confidence level. For all the arcs on DOY 062/2015, $⟨P_G+L_G⟩$ and their STD values are calculated per each arc; and the average $⟨P_G+L_G⟩$ differences between the uncorrected and corrected cases as well as average STD improvement are shown in

Table 1. Statistics of $⟨P_G+L_G⟩$ average differences between leveled carrier-to-code observables with and without multipath errors corrected and the STD improvement. All the arcs in DOY 062/2015 are taken into account and the values are averaged for GPS, BDS GEOs, BDS IGSOs and BDS MEOs.

	$⟨{\mathit{P}}_{\mathit{G}}+{\mathit{L}}_{\mathit{G}}⟩$ Difference (m)	Raw ${\mathit{P}}_{\mathit{G}}+{\mathit{L}}_{\mathit{G}}$ STD (m)	Corrected ${\mathit{P}}_{\mathit{G}}+{\mathit{L}}_{\mathit{G}}$ STD (m)
GPS	−0.025	0.518	0.499
BDS-GEO	−0.047	0.727	0.708
BDS-IGSO	−0.040	0.678	0.652
BDS-MEO	−0.066	0.553	0.493

grouped as GPS, BDS GEOs, BDS IGSOs and BDS MEOs. The $⟨P_G+L_G⟩$ values are overall smaller after multipath errors are corrected. For GPS, BDS GEOs and BDS IGSOs, the average decrement in leveled phases are 0.025, 0.047 and 0.040 m after multipath correction, respectively, and their STDs are about 0.02 m (≈0.07 ns) smaller. The BDS MEOs show most impressive improvement: it could reduce the average STD by about 0.060 m while the leveled phases are biased by −0.066 m, which is due to its large multipath errors as indicated in Figure 1.

The DCB estimation results on DOY 062/2015 for the two cases are assessed. The formal errors of the DCB estimates for the two cases are firstly examined and both can reach far better than 0.01 ns, which is too good to be realistic. This is due to huge data amount: over 330,000 observations are used in the daily estimation, which significantly reduced the posterior variance.

To validate our DCB estimations, the GPS and BDS DCB solutions from the IGS Multi GNSS Experiment (MGEX) Ionosphere ACs DLR (Deutschen Zentrums für Luft-und Raumfahrt) [17] and IGG of CAS (Institute of Geodesy and Geophysics, Chinese Academy of Science) [32] are adopted for comparison purpose. Generated using ground-based observations with different algorithms, the DLR and CAS products can be generally regarded as the best quality and are used as reference values in many studies.

The DCB estimates are compared to the DLR and CAS results and the RMS values of their differences are shown in

Table 2. RMS of GPS and BDS DCB differences between the FY3C-based results and the DLR and CAS solutions.

	GPS Difference (ns)		BDS Difference (ns)
	DLR	CAS	DLR	CAS
L4-raw	0.385	0.483	0.458	0.556
L4-cor	0.380	0.472	0.451	0.545

. For comparative analysis, the RMS values of GPS and BDS DCB differences between DLR and CAS solutions on DOY 062/2015 are also calculated, which are 0.290 ns and 0.357 ns, respectively. The precision of GPS DCBs is about 0.390 ns when compared to DLR and 0.488 ns to CAS, while the BDS precision is 0.475 ns to DLR and 0.555 ns to CAS. From this comparison, the FY3C-based DCB results show larger discrepancies than the ground-based solutions, but are still in the order of magnitude of the DCB precision. The GPS DCB precision is a little degraded than those obtained using Jason-2 satellite from [21], even though in the Jason-2 solution only pseudorange observables are used. This might be attributed to the better Jason-2 data quality and DCB estimates are mainly dependent on pseudoranges. The FY3C-based GPS DCBs are closer to the DLR solutions than CAS with smaller RMS, which is in agreement with the findings in [21]. This can be partly explained as the VTEC parametrization may be not accurate enough in the CAS DCB algorithm, which is described as the sum of a two-dimensional polynomial of geographic latitude and local time and a finite Fourier series of local time with the degrees of 2 and 4, respectively. The other possible reason is that the DLR solutions provide direct GPS C1C-C2W DCBs but CAS do not; chaining of the individual C1C-C1W and C1W-C2W DCBs from CAS products could result in precision loss through error propagation. The RMS differences between different cases are small, but improvements are still observed when multipath is corrected.

It should be noted that neither DLR nor CAS solutions considers BDS satellite-induced multipath errors, which would increase the BDS satellite DCB estimation errors; hence, the conclusions of BDS DCB evaluation from this direct comparison could be distorted when treating the DLR or CAS BDS solutions as reference truth. It is more proper to evaluate the precision of BDS satellite DCBs through stability analysis, which will be given in the subsequent section.

4. Evaluation of DCB and TEC Estimates with Multipath Correction

The FY3C-based GPS and BDS satellite DCB estimation results are first given. The estimated DCBs are then evaluated in different aspects, i.e., monthly stability and external comparison with respect to DLR and CAS solutions. Before the evaluation and comparison between different DCB sets, the DCB solutions are first aligned to a common reference of themselves using the DCB alignment algorithm proposed by [15], which is realized by constraining the sum of DCBs of a common set of satellites to a same value throughout the entire month. After alignment, the DCB stability for each satellite is calculated as the monthly standard deviation. The FY3C onboard receiver DCBs are also evaluated. Finally, the estimated VTEC is analyzed.

4.1. GPS and BDS Satellite DCB Evaluation

The GPS and BDS satellite DCBs are estimated daily based on the FY3C measurements under the designed two cases. Firstly, their daily estimated values are illustrated in Figure 3, where the L4-cor results are depicted as an example. As there are no G08, G26 or C13 measurements from onboard FY3C GNOS receiver during the entire month, their DCBs are not estimated. The estimated GPS and BDS DCBs are mainly in the range between −10 and 10 ns, except that C01 reaches a maximum around 15 ns. Their variations are all very stable with fluctuations generally within 1 ns.

The monthly stability is calculated as the STD values of the aligned daily DCBs. In Figure 4, the monthly DCB stability is shown for every GPS and BDS satellite for each case. For the GPS DCBs, their stability values vary between 0.05 ns and 0.2 ns and generally fluctuate around 0.1 ns, which is quite close to the ground-based results as shown in [32]. The results for BDS satellites are a little complex. Large STDs are found for the GEO satellites with values all greater than 0.15 ns, especially for C02 and C05. The reason for this is still under investigation. For BDS IGSOs and MEOs, their STDs are all around 0.10 ns, while C07 stability even reaches as small as 0.051 ns in L4-cor case, which might be because there are more C07 data than for the other BDS satellites during March 2015. Generally, the BDS IGSOs DCB stability is the best, followed by the MEOs, and GEOs the worst. In comparison, the typical BDS C2I–C7I DCB stability from ground-based solutions is at the order of 0.20–0.30 ns for GEOs, 0.10–0.15 ns for IGSOs and 0.15–0.20 ns for MEOs as reported in various studies (e.g., [32,33,34]). This indicates the FY3C-based BDS DCB stability shows similar performance compared to the ground-based results except C02 and C05; in particular, the FY3C-based MEOs stability is superior in both L4-raw and L4-cor cases.

The average DCB stability for GPS and BDS GEOs, IGSOs and MEOs is given in

Table 3. Stability statistics of GPS and BDS satellite DCBs from the L4-cor, L4-raw, DLR and CAS solutions. The BDS satellites are categorized as BDS GEOs, BDS IGSOs and BDS MEOs; and their respective stabilities are calculated. The L4-cor and L4-raw BDS GEOs stability values in the brackets represent the results without C02 and C05 satellites.

	L4-Cor (ns)	L4-Raw (ns)	DLR (ns)	CAS (ns)
GPS	0.106	0.112	0.058	0.112
BDS GEOs	0.200 (0.149)	0.206 (0.150)	0.148	0.169
BDS IGSOs	0.086	0.102	0.084	0.180
BDS MEOs	0.106	0.123	0.187	0.209

. As can be seen, the GPS average stability is at the order of 0.110 ns in both L4-raw and L4-cor cases. The BDS GEOs DCBs show largest variations with a monthly stability of 0.200 ns for L4-raw case and 0.206 ns for L4-cor case; by excluding C02 and C05, the GEOs stability is about 0.150 ns for both cases. The BDS IGSOs and MEOs DCB stabilities are both around 0.10 ns and show significant differences between L4-cor and L4-raw cases, which indicates the multipath effects on the DCB estimation.

By comparing L4-raw case to L4-cor solutions, the multipath effect on DCB stability is investigated. The improvement percentages are calculated and shown in Figure 5. In Figure 4 and Figure 5, it can be observed that generally the DCB stability is improved for most satellites. For GPS, over half of the GPS satellites exhibit STD reduction over 5%; however, large STD deterioration is also found for G01, G17 and G29, over 20%, indicating stability degradation. The STD of GPS L4-cor case is reduced by about 0.006 ns on average against L4-raw. For BDS, only C01 and C04 indicate stability decrement, 2.5% and 9.0%, respectively, while all others show positive improvement. The BDS GEOs DCB precision show rather small improvement. Dramatic STD decrement is found in BDS IGSOs and MEOs, reaching 0.059 ns for C08 at maximum. C08 obtains the largest stability enhancement (38.1%), followed by C07 (27.3%). The average improvement after applying multipath corrections are 5.0%, 3.1%, 16.2% and 13.6% for GPS, BDS GEOs, IGSOs and MEOs, respectively. The IGSOs and MEOs stability improves more since they suffer from larger systematic multipath errors as indicated in Figure 1. In addition, the BDS IGSOs prevail with overwhelming data amount advantage.

As the L4-cor case provides better DCB stability, it is selected as the FY3C-based solution for the inter-solution comparison hereafter. For illustration, Figure 6 depicts the DCB day-to-day variation series from FY3C-based L4-cor case, DLR and CAS solutions for G11, C03 (BDS GEOs), C07 (BDS IGSOs) and C11 (BDS MEOs) as an example. The day-to-day variations are calculated by subtracting the monthly-average DCB values after alignment. For G11, DCBs from DLR show best consistency with variations less than 0.2 ns during the whole month, while FY3C-based L4-cor and CAS DCBs give similar performance with larger scatters. All three solutions show large variations of C03 DCB. However, the FY3C-based C07 and C11 DCBs obviously outperform the DLR and CAS results with much smaller fluctuations.

For stability comparison, the DCB monthly STDs from DLR and CAS products are also calculated for each satellite during March 2015, respectively, and are depicted in Figure 4 together with the L4-cor solution. A notably superior GPS DCB estimation precision is found for DLR solution with STDs around 0.06 ns, and all less than 0.10 ns; in contrast, the CAS solution performance is comparable to L4-cor, with stability variations around 0.10 ns. Comparison results of BDS DCBs are mixed. The CAS and DLR solutions both reveal more stable estimates among different BDS GEOs than the FY3C-based solution, with STDs between 0.10 and 0.20 ns; and, additionally, their C02 and C05 stability is in the same magnitude of other BDS GEOs, indicating more realistic estimates. The FY3C-based C01, C03 and C04 stability actually achieves at a similar level with respect to the DLR and CAS solutions. The CAS solution of BDS IGSOs suffers from high fluctuations with stability above 0.15 ns; however, the DLR- and FY3C-based solutions are both below 0.10 ns and their STDs are very close. For BDS MEO DCBs, the STD is around 0.20 ns for DLR and CAS, and is improved to near 0.10 ns with the FY3C solutions.

The average stability of DLR and CAS solutions is given in Table 3. The GPS DCB stability from DLR solution outperforms the others with an average STD 0.058 ns, while that from FY3C-based L4-cor and L4-raw cases as well as the CAS solutions all shows small degraded performance around 0.11 ns. The FY3C L4-cor case provides the best stability for MEOs DCBs of 0.106 ns, while its IGSO DCBs stability (0.086 ns) is at the same level with respect to DLR solution (0.084 ns); both are superior to CAS (0.180 ns). The GEO stability is quite close between DLR and CAS solutions, which are 0.148 ns and 0.169 ns, respectively, while FY3C-based solution performs the worst. It is worth mentioning that, after excluding C02 and C05, the FY3C-based GEO DCBs stability can also achieve 0.150 ns, very close to the DLR and CAS solutions.

The above evaluations reveal that the FY3C-based GPS and BDS DCB estimates can reach similar stability level as the DLR and CAS products, or even better, except for the abnormal performance of C02 and C05 satellites. Actually, we expect more from the FY3C onboard observations as they have two main advantages over the ground-based measurements, especially for the BDS GEOs. Firstly, as the LEO satellites such as FY3C usually orbit above the Earth at the altitude of plasmasphere where electrons are much sparser than the ionosphere, the onboard observables suffer from smaller signal delays due to electronic refraction along the transmission path than the ground-based observations. Thus, the influence of STEC on DCB calibration would be alleviated. Secondly, the BDS GEOs DCBs usually show significant larger variations compared to IGSOs and MEOs as reported in previous studies utilizing ground-based measurements (e.g., [17,32,33,34]), which is usually explained as the poor geometry between BDS GEOs and ground stations. However, with an onboard receiver, the geometry constraints of GEOs observations can be greatly enhanced due to the fast movement of LEO satellite. In our experiment, the GEOs DCB stability is barely improved as expected with respect to the ground-based DLR and CAS solutions. The reason could be attributed to the FY3C data quality problem. As mentioned in [31], the tracking capacity of FY3C onboard GNOS receiver only allows for six BDS satellites for navigation at the same time at its maximum. This severely reduces the BDS data amount, especially in the Asia-Pacific Ocean regions, where most of the BDS GEOs and IGSOs are mainly visible. In addition, large BDS data gaps are found on various days during March 2015, which would also affect the BDS DCB estimation stability. However, the FY3C-based GEOs DCB stability already outperforms the results reported in several studies (e.g., [33,34]). In comparison to other LEO-based GPS DCB estimates, the FY3C-based GPS stability is superior to that from COSMIC [20], and is comparable to the Jason-2 based result [21].

4.2. Receiver DCBs Evaluation

The DCB of GNOS receiver onboard FY3C is also evaluated. The GNOS GPS and BDS DCB is illustrated in Figure 7. For comparison, the L4-cor and L4-raw cases are both shown after aligning to their own solutions. Figure 7 reveals that the GPS DCB of GNOS receiver varies between −50 and −48 ns which is actually quite large in size while BDS DCB between 10 and 11 ns. Generally, the L4-cor case exhibits smaller variations and are more consistent. The GNOS GPS and BDS average DCB values and their STDs are calculated and shown in

Table 4. Monthly-mean DCBs of FY3C GNOS receiver and stabilities from L4-raw to L4-cor solutions.

	L4-Raw (ns)		L4-Cor (ns)
	GPS	BDS	GPS	BDS
Mean (ns)	−48.784	10.715	−49.174	10.367
STD (ns)	0.254	0.233	0.242	0.209

. A bias of 0.390 ns is found between L4-cor and L4-raw cases for the GNOS GPS DCB, while 0.348 ns for BDS DCB, which is attributed to multipath errors. The DCB estimates are more stable after multipath errors are corrected. The obtained FY3C GNOS receiver DCB and its stability is −49.174 ± 0.242 ns for GPS signals and 10.367 ± 0.209 ns for BDS signals.

4.3. Plasmasphere Electron Content Estimation

During DCB estimation with the FY3C onboard GNSS data, the VTEC along the FY3C tracks is also obtained. As they refer to the electron content along the vertical path from FY3C satellite to the GNSS satellite altitude, these obtained VTEC parameters can also be considered as the plasmasphere TEC (pTEC). To investigate the realism of these pTEC estimates, they are compared to the pTEC predictions generated by the IRI 2016 model [35] along the FY3C orbit. To obtain better performance of IRI 2016 outputs, standard options of IRI 2016 are used, except that the Shubin-COSMIC model is selected for hmF2 calculation.

Figure 8 illustrates the pTEC estimates and predictions during DOY 070/2015 as a case study. The solar activity proxy F10.7 on 070/2015 is 126.7 sfu and the 3-h geomagnetic K-indices values are between 1 and 3, indicating active solar activity and quiet geomagnetic condition. Compared to the pTEC estimates in [21] from Jason-2 pseudoranges, our FY3C-based pTEC series are much smoother and less noisy, which could be mainly attributed to the use of high-precise carrier-phase observations. In comparison to other studies exhibiting larges jumps in the estimated pTEC series when averaging observations from multiple GNSS satellites [36], our pTEC series performs continuously due to the applied loose constraints on consecutive VTEC parameters. Very similar variation patterns are revealed between the FY3C and IRI pTEC series in general. At the altitude of FY3C orbit, the pTEC fluctuates mainly under 6 TECU during the selected period, and exhibits large latitude-dependent signatures as significant increment with decreasing latitudes, which is in agreement with the plasmasphere morphology. For each orbit revolution, the pTEC reaches its maximum values around the equator at around LT (local time) 1000 or LT 2200. It is also noticed that, rarely, the FY3C-based pTEC drops below 0 at high latitudes, which is unrealistic and might be due to overestimation of DCB parameters.

One major difference shown in Figure 8 is that the IRI pTEC predictions produce two significant maxima when FY3C crosses around the equator, which is identified as the equatorial ionization anomaly (EIA) phenomenon, while the FY3C-based pTEC hardly shows such features. This discrepancy, as also found in the other LEO-based plasmasphere observations [37,38], originates from the misrepresentation of the NeQuick topside profile adopted in the IRI model. Usually, the dual-peak structure of EIA is present at height below 600 km, and then gradually merges to form one peak with increasing altitude [39]. However, the NeQuick topside model tends to extrapolate the double-hump structure to higher latitudes without the merging into one peak [37,40]. pTEC enhancement in the polar zones is observed from the FY3C-based estimates while the IRI 2016 model does not produce such subtle structure. This is in accord with the CHAMP and GRACE results, which indicated underestimation of IRI model in the polar regions during March under similar solar and geomagnetic conditions, and might be due to that IRI modeling efforts have been focused on the electron density in non-auroral ionosphere [37]. The average difference between the IRI 2016 and FY3C-based pTEC on DOY 070/2015 is −0.32 ± 0.95 TECU. The above comparison indicates that the estimated VTEC during the DCB calibration is capable of capturing the real pTEC variations and is generally numerically realistic.

Another case is also investigated to examine the pTEC responses to a great geomagnetic storm occurred on DOY 076/2015 (17 March 2015). Represented by the 3-h standardized Kp as well as the 1-h Dst index, the geomagnetic activity during DOY 075–077/2015 (16–18 March 2015) period is depicted in Figure 9. During this period, the Kp indices first showed rather quiet geomagnetic activity (<4-) on 16 March; however, it increased sharply to 8- at its maximum on 17 March and then gradually reduced to a lower level around 5 on 18 March. This dramatic geomagnetic activity variation is also manifested in the Dst index which dropped to −228 nT at 23:00 on 17 March. The F10.7 proxy on 16–18 March was 116.0, 113.2 and 113.7 sfu, which indicated an active solar activity level.

The pTEC variations based on the IRI model and the FY3C estimates during this geomagnetic storm are shown in Figure 10. Only the dayside (LT 0930-1130) pTEC is displayed for better visualization of the pTEC variations. As indicated by Figure 10a,b, both the FY3C-based and IRI pTEC performed similarly with peak pTEC values about 6 TECU on 16 March when a quiet geomagnetic condition is assured. Strong responses to the geomagnetic storm could be observed in the FY3C-based pTEC estimates after UT (Universal Time) 1400 on 17 March (9 h after the storm commencement time) when the Kp indices increased to 8-. The pTEC eventually grew to about 10 TECU at its maximum and then stayed at a high level during UT 1400–2400 17 March. The largest enhancement due to this storm reached near 70%. After the main phase of the storm, the peak pTEC values reduced to less than 4 TECU, nearly half of the IRI pTEC; and then pTEC generally remained at a rather low level, even at mid and low latitudes, indicating large TEC depletions during the recovery phase of this storm. The above findings generally coincide with the results given in [41]. In contrast, the IRI 2016 model exhibits no significant variations during the geomagnetic storm. This again demonstrates the realistic physical meaning of our pTEC estimates, and implies that the estimated pTEC can contribute to plasmasphere studies.

5. Conclusions

A method for GPS and BDS DCB estimation utilizing space-borne GPS and BDS measurements from FY3C satellite is described in this study. The precise carrier-phase observations are exploited in the method through carrier-to-code leveling, while the impact of FY3C code multipath errors on the leveling algorithm as well as DCB estimation is analyzed. The FY3C pseudorange multipath errors are investigated by grid modeling in the satellite antenna reference frame, which indicates there are systematic multipath errors in both FY3C GPS and BDS code measurements, among which BDS MEOs observations suffer the largest. The FY3C code multipath grid correction model can improve the code precision by over 30% for BDS MEOs observables and 5–10% for GPS, BDS GEOs and IGSOs. The BDS IGSOs and MEOs DCB stability can be significantly improved after applying the constructed multipath grid model by 16.2% and 13.6%, respectively, while, for GPS and BDS GEOs, the improvements are comparatively smaller, reaching 5.0% and 3.1%, respectively.

The FY3C-based DCBs estimates are evaluated by stability statistics and comparison to the ground-based DLR and CAS DCB products. The monthly stability analysis manifests that the stability of FY3C-based DCBs are generally at the same order of magnitude of the conventional ground-based DCB solutions. After multipath correction, the FY3C-based GPS DCB stability archives 0.110 ns on average, which is at the same order with respect to CAS solution but is worse compared to DLR solution. The FY3C-based BDS IGSOs and MEOs DCBs show very impressive stability, reaching 0.086 ns and 0.106 ns, respectively; for the former, it is at the same level as DLR solution and better than CAS, while, for the latter, it outperforms both the DLR and CAS solutions with considerable improvement. Degraded stability of the FY3C-based C02 and C05 DCBs is observed reaching to over 0.250 ns, while the other BDS GEOs stability is in good agreement with the DLR and CAS solutions, reaching 0.150 ns level on average. The estimated GNOS receiver DCB and its stability during March 2015 is −49.174 ± 0.242 ns for GPS and 10.367 ± 0.209 ns for BDS.

The propose method can not only estimate the GNSS satellite and space-borne receiver DCBs but also the vertical plasmasphere TEC from the receiver to the altitude of GNSS satellites. Two case studies are conducted under quiet and storm geomagnetic conditions, respectively. Both cases indicate that the estimated pTEC values are generally physically and numerically realistic, and can provide insights into observing and monitoring the plasmasphere TEC variations.