Sensitivity of Receiver Differential Code Bias Estimates to Mapping Functions

Jeremy Ovadia; David Champlin; Corey Pong,; Eric Altshuler

doi:10.33012/navi.723

PDF

Abstract

Accurate estimations of differential code biases (DCBs) are critical for producing absolute ionospheric measurements using global navigation satellite system observables. DCB estimation generally requires translating slant total electronic content (TEC) measurements to the vertical domain using a mapping function. Analyzing DCB estimates from regional modeling and global ionospheric maps (GIMs) over 4.5 years reveals significant DCB differences across different mapping functions, varying by a few nanoseconds. Decompositions of receiver DCB estimates show that for some mapping functions, such as the Jet Propulsion Lab extended slab function and the 350-km thin shell function, variations in DCB estimate over time can be largely accounted for by temperature and ionospheric activity. In contrast, other mapping functions, such as the 450-km thin shell function, exhibit significantly less variation over time. Our computational and analytical results suggest the importance of selecting an appropriate mapping function for accurate DCB and TEC estimation from GIMs and shell-based spatiotemporal models.

Keywords

1 INTRODUCTION

Estimating differential code biases (DCBs), i.e., the relative hardware delays for signal pairings of global navigation satellite systems (GNSSs), for satellites and receivers is often a key task in modeling the ionosphere (Sardón et al., 1994). Geometry-free combinations, which take the difference in GNSS observables and leave frequency-dependent effects, provide measurements for ionospheric refraction that are offset by DCBs, which often show little change over several days (Coco et al., 1991; Wilson & Mannucci, 1993), along with other errors, such as wind-up and multipath. When satellite and receiver DCBs are well estimated by a global or regional network of stations and removed, such combinations can provide absolute ionospheric measurements, also referred to as “truth” or “supertruth” (Walter et al., 2018).

The performance and availability of legacy L1 satellite-based augmentation systems (SBASs) rely heavily on the ionosphere and its modeling. The Minimum Operational Performance Standards (MOPS) for the L1 SBAS outlined by the Radio Technical Commission for Aeronautics (2020) specifies the grid ionospheric vertical delays and grid ionospheric vertical errors (GIVEs) to be broadcast, enabling single-frequency GNSS users to correct errors in positioning estimates caused by ionospheric delay, the largest source of position error (Milliken & Zoller, 1978). The Wide Area Augmentation System (WAAS), i.e., the Federal Aviation Administration’s SBAS covering Global Positioning System (GPS) users in North America, employs an undersampled threat model within its GIVE monitor that accounts for deviations from planarity in worst-case situations of insufficient ionosphere sampling (Altshuler et al., 2001; Sparks et al., 2001; Blanch et al., 2002). Absolute ionospheric measurements, which require DCB estimation, populate the GIVE undersampled threat model. If DCBs are not accurately estimated, then errors will manifest in the absolute ionospheric measurements, the GIVE undersampled threat model, and, in turn, unquantized GIVEs. When these GIVEs are quantized at one of only 15 values, with the six largest possible values being 3, 3.6, 4.5, 6, 15, and 45 m, a small DCB error can prompt a significant unwarranted increase or decrease in broadcast GIVEs, strongly influencing a user’s horizontal and vertical positioning error estimates. In turn, these estimates determine the availability of SBAS corrections, according to Lateral Precision with Vertical Guidance 200 specifications (WAAS Performance Standard, 2008).

DCBs estimated jointly from the computation of global ionospheric maps (GIMs) have been widely used to benchmark methods for DCB estimation (Choi et al., 2013; Wang et al., 2021; Reid et al., 2024). Since the first computation of these maps by Mannucci et al. (1998) and Schaer (1999), GIMs have provided estimates for total ionospheric content (TEC) concentrated at a particular ionospheric height, often 450 km, on a global scale based on observation data from hundreds of GNSS receivers and models that concurrently estimate DCBs (Hernandez-Pajares et al., 2009). Five analysis centers from the International GNSS Service (IGS) have regularly published GIMs with DCB estimates over the past several years (Johnston et al., 2017): the Jet Propulsion Laboratory (JPL), Center of Orbit Determination in Europe (COD), European Space Agency (ESA), Chinese Academy of Sciences (CAS), and Natural Resources Canada (EMR). There are notable differences in the methodology and performance of TEC maps across these five GIMs (Hernandez-Pajares et al., 2009; Wu et al., 2021; Li et al., 2023); however, there is a general agreement in DCBs estimated from different GIMs, often within fractions of a nanosecond of one another, outside of some notable deviation from the EMR GIM for certain receiver DCBs.

The variation in DCBs has been well studied over both short (Li et al., 2018; Liu et al., 2019) and long (Xiang et al., 2020; Wageeh et al., 2023) timescales, with variations being attributed to temperature (Coster et al., 2013; Liu et al., 2019; Zha et al., 2019), solar activity (Xiang et al., 2020), receiver type (Choi et al., 2010), flex power (Wu et al., 2024), and Doppler shift (Liu et al., 2019), among other factors.

A notable and somewhat controversial factor in observed DCB estimate variation is ionospheric variation, for which a potential relationship was first identified by Hernandez-Pajares et al. (2009). DCB estimate variation contributions from TEC variation were later quantified by Zhang et al. (2014) and described as induced “pseudo-DCBs.” The nature of this variation has been debated among researchers, with Zhong et al. (2016, 2020) arguing that satellite DCB variation could be explained by satellite replacement and not ionospheric variation whereas Choi et al. (2019, 2020) countered with explanations that cited ionospheric activity as a source of receiver DCB variation. To our knowledge, studies have only relied on DCB estimates from GIMs; no study has explored how observed correlations between DCB estimates and ionospheric activity may be sensitive to aspects of the models that estimate DCBs.

Mapping functions that translate measured ionospheric delays to comparable vertical measurements at 90° elevation play an important role in ionospheric models in which the vertical component of the ionosphere is reduced to a thin shell or a series of thin shells. Outside of voxel models that do not reduce the ionosphere dimensionality, most ionospheric models experience obliquity errors from such slant-to-vertical mappings, which can be significant error sources for TEC and DCB estimation (Sparks, 2013). Different types of mapping functions have been used in various settings: the MOPS for L1 SBAS uses a “thin shell” model at 350 km (Radio Technical Commission for Aeronautics, 2020); some GIMs, such as the JPL and ESA GIMs, use an “extended slab” model or an approximation to it; and other studies (Zus et al., 2017; Xiang et al., 2019) have leveraged the International Reference Ionosphere (IRI) model (Bilitza, 1990) for mapping function determi-nation. Previous studies by Wang et al. (2016) and Maruyama et al. (2021) have explored how the ionospheric height used for modeling and for the associated thin shell mapping function can influence errors in joint DCB and TEC estimation through spatiotemporal modeling.

In this paper, we use DCBs estimated by regional models alongside GIM DCB estimates to study their sensitivity to the mapping function used. First, a modeling framework is outlined for DCB estimation via piecewise linear spherical harmonic basis functions in a sun-fixed geomagnetic coordinate system. For validation, our models are applied to a simulated ionosphere to recapture randomly generated DCBs. The models are also applied to actual historical GNSS observables and compared with GIM DCBs. Receiver DCB estimates for a WAAS-centric network of stations are shown for a 4.5-year period, with noticeable differences for different mapping functions used for estimation. Estimated DCBs over time are fit using measured temperature and estimated TEC as an attempt to account for observed variations in estimated DCBs. Computational modeling is supplemented with mathematical analysis that reduces an estimate for the mean receiver DCB in a network of stations to a weighted mean of observable differences by making a sequence of assumptions.

2 MODEL FORMULATION

In this section, we lay out mathematical modeling frameworks for DCB estimation and expressions that will aid in analyzing the sensitivity of receiver DCBs to mapping functions. First, a formulation of absolute ionospheric measurements in the slant domain is derived, followed by a presentation of three mapping functions that translate the measurements to the vertical domain. Next, a spatiotemporal model utilizing spherical harmonics in the vertical domain is reduced to a linear system for computational estimation of DCBs along with TEC. To supplement the computational model, analysis on the grounds of crude modeling assumptions is applied to derive mathematical expressions for mean receiver DCB estimates. Finally, a multifactor model that attempts to decompose receiver DCB estimate variation using both temperature and ionospheric activity is outlined.

2.1 Absolute Ionospheric Measurements

Consider two GNSS signals, A and B, broadcast by all satellites in a single constellation on two distinct frequencies, f_A and f_B, respectively. We will consider signal A to be the primary signal for a constellation, such as the GPS L1 coarse acquisition (C/A) signal, and signal B to be a secondary signal, such as the GPS L2(P) or L5 signal. Pseudoranges on signals A and B for the k-th observable, P_A,k and P_B,k, can be modeled for each frequency measured from a single satellite as follows:

$P_{A, k} = ρ_{k} + d_{c l o c k, k} + τ_{s l a n t, A, k} + γ_{r x, A, k} + γ_{s v, A, k} + ε_{A, k}$ 1

$P_{B, k} = ρ_{k} + d_{c l o c k, k} + τ_{s l a n t, B, k} + γ_{r x, B, k} + γ_{s v, B, k} + ε_{B, k}$ 2

Here, ρ represents the geometric range between the satellite and receiver, and d_clock is the combined clock bias of the receiver and satellite. γ_rx is the receiver code bias, and γ_sv is the satellite code bias for a given signal, denoted by its subscript. ε represents a combination of solid tide correction, ocean tide correction, observation noise along with multipath errors, and other unmodeled noise. τ_slant represents frequency-dependent slant ionospheric delays with the corresponding signal specified by the subscript.

Differencing the pseudoranges in Equations (1) and (2) yields the following:

$δ_{k} : = P_{B, k} - P_{A, k} = τ_{s l a n t, B, k} - τ_{s l a n t, A, k} - Γ_{r x, k} - Γ_{s v, k} + {\hat{ε}}_{k}$ 3

where $Γ_{r x} = γ_{r x, A} - γ_{r x, B}$ and $Γ_{s v} = γ_{s v, A} - γ_{s v, B}$ represent the DCB for the receiver and satellite, respectively, and $\hat{ε} = ε_{B} - ε_{A}$ is assumed to be zero-mean Gaussian-distributed noise.

The signal-dependent ionospheric delay can be modeled as inversely proportional to the square of the frequency of the broadcast signal:

$τ_{s l a n t, k} = k \frac{10^{16} S_{k}}{f^{2}}$ 4

