Nadir-Dependent GNSS Code Biases and their Effect on 2D and 3D Ionosphere Modeling

: Recent publications have shown that group delay variations are present in the code observables of the BeiDou system, as well as to a lesser degree in the code observables of the global positioning system (GPS). These variations could potentially affect precise point positioning, integer ambiguity resolution by the Hatch–Melbourne–Wübbena linear combination, and total electron content estimation for ionosphere modeling from global navigation satellite system (GNSS) observations. The latter is an important characteristic of the ionosphere and a prerequisite in some applications of precise positioning. By analyzing the residuals from total electron content estimation, the existence of group delay variations was confirmed by a method independent of the methods previously used. It also provides knowledge of the effects of group delay variations on ionosphere modeling. These biases were confirmed both for two-dimensional ionosphere modeling by the thin shell model, as well as for three-dimensional ionosphere modeling using tomographic inversion. BeiDou group delay variations were prominent and consistent in the residuals for both the two-dimensional and three-dimensional case of ionosphere modeling, while GPS group delay variations were smaller and could not be confirmed due to the accuracy limitations of the ionospheric models. Group delay variations were, to a larger extent, absorbed by the ionospheric model when three-dimensional ionospheric tomography was performed in comparison with two-dimensional modeling.


Introduction
Signals in the frequency bands employed by global navigation satellite systems (GNSSs) are affected by the ions and free electrons in the part of the atmosphere referred to as the ionosphere [1]. As the ionosphere is dispersive, GNSSs have been designed to employ multiple carrier frequencies. This allows almost complete elimination of the effect from GNSS code and phase observations. Conversely, availability of observations associated with multiple carrier frequencies also allow for determination of the number of free electrons and ions in the ionosphere by forming the geometryfree linear combination [1].
By utilization of the geometry-free linear combination of GNSS observations, the total electron content (TEC) along the ray paths can be estimated. Doing this for several satellites and several receivers allows for estimation of the spatial distribution of TEC. A two-dimensional (2D) vertical TEC (VTEC) distribution can be estimated and distributed as a global ionospheric map (GIM) where the modelling is based on the "thin shell" model and a mapping function between slant ionospheric delays and vertical delays [1]. A more complete three-dimensional (3D) model might be derived by computerized ionospheric tomography (CIT) [2][3][4], techniques that were inspired by and have emerged from computerized tomography (CT) for medical applications [5].

Methods
The following sections provide a summary of the methods utilized in this work for 2D and 3D ionosphere modeling, as well as for independent derivation of nadir-dependent biases from the MP linear combination.

Derivation of STECs from GNSS Measurements
GNSS signals are affected by free electrons in the ionosphere during their transmission between satellite and receiver. Slant total electron content (STEC) is given by e STEC N ds   (1) where Ne is the ionospheric electron density (IED) along the ray path of the GNSS signal. STEC is a measure of the total effect of free-charged particles of the ionosphere on the GNSS signal. It is related to the GNSS observables as a delay for code observations and an advance for carrier phase measurements. The size of the effect is dependent on the carrier frequency for the frequency bands employed. This enables either isolation or elimination of the ionospheric effect by linearly combining GNSS observations associated with several carrier frequencies. Isolation of the first-order term can be achieved by the geometry-free linear combination of code or phase observations associated with two carrier frequencies [7] In these equations, a quantity , ( ) s f r  is associated with carrier frequency f, receiver r, and satellite s; R is a code observable; L is a phase observable (in units of length); I is the ionospheric influence; DCB is the differential code bias (DCB) that is due to code biases caused by both the receiver and the satellite hardware; M is code multipath; B is combined receiver, br, and satellite phase biases, b s , lumped together with the phase ambiguities, N; m is phase multipath; and ε is noise and other unhandled errors.
In this study, STECs were estimated considering only the first-order ionospheric term, derived with Equations (2), (3), and (6). This is acceptable, as the higher-order terms only amount to one centimeter at most [29], which is well below the expected accuracy for ionosphere modeling.
Considering only the first-order term, the TEC is related to ionospheric influence for code and phase observations by where I is the ionospheric influence from Equation (2) and (3), and f is the carrier frequency [30].
Before the derivation of STECs, the DCBs were calibrated for each receiver with a GIM from the International GNSS Service (IGS) [31] over a period of 24 hours. One receiver DCB for each of the reference stations and one satellite DCB for each of the satellites were estimated by least-squares adjustment in accordance with Equation (2). Multipath and noise were assumed to average out over the 24-hour observation period. For this calibration, the DCBs were considered as constant over time. Possible nadir dependence of code biases will thereby show up in the later stage as residuals from the ionosphere modeling. To avoid the complexity in handling additional GLONASS related biases, only GNSS observations from GPS, Galileo, and BeiDou were employed in the subsequent analysis. For BeiDou, only observations from BeiDou-2 satellites were available at the time of this study.

