Next Article in Journal
Multibeam Characteristics of a Negative Refractive Index Shaped Lens
Next Article in Special Issue
A Modified TurboEdit Cycle-Slip Detection and Correction Method for Dual-Frequency Smartphone GNSS Observation
Previous Article in Journal
Defect-Induced Gas-Sensing Properties of a Flexible SnS Sensor under UV Illumination at Room Temperature
Previous Article in Special Issue
Impact of Using GPS L2 Receiver Antenna Corrections for the Galileo E5a Frequency on Position Estimates
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

GNSS-Based Non-Negative Absolute Ionosphere Total Electron Content, its Spatial Gradients, Time Derivatives and Differential Code Biases: Bounded-Variable Least-Squares and Taylor Series

Institute of Solar-Terrestrial Physics SB RAS, Irkutsk 664033, Russia
*
Author to whom correspondence should be addressed.
Sensors 2020, 20(19), 5702; https://doi.org/10.3390/s20195702
Submission received: 31 August 2020 / Revised: 2 October 2020 / Accepted: 4 October 2020 / Published: 7 October 2020
(This article belongs to the Special Issue GNSS Signals and Sensors)

Abstract

:
Global navigation satellite systems (GNSS) allow estimating total electron content (TEC). However, it is still a problem to calculate absolute ionosphere parameters from GNSS data: negative TEC values could appear, and most of existing algorithms does not enable to estimate TEC spatial gradients and TEC time derivatives. We developed an algorithm to recover the absolute non-negative vertical and slant TEC, its derivatives and its gradients, as well as the GNSS equipment differential code biases (DCBs) by using the Taylor series expansion and bounded-variable least-squares. We termed this algorithm TuRBOTEC. Bounded-variable least-squares fitting ensures non-negative values of both slant TEC and vertical TEC. The second order Taylor series expansion could provide a relevant TEC spatial gradients and TEC time derivatives. The technique validation was performed by using independent experimental data over 2014 and the IRI-2012 and IRI-plas models. As a TEC source we used Madrigal maps, CODE (the Center for Orbit Determination in Europe) global ionosphere maps (GIM), the IONOLAB software, and the SEEMALA-TEC software developed by Dr. Seemala. For the Asian mid-latitudes TuRBOTEC results agree with the GIM and IONOLAB data (root-mean-square was < 3 TECU), but they disagree with the SEEMALA-TEC and Madrigal data (root-mean-square was >10 TECU). About 9% of vertical TECs from the TuRBOTEC estimates exceed (by more than 1 TECU) those from the same algorithm but without constraints. The analysis of TEC spatial gradients showed that as far as 10–15° on latitude, TEC estimation error exceeds 10 TECU. Longitudinal gradients produce smaller error for the same distance. Experimental GLObal Navigation Satellite System (GLONASS) DCB from TuRBOTEC and CODE peaked 15 TECU difference, while GPS DCB agrees. Slant TEC series indicate that the TuRBOTEC data for GLONASS are physically more plausible.

Graphical Abstract

1. Introduction

Currently, the global navigation satellite systems (GNSS) enable the study of the ionosphere at any spot on the globe. Such studies are based on dual-frequency phase and pseudorange measurements of the total electron content (TEC) [1]. Since the phase measurements draw a lower level of noise compared to pseudorange measurements, most studies use the phase measurements to determine TEC. There are classical papers [2,3], as well as many up-to-date ones (for example, see [4] and references therein). Lanyi and Roth [2], for the first time, suggested mapping the TEC. They compared, mapped and measured total ionospheric electron content by using global positioning system (GPS) and beacon satellite observation and revealed that the difference between them is less than 1 TECU. Pioneering papers by Calais and Minster [3] showed the efficiency of dual-frequency phase measurements to investigate ionospheric responses to the strong ground displacement associated with earthquakes. Afraimovich et al. [4] presented a number of new results when studying the ionospheric irregularities caused by solar flares, solar eclipses, solar terminators, earthquakes, magnetic storms, rocket launches and tropical cyclones based on GPS/GLObal Navigation Satellite System (GLONASS) data. The studies showed the efficiency of dual-frequency phase measurements. А comprehensive study by Astafyeva et al. [5] showed a high potential of TEC phase measurements to monitor natural hazards. Nesterov and Kunitsyn [6] implemented 4D GNSS radio tomography based on phase measurements. Kunitsyn et al. [7] suggested a new approach to use L1/L5 phase measurements from satellite-based augmentation system. However, when using phase data, there is an ambiguity of phase measurements. For some applied problems, we need not relative, but absolute measurements. Forte and Aquino [8] evaluated the ionospheric effects on low frequency radio astronomy measurements. Afraimovich and Yasukevich [9] suggested an approach based on absolute TEC data for ionospheric correction of radio telescopes and radio interferometers. Ovodenko et al. [10] used absolute GNSS TEC data to correct radar data.
Pseudorange measurements are thought to be absolute, however, they are noisy. For this reason, one often uses pseudorange and phase measurements jointly to eliminate the phase ambiguity. Herewith, there is a bias related to a different time of the signal passage through the channels of a satellite and a receiver. In literature, this error is referred to as differential code biases (DCBs). This error is known to vary systematically, and it may have irregular variations [11]. Mylnikova et al. [12] showed some cases of DCB annual variations caused by receiver/antenna environment conditions. Instability in DCB can be caused by change in grounding [13]. А change in DCB can lead to significant apparent TEC variations [14]. There are some procedures to estimate the DCBs [15,16]. Hong et al. [15] suggested an efficient algorithm to estimate DCB for a network, when one of receiver DCBs is already known. Jin et al. [16] published MatLab code for DCB calculation. They reported an agreement with the global ionospheric map (GIM) data with a mean difference of less than 0.7 ns (~2 TECU) and an RMS of less than 0.4 ns (~1 TECU). Li et al. [17] analyzed triple-frequency combinations and proved that the real satellite DCBs are time-varying. Wang et al. [18] mitigated variations in DCB by separating DCB from ambiguity. As a rule, the DCB estimates are a by-product of estimating the model parameters of TEC measurements. Having the DCB estimated, one can correct the slant TEC series and then use these TEC absolute values in applications.
An important goal is to obtain absolute characteristics reflecting ionospheric conditions. A good way is to have an electron density profile that an incoherent scatter radar or an ionosonde can provide [19,20,21,22]. The cost of such equipment is very high, particularly that of radars. Therefore, the world network of such facilities is limited.
Using the GNSS data can also provide estimates for the electron density profile. This information may be obtained by using 4D spatio-temporal ionospheric tomography suggested by Mitchell and Spencer [23]. GNSS-tomography results agree with ionosphere models and ionosonde data [24]. The tomography algorithms have been developed over the last decades. Additional data (such as ionosonde and Langmuir probes [25]) and new reconstruction algorithms [26] enhance GNSS-tomography. However, these techniques require complicated calculation and large computation resources. Also, for precise measurements, it is necessary to have a very dense network, as is available in Japan or in the United States.
When using the data from sparse networks of stations, a more promising way is to build TEC maps. The International GNSS Service (IGS) Working Group on Ionosphere [27] provides the most common TEC maps. Such global TEC maps are used to analyze the regional or global electron content [28]. For regions with a small number of stations, it makes sense to estimate the vertical TEC over the station.
Currently, there are several methods to estimate the vertical TEC. These methods are based on various TEC models [29,30]. One usually uses the spherical harmonic expansion to build global vertical TEC maps. Estimation of the vertical TEC over the station by using spherical harmonics is not appropriate, it is better to use other expansions. For example, Durmas and Karslioglu [31] use a TEC model based on expansion of basic functions that are the products of univariate B-splines with different scales and variables. There are also methods based on simple polynomial expansion. For the first time, this method was offered by Lanyi and Roth [2] and by Sardón and Zarraoa [32], and then was used by Themens et al. [33], who studied the nature of seasonal and solar cycle bias variability in the polar cap region and showed the influence of user’s choice of shell-height on standard DCB estimation procedure. Also, the Taylor series expansion (as a generalized expansion) of TEC is used as a TEC model; this expansion on space was mentioned by Schaer [34] for regional vertical TEC modeling. We should note that Li et al. [35] mentioned a problem with unstable DCB at some BeiDou satellites and suggested two-step technique to reduce effects of such satellites on overall solution.
Different techniques for TEC estimation do not solve the problem of negative TEC values, especially negative slant TEC on some line of sights. Start and Parker [36] and Waterman [37] solved the mathematical problem for constrained (or bounded-variable) least squares. Zhang et al. [38] used such an approach to improve GIMs. We developed an algorithm to recover the absolute vertical TEC, its gradients, and its time derivatives from a single GNSS station data, as well as slant TEC along all lines of sight and DCB. The procedure uses the Taylor space-and-time expansion and constrained least squares. It provides both non-negative vertical TEC and non-negative slant TEC along all lines of sight. The algorithm can use data from MEO GPS/GLONASS/Galileo/BeiDou satellites. A combined solution can also use the GEO satellite data. However, here, we show only the GPS and GLONASS data. We termed the algorithm TuRBOTEC: TayloR-series and bounded-variable-least-squares-based iOnosphere TEC.