Here, S is the slant TEC in TEC units (TECU; 1 TECU = 10¹⁶ electrons/m²) of the signal delay as it passes through the ionosphere, and f represents the frequencies at which the signals are broadcast. The proportionality constant $κ = q^{2} / (8 π^{2} m_{e} \in_{0}) = c^{2} r_{e} / (2 π) \approx 40.308193 m^{3} / s^{2}$ is described in terms of the charge (q), mass (m_e), and radius (r_e) of an electron, the permittivity of free space (ϵ₀), and the speed of light (c).

The pseudorange difference for signals A and B can be expressed in terms of the slant TEC as follows:

$δ_{k} = κ 10^{16} (\frac{1}{f_{B}^{2}} - \frac{1}{f_{A}^{2}}) S_{k} - Γ_{r x, k} - Γ_{s v, k} + {\hat{ε}}_{k}$ 5

From here, the absolute slant TEC measurements, or ionospheric “truth,” can be written as follows:

$S_{k} = ξ_{A, B} (δ_{k} + Γ_{r x, k} + Γ_{s v, k} - {\hat{ε}}_{k})$ 6

where:

$ξ = \frac{f_{A}^{2} f_{B}^{2}}{κ 10^{16} (f_{A}^{2} - f_{B}^{2})}$ 7

2.2 Slant-to-Vertical Mapping Functions

Ionospheric estimates in the slant domain can be translated into the vertical domain for modeling and standardized interpretation by first defining an ionospheric height that may be fixed or a function of time and/or latitude and longitude. The location at which the line-of-sight vector from the receiver to the satellite intersects with the assumed surface height of the ionosphere is considered as the measurement’s ionospheric pierce point (IPP), as shown in Figure 1.

FIGURE 1

Schematic diagram illustrating slant and vertical domains for TEC measurements

For a given slant TEC measurement S_k, let V_k be the associated vertical TEC (VTEC) to which a user would be subject for overhead GNSS measurements passing through the IPP at an elevation of 90°. We define the slant-to-vertical mapping function, G, as follows:

$G_{k} = \frac{V_{k}}{S_{k}}$ 8

We note that mapping functions are often expressed in vertical-to-slant form (F =1/G), but the functions are expressed inversely here. Using Equation (6), we can express the VTEC as follows:

$V_{k} = ξ G_{k} (δ_{k} + Γ_{r x, k} + Γ_{s v, k} - {\hat{ε}}_{k})$ 9

The definition of G as a simple ratio is straightforward; however, accurately modeling the slant-to-vertical mapping function can be very challenging. For example, it is possible for two different slant TEC measurements from different receiver–satellite pairs to share the same IPP but have significantly different G values owing to differences in slant paths.

To study the sensitivity of DCB estimation, we will consider three different mapping function models, outlined in Sections 2.2.1–2.2.3.

2.2.1 Thin Shell Model

The most commonly used mapping function is the thin shell approximation, which assumes that all ionospheric plasma is concentrated at a fixed ionospheric height, h_I, essentially reducing the ionosphere from three dimensions to a two-dimensional surface. The thin shell slant-to-vertical mapping, G_TS, can be expressed as follows:

$G_{T S} = \sqrt{1 - {(\frac{R_{e} \cos E}{R_{e} + h_{1}})}^{2}}$ 10

where E is the line-of-sight elevation angle and R_e is the Earth’s effective radius below the associated IPP. Here, we assume a fixed R_e = 6378136.3 m, following the Radio Technical Commission for Aeronautics (2020).

2.2.2 Extended Slab Model

The extended slab model assumes an electronic density that is uniform within a “slab” defined by two ionospheric heights and decays as heights trend away from the slab. Following Coster et al. (1992) and Mannucci et al. (1998), we express the normalized electron density of the ionosphere as a function of the radial distance from the Earth’s origin, r, as follows:

$N_{e, E S} (r) = {\begin{matrix} 0, & r > r_{t o p} \\ \frac{1}{1 + \frac{r - r_{2}}{p}}, & r_{2} \leq r \leq r_{t o p} \\ 1, & r_{1} \leq r < r_{2} \\ \frac{1}{1 + \frac{r_{1} - r}{q}}, & r_{b o t t o m} \leq r < r_{1} \\ 0, & r < r_{b o t t o m} \end{matrix}$ 11

Here, the density is constant in the slab, defined by the radial heights r₁ and r₂, whereas the parameters p and q determine the decay profile outside of the slab. The radial heights r_bottom and r_top define an assumed top and bottom of the ionosphere, respectively. The turquoise curve in the top panel of Figure 2 shows the extended slab electronic density, following the parameterization outlined by Mannucci et al. (1998) in what is commonly referred to as the “JPL extended slab.”

FIGURE 2

Electronic densities as a function of atmospheric height (top) and vertical-toslant mapping functions versus elevation (bottom) for the thin shell, JPL extended slab, and IRI-based mapping functions

The extended slab mapping function can be determined by analytically integrating the density function along vertical and slant paths for the respective delay approximations, ${\hat{V}}_{E S}$ and ${\hat{S}}_{E S}$ , which yields the following:

$G_{E S} = \frac{{\hat{V}}_{E S}}{{\hat{S}}_{E S}}$ 12

${\hat{V}}_{E S} = q \log (1 + \frac{r_{1} - r_{b o t t o m}}{q}) + r_{2} - r_{1} + p \log (1 + \frac{r_{t o p} - r_{2}}{p})$

${\hat{S}}_{E S} = q g (- r_{b o t t o m}, - r_{1}, r_{F}^{2}, q + r_{1}) + \sqrt{r_{2}^{2} - r_{F}^{2}} - \sqrt{r_{1}^{2} - r_{F}^{2}} + p g (r_{2}, r_{t o p}, r_{F}^{2}, p - r_{2})$ 13

where:

$r_{F} = R_{e} \cos E$ 14

$g (y, z, b, c) = \frac{1}{2} \log (\frac{u_{1} (z, b)}{u_{1} (y, b)}) + c \log (\frac{u_{2} (z + c, b, c)}{u_{2} (y + c, b, c)})$ 15

$u_{1} (x, b) = \frac{\sqrt{x^{2} - b} + x}{\sqrt{x^{2} - b} + x}$ 16

$u_{2} (x, b, c) = \frac{c^{2} - b - c x + \sqrt{c^{2} - b} \sqrt{x^{2} - 2 c x + c^{2} - b}}{x}$ 17

The turquoise curve in the bottom panel of Figure 2 shows the inverse of G_ES as vertical-to-slant mapping using the JPL extended slab parameterization. It is important to note that the JPL extended slab curve closely follows the corresponding 550-km thin shell curve in the plot, suggesting a close similarity between the two models, as previously observed by Sparks et al. (2000).

2.2.3 IRI-Based Mapping

The last mapping function studied here is derived from electronic densities estimated by the IRI model (Bilitza, 1990). The IRI model, in development since the 1960s by the Committee on Space Research and the International Union of Radio Science, leverages measurements from ionosondes, scatter radars, topside sounders, and instruments on satellites and rockets to create a global model for the ionosphere over time. IRI electronic densities are dependent on both space and time. The dark blue curve in the top panel of Figure 2 displays a particular vertical density profile estimated above the WAAS ZLA reference station in Palmdale, California, on 2014-02-19 at 00:00.

Using an IRI density function, N_e,IRI, a slant-to-vertical mapping function can be obtained using the following ratio:

$G_{I R I} = \frac{\int_{A}^{B} N_{e, I R I} (x) d x}{\int_{C}^{D} N_{e, I R I} (x) d x}$ 18

where points A and B are along the vertical path and C and D are along the slant path outlined in Figure 1. Following Xiang et al. (2019), integration is carried out numerically while IRI model evaluation is performed using PyIRI, a Python re-build of the IRI model’s Fortran implementation (Forsythe & Burrell, 2023; Forsythe et al., 2024).

2.3 Spatiotemporal Modeling for DCB Estimation

Receiver and satellite DCBs are jointly estimated via spatiotemporal models in the vertical domain that fit piecewise linear spherical harmonic basis functions to scaled differences in GNSS observables. Our approach involves dimensional reduction of the four-dimensional ionosphere along the Earth’s radial direction, resulting in a three-dimensional model along a spherical shell above the Earth’s surface. Each ionospheric delay is attributed to an IPP at an ionospheric height of h_I.

Consider the following spatiotemporal model for VTEC along the spherical shell:

$ν_{model} (ϕ_{mag}, λ_{s}, t; a_{n}^{m}, b_{n}^{m}) = \sum_{n = 0}^{N_{\deg}} \sum_{m = 0}^{n} {\tilde{P}}_{n}^{m} (\sin ϕ_{mag}) [a_{m, n} (t) \cos (m λ_{s}) + b_{m, n} (t) \sin (m λ_{s})]$ 19

where:

$λ_{s} = \mod (λ_{mag} + \frac{π}{43200} (t - t^{*}), 2 π)$ 20

${\tilde{P}}_{n}^{m} (\sin ϕ_{mag}) = \sqrt{\frac{2 n + 1}{4 π} \frac{(n - m)!}{(n + m)!}} P_{n}^{m} (\sin ϕ_{mag})$ 21

Here, (ϕ_mag, λ_mag) represents a coordinate pair of APEX geomagnetic latitude and mean sun-fixed and APEX geomagnetic longitude of the IPP, respectively (Richmond, 1995). N_deg represents the spherical harmonic degree of the spherical harmonic model. $P_{n}^{m}$ and ${\tilde{P}}_{n}^{m}$ are the unnormalized and normalized associated Legendre functions, respectively, whereas a_m,n and b_m,n are piecewise linear spherical harmonic coefficients with a temporal discretization of 2 h. t represents GPS time on some time interval t₀ ≤ t ≤ t_f. The offset t^∗ = 50,400 s is a solar shift to align λ_s = 0 with 14:00 in local time, which is approximately when the ionosphere is most active throughout the day.

By replacing the vertical delay in Equation (9) with the spatiotemporal model for TEC and the pseudorange difference with more accurate leveled carrier-phase differences, $\hat{δ}$ , and rearranging terms, we can write the following:

$ν_{model, k} - ξ {\tilde{G}}_{k} Γ_{r x, k} - ξ {\tilde{G}}_{k} Γ_{s v, k} = ξ {\tilde{G}}_{k} {\tilde{δ}}_{k}$ 22

where $\tilde{G}$ is an estimated slant-to-vertical mapping function approximating the value in Equation (8). Terms with unknown variables (a_m,n, b_m,n, Γ_rx, Γ_sv) are grouped on the left side of Equation (22) whereas scaled measurements are on the right side, forming an organized linear system. Unlike spherical harmonic coefficients, it is assumed that the DCBs Γ_rx and Γ_sv are constant throughout the day.

Owing to the ambiguity in determining satellite and receiver DCBs, a common but important condition requiring all satellite DCBs for each constellation to sum to zero is assumed:

$\sum_{i = 1}^{N_{s v}} {\tilde{Γ}}_{r x}^{(i)} = 0$ 23

where N_sv is the number of satellites for the constellation. This zero-sum satellite DCB condition is important during our analysis because any spatiotemporally constant errors in DCB estimation will manifest in receiver DCBs and not satellite DCBs. As a result, the analysis provided here will focus on receiver DCBs and not satellite DCBs. We note that GIMs also make the zero-sum satellite DCB assumption.

Computational details for spatiotemporal modeling are further outlined in Section 7.1 of the appendix.

2.4 Analytical Approximation for Mean Receiver DCBs

Beyond computational studies, a sequence of assumptions in this section will enable the derivation of a closed-form approximation for the mean receiver DCB for a ground network, allowing for an analytical understanding of how receiver DCB estimates may depend on the slant-to-vertical mapping.

Our first assumption will be the crudest one: that the VTEC is constant across all measurements, meaning that v_model,k = V for all k and for some value V. The ionosphere is widely known to have significant diurnal and geomagnetic variation, but this crude assumption will yield simplicity in the modeling. We define the following:

$f (V, Γ_{r x}) = \frac{1}{2} {\sum_{k = 1}^{M} [\frac{1}{ξ} V - G_{k} ({\hat{δ}}_{k} + Γ_{r x, k} + Γ_{s v, k})]}^{2}$ 24

Let $\tilde{V}$ and ${\tilde{Γ}}_{r x}$ be the minimizers for f, or the ordinary least-squares approximation, as follows:

$\begin{matrix} 0 = {ξ \frac{\partial f}{\partial V} |}_{(\tilde{V}, {\tilde{Γ}}_{r x})} & = \sum_{k = 1}^{M} [\frac{1}{ξ} \tilde{V} - G_{k} ({\hat{δ}}_{k} + {\tilde{Γ}}_{r x, k} + Γ_{s v, k})] \\ = \frac{M}{ξ} \tilde{V} - \sum_{k = 1}^{M} G_{k} ({\hat{δ}}_{k} + {\tilde{Γ}}_{r x, k} + Γ_{s v, k}) \end{matrix}$ 25

$\begin{matrix} 0 = {\frac{\partial f}{\partial Γ_{r x}} |}_{(\tilde{V}, {\tilde{Γ}}_{r x})} & = - \sum_{k = 1}^{M} G_{k} [\frac{1}{ξ} \tilde{V} - G_{k} (δ_{k} + {\tilde{Γ}}_{r x, k} + Γ_{s v, k})] \\ = - \frac{M}{ξ} \tilde{V} \bar{G} + \sum_{k = 1}^{M} G_{k}^{2} ({\hat{δ}}_{k} + {\tilde{Γ}}_{r x, k} + Γ_{s v, k}) \end{matrix}$ 26

where $\bar{G} = \frac{1}{M} Σ_{i = 1}^{M} G_{k}$ denotes the mean slant-to-vertical mapping function across all observations.

Combining Equation (25) multiplied by $\bar{G}$ with Equation (26) yields the following:

$0 = \sum_{k = 1}^{M} G_{k} (G_{k} - \bar{G}) ({\hat{δ}}_{k} + {\tilde{Γ}}_{r x, k} + Γ_{s v, k})$ 27

Next, if we assume that all stations have similar look angles to satellites and similar numbers of measurements, then it is natural to further assume the following approximation:

$\sum_{k = 1}^{M} G_{k} (G_{k} - \bar{G}) {\tilde{Γ}}_{r x, k} \approx \bar{{\tilde{Γ}}_{r x}} \sum_{k = 1}^{M} G_{k} (G_{k} - \bar{G})$ 28

where:

$\bar{{\tilde{Γ}}_{r x}} = \frac{1}{N_{r x}} \sum_{n = 1}^{N_{r x}} {\tilde{Γ}}_{r x}^{(n)}$ 29

Here, ${\tilde{Γ}}_{r x}^{(n)}$ denotes the minimizing n-th receiver bias with N_rx counting the number of receivers. The mean minimizing receiver bias can then be approximated as follows:

$\bar{{\tilde{Γ}}_{r x}} \approx - \frac{\sum_{k = 1}^{M} ({\hat{δ}}_{k} + Γ_{s v, k}) w_{k}}{\sum_{k = 1}^{M} w_{k}}$ 30

where:

$w_{k} = G_{k} (G_{k} - \bar{G})$ 31

We further apply the assumption of a zero-mean satellite DCB condition in Equation (23) by taking the following approximation:

$\sum_{k = 1}^{M} w_{k} Γ_{s v, k} \approx \bar{Γ_{s v}} \sum_{k = 1}^{M} w_{k} \approx 0$ 32

which yields a closed-form approximation for the mean receiver DCB estimated by a least-squares fit to Equation (22) as follows:

$\bar{{\tilde{Γ}}_{r x}} \approx - \frac{\sum_{k = 1}^{M} w_{k} {\hat{δ}}_{k}}{\sum_{k = 1}^{M} w_{k}}$ 33

The simplest interpretation of the approximation in Equation (33) is a weighted average of leveled carrier-phase differences. The largest weights are attributed to high-elevation measurements, and the lowest weights are attributed to either low-elevation measurements or measurements with $G_{k} \approx \bar{G}$ . The form in Equation (33) can be equivalently rewritten in terms of statistical functions as follows:

$\begin{matrix} \bar{{\tilde{Γ}}_{r x}} & \approx \frac{cov ({\hat{δ}}_{k} G_{k}, G_{k})}{var (G_{k})} = corr ({\hat{δ}}_{k} G_{k}, G_{k}) \sqrt{\frac{var ({\hat{δ}}_{k} G_{k})}{var (G_{k})}} \\ = corr ({\hat{δ}}_{k} G_{k}, G_{k}) \sqrt{{\bar{\hat{δ}}}^{2} + var ({\hat{δ}}_{k}) + \frac{{\bar{G}}^{2} var ({\hat{δ}}_{k})}{var (G_{k})}} \end{matrix}$ 34

where cov, var, and corr denote the covariance, variance, and Pearson correlation coefficient functions, respectively, and $\bar{\hat{δ}} = \frac{1}{M} Σ_{i = 1}^{M} δ_{k}$ .

The degree of sensitivity of receiver DCB estimates to the mapping function, G, is revealed by the approximation to the mean receiver DCB in Equation (34). Slant-to-vertical mapping functions with larger means, $\bar{G}$ , and smaller variances, var (G_k), across all collected measurements, such as the JPL extended slab model and the 550-km thin shell model, whose 1/ G representations are shown in Figure 2, are likely to yield larger $\bar{{\tilde{Γ}}_{r x}}$ approximations than those with lower respective means and larger respective variances, such as the 350-km thin shell, because of an inflation of the last term in the radical of Equation (34).

2.5 Decomposition of Receiver DCB Estimate Variation

Whether computed as in this study or estimated by GIMs, receiver DCB estimates will vary over time owing to estimation error or actual physical changes in hardware. Here, a multifactor linear model is proposed to explain receiver DCB variations over time:

$Γ_{r x, i} = α_{0} + α_{T} T_{i} + α_{V} V_{i} + Δ_{hardware, i} + \in_{i}$ 35

When T represents the local temperature, V represents the ionospheric activity as captured by VTEC estimates, Δ_hardware,i is the offset due to known events in which receiver hardware was replaced or upgraded, and ∈ represents noise and unexplained variation. α variables are unknown coefficients to be fit to daily receiver DCB estimates.

Receiver DCB sensitivity to temperature has been well-studied and, seemingly, changes the physical properties of signal reception (Coster et al., 2013; Liu et al., 2019; Zha et al., 2019). As a result, it is natural to include temperature sensitivity in the multifactor model. In contrast, the vertical TEC may not be expected to have a significant effect on the physical aspects of signal reception, if any at all, and has been considered as a controversial factor (Zhong et al., 2016; Choi et al., 2019; Zhong et al., 2020; Choi et al., 2020). Instead, it is possible that any explanatory behavior of VTEC in receiver DCB variation is suggestive of modeling errors in joint estimation of DCBs and VTEC.

3 RESULTS

This section presents computational tests of our modeling framework, first validating the framework in two ways: (1) by accurately estimating known DCBs in a simulated ionosphere and (2) by producing receiver DCB estimates that are consistent with GIMs. Following validation, computational and analytical receiver DCB estimates over a 4.5-year period are studied via an array of mapping functions. Finally, the variations of the receiver DCB estimates produced using each mapping function are explained by the measured temperature and estimated VTEC, revealing potential obliquity errors.

3.1 Validation and Testing with an IRI-Simulated Ionosphere

To validate the computational approach in this study, we tested the DCB estimation methodology presented here on a simulated ionosphere using the IRI model. Synthetic pseudorange and carrier-phase measurements were computed for actual GPS satellite positions with respect to actual WAAS reference station (WRS) positions, denoted by orange markers in Figure 3, by adding the following measurements:

True range from satellite to receiver
Simulated slant delays obtained by numerically integrating IRI electronic densities along slant paths, as described in Section 2.2.3
Satellite and receiver DCBs randomly determined daily

FIGURE 3

Map of the WRS network along with 18 IGS stations with regularly published GIM DCBs

All other errors and noise are ignored in simulations in order to allow for the majority of DCB estimation inaccuracy to stem from dimensional reduction and obliquity errors. Simulations were conducted during the peak of Solar Cycle 24 from 2014-02-01 to 2014-03-31, a timeframe for which the IRI was calibrated with historical measurements. Further details on ionospheric simulations are outlined in Section 7.2 in the appendix.

The known randomly generated DCBs were recaptured by evaluating the spatiotemporal models for each of the 59 days in the testing set individually, each with a different set of DCBs, with different slant-to-vertical mapping functions. In this process, simulated daily DCBs were nearly uncorrelated across consecutive days in the testing set, yielding a greater independence of sampling from day to day versus assuming relatively constant simulated DCBs over time. Eight different mapping functions were used for testing:

Thin shells at 300 km, 350 km, 400 km, 450 km, 500 km, and 550 km
The JPL extended slab mapping function with IPPs determined at 450 km
The IRI-based mapping function with IPPs determined at 350 km

Errors in receiver DCB estimation were computed on a daily basis for all 38 WRSs for all eight mapping functions. Table 1 lists the mean of all such errors for each mapping function in the second column and the mean of each WRS’s root-mean-square error computed over time. Although the IRI-based mapping function was nearly exact and yielded an overall mean error of only 0.089 ns, this function was not the most accurate in recapturing receiver DCBs; rather, the 350-km thin shell function was the most accurate. One potential reason for this finding is that the mathematical approach for DCB estimation assumes that the ionosphere is a thin shell when fitting combinations of piecewise linear spherical harmonic functions; thus, using a mapping function that aligns with this approach to dimensional reduction could provide the best fit. Another reason may be that the assumed ionospheric height of 350 km for IPP determination in the IRI-based mapping function may not be the most appropriate when considering ionospheric dynamics.

View this table:

Table 1 Accuracy in Estimating Randomly Generated DCBs in an IRI-Simulated Ionosphere Using Eight Different Mapping Functions From 2014-02-01 to 2014-03-31

The low mean receiver DCB estimation errors for the 350-km thin shell and the IRI-based mapping function validate the methodology outlined in Section 2.3. However, these errors increase to over 1.5 ns when the 550-km thin shell function or the JPL extended slab function is used, suggesting that there may be a strong sensitivity to the choice of mapping function.

3.2 Estimation of Solar Cycle 25 Receiver DCB

Following our validation with simulated data, actual GNSS measurements fed the spatiotemporal modeling framework outlined in Section 2.3 to estimate receiver DCBs over a sufficiently long timeframe to study their dispersion and variation. Specifically, receiver DCBs were computed for the thin shell, JPL extended slab, and IRI-based mapping functions over a 4.5-year timeframe from January 2020 to June 2024 and compared with receiver DCB estimates computed by GIMs. For all model evaluations, only the mapping function was varied; all other components were left unchanged.

3.2.1 Construction of a WRS Proxy Network

Receiver DCB estimates were computed for a WAAS-centric network of 38 stations to enable a comparison with GIM receiver DCB estimates. To begin network construction, we identified 18 stations from the IGS, shown as blue triangles with 4-letter station codes in Figure 3, that were in close proximity to distinct WRSs and were subject to regular GPS L1(P)–L2(P) DCB estimates by five IGS GIM products (JPL, COD, ESA, CAS, EMR). For a given day, each of these stations was assigned as a “proxy” to their nearby WRS if low-rate observables (30 s/sample) had a sufficient amount of data (an average of at least five usable observables above a 5° elevation mask per epoch) and had data available the previous and following days for proper carrier-phase difference leveling to code difference measurements. One by one, continuously operating reference stations (CORSs) from the IGS, University Navstar Consortium (UNAVCO), and the National Oceanic and Atmospheric Administration’s (NOAA) National Geodetic Survey (NGS) networks were added to complete a 38-station WRS proxy network for a given day; here, the nearest WRS-CORS pairings that were not within a prescribed 200-km radius of any other proxy site were added. It is important to note that the ‘A’ thread receivers in the WRS network provide NGS with observation data; thus, WRS receivers were often used within the proxy networks created. Figure 4 shows which stations comprised each day’s proxy network over the 4.5-year testing period, with blue markers denoting the 18 IGS stations with regular DCBs published, green markers representing WAAS data provided by NGS, and purple markers representing all other CORS sites.

FIGURE 4

Stations used for spatiotemporal modeling and DCB estimation from 2020-01-01 to 2024-06-30 including IGS stations with DCBs published by GIMs (blue), WRSs published by NGS (green), and other CORSs from the IGS, UNAVCO, and NGS networks (purple)

3.2.2 Receiver DCB Estimates

Receiver DCBs for WRS proxy stations were estimated each day for all eight mapping functions listed in Section 3.1 from January 2020 to June 2024 using our spatiotemporal modeling framework. Only GPS DCBs were estimated, with a preference for the L1(P)–L2(P) signal pairing. When unavailable, L1(P) measurements were replaced with L1-C/A measurements with the use of L1-C/A–L1(P) DCBs that were estimated by solving a linear system for L1-pseudorange differences that assumed a zero-sum satellite bias condition similar to Equation (23).

With the knowledge that some GIMs use the JPL extended slab model or approximations to it, we compared our receiver DCBs estimated by a regional network using the JPL extended slab function with those from GIMs that were estimated on a global scale; we observed a close agreement for most stations. Figure 5 displays the difference between our JPL extended slab receiver DCBs and those estimated by the two longest established GIMs, JPL and COD, shown in the top and bottom panels, respectively, over the entire testing timeframe. Except for CRO1 and KOKB, almost all DCB residual differences lie within a steady band near zero, with a width that appears to be primarily no larger than 1 ns, suggesting a strong agreement among DCB estimates. CRO1, the proxy for the WRS in San Juan, Puerto Rico, and KOKB, the proxy for the WRS in Honolulu, Hawaii, are both on the periphery of the WRS network and near the geomagnetic equator, where the ionosphere is most active, making their DCBs the most challenging to estimate. As a result, it is unsurprising to observe discrepancies between the GIM DCBs and those from our modeling using the JPL extended slab function. It is also worth noting that the set of GPS satellites used for computing the JPL GIM was different from that used in our models for 12.2% of the days tested, which led to an average contribution of receiver DCB discrepancies of 0.10 ns on those days; in contrast, the COD GIM used a different set of GPS satellites 1.3% of the time, which contributed an average receiver discrepancy of 0.11 ns. Overall, the consistency of the estimation from our model with the JPL extended slab and the JPL and COD GIMs for 16 other IGS station DCBs further validates our modeling framework beyond the simulated testing in Section 3.1.

FIGURE 5

Differences between receiver DCBs modeled using the JPL extended slab mapping function and GIM DCBs estimated by JPL (top) and COD (bottom) for 18 IGS stations, along with respective means (black) from 2020-01-01 to 2024-06-30

Despite the consistency between the GIM receiver DCBs and our regional modeling using the JPL extended slab function, significant DCB estimate differences were noted within our own modeling across different mapping functions. Figure 6 plots receiver DCBs for 4 of the 18 IGS stations of interest (DUBO, GODE, MDO1, FAIR) over time as estimated by (A) the thin shell mapping function at six different ionospheric heights and (B) five GIMs. Across the eight mapping functions, DCB estimates differed largely by no more than 1 ns during the early portion of Solar Cycle 25 in 2020 and 2021, with the largest difference observed for shell heights of 300 km and 550 km; however, when ionospheric activity increased in 2023 and 2024, these discrepancies increased to 2 or 3 ns on a regular basis. Furthermore, the 550-km thin shell DCB estimates displayed general increases over lengthy periods of time whereas the 300-km thin shell DCB estimates showed decreasing trends. Receiver DCBs for the GIMs, which, except for the EMR GIM, showed good agreement with each other, largely follow the DCBs that we modeled at 550 km, as expected, and similarly trended upwards over time. It is also worth noting that DCB estimates for FAIR (located in Fairbanks, Alaska) displayed significant seasonality, with higher values during the warmer summers and lower values during colder winters, likely due to the effects of temperature on the receiver and antenna hardware (Coster et al., 2013; Liu et al., 2019; Zha et al., 2019). We note that the FAIR station also experienced a hardware change around 2021-04-14, which was accounted for in the figure using the methods outlined in Section 2.5 and is further studied in Section 3.2.4.

FIGURE 6

Receiver DCBs for (clockwise from top left) DUBO, GODE, MDO1, and FAIR estimated by six thin shell mapping functions (left) and five GIMs (right) from 2020-01-01 to 2024-06-30

Aggregated outputs from our spatiotemporal modeling and GIMs for each of the 18 notable IGS receivers over the 4.5-year timeframe are shown in Figure 7. These results were obtained after modeling out the hardware changes listed in Table 2 using the approach described in Section 2.5 and further specified in Section 3.2.4. The top panel of Figure 7 displays the availability of each DCB estimate over the testing period, with GIMs shaded in gray on the left and our modeled results unshaded on the right, showing high overall availability. For the GIMs, the mean VTEC plotted in the second panel represents the mean GIM VTEC interpolated at the specified 450-km height above the receiver over the entire 4.5-year period, whereas our modeled estimates for each of the eight mapping functions present the mean of all daily mean spatiotemporal modeled TEC values for all measurements for each station. VTEC estimates for the GIMs, the modeled JPL extended slab, and the 550 km thin shell were noticeably higher than the modeled VTEC at lower ionospheric heights, such as the 450-km and 350-km thin shell results, which were roughly 1.7 and 3.1 TECU lower on average than those for the JPL extended slab function, respectively. Similar trends were observed for the mean receiver DCBs, which are shown relative to the respective JPL extended slab estimates in the third panel of the figure; receiver DCBs for the GIMs, JPL extended slab, and 550-km thin shell exceeded estimates at lower ionospheric heights, with the JPL extended slab results exceeding the DCBs for the 450-km and 350-km thin shell function by approximately 0.63 and 1.17 ns on average, respectively.

FIGURE 7

Aggregate statistics over time for receiver DCBs estimated from GIMs (left, shaded) and spatiotemporal modeling using estimates with eight different mapping functions for 18 IGS stations from 2020-01-01 to 2024-06-30 (right, unshaded), including the proportion of time during which the receiver DCBs were estimated (top panel), the mean of locally estimated VTEC (second panel), the mean receiver DCB differenced with the corresponding JPL extended slab mean (third panel), and the receiver DCB root-mean-square deviation (RMSD) (bottom panel)

View this table:

Table 2 Effective Dates for Which Hardware Changes for Seven IGS Stations Are Accounted in Long-Term Receiver DCB Analysis

Statistical dispersions for the 18 notable IGS receiver DCBs estimated by the eight different mapping functions and the five GIMs are shown in the bottom panel of Figure 7. Dispersions are calculated as root-mean-square deviations of DCB estimates and can be interpreted as the expected accuracy that can be achieved in the calibration of the receiver DCBs. Among the GIM estimates shown on the left side of the figure, the EMR estimates had the highest dispersions whereas COD and JPL seemed to have the overall lowest values. Of the DCBs we modeled, the 450-km thin shell DCBs had noticeably lower dispersions when compared with the other mapping functions and GIMs, indicating that these estimates were the most consistent over time. It is worth noting that receiver DCB estimates for CRO1 and KOKB, the two stations that displayed the highest VTEC values in the second panel of Figure 7 and the least consistency in Figure 5, had the highest dispersions amongst all stations.

Although its estimates showed the highest consistency over time, the 450-km thin shell DCB did not always have the best goodness-of-fit statistics among the eight functions tested. Figure 8 shows rolling 30-day means for the coefficient of determination, R², and reduced χ² statistics, as described in Section 7.3 in the appendix, for the spatiotemporal modeling involved in DCB estimation in the top left and top right panels, respectively, with their ranks across all mapping functions displayed underneath. Different mapping functions achieved the best goodness-of-fit statistics under different ionospheric conditions and during different stretches of time, similar to observations from Wang et al. (2016). The 450-km thin shell function did yield the best overall R² values, but it appears as though a height of approximately 350 km yields the best reduced χ² values. The primary explanation behind these discrepancies is that the reduced χ² does not normalize by residual magnitude whereas R² does, and the lowest and highest modeled TEC values were observed for shell heights of 300–350 km and 500–550 km, respectively.

FIGURE 8

Goodness-of-fit statistics for spatiotemporal modeling using eight different mapping functions from 2020-01-01 to 2024-06-30

Top panels show the 30-day rolling means for the coefficient of determination R² (left) and the reduced χ² statistic (right); the corresponding plots underneath rank the mapping functions by magnitude of the respective statistics.

3.2.3 Comparison of Analytical Approximate Mean Receiver DCBs With Spatiotemporally Modeled Estimates

In contrast to the complexity of our spatiotemporal modeling using piecewise linear spherical harmonic functions, analytical estimation of mean receiver DCBs equivalently expressed in Equations (33) and (34) requires a mean of leveled carrier-phase differences, using mapping-function-dependent weights or a simple covariance/variance ratio. The mean of the DCBs of the 18 IGS stations of interest, as estimated by spatiotemporal modeling for the eight different mapping functions studied, shown in the top left panel of Figure 9, is not well approximated by its analytically estimated counterpart from Equation (33), shown in the top right panel, with differences of over 2 ns. However, when comparing mean DCB discrepancies across mapping functions, the analytical estimates provide good insight into the spatiotemporally modeled results. Notably, the difference of each mean spatiotemporally modeled DCB estimate with respect to the JPL extended slab estimate, plotted in the bottom left panel of Figure 9, looks similar to the differences for the analytically estimated mean DCBs, shown in the bottom right panel. Thus, the estimate in Equation (33) appears to be poor in an absolute sense but seemingly provides a simple insight into the sensitivity of DCBs to the mapping function used for estimation.

FIGURE 9

Spatiotemporally modeled (top left) and analytically estimated (top right) mean receiver biases for 18 IGS stations with their respective differences from the JPL extended slab estimates (bottom) from 2020-01-01 to 2024-06-30

3.2.4 Fits to Receiver DCB Estimates

With significant differences in the variation of DCB estimates based on the mapping function, as shown in Figures 6, 7, and 9, our attention turns to modeling such variation to better understand its nature. The multifactor model described in Equation (35) was fit to receiver DCBs estimated by our spatiotemporal models with eight different mapping functions and for the five GIMs on a daily basis using least-squares regression. Local temperature measurements, T, were daily averages sourced from nearby weather station data provided by NOAA’s National Centers for Environmental Information (NCEI); no reference weather station was further than 100 km from the modeled GNSS receiver. For spatiotemporal modeling, daily VTEC values, V, were determined by the mean of all spatiotemporal model values for all measurements for each station and each mapping function. VTEC values for GIMs were determined as the average of GIM values interpolated at the shell position above the receiver over each entire day. Seven hardware changes modeled by Δ_hardware are listed in Table 2, with no other potential changes accounted for. We note that none of the 18 IGS stations of interest specified any type of temperature stabilization in their most recent IGS receiver log files.

Figure 10 illustrates the components and predictions for multifactor modeling of receiver DCB estimates using our JPL extended slab (blue), 450-km thin shell (red), and 350-km thin shell (green) mapping functions for (clockwise from top left) DUBO, GODE, MDO1, and FAIR. For each of the four station’s subfigures, the top panel shows the spatiotemporally modeled receiver DCB in dark coloring, with the predicted DCB from fits using Equation (35) shown in a fainter color. The middle and bottom panels show the timeseries used for V and T, respectively. First, the spatiotemporally modeled receiver DCBs appear to be well shadowed by the fitted DCBs, with notable fluctuations captured by the model. Second, the DUBO, GODE, and MDO1 DCBs estimated using the JPL extended slab function appear to be highly correlated with the daily VTEC beneath them, whereas the 350-km thin shell DCBs appear to be strongly negatively correlated with the daily VTEC. For FAIR, all DCB estimates appear to display a seasonality that is captured by not only temperature but also by daily VTEC values.

FIGURE 10

Spatiotemporally modeled DCBs (dark markers) plotted above fits (light markers) based on daily local VTEC and temperature estimates (below) for (clockwise from top left) DUBO, GODE, MDO1, and FAIR stations from 2020-01-01 to 2024-06-30

Aggregate statistics from the multifactor fitting for each GIM (left side of each subfigure, shaded gray) and our receiver DCB estimation using eight different mapping functions (right side of each subfigure, unshaded) for each of the 18 IGS stations of interest are shown in Figure 11. Outcomes for fits with less than 50% data availability, as shown in the top panel of Figure 7, are excluded from the plots. The top left panel of Figure 11 displays the Pearson correlation coefficient between Γ_rx and V for the fitting whereas the top right panel shows the regression coefficient α_V. GIM DCB–VTEC correlations range primarily between 0.5 and 0.8 for almost all stations, indicating a high correlation between two quantities whose relationship has been debated: receiver hardware and space weather (Zhong et al., 2016; Choi et al., 2019; Zhong et al., 2020; Choi et al., 2020). Corresponding α_V values for GIMs were often between 0.03 and 0.06 ns/TECU, suggesting that a VTEC of only 10 TECU could account for 0.3–0.6 ns of DCB variation. Modeled 550-km thin shell and 350-km thin shell correlations and α_V values appear to have similar magnitudes, although the values for the 350-km thin shell function were negatively oriented. As thin shell heights dropped from 550 km to 350 km, DCB–VTEC correlations and α_V values reached zero at approximately 450 km, suggesting that ionospheric activity had little effect on DCB estimation.

FIGURE 11

Aggregate correlations (left top and middle) and regression coefficients (right top and middle) for daily local VTEC (top row) and temperature (middle row) estimates along with the coefficient of determination R² (bottom left) and the reduced χ² statistic (bottom right) for regressions using estimates for 18 IGS stations from 2020-01-01 to 2024-06-30

The middle left panel of Figure 11 displays the Pearson correlation coefficient between Γ_rx and T for the fitting whereas the middle right panel shows the regression coefficient α_T. Stations at high latitude, such as FAIR and IQAL, and at high altitude, such as QUIN and PRDS, had the largest magnitudes for DCB–temperature correlations and α_T values, suggesting the potential influence of temperature on actual receiver hardware. Notably, these correlations and α_T values were much less sensitive to changes in mapping function than the DCB–VTEC correlations and α_V values. It is worth noting that CRO1 and KOKB exhibited some of the highest α_T magnitudes and were most sensitive to the mapping functions used; however, goodness-of-fit statistics for the multifactor fitting highlighted relatively poor overall fits for both of those stations, as shown by the coefficient of determination, R², and reduced χ² statistics shown in the bottom panels of Figure 11, which displayed very high dispersion in the estimates in Figure 7. Further examination of the R² plot suggests relatively poor fits for several receiver DCBs estimated by the 450-km thin shell mapping function; however, the R² statistic normalizes by dispersion, unlike the reduced χ² statistic, and the 450-km thin shell results generally displayed the lowest dispersions in Figure 7.

4 DISCUSSION AND CONCLUSIONS

Analyses of receiver DCB estimates from regional spatiotemporal modeling and GIM receiver DCBs over a 4.5-year period revealed a sensitivity of DCBs to slant-to-vertical mapping functions. DCBs estimated via the JPL extended slab function, which is effectively used by GIMs, and the 350-km thin shell function, which is outlined in the Radio Technical Commission for Aeronautics (2020), showed significant variation over time with respective root-mean-square deviations of 0.49 ns and 0.47 ns, which were notably higher than the 0.34-ns variation observed for the 450-km thin shell function. Some of the variation from the JPL extended slab and 350-km thin shell receiver DCBs can likely be attributed to modeling errors because the variation aligned with ionospheric activity, as evidenced by respective mean correlation coefficients of 0.71 and -0.63 across the 18 North American IGS stations tested. In contrast, receiver DCBs estimated by the 450-km thin shell function were nearly uncorrelated with ionospheric activity, with a mean correlation coefficient of -0.01 across the same set of stations, suggesting that the 450-km thin shell function was a more appropriate choice than other mapping functions for North America during the 2020-2024 timeframe studied. Accordingly, our analysis suggests that GIM receiver DCBs, which are effectively estimated by the JPL extended slab mapping function, carry obliquity errors that cause overestimation, consistent with findings from Reid et al. (2024). Furthermore, the observed combination of obliquity errors and modeled TEC suggest that GIM TEC values are likely to be overestimated as well.

This study focused on the variation of receiver DCB estimates and not satellite DCB estimates; however, a primary factor in this study is the zero-sum satellite DCB constraint in Equation (23). In the estimation framework outlined in Section 2.3, any modeling errors that are spatiotemporally constant will be attributed to all receiver DCB estimates, as attributing the errors to all satellite DCB estimates would result in a non-zero sum of satellite DCBs. As these modeling errors vary with ionospheric activity, receiver DCBs exhibit more variation due to modeling errors than their satellite counterparts. If the zero-sum satellite DCB constraint in Equation (23) was modified so that all receiver DCBs were required to sum to zero, these spatiotemporally constant modeling errors would be captured by satellite DCBs instead of receiver DCBs. It is important to note that the actual, not estimated, DCBs of satellites are likely more stable than actual receiver DCBs owing to the onboard oscillators, which are more stable than those used with ground-based receivers, resulting in less actual DCB variation.

From an aviation-focused perspective, the choice of mapping function for DCB estimation is crucial for computing the absolute ionospheric measurements, or “truth,” that seed the GIVE undersampled threat model. Any DCB errors that are removed from geometry-free combinations will manifest as errors in absolute ionospheric measurements in the vertical domain that are directly proportional to an elevation-dependent slant-to-vertical mapping function. For threat estimation occurring in the vertical domain (Altshuler et al., 2001; Sparks et al., 2001; Blanch et al., 2002), GIVE undersampled threat models and overall L1 SBAS performance can, ultimately, be sensitive to the mapping functions used for DCB estimation.

Future studies could perform similar analyses of mapping function sensitivity by using global models, incorporating multiple constellations, and accounting for the diurnal variation of ionospheric heights (Roble, 1975). Furthermore, the modeling framework presented here could be extended to account for multiple shells, similar to the work by Maruyama et al. (2021); however, our results using receiver DCBs computed from the JPL GIM that uses three shells (Komjathy et al., 2002) and findings from Maruyama et al. (2021) suggest that the shell heights and mapping functions used with multi-shell models still have a significant effect on DCB estimation. To better understand DCB accuracy, spatiotemporal models computed or planar fits to absolute truth computed via different mapping functions can be used to account for ionospheric delays in a single-frequency precise point positioning (PPP) model. In this case, lower PPP errors would be suggestive of more accurate mapping functions for DCB estimation.

HOW TO CITE THIS ARTICLE

Ovadia, J., Champlin, D., Pong, C., & Altshuler, E. (2025). Sensitivity of receiver differential code bias estimates to mapping functions. NAVIGATION, 72(4). https://doi.org/10.33012/navi.723

5 ACKNOWLEDGMENTS

The contents of this paper are the views of the authors and do not represent those of the Federal Aviation Administration.

The authors acknowledge helpful conversations with Luke Bucklew and Julia Garver from Sequoia Research Corporation and Larry Sparks, Attila Komjathy, Anthony Mannucci, Siddarth Krishnamoorthy, Panagiotis Vergados, and Olga Verkhoglyadova from JPL.

The authors appreciate the data provided by IGS, UNAVCO, Natural Resources Canada (NRCan), and NOAA’s NGS and NCEI and the modeling provided by IRI and PyIRI.

Appendix

7 APPENDIX

GPS satellite positions were determined using precise orbit data from COD (Dach et al., 2024), propagated using Everett’s formula implemented with a 12-point interpolation stencil (Nykiel, 2015).

IPPs in geodetic coordinates $(ϕ_{geo}, λ_{geo})$ from a user’s geodetic position $(ϕ_{u}, λ_{u})$ are determined in radians using the azimuth, α, and elevation, E, look angles as follows:

$ϕ_{geo} = {sin}^{- 1} (sin ϕ_{u} cos ψ_{p p} + cos ϕ_{u} sin ψ_{p p} cos α)$ 36

$λ_{geo} = λ_{u} + {sin}^{- 1} (\frac{sin ψ_{p p} sin α}{cos ϕ_{geo}})$ 37

where:

$ψ_{p p} = \frac{π}{2} - E - \sin^{- 1} (\frac{R_{e}}{R_{e} + h_{I}} \cos E)$ 38

We note that the expression for geodetic longitude, λ_geo, in Equation (37) is used for all latitudes, unlike the Radio Technical Commission for Aeronautics (2020) approach, which uses the geodetic longitude only for the latitude range of $- 70^{\underline{o}} < ϕ_{geo} < 70^{\underline{o}}$ .

7.1 Computational Details for Spatiotemporal Modeling

The APEX geomagnetic coordinates in Section 2.3 are determined from geodetic coordinates using the Python package “apexpy,” version 2.0.1 (Burrell, 2023).

7.1.1 Cycle Slip Detection

We denote carrier measurements for signals A and B using Φ_A and Φ_B, respectively, and denote code measurements for signals A and B as P_A and P_B, respectively.

The algorithm applied for cycle slip detection uses a difference between the Melbourne-Wubbena combination and its rolling mean as a signal for a threshold-based approach to cycle slip detection:

The Melbourne-Wubbena combination is evaluated as follows:
$C_{M W}^{(i)} = λ_{W L} (Φ_{B}^{(i)} - Φ_{A}^{(i)}) - \frac{1}{f_{A} + f_{B}} (f_{B} P_{B}^{(i)} + f_{A} P_{A}^{(i)})$ 39
where $λ_{W L} = c / (f_{B} - f_{A})$ is the wide-lane λ and P_s is the code measurement for signal s.
A rolling mean is computed over the N_MW-sized window of $[t^{(i)}, t^{(i + N_{M W}^{- 1})}]$ for all possible i:
$C_{r o l l}^{(i)} = \frac{1}{N_{M W}} \sum_{k = 0}^{^{N} M W - 1} C {}_{M W}^{(i - k)}$ 40
A dynamic estimate of the magnitude of the short-term deviation of the Melbourne-Wubbena combination is computed as follows:
$C_{n o i s e}^{(i)} = C_{M W}^{(i)} - C_{r o l l}^{(i)}$ 41

A set of outliers amongst ${C_{n o i s e}^{(i)}}$ measurements is iteratively determined using the following algorithmic pseudocode:
- Initialize an empty outlier set.
- While (new outliers are found):
  - - Identify new outliers as any C_noise values such that:
    $\frac{| C_{n o i s e} - mean (C_{n o i s e}) |}{std (C_{n o i s e})} > Z_{n o i s e}$ 42
  - where mean(⋅) and std(⋅) denote the mean and standard deviation, respectively, of the non-outlier set of values.
- - If new outliers are found, add them to the outlier set.
- - If no new outliers are found, then exit the while loop.
The outlier-removed standard deviation of the Melbourne-Wubbena noise is computed as follows:
$σ_{M W} = std (C_{n o i s e}^{i \notin {outliers}})$ 43
The Melbourne-Wubbena signal threshold is set as a thresholded multiple of the outlier-removed standard deviation:
$C_{t h r e s h} = \min (C_{h i g h t h r e s h}, \max (C_{l o w t h r e s h}, 2 Z_{n o i s e} σ_{M W}))$ 44
The magnitude of differences of rolling measurements, scaled by the widelane λ, is taken as a signal for carrier-phase cycle slip and run-off detection:
$C_{diff}^{(i)} = \frac{1}{λ_{W L}} | C_{r o l l}^{(i + N_{M W})} - C_{r o l l}^{(i)} |$ 45
We note that there are quantitative differences in the average Melbourne-Wubbena combination based on N_MW measurements up to time t⁽ⁱ⁾ and the average based on N_MW measurements after time t⁽ⁱ⁾.
The Melbourne-Wubbena-detected anomalies are characterized as any k such that $C_{d i f f}^{(k)} > C_{t h r e s h}$ .

Going forward, a single satellite’s observations from a single receiver is re-partitioned into intervals based on discontinuities determined by either the presence of a δ_{discontinuity} -second data anomaly or a Melbourne-Wubbena-detected anomaly. Then, each of these intervals is used to determine carrier-to-code leveling constants, as specified further in Section 7.1.2.

Our parameter choices for cycle slip detection are as follows:

$\begin{matrix} N_{M W} = 60 s (2 epochs) \\ Z_{n o i s e} = 3 \\ C_{h i g h t h r e s h} = 1.25 \\ C_{l o w t h r e s h} = 5 \end{matrix}$

7.1.2 Carrier-to-Code Leveling

Let us continue to denote carrier measurements for signals A and B using Φ_A and Φ_B, respectively, and denote code measurements for signals A and B as P_A and P_B, respectively. These measurements are input from receiver independent exchange format (RINEX) observation files; when either code or carrier measurements are not available for a single observable, the measurement is considered unusable and removed from the set of observables. Furthermore, data records in which the magnitude of the code pseudorange difference exceeds a prescribed threshold of D_outlier are considered statistical outliers and are removed from the set of observables. Accordingly, outliers are defined by |P_B – P_A| > D_outlier.

The time series that can be used to estimate a carrier-to-code leveling constant, denoted by ζ, is given by the following:

$ζ_{A, B} = (P_{B} - P_{A}) + c (\frac{Φ_{B}}{f_{B}} - \frac{Φ_{A}}{f_{A}})$ 46

where c = 299,792,458 m/s and f is the corresponding signal’s frequency in Hz.

The carrier-phase difference is then leveled to the code pseudorange difference via the following dispersion-based weighting scheme:

The observations of a single satellite are parsed from a single receiver into intervals in which no gaps in collection exceed a time interval length of δ_{discontinuity} and no cycle slips are detected.
For each interval, which we assume spans N_T seconds:
- If N_T < max (T_thresh, T_roll), then omit the measurements in the interval.
- For measurements between T_roll seconds and N_T - T_roll seconds, perform an Ornstein-Uhlenbeck model fit over a T_roll-second rolling window to obtain a time series of noise measurements, denoted by $σ_{O U}^{(i)}$ , with the superscript denoting the time index i.
- Extend the $σ_{O U}^{(i)}$ time series over the entire interval by assigning its value at T_roll seconds to all prior T_roll-1 epochs and assigning its value at N_T - T_roll to the last T_roll-1 epoch.
- º Compute the time-dependent weight vector over the entire interval of N_T epochs:
  $w_{O U}^{(i)} = \frac{1}{β_{O U, E C 50} + σ_{O U}^{(i)}}$ 47
  where β_OU,EC50 is the half-maximal effective concentration parameter for the Hill function for w_OU.
- Iteratively determine a set of outliers amongst ${ζ_{A, B}^{(i)}}$ measurements using the following algorithmic pseudocode:
  - - Initialize an empty outlier set.
  - - While (new outliers are found):
    - - Identify new outliers as any ζ_A,B values such that:
      $\frac{| ζ_{A, B} - mean (ζ_{A, B}) |}{std (ζ_{A, B})} > Z_{o u t l i e r}$ 48
      where mean(⋅) and std(⋅) denote the mean and standard deviation, respectively, of the non-outlier set of values.
    - - If new outliers are found, add them to the outlier set.
    - - If no new outliers are found, then exit the while loop.
Normalize the weight vector:
${\tilde{w}}_{O U}^{(i)} = \frac{w_{O U}^{(i)}}{\sum_{k \neq o u t l i e r} w_{O U}^{(k)}}$ 49
where the superscript denotes the time index corresponding to non-outlier measurements.
Determine the leveling constant, C_level, of the interval by the weighted mean of non-outlier ζ_A,B measurements using the weight vector ${\tilde{w}}_{O U}$ :
$C_{level} = \sum_{k \neq outlier} w_{O U}^{(k)} ζ_{A, B}^{(k)}$ 50
Compute the leveled carrier-phase difference, $\hat{δ}$ , on the interval as follows:

$\hat{δ} = - c (\frac{Φ_{B}}{f_{B}} - \frac{Φ_{A}}{f_{A}}) + C_{level}$ 51

Our parameter choices for carrier-to-code leveling are as follows:

$\begin{matrix} δ_{d i s c o n t i n u i t y} = 5 s \\ D_{o u t l i e r} = 300 m \\ T_{r o l l} = 180 s \\ T_{t h r e s h} = 300 s \\ β_{O U, E C 50} = 0.01 m \\ Z_{o u t l i e r} = 3 \end{matrix}$

7.1.3 QR-Decomposition Recursive Least-Squares (QR-RLS) Solver

Consider the overdetermined linear system in Equation (22), rewritten as follows:

$A_{M \times N} x_{N \times 1} = b_{M \times 1}$ 52

where A is the coefficient matrix from the left-hand size, x is a vector of unknowns, b represents the right-hand side, and the subscript denotes matrix size.

Suppose that the interval of integers [1, M] is partitioned by ${n_{k}}_{k = 0}^{k}$ , with n₀ = 1 and n_K = M for some 1 ≤ K ≤ M - 1. Let A_(i,j) for some i ≤ j denote the block matrix of size (j-i+1)×N formed by rows i through j of A and b_(ij) denote the vector of length j–i+1 formed by entries i through j of b.

The QR-RLS algorithm can be used to solve for x_N×1 via an approach that iteratively processes partitions of A:

Initialize R_0×N = 0_0×N as an empty matrix with no rows.
Initialize B_0×1 = 0_0×1 as an empty vector with no rows.
For k = 1 to K:
- Construct the block matrix:
  $R_{n_{k} \times N} = [\begin{matrix} R_{n_{k - 1} \times N} \\ A_{(n_{k - 1}, n_{k})} \end{matrix}]$ 53
- Q, $R_{\min} (N, n_{k}) \times N \leftarrow$ QR_Decomposition $(R_{_{n_{k} \times N}})$
- Apply matrix multiplication:
  $B_{\min (N, n_{k}) \times N} = Q^{T} [\begin{matrix} B_{n_{k - 1} \times 1} \\ b_{(n_{k - 1}, n_{k})} \end{matrix}]$ 54
Solve the upper triangular system using back substitution:

$R_{N \times N} x_{N \times 1} = B_{N \times 1}$ 55

Here, the routine QR_Decomposition(X) computes the matrix decomposition X=QR, where R is an upper triangular and Q is orthogonal. Various algorithms, such as those in which the lower triangular matrix elements are annihilated by a sequence of Givens rotations, can be used for this decomposition.

One caveat to the outlined QR-RLS algorithm is that the order of processing of the first N rows can be important for A matrices with zero entries, specifically sparse A matrices. For the m-th row for some m < N, it is necessary that A_m,m be non-zero.

In practice, the matrix A is partitioned by using chunks of 100,000 rows at a time and modified based on memory constraints.

7.2 IRI Simulation Details

For all IRI electronic densities, the “fluxadjflux” F10.7 cm data set provided by NRCan is used as a solar flux input for PyIRI, and the Consultative Committee on International Radio model (Bilitza et al., 2022) is used to determine the F2 critical frequency for the electronic density at each epoch.

For the IRI-based mapping function in Section 2.2.3, numerical integration is performed for ionospheric heights of 50-1000 km using trapezoidal approximations with a step size of 5 km.

For the simulated slant delays in Section 3.1, numerical integration is performed for ionospheric heights of 25-750 km using trapezoidal approximations with a step size of 25 km. Randomly generated L1-L2 DCBs are determined from the difference of L1 and L2 code biases sampled from zero-mean Gaussian distributions with a standard deviation of 1 m.

7.3 Goodness-of-Fit Statistics

Given a data sequence of scalar values ${y_{i}}_{i = 1}^{M}$ fit by ${f_{i}}_{i = 1}^{M}$ , a model with N_free free parameters, the coefficient of determination R² and the reduced χ² statistic are given as follows:

$R^{2} = \frac{\sum_{i = 1}^{M} {(y_{i} - f_{i})}^{2}}{\sum_{i = 1}^{M} {(y_{i} - \bar{y})}^{2}}$ 56

$χ^{2} = \frac{\sum_{i = 1}^{M} {(y_{i} - f_{i})}^{2}}{M - N_{f r e e}}$ 57

where $\bar{y} = Σ_{j = 1}^{M} y_{j}$ is the mean of {y_i}. We note that M – N_free quantifies the degrees of freedom of the fit, and the numerators of Equations (56) and (57) represent the residual sum of squares of the fit. Higher R₂ values and lower reduced χ² statistics indicate a better goodness of fit.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

6 REFERENCES

↵
Altshuler, E., Fries, R., & Sparks, L. (2001). The WAAS ionospheric spatial threat model. Proc. of the 14th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 2001), Salt Lake City, UT2463–2467. https://www.ion.org/publications/abstract.cfm?articleID=1921
Google Scholar
↵
Bilitza, D. (1990). International Reference Ionosphere 1990. Report 90-22. National Space Science Data Center. https://ntrs.nasa.gov/citations/19910021307
Google Scholar
↵
Bilitza, D., Pezzopane, M., Truhlik, V., Altadill, D., Reinisch, B. W., & Pignalberi, A. (2022). The International Reference Ionosphere model: A review and description of an ionospheric benchmark. Reviews of Geophysics, 60, e2022RG000792. https://doi.org/10.1029/2022RG000792
Google Scholar
↵
Blanch, J., Walter T., & Enge P. (2002). Ionospheric threat model methodology for WAAS. NAVIGATION, 49(2), 103–108. https://doi.Org/10.1002/j.2161-4296.2002.tb00259.x
Google Scholar
↵
Burrell, A. (2023). aburrell/PyIRI: v2.0.1. GitHub.
Google Scholar
↵
Choi, B., Chung, J., & Cho, J. (2010). Receiver DCB estimation and analysis by types of GPS receiver. Journal of Astronomy and Space Sciences, 27(2), 123–128. https://doi.org/10.5140/JASS.2010.27.2.123
Google Scholar
↵
Choi, B., Park, J., Min Rho, K., & Lee, S. (2013). Comparison of GPS receiver DCB estimation methods using a GPS network. Earth, Planet and Space, 65, 707–711. https://doi.org/10.5047/eps.2012.10.003
Google Scholar
↵
Choi, B., Sohn, D., & Lee, S. (2019). Correlation between ionospheric TEC and the DCB stability of GNSS receivers from 2014 to 2016. Remote Sensing, 11(22), 2657. https://doi.org/10.3390/rs11222657
Google Scholar
↵
Choi, B., Sohn, D., & Lee, S. (2020). Reply to comment on Choi et al. correlation between ionospheric TEC and the DCB stability of GNSS receivers from 2014 to 2016. Remote Sens. 2019, 11, 2657. Remote Sensing, 12(21), 3510. https://doi.org/10.3390/rs12213510
Google Scholar
↵
Coco, D., Coker, C., Dahlke S., & Clynch, J. (1991). Variability of GPS satellite differential group delay biases. IEEE Transactions on Aerospace and Electronic Systems, 27(6), 931–938. https://doi.org/10.1109/7.104264
Google Scholar
↵
Coster, A., Gaposchkin, E., & Thorton, L. (1992). Real-time ionospheric monitoring system using GPS. NAVIGATION, 39(2), 191–204. https://doi.org/10.1002/j.2161-4296.1992.tb01874.x
Google Scholar
↵
Coster, A., Williams, J., Weatherwax, A., Rideout, W., & Herne, D. (2013). Accuracy of GPS total electron content: GPS receiver bias temperature dependence. Radio Science, 48, 190–196. https://doi.org/10.1002/rds.20011
Google Scholar
↵
Dach, R., Schaer, S., Arnold, D., Brockmann, E., Kalarus, M., Lasser, M., Stebler, P., & Jäggi, A. (2024). CODE final product series for the IGS. Astronomical Institute, University of Bern. https://doi.org/10.48350/197025
Google Scholar
↵
Federal Aviation Administration (2008). Global Positioning System Wide Area Augmentation System (WAAS) performance standard (1st ed.). https://www.gps.gov/technical/ps/2008-WAAS-performance-standard.pdf
Google Scholar
↵
Forsythe, V., Bilitza, D., Burrell, A., Dymond, K., Fritz, B. A., & McDonald, S. (2024). PyIRI: Whole-globe approach to the International Reference Ionosphere modeling implemented in Python. Space Weather, 22(4), e2023SW003739. https://doi.org/10.1029/2023SW003739
Google Scholar
↵
Forsythe, V., & Burrell, A. (2023). victoriyaforsythe/PyIRI: v0.0.1. GitHub. https://doi.org/10.5281/zenodo.8235173
Google Scholar
↵
Hernández-Pajares, M., Juan, J., Sanz, J., Orus, R., Garcia-Rigo, A., Feltens, J., Komjathy, A., Schaer, S.C., & Krankowski, A. (2009). The IGS VTEC maps: A reliable source of ionospheric information since 1998. Journal of Geodesy, 83, 263–275. https://doi.org/10.1007/s00190-008-0266-1
Google Scholar
↵
Johnston, G., Riddell, A., &Hausler, G. (2017). The International GNSS Service. In Peter J. G.Teunissen & O. Montenbruck (Eds.), Springer handbook of global navigation satellite systems (1st ed., pp. 967–982). Springer International Publishing. https://doi.org/10.1007/978-3-319-42928-1
Google Scholar
↵
Komjathy, A., Wilson, B., Runge, T., Boulat, B., Mannucci, A., Sparks, L., & Reyes, M. (2002). A new ionospheric model for wide area differential GPS: The multiple shell approach. Proc. of the 2002 National Technical Meeting of the Institute of Navigation, San Diego, CA, 460–466. https://www.ion.org/publications/abstract.cfm?articleID=238
Google Scholar
↵
Li, M., Yuan, Y., Wang, N., Liu, T., & Chen, Y. (2018). Estimation and analysis of the short-term variations of multi-GNSS receiver differential code biases using global ionosphere maps. Journal of Geodesy, 92, 889–903. https://doi.org/10.1007/s00190-017-1101-3
Google Scholar
↵
Li, W., Wang, K., & Yuan, K. (2023). Performance and consistency of final global ionospheric maps from different IGS analysis centers. Remote Sensing, 15(4), 1010. https://doi.org/10.3390/rs15041010
Google Scholar
↵
Liu, A., Li, Z., Wang, N., Yuan, C., & Yuan, H. (2019). Analysis of the short-term temporal variation of differential code bias in GNSS receiver. Measurement, 153, 107448. https://doi.org/10.1016/j.measurement.2019.107448
Google Scholar
↵
Mannucci, A., Wilson, B., Yuan, D., Ho, C., Lindqwister, U., & Runge, T. (1998). A global mapping technique for GPS-derived ionospheric total electron content measurements. Radio Science, 33(3), 565–582. https://doi.org/10.1029/97RS02707
Google Scholar CrossRef Web of Science
↵
Maruyama, T., Hozumi, K., Ma, G., Supnithi, P., Tongkasem, N., & Wan, Q. (2021). Double-thin-shell approach to deriving total electron content from GNSS signals and implications for ionospheric dynamics near the magnetic equator. Earth, Planets and Space, 73, 109. https://doi.org/10.1186/s40623-021-01427-y
Google Scholar
↵
Milliken, R., & Zoller, C. (1978). Principle of operation of NAVSTAR and system characteristics. NAVIGATION, 25(2), 95–106. https://doi.org/10.1002/j.2161-4296.1978.tb01322.x
Google Scholar
↵
Nykiel, G., Figurski, M., & Baldysz, Z. (2015). Comparison of GPS precise ephemerides interpolation methods. In Conference: 15th International Multidisciplinary Scientific GeoConference SGEM 2015 (Vol. 2, pp. 161–170). https://doi.org/10.1016/j.measurement.2019.10744810.5593/SGEM2015/B22/S9.020
Google Scholar
↵
Radio Technical Commission for Aeronautics. (2020). Minimum operational performance standards for Global Positioning System/satellite-based augmentation system airborne equipment, technical document RTCA DO-229F. Special Committee SC-159. https://my.rtca.org/nc__store?search=229f
Google Scholar
↵
Reid, B., Themens, D. R., McCaffrey, A., Jayachandran, P. T., Hoque, M., & Mazzella, A. J., Jr. (2024). GNSS differential code bias determination using Rao-Blackwellized particle filtering. Space Weather, 22(5), e2023SW003611. https://doi.org/10.1029/2023SW003611
Google Scholar
↵
Richmond, A. (1995). Ionospheric electrodynamics using magnetic apex coordinates. Journal of Geomagnetism and Geoelectricity, 47(2), 191–212, https://doi.org/10.5636/jgg.47.191
Google Scholar
↵
Roble, R. (1975). The calculated and observed diurnal variation of the ionosphere over Millstone Hill on 23-24 Mar 1970. Planetary and Space Science, 23(7), 1017–1033. https://doi.org/10.1016/0032-0633(75)90192-0
Google Scholar
↵
Sardón, E., Rius, A., & Zarraoa, N. (1994). Estimation of the transmitter and receiver differential biases and the ionospheric total electron content from Global Positioning System observations. Radio Science, 29(3), 577–586. https://doi.org/10.1029/94RS00449
Google Scholar CrossRef Web of Science
↵
Schaer, S. (1999). Mapping and predicting the Earth’s ionosphere using the Global Positioning System. [Thesis, University of Berne]. https://www.scirp.org/reference/referencespapers?referenceid=2405487
Google Scholar
↵
Sparks, L. (2013). Eliminating obliquity error from the estimation of ionospheric delay in a satellite-based augmentation system. Proc. of the ION 2013 Pacific PNT Meeting, Honolulu, HI, 307–318. https://www.ion.org/publications/abstract.cfm?articleID=10983
Google Scholar
↵
Sparks, L., Iijima, B., Mannucci, A., Pi, X., & Wilson, B. (2000). A new model for retrieving slant TEC corrections for wide area differential GPS. Proc. of the 2000 National Technical Meeting of the Institute of Navigation, Anaheim, CA, 464–473. https://www.ion.org/publications/abstract.cfm?articleID=56
Google Scholar
↵
Sparks, L., Mannucci, A.J., Walter, T., Blanch, J., Hansen, A., Enge, P., Altshuler, E., & Fries, R. (2001). The WAAS ionospheric threat model. Proc. of the International Beacon Satellite Symposium 2001, Boston, MA. https://dataverse.jpl.nasa.gov/dataset.xhtml?persistentId=hdl:2014/12868
Google Scholar
↵
Wageeh, A., Doma, M., Sedeek, A., & Elghazouly, A. (2023). DCB estimation and analysis using the single receiver GPS/GLONASS observations under various seasons and geomagnetic activities. Advances in Space Research, 72(9), 3933–3945. https://doi.org/10.1016/j.asr.2023.07.063
Google Scholar
↵
Walter, T., Blanch, J., DeGroot, L., Norman, L., & Joerger, M. (2018). Ionospheric rates of change. In Proc. of the 31st International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2018), Miami, FL, 4158–4170. https://doi.org/10.33012/2018.16112
Google Scholar
↵
Wang, X., Wan, Q., Ma, G., Li, J., & Fan, J. (2016). The influence of ionospheric thin shell height on TEC retrieval from GPS observation. Research in Astronomy and Astrophysics, 16, 016. https://doi.org/10.1088/1674-4527/16/7/116
Google Scholar
↵
Wang, Y., Zhao, L., & Gao, Y. (2021). Estimation and analysis of GNSS differential code biases (DCBs) using a multi-spacing software receiver. Sensors, 21(2), 443. https://doi.org/10.3390/s21020443
Google Scholar
↵
Wilson, B., & Mannucci, A. (1993). Instrumental biases in ionospheric measurements derived from GPS data. Proc. of the 6th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 1993), Salt Lake City, UT, 1343–1351. https://www.ion.org/publications/abstract.cfm?articleID=4317
Google Scholar
↵
Wu, Q., Zhang, P., Sun, M., Liu, S. , Wang, H., & Chen, S. (2021). Performance evaluation of GIMs released by different IGS ionosphere associate analysis centers in ionospheric constrained single-frequency precise point positioning. Advances in Space Research, 68(12), 4834–4856. https://doi.org/10.1016/j.asr.2020.12.006
Google Scholar
↵
Wu, Z., Li, S., Wan, H., Ji, M., Mao, P., & Xiong, S. (2024). Effects of BDS flex power on DCB estimation and PPP convergence. GPS Solutions, 28, 41. https://doi.org/10.1007/s10291-023-01581-8
Google Scholar
↵
Xiang Y., & Gao Y. (2019). An enhanced mapping function with ionospheric varying height. Remote Sensing, 11(12), 1497. https://doi.org/10.3390/rs11121497
Google Scholar
↵
Xiang, Y., Xu, Z., Gao, Y., & Yu, W. (2020). Understanding long-term variations in GPS differential code biases. GPS Solutions, 24, 118. https://doi.org/10.1007/s10291-020-01034-6
Google Scholar
↵
Zha, J., Zhang, B., Yuan, Y., Zhang, X., & Li, M. (2019). Use of modified carrier-to-code leveling to analyze temperature dependence of multi-GNSS receiver DCB and to retrieve ionospheric TEC. GPS Solutions, 23, 1–12. https://doi.org/10.1007/s10291-019-0895-2
Google Scholar
↵
Zhang, D., Shi, H., Jin, Y., Zhang, W., Hao, Y., & Xiao, Z. (2014). The variation of the estimated GPS instrumental bias and its possible connection with ionospheric variability. Science China Technological Sciences, 57, 67–79. https://doi.org/10.1007/s11431-013-5419-7
Google Scholar
↵
Zhong, J., Lei, J., Dou, X., & Yue, X. (2016). Is the long-term variation of the estimated GPS differential code biases associated with ionospheric variability? GPS Solutions, 20, 313–319. https://doi.org/10.1007/s10291-015-0437-5
Google Scholar
↵
Zhong J., Lei J., & Yue, X. (2020). Comment on Choi et al. correlation between ionospheric TEC and the DCB stability of GNSS receivers from 2014 to 2016. Remote Sens. 2019, 11, 2657. Remote Sensing, 12(21), 3496. https://doi.org/10.3390/rs12213496
Google Scholar
↵
Zus, F., Deng, Z., Heise, S., & Wickert, J. (2017). Ionospheric mapping functions based on electron density fields. GPS Solutions, 21, 873–885. https://doi.org/10.1007/s10291-016-0574-5
Google Scholar

