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

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.


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 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.

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 (  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 Sensors 2020, 20, 5702 4 of 20 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.
Sensors 2020, 20, x FOR PEER REVIEW 4 of 24 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.

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:

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: where I M 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, I V 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 I DCB 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 I V =I V (φ 0 , l 0 , t 0 ), where φ 0 , l 0 , and t 0 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 I V 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 I V Taylor series expansion in (1) only on space. In our research, we expand the vertical TEC function I V (φ, l, t) into the Taylor series on space and time: where ∆φ, ∆l are the difference between latitude, longitude of the ionospheric pierce point and the point of station (φ 0 , l 0 ), ∆t is the difference between time of measurement and t 0 .
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 I P 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): 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 I Exp . (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 (I V ) 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.     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: 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_ϕ = ∂ 2 IV/∂ϕ 2 , and Gq_l = ∂ 2 IV/∂l 2 are the linear and quadratic spatial TEC gradients, and Gt = ∂IV/∂t and Gq_t = ∂ 2 IV/∂t 2 are the first and second time derivatives. Equation (3) represents the secondorder Taylor series expansion 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 Meanwhile, using the second order provides a sufficient accuracy to be able to neglect higher orders. Therefore, we use the expansion I V into the Taylor series up to the second order: where I V is the absolute vertical TEC value, ∆φ (∆l) is the latitude (longitude) difference between the ionospheric point coordinate φ(l) and that of the station φ 0 (l 0 ), ∆t is the difference between the measurement time t and the time t 0 for which the calculation is performed. Further, G φ = ∂I V /∂φ, G l = ∂I V /∂l, G q_φ = ∂ 2 I V /∂φ 2 , and G q_l = ∂ 2 I V /∂l 2 are the linear and quadratic spatial TEC gradients, and G t = ∂I V /∂t and G q_t = ∂ 2 I V /∂t 2 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]: where R E = 6371 km is the Earth radius, h max = 450 km is the ionospheric point height, and α is the correcting coefficient. Figure 4 presents the slant TEC dependence, I S , on the elevation when converting I S = S·I V 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 latitudedependent. 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 = ∑U k (5) for the set of selected instant t k , for which we estimate the parameters. For computation, we represent (5) in the form (7): 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, is the i th instant of measurement for the j th satellite, ∆ , is the time difference between the current measurement and 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 = U k (5) for the set of selected instant t k , for which we estimate the parameters. For computation, we represent (5) in the form (7): where I Exp is the experimental phase TEC measurements with the eliminated phase ambiguity obtained after Stage 4 of the algorithm, Θ is the Heaviside step function, t i j is the ith instant of measurement for the jth satellite, ∆t i,k j is the time difference between the current measurement t i j and the time t k , 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 t k -estimate). The 1/S factor in Equation (6) causes the measurements at high elevations to produce the greatest contribution. For each time instant G k q_t , and I DCB, 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: 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). 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 I DCB , we calculate the parameters simultaneously for different time instants over 24 h, thus solving a consistent set of Equations (7). The temporal resolution for t k may vary from 2 h to 5 min.
To analyze the satellite DCB separately, we apply the often-used zero-mean condition [34]: 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]. 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 t k may vary from 2 h to 5 min. To analyze the satellite DCB separately, we apply the often-used zero-mean condition [34]: 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·S 10 (θ + 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 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·S 10 (θ + 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.

Technique Validation and Discussion
We validated data based on ionosphere modeling and TEC products from alternative software mentioned in Section 2. 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).

Absolute Total Electron Content
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. 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).

Absolute Total Electron Content
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.   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 March 17, 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.  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 ∆I V difference distribution between the vertical TEC data in the Irkutsk region.  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. 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.

Spatial Gradients. Accuracy of Determining TEC at a Growing Distance from a Station
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). 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. 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. 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 Figures 6 and 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. 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. 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 Figures 6 and 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 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 April 10, 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 March 1, 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- 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 January 5, 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.

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): 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

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): where, I const 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.

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 highlatitude 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

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.