2. Data and Background

We use the GNSS data from the IGS network [39] as the TuRBOTEC input data. We use the data from three stations: IRKJ, NTUS, and THU2. Table 1 shows the stations’ coordinates: geographic latitude (Lat), geographic longitude (Lon), geomagnetic latitude (MLat), geomagnetic longitude (MLon). The coefficient α is a parameter for the modified single-layer model mapping function (see below). To compare with the obtained results, we used data from the IONOLAB software (http://www.ionolab.org/) [40] and software by Seemala [41], as well as data from ionospheric maps: Madrigal TEC [42] and CODE GIM [34]. It is free-to-use software and data. Figure 1 shows TEC from different techniques for the IRKJ region (52° N, 104° E). The upper panels display the IONOLAB-TEC (left, red line for vertical TEC) and SEEMALA-TEC (right, blue line for vertical TEC). The left bottom panel shows Madrigal TEC (orange line for vertical TEC). The right bottom panel shows TEC from all the above methods and TEC from CODE GIM.
The IONOLAB software is based on a single-station solution [40]. The software uses only GPS data. The input is a pseudorange slant TEC corrected by CODE DCB. Slant TEC is projected to the local zenith to obtain the vertical TEC. Actually it is a zero order expansion. The least-squares method is used to obtain TEC estimates. Our experience showed the problem when negative or zero values appear: negative TEC values result in rejecting data for the whole day.
Seemala developed a software for TEC calculations (referred to as SEEMALA-TEC) [41] based on a similar approach as that in the IONOLAB. However, it uses phase TEC leveled by pseudorange TEC. The vertical TEC is estimated as a mean from all satellites’ vertical TECs. The software can estimate its own DCB, but our experience showed that it is better to use CODE DCB.
Madrigal TEC is also based on projection of slant TEC to vertical TEC, but for different cells (1° × 1°) on the map [42]. They use a novel 3-steps technique to estimate DCB. Averaging data in cells improves the data. However, in the regions, where there are few stations, it is difficult to have continuous and reliable data. Thus, we use a 5° × 5° region around a GNSS station and choose the median value as a Madrigal vertical TEC. Spline is used to smooth median TEC series.
Global ionosphere maps are used for a long time to monitor and model the ionosphere [27]. We used CODE GIMs [30] that are based on spherical harmonics and modified single layer mapping function. We also used JPL GIMs based on multi-layer model. The GIM time resolution was 2 h for the analyzed data. We performed an algorithm operation simulation based on the international reference ionosphere model IRI-2012 [43] and the IRI-plas model [44] as a test.

3. TuRBOTEC Algorithm

To obtain absolute TEC, one estimates the preset measurement model parameters by using experimental data. As a rule, one uses the following model of TEC:
IM = IS + IDCB = IV∙S + IDCB
where IM is the model (expected) value of the recorded slant TEC, Is is the slant TEC along the line of sight without a differential code bias, IV is the vertical TEC estimate at a point corresponds with the measurement (ionospheric pierce point), S is the function that converts the slant TEC into the vertical one (mapping function), and IDCB is an error related to the DCB of the satellite and the receiver.
For geostationary Earth orbit (GEO) satellites it is difficult to estimate DCB (or similar error in single frequency data) from (1), because S in Equation (1) is a constant. The corresponding error can be regarded as a TEC offset [45]. However, GEO satellites’ DCB could be estimated along with DCB of medium Earth orbit (MEO) satellites [46].
Ma and Maruyama [47] used de facto the approximation IV=IV0, l0, t0), where ϕ0, l0, and t0 are the station latitude, longitude, and the time, for which the calculation is performed. This implies that all the TEC measurements were supposed to be performed over the station. However, such an approach does not allow one to obtain spatial gradients. Such an approximation is correct in a limited number of cases when the spatial gradients and the time derivative can be neglected. In the most studies, a spatial-temporal expansion of the IV function is used. The spherical harmonic expansion is often used. Komjathy et al. [48] increased the capacity of GIM technique to simultaneously treat data from more than 1000 dual-frequency GPS receivers. Zhang et al. [49] used the DCB estimation technique based on spherical expansion to evaluate BeiDou DCBs. According to Schaer [34], the Taylor expansion is more appropriate than the spherical harmonic expansion for regional modeling. Schaer [34] mentioned the IV Taylor series expansion in (1) only on space. In our research, we expand the vertical TEC function IV (ϕ, l, t) into the Taylor series on space and time:
I V ( ϕ , l , t ) = m = 0 n = 0 k = 0 ( Δ ϕ ) m ( Δ l ) n ( Δ t ) k m ! n ! k ! m + n + k I V ϕ m l n t k
where Δϕ, Δl are the difference between latitude, longitude of the ionospheric pierce point and the point of station (ϕ0, l0), Δt is the difference between time of measurement and t0.
We developed the following algorithm to estimate the vertical TEC, the TEC gradients, the time derivative, and DCBs. The algorithm involves:
(1)
Calculating TEC based on the pseudorange IP and phase Iφ measurements. For the analysis, we use the data with elevations greater than 10°.
(2)
Dividing the data into continuous samples.
(3)
Detecting and eliminating outliers and cycle slips in the TEC data [50] (Figure 2).
(4)
Eliminating the phase measurement ambiguity (“leveling”, see Figure 2a):
c o n s t = 1 N i = 1 N S i 1 · i = 1 N ( I P I φ ) S i , where N is the number of measurements over a continuous interval, S is the mapping function (see below). At this state we obtain experimental slant TEC IExp.
(5)
Estimating DCBs by a simple measurement model and determining the model parameters based on minimizing the model data root-mean-square deviation.
Although in some research the zero-order expansion of (2) is used (only term with m = 0, n = 0, k = 0 in Equation (2)), indeed, it is necessary to substantiate the selection of the expansion order in greater detail. We performed an analysis by using the IRI-2012 model. The TEC slant values, corresponding to real observed elevations and azimuths, as well as the vertical TEC, the time derivative, and the spatial gradients were modeled by latitude and by longitude. Then, we recovered the values of vertical TEC (IV) from slant measurements by using expansions to the specified order of Taylor-expansion.
At the first step, this recovery was performed based on the zero-order Taylor expansion (2) (leaving term with m = 0, n = 0, k = 0 in (2)). Next, TEC was recovered by using the first-order (leaving terms with n + m + k ≤ 1 in (2)) and second-order (leaving terms with n + m + k ≤ 2 in (2)) expansions.
Figure 3 provides the results for recovery of the vertical TEC, the time derivative, and the latitudinal gradient. The solid black line denotes the model values on the plots. The recovered parameters obtained based on of the zero order are denoted by the black dash-dotted line, those based on the first order―by the blue dotted line, and those based on the second order―by the red solid line. From the plots, one can see that using the second order of the IV expansion into the Taylor series, viz., using of time derivatives and square gradients on space, considerably improves the absolute vertical TEC recovery accuracy.
Meanwhile, using the second order provides a sufficient accuracy to be able to neglect higher orders. Therefore, we use the expansion IV into the Taylor series up to the second order:
I M = S j i [ I V ( ϕ 0 , l 0 , t 0 ) + G ϕ Δ ϕ j i + G q _ ϕ ( Δ ϕ j i ) 2 + + G l Δ l j i + G q _ l ( Δ l j i ) 2 + G t Δ t j i + G q _ t ( Δ t j i ) 2 ] + I D C B , j ,
where IV is the absolute vertical TEC value, Δϕ (Δl) is the latitude (longitude) difference between the ionospheric point coordinate ϕ(l) and that of the station ϕ0(l0), Δt is the difference between the measurement time t and the time t0 for which the calculation is performed. Further, Gϕ = ∂IV/∂ϕ, Gl = ∂IV/∂l, Gq_ϕ = ∂2IV/∂ϕ2, and Gq_l = ∂2IV/∂l2 are the linear and quadratic spatial TEC gradients, and Gt = ∂IV/∂t and Gq_t = ∂2IV/∂t2 are the first and second time derivatives. Equation (3) represents the second-order Taylor series expsion of (1). We neglect the mixed derivatives.
In the literature different mapping functions are used. Selecting the most correct one is a challenge. In most papers, one uses the one-layer approximation of the ionosphere as a thin layer. This approach was suggested by Klobuchar [51]. Hernández-Pajares et al. [52] used the two-layer approximation of the ionosphere. Based on a tomographic procedure, Lyu et al. [53] proposed a climatological mapping function―the Barcelona Ionospheric Mapping Function. There is another approach: Schüler and Oladipo [54] suggested one use the NeQuick model instead of the thin layer approximation to convert the slant TEC into the vertical one. However, such an approach did not provide an essential accuracy increase [54]. We use the modified single-layer model mapping function, and coefficient α is introduced to more correctly account for the ionospheric layer height [34]:
S j i = [ c o s { arcsin ( R E R E + h m a x · sin [ α ( 90 θ j i ) ] ) } ] 1 ,
where RE = 6371 km is the Earth radius, hmax = 450 km is the ionospheric point height, and α is the correcting coefficient.
Figure 4 presents the slant TEC dependence, IS, on the elevation when converting IS = S∙IV for different values α. Simulation used the IRI-2012. By its physical implication, the coefficient α adjusts a sharp non-physical TEC growth at low elevations, less than 30°.
It is worth noting that due to the ionosphere peculiarities, the coefficient α should be latitude-dependent. Also, the α parameter should affect the ionosphere disturbance level that may dramatically vary the ionization global distribution [55]. We estimated the α values for several points on the earth surface. Table 1 presents the results.
We obtain the set of equations by minimizing the functional U = ∑Uk (5) for the set of selected instant tk, for which we estimate the parameters. For computation, we represent (5) in the form (7):
U k = j = 1 N k i = 1 N j k ω k , j i ( I M j i I E x p j i ) 2 ,
ω k , j i ω k ( t j i ) = Θ ( t k t j i + Δ t ) Θ ( t j i + Δ t t k ) 1 S j i [ 1 + ( Δ t j i , k 1 h o u r ) 2 ] 1 ,
U = ∑Uk = ||Ax –b ||2 -> min,
where IExp is the experimental phase TEC measurements with the eliminated phase ambiguity obtained after Stage 4 of the algorithm, Θ is the Heaviside step function, t j i is the ith instant of measurement for the jth satellite, Δ t j i , k is the time difference between the current measurement t j i and the time tk, for which calculating is performed, and Δt = 1 h is the maximal time difference, for which the data are still used for analysis (Figure 5 shows which measurements affect the tk-estimate). The 1/S factor in Equation (6) causes the measurements at high elevations to produce the greatest contribution. For each time instant k, we have 7∙J + N variables (IkV, Gk ϕ, Gkl, Gkq _ ϕ, Gkq_l, Gkt, Gkq_t, and IDCB, j) in Equation (5), where J is the number of instants over the investigated interval, for which calculation is performed, and N is the number of the satellites observed.
After that, minimization should be applied, which is based on the least square technique. A typical problem for TEC estimation is emerging negative or zero values. For GIM data, this leads to zero values in some GIM cells. Actually, the solution for such problems was suggested a long time ago. We need to restrict the estimated values [36,37]. First results to apply restriction was used to improve GIMs [38]. For our task, we need to obtain both non-negative vertical TEC and non-negative slant TEC. The slant TEC after DCB removal should at least exceed some value C. So we introduce the next boundaries:
IV(t0) > C, ∀t0
(IDCB)j < (IExp)min,jC, ∀ satellite j
where C is a non-negative value of minimal TEC, which can be observed in principle. We chose C = 0.5 TECU.
We used the Python library scipy.optimize.lsq_linear based on algorithm suggested by Start and Parker [36]. The main steps of the algorithm are as follows:
(1)
The algorithm first computes the usual least-squares solution. This solution is returned as optimal, if it lies within the bounds. If not, the algorithm finds all variables within the bounds (free set) and beyond (active set).
(2)
At each iteration the algorithm chooses a new variable (which has maximal gradient of the squared objective) to move from the active set to the free set.
(3)
New equation system for free set is created where b in (7) is changed by active set. Least-squares solution for new equation system contains variables beyond the bounds, the gradient correction is applied to all the free set (see [36] for details).
(4)
The iterations continue until all the variables are in the free set.
This algorithm ensures an accurate solution eventually, but may require about n iterations for a problem with n variables. As a result we obtain non-negative (positive) vertical TEC and slant TEC at all the lines of sight. To obtain robust DCB estimates IDCB, we calculate the parameters simultaneously for different time instants over 24 h, thus solving a consistent set of Equations (7). The temporal resolution for tk may vary from 2 h to 5 min.
To analyze the satellite DCB separately, we apply the often-used zero-mean condition [34]:
i = 1 N I D C B i = 0 ,  
where N is the number of satellites in the constellation. We applied condition (9) for GPS and GLONASS (Galileo and BeiDou) separately, following Schaer [56].
To simulate the algorithm operation, we used the IRI-2012 model with a set of the International Union of Radio Science (URSI) coefficients, recommended by the URSI and the IRIcorr topside ionosphere profile option [43]. To check the influence of the plasmasphere, we performed calculations based on the IRI-plas model [44]. Currently IRI-plas is a standard for the plasmasphere calculation.
Simulation was performed for a selected real station. For each recorded satellite at each time instant, we calculated the electron density and TEC along the line of sight. Further, we introduced the DCB-related error, as well as random noise to the phase TEC and the pseudorange TEC. This value was 0.01 TECU for the phase TEC. For the pseudorange TEC, the noise value was assigned depending on the elevation: five TECU at θ > 60°; 5 + 1·S10(θ + 30) TECU at θ < 60°. Also, we introduced outliers for the pseudorange and losses of phase lock for the phase, which could occur with a probability. Thus, at the output, we obtained a series of the phase TEC and the pseudorange TEC corresponding to a certain time instant, elevation, and azimuth. Further, we used these values and calculated ionospheric parameters based on the above algorithm. Finally, the parameters obtained as a result of the algorithm operation were compared with the modeled vertical TEC, time derivatives, and gradients. The obtained DCBs were compared with those specified as the simulation input.

