Estimating satellite and receiver differential code bias using a relative Global Positioning System network

Precise total electron content (TEC) is required to produce accurate spatial and temporal resolution of global ionosphere maps (GIMs). Receivers and satellite differential code biases (DCBs) are one of the main error sources in estimating precise TEC from Global Positioning System (GPS) data. Recently, researchers have been interested in developing models and algorithms to compute DCBs of receivers and satellites close to those computed from the Ionosphere Associated Analysis Centers (IAACs). Here we introduce a MATLAB code called Multi Station DCB Estimation (MSDCBE) to calculate satellite and receiver DCBs from GPS data. MSDCBE based on a spherical harmonic function and a geometry-free combination of GPS carrier-phase, pseudo-range code observations, and weighted least squares was applied to solve observation equations and to improve estimation of DCB values. There are many factors affecting the estimated values of DCBs. The first one is the observation weighting function which depends on the satellite elevation angle. The second factor is concerned with estimating DCBs using a single GPS station using the Zero Difference DCB Estimation (ZDDCBE) code or using the GPS network used by the MSDCBE code. The third factor is the number of GPS receivers in the network. Results from MSDCBE were evaluated and compared with data from IAACs and other codes like M_DCB and ZDDCBE. The results of weighted (MSDCBE) least squares show an improvement for estimated DCBs, where mean differences from the Center for Orbit Determination in Europe (CODE) (University of Bern, Switzerland) are less than 0.746 ns. DCBs estimated from the GPS network show better agreement with IAAC than DCBs estimated from precise point positioning (PPP), where the mean differences are less than 0.1477 and 1.1866 ns, respectively. The mean differences of computed DCBs improved by increasing the number of GPS stations in the network.

1 Introduction

Total electron content (TEC) is an important parameter in the study of ionospheric dynamics, structures, and variabilities. The ionosphere is a dispersive medium for space geodetic techniques operating in the microwave band (Böhm and Schuh, 2013) that allows calculation of TEC using Global Positioning System (GPS) dual-frequency radio transmissions. The global availability of GPS has made it a valuable tool for monitoring regional and global Earth ionospheric activity (Hernández-Pajares et al., 1999; Komjathy et al., 2005; Li et al., 2015; Liu and Gao, 2004; Mannucci et al., 1993). Unfortunately, GPS-derived TEC measurements are adversely affected by an inherent inter-frequency bias within the receiver and satellite hardware, typically referred to as the differential code biases (DCBs). Careful estimation of the DCBs is required to obtain accurate TEC, which is used in several applications, such as in several ionospheric prediction models, and in the correction of GPS positioning measurements (McCaffrey et al., 2017). A number of methods have been proposed for the estimation of GPS receiver DCBs, each with varying requirements and limitations, including making assumptions about the ionospheric structure, the use of internal calibration (Arikan et al., 2008; Themens et al., 2013, 2015), or the use of a reference instrument or model. Estimating DCBs for receivers and satellites from GPS observations depends on two approaches, the relative and absolute methods. The relative method utilizes a GPS network, while the absolute method determines DCBs from a single station (Sedeek et al., 2017). In the current study, we applied a relative method to calculate DCBs of satellites and GPS receivers.

Figure 1 The flowchart shows how the code works.

There has also been growing interest in measuring the accuracy of these methods and how different factors, e.g., ionospheric activity, play a role in these methods (McCaffrey et al., 2017). Nowadays, reliable GIMs and accurate DCBs of satellites and the International GNSS Service (IGS) stations can be obtained from IAACs like CODE (Schaer, 1999), the European Space Agency (ESA, Germany) (Feltens and Schaer, 1998), the Jet Propulsion Laboratory (JPL, USA) (Mannucci et al., 1998), and UPC (Technical University of Catalonia, Spain) (Hernández-Pajares et al., 1999; Orús et al., 2005). However, the IAAC DCB receivers' values are only available for IGS stations. Furthermore, some of the IGS ground receivers' DCB estimates are not available from all analysis centers. Also, some regions do not have any IGS ground stations, like our country Egypt, which means the TEC values over them would be interpolated from the nearest calculated values. As TEC values are dependent on DCB values, a mathematical model is required to calculate DCBs from GPS data.

In this study we introduce a mathematical model estimating satellites and receiver DCBs for a GPS network based on a spherical harmonic function (SHF) written under the MATLAB environment; the developed mathematical model uses a geometry-free combination of pseudo-range observables (P code). Weighted least squares were used to consider variation of the satellite elevation angle. The code was evaluated and compared with other researchers' codes in the “Results and analysis” section. In the “Conclusion” section we summarize the overall paper results.

2 GPS observation model

For a GPS satellite, the pseudo-range and carrier-phase observations between a receiver and a satellite can be expressed as (Jin et al., 2008; Leandro, 2009; Leick et al., 2015; Zhang et al., 2018)

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{P}_{r,j}^s\left(i\right)& ={\mathit{\rho }}_r^s\left(i\right)+c(\mathrm{d}{t}_r-\mathrm{d}{t}^s)+T_r^s+I_{r,j,P}^s\\ & +{\mathrm{DCB}}_r^P-\mathrm{DC}{B}_s^P+M_j+E_j,\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathrm{\Phi }}_{r,j}^s\left(i\right)& ={\mathit{\rho }}_r^s\left(i\right)+c(\mathrm{d}{t}_r-\mathrm{d}{t}^s)+T_r^s-I_{r,j,\mathrm{\Phi }}^s+{\mathit{\lambda }}_j{N}_j\\ & +p{b}_{r,j}-\mathrm{p}{b}_{\mathrm{s},j}+{\mathrm{DCB}}_r^{\mathrm{\Phi }}-\mathrm{DC}{B}_s^{\mathrm{\Phi }}+m_j+e_j,\end{array}\end{array}$$

with r, s, j, and i the receiver, satellite, frequency, and epoch indices, and where $P_{r,j}^s\left(i\right)$ are pseudo-range measurements, in meters, ${\mathrm{\Phi }}_{r,j}^s\left(i\right)$ carrier-phase measurements, in meters, ${\mathit{\rho }}_r^s\left(i\right)$ the geometric distance between satellite and receiver antennas, in meters, c the speed of light, in meters per second, dt_r and dt^s receiver and satellite clock errors, respectively, in seconds, $T_r^s$ the neutral troposphere delay, in meters, $I_{r,j,P}^s$ and $I_{\mathrm{r},j,\mathrm{\Phi }}^s$ the ionosphere delay of pseudo-range and carrier-phase observations, in meters, N_j carrier-phase integer ambiguities, in cycles, λ_j carrier-phase wavelength, in meters, ${\mathrm{DCB}}_{\mathrm{r},}^p$ and ${\mathrm{DCB}}_{\mathrm{s}}^p$ receiver and satellite pseudo-range hardware delays, respectively, in metric units, ${\mathrm{DCB}}_{\mathrm{r}}^{\mathrm{\Phi }}$ and DCB${}_{\mathrm{s}}^{\mathrm{\Phi }}$ receiver and satellite carrier-phase hardware delays, respectively, in metric units, M_j the pseudo-range multipath, in meters, E_j other un-modeled errors of pseudo-range measurements, in meters, pb_r,i and pb_s,i receiver and satellite carrier-phase initial phase bias, respectively, in metric units, m_j the carrier-phase multipath, in meters, and e_j other un-modeled errors of carrier-phase measurements, in meters.

Figure 2 IGS station locations.

Here, we consider a measurement scenario where one GPS receiver tracks dual-frequency code and phase data from a total of m satellites over t epochs, thereby implying that r=1, s=1, … m, j=1, 2, and i=1, …, t.

Firstly, the code read the RINEX files and extracted the pseudo-range and carrier-phase observations which are the range distances between the receivers and satellites measured using L₁ and L₂ frequencies. The “geometry-free” linear combination of GPS observations is used to derive the observation equations. The geometric range, clock offsets, and tropospheric delay are frequency independent and can be eliminated using this combination. The “geometry-free” linear combinations for pseudo-range and carrier-phase observations are given as (Al-Fanek, 2013)

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{P}_{\mathrm{4}}& =P_{r,\mathrm{1}}^s\left(i\right)-P_{r,\mathrm{2}}^s(i{)}_{=}{I}_{r,\mathrm{1},p}^s-I_{r,\mathrm{2},p}^s\\ & +{\mathrm{DCB}}_r^p+{\mathrm{DCB}}_s^p+E_{\mathrm{12}},\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{4}}& ={\mathrm{\Phi }}_{r,\mathrm{1}}^s\left(i\right)-{\mathrm{\Phi }}_{r,\mathrm{2}}^s(i{)}_{=}{I}_{r,\mathrm{2},\mathrm{\Phi }}^s-I_{r,\mathrm{1},\mathrm{\Phi }}^s+{\mathit{\lambda }}_{\mathrm{1}}{N}_{\mathrm{1}}-{\mathit{\lambda }}_{\mathrm{2}}{N}_{\mathrm{2}}\\ & +{\mathrm{DCB}}_r^{\mathrm{\Phi }}+{\mathrm{DCB}}_s^{\mathrm{\Phi }}+e_{\mathrm{12}}.\end{array}\end{array}$$

$E_{\mathrm{12}}=\sqrt{(E_{\mathrm{1}}{)}^{\mathrm{2}}+(E_{\mathrm{2}}{)}^{\mathrm{2}}}$ is the combination of multipath and measurement noise on $P_{r,\mathrm{1}}^s\left(i\right)$ and $P_{r,\mathrm{2}}^s\left(i\right)$ (m), and $e_{\mathrm{12}}=\sqrt{(e_{\mathrm{1}}{)}^{\mathrm{2}}+(e_{\mathrm{2}}{)}^{\mathrm{2}}}$ is the combination of multipath and measurement noise on ${\mathrm{\Phi }}_{r,\mathrm{1}}^s\left(i\right)$ and ${\mathrm{\Phi }}_{r,\mathrm{2}}^s\left(i\right)$ (m).

Figure 3 Mean difference between the receiver DCB values of CODE and the computed values by each of MSDCBE and M_DCB estimated from 1 to 31 January 2010.

Figure 4 Mean differences between receiver DCB values estimated by MSDCBE and those released by CODE, JPL, and IGS combined from 1 to 31 January 2010.

To reduce the multipath and noise level in the pseudo-range observables, the carrier-phase measurements are used to compute a more precise relative smoothed range. Although the carrier-phase observables are more precise than the code derived, they are ambiguous due to the presence of integer-phase ambiguities in the carrier-phase measurements. To take advantage of the low-noise carrier-phase-derived and unambiguous nature of the pseudo range, both measurements are combined to collect the best of both observations.

Smoothed P_4,sm observations can be expressed as follows (Jin et al., 2012):

$$\begin{array}{}\text{(5)}& P_{\mathrm{4},\mathrm{sm}}={\mathit{\omega }}_t{P}_{\mathrm{4}}\left(t\right)+(\mathrm{1}-{\mathit{\omega }}_t)P_{\mathrm{4},\mathrm{prd}}\left(t\right)\left(t\mathit{>}\mathrm{1}\right),\end{array}$$

where t stands for the epoch number and ω_t is the weight factor related to epoch t, and

$$\begin{array}{}\text{(6)}& P_{\mathrm{4},\mathrm{prd}}\left(t\right)=P_{\mathrm{4},\mathrm{sm}}(t-\mathrm{1})+\left[L_{\mathrm{4}}\left(t\right)-L_{\mathrm{4}}(t-\mathrm{1})\right]\left(t\mathit{>}\mathrm{1}\right),\end{array}$$

when t is equal to 1, which means the first epoch of one observation arc, and P_4,sm is equal to P₄.

Figure 5 Mean difference between the receiver DCB values of the IGS and the computed values by MSDCBE estimated from 1 to 5 January 2010.

Figure 6 Mean difference between the satellite DCB values of the IGS and the computed values by MSDCBE estimated from 1 to 5 January 2010.

3 Spherical harmonic model

To determine the receiver DCB, there are two different methods. The first one is to calibrate the receiver device and obtain the DCB directly. This method calculates the DCB of the receiver device, ignoring that from the antenna cabling used during observation (Hansen, 2002). The second method calculates the receiver DCB as a part of the GPS signal time delay which is independent of the type of antenna. MSDCBE code works like the second method (Fig. 1).

The ionosphere delay can be expressed as follows (Abid et al., 2016):

$$\begin{array}{}\text{(7)}& d_{\mathrm{ion}}={\displaystyle \frac{\mathrm{40.3}}{f^{\mathrm{2}}}}\mathrm{STEC},\end{array}$$

where f stands for the frequency of the carrier and slant total electron content (STEC) is the total electron content along the path of the signal. The observation equation can be formed by substituting Eq. (7) into Eq. (3) and replacing P₄ by smoothed P_4,sm; we get (Abid et al., 2016)

$$\begin{array}{}\text{(8)}& P_{\mathrm{4},\mathrm{sm}}=\mathrm{40.3}({\displaystyle \frac{\mathrm{1}}{f_{\mathrm{1}}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{f_{\mathrm{2}}^{\mathrm{2}}}})\mathrm{STEC}+c\cdot {\mathrm{DCB}}_r+c\cdot {\mathrm{DCB}}_s,\end{array}$$

where c is the speed of light and DCB_r and DCB_s are differential code biases for receivers and satellites in seconds.

Figure 7 Mean difference between the receiver DCB values of the IGS and the computed values by each of M_DCB, ZDDCBE, and MSDCBE estimated from 1 to 5 January 2010.

Figure 8 Mean difference between the satellite DCB values of the IGS and the computed values by M_DCB, ZDDCBE, and MSDCBE estimated from 1 to 5 January 2010.

STEC can be translated into vertical total electron content (VTEC) using the modified single-layer model (MSLM) (Haines, 1985; Jin et al., 2012):

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathrm{VTEC}=\mathrm{MF}\left(z\right)\mathrm{STEC},\text{(10)}& {\displaystyle }\mathrm{MF}=\cos \left(\arcsin \left({\displaystyle \frac{R}{R+H}}\sin \left(\mathit{\alpha }z\right)\right)\right),\end{array}$$

where MF is the mapping function, z is the satellite elevation angle, R is the radius of the Earth (6371 km), H is the attitude of the ionosphere thin shell (assumed as used by CODE = 506.7 km), and α=0.9782 (Jin et al., 2012).

To estimate the satellite and receiver HDs, the current study applies a model based on a spherical harmonic function to calculate them using zero-difference observations. The used model is expressed as follows (Schaer, 1999; Li et al., 2015; Elghazouly et al., 2019):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{VTEC}(\mathit{\beta },s)& ={\sum }_{n=\mathrm{0}}^N{\sum }_{m=\mathrm{0}}^n{P}_n^m\left(\sin \right(\mathit{\beta }\left)\right)\\ & \left(A_n^m\cos \left(m\mathit{\lambda }\right)+B_n^m\sin \left(m\mathit{\lambda }\right)\right),\end{array}\end{array}$$

where β is the geocentric latitude of IPPs (ionosphere Pierce points), s is the solar fixed longitude of IPPs, N is the degree of the spherical function, and M is the order of the spherical harmonic function; fourth order is used. P_mn is the regularization Legendre series and A_mn and B_mn are the estimated spherical harmonics coefficients.

By substituting Eqs. (8), (10), and (11) into Eq. (9), we get

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \sum _{n=\mathrm{0}}^N\sum _{m=\mathrm{0}}^n{P}_n^m\left(\sin \right(\mathit{\beta }\left)\right)\left(A_n^m\cos \left(m\mathit{\lambda }\right)+B_n^m\sin \left(m\mathit{\lambda }\right)\right)\\ & =\cos \left(\arcsin \left({\displaystyle \frac{R}{R+H}}\sin (\mathit{\alpha }z\right)\right)\\ & \left[-{\displaystyle \frac{f_{\mathrm{1}}^{\mathrm{2}}{f}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{40.3}\left(f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}}\right)}}\left(P_{\mathrm{4},sm}-c\cdot {\mathrm{DCB}}_r-c\cdot {\mathrm{DCB}}_s\right)\right].\end{array}\end{array}$$

Only one GPS station has more than 20 000 observations per day. When applying Eq. (12) using station observation data, the number of equations is much more than the number of unknown coefficients. These coefficients were determined using the weighted least squares method. The general form of the weighted least square function can be expressed as (Ghilani and Wolf, 2012)

$$\begin{array}{}\text{(13)}& \mathit{X}=({\mathbf{A}}^{T}PA{)}^{-\mathrm{1}}{\mathbf{A}}^{T}PL,\end{array}$$

where X is the unknown parameter vector, namely ${\mathbf{A}}_n^m$, $B_n^m$, DCB_r, and DCB_s, A is the coefficient (design) matrix (coefficients of ${\mathbf{A}}_n^m$, $B_n^m$, DCB_r, and DCB_s), L is the observation vector (values of P_4,sm), and P is the weight matrix.

As is known, the quality of observations is affected by satellite elevation angle, and each observation has a weight value dependent on its satellite elevation angle. The weight value can be computed from the following Eqs. (14), (15), and (16) (Luo, 2013):

$$\begin{array}{}\text{(14)}& {\displaystyle }w={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}},\text{(15)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}={\left[\mathrm{0.05}+{\displaystyle \frac{\mathrm{0.02}}{\sin (z{)}^{\mathrm{2}}}}\right]}^{\mathrm{2}},\text{(16)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}=(f+d{)}^{\mathrm{2}},\end{array}$$

where f and d are two constants equal to 5 and 2 cm, respectively (Ray and Griffiths, 2008).

4 Mathematical model evaluation

The MSDCBE software was written in MATLAB (version 2016a). The first input is GPS observations in Receiver Independent Exchange (RINEX) format according to the selected stations (Fig. 2) downloaded from ftp://garner.ucsd.edu/rinex (last access: 8 August 2018) and precise ephemeride (SP3) files of test days downloaded from http://www.GPScalendar.com/index.html?year=2010 (last access: 8 August 2018). In addition, IONosphere Map EXchange Format (IONEX) files of the IGS, CODE, and JPL are downloaded – as threshold values – from ftp://cddis.gsfc.nasa.gov/GPS/products/ionex/ (last access: 8 August 2018).

In the present contribution, to evaluate the performance of the developed model, numerical case studies were performed. The main goals of the numerical case studies are to investigate three issues.

The first issue is to investigate the effect of applying weighted least squares instead of least squares on satellites and GPS receiver DCBs, and this is done by comparing results from MSDCBE which apply weighted least squares with the published results of M_DCB by Jin et al. (2012) and with those of IAAC.

BOGO, BRUS, GOPE, GRAS, ONSA, POTS, PTBB, SOFI, and WTZA IGS station data from 1 to 31 January 2010 were applied as it was the same network used by Jin et al. (2012).

The second issue is to investigate the correlation between the size (number of receivers) of the GPS network and estimated DCBs for satellite and GPS receivers, and this is done by comparing DCB values of three stations, namely GOPE, GRAS, and ONSA, estimated from a network consisting of a 3 GPS receiver and a network consisting of a 9 GPS receiver.

This study was applied using IGS station data from 1 to 5 January 2010 of six stations, namely BOGO, BRUS, GOPE, GRAS, ONSA, POTS, PTBB, SOFI, and WTZA.

The third issue is to investigate the congruence of DCBs estimated from absolute and relative methods with other IAAC, and this is done by comparing results from MSDCBE with the published results of ZDDCBE by Sedeek et al. (2017).

This study was applied using data from 1 to 5 January 2010 of six stations, namely GOPE, GRAS, ONSA, MADR, PTBB, and SOFI, which was the same network used by Jin et al. (2012) and Sedeek et al. (2017).

4.1 Comparison of multi-station test results from MSDCBE and M_DCB

The first evaluation made by this paper is the evaluation of a weight function. MSDCBE used a weight function depending on the satellite elevation angle as mentioned before. Table 1 shows the differences and rms between satellites and receivers estimated from 1 to 31 January 2010 using multiple GPS stations of both MSDCBE (weighted) and M_DCB (unweighted).

Table 1 The differences and rms between satellites and receivers estimated from 1 to 31 January 2010 using multiple GPS stations (MSDCBE and M_DCB minus CODE).

Table 2 Mean difference between the receiver DCB values computed by the IGS and the computed values by using single-station M_DCB, ZDDCBE, and multi-station MSDCBE estimated from 1 to 5 January 2010.

From the table one can see that the differences of MSDCBE-estimated satellite DCBs are less than 0.302 ns and the rms of all satellite DCB differences is less than 0.128, except G1, whose rms = 0.250. The maximum difference of MSDCBE-estimated receiver DCBs is 0.150 ns of receiver GOPE and the minimum is 0.045 ns of receiver SOFI (Fig. 3). The maximum rms of MSDCBE-estimated receiver DCBs is 0.125. On the other hand, M_DCB results show that receiver DCB biases are slightly larger than those for satellites, but most of them are less than 0.4 ns, except G1, whose DCB bias reaches 0.746 ns. The rms of all differences is lower than 0.3 ns (Jin et al., 2012). Figure 4 shows the mean differences between receiver DCB values estimated by MSDCBE and those released by CODE, IGS, and JPL combined from 1 to 31 January 2010. The figure shows that the results of MSDCBE are mostly closer to those of CODE than IGS and JPL. By comparing Fig. 4 with the corresponding chart published by Jin et al. (2012), it is clear that all differences between MSDCBE receivers' DCB results and between CODE, IGS, and JPL are less than those from M_DCB, except station GOPE, which is almost equal.

4.2 Effect of network size factor on DCB estimation

By using multi-station DCB estimation, the number of stations used will appear as a factor influencing DCB estimation. This test was done by comparing DCBs computed by MSDCBE of a network of three receivers, namely GOPE, GRAS, ONSA, and DCBs of the same receivers, but this time as a part of a network of nine receivers, namely BOGO, BRUS, GOPE, GRAS, ONSA, PTBB, SOFI, and WTZA. Figure 5 shows these results which demonstrate that using nine receivers gives more accurate DCBs. Also, the satellite DCB differences (Fig. 6) almost improved, but not like receiver DCBs, because satellite DCBs are small values compared with those of receivers.

4.3 Comparison of multi station from MSDCBE and single station from ZDDCBE and M_DCB test results

In this section the performance of a multi-station network against single-station DCB estimation will be evaluated. Table 2 shows the mean difference between the receiver DCB values computed by IGS and the computed values by each of M_DCB, ZDDCBE, and MSDCBE estimated from 1 to 5 January 2010. Figure 7 shows these results graphically and Fig. 8 shows the mean differences computed from M_DCB, ZDDCBE, and MSDCBE for GPS satellites. The results show a significant difference between multi-station network and single-station DCB estimation. The maximum difference between receiver DCB estimation using IGS and MSDCBE is 0.1477 ns of the MADR station, but it is 1.1866 and 0.7982 ns for M_DCB and ZDDCBE, respectively.

5 Conclusions

The current study proposes a new MATLAB code called MSDCBE able to calculate DCBs of GPS satellites and receivers. This code was compared with two other codes and evaluated using IAAC data and, from all the above, we can conclude the following.

1. The estimated DCB results are affected by using a weight function according to satellite elevation angle observations. In addition, results show good agreement with IGS, CODE, and JPL results than using multi-station estimation DCB without a weight function.

2. When using multi-station DCB estimation, the number of input stations influences the DCB results. However, it is recommended to enlarge the size of the used network, but it needs high computer requirements and much more analysis time (only one station has more than 20 000 observations per a day).

3. The most effective factor in DCB estimation is using a multi-station network instead of a single station that appeared from results which improved from 1.1866 and 0.7982 ns maximum DCB mean differences for M_DCB and ZDDCBE single-station analysis to 0.1477 ns for MSDCBE. So, using multi-station network DCB estimation – if available – is strongly recommended.