Previous Next

Sensitivity of Receiver Differential Code Bias Estimates to Mapping Functions

Abstract

1 INTRODUCTION

2 MODEL FORMULATION

2.1 Absolute Ionospheric Measurements

2.2 Slant-to-Vertical Mapping Functions

2.2.1 Thin Shell Model

2.2.2 Extended Slab Model

2.2.3 IRI-Based Mapping

2.3 Spatiotemporal Modeling for DCB Estimation

2.4 Analytical Approximation for Mean Receiver DCBs

2.5 Decomposition of Receiver DCB Estimate Variation

3 RESULTS

3.1 Validation and Testing with an IRI-Simulated Ionosphere

3.2 Estimation of Solar Cycle 25 Receiver DCB

3.2.1 Construction of a WRS Proxy Network

3.2.2 Receiver DCB Estimates

3.2.3 Comparison of Analytical Approximate Mean Receiver DCBs With Spatiotemporally Modeled Estimates

3.2.4 Fits to Receiver DCB Estimates

4 DISCUSSION AND CONCLUSIONS

HOW TO CITE THIS ARTICLE

5 ACKNOWLEDGMENTS

Appendix

7 APPENDIX

7.1 Computational Details for Spatiotemporal Modeling

7.1.1 Cycle Slip Detection