4. Technique Validation and Discussion

We validated data based on ionosphere modeling and TEC products from alternative software mentioned in Section 2.

4.1. Absolute Total Electron Content

Figure 6 presents the results of the IRI-2012 simulating the TuRBOTEC algorithm operation (black dotted lines) for different stations: the mid-latitude IRKJ, equatorial NTUS, high-latitude THU2. The results were obtained for the 10 April 2012. One can see good qualitative and quantitative agreement of the vertical TEC recovery results with the IRI-2012 input data (red line). The average deviation for IRKJ by the absolute value is 0.1 TECU, maximal is 0.3 TECU, and RMS being 0.09 TECU. For NTUS, the average deviation is 0.4 TECU, maximal 0.95 TECU, RMS being 0.35 TECU. For THU2, the average deviation is 0.23 TECU, maximal −0.6 TECU, RMS being 0.17 TECU.
The differences between TuRBOTEC simulation and the IRI-plas vertical TEC are shown in green. The results show that TuRBOTEC provides even better performance for modeling based on IRI-plas at mid-latitudes, and comparable performance at low latitudes. Higher difference at high latitudes can be due to using α obtained through the IRI-2012 modeling, so we used α=0.96 for THU2. We did not find significant influence of the plasmosphere on TuRBOTEC estimations, except high latitudes, where we should to use another α in (4).
The greatest deviation occurs in the equatorial region. This is related to a substantially inhomogeneous ionosphere structure in this region, particularly, during daylight hours. The error is less at high latitudes. Although, as compared with the mid-latitude region, the error is higher. This is related to the absence of satellite observations at high elevations in high latitude regions, and to the dominance of the southward contribution to the total measurement statistics.
Figure 7 presents the results of estimating the vertical TEC for IRKJ by real dual-frequency measurements from GPS and GLONASS (red line). The dot-and-dash blue line and the dotted black line show the data on GIMs from the JPL and CODE laboratories, respectively. The data are presented for the 5 March 2015 quiet conditions (left panel, Kpmax = 2.3, maximal Kp for the day) and for the 17 March 2015 magnetic storm (right panel, Kpmax = 7.7).
One can see well that the TEC curves reproduce the diurnal variation similarly, but they can quantitatively differ. Such systematic differences are well-known and were repeatedly addressed in literature [57]. Like discussed earlier, the JPL data, as a rule, surpass the data from other laboratories.
One can see that, at individual instants, synchronous variations with TEC are noted in the CODE and TuRBOTEC data. For example, for the 17 March 2015 strong magnetic storm, one can see a slight increase in the vertical TEC around 18 UT. In general, the estimates for the vertical TEC appear plausible for all the considered cases. Deviations in the data from other laboratories are within the variance interval of different laboratories among themselves. To analyze systematic variances, we built a histogram for the ΔIV difference distribution between the vertical TEC data in the Irkutsk region.
Figure 8 shows distribution of differences between the vertical TEC from alternative and TuRBOTEC data. Data for IRKJ station over 2014. Panel (a) shows GIM CODE vs. TuRBOTEC; (b) Madrigal vs. TuRBOTEC; (c) IONOLAB vs. TuRBOTEC; (d) SEEMALA-TEC vs. TuRBOTEC. Panel (e) shows the distribution between the difference in JPL GIM TEC and CODE GIM TEC. (JPL vs. CODE)-distribution (e) features similar against (CODE vs. TuRBOTEC) or (IONOLAB vs. TuRBOTEC) root-mean-squares, that are 1.5 TECU vs. 1.7 or 2.1 TECU, and similar mean - ~2.7 TECU vs. 2.3 or 2.1 TECU. The Madrigal data and the SEEMALA-TEC data show a high discrepancy with the TuRBOTEC: the root-mean-squares were 10–12 TECU. Panel (f) compares the techniques with (TuRBOTEC) and without constraint. It shows that most vertical TECs are the same with and without constraint. However, about 9% of vertical TECs from the TuRBOTEC estimates exceed (by more than 1 TECU) those from the same algorithm but without constraints.
The vertical TEC from the previous version of the suggested technique (without constrains and several other issues) was experimentally checked in [58,59]. Those results also show relevant vertical TEC dynamics during space weather events.