Ionosphere Modeling in 2D
One common way to model the ionosphere is by an assumption that all the electron content of the ionosphere is distributed over an infinitely thin layer in the so-called "thin shell" model. The height of the shell is preferably chosen to the height of the largest concentration level of free-charged particles. STECs that refer to the slant signal path are in this model converted to vertical TEC (VTEC). This is done with the ionospheric mapping function for a thin shell [1] 1 ( ) cos where ( ) STEC F z VTEC  and z is the zenith angle at the ionospheric pierce point (IPP), which is the intersection between the GNSS signal path and the thin shell. STECs derived from GNSS measurements are thereby mapped as VTECs on the IPPs. In the next step, a 2D model of the VTEC distribution over the shell is created with a spherical harmonics expansion of the VTECs on the IPPs. The spherical harmonics expansion is achieved through where ϕ is the latitude, λ is the longitude, nm P  is the normalized associated Legendre function of degree n and order m, nm C  and nm S  are spherical harmonics coefficients to be determined.
Equation (8) is a slightly modified version of the one found in Schaer [1]. In this study, a geographic longitude was used instead of sun-fixed longitude, as each observation epoch was handled independently and no interpolation was done between epochs.

Computerized Iionospheric Tomography (CIT)
The 3D TEC distribution can be estimated with tomographic inversion. For an inversion volume defined regionally, this is often done by representing the ionospheric electron density (IED) distribution over a finite number of voxels, i.e., volume elements, with constant IED in each voxel. This can be represented as a linear observation model expressed in vector form as where l is a vector with M STEC observations along the ray paths, x is the vector of IEDs for each of the N voxels, and A is a matrix with distances that the signal travels through each of the voxels.
The problem of finding the IEDs in x is typically ill-posed for a network of receivers placed on the ground, and a unique least-squares solution to Equation (9) is therefore generally not possible to find. Several strategies that address this issue have been developed over the years. These strategies can be divided into iterative and noniterative inversion methods. The basic iterative methods include the algebraic reconstruction technique (ART) [4] and the multiplicative algebraic reconstruction technique (MART) [32]. Improvements of these, as well as new methods, have been suggested and applied in ionospheric tomography, for instance, the simultaneous iteration reconstruction technique (SIRT) [33], the decomposed algebraic reconstruction technique (DART) [2], and the improved algebraic reconstruction technique (IART) [34].
Noniterative inversion methods, including for instance singular value decomposition (SVD) [35] and truncated singular value decomposition (TSVD) [36], are worth mentioning as well. These methods achieve solutions in a least-squares-least-norm sense and share similarities with Tikhonov regularization. Generally, an ill-posed and unsolvable problem might be regularized to a solvable problem by adding additional cost terms to the objective function. In the case of Tikhonov regularization, an objective to minimize the norm is added to the problem formulation. In other noniterative inversion methods ,regularization is demonstrated by adding smoothness objectives, minimizing the difference of IEDs between adjacent voxels [37][38][39].
Successful tomographic inversion is dependent on the geometric distribution of the observations available. Figure 1 depicts the average number of signals that travelled through each voxel for an inversion volume defined over Sweden where observations from the three GNSSs GPS, Galileo, and BeiDou were employed. As can be observed, voxels on high latitudes and high altitudes generally suffer by a lower availability of GNSS observations. The relatively low availability of observations on higher latitudes means that the tomographic inversion will be more challenging than for regions close to the equator with more favorable GNSS satellite geometries. To achieve optimal usage of the available GNSS measurements, side rays were also utilized in this study. Side rays are rays that partially pass through the ionosphere outside of the inversion volume. As they therefore only partially coincide with the inversion volume, they are traditionally not employed in ionospheric tomography, since the part of the STEC originating for the inversion volume is unknown. According to the method described by Yao et al. [40], partial STECs (PSTECs) can be determined for the side rays by where GNSS STEC are the total STECs determined from GNSS measurements by Equation (2) and (3).
Even though additional information might be incorporated into the tomographic inversion by regularization, it is not possible to include inequality constraints with the closed form least squares adjustment formulas derived from the normal equations. One drawback is thus that, even with regularization, the solutions produced by such approaches might still, for instance, include negative IEDs.
In this study, we imposed smoothness and constrained the solution to only positive IEDs by applying convex optimization. This optimization technique finds the optimal value of a convex objective function that is constrained over a convex set of possible solutions. Optimality of the final solution is guaranteed by strong duality, which applies generally for convex optimization problems [42].
With this method, the inversion was done by solving the convex optimization problem where 1  and 2  are the L1 and L2 norms, respectively. B is a matrix for estimation of differences between adjacent voxels belonging to the same layer, i.e., for the row vector i for two adjacent voxels x . This constraint is only active above 350 km and prevents unreasonably large IED values for high altitudes. The C matrix identifies voxels in the top and bottom layer, and thus according to the problem, the formulation sets these to zero.
The convexity of Equation (11) follows from convexity of the L1 and L2 norms together with the convexity of linear constraints [42]. This inversion method was inspired by methods employing regularization achieving smoothness by adding the total variation (TV) norm as a cost term to the objective function [38,43]. A slightly simplified approach is to minimize the L1 norm of differences between adjacent voxels instead. Like with the TV norm, some edge preservation is achieved as minimization of the L1 norm tends to produce sparse solutions. This approach is beneficial, as the problem can be solved with a convex solver. In this study, this was done in Julia [44] with the embedded conic solver (ECOS) [45]. The suggested inversion method described in Equation (11) was evaluated against the NeQuick 2 model and some of the previously mentioned inversion methods, including ART, MART, SIRT, SVD, and regularization with minimization of IED differences between adjacent voxels. It proved to be superior regarding the challenging satellite geometries for the high latitudes involved in this study.

Alternative Method to Determine Nadir-Dependent Biases
The existence of nadir-dependent GDV was demonstrated in the investigation of signal anomalies associated with the GPS SVN 49 satellite [46]. Later studies also proved the existence of prominent nadir-dependent GDV of BeiDou second generation medium Earth orbiting (MEO) and inclined geosynchronous orbiting (IGSO) satellites [11]. Less prominent nadir-dependent GDV have also been demonstrated to exist for GPS observations [12].
Earlier methods have proved the existence of GDV by employing the MP linear combination of GNSS observations. By forming this combination, it is possible to isolate variations of a selected code observable, given that dual frequency carrier phase observations are available from the same satellite. The MP linear combination is formed by where 1 1 2 where P is a code observation, L is a phase observation in units of length, λ is the carrier wavelength, and B is phase ambiguities lumped together with hardware code and phase biases [29].
To separate the nadir-dependent part of the GDV from noise and local multipath, a spherical harmonics expansion of degree 5 and order 0 was used as regression model. This regression model for determination of nadir-dependent GDV was demonstrated by the authors of [12]. With this technique, high-frequency variations, probably consisting of noise and local high-frequency multipath, are filtered away, while low-frequency variations persist.

Data Processing
In this study, data processing was done according to the flowchart presented in Figure 2. First, DCBs were calibrated with GNSS observations from SWEPOS reference stations and a GIM from the IGS. This was done with Equation (2), as described in Montenbruck and Hauschild [47], with data covering a period of 24 hours. An elevation mask of 15 degrees was applied for the GNSS observation data. In the next step, STECs were estimated from calibrated DCBs and reference station GNSS observation data with Equations (2) and (3). This was done from the GNSS observations of about 400 SWEPOS reference stations, which are distributed over Sweden with a maximum distance between stations of 70 km. The estimated STECs were used as observations in both 2D and 3D ionosphere modeling as described by Equations (8) and (11). In the later analysis, residuals from the ionosphere modeling were compared with GDV derived from the MP linear combination described in Equation (14).

Handling of Remaining Error Sources
Local code multipath is a considerable error source that is not cancelled when forming the geometry-free linear combination. To mitigate this error source without removing the long-term variations of the code observable, the estimated code STECs were fitted to carrier phase STECs over 15-minute intervals. Of the fitted STECs, only the middle epoch in each interval was used in the subsequent ionosphere modeling. By doing this, a considerable part of the TEC estimation errors due to multipath and code noise were removed. Multipath and code noise amounted to a few decimeters at the SWEPOS sites, which might have employed choke ring antennas, and were located at places with favorable multipath conditions. After the fitting of code STECs to carrier phase STECS, all data processing was performed independently between epochs. This approach was chosen to avoid introducing extra time dependencies in the end results due to assumptions of the dynamics made by the estimation method.
Error sources with contributions of only a few centimeters were ignored in the data processing because uncertainties at this level did not affect the ionosphere modeling performed in this study, as the modeling errors themselves amounted to the decimeter level. The error sources ignored included phase wind-up [48], which only amounted to a few centimeters, and phase center offsets at the reference stations and at the satellites. For phase center offsets, only the relative offsets associated with different carrier frequencies remained for the geometry-free linear combination. As only variations over the nadir and elevation angle intervals were of interest in this study, this, in practice, also only amounted to a few centimeters. For the receiver antennas, the relative phase center offsets for the Dorne-Margolin antennas used within SWEPOS were a few centimeters, which was negligible. At the satellite side, a nadir angle interval of 0-14 degrees amounts to a variation of only 0.03 times the phase center offset [12]. In the case of the geometry-free linear combination, the relative phase center offsets were assumed to be at the meter level, which meant an effect of a few centimeters for the small nadir angle interval.

Results
Ionosphere modeling was performed in two and three dimensions from using three days of GNSS observation data, day of year (DOY) 180-182 in 2018. This is a period almost free from ionospheric disturbances, which means that the smoothness assumed in Equations (8) and (11) was more likely to be valid. STECs were estimated from the geometry-free linear combination according to Equations (2) and (3) from the GNSS signal combinations C1C-C2W, C6X-C1X, and C1X-C5X for GPS, BeiDou, and Galileo, respectively with signals denoted in accordance with the receiver independent exchange format RINEX 3 [49]. The residuals of the ionosphere modeling were analyzed and compared with GDV derived from the MP linear combination. As only regional ionosphere modeling was performed, only some of the satellites were observed over longer nadir angle intervals. It was thereby only possible to analyze the nadir angle dependence for a subset of all satellites belonging to the constellations of GPS, BeiDou, and Galileo.

Ionosphere Modeling in 2D
A 2D ionosphere model was estimated over Sweden with spherical harmonics expansion by Equation (8) and the assumption that the ionosphere was distributed over an infinitely thin layer according to the "thin shell" model. The height of the layer was set to 450 km, which was the same as for the GIMs used for the DCB calibration. Both the degree and order of the spherical expansion were conservatively set to 20 based on the reference station density of the network. Figures 3-5 present residuals for a sample of three satellites of 2D ionosphere modeling compared with GDV derived by the MP linear combination for GPS, BeiDou, and Galileo, respectively. Due to the estimation method, with regional ionosphere modeling, only a subset of all satellites could be processed over larger nadir intervals. The satellites with residuals available for the longest nadir angle intervals were selected, but the results are representable also for other visible satellites. In these figures, residuals are marked with blue stars and connected with line segments for a complete satellite pass, which includes angle dependencies both for the satellite rising and for the satellite setting. In these figures, the MP combination values themselves (orange dots), as well as the spherical harmonics expansion regression model values (green line), were shifted to the average level of the residuals. The residuals and the regression model values followed the same trend, especially for BeiDou-2 satellites, which were also expected to have larger GDV. The sizes of the BeiDou GDV, with a range of about 10 TECU, are in good agreement with previous studies by Hauschild, Montenbruck, Sleewaegen, Huisman, and Teunissen [10] and Wanninger and Beer [11]. For GPS and Galileo, the variations were considerably smaller, and it was thus harder to draw any conclusion as these were at, or below, the expected accuracy of the ionosphere model. The variations that can be observed were at a level that agrees with the findings of Wanninger, Sumaya, and Beer [12], which also confirmed that the GDV are less prominent for GPS. However, the receiver antenna also induced GDV of similar sizes, which could not be separated in the results.

Ionospheric Tomography
Ionospheric tomography with the inversion method described by Equation (11) was utilized to determine a 3D ionospheric model of IEDs over Sweden. The inversion volume with voxels and reference stations is depicted in Figure 1 and included all voxel columns coincident with the colored squares of the figure.
Figures 6-8 present the residuals from ionospheric tomography compared with GDV derived by the MP linear combination. Similar to the 2D case, the GDV for BeiDou, with a range of about 7 TECU, were considerably larger than the GDV for GPS and Galileo, but still not as large as the MP linear combination suggests. This might be an indication that some of the GDV were absorbed by the ionospheric model. The absorption of GDV by the model could be caused by the correlation between the GDV and the model parameters. This was more likely for the tomographic inversion problem, as this model was ill-conditioned and had considerably more degrees of freedom than the model in the 2D case. Similar to the 2D case, the GDV were small for GPS and Galileo, and the values fell below the expected accuracy of the ionosphere model.

Discussion
The results show large nadir-dependent GDV for the BeiDou-2 satellite for both 2D and 3D ionosphere modeling, making the results consistent for both types of ionosphere modeling even though they fundamentally differed in the way they modeled the ionosphere. They are also consistent in size with the results of previous studies [10,11]. Furthermore, the shown BeiDou-2 variations of the residuals, of about 10 TECU for 2D ionosphere modeling and about 7 TECU for 3D ionosphere modeling, were too large to be caused by remaining possible error sources, which included TEC estimation errors due to, for instance, low-frequency multipath, modeling errors, unmodeled phase wind-up and phase center offsets, and GDV of the receiver antenna. These error sources did not contribute more than on the decimeter level. Using the same reasoning, nadir-dependent GDV cannot be inferred for GPS and Galileo by analyzing residuals from the two types of ionosphere modeling, as the variations for these systems lie below the expected sizes of remaining not considered error sources.

Conclusions
This paper describes developments and testing of a new approach for estimation, as well as a characterization of the nadir-dependent GDV. This is done by analyzing the residuals from ionosphere modeling based on the geometry-free linear combination of GNSS code and phase observations. The approach was investigated for ionosphere modeling both in 2D, with the "thin shell" model, as well as 3D with tomographic inversion.
Nadir-dependent GDV are known from previous studies to be present in the code observable, especially for BeiDou-2 IGSO and MEO satellites, but also for GPS to a lesser degree. It was verified that these GDV indeed were present and affects both 2D and 3D ionosphere modeling by spherical harmonics expansion of VTECs and ionospheric tomography. Investigating the residuals of ionosphere modeling shows that these effects are especially prominent for BeiDou.
The values of GDV estimated for GPS and Galileo could not be verified due to the limited accuracy of the ionosphere modeling. It was also concluded that these GDV were probably absorbed by the ionospheric model in the 3D case with ionospheric tomography to a greater extent. This might be due to the higher number of degrees of freedom in the tomographic problem formulation, and it suggests that ionospheric tomography is more sensitive to GDV than the simpler 2D VTEC modeling. This result also implies that 2D ionosphere modeling with the "thin shell" model might be better suited for GDV characterization than for 3D ionosphere modeling.
It is important to be able to characterize GDV because, as earlier studies demonstrated, they affect PPP, integer ambiguity resolution by the Hatch-Melbourne-Wübbena linear combination, as well as TEC estimation. The latter is important for various precise positioning applications where the ionospheric influence explicitly must be considered.