A Least Squares Solution to Regionalize VTEC Estimates for Positioning Applications

Abstract

A new approach is presented to improve the spatial and temporal resolution of the Vertical Total Electron Content (VTEC) estimates for regional positioning applications. The proposed technique utilises a priori information from the Global Ionosphere Maps (GIMs) of the Center for Orbit Determination in Europe (CODE), provided in terms of Spherical Harmonic (SH) coefficients of up to degree and order 15. Then, it updates the VTEC estimates using a new set of base-functions (with better resolution than SHs) while using the measurements of a regional GNSS network. To achieve the highest accuracy possible, our implementation is based on a transformation of the GIM/CODE VTECs to their equivalent coefficients in terms of (spherical) Slepian functions. These functions are band-limited and reflect the majority of signal energy inside an arbitrarily defined region, yet their orthogonal property is remained. Then, new dual-frequency GNSS measurements are introduced to a Least Squares (LS) updating step that modifies the Slepian VTEC coefficients within the region of interest. Numerical application of this study is demonstrated using a synthetic example and ground-based GPS data in South America. The results are also validated against the VTEC estimations derived from independent GPS stations (that are not used in the modelling), and the VTEC products of international centres. Our results indicate that, by using 62 GPS stations in South America, the ionospheric delay estimation can be considerably improved. For example, using the new VTEC estimates in a Precise Point Positioning (PPP) experiment improved the positioning accuracy compared to the usage of GIM/CODE and Klobuchar models. The reductions in the root mean squared of errors were ∼23% and 25% for a day with moderate solar activity while 26% and ∼35% for a day with high solar activity, respectively.

1. Introduction

The Global Navigation Satellite System (GNSS) technique has become an integral part of all applications, where mobility plays an important role. The basic observable in the GNSS positioning is the time required for electromagnetic signals to travel from the GNSS satellite (transmitter) to a GNSS receiver. This travelling time, multiplied by the speed of light, provides a measure of the apparent distance (pseudo-range) between them. By knowing the position of GNSS satellites from Precise Orbit Determination (POD), the unknown position of the GNSS receiver and its uncertainty can be computed when at least four range measurements exist.

From the mostly used GNSS constellations, GPS and GLONASS satellites orbit at altitudes of around 20,000 km, while BeiDou and Galileo satellites orbit a bit higher, i.e., around 21,500 km and 23,000 km, respectively. The signals from GNSS satellites must transit the ionosphere (i.e., the part of atmosphere between 60 and 2000 km containing ionized plasma of different gas components) on their way to receivers. These free electrons add delay on the code-derived pseudo-range and advance the career phase signals. These effects must be eliminated in some way to achieve high accuracy in GNSS positioning, navigation and timing applications. As a result, the ionospheric modelling has received ever increasing attention in various fields including radio communication, navigation, satellite positioning and other space technologies [1].

Ionospheric models can be divided into three main categories: group (i) includes physical models, in which ionospheric changes are simulated based on the physical laws or the assumptions concerning the structure and variations of ionosphere such as the Global Assimilative Ionospheric Model GAIM [2]; group (ii) is known as empirical models that use deterministic functions to describe periodic and sudden changes in ionosphere, see, e.g., [3], including the International Reference Ionosphere IRI [4] and the NeQuick model [5]; and finally group (iii) consists of mathematical models, which are estimated in terms of mathematical (base) functions using observation techniques that provide ionospheric variables such as the Total Electron Content (TEC), Vertical TEC (VTEC) and Ionospheric Electron Density (IED).

This study follows the computational strategy of (iii), for which the spatial changes of ionospheric density can be formulated as two-dimensional (2D) or three-dimensional (3D) models. The 2D technique is formulated as a Single-Layer Model (SLM), where all free electrons within the ionosphere are concentrated in an extremely thin layer at a constant height. Thus, in the 2D approach, the vertical gradient of electron density changes is not considered [6,7,8,9,10]. The 2D models are often used for estimating total ionospheric delay in (precise) positioning applications. A comprehensive summary of these models can be found in El-Arini et al. [11].

In the 3D modelling techniques, horizontal changes in the ionospheric electron contents (e.g., with respect to the geodetic latitude and longitude), as well as their vertical variations (i.e., along the altitude that is measured from surface of the Earth) are described. Therefore, the 3D models require more observations (than the 2D techniques) and a careful parameterization to obtain a relatively stable relationships between observations and model’s unknown parameters. To gain reasonable horizontal and vertical coverage for the 3D ionospheric tomography, a combination of various TEC observations is considered in previous studies, see e.g., [12,13,14]. This includes STEC (e.g., from GNSS observation), IED (e.g., from Radio Occultation (RO) of Low-Earth-Orbiting (LEO) satellites) and VTEC measurements (e.g., from satellite altimetry) for various techniques of ionospheric density modelling, see, e.g., [15]. While the 3D techniques are often used to study the structure of ionosphere, the 2D models (SLMs) are applied for computing total ionospheric delays in Precise Point Positioning (PPP) applications. It is worth mentioning here that the temporal variation of above mentioned techniques can be achieved by modelling the ionosphere as a short-time (e.g., 2-hourly) static field, see, e.g., [16] or dynamic as in [17]. In this study, we will show how our regionalized 2D (spatial) TEC model can improve the performance of PPP applications.

The single-layer TEC models, see, e.g., [11], are often described based on a set of (mathematical) base functions. The choice of these functions is arbitrary and application-dependent. For example, polynomial functions were chosen in [18,19], while Spherical Harmonics (i.e., often selected to be up to degree and order 15) are considered in, e.g., [20] for global TEC inversions. The Center for Orbit Determination in Europe (CODE) is one of the IGS analysis centers, which determines the precise GPS orbits using the IGS network data and provides the orbit information to GPS users worldwide. The CODE has also produced daily maps of the Earth’s ionosphere on a regular basis since 1 January 1996. The GIM/CODE is modeled by 256 coefficients of the Spherical Harmonics (SHs) expansion up to degree and order 15. Principles of the TEC mapping technique by CODE are described in, e.g., [6,21].

Though modelling based on SHs is well-understood and it is convenient to be applied for the global representation and analyzing the ionosphere, it is not well suited for regional applications [22]. Therefore, regional base functions are introduced to the regional TEC modelling applications. Most of these studies have taken advantage of the orthogonal but strictly band-limited functions that can be concentrated within a region of interest, or by considering an appropriate orthogonal family of the strictly space-limited functions [23]. For example, the B-Spline functions were applied in [14,24,25,26,27], which are based on the Euclidean quadratic B-Spline wavelets. The useful properties of these base functions, i.e., continuity, smoothness and computational efficiency, provide a great advantage for regional modelling of the ionosphere or when the observations are unevenly distributed over the globe. The Spherical-cap Harmonic were applied in [28,29,30] to reduce the lack of orthogonality of the global SHs for regional applications. To mitigate a limitation of these functions that can only be built in regions with smooth boundaries, the regional TEC inversion was formulated in [10] based on the spherical Slepian functions [23] that optimizes field separation over arbitrary regions with irregular boundaries. The other benefit of this formulation is the direct relationship between spherical Slepian functions and the global SHs, which will be discussed in Section 3.

Building on the methodology of Farzaneh and Forootan [10], the main focus of this study is to present an efficient approach to use already available TEC information or maps as a priori information and update them when new TEC observations become available. This update is implemented through a Least Squares (LS) approximation of the Bayesian formulation that considers the uncertainties of both a priori information and observations. We will demonstrate the efficiency of this approach through a synthetic example, where the ‘true’ solutions of the regional VTEC is known by definition. We also test the proposed algorithm by providing an accurate and fast estimation of the regional TECs for a Precise Point Positioning (PPP) application. Therefore, the methodology of this study is formulated as a 2D (spatial) mathematical approach followed by an LS update to regionalize the available TEC maps in the region of interest. Our main assumption is that the local TEC changes are well presented in the dual frequency GNSS observations of the local network. The presented method, after some modifications, can also be extended to the 3D formulation. Theoretically, the presented method can also accept other TEC observations (e.g., from RO, altimetry and the Swarm mission) besides the GNSS observations. These extensions, however, will be a subject of our future investigations.

Recently, Erdogan et al. [17] developed a near real-time processing framework to model spatial and temporal variations within the ionosphere in terms of compactly supported B-spline functions and the recursive filtering using GNSS measurements. The DGFI-TUM’s high-resolution ionospheric product was compared by Goss et al. [31] with the GIM/CODE and the voxel solution from the Universitat Politècnica de Catalunya (UPC). The authors found a better ionosphere representation using DGFI-TUM estimations during the test period of September 2017. Olivares-Pulido et al. [26] presented a real-time TEC modelling approach using B-Spline base functions and a sequential Kalman filter updating scheme, with the goal of providing ionospheric corrections for Precise Point Positioning-Real Time Kinematic (PPP-RTK) applications. Compared to these techniques, the proposed approach of this study addresses a regional estimation of TEC by taking advantage of available a priori information and models the localized anomalies based on the spherical Slepian functions.

In what follows, the data sources of this study are described in Section 2. Estimation of the Slant TEC (STEC) using GNSS observations and the description of global Vertical TEC (VTEC) maps used as a priori information are described in Section 2.1 and Section 2.2, respectively. In Section 2.3, the VTEC products of the NASA’s Jet Propulsion Laboratory (JPL), European Space Agency (ESA) and the International GNSS Service (IGS) are introduced, where they are later compared with the regionalized VTEC estimates of this study and those of GIM/CODE. TEC modelling is introduced in Section 3, where, in Section 3.1, the (spherical) Slepian functions are introduced, and the LS formulation to estimate TECs in the presence of a priori data are described in Section 3.2. Results of this study are presented in Section 4, where the simulation is presented in Section 4.1 and Section 4.2, the Independent Component Analysis (ICA), as in [32,33,34], is applied to compare one year of the regionalized TEC maps with those of GIM/CODE. The regionalized VTEC estimates are compared with those of CODE, JPL, ESA and IGS in Section 4.3, where the VTEC estimates from three GPS stations are used as an independent validation. An assessment of the regionalized TEC estimations in a PPP application is performed in Section 4.4, and, finally, this study is concluded in Section 5.

2. Data

2.1. TEC and VTEC Determination Using GNSS Observations

In theory, the Slant Total Electron Contents (STECs) can be directly computed from the differential code or carrier phase measurements received by dual frequency receivers. Here, the formulation is presented based on the GPS measurements on both L1 (1575.420 MHz) and L2 (1227.600 MHz) frequencies. The noise level of the carrier phase measurements is significantly lower than that of the pseudo-range ones. However, for estimating STECs from the carrier phase measurements, one must account for the (float/integer) full cycle ambiguity (see [35], which is often estimated at the pre-processing step). In order to benefit from the ambiguity-independent estimates of STECs derived from the code pseudo-ranges and the high precision of the carrier phase measurement, the pseudo-range ionospheric observations are smoothed using the “carrier to code leveling process” method [25,36], i.e.,

$$STEC=({\tilde{P}}_4-b_r-b_s-{〈{\epsilon }_p〉}_{arc}+{\epsilon }_L)\left(\frac{f_1^2{f}_2^2}{40.3(f_2^2-f_1^2)}\right),$$

(1)

where ${\tilde{P}}_4$ is the pseudo-range ionospheric observable smoothed by the carrier-phase ionospheric one, i.e.,

$${\tilde{P}}_4={〈P_4+{\Phi }_4〉}_{arc}-{\Phi }_4\approx I_1-I_2+br+bs+{〈{\epsilon }_p〉}_{arc}-{\epsilon }_L.$$

(2)

In these equations, $b_r$ and $b_s$ are the code Inter-Frequency Biases (IFBs) for the receiver and satellite, $f_1$ and $f_2$ are the L1 (1575.420 MHz) and L2 (1227.600 MHz) frequencies, $I_1$ and $I_2$ are the ionospheric refraction delays at L1 and L2, and $P_4$ and ${\Phi }_4$ are the geometry-free linear combination of pseudorange and carrier phase measurements in the continuous observational arc (the interval at which no cycle-slip error has occurred). Finally, ${\epsilon }_p$ and ${\epsilon }_L$ represent the effects of multi-path and measurement noise on the pseudo-range and carrier phase, respectively. The STEC values can be converted into the height-independent VTEC by introducing the single layer mapping function [20]:

$$MF=\frac{STEC}{VTEC}=\frac{1}{\cos z^{\prime }},$$

(3)

with

$$\sin z^{\prime }=\frac{R_E}{R_E+H}\sin z,$$

(4)

where R is the mean Earth radius, z and $z^{\prime }$ are respectively the zenith angles of GPS satellite at the user position and the ionospheric pierce point, and H is the mean altitude (that approximately corresponds to the altitude of the maximum electron density and its height, i.e., varying between 250 and 500 km) depending on the latitude, season, solar and geomagnetic activity conditions [37,38].

The uncertainty of GPS-derived VTECs is related to the quality of code pseudo-range and carrier phase measurements at $L1$ and $L2$ frequencies, and phase and code inter-frequency biases for the receiver and satellites. The usage of the mapping function to covert STEC to VTEC (Equation (3)) can also be considered as an additional source of uncertainties, but this has not been considered in this study. According to covariance propagation theory, these uncertainties can be estimated as follows:

$${\sigma }^{VTEC-GNSS}=\frac{\alpha }{MF}\times \sqrt{{\sigma }_{{\tilde{P}}_4}^2+{\sigma }_{br}^2+{\sigma }_{bs}^2},\alpha =\frac{f_1^2{f}_2^2}{40.3(f_2^2-f_1^2)},$$

(5)

where ${\sigma }_P={\sigma }_{P_1}={\sigma }_{P_2}=0.2$ m is standard deviation of code pseudo ranges observation and ${\sigma }_{\varphi }={\sigma }_{{\varphi }_1}={\sigma }_{{\varphi }_2}=0.02$ cycle is the standard deviation of carrier phase pseudo ranges observation. Assuming that the code pseudo-range and carrier phase derived TECs are uncorrelated, uncertainties of the pseudo-range ionospheric observable smoothed ${\tilde{P}}_4$ can be derived using the formal error propagation law, which gives:

$${\sigma }_{{\tilde{P}}_4}^2={\sigma }_{\varphi }^2\times {\lambda }_1^2+{\sigma }_{\varphi }^2\times {\lambda }_2^2+\frac{1}{n}({\sigma }_{\varphi }^2\times {\lambda }_1^2+{\sigma }_{\varphi }^2\times {\lambda }_2^2+2\times {\sigma }_P^2),$$

(6)

where ${\lambda }_1$ and ${\lambda }_2$ are the wavelength of carrier phase observations (i.e., $19.03$ cm and $24.42$ cm), n is the number of measurements in the continuous arc, while ${\sigma }_{br}$ and ${\sigma }_{bs}$ are reachable from the IONEX files, which are produced by the IGS Analysis Centers (ftp://ftp.unibe.ch/aiub/CODE/) and contain the GPS-derived TEC maps with their uncertainties.

2.2. Global Ionospheric Model (GIM/CODE)

The Global Ionosphere Maps (GIMs) from the Center for Orbit Determination in Europe (CODE) are modeled by 256 coefficients of the Spherical Harmonics expansion up to degree and order 15. The 2-hourly GIM/CODE coefficients (${\tilde{C}}_{nm}$ and ${\tilde{S}}_{nm}$) and their errors ($\sigma \left({\tilde{C}}_{nm}\right)$ and $\sigma \left({\tilde{S}}_{nm}\right)$) during 2013 are downloaded from ftp://ftp.unibe.ch/aiub/CODE/. We compute the VTEC values and their standard deviations form the GIM/CODE data ($VTE{C}^{GIM/CODE}$ and ${\sigma }_i^{VTEC-GIM/CODE}$) in the geomagnetic frame (i.e., latitude ${\phi }_i$ and $S_{Si}={\lambda }_i-{\lambda }_s$ with ${\lambda }_i$ being the longitude of the point of interest and ${\lambda }_s$ representing the geomagnetic longitude of the Sun) as:

$$VTE{C}^{GIM/CODE}=\sum _{n=0}^{n_{max}}\sum _{m=0}^n{\tilde{P}}_{nm}\left(\sin {\phi }_i\right)({\tilde{C}}_{nm}\cos \left(m{s}_{si}\right)+{\tilde{S}}_{nm}\sin \left(m{s}_{si}\right)),\mathrm{and}$$

$${\sigma }_i^{VTEC-GIM/CODE}=\sum _{n=0}^{n_{max}}\sum _{m=0}^n{\tilde{P}}_{nm}\left(\sin {\phi }_i\right)(\sigma \left({\tilde{C}}_{nm}\right)\cos \left(m{s}_{si}\right)+\sigma \left({\tilde{S}}_{nm}\right)\sin \left(m{s}_{si}\right)),$$

(7)

where ${\tilde{P}}_{nm}$ are normalized Legendre polynomials. Equation (8) can be written as a linear transformation in the form of ${\sigma }_i^{VTEC-GIM/CODE}=\mathbf{M}{\Sigma }^X$, where ${\Sigma }^X$ contains the errors ($\sigma \left({\tilde{C}}_{nm}\right)$, and $\sigma \left({\tilde{S}}_{nm}\right)$) and $\mathbf{M}$ contains the known values of ${\tilde{P}}_{nm}$, $\cos \left(m{s}_{si}\right)$ and $\sin \left(m{s}_{si}\right)$. The covariance matrix of GIM/CODE data in the grid domain can be derived as:

$${\Sigma }^{VTEC-GIM/CODE}=\mathbf{M}{\Sigma }^X{\mathbf{M}}^T.$$

(8)

2.3. Ionospheric Models Used for Comparisons

VTEC estimates from the official products of the NASA’s Jet Propulsion Laboratory (JPL), the European Space Agency (ESA) and the International GNSS Service (IGS) are used here for comparison. The IGS ionosphere working group was established since 1998, and has developed different techniques to provide VTEC maps using the IONosphere EXchange (IONEX) format [39,40,41] since then. IONEX provides TEC estimates with a spatial resolution of 5${}^{\circ }$ and 2.5${}^{\circ }$ in longitude and latitude, respectively, and a temporal resolution of few minutes to several hours in real-time, rapid and final modes. Although the real-time TEC products have been proposed by the IGS, users now can only access the rapid and final products with a latency of a few days. Various VTEC products from CODE, Universitat Politècnica de Catalunya/IonSAT (UPC), JPL, ESA are accessible from ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex/ with 2 h steps [20,21,42,43]. Comparing TEC estimates from different centers can reflect the consistency of our estimates with respect to the existing modelling methods.

3. Method

In what follows, in Section 3.1, we introduce the (spherical) Slepian functions that can be used for efficient regional TEC modelling, see, e.g., [10]. The relationships between these functions and SHs are also presented. In Section 3.2, a Least Squares (LS) approximation of the Bayesian-type update is provided to compute TEC maps using GNSS observations, while taking the GIM/CODE maps (in terms of SHs) as a priori information.

3.1. Spherical Slepian Functions

A global field can be expanded by SH functions and their coefficients as:

$$f\left(r\right)=\sum _{l=0}^{\infty }\sum _{m=-l}^l{f}_{lm}{Y}_{lm}\left(r\right),$$

(9)

where $Y_{lm}\left(r\right)$ represents SHs of degree l and order m, r is the location of a point on the surface of a unit sphere $\Omega$. The corresponding coefficients ($f_{lm}$) can be estimated by solving, e.g., an integral over the entire sphere (unit sphere $\Omega$), i.e.,

$$f_{lm}=\frac{1}{4\pi }{\int }_{\Omega }f\left(r\right)Y_{lm}\left(r\right)\mathrm{d}\Omega .$$

(10)

In order to localize these functions into a region of interest (target region), the optimization of a local energy criterion can be utilized. This will give a new set of functions in the sense of [44]. The spherical Slepian function can be presented as a band-limited spherical harmonic expansion as:

$$g\left(r\right)=\sum _{l=0}^L\sum _{l=-m}^m{g}_{lm}{Y}_{lm}\left(r\right),$$

(11)

with

$$g_{lm}={\int }_{\Omega }g\left(r\right)Y_{lm}\left(r\right)\mathrm{d}\Omega .$$

(12)

To maximize the spatial concentration of the band-limited function $g\left(r\right)$ within the target region R, the ratio of the norms should be maximized as:

$$\lambda =max\frac{{\int }_R{g}^2(\Omega )d\Omega }{{\int }_{\Omega }{g}^2(\Omega )d\Omega },$$

(13)

where $0\le \lambda \le 1$ is a measure of spatial concentration. The maximization of this concentration criterion can be achieved in the spectral domain by solving the algebraic eigenvalue problem as:

$$\sum _{l^{\prime }=0}^L\sum _{m^{\prime }=-l^{\prime }}^{l^{\prime }}{D}_{lm,l^{\prime }{m}^{\prime }}{g^{\prime }}_{l^{\prime }{m}^{\prime }}=\lambda g_{lm},$$

(14)

with $D_{lm,l^{\prime }{m}^{\prime }}$ being the Gram matrix of energy within the target region R, i.e.,

$$D_{lm,l^{\prime }{m}^{\prime }}=\underset{R}{\int }{Y}_{lm}{Y}_{l^{\prime }{m}^{\prime }}d\Omega .$$

(15)

Therefore, the desired signal within the region of interest R can be efficiently approximated by:

$$d\left(\widehat{r}\right)\approx \sum _{n=1}^N{d}_n{g}_n\left(\widehat{r}\right),$$

(16)

where $g_n\left(\widehat{r}\right)$, $d_n$ and N are the spherical Slepian functions, corresponding unknown coefficients and Shannon number, respectively; see more details in [45,46]. A linear transformation can be used to convert SH coefficients into spherical Slepian ones, whose energy is concentrated onto specific patches of the sphere [22,23,47]. Considering f to be the coefficients of SHs, which can be, for example, those of VTEC coefficients from the GIM/CODE, and a given variance-covariance matrix ${\Sigma }^X$; this global field can be localized inside the region of interest using:

$$d_n^{\mathrm{localized}}=\left(g{g}^T\right)f,$$

(17)

where $d_n^{\mathrm{localized}}$ is the localized field and $g={\left(g_{00}\dots g_{lm}\dots g_{LL}\right)}^T$ is a localization matrix from Equation (10). The covariance matrix of the localized GIM/CODE coefficients, which is used as the prior covariance matrix in the update step (Section 3.2), can be estimated through a covariance propagation as:

$${\Sigma }_{d_n^{\mathrm{localized}}}=\left(g{g}^T\right){\Sigma }^X\left(g^{T}g\right).$$

(18)

3.2. A Least Squares (LS) Approximation of the Bayesian Update

In Equation (17), we showed how to convert SH coefficients to their corresponding spherical Slepian coefficients. These values and their errors can then be used as a priori information in a Bayesian formulation to be updated by the VTEC values that are estimated from GNSS observations (or any other techniques that measure TEC or VTEC). To estimate these updates, suppose L is a vector of VTEC observations (e.g., computed by rearranging Equation (3), i.e., $VTEC=STEC\cos z^{\prime }$. The distribution of L (i.e., $P\left(L\right)$) conditionally relates to the distribution of the unknown Slepian coefficients $d_n$ (i.e., shown by $P\left(d_n\right)$). Relationships between observation and unknowns are introduced by $P\left(L\right|d_n)$, which is known as the likelihood function and the distribution of unknown parameters in the presence of observation L (i.e., shown by $P\left(d_n\right|L)$) is derived by the Bayesian theory, i.e.,

$$P\left(\left.d_n\right|L\right)\propto P\left(d_n\right)P\left(\left.L\right|d_n\right).$$

(19)

The distribution of unknown parameters ($P\left(d_n\right)$) is already known before the observations (L) were taken. Once the observations (L) are introduced, $P\left(d_n\right|L)$ represents the posterior distribution of the parameter vector ($d_n$). Thus, this is an update of a priori information by the introduced observations (L). By knowing the distribution of parameters ($P\left(d_n\right)$), one can compute the mathematical expectation of parameters, i.e., $d_n^{\mathrm{localized}}$ and its covariance matrix, i.e., ${\Sigma }_{d_n^{\mathrm{localized}}}$. Here, for simplicity, we suppose that $P\left(d_n\right)$ is normal. Thus, a priori distribution of the unknown parameters $d_n$ is $P\left(d_n\right)\sim N(d_n^{\mathrm{localized}},{\Sigma }_{d_n^{\mathrm{localized}}})$. Moreover, if we assume that the observation vector (L) is also normally distributed and the variance factor to be ${\sigma }^2$, a priori distribution of the unknown parameters ($d_n$) is a conjugate prior. Hence, the posterior distribution of $d_n$ is also normal and can be computed as: $d_n|L\sim N(d_0,{\sigma }^2{(A^{T}PA+{\Sigma }_{d_n^{\mathrm{localized}}}^{-1})}^{-1})$, where A is the design matrix containing the Slepian functions, i.e., g in Equation (12), P is the weight matrix of the observations, and $d_0$ is the posterior mathematical expectation. Therefore, the LS approximation of the Bayesian update that provides $d_n$ is given by:

$${\widehat{d}}_n={(A^{T}PA+{\Sigma }_{d_n^{\mathrm{localized}}}^{-1})}^{-1}(A^{T}PL+{\Sigma }_{d_n^{\mathrm{localized}}}^{-1}{d}_n^{\mathrm{localized}}),$$

(20)

with A being a full column rank $n\times u$ matrix. Here, n is the total number of observations (length of the observation vector L); u is the number of unknowns (i.e., the total number of unknown Slepian coefficients); P is the known positive definite $n\times n$ weight matrix of the observations. The variance of the unknown Slepian coefficients can be computed as:

$${\widehat{\sigma }}^2=\frac{V^{T}PV+V_X^T{\Sigma }_{\overline{X}}^{-1}{V}_X}{n},$$

(21)

where $V=A\widehat{d_n}-L$, $V_X=\widehat{d_n}-d_n^{\mathrm{localized}}$, ${\Sigma }_{\overline{X}}^{-1}$ is the prior covariance matrix and n is the total number of observations. An overview of VTEC estimation using the proposed approach is presented in Figure 1.

4. Results and Discussion

The regional VTEC modelling of this study is based on the ground-based GNSS observations collected across South America and few stations in North America. GNSS observations of 62 stations belong to IGS and the Brazilian Network for Continuous Monitoring (RBMC). The data are obtained from www.ibge.gov.br with the sampling rate of 30 s. In order to solve STEC from these observations, receivers’ Differential Code Bias (DCB) and the Inter-Frequency Bias (IFB) values for the GNSS satellites were calculated using the regional formulation as in [48,49,50].

The STEC and VTEC values from each GNSS observation were computed using Equations (1)–(3). The altitude of the Single Layer Model (SLM) was set to 450 km (to be consistent with those of GIM/CODE), and the elevation cut-off angle of ${15}^{\circ }$ was used to select valid GNSS satellites. The precise orbit files, which are provided the IGS agencies, were downloaded from ftp://cddis.gsfc.nasa.gov/pub/gps/products and interpolated to determine satellite positions. An overview of the input data used for the ionospheric modelling of this study is shown in Figure 2, where one can see that our VTEC modelling domain covers an extended region above the the GNSS network.

In what follows, we validate the proposed approach by a synthetic example (Section 4.1). The VTEC estimates are then assessed in three different ways, in Section 4.2, the 2-hourly regionalized maps of this study are compared with those of GIM/CODE to see whether the new model represents expected spatial-temporal as reflected in the global model. In Section 4.3, the regionalized vTEC estimates are compared with the predicted values provided by the NASA’s JPL, European Space Agency (ESA), IGS and GIM/CODE. Finally, in Section 4.4, the VTEC estimations are assessed in a Precise Point Positioning (PPP) application.

4.1. Simulation

In order to validate our modelling approach, a synthetic example is designed to assess the ability of the Bayesian approach in recovering the regional signals. To introduce the ‘true’ VTEC patterns, the spherical harmonic coefficients of the GIM/CODE is used to produce a smooth grid with 0.5 degree resolution that covers South America. Then, we added periodic oscillations with the magnitude of 10 TECU $\left(F\right(longitude,latitude)=10\ast \sin (20longitude)\ast \cos (20latitude\left)\right)$ on the top of GIM (see Figure 3, plot on the left). For comparison, we show a simple synthesizing of the true signal using the spherical harmonics expansions of degree and order 15 and 90 in Figure 3A1,A2, respectively. Their differences with the true pattern are shown in B1 and B2, respectively. The regionalize VTEC estimation method is implemented by considering the low-degree VTEC estimates of Figure 3A1 as a priori information.

Two regionalization experiments are considered, where the first is shown in Figure 3A3. Here, the ‘true’ VTECs are interpolated at 2140 points that are located at the pierce points (see the locations in Figure 3C1) that are derived by connecting the 62 stations of Figure 2 to the GPS satellites during 30 s. In Figure 3A4, 1${}^{\circ }$ VTEC values of the true VTEC signal (i.e., 8191 points of Figure 3C2) are considered as observation. The differences with the true signal are shown in Figure 3B3,B4, respectively. Though the VTEC recovery by spherical harmonics is a simple synthesizing the introduced VTEC field, its accuracy is found to be limited due to the truncation that is dictated by the selected spectral resolution. This is demonstrated in the residual plots (Figure 3B1,B2, where the differences of 10 TECU can be detected that include mismatch anomalies compared with the VTEC estimates from the GIM/CODE, as well as the artificial regional anomalies. The magnitude of residuals derived from the Bayesian update using only 30 s of the local GPS network is one level of magnitude better than Figure 3A2,B2. In Figure 3A4, it can be seen that using a well covered VTEC observations results in an accurate recovery of the truth with very negligible error magnitude (i.e., ${10}^{-4}$ in Figure 3B4). This simulation indicates that indeed the regional anomalies are able to modify a priori information. Therefore, the method can will be applied to real case studies.

4.2. A Comparison between the Regionalized VTEC Maps and GIM/CODE Products

Two-hourly snapshots of the regionalized VTEC maps for 17 March (DOY 76) and 21 December 2013 (DOY 355) are presented in 12 maps presented on the top panel of Figure 4 and Figure 5, respectively. These values are presented as TEC Unite (TECU). Considering these maps, the ionosphere maximum appears around the local noon as travelling along with the Sun. Relatively bigger values estimated in DOY 76 (compared to DOY 355) are related to the higher magnetic activity during this day (i.e., the $K_p$ values of these days were $+6$ and $-2$, respectively). Equatorial TEC anomalies can be detected in both days, where a sunrise enhancement is seen in VTEC estimates at 12:00 UT during both days. The high values of TEC at 19:00 and 23:00 UT are related to the physical structures of the diurnal equatorial ionization anomaly and its resurgence after sunset, respectively [51,52,53].

The differences between the regionalized estimates of VTECs and those of GIM/CODE are shown on the bottom panel of Figure 4 and Figure 5. On 17 March (DOY 76) at 00:00, 04:00, 06:00, 20:00 and 22:00 UT large differences in VTEC estimates are found around ${10}^{\circ }\mathrm{S}$–${30}^{\circ }\mathrm{N}$. On 21 December 2013 (DOY 355), the same magnitude of differences is detected at 00:00 and 22:00 UT. The results on 17 March 2013 indicate a minimum VTEC of ∼16 around 8:00 a.m. and ∼95 around 8:00 p.m. (see Figure 4). An average value of VTEC during March (with high magnetic activity) is found to be ∼60. On 21 December 2013, due to its lower magnetic activity, the values VTEC are found to be smaller, i.e., the minimum VTEC of ∼10 around 8:00 AM and ∼68 around 8:00 PM (see Figure 5). Considering the differences between the regionalized solutions and the GIM/CODE (12 plots on the bottom of Figure 4 and Figure 5), the GIM/CODE products are found to underestimate VTEC variations, mostly when the magnitude of VTEC is higher. The maximum differences are found to be ∼15 TECU on 17 March 2013 and ∼10 TECU on 21 December 2013.

Here, we extend the assessment of the new VTEC estimates for the entire year 2013, by applying the Independent Component Analysis ICA [32,33] on the differences between the two-hourly GIM/CODE products and regionalized results. The first two dominant independent modes are shown in Figure 6 that correspond to $65\%$ of the total variance of VTEC differences. Signals on the oceans are masked to highlight the changes over the land and to provide an indication for possible impacts on positioning applications. In each mode, spatial functions (Figure 6-left) are anomaly maps in terms of TECU, which can be multiplied by the unit-less time series (Figure 6-right) known as normalized Independent Components (ICs) to derive independent modes of variability. The results indicate that the magnitude of differences in the year 2013 reach up to 6 TECU, and they are dominated by the diurnal and semi-diurnal frequencies. By analysing the temporal ICs, it became clear that the magnitude of differences during May to middle September (DOY ∼120–258) is almost half of the VTEC differences during the rest of 2013 (see IC1 and IC2 of Figure 6). To assess whether there is a relationship between the VTEC differences and geomagnetic activity, the ICs are smoothed by “loess” and “rloess” methods [54], while using $5\%$ of the data and they are compared with the loess-smoothed geomagnetic $K_p$ index. The numerical results indicate a considerably high correlation coefficient ($0.7$) between them during January to middle May and middle September to December, where the magnitude of VTEC changes was higher than rest of the year.

In Figure 7, the amplitude of diurnal and semi-diurnal VTEC differences is shown. These two frequencies are selected because of their dominance as it was shown in Figure 6. The amplitude of diurnal differences reach up to 5 TECU (Figure 7-left), while that of semi-diurnal is found to be 3 TECU (Figure 7-right). These differences represent considerable impact on the accuracy of positioning applications, where, roughly speaking, 1 TECU and corresponds to 16 cm positioning error at the $f_1$ frequency (1575.420 MHz).

4.3. Comparing Regionalized VTEC Maps with the JPL, ESA and IGS Products

In this section, we evaluate the regionalized VTEC estimates by comparing them with those derived from other groups. For this, TEC values along the line of sights between GNSS stations and satellites are computed (Equation (1)). These values are then converted to VTEC by implementing the single layer mapping function (Equation (3)) [20], and the results are compared with other models. As reference stations, we used the dual frequency GPS observations of Bogota (lat = 4.64007 and lon = 285.91906), Unsa (lat = −24.72746 and lon = 294.59236) and Punta Arenas (lat = −53.13695 and lon = 289.12011). These observations were not used in the regionalized VTEC estimates of this study, but they are used for validations (locations of the stations are shown by three black dots in Figure 1).

Plots in Figure 8 top-left and top-right show the observed VTEC estimates (from dual frequency GPS observations) of the three stations during 17 March 2013 (DOY 76 with $K_p=+6$) and 21 December 2013 (DOY 355 with $K_p=-2$), respectively. For comparison, the differences between these values and the regionalized VTEC estimates, as well as those of GIM/CODE, JPL, ESA and IGS centers [21] are shown in separate plots. The regionalized VTEC estimates are computed in near-real time mode whenever the observations are available. Therefore, they are generated with the same sampling rats of the GPS observation (i.e., every 30 s), while the official products are provided 2-hourly (120 min) and with the latency of several days. By comparing the results of Figure 8, the estimated VTEC residuals of the regionalized model are found to be smaller than the other products.

Table 1. Deviation of the regionalized VTEC estimates, as well as those of GIM/CODE, ESA, IGS and JPL from the GPS derived VTEC values.

Epoch	Regionalized VTEC	GIM/CODE	ESA	IGS	JPL
2013 DOY 76 (Bogota)	4.36	10.13	13.46	7.52	13.70
2013 DOY 76 (Punta Arenas)	3.53	4.81	15.15	9.50	30.76
2013 DOY 76 (Unsa)	6.37	10.14	24.55	20.25	27.27
2013 DOY 355 (Bogota)	3.30	14.96	9.27	10.33	6.94
2013 DOY 355 (Punta Arenas)	2.34	2.56	11.34	4.70	9.34
2013 DOY 355 (Unsa)	2.95	4.56	5.03	4.57	11.57

provides the Average of the Absolute Percentage Deviations ($AAPD=\frac{{\sum }_{i=1}^{size\left(obs\right)}\left(100\times \left|\frac{ob{s}_i-mode{l}_i}{ob{s}_i}\right|\right)}{\mathrm{size}\left(obs\right)}$) of the regionalized VTEC estimates and those of GIM/CODE, ESA, IGS and JPL values from the GNSS estimates.

4.4. A PPP Assessment of the Regionalized VTEC Estimates

Three GNSS stations (i.e., Bogota, Unsa and Punta Arenas) of Section 4.3 are chosen to assess the impact of VTEC estimates on the positioning accuracy. The dates of assessments are also chosen to be the same as previous section to make the interpretation easier. As a measure of accuracy, the Root Mean Squares Error (RMSE) of the positioning residuals is calculated. To compute an accurate position (to be used as reference for estimating the position accuracy), the Precise Point Positioning (PPP) solution as in [55] is chosen as our computation technique. Ionospheric-free combination is created using GPS observations including both L1 and L2 signals. Based on this, station coordinates, receiver clock error, systems time difference parameters with the GPS system, troposphere parameter and phase ambiguity are estimated during 24 h of 17 March 2013 (DOY 76 with $K_p=+6$) and 21 December 2013 (DOY 355 with $K_p=-2$).

Table 2. An overview of the observations, models and processing strategy used for the precise Point Positioning (PPP) experiment of this study.

Items	Setting
Observations	Undifferenced phase and code ionosphere-free combination of measurements
Frequency	GPS: L1/L2
Elevation Cutoff	10 (Deg)
Estimator	The recursive least squares with the unknowns added [56]
Sampling rate	5 (Sec)
Satellite orbit and clock	GFZ multi-GNSS (GBM) [57]
Satellite antenna PCO and PC	GPS: IGS14.atx [58]
Receiver antenna phase center and variation	Corrected by IGS14.atx [58]
Troposphere Dry component	Saastamoinen model and Global Mapping Function GMF, [59]
Relativistic effects	Corrected [60]
Phase wind-up	Corrected [61]
Station displacement	Solid Earth tides, ocean tide loading and pole tides [62]
Observation weight	Using the Recursive Least Square Estimation- Variance Components Estimation (RLS-VCE) method [56] to estimate the accuracy of observations and the correlations between them
Estimated parameters	Receiver position, tropospheric wet delay, ambiguity parameters and system time difference parameters

presents the processing strategy and the error modelling for the performed PPP experiment. It is worth mentioning that the convergence period of the PPP experiments was not considered in the computation of the quality measures.

To assess the impact of VTEC modelling on the position accuracy, the above PPP experiment is repeated with the same setup but, for estimating the ionospheric delays, the regionalized VTEC estimates, and those of Klobuchar [63] and GIM/CODE are replaced. The ionosphere-free combination is therefore replaced by the single-frequency PPP [64]. The position estimates of these experiments are compared with those of the reference (computed by the ionosphere-free combination of L1 and L2 measurement as described before).

Figure 9 and Figure 10 show the computed position errors on 17 March 2013 (DOY $=76$) and 21 December 2013 (DOY $=355$), respectively. Plots indicate that indeed the use of regionalized VTECs improves the positioning accuracy in comparison with the GIM/CODE and Klobuchar models. The error plots that correspond to the regionalized VTEC estimates are found to be very smooth similar to those of the ionosphere-free combinations.

Figure 11 summarizes the RMSE of the positioning residuals for each station compared with errors from the dual frequency PPP estimates. In comparison with GIM/CODE and Klobuchar models, the use of regionalized VTEC estimates improves the positioning accuracy by 30% and 33% for Bogota, 25% and 38% for Punta Arenas, as well as 24% and 42% for Unsa during 17 March 2013. These values are found to be 27% and 24% for Bogota, 15% and $23\%$ for Punta Arenas, as well as 28% and 29% for Unsa during 21 December 2013, respectively. The differences in the magnitude of improvements are related to the differences in geomagnetic activity of these two days.

5. Summary and Conclusions

The ionosphere is the major error source that affects the positioning accuracy of GNSS positioning. In this study, a Least Squares (LS) approximation of the Bayesian formulation is introduced to use a priori information from, e.g., already available VTEC maps from the Global Ionosphere Map (GIM) of the Center for Orbit Determination in Europe (CODE). Then, we use TEC estimates from local GNSS networks to update a priori values within the region of interest. The presented VTEC estimation follows a 2-step algorithm, where, in step-1, the GIM/CODE’s VTEC values are transformed from the global spherical harmonic coefficients to an optimum band-limited local spherical Slepian coefficients in the region of interest. In step-2, we use new VTEC observations from GNSS observations in a Bayesian equation to update the spherical Slepian VTEC coefficients of step-1. The numerical assessments are performed on a network in South America including 62 stations. Comparisons are performed with the VTEC products of GIM/CODE, and the external data of JPL, ESA and IGS. A Precise Point Positioning (PPP) experiment is implemented to compare the impact of VTEC models in terms of position errors with the position derived from 24-hour double-frequency measurements of three IGS stations. The main findings of this study include:

• Comparisons with the GIM/CODE confirm that the regionalized model estimates TEC without unexpected oscillations, though the range of variations from the IGS models is found to be underestimated.
• Comparing the regionalized estimates of this study with the VTEC estimations using the dual-frequency measurements of three GPS stations indicates that the average of absolute differences is less than 2 TECU, which indicates an accepted performance of the presented technique.
• Performing VTEC analyses for the entire 2013 shows that the presented regionalization technique is appropriate for VTEC modelling under normal and high geomagnetic conditions.
• Comparison with various VTEC models, the quality of the new estimates, and hence the ionospheric corrections, is found to be better within near real-time PPP applications.
• The results showed that the positioning accuracy of single-frequency positioning with the external ionospheric model correction can obtain meter-level accuracy, and the vertical error is found to be relatively larger than the horizontal components.
• Results indicate that the regionalized model is better suited to correct ionospheric impact of GPS positioning compared with the usage of Klobuchar and GIM/CODE in a precise point positioning setup.

Some ideas for further development and improvement of the results for future work are:

• The new European satellite navigation system, Galileo and the restored Russian system, GLONASS, are examples of other constellations that can double the quantity of (V)TEC data. The multi-constellation observations will be used for future studies to improve the quality of TEC observations.
• Combining different data sources, e.g., from radio occultation and satellite altimetry, will be considered to improve the spatial coverage of TEC observations.
• Further investigations need to be conducted for other GNSS networks at different latitudes with higher or lower reference station density.
• The LS approximation of the Bayesian update can be replaced by a more efficient Markov Chain Monte Carlo optimization to avoid the assumption of Gaussian distribution for a priori information and observations.