Adaptive Modeling of the Global Ionosphere Vertical Total Electron Content

The Kalman filter (KF) is widely applied in (ultra) rapid and (near) real-time ionosphere modeling to meet the demand on ionosphere products required in many applications extending from navigation and positioning to monitoring space weather events and naturals disasters. The requirement of a prior definition of the stochastic models attached to the measurements and the dynamic models of the KF is a drawback associated with its standard implementation since model uncertainties can exhibit temporal variations or the time span of a given test data set would not be large enough. Adaptive methods can mitigate these problems by tuning the stochastic model parameters during the filter run-time. Accordingly, one of the primary objectives of our study is to apply an adaptive KF based on variance component estimation to compute the global Vertical Total Electron Content (VTEC) of the ionosphere by assimilating different ionospheric GNSS measurements. Secondly, the derived VTEC representation is based on a series expansion in terms of compactly supported B-spline functions. We highlight the morphological similarity of the spatial distributions and the magnitudes between VTEC values and the corresponding estimated B-spline coefficients. This similarity allows for deducing physical interpretations from the coefficients. In this context, an empirical adaptive model to account for the dynamic model uncertainties, representing the temporal variations of VTEC errors, is developed in this work according to the structure of B-spline coefficients. For the validation, the differential slant total electron content (dSTEC) analysis and a comparison with Jason-2/3 altimetry data are performed. Assessments show that the quality of the VTEC products derived by the presented algorithm is in good agreement, or even more accurate, with the products provided by IGS ionosphere analysis centers within the selected periods in 2015 and 2017. Furthermore, we show that the presented approach can be applied to different ionosphere conditions ranging from very high to low solar activity without concerning time-variable model uncertainties, including measurement error and process noise of the KF because the associated covariance matrices are computed in a self-adaptive manner during run-time.


Introduction
The selection of an appropriate parameter estimation strategy, which allows for handling of the large data sets from various space geodetic observation techniques, e.g., the continuously operating IGS network of GNSS receivers, is vital for (ultra) rapid and (near) real-time VTEC modeling. Recursive filtering methods, e.g., the Kalman Filter [1] and the Ensemble Kalman Filter [2], provide mathematical tools in terms of assimilating new data immediately once they are available [3,4].
These conventional filters are generally extended by adaptive approaches to cope with time-varying model uncertainties in many applications. In this sense, the focus of this paper is on the application of the adaptive discrete Kalman Filter (KF) to ultra-rapid VTEC modeling based on B-splines. The presented approach is a comprehensive extension to the study by Erdogan et al. [5] employing a conventional KF.
The discrete Kalman filter (KF) relies on consecutively performed steps, namely the time update and the measurement update; the former represents the propagation of the unknown parameters with time and the latter step is responsible for the incorporation of new allocated measurements into the filter [3]. A stochastic model attached to the deterministic models describes the corresponding covariance matrices of process and measurement noise. The conventional KF is very sensitive to these matrices which are required to be introduced beforehand. The common practice to manually define the parameters of these matrices is to perform multiple tests by executing the filter for different values of the parameters; for further details, we refer to Maybeck [6]. However, the data used in tests may not be appropriate to define these matrices. Furthermore, the covariance matrices may change with time. An inappropriate definition of these matrices can lead to an estimation of poor quality, or even worse, the filter may diverge. A self-learning filter which is capable of adapting itself to changes during run-time can cope with such situations [4]. Accordingly, several approaches, classified under the name of adaptive modeling, were proposed after Kalman's seminal study from 1960 to tune the parameters in run-time for achieving optimal results and avoiding a filter divergence.
According to Hide et al. [7], the approaches for adaptive modeling can be categorized as (1) covariance scaling, (2) multiple-model adaptive estimation (MMAE) and (3) stochastic modeling. The first method means a simple and effective algorithm in which the predicted or process noise covariance matrix of the KF is multiplied by a factor which can be deduced from the innovations or the residuals of the filter [8][9][10]. The MMAE was introduced by Magill [11] and relies on the combination of multiple outputs derived from multiple KF instances running in parallel with different values of measurement and process noise parameters. High computational cost is the drawback of the MMAE method. Furthermore, methods based on a stochastic modeling attempt to estimate the overall covariance matrices using the innovation or residual sequences collected and stored at previous steps of the filter, see e.g., [12]. This approach makes use of data sequences collected in a time window of past epochs and suffer mostly from the fact that the size of innovations, residuals or state vectors of the filter at the consecutive epochs can change over time [13]. Furthermore, the method of variance component estimation (VCE), which goes back to the study by Helmert [14], was extensively investigated using different types of estimators to obtain a realistic stochastic model of measurement uncertainties in the sense of batch parameter estimation. We refer to, e.g., [15][16][17][18] and references therein for a comprehensive review. In recent decades, many studies were carried out in an attempt to incorporate the VCE approach into the KF to achieve a recursive and an adaptive estimator. A VCE-based adaptive KF does not require the storage of data from previous epochs of the filter and may be categorized as a covariance scaling approach. For example, Yang and Xu [13] introduced an adaptive robust Kalman filter resisting to outliers as well as determining an adaptive factor using the ratio of variance components between the measurements and the predicted state vector computed by referring to the study by [19]. Hu et al. [9] developed an adaptive estimator derived from innovations. Gao et al. [20] presented an adaptive KF based on the Helmert VCE method for the positioning and the attitude determination by integrating measurements from the Multi-Constellation GNSS and the Inertial Navigation System. Chang et al. [21] incorporated the least-squares VCE approach into a KF to estimate the epoch differenced ionospheric delay for repairing cycle slips of GNSS signals. In the framework of GNSS meteorology, Yang et al. [22] presented an adaptive KF algorithm to tune the process noise in real time by inserting the least-squares VCE method into the filter to compute precipitable water vapor content derived from estimated zenith tropospheric delays of GNSS signals.
VCE allows for the determination of proper weights for different types of observation groups with varying model uncertainties, which can exhibit limited knowledge of noise characteristics [23].
A covariance matrix of an observation type can generally be decomposed into two parts, namely a known co-factor matrix and an unknown scaling factor, i.e., the variance component. The known matrix refers to a weighting matrix determining the contributions of observations according to pre-defined quality criteria. For instance, in this study it is defined according to the elevation angle of GNSS observations. Run-time computations of the unknown variance components in a Kalman filter for each observation group allows the implementation of a computationally efficient adaptive filtering approach. When the VCE approach is chosen, it decreases the number of unknown parameters compared to the stochastic modeling approach attempting to estimate an entire covariance matrix which generally leads to a less stable filter implementation due to an increased number of unknowns as well as lack of enough observations. In this paper, we use the VCE approach for the adaptive Kalman filtering implementation to separately handle different observation groups. The estimation of variance components for each type of observations acquired from GPS and GLONASS is iteratively carried out and can be driven according to the Maximum-Likelihood method or Förstner's iterative method, see, e.g., [16,18,[23][24][25].
Contrary to physical models, empirical models are mostly data driven approaches and do not rely on physical equations and quantities. The selection of a proper empirical model which allows analyzing local structures and deducing physical interpretations from the model parameters is crucial. Ionosphere models relying on localized basis functions have been increasingly gaining popularity, for instance, applications of B-spline functions [26][27][28][29] and Slepian functions [30]. The B-spline representation of VTEC provides an empirical approximation for evaluating temporal and spatial variations. Whereas the VTEC products provided by the Ionosphere Associated Analysis Centers (IAAC) of the International GNSS Service (IGS) are mostly based on spherical harmonics (SH), i.e., basis function of global support, DGFI-TUM aims at VTEC representations using series expansions in terms of B-spline functions. B-splines are an appropriate mathematical tool for handling the heterogeneous input data distribution including data gaps [31] due to their compact support. Moreover, they allow the generation of a multi-scale signal analysis of VTEC [26,[32][33][34] and can easily be extended for multi-dimensional model definitions, e.g., for 3D/4D electron density modeling [27,35,36].
By exploiting the advantages of B-splines, we show that the estimated B-spline coefficients resemble the global VTEC structure in terms of shape and magnitude, which paves the way of assigning a physical meaning to the B-spline coefficients. For instance, we found a very high correlation between a time series of mean B-spline coefficients and corresponding mean VTEC values providing a valuable indicator about temporal VTEC activity as well as space weather events. The physical interpretation of the coefficients leads to extending the definition of the stochastic model attached to the dynamic model, which is also called the prediction model of the KF. In this context, an exponential empirical model taking the mean value of the B-spline coefficients and their relative variations into account is introduced to define the uncertainty of the prediction model. Moreover, the compact support property of the B-spline representation provides additional flexibility to locally tune the sensitivity of the filter to observations without causing a global effect. This is unlikely for models making use of globally defined basis functions e.g., spherical harmonics. The empirical model is extended to consider the number of observations to keep the filter more sensitive to the observations at regions including very high data density.
In summary, the main goal of the present study is to introduce an adaptive KF for VTEC modeling using a B-splines series expansion. The measurement covariance matrices of the GNSS observations are adaptively obtained using the VCE approach. We highlight structural similarities between the spatial distributions of VTEC over the globe and the estimated B-spline coefficients that allow for deducing physical interpretations from the coefficients. Accordingly, we introduce an adaptive empirical model for the process noise covariance matrix of the B-spline coefficients regarding the temporal and spatial variation of VTEC. The presented adaptive filtering approach can significantly reduce the efforts to define or to set up the model uncertainties. Moreover, exploiting the localizing property of the B-spline basis functions allows for tuning the filter sensitivity to observations in the vicinity of a point without causing a significant global effect, which means a huge advantage w.r.t. methods relying on globally defined basis functions such as spherical harmonics. The approach is applied to ultra-rapid VTEC modeling by assimilating carrier phase and pseudo-range observations from GPS and GLONASS, which can be extended for additional GNSS constellations such as GALILEO or measurement techniques such as the DORIS (Doppler Orbitography and Radiopositioning Integrated by Satellite) technique. The presented adaptive algorithm was recently employed in [37] to generate high-resolution VTEC maps of Europe by means of a series expansion in B-spline functions.
The paper is outlined as follows. The first section shows the B-spline representation for global VTEC modeling and the definition of the selected coordinate system. In Section 3, the extraction of ionosphere observations from raw hourly GNSS data is explained. Section 4 gives the background for the sequential estimation algorithm using the adaptive KF. This section comprises the introduction of the measurement and prediction models of the filter, the fundamentals of the adaptive Kalman filtering using variance components, and the handling of the model constraints in the filter. Section 5 explains the validation techniques, accuracy of the presented approach, and discussions. Finally, Section 6 provides the conclusion and future work.

Global VTEC Representation Based On B-Splines
In this study, the B-spline representation of the global VTEC distribution reads VTEC(ϕ, λ, t) = increased further. A higher resolution B-spline representation with level values J 1 = 5 and J 2 = 3 in comparison to a spherical harmonic representation and the VTEC products of the IGS analysis centers (CODE and UPC) were studied by [32]. By taking into account the trade-off between the resolution level and the computational burden when running the filter in an operational mode, B-spline level values of J 1 = 5 and J 2 = 3 were selected in this study.

Coordinate System
The Earth's magnetic field plays a crucial role on dynamics, structure and formation of the upper atmosphere and space weather phenomena such as ionospheric plasma motion, polar lights, the electron density distribution and the equatorial ionization anomaly. A proper selection of the coordinate system can become significant in terms of the strength of a dominating driver. For instance, at the altitudes close to the Earth's surface, the geomagnetic field is stronger than the magnetic disturbances derived by solar winds, so that, a coordinate system aligned with the Earth's magnetic dipole axis would be more appropriate for ionospheric phenomena [38]. Solar magnetic coordinates, which take into account the shift of the dipole axis from the Earth's rotation axis, are usually adapted for ionospheric electron content modeling (e.g., Mannucci et al. [39]). In opposite to this simple dipole approximation, a better representation of the Earth's magnetic field can be obtained by considering the non-dipole features of the field which can show significant deviation from the center dipole approximation at the ionospheric altitudes. Corrected geomagnetic coordinates [40,41], modified apex coordinates and quasi-dipole (QD) coordinates [42] are some examples which come at the expense of, e.g., a high computational load, complexity and non-orthogonality. For example, the apex coordinates are employed in physics-based models simulating the thermospheric and ionospheric behavior, such as the Thermosphere Ionosphere Electrodynamics General Circulation Model (TIEGCM) [43][44][45]. Among the ionosphere modeling community, the modified dip latitude (MODIP) coordinates [46], which takes the Earth's magnetic field into account to compute the latitude of a given point in terms of the dip latitude, are also in use. Azpilicueta et al. [47] applied the MODIP coordinates to the VTEC representation using spherical harmonics and declared a significant improvement in their solutions compared to the VTEC solution obtained using geographic latitude.
Since the Sun is the primary driver of the photo-ionization process, it causes spatio-temporal variations in the ionosphere. In VTEC modeling, keeping the VTEC structure aligned with the position of the Sun leads to a Sun-fixed coordinate system definition. This mitigates the effect of the Earth's diurnal motion that leads to a much slower temporal variation of ionospheric VTEC (see e.g., Mannucci et al. [39], Komjathy and Langley [48]).
In this study, in order to consider the effects of the magnetic field and the Sun on the formation of VTEC structure, a solar magnetic coordinate system (GSM) (e.g., [49]) is used. The z-axis of the GSM coordinates aligns with the centered dipole axis and points to the Northern Hemisphere. The x-axis is defined such that the x-z plane contains the Sun-Earth line and the y-axis is perpendicular to the Sun-Earth line.

GNSS Ionospheric Observables
Due to the dispersive characteristics of the ionosphere, electromagnetic signals transmitted by the GNSS satellites are refracted while traveling through the medium. The rate of the ionospheric refraction varies with respect to the frequencies of the transmitted signals. The geometry-free linear combination of the code pseudo-ranges observations P s r, f 1 and P s r, f 2 as well as the carrier-phases observations Φ s r, f 1 and Φ s r, f 2 , acquired from GPS and GLONASS signals at the distinct carrier frequencies f 1 and f 2 , lead to the ionospheric observables where P s I,r and L s I,r are the pseudo-range and the carrier-phase ionosphere combinations, respectively [50]. B r and B s refer to the receiver and the satellite inter-frequency biases (IFB) of the carrier phase observations, b r and b s are the pseudo-range IFB of the receiver and the satellite, respectively, which are also commonly called differential code biases (DCB). The combined ambiguity bias of the carrier-phase observations is denoted by B s A,r . The terms e P s I,r and e L s I,r stand for the measurement errors. Furthermore, considering P s I,r and L s I,r are given in units of meters, α s r is a frequency-dependent factor defined as α s r = 40.3 16 . The pseudo-range ionosphere combination (2a) is rather noisy but unambiguous. Although the carrier-phase ionosphere combination is around two orders of magnitude more precise, it is biased by an unknown integer number of cycles [50,51]. A widely common approach for mitigating the ambiguity and exploiting the precision of the carrier-phase observation is the so-called leveling bias computation [39,52]. The leveling bias C s r includes the IFBs and the carrier-phase ambiguity term. It can be determined by the weighted averaging of the epoch-wise differencing of L s I,r and P s I,r and reads where w i is the weight of the i th observation referring to the precision of the differenced observations on a phase continuous arc (see e.g., [39]). The term e C s r is the error of the leveling bias. The value for the weight w i can be determined by taking into account, e.g., the elevation angle of the observations, the receiver tracking modes and the signal power. Prior to the computation of the leveling bias, a two step cycle slip detection method, which is composed of a double differences of ionospheric phase observations and the Melbourne-Wübbena combination, is applied for the detection of jumps and the construction of phase continuous arcs between receivers and satellite pairs [53][54][55].
Once the leveling bias is obtained, the leveled carrier-phase ionospheric observations L s I,r along a continuous arc are computed by where e L s I,r refer to the error of the leveled carrier-phase ionosphere combination which is dominated by the error e C s r of the leveling bias. The accuracy and the precision of the leveling bias in Equation (3) depends on different factors which are primarily measurement noise and multipath effects. The measurement noise on carrier-phase measurements can be neglected, since those effects are much smaller compared to the ones on the code measurements. Ciraolo et al. [50] showed that the multipath effect on code measurements is a primary error source that cannot totally be removed in Equation (3) by the averaging procedure and can vary for different receivers and antennas types. Besides, depending on the environmental conditions around the receivers, the IFBs of code measurements can present intra-day variations that violate the long term stability assumption on the IFBs, and consequently result in a systematic error in the leveling bias.
To reduce the effects of the systematic errors, the following measures are considered in the computation of the leveled carrier-phase observable:

•
Observations with an elevation angle of less than 10 • were rejected from the data to avoid contributions from likely very noisy measurements.
• Since the quality of the averaging procedure in Equation (3) relies on the length of a given phase continuous arc [39], a minimum threshold of 30 min arc length is applied to mitigate the pseudo-range noise.

•
To mitigate the errors due to multi-path effects on pseudo-range measurements, which are generally inversely proportional to the satellite elevation angle, only the observations with elevation angles larger than 20 • are included in the determination of the leveling bias (3).

•
Moreover, an elevation dependent weighting function, e.g., w i = sin(e i ), is introduced to leverage the influence of more precise observations. e i means the elevation angle of the i th observation along the arc. • A pre-processing algorithm was developed for the recursive processing of GNSS data acquired as hourly data blocks broadcasted, e.g., by the IGS global data centers. To reach a maximum number of observations and a high pre-processing quality, data blocks with a moving window of 3 h length are considered. For instance, to apply the pre-processing procedures to a new acquired data set between t k and t k + 1h, a data set with a window length of 3 h extending from t k − 2h to t k + 1h is considered for a more accurate computation of the leveling bias.
Assuming that the code and carrier-phase ionosphere observations are uncorrelated throughout a continuous arc, the variance of the leveled carrier-phase observations can be derived by applying the law of error propagation to Equation (4). Consequently, the variance σ 2 where k = 1, . . . , N is the time index of the observations along an arc with a total number of N observations. l b and l e , respectively, refer to the beginning and the end points on the arc that bind the data used for the determination of the leveling bias in Equation (3)

Model Definition
The discrete linear system of equations describing a dynamical process is given by where k is the time stamp, F k is the transition matrix, β k the vector of the unknown parameters, y k is the vector of the observations and X k is the corresponding design matrix. The measurement error vector e k and the vector of the process noise w k are normal distributed, i.e., w k ∼ N(0, Σ w,k ) and e k ∼ N(0, Σ y,k ) and fulfill the assumptions and The covariance matrices referring to the vectors e k and w k are automatically computed by the introduced algorithms in a self-learning manner. Accordingly, the adaptive algorithm defining the covariance matrix of the measurement error vector is explained in Section 4.4. The developed time-varying adaptive model of the process noise covariance matrix for the ionospheric target parameters is presented and discussed in Section 4.3.

Measurement Model
The measurement equations following from Equation (4) are re-written as where y GPS s r ,k and y GLO s r ,k are the leveled ionospheric observations in TECU at time stamp k for the GPS and GLONASS constellations. The quantity STEC in Equation (4) at time t k is represented by the so-called Single Layer Model (SLM) (see e.g., [56]) as STEC k = m(z s r,k ) · VTEC k where m(·) stands for the Modified Single Layer Mapping Function depending on the satellite elevation angle z s r,k [56,57]. Moreover, VTEC k at time t k is given by Equation (1). The values b r,GPS , b r,GLO and b s GPS , b s GLO stand for the receiver and the satellite DCBs in nanosecond for GPS and GLONASS. The speed of light is denoted as c in m/s. The state vector β k of the linear models given by the Equations (6) and (7) consists of the sub-vector Consequently, the design matrix X k is given as Herein, X d GPS ,k and X d GLO ,k are the design sub-matrices of GPS and GLONASS observations in connection with the vector d k of unknown B-spline coefficients. The design sub-matrices for the GPS and GLONASS receiver DCBs are X b GPS ,k and X b GLO ,k . In the same way, the sub-matrices X b GPS ,k and X b GLO ,k refer to the design matrices of the satellite DCBs.
Assuming the vectors y GPS,k , and y GLO,k are uncorrelated, the measurement covariance matrix Σ y,k consisting of the covariance matrices Σ y GPS ,k and Σ y GLO ,k for GPS and GLONASS reads This study aims on the computation of the measurement covariance matrix Σ y,k in run-time by employing an appropriate adaptive estimation method. For problems with a large number of unknowns and/or measurements, the computation of each component in the measurement covariance matrix results in very high complexity and computational burden. Additionally, an unstable filter implementation would emerge due to a lack of a sufficient number of observations. To cope with these problems, a common simplified approach is to assume that the covariance matrix is known up to a scaling factor. Taking VCE into account, this factor can be handled as an unknown variance component. Accordingly, the measurement covariance matrix Σ y j ,k for each observation technique j ∈ {1, . . . , p} with the associated variance component σ 2 y j ,k is expressed by (see [19] ) where P y j ,k is the appropriate weight matrix. Consequently, considering that y 1 = y GPS , σ 2 y 1 ,k = σ 2 y GPS ,k and y 2 = y GLO , σ 2 y 2 ,k = σ 2 y GLO ,k , Equation (13) can be written as where σ 2 y GPS ,k and σ 2 y GLO ,k are the unknown variance components for GPS and GLONASS. The matrices P y GPS ,k and P y GLO,k are given positive definite matrices which can be determined to ensure more contributions from the precise observations by introducing high relative weights. In this sense, the weight matrix can refer to numerous quality factors for GNSS observations, such as the satellite elevation angle, the receiver-antenna configuration and the signal strength. A common and simple strategy is to adapt an elevation-dependent weighting scheme for each individual observation. Assuming that the observations are uncorrelated, the weight matrix P y j ,k becomes a diagonal matrix and its diagonal element p nn can be defined by where z m is the zenith angle of the mth measurement. The term σ 2 L s I,r ,k refers to the prior observation variance described by Equation (5).

Prediction Model
Since VTEC is a time varying phenomenon, a proper model is required to represent its propagation between consecutive epochs t k−1 and t k in accordance with Equation (6).
In the current implementation of the approach, VTEC is represented in the GSM, and thus, leading to slow temporal variations in the corresponding B-spline coefficients. Accordingly, a simple model based on a random walk process is set to represent the temporal variations. Moreover, it should be noted that the unknown system biases, precisely the GNSS DCBs are also quite stable over the course of time, which allows the use of simple prediction models, too. Consequently, a random walk process is introduced as a prediction model for these biases. The entire prediction model with F k = I of the state vector reads, therefore, following from Equation (6) where w k is the process noise vector with the covariance matrix Σ w,k defined as where Σ w d ,k is the diagonal matrix consisting of the components σ 2 d i ,k with i = {1, . . . , K J 1 · K J 2 } corresponding to the variance of the random walk process for each B-spline coefficients. The quantities Process noise variances should be set properly, since these values effect the filter performance. The improper selection of the process noise could lead to a divergence of the filter. Moreover, the process noise plays an important role to tune the sensitivity of the filter to the measurements. The higher the value of the process noise the more sensitive is the filter not only to the measurements but also to, e.g., outliers and undetected biases in the measurements. For instance, increasing the process noise results in larger values of the predicted covariance matrix Σ − β,k in Equation (31) which leads to a relatively high weighting of the measurements.
Considering the changing environmental conditions, such as seasonal solar activity or solar storms, the error of the prediction model for the B-spline coefficients is time dependent. Therefore, each process noise variance component σ 2 d i ,k of Σ w d ,k in Equation (18)   In this regard, the mean VTEC value at a given time provides an appropriate indicator representing the overall VTEC activity sensitive to seasonal, monthly and daily variations (see e.g., [56]). Moreover, abnormal variations due to high solar activities can be detected from the mean VTEC value (see e.g., [58] in case of coronal mass ejections). The global mean VTEC values of the IGS final GIMs during the year 2015 with a temporal sampling of 2 h are plotted (red dots) in Figure 2. In the same figure, the corresponding mean values of the B-spline coefficients are illustrated by blue dots. There is a very high correlation between the mean VTEC values and the corresponding mean B-spline coefficients, namely ρ = 0.999.
In the general case, the variance of a random walk process can be defined as where ∆t is the length of the time step between two consecutive epochs (i.e., ∆t = t k − t k−1 ) andq k is the variance rate referring to the change of the variance with time [39,59]. Considering that the B-spline coefficients resemble the structure of the VTEC signal, the variance rateq d i ,k of the process noise σ 2 d i ,k for the coefficient d i,k is defined by the producṫ The coefficient C 0,k is given as where m s is a constant value determined experimentally. The mean absolute valued k of the absolute value of the B-spline coefficients for global VTEC modeling is given as C 1,d i,k and C 2,d i,k are coefficients to tune the filter sensitivity according to spatial variations of the B-spline coefficients and data distribution. The model error of the random walk approach for a B-spline coefficient is proportional to the ionospheric activity, since the higher the VTEC values, the larger the associated B-spline coefficients as well as the higher is the ionospheric activity. In this context, the coefficients C 1,d i,k are introduced to represent the relative variations of the process noise in magnitude between the B-spline coefficients. The coefficients C 1,d i,k are defined by using an exponential empirical model as where the e-function allows the generation of bounded values approaching to a limit value. From the Equations (19), (20) and (23), a relatively high process noise is assigned to a coefficient which has a large absolute value. Additionally, a further control on the filter behavior with respect to the data distribution can locally be accomplished using the B-spline representation. The sensitivity of the filter to ionosphere observations can be tuned regionally. For instance, the European and the North American regions have a very dense data distribution. By a proper setting of the process noise of the B-spline coefficients corresponding to these regions, the filter can be forced to be more sensitive to these observations without causing a global effect. To this end, the coefficient C 2,d i,k is introduced and defined by where m w is a constant value determined experimentally. The term N d i ,k stands for the number of measurements supporting the B-spline coefficient d i and it can be obtained, e.g., by counting the total number of non-zero elements through the corresponding column of the design matrix in Equation (12). Furthermore, the term N d total ,k is the total number of observations at the epoch t k . Exemplary illustrations for the maps of the coefficients, C 1,d i ,k and C 2,d i ,k with C 0,k = 0.0006 are shown in Figure where the constants C b r,GPS , C b r,GLO , C b s GPS and C b s GLO are defined by conducting tests.

Adaptive Filtering Using Method of Variance-Components (VC) Estimation
Kalman [1] proposed a recursive solution to the problem defined in the Equations (6) and (7). Following this pioneering study, various authors have re-interpreted the filtering problem from different perspectives leading to the same set of KF equations, for instance, the maximum likelihood estimation [60], as a sub-class of the Bayesian estimation [61] and the minimum variance estimation (e.g., Jazwinski [62]). The prediction and correction equations that are solutions of the filtering problem read where K k is the gain matrix and the index "· − " refers to "predicted" whereas " · " stands for "estimated" [3,59,62]. The Equations (28) and (29) of the measurement update step can be re-written as by taking into account the matrix identities g., Koch [17], Koch [18] and Jazwinski [62]).
Observations from different techniques generally exhibit different characteristics in terms of precision and accuracy. Therefore, the determination of the measurement covariance matrix should be handled carefully to provide a reliable estimation of the unknown parameters. Assuming that intra-technique observations are uncorrelated, Equation (7) can be partitioned into blocks for each of the related techniques (see e.g., [23,26]). This leads to the solution with where each group of observation vector y y j ,k (GPS and GLONASS observations in this study) is defined by the index j ∈ {1, . . . , p}. Each covariance matrix Σ y j ,k of observations of the group j is replaced by the expression given in Equation (14) to estimate the associated variance components. Accordingly, the overall set of adaptive filtering equations is given by with Herein the unknown variance component σ 2 y j ,k for each group j is determined via VCE as σ 2 y j ,k = e T y j ,k P y j ,k e y j ,k r y j ,k where e y j ,k = y j,k − X y j,k β k is the residual vector for the observation group j and r y j ,k is the corresponding partial redundancy defined as r y j ,k = n j − trace N y j ,k N −1 k (38) where N k is the overall normal equation matrix, i.e., N −1 k = Σ β,k and n j is the total number of observations of group j. The normal equation matrix N y j ,k for the observation group j is given by An iterative estimation of the unknown variance components σ 2 y j ,k is carried out until a point of convergence is reached. The iteration starts from the approximate values σ 2,o y j ,k of the variance components. The estimated variance components at the previous epoch t k−1 of the filter are used as initial values at the current epoch t k , i.e., σ 2,o y j ,k = σ 2 y j ,k−1 . At each iteration step, the Equations (35)-(37) are re-executed for the updated values of the variance components.

Constrained Filtering
Numerous ways to incorporate constraints into the KF have been studied. For a comprehensive review the reader is referred to [63] and references therein. Here, the focus is on perfect measurements and the estimate projection methods. The main distinction between the two methods is that the former one is performed during the measurement update step of the filter whereas the second method uses the constraints following the measurement update step. The estimate projection method was employed by [5] in the context of ionosphere VTEC modeling relying on a standard KF implementation referring to the Equations (28) and (29). Contrarily, in this study, the "perfect measurement" method is carried out to ensure a more stable KF implementation.
Although the alternative implementation of KF given by the Equations (35) and (36) offers a faster execution time and a lower memory consumption via the normal equation representation for large problems, i.e., for a design matrix X mn with m n, the filter becomes more vulnerable to numerical problems such as ill-conditioning and rank deficiency. For instance, data gaps and an inhomogeneous GNSS data distribution result in unobserved unknown parameters. Moreover, the observation equations have a linear dependency problem due to the terms including satellite DCBs. However, the perfect measurement method increases the numerical stability of the filter, since the method considers the constraints as additional information exploited during the measurement update step of the filter. It threats constraints as fictitious observations and augments the measurement model as where the augmented measurement vector, the design matrix and the error vector are, respectively, given by the new augmented measurement covariance matrix reads where Σ y c ,k is a zero matrix by definition, since the constraints are handled as perfect measurements. However, in order to avoid singularity in practice, Σ y c ,k is generally replaced by a matrix of very small numbers, i.e., Σ y c ,k = σ 2 c I [64]. In this study, σ 2 c is set to 10 −8 which is defined as small as possible by considering the computational precision. The equality constraint equations can be handled in the measurement model described by the Equations (35) and (36) via considering them as a block of measurements with a known covariance matrix (see, e.g., [17,65]). In the scope of this study, the constraint equations introduced for preserving the spherical geometry read where X c d ,k is a known matrix [28]. Furthermore, a zero mean constraint, which is usually adapted in the parameter estimation methods employed by the IAACs, is introduced for the DCBs to vanish the linear dependency in the observation equations. It is given as where X c b GPS ,k and X c b GLO ,k are known matrices. A simplified notation for a general definition of the equality constrained equations is defined as From the Equations (43) and (44) the design matrix X y c ,k of the constraint equations for the current problem reads 46) and the vector y c,k is given by y c,k = 0.

Filter Settings
Analysis centers rely on different sets of GNSS stations selected according to the requirements for their modeling approaches and the chosen performance criteria. In the presented approach, the list of available GNSS stations is downloaded and updated every hour. The selected stations are mostly from the IGS network, which is extended by a few additional stations from the UNAVCO and the EUREF network to obtain a better geographical coverage.
Data from the GNSS receivers is based on the hourly data blocks provided from GNSS data centers with a latency of about 1 h. The temporal sampling of the GNSS data is mostly 30 s or less. This allows for running the filter theoretically with a step size of 30 s. However, high temporal sampling of the filter results in an increased computational burden and a storage requirement, which can lead to difficulties for efficiently managing the resources when the model runs in operational mode. Moreover, VTEC varies slowly in time in a Sun-fixed coordinate system during low and moderate ionospheric activity. However, a very long filter step size should also be avoided since the filter prediction error generally grows in line with the filter step size. Accordingly, the filter step size is set to 5 min by taking the computational load and the Kalman Filter prediction error into account. The estimated target parameters, related auxiliary data, and the VTEC products generated from the estimated parameters are stored every 10 min by considering the storage requirements.
The convergence rate of the filter depends on the accuracy of the initial values of the unknown ionospheric target parameters, their initial covariance matrix and the size of the filter state vector. The initial values of the B-spline coefficients are computed using the IGS analysis center products to decrease the convergence rate. The analysis of the estimated covariance matrices and the residuals over time shows that the filter state vector is stabilized around 1 h (12 epochs with a filter step size of 5 min) after the initialization of the filter.

Validation Methods and Data Sets
The quality of the VTEC products derived by executing the developed adaptive filtering approach is assessed using (1) the dSTEC analysis and (2) altimeter VTEC comparisons. The dSTEC analysis is a very precise and appropriate method for the validation of VTEC products on areas including GNSS receivers. The method has been widely preferred for assessments of ionospheric product quality, for example, see [66][67][68][69].
The comparison of the results with altimeter data acquired from the Jason-2/3 missions allows for a validation of the VTEC maps over oceans. Although the VTEC values obtained from the altimetry missions are noisy and biased, the method provides an appropriate external validation approach for qualifying VTEC values over the oceans, and it is one of the fundamental approaches for the comparison of VTEC products provided by the IGS analysis centers, see e.g., [70,71].
The VTEC products in IONEX format with the labels "igsg", "uqrg", "esag", "jplg" and "codg" provided by the IAACs, which are respectively IGS, UPC, ESA, JPL and CODE, are used for the validation of the results. The VTEC maps of IGS, ESA and JPL refer to their final products with a temporal sampling of 2 h, whereas the product of CODE has a sampling interval of 1 h. Although UPC has final VTEC products, we consider their rapid solution labeled "uqrg" which has a very high accuracy and a temporal resolution of 15 min and is one of the best products within the IGS analysis centers (see, e.g., [70]). The estimated VTEC maps generated using the developed method will be labeled "othg" in the remainder of the paper. The letters "o", "t", "h" and "g" stands for the following expressions; project "O"PTIMAP, temporal resolution of 10 ("t"en) minutes, "h"igh resolution with J 1 = 5, J 2 = 3 and "G"lobal VTEC Maps, respectively. For further information on product generation at different resolution levels and the corresponding label definitions, we refer to [32,37].
Two data sets referring to different time intervals are selected for the VTEC product generation and validation. The first data set comprises GNSS data collected during a considerably high ionospheric activity between 10

dSTEC Analysis
The results of the dSTEC analysis are shown in Figure 4 to reveal the performance of the products according to the varying VTEC activity in different geographical regions. The geographical distribution of the selected GNSS receivers is illustrated in the panel (a) of Figure 4. The stations used for validation were selected amongst those which exist in both the 2015 and the 2017 time periods. Moreover, we searched in all VTEC products of the analysis centers for a set of geographically well-distributed common stations. Next, we selected GNSS stations either used by all IAACs or are in a very close distance to a used one. The mean RMS values of the dSTEC variations covering the test periods are separately computed for each of the sites and the analysis centers and depicted in the panels (b) and (c) of Figure 4 corresponding to the data sets of 2015 and 2017, respectively. The numbers within the brackets given in the legends show the average values of the RMS deviations computed from all GPS receivers for each analysis center; for further information on the computations, we refer to [5].
The average RMS values of "othg" for 2015 and 2017 are 1.62 and 0.85 for the dSTEC analysis, respectively, whereas the RMS values of the IAACs are varying between 1.64 and 2.33 TECU for 2015 and 0.74 and 1.17 TECU for 2017. The RMS values indicate that the quality of the VTEC product "othg" is consistent with those from the other analysis centers. In particular, the average RMS values of "othg" are in very close agreement with the values of "uqrg" product which has a high temporal resolution.

Altimetry Comparisons
The VTEC measurements acquired from the Jason-2 and Jason-3 missions equipped with dual frequency altimeters were used to assess the quality of the "othg" solution in comparison with the products of the other analysis centers. Prior to the performing analysis, the VTEC data from the altimeter missions were pre-processed to remove outliers and to reduce noise by applying filtering and smoothing steps. Furthermore, VTEC maps of the analysis centers were interpolated to obtain VTEC values at points of altimeter observations as explained by [5]. In the comparisons, data from Jason-2 were used to evaluate the data set of 2015 whereas VTEC data from the Jason-3 mission was considered for the year 2017.
The RMS deviations with respect to altimeter VTEC values are shown in the panels (a) and (b) of Figure 6

Evaluation of the Proposed Approach
This section is devoted to scrutinizing individual contributions of each developed method, including the VCE-based adaptation and the process noise definition, to the performance of the overall KF algorithm. The data set of 2015, which comes with challenging ionosphere conditions, is considered in the evaluations. Several test cases with different filter settings, summarized in Table 1, were designed and the KF was carried out for each case to reveal the model performance. The test case with the label "othg" in Table 1 refers to the final version of the presented adaptive approach.  Firstly, one drawback of the standard KF is the requirement for the prior knowledge of the process noise and the measurement error covariance matrices. In common practice, these matrices are defined by conducting multiple test runs of the filter with different candidate values of the variances defining the covariance matrices. To simulate this situation, the nominal values σ 2 y GPS,NOM , σ 2 y GLO,NOM for the corresponding variance components in Equation (15) and the constant nominal value σ 2 d NOM for the process noise variance of the B-spline coefficients in Equation (19) were manually selected. As it is expected, the results of the altimetry and the dSTEC comparisons for the products derived by the proposed approach with label "othg" and the nominal solution, labeled "TC1" in Table 1, are approximately in line. The RMS values for the altimetry and the dSTEC comparisons of the nominal solution are respectively RMS ALT,TC1 = 5.4 and RMS dSTEC,TC1 = 1.66 TECU. These results show that the "othg" solution exhibits about 2% better performance compared to the nominal solution "TC1". The manual computation procedure, which is what the nominal solution relies on, is time-consuming since it requires executing the filter multiple times to find out the optimal values defining the measurement and process noise covariance matrices. Moreover, it is not guaranteed that the nominal values are valid for any given time since they are computed using a small data set in 2015. However, the presented approach "othg" takes the variances into account in an adaptive way by computing them during the filter run-time. This assures that the values of the variances are always up-to date and sensitive to the changes in errors of GNSS measurements as well as temporal ionospheric variations. For example, the nominal values of the GPS and GLONASS variance components computed using the data set of 2015 are approximately two times larger than the ones computed using the data set of 2017. The data set of 2015 includes considerably larger VTEC variations in magnitude due to high VTEC activity which can significantly degrade the quality of GNSS measurements. The adaptive filter approach successfully handles this issue by assigning larger values to the error covariance matrix of GNSS measurements for the data set of 2015 compared to those computed for the data set of 2017.
Adaptive filters generally allow flexible updates of the underlying model parameters. For instance, in case that the computation of the weight matrices in Equation (15) defined by (16) is replaced by another model, the corresponding value of the nominal variances have to be re-computed in a standard KF. However, the adaptive filtering approach does not require such a time consuming effort. The effect of an inappropriate variance definition to the filter performance is simulated in the test scenario labeled "TC2" in Table 1. The weight matrices P y GPS ,k and P y GLO,k in Equation (15) were set to identity matrices but the nominal values of the corresponding variance components are kept as they were computed in "TC1". The RMS values for the altimetry and dSTEC comparisons for the test case "TC2" are RMS ALT,TC2 = 5.9 TECU and RMS dSTEC,TC2 = 1.94 TECU. The values respectively exhibit about 11% and 17% performance degradation compared to those of the "othg".
To compute ionosphere target parameters reliably, data sets from different sources (e.g., GPS and GLONASS in this study) are required to be handled carefully by means of an appropriate selection of relative weights. The presented adaptive approach deals with this issue by computing different variance components for the observations from each of the constellations. The test case labeled "TC3" is designed to show a performance degradation in case of improper assignment of weights to the observation groups. A unique value, σ 2 y GLO,NOM obtained from the test case "TC1", was assigned to the variance components of GPS and GLONASS observations to conduct the simulation labeled "TC3" in Table 1. The RMS values for the altimetry and dSTEC comparisons for the test case are RMS ALT,TC3 = 5.5 TECU and RMS dSTEC,TC3 = 1.75 TECU. The RMS values show about 4% and 8% performance degradation compared to the corresponding RMS values of the "othg" solution.
The final evaluations are dedicated to assess the effectiveness of the process noise definition given by the Equations (19) and (20) for the B-spline coefficients. In this context, the test case labeled "TC4" in Table 1  The RMS values exhibit about 4% and 2% performance degradation compared to those of the "othg" solution. The RMS values of the dSTEC analysis at the test stations depicted in Figure 7a clearly show that the "othg" solution relying on the signal and data-adaptive process noise definition outperforms the nominal solution over the continental area. A closer look at the RMS values of the dSTEC analysis reveals that the improvements are considerably larger for the receivers at the sites NIUM, KOKB, BOGT, YKRO, COCO, GUAM which are located at low latitudes closer to the geomagnetic equator. These regions generally exhibit very high VTEC activity during the day due to the Equatorial Ionization Anomaly (EIA). The daily RMS values of the altimetry analysis illustrated in Figure 7b indicate that the "othg" solution performs slightly better than the nominal solution for the entire test period of 2015 over oceanic areas. Notable differences between the RMS values of the two solutions can be observed for the days ranging from DOY-72 to DOY-76, which are characterized by very high VTEC activity. During these days, especially the Saint Patrick storm day (DOY-76), the "othg" solution significantly outperforms the nominal solution.
Furthermore, in order to test the sub-parameters of the empirical model defining the process noise variances, a test run of the filter, which is labeled "TC5" in Table 1, was conducted by only considering the C 0,k coefficient in Equation (20). The coefficients C 1,k and C 2,k are set to 1, which means omitting the assumption that the temporal uncertainties of the B-spline coefficients are correlated to their relative variations. The RMS values for the altimetry and dSTEC comparisons are RMS ALT,TC5 = 5.6 TECU and RMS dSTEC,TC5 = 1.76 TECU, where the values show about 6% and 9% performance degradation, respectively, compared to the RMS values of the "othg" solution. Considering that the process noise variance behaves like a weighting factor, the results indicate that the performance of the estimator is significantly sensitive to the relative variations of the B-spline coefficients, which are effectively taken into account by the "othg" solution.

Conclusions
Adaptive approaches allowing for tuning of model uncertainties in run-time are introduced into the Kalman Filter for global VTEC modeling relying on the B-spline representation. The representation of the global VTEC signal using a series expansion of B-spline functions allows proper handling of regional signal variations as well as the heterogeneous data distribution. Incorporating B-splines into a recursive filter results in a very powerful data assimilation framework for ultra-rapid and (near)real-time VTEC modeling. In this regard, the study by [5], which includes the basic recursive algorithms for VTEC product generation using B-splines, was comprehensively extended for adaptive estimation of the unknown ionospheric target parameters.
The adaptive algorithm, incorporated into the KF implementation, makes use of the VCE estimation approach to tune the quality of GPS and GLONASS observations by estimating proper scaling factors (variance components) for each type of observation group during the filter run-time. The adaptive KF approach leads to the construction of the measurement covariance matrix automatically, whereas the standard KF requires manually conducting multiple tests for the identification of the covariance matrix. The equality constraint equations, which are introduced to preserve the spherical geometry of the global VTEC distribution and DCB related restrictions, were inserted into the filter by considering the constraints as perfect measurements.
We also showed that there are very high morphological similarities between the global VTEC distribution and the corresponding B-spline coefficients, which allows deducing physical interpretation directly from the magnitude and the distribution of B-spline coefficients. By exploiting this advantage of the B-spline representation, an empirical model to define the process noise of the filter was developed. The process noise represents the uncertainty of the KF time update (prediction) step; therefore, it plays an essential role in the performance of the filter. In this way, value of the process noise for each coefficient is tuned in an adaptive manner during the course of time according to the structure of the B-spline coefficients.
Furthermore, the hourly block of raw measurements acquired from GPS and GLONASS receivers are provided by IGS's data centers with a latency of 1 h. The ionosphere measurements are obtained from the combination of raw GNSS pseudorange and carrier-phase measurements by applying the method of geometry-free linear combinations.
The performance of the implemented algorithms on VTEC modeling is verified by carrying out the dSTEC analysis and altimetry VTEC comparisons, which shows the generated VTEC products are consistent and in good agreement with the products of the other analysis centers.
The proposed adaptive KF is applied to groups of measurements from the GPS and GLONASS constellations. The approach as well as the data pre-processing module will be extended in the near future to handle the additional constellations, including, e.g., GALILEO and BeiDou. VTEC maps with label "othg", derived using the presented approach, are referred to as ultra-rapid products due to latency in the acquisition of GNSS hourly data from the IGS data centers. However, the adaptive filtering algorithm can be directly used in a real-time modeling study. Therefore, the modeling approach including the data pre-processing will be extended to compute real-time global VTEC as a further step. Moreover, the VCE approach was applied to the measurement covariance matrix. Further studies will be conducted to extend the algorithm for self-tuning of the predicted covariance matrix of the unknown state vector in the sense of VCE.
Author Contributions: The design and concept of the paper were proposed by E.E. The initial full draft of the paper, including preparation of the figures, selection of data set and evaluation of the proposed model, was prepared by E.E. with assistance from M.S and A.G. The software used in the analysis partly developed by E.E. and A.G. The authors M.S., F.S., and B.G. have the responsibility for the project OPTIMAP for which significant contributions entered into this article. All co-authors reviewed and commented on the original manuscript. All authors have read and agreed to the published version of the manuscript.