7.1.2 Carrier-to-Code Leveling

7.1.3 QR-Decomposition Recursive Least-Squares (QR-RLS) Solver

7.2 IRI Simulation Details

7.3 Goodness-of-Fit Statistics

6 REFERENCES

In this Issue

NAVIGATION: Journal of the Institute of Navigation

Jump to section

Keywords

Sensitivity of Receiver Differential Code Bias Estimates to Mapping Functions

Abstract

1 INTRODUCTION

2 MODEL FORMULATION

2.1 Absolute Ionospheric Measurements

2.2 Slant-to-Vertical Mapping Functions

2.2.1 Thin Shell Model

2.2.2 Extended Slab Model

2.2.3 IRI-Based Mapping

2.3 Spatiotemporal Modeling for DCB Estimation

2.4 Analytical Approximation for Mean Receiver DCBs

2.5 Decomposition of Receiver DCB Estimate Variation

3 RESULTS

3.1 Validation and Testing with an IRI-Simulated Ionosphere

3.2 Estimation of Solar Cycle 25 Receiver DCB

3.2.1 Construction of a WRS Proxy Network

3.2.2 Receiver DCB Estimates

3.2.3 Comparison of Analytical Approximate Mean Receiver DCBs With Spatiotemporally Modeled Estimates

3.2.4 Fits to Receiver DCB Estimates

4 DISCUSSION AND CONCLUSIONS

HOW TO CITE THIS ARTICLE

5 ACKNOWLEDGMENTS

Appendix

7 APPENDIX

7.1 Computational Details for Spatiotemporal Modeling

7.1.1 Cycle Slip Detection

7.1.2 Carrier-to-Code Leveling

7.1.3 QR-Decomposition Recursive Least-Squares (QR-RLS) Solver

7.2 IRI Simulation Details

7.3 Goodness-of-Fit Statistics

6 REFERENCES

In this Issue

NAVIGATION: Journal of the Institute of Navigation

Alerts for this Article

Email This Article

Citation Tools

Sensitivity of Receiver Differential Code Bias Estimates to Mapping Functions

Citation Manager Formats

Share This Article

Sensitivity of Receiver Differential Code Bias Estimates to Mapping Functions

SAVE ITEM IN SELECTED FOLDER.

Jump to section

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords