Near real-time estimation of ionosphere vertical total electron content from GNSS satellites using B-splines in a Kalman filter

Although the number of terrestrial global navigation satellite system (GNSS) receivers supported by the International GNSS Service (IGS) is rapidly growing, the worldwide rather inhomogeneously distributed observation sites do not allow the generation of high-resolution global ionosphere products. Conversely, with the regionally enormous increase in highly precise GNSS data, the demands on (near) real-time ionosphere products, necessary in many applications such as navigation, are growing very fast. Consequently, many analysis centers accepted the responsibility of generating such products. In this regard, the primary objective of our work is to develop a near real-time processing framework for the estimation of the vertical total electron content (VTEC) of the ionosphere using proper models that are capable of a global representation adapted to the real data distribution. The global VTEC representation developed in this work is based on a series expansion in terms of compactly supported B-spline functions, which allow for an appropriate handling of the heterogeneous data distribution, including data gaps. The corresponding series coefficients and additional parameters such as differential code biases of the GNSS satellites and receivers constitute the set of unknown parameters. The Kalman filter (KF), as a popular recursive estimator, allows processing of the data immediately after acquisition and paves the way of sequential (near) real-time estimation of the unknown parameters. To exploit the advantages of the chosen data representation and the estimation procedure, the Bspline model is incorporated into the KF under the consideration of necessary constraints. Based on a preprocessing strategy, the developed approach utilizes hourly batches of GPS and GLONASS observations provided by the IGS data centers with a latency of 1 h in its current realization. Two methods for validation of the results are performed, namely the self consistency analysis and a comparison with Jason-2 altimetry data. The highly promising validation results allow the conclusion that under the investigated conditions our derived near real-time product is of the same accuracy level as the so-called final post-processed products provided by the IGS with a latency of several days or even weeks.


Introduction
The ionosphere constitutes the upper part of the atmosphere, extending from approximately 60 to 1500 km above the Earth's surface, enriched with free electrons and ions (Schaer, 1999).The knowledge of the structure and dynamics of the ionospheric plasma has great importance for various scientific applications and services, such as telecommunication through radio signals, point positioning based on global navigation satellite systems (GNSSs) (Gao et al., 2006;Le et al., 2009), the monitoring of space weather events such as solar flares and coronal mass ejections (Monte-Moreno and Hernández-Pajares, 2014; Wang et al., 2016), the investigations of ionospheric anomalies preceding or following a natural hazard such as earthquakes (Liu, 2004), and investigation of the ionospheric effects on thermospheric mass density, wind, and temperature (Jee et al., 2008;Lee et al., 2012).
Published by Copernicus Publications on behalf of the European Geosciences Union.

E. Erdogan et al.: Near real-time estimation of ionosphere VTEC
The ionospheric plasma density varies with time and location and exhibits a coupled system with its environment: the Sun, the Earth's lower atmosphere, the thermosphere and the magnetosphere (Heelis, 2014).Interactions between the thermospheric neutral winds and ionized plasma drive the ionospheric charged particles in motion and lead to separation of charges, resulting in the creation of a polarized electrical field E. During daytime, the combination of this eastward-oriented electrical field E with the Earth's horizontally northward-oriented magnetic field B causes an upward E × B drift of the ionized plasma.Following the upward drift, the plasma diffuses through the magnetic field lines and creates two crests with high ionization at both sides of the magnetic equator, which is also known as fountain effect (Schunk and Nagy, 2009).In addition, at high latitudes, the interaction of the Earth's magnetosphere with the interplanetary magnetic field attached to the solar winds as well as to the space weather events increases the complexity of the system.Today, the advances in space-geodetic techniques, such as terrestrial GNSS, spaceborne radio occultations to low-Earth-orbiting (LEO) satellites as well as satellite altimetry, facilitates the monitoring of the structure of the ionosphere with an improved spatial and temporal resolution.Particularly, GNSS offers an attractive alternative to traditional methods, such as ionosondes, for monitoring the electron content within the ionosphere in terms of volume and global data distribution.
The International GNSS Service (IGS) delivers large volumes of GNSS data with different latencies (e.g., real time, hourly) acquired from continuously operating terrestrial GNSS receivers distributed worldwide.The four IGS Ionosphere Associate Analysis Centers (IAACs), namely the Jet Propulsion Laboratory (JPL), the Center for Orbit Determination in Europe (CODE), the European Space Operations Center of the European Space Agency (ESOC) and the Universitat Politècnica de Catalunya (UPC), monitor the ionosphere and evaluate relevant parameters using dual frequency GNSS receivers.Several modeling approaches for ionospheric parameters have been proposed.A common approach, generally denoted as single-layer model (SLM), is based on the assumption that the electrons in the ionosphere are concentrated within a thin shell at a fixed altitude above the Earth (Mannucci et al., 1998;Komjathy and Langley, 1996).The spatial variations of electron content in this single layer are represented by a proper mathematical model such as spherical harmonics (Schaer, 1999), B-splines (Schmidt, 2007;Schmidt et al., 2008), spatially defined total electron content (TEC) grids (Skone, 1999), polynomials (Komjathy and Langley, 1996;Komjathy, 1997), wavelets (Schmidt, 2007), or the MARS (Durmaz et al., 2010) and BMARS approaches (Durmaz and Karslioglu, 2015).Furthermore, three-dimensional models that consider the variation with altitude were also studied (Hernández-Pajares et al., 1999;Mitchell and Spencer, 2003;Liu, 2004;Zeilhofer et al., 2009;Limberger et al., 2013).The IGS delivers post-processed global VTEC products by combining VTEC maps from the different analysis centers mentioned before (Hernández-Pajares et al., 2011).The combination of these products computed with independent algorithms performed by the analysis centers contributes to the accuracy, the integrity and the availability of the global VTEC.
In this context, one of the main goals of this study is to develop an alternative approach that contributes to these modeling efforts by generating near real-time products.To be more specific, series expansions in terms of tensor products of compactly supported B-spline functions are used to represent the spatial variation of the global VTEC; introductory studies on this topic were published by Schumaker and Traas (1991), Schmidt (2007), Schmidt et al. (2008) and Zeilhofer (2008).The advantages of B-spline representations for VTEC modeling in comparison to spherical harmonics were also revealed by Schmidt et al. (2011).The results of this study show that the B-spline series approach features in considerable improvements if the input data is distributed heterogeneously.
Considering the increasing demands on high-precision (near) real-time global ionosphere products, including global VTEC maps, an estimation strategy becomes obviously important to appropriately handle the large amount of ionosphere data as well as processing the data once it is available.The Kalman filter (KF) is a popular filtering technique and does not require storing measurements of the past since it is in the batch filtering for estimating the current state of a system and data can be processed immediately after acquisition (Gelb, 1974;Grewal and Andrews, 2008).The KF has been successfully applied in (near) real-time ionosphere modeling by many authors; see Komjathy (1997), Liu (2004) and Anghel et al. (2009).Here, apart from the other studies, we focus on the implementation of a B-spline representation into the KF for ionospheric parameter estimation.
In summary, the goal of the present study is to develop a near real-time processing framework to globally monitor the spatial and temporal variations within the ionosphere by exploiting the advantages of the B-spline series expansion and the recursive filtering using GPS and GLONASS measurements.Finally, it shall be discussed how the quality of this near real-time product is in comparison to the so-called final, i.e., the post-processed high-quality VTEC products of the IGS and its IAACs provided with a latency of days or even weeks.
The paper is outlined as follows.Section 2 explains the theoretical foundations for obtaining ionosphere observations from raw hourly GNSS data.In Sect. 3 the B-spline representation for global VTEC modeling is explained.Section 4 gives the background for the sequential estimation algorithm using the KF.The section comprises subparts, introducing the measurement model of the filter, the definition of model constraints, the prediction model of the filter, the foundation of the Kalman filtering, and the handling of the model constraints in the filter and data editing concepts.Sec-tion 5 summarizes the entire near real-time estimation procedure from an application point of view.Note that the exemplified generation of VTEC maps in this paper is based just on GNSS data.Validation concepts -partly based on satellite altimetry data -and the results are given in Sect.6.Finally, Sect.7 provides the conclusion and future work.
2 Ionosphere TEC measurements from GNSS The free electrons within the ionosphere affect the propagation of electromagnetic waves, i.e., they cause a frequencydependent delay in the transmitted radio signals.Although the ionospheric effect on GNSS signals is not desirable in positioning and navigation, it provides valuable information for the investigation of the electron content of the ionosphere.The magnitude of the delay depends on the electron density N e along the ray path between the satellite s and the receiver r and is proportional to the slant total electron content (STEC) defined as N e is measured in electrons per cubic meter.Following Eq. ( 1), the vertical total electron content (VTEC) is defined as the integration of the electron density along the vertical, i.e., in height direction.
In order to extract ionospheric information from dualfrequency GNSS measurements, the geometry-free linear combination can be used (Ciraolo et al., 2007).Firstly, the equations leading to the ionospheric observables P s r,I and L s r,I derived from combinations of the pseudo-range (code) and carrier-phase measurements, respectively, are defined as where P s r,f i and s r,f i with i ∈ {1, 2} stand for the pseudorange and the carrier-phase measurements observed simultaneously by the receiver r from satellite s.The subscripts 1 and 2 denote the two carrier frequencies f 1 and f 2 with the corresponding wavelengths λ 1 and λ 2 .B r and B s are the receiver and satellite inter-frequency biases (IFBs) for the carrier phase observations; similarly, b r and b s stand for the differential code biases (DCBs) of the receiver and the satellite.α is a frequency-dependent constant factor.Furthermore, C s arc,r is the ambiguity bias of the carrier-phase, and P and L account for the measurement errors.The pseudo-range measurements are rather noisy but unambiguous, while the carrier-phase data are significantly more precise but biased.To exploit the precision of the phase measurements, an offset CPB s r , including the terms IFB, DCB and C s arc,r , is computed by averaging the differences between L s r,I and P s r,I for every continuous arc that shares a common phase bias (Ciraolo et al., 2007;Mannucci et al., 1998) where N is the number of observations continuously measured along the arc.An elevation-dependent threshold is established to select the more precise observations with higher elevation angle for the computation of CPB s r .Then, a leveled geometry-free phase observation L s r,4 for a continuous arc becomes The leveling technique is applied to the hourly data sets of GPS and GLONASS observations obtained from the IGS data servers.

Global VTEC representation with B-splines
As already mentioned in the introduction, appropriate approaches for representing VTEC are two-dimensional series expansions in terms of spherical harmonics or B-spline functions.In the latter case the basis functions can be set up by tensor products of polynomial and trigonometric B-splines.
To be more specific, the B-spline representation of the global VTEC reads VTEC(ϕ, λ) = (Schmidt et al., 2011), where the quantities d J 1 ,J 2 k 1 ,k 2 mean the unknown series coefficients, N 2 J 1 ,k 1 (ϕ) are the end-point interpolating polynomial B-spline functions of order three depending on the latitude ϕ, and T 2 J 2 ,k 2 (λ) are the trigonometric B-splines of order three depending on the longitude λ.The values J 1 and J 2 are called levels, the numbers k 1 and k 2 define the geometrical positions of the two-dimensional basis functions on the sphere.The value K J 1 stands for the number of polynomial B-spline functions according to the associated level J 1 ; its numerical value is given by K J 1 = 2 J 1 + 2. Similarly, the number of trigonometric Bspline functions K J 2 for the level J 2 is defined by and the trigonometric B-spline functions T 2 J 2 ,k 2 (λ) are compactly supported, which means the functions are different from zero only within a small subinterval (see Schumaker and Traas, 1991;Schmidt et al., 2011Schmidt et al., , 2015, and the references therein).
The selection of appropriate resolution levels J 1 and J 2 requires the consideration of different criteria, namely the distribution of the input data, the computational burden and the The circles and triangles represent the data locations for GLONASS and GPS, respectively, and the colors indicate the corresponding VTEC magnitudes related to the ionospheric pierce points (IPP) obtained from a reference map.In all the four panels at the bottom and on the left side the K J 2 trigonometric B-spline functions T 2 The straight black lines on the maps show the corresponding knot point locations.The red line shows the prime meridian at Greenwich.desired level of smoothness.Figure 1 clearly shows that the higher the level chosen, the larger the number of B-spline basis functions, i.e., the higher the resolution of the VTEC representation.Appropriate levels of the B-splines can be obtained from the distribution of the input data (Schmidt et al., 2011(Schmidt et al., , 2015)).This procedure has already been applied successfully for ionosphere modeling from GPS occultation data (Limberger, 2015).In case of near real-time applications, the level values may have to fulfill computation time constraints.A trade-off between the resolution level and the computational burden becomes another necessary issue since the higher the resolution, the larger the number of unknown coefficients and the higher the time consumption the parameter estimation procedure based on Kalman filtering requires.The GNSS data distribution is generally rather inhomogeneous, with large data gaps, especially over the oceans.Bspline representations with high level values can result in ba-sis functions without any data support.For instance, all the basis functions in Fig. 1a are supported by GNSS observations.However, this is not the case for the different levels of representations shown in the other three panels of Fig. 1.Conversely, a representation with level values chosen too low may not describe the spatial VTEC variations sufficiently and may cause an undesired smoothing, i.e., a loss of information.For example, a reference VTEC map used as input observation is illustrated in Fig. 2a. Figure 2b depicts the reconstructed map using the B-spline levels J 1 = 3 and J 2 = 1: it smooths the details too much.An increase in the level values to J 1 = 4 and J 2 = 3 results in a better reconstruction of the reference map, as is shown in Fig. 2c.Consequently, in order to achieve a representation quality comparable to IGS products while considering the aforementioned criteria, the B-spline levels are set to J 1 = 4 and J 2 = 3 in this study.6) with Bspline levels J 1 = 3 and J 2 = 1; (c) same as (b) but with B-spline levels J 1 = 4 and J 2 = 3.

Sequential estimation of global VTEC
The KF (Kalman, 1960) was used in this work as a sequential estimator to compute the ionospheric parameters in near real time.Since the KF is of recursive nature, measurements from the past do not need to be stored (Gelb, 1974).The current state is updated as soon as new observations are available.This means a crucial advantage for (near) real-time applications because this way the filter allows assimilation of observations as soon as possible without waiting for another group of observations.
Here, for the estimation of the ionospheric target parameters, the system equations, including the measurement model and the prediction model, are linear and the observations are assumed to have a Gaussian distribution.Then the KF provides an optimal recursive estimator in terms of minimum variance estimation (see Grewal and Andrews, 2008;Simon, 2006;Gelb, 1974).The linear system of equations in a discrete form is defined as where k is the time stamp, F k is the transition matrix, β k is the vector of the unknown parameters, y k is the vector of the measurements and X k is the corresponding design matrix.The measurement error vector e k and the vector w k of the process noise are assumed to be white noise vectors with the expectation values E(e k ) = 0 and E(w k ) = 0; furthermore, the covariance matrices y and w , respectively, fulfill the assumptions and where δ k,l is the delta symbol with δ k,l = 1 if k = l and δ k,l = 0 for k = l.Equation ( 10) means that the vectors w k and e l are assumed to be mutually independent.

Measurement model
For the sake of clarity, the GNSS STEC measurement L s r,4 from Eq. ( 5) is redefined to express both GPS and GLONASS measurement models separately as A widely accepted mapping function, namely the modified single-layer mapping (MSLM) function is defined as (Schaer, 1999) with the single layer is denoted as the ionospheric pierce point (IPP).The angles z and z are the zenith angles of the satellite at the receiver position and the IPP.The two observation vectors y GPS and y GLO for the GPS and the GLONASS measurements, respectively, build the measurement vector y as given in Eq. ( 8).The state vector β consists of the sub-vector d = (d J 1 ,J 2 k 1 ,k 2 ) of the unknown Bspline coefficients d J 1 ,J 2 k 1 ,k 2 as defined in Eq. ( 6), and the subvectors b r,GPS , b r,GLO , b s GPS and b s GLO of the receiver and satellite DCBs.Consequently, the vectors β and y read Although the satellite systems GPS and GLONASS refer to the same space geodetic technique, they are operated by different agencies with a different design, constellation and signal structure, which can lead to different sensitivities within the parameter estimation.To account for this fact, instead of assigning one variance factor, an individual variance factor for each observation group is introduced.Furthermore, it is assumed that the vectors y GPS and y GLO are uncorrelated.Therefore, the measurement covariance matrix y reads herein σ 2 GPS and σ 2 GLO are the two unknown variance factors mentioned before, P GPS and P GLO are given positive definite weight matrices.Since the precision of the GNSS measurement depends on the elevation angle, a common strategy is to adapt an elevation-dependent weighting scheme for each individual observation.Therefore, a diagonal element p ii of the weight matrix P = (p ij ) is defined by (adapted from Wang et al., 2015) where z i is the zenith angle of the ith measurement at a given GNSS receiver location.The non-diagonal elements p ij with i = j are usually set to zero.

Model constraints
The modeling approach includes constraints defined for the different groups of unknown parameters.One group of constraints preserves the spherical geometry (Schumaker and Traas, 1991) since the two-dimensional B-spline model as introduced in Eq. ( 6) is applied to the representation of a function defined on a sphere.The constraints to be applied can be summarized as pole constraints, pole continuity constraints and longitude periodicity constraints.The first group of constraints forces the VTEC values to be the same at the poles.
The second group ensures the equality of tangent planes of VTEC values at the poles, and the last group is to preserve continuity of the VTEC values at the longitudinal boundaries.Considering these requirements the constraint equation reads with the matrix H d of given coefficients (Schumaker and Traas, 1991).To handle the rank deficiency problem related to the satellite and receiver DCBs, a zero mean condition is usually applied by the IAACs within the adjustment procedure.In the following, we also rely on this assumption and introduce the constraint equations (Schaer, 1999).Herein H GPS and H GLO are matrices of given coefficients, respectively.

Prediction model
Many time-varying models for representing the ionospheric dynamics are based on a KF approach, for instance, the physics-based model developed by Scherliess et al. (2006) for electron density modeling.The dynamic system within a KF is often realized by empirical models such as the random walk (Komjathy and Langley, 1996) or the Gauss-Markov process (Liu et al., 2004;Schunk et al., 2004).
The selection of a proper coordinate system is probably one of the most important issues in monitoring the temporal variations of the ionosphere.The KF problem can be solved both in a Sun-fixed (Komjathy and Langley, 1996) or in an Earth-fixed coordinate system (Dettmering et al., 2011).Since the effect of the Earth's diurnal motion is mitigated, the ionosphere varies much slower in a Sun-fixed system and could be assumed as static for a certain time interval.Thus, we follow this argumentation and handle the global VTEC representation according to Eq. ( 6) in a Sun-fixed coordinate system spanned by the solar geomagnetic longitude, the geomagnetic latitude and the radial distance from the origin, i.e., the geocenter.In this coordinate system not only VTEC itself but also the B-spline coefficients vary rather slowly.Hence, a random walk approach (Gelb, 1974) could be performed for the prediction of the coefficients from one measurement epoch to the next.Moreover, the satellite and receiver DCBs are quite stable over a long period, even over a few days (Schaer, 1999).Consequently, they can also be modeled with a random walk approach.In summary, all parameters of the ionospheric state vector are treated as random walk processes in the time domain, which leads to a transition matrix F k as introduced in Eq. ( 7) equal to the identity matrix I.
The covariance matrix w of the process noise stands for uncertainties introduced by deficiencies in the model.Here, w is defined as a diagonal matrix with a variable but unknown precision factor for each component of the ionospheric state vector such that where I n , I m , I k , I l and I h are identity matrices of different sizes.The quantities GLO are the unknown variance factors of the process noise covariance matrix referring to the corresponding B-spline coefficients as well as to the receiver and satellite DCBs for GPS and GLONASS, respectively.If model uncertainties cannot be described precisely, a common practice is manually conducting tests via multiple runs of the filter (Maybeck, 1979).However, alternatively, adaptive methods can be applied to compute the process noise covariance matrix as well as the measurement covariance matrix in Eq. ( 15) during the running time (see Yang, 2010).In this research, the process noise covariance matrix was constructed from an extensive data analysis.

Kalman filtering
The solution of the estimation problem as defined in Eqs. ( 7) and ( 8) by using a KF consists generally of the sequential application of a prediction step -also known as time update -and a correction step -known as the measurement update.In the prediction step, the current ionospheric state vector β k−1 and its covariance matrix β,k−1 are propagated from the time epoch k − 1 to the next time epoch k using a proper prediction model in order to obtain a predicted state vector β − k and the related predicted covariance matrix − β,k of the next step k as where the symbol "−" indicates predicted values.In our application, the sampling interval, i.e., the step size from one epoch k − 1 to the next epoch k is set to 5 min.
Once the prediction step is performed, the corrected state vector and its covariance matrix are computed by incorporating the new allocated ionospheric measurements as where β k and β,k are the updated (or corrected) state vector and the covariance matrix.The so-called gain matrix behaves like a weighting factor between the new measurements and the predicted state and is defined as The covariance matrix computed by Eq. ( 23) does not guarantee symmetry and positive definiteness due to computer rounding errors or an ill-posed character of the problem (Grewal and Andrews, 2008).In order to preserve the symmetry and the positive definiteness, an alternative form of Eq. ( 23) can be obtained as where both terms are symmetric; the first term in Eq. ( 25) is additionally positive definite, and the second one is positive semi-definitive.

Kalman filtering with constraints
Two types of equality constraints were defined by Eqs. ( 17) and ( 18).In order to simplify the notation we define the more general equation of constraints -matrix H and vector q are given -which has to be incorporated into the KF.To solve such a problem, numerous approaches have been proposed, for instance, a model reduction or the Kalman gain restriction method (see Simon, 2010, and references therein).In our study, the socalled estimate projection method is applied and realized as an additional step following the measurement update of the KF.The idea behind this method is to project the estimated state vector β k onto the constraint surface.In other words, it essentially solves the constrained minimization problem min (Simon, 2002), where P k is a symmetric positive definitive weighting matrix and β k means the projected state vector.
The former can be selected as The solution of the minimization problem (Eq.27) is given as wherein the gain matrix reads 4.6 Data editing: pre and post-processing of the filter The composition of the ionospheric state vector β as defined in Eq. ( 14) can change in time.For example, the current ob-  k − 1.Thus, an additional unknown DCB has to be considered in the state vector β k , and an additional preprocessing step prior to the measurement update step has to be carried out for diagnosing the state vector.To be more specific, the composition of the state vector has to be checked at every cycle of the filter to detect a new allocation or a loss of a GNSS receiver as well as the status of the GNSS satellites.If a DCB of the predicted state vector has to be deleted since a receiver is lost, the corresponding row and column also have to be removed from the predicted covariance matrix.Conversely, if a new DCB has to be added, the predicted state vector and the covariance matrix are extended and filled up with predefined values.
The filter post-processing step refers to secondary tasks that are not directly related to the filter but to the generation of products.For instance, the estimated ionospheric parameters and their covariance matrix in combination with other relevant data that would be informative for diagnosing and monitoring the filter results are stored in a database.

Near real-time estimation and product generation
Although the presented filter is capable of running in real time, measurements can only be assimilated with a time delay due to the latency arising from the availability of hourly GNSS observations that have been provided by the IGS data servers with at least a 1 h delay.Furthermore, downloading and processing of the raw GNSS data as well as filtering introduce additional delays.Therefore, the presented approach is called near real time and is capable of generating global VTEC products usually with less than a 1.5 h delay.
Figure 3 shows the overall steps of the presented approach: the hourly data processing, the filtering of the hourly data and the product generation.The routines for downloading raw GNSS data, the computation of the ionospheric observable from the raw data and storing of the relevant data into a database are accomplished within the hourly data processing step.A data preprocessing module, including the following steps, runs sequentially at the beginning of every hour within the implemented software.After the acquisition of the observations, a cutoff elevation angle of 10 • is applied to eliminate the very noisy measurements.Data arcs containing cycle slips are detected and then split into parts.The arcs with a number of observations less then a given threshold are rejected since the leveling accuracy depends on the arc length (Mannucci et al., 1998).Furthermore, an algorithm to detect and remove degraded observations is performed using a socalled "3σ " outlier test, which is carried out by screening the differences of the leveled carrier phase and the pseudo-range measurements for each arc separately, i.e., d s r = P s r,I − L s r,4 .Thus, the observations exceeding the threshold of 3σ are rejected from the data set.
The next step, the filtering of hourly data, includes the parameter estimation procedures driven by the implemented KF.Finally, the estimated ionospheric parameters are stored and utilized to generate the ionospheric products, for instance, IONEX-formatted files including global VTEC maps.

Results and validation
To asses the quality of the VTEC products generated by the modeling approach described before, two different evaluation methods are considered.The first validation method, called self-consistency analysis, performs a very precise sensitivity analysis using the differences in STEC observations to locally detect temporal and spatial variations around a given GNSS receiver.The second validation procedure shows the estimation quality of the VTEC maps on water surfaces; the results are compared with altimeter data acquired from the Jason-2 mission (Fu and Cazenave, 2001).
For the validation of our results we use the final, i.e., the post-processed, products of the IGS and its IAACs because they are widely accepted as well-established standards.Note that these final products are available with a latency of days or even weeks, whereas our results are evaluated using preprocessing strategies in near real time.We use statistical metrics, namely the RMS value, the error mean value and the standard deviation, to evaluate variations of the VTEC products with respect to reference values derived from the selfconsistency analysis and the Jason-2 altimetry.
In this analysis, the VTEC products are labeled by the following standard convention for the IONEX files, including VTEC maps, as "igsg", "codg", "jplg", "esag" and "upcg",  which are provided by the IGS and its IAACs, namely CODE, JPL, ESOC and UPC.In this sense, the label "dfrg" throughout this work refers to the estimated VTEC maps of the German Geodetic Research Institute -Technical University of Munich (DGFI-TUM) with a temporal resolution of the KF step size set to 5 min, whereas "d1rg" is generated from "dfrg" and comprises VTEC maps with a temporal resolution of 1 h.
A time interval between 11 August 2016 (day of year, DOY 224) and 25 August 2016 (DOY 238) covering 2 weeks of data was considered with the following geomagnetic and ionospheric conditions: the 3 h Kp index data 1 , which quantifies the intensity of the planetary geomagnetic activity, shows low variability with Kp ≤ 5.The sunspot number 2 , which is a good indicator of solar activity, shows a considerably high peak value of 80 for 16 August and the lowest value of 14 for 21 August.The characteristics of these data sets are -as already mentioned -downloaded and preprocessed in near real E. Erdogan et al.: Near real-time estimation of ionosphere VTEC time.The following subsections are dedicated to the results of these assessments.

Self-consistency analysis
The derivation of very accurate absolute STEC values from GNSS measurements may be a challenging procedure since the observations include the DCBs of the receivers and the transmitting satellites.Several research groups have provided GNSS-based solutions regarding the TEC modeling with appropriate approaches for quality assessment; for example, see Rovira-Garcia et al. (2015), Li et al. (2015), Brunini et al. (2011) and Orus et al. (2007).In the context of quality assessment, through a GPS phase-continuous arc, differential STEC, i.e., dSTEC values, can be obtained with an accuracy of less than 0.1 TECU (Feltens et al., 2011).Note that 1 TECU is equivalent to 10 16 electrons m −2 .A test value dSTEC k for assessing the quality of the products at an epoch k can be obtained by (Orus et al., 2007), where dSTEC obs,k is the difference of the GPS geometry-free linear combination at the epoch k with another linear combination computed on the same continuous arc at a reference epoch characterized by the highest elevation angle.The computed dSTEC values from the VTEC maps denoted as dSTEC map,k at the same epoch k and the reference epoch are obtained by multiplying the VTEC values with the elevation-dependent mapping function (Eq.13).
The geographical locations and the names of the receiver stations selected for the dSTEC evaluation are depicted in Fig. 4. The test receivers are chosen globally and located at low and high latitudes, which can reveal the VTEC model accuracy at the regions characterized by varying VTEC activity.The receivers at the sites MKEA, ASPA, BOGT, CHPI, YKRO, DGAR and PIMO are located at middle and low latitudes, whereas DUBO, PENC, URUM, SYOG and MAC1 are established at higher latitudes.As an exemplified analysis, the mean, the standard deviation and the RMS values of daily dSTEC variations are computed using the data from the observation site PIMO as presented in Fig. 5.The numbers within the parentheses shown on the legends give the average values of the corresponding statistical measures for the entire test period, which are also summarized in Fig. 6.The biases of the VTEC maps at the PIMO station show day-to-day variations from 0.1 to 1.2 TECU.The average RMS errors of the dfrg and the d1rg solutions are 2.34 and 2.39 TECU, respectively, which are in close agreement with the RMS values of the analysis centers ranging from 1.99 to 2.72 TECU.Taking into account the daily RMS variations plotted in Fig. 5, the VTEC solutions dfrg and d1rg show larger values only at DOY 230 but smaller deviations for the rest of the test period.In Fig. 6, as a summary of statistical measures, the average mean values and standard deviations as well as the average RMS errors for the entire test period are visualized for each of the 12 receiver sites and the seven VTEC models.As is expected, the dSTEC error for the sites at low latitudes (ASPA, BOGT, YKRO, DGAR, PIMO) are general higher than those obtained at high latitudes.Especially, the ASPA site has rather high errors for each of the VTEC maps in terms of RMS error.The site is very close to the geomagnetic equator, and it can suffer from poor nearby data coverage due to its geographic location on an oceanic island.The result of our solution dfrg has an average bias of 0.04 TECU.The biases with respect to the other maps vary between −0.02 and 0.36 TECU.The average standard deviation for dfrg is 1.78 TECU, which shows a moderate accuracy compared to those of the other maps, which range from 1.66 to 1.96 TECU.The RMS error of dfrg has 1.81 TECU, whereas the RMS errors of the other VTEC solutions vary between 1.70 and 2.00 TECU.A closer look into the average RMS errors for each of the sites reveals that although the presented approach shows a slightly higher deviation only for the ASPA station, the accuracy at the remaining sites is compatible with that of the analysis centers.Moreover, concerning the temporal resolution, the RMS error of dfrg amounts to 1.81 TECU and is thus slightly better than d1rg with 1.85 TECU.However, during high solar activity, a considerable difference might be expected.

Validation by the Jason-2 altimetry
The dual-frequency altimeter onboard the Jason-2 satellite can directly measure in the nadir direction using two different frequencies, which allow the extraction of VTEC data with less effort and without applying any mapping between STEC and VTEC (Dettmering et al., 2011) according to Eq. ( 12).The VTECs provided by Jason-2 are taken into account as reference since they allow the evaluation of the errors of the global VTEC maps over the oceans as well as over regions that exhibit poor estimation quality due to the availability of only a few GNSS measurements.It should be noted that although the Jason-2 satellite provides accurate, direct and independent VTEC data, several studies reported that the measurements are contaminated by an offset of around 3 TECU compared to the GNSS-derived VTEC products (Tseng et al., 2010;Chen et al., 2016).
Before using the altimeter data for the comparisons, a median filter with a window size of 20 s was applied to smooth the data.It is worth mentioning that Jason-2 radar altimetry provides data with a higher spatial and temporal resolution compared to the VTEC maps.Therefore, a linear interpolation in the spatial and time domains was applied to obtain VTEC values for the corresponding time and location of the altimetry observations.The interpolation is performed in the Sun-fixed coordinate system between consecutive epochs.Figure 7 shows exemplified ground tracks of Jason-2 with the associated VTEC values for the entire day of 16 August 2016 (DOY 229).High VTEC variations around the equatorial re-    gions can be clearly seen.This day has the highest sun spot number during the test period.
For the sake of clarity, a selected data set between 00:00 and 01:00 UTC from Fig. 7 is depicted in Fig. 8.The equatorial ionization anomaly (EIA), which is characterized by two crests at both sides of the geomagnetic equator and a trough around the equator (Appleton, 1946), can be seen in Figs.7 and 8 by carefully interpreting the colors along the ground tracks.The VTEC values obtained from Jason-2 and the seven VTEC maps, using the aforementioned interpolation method, are visualized in Fig. 9 for the ground track displayed in Fig. 8.The camelback-shaped EIA in terms of VTEC is illustrated by two regions, including the peak VTEC values between 0.0 and 0.3 h of day (hod).The morphology of the EIA defined by the Jason-2 VTEC is clearly represented by the VTEC maps.The DGFI-TUM solutions are in good agreement with all the other VTEC values.For the For the entire test period, the results of the comparisons in terms of daily mean, standard deviation and RMS values as well as their overall averaged values for the entire test period are illustrated in Fig. 10.The average relative biases of the VTEC maps show variations between 0.5 and −2.0 TECU; our two solutions labeled as dfrg and d1rg both have an average relative bias of −1.6 TECU.The daily standard deviations vary between 3.2 and 5.0 TECU; their averages deviate around 4.0 TECU and are also in accordance with those derived from a recent comparison study by Hernández-Pajares et al. (2016) stating that the daily deviations for the results of the IAACs can vary from a few TECU to 10 TECU.Our solutions have an average RMS error of 4.7 TECU, which does not exceed that computed from the other analysis centers ranging between 4.0 and 4.7 TECU.A detailed look into the daily RMS values indicates that the estimated VTEC values of our solutions are generally in good agreement with the results of the analysis centers.

Conclusions and future improvements
A near real-time processing framework that is capable of automated data downloading, data preprocessing, Kalman fil-tering and formatted product generation is presented to provide VTEC maps as well as satellite and receiver DCBs of GPS and GLONASS in near real time.The B-spline representation of global VTEC is incorporated into the Kalman filter procedure.The filter was also extended to integrate the equality constraint equations comprising the spherical and DCB-related restrictions.Coefficients of the B-spline model and the DCBs, which constitute the unknown parameters, are recursively estimated by exploiting hourly GNSS observations acquired from the IGS data centers with 1 h latency.The ionosphere observable is derived from raw GNSS code and phase measurements using the geometry-free linear combinations.
The validation of the proposed approach is carried out using GNSS data downloaded in near real time covering a time span of 2 weeks.To summarize, according to the self-consistency analysis, an RMS value of 1.81 TECU was found.The four IAACs, CODE, JPL, ESOC and UPC, as well as the IGS combination product exhibit comparable RMS errors between 1.70 and 2.00 TECU.Moreover, the Jason-2 validation shows that the RMS error achieved by the proposed method fits well with the results of the IAACs.Considering the comparisons, specific to the test period it might be concluded that the estimated VTEC products using the presented near real-time strategy shows promising initial results in terms of accuracy and overall agreement with the post-processed final products of IGS and its analysis centers, which are publicly available with several days of latency.

Figure 1 .
Figure1.Global VTEC representation in a solar geomagnetic coordinate system using different resolution levels; (a)J 1 = 3, J 2 = 1; (b) J 1 = 4, J 2 = 2; (c) J 1 = 4, J 2 = 3; (d) J 1 = 5, J 2 = 3.The circles and triangles represent the data locations for GLONASS and GPS, respectively, and the colors indicate the corresponding VTEC magnitudes related to the ionospheric pierce points (IPP) obtained from a reference map.In all the four panels at the bottom and on the left side the K J 2 trigonometric B-spline functions T 2
11) i.e., L s r,4 ∈ {y GPS , y GLO }.The DCB values b r,GPS and b r,GLO stand for the GPS and GLONASS receiver biases.Furthermore, b s GPS and b s GLO represent the unknown DCB values of the respective satellites.The mapping function m(z) depending on the zenith angle z projects STEC into the vertical by VTEC = 1 m(z) STEC .

Figure 3 .
Figure 3. Overall scheme for the global VTEC estimation and the product generation.

Figure 4 .
Figure 4.The geographical locations and the identifiers of the receiver sites used in the dSTEC analysis.

Figure 5 .
Figure 5. Results of the statistical evaluations presenting the differences between the observed and computed dSTEC values at the station PIMO.The investigation covers the days DOY 224 to DOY 238, 2016: (upper panel) mean of deviations, (middle panel) standard deviations and (bottom panel) RMS values of errors in terms of TECU.

Figure 6 .
Figure 6.Results of the statistical evaluations presenting the differences between the observed and computed dSTEC values, which cover the days between DOY 224 and DOY 238, 2016: (upper panel) mean of deviations, (middle panel) standard deviations and (bottom panel) RMS values of deviation in terms of TECU.

Figure 7 .
Figure 7. Ground tracks of the Jason-2 altimetry mission for 16 August 2016 (DOY 229).The colors show the magnitude of VTEC as acquired from the satellite measurement system.

Figure 8 .
Figure 8. Ground track of the Jason-2 altimetry satellite between 00:00 and 01:00 UTC on DOY 229, 2016.The colors show the magnitude of VTEC in TECU acquired from the satellite.

Figure 10 .
Figure 10.Comparison of VTEC values acquired from the analysis centers and the DGFI-TUM solutions with Jason-2 altimetry VTEC data between DOY 224 and DOY 238, 2016; (upper panel) mean deviations, (middle panel) standard deviations and (bottom panel) RMS values of deviations in terms of TECU.