4.2. Spatial Gradients. Accuracy of Determining TEC at a Growing Distance from a Station

Figure 9 presents the TEC latitude and longitude gradients for IRKJ, NTUS, and THU2. The data were obtained as a result of simulation based on the IRI-2012. The initial measurements are the same, as those for Figure 6. The parameters presented in the figure are a solution for a consistent set of Equations (7). The black dotted line marks the results of the TuRBOTEC operation; the red line shows the vertical TEC data and gradients from IRI-2012 model.
Absolute gradients differ from <0.1 TECU/deg. for longitude gradients in the high-latitude regions to 2.5 TECU/deg. for latitude gradients in the equatorial regions during the equatorial anomaly evolution.
The IRI-2012 gradients and their TuRBOTEC estimates are qualitatively similar. However, as compared with the vertical TEC recovery accuracy, the accuracy of determining the spatial gradients is substantially lower. Notably, for equatorial stations, the divergence in the latitude gradient may reach 1 TECU/deg. during the equatorial anomaly evolution. This is related to that model (3) involves S function that converts the slant TEC into the vertical TEC. In the equatorial anomaly region, such an approximation works worst of all, and, whereas the vertical TEC estimate is quite reasonable (see Figure 6), the spatial gradient estimate may not be satisfactory.
When estimating the vertical TEC at a growing distance from a station, this leads to an error. In Figure 10, we present the value of this error at different UT for various distances from a station. The figure presents the simulation data for the mid-latitude IRKJ station. Panel (a) shows the error when moving away from the station along the latitude; panel (b) does the same when moving away from the station along the longitude. The arrowlets mark the local midday and the solar terminator time at 200 km: SR is the morning (sunrise) terminator, SS being the night (sunset) one. On panel (b), the time of local midday and the terminator is marked for the longitude extremities: for 124.3° E (arrowlets at the plot top) and for 84.3° E (arrowlets at the plot bottom).
One can see that the latitude errors grow more than when moving at a similar distance by longitude. This may be determined by that the S mapping function has an essential latitude dependence. Therefore, at a significant deviation, the vertical TEC value converted from the slant TEC is not precisely determined. Hence, the gradients are not precisely determined, either. It is worth noting that, in this case, this error would grow almost linearly. From Figure 10, one can see that this is not absolutely so. The real and the model latitude gradients have a dramatically spatially inhomogeneous character as one moves away from the station, and, at distances more than 10°–15°, the estimates start to surpass 10 TECU.

4.3. TEC Time Derivative

Figure 11 shows the results of simulating the algorithm operation to recover the TEC time derivative. We obtained TEC time derivatives when solving the consistent set of equations, whose results are presented in Figure 6 and Figure 9. The designations are the same as those in Figure 6. Similar to the results for the absolute vertical TEC (Figure 6), there is agreement of both general dynamics and the quantitative values of the TEC time derivative. This indicates that such data may be used to efficiently predict the vertical TEC in the station region.
Figure 12 presents the results of the TEC time derivative experimental recovery through the TuRBOTEC algorithm from the IRKJ data (red solid curve) and the GIM-obtained time derivative from CODE (black dotted line), from JPL (blue dot-and-dash line) for 5 March 2015, Kpmax = 2.3 (left panel), and for 17 March 2015, Kpmax = 7.7 (right panel). Like for the absolute vertical TEC data (Figure 7), there is qualitative and quantitative agreement with the GIM data, the divergence values from other laboratories are within the divergence interval of different laboratories among themselves.

4.4. Differential Code Biases and Absolute Slant TEC

Figure 13 presents the modeled results of the DCB (red circles) recovery. The data are in TECU. A random number generator assigned the initial biases (black triangles). Then, there was an IRI-2012 simulation and determining DCBs by a set of integral TEC data through the above algorithm. For simulation, we used the geometry of a real station, IRKJ, as of 10 April 2012.
One can see a good performance of the algorithm both for the GPS data and for the GLONASS data. The maximal error is 0.6 TECU, RMS being 0.3 TECU. The obtained estimates are sufficiently robust both for quiet and for disturbed conditions.
Figure 14 presents the DCB experimental estimates (in TECU, red circles) obtained from real GPS and GLONASS measurements through the algorithm operation for all the GPS (a) and GLONASS (b) satellites. The estimates were obtained from the IRKJ 1 March 2015 measurements. For comparison, we show the CODE estimates (black diamonds). For the GPS DCBs, the divergence of estimates is sufficiently small. The maximal error for the GPS satellites is 5 TECU, RMS being 2.4 TECU. For the GLONASS satellites, this error is higher. The GLONASS maximal error is 17 TECU, RMS being 8.8 TECU. The results for GPS quite correlate with the estimates by Jin et al. [16]. At the same time, we observe a significant deviation of the CODE and TuRBOTEC estimates for GLONASS. There is no such difference between GPS and GLONASS in simulation. Slant TEC demonstrates non-negative values for GLONASS, when we use TuRBOTEC DCB instead of CODE DCB for GLONASS.
Figure 15 shows absolute slant TEC for both techniques with constraint (TuRBOTEC, blue dots) and without (black dots). The data are for 5 January 2014. We show the data for all satellites in view. After least squares without constraint, we note negative TEC values after 12 UT for several satellites. Constraint provides non-negative values for all the satellites.

4.5. Influence of Intra-Day DCB Variations

Nie et al. [60] revealed that the intra-day DCB variations result in a mis-modeling error of several tenths of TECU. Grounding change could lead to such variations [13]. Figure 16a shows an example of pseudorange TEC that features by rapid step-like changes (jumps). Such huge errors destroy the solution of equation system based on (1). Black line in Figure 16b shows the obtained vertical TEC dynamics without typical daily TEC variation.
We can solve the emerged problem by excluding DCB from Equation (1):
IM = IS + IDCB = IV∙S + Iconst,
where, Iconst is a term corresponding to the phase ambiguity. The term is constant for a continuous TEC series. We change (3) in the same way. Thus, for Equation (7), we exchange N variables (where N is a number of satellites) by N’ those (where N’ is a number of continuous series). This improves the solution (see red line in Figure 16b), but does not provide DCB anymore. It is the users, who should choose, if they need DCB, and if the pseudorange data quality is sufficient. We would note additionally, that short series should be removed from the treated data.

5. Conclusions

We have developed an algorithm to recover the absolute TEC, its gradients, its time derivative, and DCBs. The procedure is based on the space-and-time Taylor expansion and bounded-variable least squares. We termed it TayloR-series and Bounded-variable-least-squares based iOnosphere TEC (TuRBOTEC). We simulated the algorithm operation by using the IRI-2012 and the IRI-plas models. The absolute TEC values recovered through the developed algorithm were established to agree (qualitatively and quantitatively) with the IRI-2012 model-set values. The mean standard deviation for the mid-latitude IRKJ station is 0.09 TECU, for the equatorial NTUS is 0.35 TECU, for the high-latitude THU2 is 0.17 TECU. We did not find significant influence of the plasmosphere on TuRBOTEC estimations, except for high latitudes, where we should use another α for the mapping function. About 9% of experimental vertical TECs from the TuRBOTEC estimates exceed (by more than one TECU) those from the same algorithm but without constraints.
The recovered values of the TEC spatial gradients and of the TEC time derivative agree qualitatively with the model-set values. Also, we studied the accuracy of TEC estimates by means of the latitudinal and longitudinal gradients, for the ionosphere at a distance from a station. The latitude errors were established to grow more dramatically, than those of longitude. This could happen, because the S mapping function has an essential latitude dependence. Therefore, at a significant distance, the vertical TEC value converted from the slant TEC is not precisely determined. Hence, the gradients are not precisely determined, either. One should note that, in this case, this error would grow almost linearly, but this is not absolutely so. The real (including the model ones) latitude gradients have a considerably spatially non-uniform character at a distance from a station, and, at distances more than 10°–15°, the TEC estimates based on gradients, start to surpass 10 TECU.
The DCB values obtained through the developed algorithm for GPS satellites agree with the GIM and CODE data, but, for the GLONASS DCB values, the deviation from the CODE data is up to 17 TECU. At the same time, the recovered DCBs (at the IRI-2012 simulation) agree well with the initial data. At large errors in determining DCBs after correcting the slant TEC series, one may observe negative unphysical TEC values.
The developed software may be used to calculate the vertical TEC from the local networks or to locally update ionosphere models [61]. The vertical TEC obtained through this procedure generally agrees with the TEC from the IONOLAB and CODE GIM. The TuRBOTEC data disagree with the Madrigal for a region with few stations and with the SEEMALA-TEC data.

Author Contributions

Conceptualization, Y.Y.; methodology, Y.Y., A.M., and A.V.; software, A.M., and A.V.; validation, Y.Y., A.M. and A.V.; formal analysis, A.M.; investigation, Y.Y., A.M., and A.V.; writing—original draft preparation, Y.Y., A.M., and A.V.; writing—review and editing, Y.Y., A.V.; visualization, A.M., Y.Y.; supervision, Y.Y.; project administration, Y.Y.; funding acquisition, Y.Y. All authors have read and agreed to the published version of the manuscript.

Funding

The study was supported by the Russian Foundation for Basic Research (grant N 18-35-20038) and the Ministry of Education and Science (Basic Research program II.16).

Acknowledgments

We acknowledge the IGS network for the GNSS data, as well as we do the IRI community and T.L. Gulyaeva for the IRI-2012 and IRI-plas models’ source code. The authors thank scipy.org for python library used in the study. The authors give a credit Alexandra Kustavinova for the graphical abstract made for the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hofmann-Wellenhof, B.; Lichtenegger, H.; Collins, J. Global Positioning System: Theory and Practice; Springer-V Sciernce and Business Media LLC: New York, NY, USA, 2001. [Google Scholar]
  2. Lanyi, G.E.; Roth, T. A comparison of mapped and measured total ionospheric electron content using global positioning system and beacon satellite observations. Radio Sci. 1988, 23, 483–492. [Google Scholar] [CrossRef]
  3. Calais, E.; Minster, J.B. GPS detection of an ionospheric perturbation following the January 17, 1994, Northridge earthquake. Geophys. Res. Lett. 1995, 22, 1045–1048. [Google Scholar] [CrossRef]
  4. Afraimovich, E.L.; Astafyeva, E.I.; Demyanov, V.V.; Edemskiy, I.K.; Gavrilyuk, N.S.; Ishin, A.B.; Kosogorov, E.A.; Leonovich, L.A.; Lesyuta, O.S.; Palamartchouk, K.S.; et al. A review of GPS/GLONASS studies of the ionospheric response to natural and anthropogenic processes and phenomena. J. Space Weather. Space Clim. 2013, 3, A27. [Google Scholar] [CrossRef] [Green Version]
  5. Astafyeva, E.; Rolland, L.; Lognonné, P.; Khelfi, K.; Yahagi, T. Parameters of seismic source as deduced from 1 Hz ionospheric GPS data: Case study of the 2011 Tohoku-oki event. J. Geophys. Res. Space Phys. 2013, 118, 5942–5950. [Google Scholar] [CrossRef]
  6. Nesterov, I.; Kunitsyn, V. GNSS radio tomography of the ionosphere: The problem with essentially incomplete data. Adv. Space Res. 2011, 47, 1789–1803. [Google Scholar] [CrossRef]
  7. Kunitsyn, V.E.; Padokhin, A.; Kurbatov, G.A.; Yasyukevich, Y.V.; Morozov, Y.V. Ionospheric TEC estimation with the signals of various geostationary navigational satellites. GPS Solut. 2015, 20, 877–884. [Google Scholar] [CrossRef]
  8. Forte, B.; Aquino, M. On the estimate and assessment of the ionospheric effects affecting low frequency radio astronomy measurements. In Proceedings of the 2011 XXXth URSI General Assembly and Scientific Symposium, Istanbul, Turkey, 13–20 August 2011. [Google Scholar] [CrossRef]
  9. Afraimovich, E.; Yasukevich, Y. Using GPS–GLONASS–GALILEO data and IRI modeling for ionospheric calibration of radio telescopes and radio interferometers. J. Atmos. Sol.-Terr. Phys. 2008, 70, 1949–1962. [Google Scholar] [CrossRef]
  10. Ovodenko, V.; Trekin, V.; Korenkova, N.; Klimenko, M. Investigating range error compensation in UHF radar through IRI-2007 real-time updating: Preliminary results. Adv. Space Res. 2015, 56, 900–906. [Google Scholar] [CrossRef]
  11. Schaer, S.; Overview of GNSS biases. International GNSS Service Workshop on GNSS Biases. 2012. Available online: http://www.biasws2012.unibe.ch/pdf/bws12_1.3.1.pdf (accessed on 31 August 2020).
  12. Mylnikova, A.; Yasyukevich, Y.V.; Kunitsyn, V.; Padokhin, A. Variability of GPS/GLONASS differential code biases. Results Phys. 2015, 5, 9–10. [Google Scholar] [CrossRef] [Green Version]
  13. Choi, B.-K.; Lee, S.-J. The influence of grounding on GPS receiver differential code biases. Adv. Space Res. 2018, 62, 457–463. [Google Scholar] [CrossRef]
  14. Yasyukevich, Y.V.; Mylnikova, A.A.; Kunitsyn, V.E.; Padokhin, A.M. Influence of GPS/GLONASS differential code biases on the determination accuracy of the absolute total electron content in the ionosphere. Geomag. Aeron. 2015, 55, 790–796. [Google Scholar] [CrossRef]
  15. Hong, C.-K.; Grejner-Brzezinska, D.A.; Kwon, J.H. Efficient GPS receiver DCB estimation for ionosphere modeling using satellite-receiver geometry changes. Earth Planets Space 2008, 60, e25–e28. [Google Scholar] [CrossRef] [Green Version]
  16. Jin, R.; Jin, S.; Feng, G. M_DCB: Matlab code for estimating GNSS satellite and receiver differential code biases. GPS Solut. 2012, 16, 541–548. [Google Scholar] [CrossRef]
  17. Li, H.; Xiao, J.; Zhu, W. Investigation and Validation of the Time-Varying Characteristic for the GPS Differential Code Bias. Remote. Sens. 2019, 11, 428. [Google Scholar] [CrossRef] [Green Version]
  18. Wang, J.; Huang, G.; Yang, Y.; Zhang, Q.; Gao, Y.; Zhou, P. Mitigation of Short-Term Temporal Variations of Receiver Code Bias to Achieve Increased Success Rate of Ambiguity Resolution in PPP. Remote. Sens. 2020, 12, 796. [Google Scholar] [CrossRef] [Green Version]
  19. Gordon, W. Incoherent Scattering of Radio Waves by Free Electrons with Applications to Space Exploration by Radar. IEEE Proc. IRE 1958, 46, 1824–1829. [Google Scholar] [CrossRef]
  20. Potekhin, A.P.; Medvedev, A.; Zavorin, A.V.; Kushnarev, D.S.; Lebedev, V.P.; Lepetaev, V.V.; Shpynev, B.G. Recording and control digital systems of the Irkutsk Incoherent Scatter Radar. Geomagn. Aeron. 2009, 49, 1011–1021. [Google Scholar] [CrossRef]
  21. Reinisch, B.; Galkin, I.; Khmyrov, G.M.; Kozlov, A.V.; Bibl, K.; Lisysyan, I.A.; Cheney, G.P.; Huang, X.; Kitrosser, D.F.; Paznukhov, V.V.; et al. New Digisonde for research and monitoring applications. Radio Sci. 2009, 44, 1–15. [Google Scholar] [CrossRef]
  22. Reinisch, B.W.; Galkin, I.A. Global Ionospheric Radio Observatory (GIRO). Earth Planets Space 2011, 63, 377–381. [Google Scholar] [CrossRef] [Green Version]
  23. Mitchell, C.N.; Spencer, P.S. A three-dimensional time-dependent algorithm for ionospheric imaging using GPS. Ann. Geophys. 2003, 46, 687–696. [Google Scholar] [CrossRef]
  24. Jin, S.; Park, J.-U. GPS ionospheric tomography: A comparison with the IRI-2001 model over South Korea. Earth Planets Space 2007, 59, 287–292. [Google Scholar] [CrossRef] [Green Version]
  25. Bust, G.; Garner, T.W.; Ii, T.L.G. Ionospheric Data Assimilation Three-Dimensional (IDA3D): A global, multisensor, electron density specification algorithm. J. Geophys. Res. Space Phys. 2004, 109. [Google Scholar] [CrossRef] [Green Version]
  26. Kunitsyn, V.E.; Nesterov, I.A.; Padokhin, A.; Tumanova, Y.S. Ionospheric radio tomography based on the GPS/GLONASS navigation systems. J. Commun. Technol. Electron. 2011, 56, 1269–1281. [Google Scholar] [CrossRef]
  27. Hernández-Pajares, M.; Juan, J.M.; Sanz, J.; Orús, R.; García-Rigo, A.; Feltens, J.; Komjathy, A.; Schaer, S.C.; Krankowski, A.; Zornoza, J. The IGS VTEC maps: A reliable source of ionospheric information since 1998. J. Geod. 2009, 83, 263–275. [Google Scholar] [CrossRef]
  28. Afraimovich, E.L.; Astafyeva, E.I.; Oinats, A.V.; Yasukevich, Y.V.; Zhivetiev, I. Global electron content: A new conception to track solar activity. Ann. Geophys. 2008, 26, 335–344. [Google Scholar] [CrossRef] [Green Version]
  29. Mannucci, A.J.; Wilson, B.D.; Yuan, D.N.; Ho, C.H.; Lindqwister, U.J.; Runge, T.F. A global mapping technique for GPS-derived ionospheric TEC measurements. Radio Sci. 1998, 33, 565–582. [Google Scholar] [CrossRef]
  30. Schaer, S.; Beutler, G.; Rothacher, M. Mapping and predicting the ionosphere. In Proceedings of the 1998 IGS Analysis Center Workshop, Darmstadt, Germany, 9–11 February 1998; pp. 307–320. [Google Scholar]
  31. Durmaz, M.; Karslioglu, M.O. Regional vertical total electron content (VTEC) modeling together with satellite and receiver differential code biases (DCBs) using semi-parametric multivariate adaptive regression B-splines (SP-BMARS). J. Geod. 2014, 89, 347–360. [Google Scholar] [CrossRef]
  32. Sardón, E.; Zarraoa, N. Estimation of total electron content using GPS data: How stable are the differential satellite and receiver instrumental biases? Radio Sci. 1997, 32, 1899–1910. [Google Scholar] [CrossRef]
  33. Themens, D.R.; Jayachandran, P.T.; Langley, R.B. The nature of GPS differential receiver bias variability: An examination in the polar cap region. J. Geophys. Res. Space Phys. 2015, 120, 8155–8175. [Google Scholar] [CrossRef]
  34. Schaer, S. Mapping and predicting the Earth’s ionosphere using the global positioning system. Ph.D. Thesis, University of Berne, Berne, Switzerland, 1999. Available online: http://ftp.aiub.unibe.ch/papers/ionodiss.ps (accessed on 31 August 2020).
  35. Li, Z.; Yuan, Y.; Li, H.; Ou, J.; Huo, X. Two-step method for the determination of the differential code biases of COMPASS satellites. J. Geod. 2012, 86, 1059–1076. [Google Scholar] [CrossRef]
  36. Start, P.B.; Parker, R.L. Bounded-Variable Least-Squares: An Algorithm and Applications. Comput. Stat. 1995, 10, 129–141. [Google Scholar]
  37. Waterman, M.S. A restricted least squares problem. Technometrics 1974, 16, 135–136. [Google Scholar] [CrossRef]
  38. Zhang, H.; Xu, P.; Han, W.; Ge, M.; Shi, C. Eliminating negative VTEC in global ionosphere maps using inequality-constrained least squares. Adv. Space Res. 2013, 51, 988–1000. [Google Scholar] [CrossRef]
  39. Dow, J.M.; Neilan, R.E.; Rizos, C. The International GNSS Service in a Changing Landscape of Global Navigation Satellite Systems. J. Geod 2009, 83, 191–198. [Google Scholar] [CrossRef]
  40. Arikan, F.; Arikan, O.; Erol, C.B. Regularized estimation of TEC from GPS data for certain midlatitude stations and comparison with the IRI model. Adv. Space Res. 2007, 39, 867–874. [Google Scholar] [CrossRef]
  41. Seemala, G.K. GPS-TEC analysis application. 2017. Available online: https://seemala.blogspot.com/ (accessed on 31 August 2020).
  42. Rideout, W.C.; Coster, A. Automated GPS processing for global total electron content data. GPS Solut. 2006, 10, 219–228. [Google Scholar] [CrossRef]
  43. Bilitza, D.; Altadill, D.; Zhang, Y.; Mertens, C.; Truhlik, V.; Richards, P.; McKinnell, L.-A.; Reinisch, B. The International Reference Ionosphere 2012―A model of international collaboration. J. Space Weather. Space Clim. 2014, 4, A07. [Google Scholar] [CrossRef]
  44. Gulyaeva, T.L.; Bilitza, D. Towards ISO Standard Earth Ionosphere and Plasmasphere Model. In Proceedings of the 39th COSPAR Scientific Assembly, Mysore, India, 14–22 July 2012; pp. 1–39. [Google Scholar]
  45. Cooper, C.; Mitchell, C.N.; Wright, C.J.; Jackson, D.R.; Witvliet, B.A. Measurement of Ionospheric Total Electron Content Using Single-Frequency Geostationary Satellite Observations. Radio Sci. 2019, 54, 10–19. [Google Scholar] [CrossRef]
  46. Krishna, K.S.; Ratnam, D.V. Determination of NavIC differential code biases using GPS and NavIC observations. Geod. Geodyn. 2020, 11, 97–105. [Google Scholar] [CrossRef]
  47. Ma, G.; Maruyama, T. Derivation of TEC and estimation of instrumental biases from GEONET in Japan. Ann. Geophys. 2003, 21, 2083–2093. [Google Scholar] [CrossRef] [Green Version]
  48. Komjathy, A.; Sparks, L.; Wilson, B.D.; Mannucci, A.J. Automated daily processing of more than 1000 ground-based GPS receivers for studying intense ionospheric storms. Radio Sci. 2005, 40, 1–11. [Google Scholar] [CrossRef] [Green Version]
  49. Zhang, Q.; Zhao, Q.; Zhang, H.; Chen, G. BDS Satellites and Receivers DCB Resolution. In China Satellite Navigation Conference (CSNC) 2014 Proceedings: Volume III; Springer: Berlin, Heidelberg, 2014; Volume 305, pp. 187–197. [Google Scholar]
  50. Blewitt, G. An Automatic Editing Algorithm for GPS data. Geophys. Res. Lett. 1990, 17, 199–202. [Google Scholar] [CrossRef] [Green Version]
  51. Klobuchar, J.A. Ionospheric Time-Delay Algorithm for Single-Frequency GPS Users. IEEE Trans. Aerosp. Electron. Syst. 1987, 23, 325–331. [Google Scholar] [CrossRef]
  52. Hernández-Pajares, M.; Juan, J.; Sanz, J.; Zornoza, J. New approaches in global ionospheric determination using ground GPS data. J. Atmos. Sol.-Terr. Phys. 1999, 61, 1237–1247. [Google Scholar] [CrossRef]
  53. Lyu, H.; Hernández-Pajares, M.; Nohutcu, M.; García-Rigo, A.; Zhang, H.; Liu, J. The Barcelona ionospheric mapping function (BIMF) and its application to northern mid-latitudes. GPS Solut. 2018, 22, 67. [Google Scholar] [CrossRef] [Green Version]
  54. Schüler, T.; Oladipo, O.A. Single-frequency single-site VTEC retrieval using the NeQuick2 ray tracer for obliquity factor determination. GPS Solut. 2013, 18, 115–122. [Google Scholar] [CrossRef]
  55. Astafyeva, E.; Zakharenkova, I.; Forster, M. Ionospheric response to the 2015 St. Patrick’s Day storm: A global multi-instrumental overview. J. Geophys. Res. Space Phys. 2015, 120, 9023–9037. [Google Scholar] [CrossRef] [Green Version]
  56. Schaer, S. SINEX BIAS-Solution (Software/technique) INdependent EXchange Format for GNSS Biases Version 1.00. Available online: http://ftp.aiub.unibe.ch/bcwg/format/draft/sinex_bias_100_dec07.pdf (accessed on 7 December 2019).
  57. Afraimovich, E.L.; Astafyeva, E.I.; Yasukevich, Y.V.; Oinats, A.V.; Zhivetiev, I.V. Response of global and regional ionosphere electron content to solar activity changes. Geomag. Aeron. 2008, 48, 187–200. [Google Scholar] [CrossRef]
  58. De La Luz, V.; Gonzalez-Esparza, J.A.; Sergeeva, M.A.; Corona-Romero, P.; González, L.X.; Mejia-Ambriz, J.C.; Valdés-Galicia, J.F.; Aguilar-Rodriguez, E.; Rodriguez-Martinez, M.; Romero-Hernandez, E.; et al. First joint observations of space weather events over Mexico. Ann. Geophys. 2018, 36, 1347–1360. [Google Scholar] [CrossRef] [Green Version]
  59. Sergeeva, M.; Maltseva, O.; Gonzalez-Esparza, J.A.; De La Luz, V.; Corona-Romero, P. Features of TEC behaviour over the low-latitude North-American region during the period of medium solar activity. Adv. Space Res. 2017, 60, 1594–1605. [Google Scholar] [CrossRef]
  60. Nie, W.; Xu, T.; Rovira-Garcia, A.; Zornoza, J.; Sanz, J.; González-Casado, G.; Chen, W.; Xu, G. Revisit the calibration errors on experimental slant total electron content (TEC) determined with GPS. GPS Solut. 2018, 22, 85. [Google Scholar] [CrossRef] [Green Version]
  61. Kotova, D.S.; Ovodenko, V.B.; Yasyukevich, Y.V.; Klimenko, M.V.; Ratovsky, K.G.; Mylnikova, A.A.; Andreeva, E.S.; Kozlovsky, A.E.; Korenkova, N.A.; Nesterov, I.A.; et al. Efficiency of updating the ionospheric models using total electron content at mid- and sub-auroral latitudes. GPS Solut. 2019, 24, 25. [Google Scholar] [CrossRef]
Figure 1. IONOLAB TEC (a), SEEMALA-TEC (b), Madrigal TEC (c). Grey dots on panels (a,b) show absolute slant TEC from different satellites. Grey dots on panel (c) show the vertical TEC for different cells (50°–54°N, 102°–106°E). Orange dots show the median along the 5° × 5° region, while the orange line shows the spline interpolation regarded as the absolute vertical TEC. Panel (d) shows all the mentioned vertical TECs along with the TEC from CODE GIMs (black line). The data are for 1 January 2014.
Figure 1. IONOLAB TEC (a), SEEMALA-TEC (b), Madrigal TEC (c). Grey dots on panels (a,b) show absolute slant TEC from different satellites. Grey dots on panel (c) show the vertical TEC for different cells (50°–54°N, 102°–106°E). Orange dots show the median along the 5° × 5° region, while the orange line shows the spline interpolation regarded as the absolute vertical TEC. Panel (d) shows all the mentioned vertical TECs along with the TEC from CODE GIMs (black line). The data are for 1 January 2014.
Sensors 20 05702 g001
Figure 2. Eliminating the phase ambiguity (a), outliers (b), and cycle slips (c). The panels show the numbers for the corresponding GNSS satellite. Red dots mark the uncorrected phase TEC, blue dots-corrected phase TEC, orange dots-uncorrected pseudorange TEC, grey dots-corrected pseudorange TEC.
Figure 2. Eliminating the phase ambiguity (a), outliers (b), and cycle slips (c). The panels show the numbers for the corresponding GNSS satellite. Red dots mark the uncorrected phase TEC, blue dots-corrected phase TEC, orange dots-uncorrected pseudorange TEC, grey dots-corrected pseudorange TEC.
Sensors 20 05702 g002
Figure 3. Results of simulating the recovery of the absolute vertical TEC (a), the time derivative (b), and the latitudinal gradient (c) for 10 April 2012. Red solid line indicates values (vertical TEC, time derivative, and spatial gradient) recovered by using the second order IV expansion; blue dotted line indicates values recovered by using the first order IV expansion, dot-and-dash black line shows values recovered by using the zero order IV expansion). The solid black exhibits the IRI-2012 data.
Figure 3. Results of simulating the recovery of the absolute vertical TEC (a), the time derivative (b), and the latitudinal gradient (c) for 10 April 2012. Red solid line indicates values (vertical TEC, time derivative, and spatial gradient) recovered by using the second order IV expansion; blue dotted line indicates values recovered by using the first order IV expansion, dot-and-dash black line shows values recovered by using the zero order IV expansion). The solid black exhibits the IRI-2012 data.
Sensors 20 05702 g003
Figure 4. TEC dependence on the elevation for different α. The simulation used the IRI-2012 model. The red solid line shows the IRI-2012 slant TEC data; the red dotted line exhibits the slant TEC converted from the vertical one by using α=1; the blue solid line shows the slant TEC converted from the vertical one by using α = 0.97; the green dotted line shows the same, but with the use of α = 0.94; the black dotted line exhibits the same, but with the use of α 0.87. The blue dot-and-dash line shows the elevation of the satellite.
Figure 4. TEC dependence on the elevation for different α. The simulation used the IRI-2012 model. The red solid line shows the IRI-2012 slant TEC data; the red dotted line exhibits the slant TEC converted from the vertical one by using α=1; the blue solid line shows the slant TEC converted from the vertical one by using α = 0.97; the green dotted line shows the same, but with the use of α = 0.94; the black dotted line exhibits the same, but with the use of α 0.87. The blue dot-and-dash line shows the elevation of the satellite.
Sensors 20 05702 g004
Figure 5. Contribution of measurements from individual satellites to the estimates at different instants.
Figure 5. Contribution of measurements from individual satellites to the estimates at different instants.
Sensors 20 05702 g005
Figure 6. Results of simulating the absolute vertical TEC recovery for different stations as 10 April 2012: (a) mid-latitude IRKJ, α = 0.97; (b) equatorial NTUS, α = 0.87; (c) high-latitude THU2, α = 0.94. The black dotted line shows the results of the TuRBOTEC operation, the red line exhibits the IRI-2012 vertical TEC data, and the blue line shows differences between the results of the TuRBOTEC operation and the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Figure 6. Results of simulating the absolute vertical TEC recovery for different stations as 10 April 2012: (a) mid-latitude IRKJ, α = 0.97; (b) equatorial NTUS, α = 0.87; (c) high-latitude THU2, α = 0.94. The black dotted line shows the results of the TuRBOTEC operation, the red line exhibits the IRI-2012 vertical TEC data, and the blue line shows differences between the results of the TuRBOTEC operation and the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Sensors 20 05702 g006
Figure 7. Absolute vertical TEC recovered from real the GPS/GLONASS data for IRKJ (red line). Panel (a) is for 5 March 2015 (Крmax = 2.3), panel (b) 17 March 2015 (Крmax = 7.7). The dot-and-dash blue line presents the JPL data, the black dotted line presents the CODE data, and the red solid line shows the TuRBOTEC values.
Figure 7. Absolute vertical TEC recovered from real the GPS/GLONASS data for IRKJ (red line). Panel (a) is for 5 March 2015 (Крmax = 2.3), panel (b) 17 March 2015 (Крmax = 7.7). The dot-and-dash blue line presents the JPL data, the black dotted line presents the CODE data, and the red solid line shows the TuRBOTEC values.
Sensors 20 05702 g007
Figure 8. Normalized ΔIV histograms (difference between IV values from alternative and TuRBOTEC data). (a) GIM CODE vs. TuRBOTEC; (b) Madrigal vs. TuRBOTEC; (c) IONOLAB vs. TuRBOTEC; (d) SEEMALA-TEC vs. TuRBOTEC. Panel (e) shows the distribution between the difference in JPL GIM TEC and CODE GIM TEC; panel (f) compares techniques with (TuRBOTEC) and without constraint. The data are for the IRKJ over 2014.
Figure 8. Normalized ΔIV histograms (difference between IV values from alternative and TuRBOTEC data). (a) GIM CODE vs. TuRBOTEC; (b) Madrigal vs. TuRBOTEC; (c) IONOLAB vs. TuRBOTEC; (d) SEEMALA-TEC vs. TuRBOTEC. Panel (e) shows the distribution between the difference in JPL GIM TEC and CODE GIM TEC; panel (f) compares techniques with (TuRBOTEC) and without constraint. The data are for the IRKJ over 2014.
Sensors 20 05702 g008
Figure 9. Simulation for the TEC latitude (a–c) and longitude (d–f) gradients recovery for IRKJ (a,d), NTUS (b,e), and THU2 (c,f). The black dotted line presents the results of the TuRBOTEC operation; the red line exhibits the IRI-2012 TEC gradient data, and the blue line shows the difference between the results of the TuRBOTEC operation and those from the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Figure 9. Simulation for the TEC latitude (a–c) and longitude (d–f) gradients recovery for IRKJ (a,d), NTUS (b,e), and THU2 (c,f). The black dotted line presents the results of the TuRBOTEC operation; the red line exhibits the IRI-2012 TEC gradient data, and the blue line shows the difference between the results of the TuRBOTEC operation and those from the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Sensors 20 05702 g009
Figure 10. Errors in TEC estimations at different UT for various distances from IRKJ. ΔIV, is the difference between the model-specified values and the model-recovered data, panel (a) and panel (b) are for the distance from the station by latitude and longitude, respectively. The IRI-2012 simulation as of 10 April 2012 for the IRKJ geometry.
Figure 10. Errors in TEC estimations at different UT for various distances from IRKJ. ΔIV, is the difference between the model-specified values and the model-recovered data, panel (a) and panel (b) are for the distance from the station by latitude and longitude, respectively. The IRI-2012 simulation as of 10 April 2012 for the IRKJ geometry.
Sensors 20 05702 g010
Figure 11. Simulating the TEC time derivative recovery for IRKJ (a), NTUS (b), and THU2 (c) stations. The black dotted line presents the results of the TuRBOTEC operation, the red line exhibits the IRI-2012 TEC time derivative data, and the blue line shows the difference between the results of the TuRBOTEC operation and those from the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Figure 11. Simulating the TEC time derivative recovery for IRKJ (a), NTUS (b), and THU2 (c) stations. The black dotted line presents the results of the TuRBOTEC operation, the red line exhibits the IRI-2012 TEC time derivative data, and the blue line shows the difference between the results of the TuRBOTEC operation and those from the IRI-2012 model. The green line shows the similar differences, but for the IRI-plas simulation. The axis for the blue and green lines are on the right, for the red and black lines on the left.
Sensors 20 05702 g011
Figure 12. GNSS-based TEC time derivative for IRKJ. (a) 5 March 2015 (Kpmax = 2.3); (b) 17 March 2015 (Kpmax = 7.7). The blue dot-and-dash line presents the JPL data, the black dotted line denotes the CODE data, and the red solid line displays the TuRBOTEC-obtained values.
Figure 12. GNSS-based TEC time derivative for IRKJ. (a) 5 March 2015 (Kpmax = 2.3); (b) 17 March 2015 (Kpmax = 7.7). The blue dot-and-dash line presents the JPL data, the black dotted line denotes the CODE data, and the red solid line displays the TuRBOTEC-obtained values.
Sensors 20 05702 g012
Figure 13. Simulating the DCB recovery from GPS (a) and GLONASS (b). The black triangles are the initial values; the red circles present the recovered values; the blue line shows the difference between the results of the TuRBOTEC operation and the initial values. The data are in TECU.
Figure 13. Simulating the DCB recovery from GPS (a) and GLONASS (b). The black triangles are the initial values; the red circles present the recovered values; the blue line shows the difference between the results of the TuRBOTEC operation and the initial values. The data are in TECU.
Sensors 20 05702 g013
Figure 14. DCBs estimate for GPS (a) and GLONASS (b) satellites from the IRKJ 1 March 2015 data. The red circles show the results of the TuRBOTEC algorithm operation; the black diamonds present the CODE data.
Figure 14. DCBs estimate for GPS (a) and GLONASS (b) satellites from the IRKJ 1 March 2015 data. The red circles show the results of the TuRBOTEC algorithm operation; the black diamonds present the CODE data.
Sensors 20 05702 g014
Figure 15. GNSS absolute slant TEC after least squares without constraint (black dots) and TuRBOTEC techniques (blue dots).
Figure 15. GNSS absolute slant TEC after least squares without constraint (black dots) and TuRBOTEC techniques (blue dots).
Sensors 20 05702 g015
Figure 16. Influence of intra-day DCB variations on pseudorange TEC (a, blue dots) and vertical TEC estimates (b, black line). Red line on panel (b) shows the phase-without-pseudorange solution. The data are for IRKJ 12 July 2009.
Figure 16. Influence of intra-day DCB variations on pseudorange TEC (a, blue dots) and vertical TEC estimates (b, black line). Red line on panel (b) shows the phase-without-pseudorange solution. The data are for IRKJ 12 July 2009.
Sensors 20 05702 g016
Table 1. Position of GNSS stations.
Table 1. Position of GNSS stations.
StationLat, °Lon, °MLat, °MLon, °α
IRKJ52.2104.347.7178.30.97
NTUS1.3103.7−7.2176.30.87
THU276.5291.283.827.10.94

Share and Cite

MDPI and ACS Style

Yasyukevich, Y.; Mylnikova, A.; Vesnin, A. GNSS-Based Non-Negative Absolute Ionosphere Total Electron Content, its Spatial Gradients, Time Derivatives and Differential Code Biases: Bounded-Variable Least-Squares and Taylor Series. Sensors 2020, 20, 5702. https://doi.org/10.3390/s20195702

AMA Style

Yasyukevich Y, Mylnikova A, Vesnin A. GNSS-Based Non-Negative Absolute Ionosphere Total Electron Content, its Spatial Gradients, Time Derivatives and Differential Code Biases: Bounded-Variable Least-Squares and Taylor Series. Sensors. 2020; 20(19):5702. https://doi.org/10.3390/s20195702

Chicago/Turabian Style

Yasyukevich, Yury, Anna Mylnikova, and Artem Vesnin. 2020. "GNSS-Based Non-Negative Absolute Ionosphere Total Electron Content, its Spatial Gradients, Time Derivatives and Differential Code Biases: Bounded-Variable Least-Squares and Taylor Series" Sensors 20, no. 19: 5702. https://doi.org/10.3390/s20195702

APA Style

Yasyukevich, Y., Mylnikova, A., & Vesnin, A. (2020). GNSS-Based Non-Negative Absolute Ionosphere Total Electron Content, its Spatial Gradients, Time Derivatives and Differential Code Biases: Bounded-Variable Least-Squares and Taylor Series. Sensors, 20(19), 5702. https://doi.org/10.3390/s20195702

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop