High-resolution vertical total electron content maps based on multi-scale B-spline representations

For more than 2 decades the IGS (International GNSS Service) ionosphere associated analysis centers (IAACs) have provided global maps of the vertical total electron content (VTEC). In general, the representation of a 2-D or 3-D function can be performed by means of a series expansion or by using a discretization technique. While in the latter case, pixels or voxels are usually chosen for a spherical function such as VTEC, for a series expansion spherical harmonics (SH) are primarily used as basis functions. The selection of the best suited approach for ionosphere modeling means a trade-off between the distribution of available data and their possibility of representing ionospheric variations with high resolution and high accuracy.

Most of the IAACs generate global ionosphere maps (GIMs) based on SH expansions up to the spectral degree n=15 and provide them with a spatial resolution of $\mathrm{2.5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ with respect to the latitudinal and longitudinal directions, respectively, and a temporal sampling interval of 2 h. In recent years, it has frequently been claimed that the spatial resolution of the VTEC GIMs has to be increased to a spatial resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ and to a temporal sampling interval of about 15 min. Enhancing the grid resolution means an interpolation of VTEC values for intermediate points but with no further information about variations in the signal. n=15 in the SH case, for instance, corresponds to a spatial sampling of $\mathrm{12}{}^{\circ }\times \mathrm{12}{}^{\circ }$. Consequently, increasing the grid resolution concurrently requires an extension of the spectral content, i.e., to choose a higher SH degree value than 15.

Unlike most of the IAACs, the VTEC modeling approach at Deutsches Geodätisches Forschungsinstitut der Technischen Universität München (DGFI-TUM) is based on localizing basis functions, namely tensor products of polynomial and trigonometric B-splines. In this way, not only can data gaps be handled appropriately and sparse normal equation systems be established for the parameter estimation procedure, a multi-scale representation (MSR) can also be set up to determine GIMs of different spectral content directly, by applying the so-called pyramid algorithm, and to perform highly effective data compression techniques. The estimation of the MSR model parameters is finally performed by a Kalman filter driven by near real-time (NRT) GNSS data.

Within this paper, we realize the MSR and create multi-scale products based on B-spline scaling, wavelet coefficients and VTEC grid values. We compare these products with different final and rapid products from the IAACs, e.g., the SH model from CODE (Berne) and the voxel solution from UPC (Barcelona). In contrast to the abovementioned products, DGFI-TUM's products are based solely on NRT GNSS observations and ultra-rapid orbits. Nevertheless, we can conclude that the DGFI-TUM's high-resolution product (“othg”) outperforms all products used within the selected time span of investigation, namely September 2017.

1 Introduction

The properties of the atmosphere can be described by means of different variables, e.g., the temperature or the charge state. In the case of temperature, we distinguish between the troposphere (up to a height of about 15 km), the stratosphere (about 15 to 50 km), the mesosphere (about 50 to 90 km), the thermosphere (about 90 to 800 km) and the exosphere (above 800 km) using increasing height above the Earth's surface. In the case of the charge state, the atmosphere is split into the neutral atmosphere (up to a height of 80 km), the ionosphere (about 80 to 1000 km) and the plasmasphere (above 1000 km) (see, e.g., Limberger, 2015).

The ionosphere is mostly driven by the Sun; extreme UV (EUV), X-ray and solar particle radiation cause ionization processes. In geodesy, the main ionospheric impact is the influence of free electrons on radio wave propagation. This effect mainly depends on the signal frequency, i.e., the ionosphere is a dispersive medium (Schaer, 1999). Signals with frequencies lower than 30 MHz will be blocked and reflected by the ionosphere, whereas signals with shorter wavelengths penetrate the ionosphere but are affected with respect to speed and direction. The ionospheric influence on radio waves is twofold: the signal travel times are changed (delay) and the signal paths are modified (bending). While the latter effect can be neglected for most applications, the ionospheric delay

$$\begin{array}{}\text{(1)}& d_{\mathrm{ion}}=\pm {\displaystyle \frac{\mathrm{40.3}}{f^{\mathrm{2}}}}\underset{S}{\overset{R}{\int }}{N}_{\mathrm{e}}\mathrm{d}s\end{array}$$

depends directly on the electron density N_e along the signal path s between satellite S and receiver R and inversely on the carrier frequency f. Equation (1), which can be derived from dual-frequency measurements, is only an approximation as higher-order effects are neglected. These terms depend on the magnetic field, signal frequency, signal elevation, and ionospheric conditions and reach about 0.2 cm in zenith for GPS signals (Bassiri and Hajj, 1993). The sign on the right-hand side changes depending on whether it is applied for a carrier phase observation (“−”) or for a pseudorange measurement (“+”) (see, e.g., Langley, 1998).

Observations of space geodetic techniques, such as the global navigation satellite systems (GNSS) and the Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) tracking system as well as satellite altimetry and ionospheric radio occultation (IRO) are based on electromagnetic signal propagation; thus, they are disturbed by the ionosphere. Most of the techniques are not directly sensitive to the electron density, but to the integrated effect along the ray path. In Eq. (1) the integral

$$\begin{array}{}\text{(2)}& \mathrm{STEC}({\mathit{x}}^S,{\mathit{x}}_R,t)=\underset{S}{\overset{R}{\int }}{N}_{\mathrm{e}}(\mathit{x},t)\mathrm{d}s\end{array}$$

is called the slant total electron content (STEC). In Eq. (2), in addition to the time t we introduce the position vectors x^S x_R and

$$\begin{array}{}\text{(3)}& \mathit{x}=r[\cos \mathit{\phi }\cos \mathit{\lambda },\cos \mathit{\phi }\sin \mathit{\lambda },\sin \mathit{\phi }{]}^T\end{array}$$

of the satellite S, the receiver R and an arbitrary point P moving along the signal path s; the coordinate triple $(\mathit{\phi },\mathit{\lambda },r)$ comprises latitude φ, longitude λ and radial distance r within a geocentric coordinate system Σ_E.

The vertical total electron content (VTEC),

$$\begin{array}{}\text{(4)}& \mathrm{VTEC}(\mathit{\phi },\mathit{\lambda },t)=\underset{h_{\mathrm{l}}}{\overset{h_{\mathrm{u}}}{\int }}{N}_{\mathrm{e}}(\mathit{\phi },\mathit{\lambda },h,t)\mathrm{d}h,\end{array}$$

is defined as the integration of the electron density in the vertical direction, i.e., along the height h above the Earth's surface, defined as $h=r-R_{\mathrm{E}}$; R_E refers to the radius of a spherical Earth. Furthermore, in Eq. (4) h_l and h_u are the respective heights of the lower and the upper boundary of the ionosphere (see, e.g., Dettmering et al., 2011, 2014; Limberger, 2015).

Equations (2) and (4) require a 3-D integration of the electron density. Often a simplification is preferred which is based on the so-called single-layer model (SLM). It assumes that all free electrons are concentrated in an infinitesimally thin shell, i.e., the sphere Ω_H with radius $R_{\mathrm{H}}=R_{\mathrm{E}}+H$ (Schaer, 1999), where H is the single-layer height. As a consequence of this assumption and according to

$$\begin{array}{}\text{(5)}& \mathrm{STEC}({\mathit{x}}^S,{\mathit{x}}_R,t)=M\left(z\right)\cdot \mathrm{VTEC}({\mathit{x}}_{\mathrm{IPP}},t),\end{array}$$

VTEC can be transformed into STEC by introducing a mapping function M(z) depending on the zenith angle z of the ray path between satellite S and receiver R. In Eq. (5) the position vector x_IPP, i.e., the spherical coordinates (φ_IPP, λ_IPP, R_H) define the ionospheric pierce point (IPP), which refers to the geometrical piercing point of the straight line between S and R with the sphere Ω_H of the SLM. This point denotes the reference point of an observation including the STEC, such as a GNSS measurement (see, e.g., Erdogan et al., 2017). Figure 1 shows, for instance, the global distribution of the IPPs from GNSS observations on 6 September 2017 between 12:55 and 13:05 UT. However, it must be pointed out that the introduction of a simple isotropic mapping function M(z), depending solely on the zenith angle z, can only generate an approximation of STEC. Recently, more sophisticated approaches, e.g., the Barcelona ionospheric mapping function (BIMF), have been developed to improve the projection of VTEC on STEC (see, e.g., Lyu et al., 2018).

Combining Eqs. (1), (2) and (5) yields the relation

$$\begin{array}{}\text{(6)}& d_{\mathrm{ion}}({\mathit{x}}^S,{\mathit{x}}_R,t)=-{\displaystyle \frac{\mathrm{40.3}}{f^{\mathrm{2}}}}\cdot M\left(z\right)\cdot \mathrm{VTEC}({\mathit{x}}_{\mathrm{IPP}},t)\end{array}$$

between VTEC and the ionospheric delay d_ion in the case of a phase observation. Equation (6) can be interpreted and applied in two ways:

• if VTEC is given from an ionospheric model, the delay d_ion can be computed and used as a correction for GNSS observations;

• or if the delay d_ion can be derived from double-frequency GNSS measurements, it can be used as an observation to develop or improve VTEC models.

Applications, such as satellite navigation and positioning require high-precision and high-resolution VTEC models. For this purpose the correction d_ion could, according to Eq. (6), be derived from VTEC maps, usually the so-called global ionosphere maps (GIM). The most prominent GIM is provided by the International GNSS Service (IGS) (Feltens and Schaer, 1998; Hernández-Pajares et al., 2011) as a weighted combination product of VTEC maps from various IGS ionosphere associated analysis centers (IAACs), namely (1) the Jet Propulsion Laboratory (JPL), (2) the Center for Orbit Determination in Europe (CODE), (3) the European Space Operations Center of the European Space Agency (ESOC), (4) the Universitat Politècnica de Catalunya (UPC), (5) the Canadian Geodetic Survey of Natural Resources Canada (NRCan), (6) the Wuhan University (WHU) and (7) the Chinese Academy of Sciences (CAS). Recently, Roma-Dollase et al. (2017) published a review paper on these seven GIMs concerning their mapping techniques and their consistency during one solar cycle.

There are several modeling strategies for generating GIMs; the most prominent approach is based on spherical harmonics (SH) and was introduced by Schaer (1999). Further approaches include the tomographic approach based on voxels (Hernández-Pajares et al., 1999) and other approaches based on B-spline scaling functions and wavelets (Schmidt, 2007, 2015b; Schmidt et al., 2011), multivariate adaptive regression splines (MARS) and adaptive regression B-splines (BMARS) (Durmaz et al., 2010; Durmaz and Karslioglu, 2015), and polynomials (Komjathy and Langley, 1996).

Generally, we distinguish between GIMs provided as “final”, “rapid”, “near real-time” (NRT) or “real-time” (RT) products. This classification is based on the latency of the underlying input data. For final products, for instance, only post-processed observations and orbits are used, whereas NRT products are based on rapid orbits and observations with a latency of some minutes up to a few hours. GIMs are typically provided with a temporal resolution of 1 or 2 h and with a spatial resolution of $\mathrm{2.5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ with respect to geographical latitude and longitude (Hernández-Pajares et al., 2017).

VTEC variations basically follow annual, seasonal, diurnal and semidiurnal periods. Earthquakes or incidental natural hazards can also cause small but visible signatures (Liu et al., 2004; Zhu et al., 2013). However, during space weather events, such as solar flares or coronal mass ejections (CME), the number of free electrons may drastically increase. In the latter case, solar plasma consisting of electrons, ions and photons may enter the Earth's atmosphere and cause short period variations within the electron density distribution (see, e.g., Monte-Moreno and Hernández-Pajares, 2014; Wang et al., 2016; Tsurutani et al., 2006, 2009). As a consequence, the modeling of the disturbed ionosphere requires both a high temporal and a high spatial resolution. In 2012, during the IGS 2012 workshop in Olsztyn, Poland, it was recommended that high-resolution IGS combined GIMs be provided. The UPC and JPL IAACs agreed on disseminating GIMs with a temporal resolution of 15 min and a spatial resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ with respect to latitude and longitude, respectively (Dach and Jean, 2013).

As already confirmed by Roma-Dollase et al. (2017), an increase in temporal resolution allows for an improvement in the overall accuracy of the GIMs. The authors compared the final products with a temporal resolution of 2 h with rapid products with a temporal resolution of 15 min using the dSTEC analysis, which is the most reliable method of assessing the accuracy of VTEC products (Hernández-Pajares et al., 2017). Following the results of their investigations, it can be stated that an increase in the temporal resolution yields better results in the dSTEC analysis.

To the knowledge of the authors, the spatial resolution of GIMs has not been investigated in detail to date. Most of the GIMs are based on series expansions in terms of SHs with a maximum degree of n_max=15. This value fits to a block size of about $\mathrm{12}{}^{\circ }\times \mathrm{12}{}^{\circ }$ on the sphere Ω_H. In contrast, a grid spacing of $\mathrm{2.5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ corresponds to a maximum SH degree of around n=36; a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid spacing, i.e., a spatial resolution of around 110 km along the Equator fits to a SH expansion up to degree n=180. As a matter of fact, a reliable computation of the corresponding SH series coefficients requires a global input data coverage of the same spatial sampling. As the IAAC VTEC maps are based solely on GNSS observations with a rather inhomogeneous distribution (cf. Fig. 1 showing the IPPs of NRT observations with dense clusters over continents and large data gaps over oceans), finer ionospheric structures can only be monitored and modeled where high-resolution input data are available.

Figure 1 Global distribution of the IPPs from GPS (red dots) and GLONASS (blue stars) measurements for 6 September 2017 collected within a 10 min interval between 12:55 and 13:05 UT. The regional maps at the top are “zoom-ins” of Europe and Indonesia.

By increasing the temporal resolution of the GIMs, the number of observations supporting the individual maps decreases. The two “zoom-in” maps at the top of Fig. 1 show the strong incongruity between the data distribution and the signal structure (cf. Fig. 9a, c). In areas with high-resolution data, such as Europe, the US or Australia, the VTEC signal is usually rather smooth. In areas with highly variable spatial and temporal signal structures such as in the equatorial belt, a much smaller number of observations is generally given. As a consequence, for global modeling we have to deal with a trade-off between signal structure and data resolution.

It is a well-known fact that SHs as global basis functions are not suitable for representing unevenly, globally distributed data. Consequently, in such cases, a series expansion in terms of localizing basis functions is more appropriate. In the following, we apply tensor products of polynomial and trigonometric B-splines as localizing 2-D basis functions. Besides the localizing features, B-splines additionally generate a multi-scale representation (MSR), also known as multi-resolution representation (MRR). The basic feature of a MSR is to split a target function into a smoothed, i.e., low-pass-filtered version, and a number of detail signals, i.e., band-pass-filtered versions via successive low-pass filtering (Mertins, 1999). Hence, a spatial MSR of VTEC adapts the model resolution to the data distribution and, thus, fulfills IGS' requirement of high-resolution VTEC modeling.

In this study, we compare global VTEC maps based on series expansions in terms of both globally defined SHs and localizing B-spline functions, including the MSR with respect to the spectral content. For this purpose, we use the SH degree as the common measure for the spectral content of a spherical signal. In detail, we study the interrelations between the SH degree, the spatial sampling intervals of the input data and the resolution levels of B-spline expansions. In addition, we discuss the influence of different temporal resolutions of the GIMs. For the estimation of the unknown series coefficients of the B-spline expansion, we use a Kalman filter (KF) procedure as explained by Erdogan et al. (2017). In order to assess the quality of the approach, we perform a dSTEC analysis (Hernández-Pajares et al., 2017).

The paper is outlined as follows: in Sect. 2 a description of VTEC modeling procedures based on both, SH and B-spline expansions are presented. In Sect. 2.3 we study the spectral resolution of global VTEC maps. Section 3 comprises a detailed description of the MSR and the estimation procedure. In Sect. 4 case studies are set up to verify the results of the previous sections numerically. Furthermore, this section provides a final assessment by means of a dSTEC validation. The final section provides conclusions and an outlook for future work.

2 VTEC modeling approaches

The 3-D signal VTEC$(\mathit{\phi },\mathit{\lambda },t)=f(\mathit{x},t)$, introduced in Eqs. (4) and (5), can be modeled as series expansion

$$\begin{array}{}\text{(7)}& f(\mathit{x},t)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{c}_k\left(t\right){\mathit{\varphi }}_k\left(\mathit{x}\right)\end{array}$$

in terms of given space-dependent basis functions ϕ_k(x) and unknown time-dependent series coefficients c_k(t). ¹ Assuming that at discrete times $t_s=t_{\mathrm{0}}+s\cdot \mathrm{\Delta }t$ with s∈ℕ₀ and sampling interval Δt the total number of I_s observations y(${\mathit{x}}_{i_s}$, t_s) of VTEC at IPP position $P_{i_s}\in {\mathrm{\Omega }}_{\mathrm{H}}$ with i_s=1, 2, …, I_s are given. Considering the measurement errors e(${\mathit{x}}_{i_s}$, t_s), the observation equation follows from Eq. (7) and reads

$$\begin{array}{}\text{(8)}& \begin{array}{rl}y({\mathit{x}}_{i_s},t_s)+e({\mathit{x}}_{i_s},t_s)& =f_N({\mathit{x}}_{i_s},t_s)=\\ & =\sum _{k=\mathrm{0}}^N{c}_k\left(t_s\right){\mathit{\varphi }}_k\left({\mathit{x}}_{i_s}\right).\end{array}\end{array}$$

Note, that we neglect the truncation error in the following

$$\begin{array}{}\text{(9)}& r_N({\mathit{x}}_{i_s},t_s)=\sum _{k=N+\mathrm{1}}^{\mathrm{\infty }}{c}_k\left(t_s\right){\mathit{\varphi }}_k\left({\mathit{x}}_{i_s}\right)\end{array}$$

and omit other unknown parameters such as the satellite and receiver differential code biases (DCB) for GNSS geometry-free observations on the right-hand side of Eq. (8) (see, e.g., Erdogan et al., 2017).

In the following (Sect. 2.1 and 2.2), the SH expansion – as likely the most frequently used approach in ionosphere modeling – and the 2-D B-spline tensor product approach are described.

2.1 Spherical harmonic expansion

In the SH approach, the observation equation, Eq. (8), can be rewritten as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}y({\mathit{x}}_{i_s},t_s)& +e({\mathit{x}}_{i_s},t_s)=f_{n_{\max }}({\mathit{x}}_{i_s},t_s)=\\ & =\sum _{n=\mathrm{0}}^{n_{\max }}\sum _{m=-n}^n{c}_{n,m}\left(t_s\right)Y_{n,m}\left({\mathit{x}}_{i_s}\right),\end{array}\end{array}$$

where the functions Y_n,m(x), i.e., the SHs of degree n=0, …, n_max and order $m=-n$, …, n, are defined as

$$\begin{array}{}\text{(11)}& Y_{n,m}\left(\mathit{x}\right)=P_{n,\left|m\right|}\left(\sin \mathit{\phi }\right)\cdot \left\{\begin{array}{l}\cos m\mathit{\lambda }\text{if}m\ge \mathrm{0}\\ \sin \left|m\right|\mathit{\lambda }\text{if}m<\mathrm{0}\end{array}\right.\end{array}$$

with $P_{n,\left|m\right|}$ being the normalized associated Legendre functions of degree n and order m. The total (n_max+1)² quantities c_n,m(t) in Eq. (10) are the time-dependent SH coefficients. According to the sampling theorem on the sphere, the maximum degree n_max is related to the sampling intervals Δφ and Δλ of the input data with respect to latitude φ and longitude λ, namely

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }\mathit{\phi }<{\displaystyle \frac{\mathrm{180}{}^{\circ }}{n_{\max }}}\text{and}\mathrm{\Delta }\mathit{\lambda }<{\displaystyle \frac{\mathrm{180}{}^{\circ }}{n_{\max }}}.\end{array}$$

As can be seen from Eq. (11), SHs are basis functions of global support. This implies that each single SH function is different from zero almost everywhere on the sphere Ω_H. Consequently, each coefficient c_n,m has to be recomputed if only one additional observation is considered in the set of observation Eq. (10).

As the VTEC observations $y({\mathit{x}}_{i_s},t_s)$ at the IPP positions will usually not be given on a spatial grid with constant mesh size, the sampling intervals Δφ and Δλ in the formulae of Eq. (12) have to be interpreted as global average values.

2.2 B-spline expansion

At DGFI-TUM we rely on B-splines as basis functions for ionosphere modeling, as they are (1) characterized by their localizing feature and (2) they can be used to generate a MSR. For VTEC modeling we rewrite Eq. (8) as

$$\begin{array}{}\text{(13)}& \begin{array}{rl}y({\mathit{x}}_{i_s},t_s)& +e({\mathit{x}}_{i_s},t_s)=f_{J_{\mathrm{1}},J_{\mathrm{2}}}({\mathit{x}}_{i_s},t_s)=\\ & =\sum _{k_{\mathrm{1}}=\mathrm{0}}^{K_{J_{\mathrm{1}}}-\mathrm{1}}\sum _{k_{\mathrm{2}}=\mathrm{0}}^{K_{J_{\mathrm{2}}}-\mathrm{1}}{d}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right){\mathit{\varphi }}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}({\mathit{\phi }}_{i_s},{\mathit{\lambda }}_{i_s})\end{array}\end{array}$$

with initially unknown time-dependent scaling coefficients $d_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right)$ and the 2-D scaling functions ${\mathit{\varphi }}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}({\mathit{\phi }}_{i_s},{\mathit{\lambda }}_{i_s})$ of levels J₁ and J₂ with respect to φ and λ. The latter are defined as tensor products

$$\begin{array}{}\text{(14)}& {\mathit{\varphi }}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })={\mathit{\varphi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}\left(\mathit{\phi }\right){\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)\end{array}$$

of the 1-D scaling functions ${\mathit{\varphi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}\left(\mathit{\phi }\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)$ depending on latitude φ and longitude λ, respectively. As the B-spline approach is not as well known as the SH approach, it will be described in more detail in the following; we further refer to Dierckx (1984); Stollnitz et al. (1995a, b); Lyche and Schumaker (2001); Jekeli (2005); Schmidt (2015b) and citations therein.

To decompose VTEC into its spectral components via the MSR in Sect. 3, Eqs. (13) and (14) need to be rewritten in vector and matrix notation. For this purpose we introduce the $K_{J_{\mathrm{1}}}\times \mathrm{1}$ vector

$$\begin{array}{}\text{(15)}& {\mathit{\varphi }}_{J_{\mathrm{1}}}\left(\mathit{\phi }\right)=\left[{\mathit{\varphi }}_{\mathrm{0}}^{J_{\mathrm{1}}}\right(\mathit{\phi }),{\mathit{\varphi }}_{\mathrm{1}}^{J_{\mathrm{1}}}(\mathit{\phi }),\mathrm{\dots },{\mathit{\varphi }}_{K_{J_{\mathrm{1}}}-\mathrm{1}}^{J_{\mathrm{1}}}(\mathit{\phi }){]}^T,\end{array}$$

the $K_{J_{\mathrm{2}}}\times \mathrm{1}$ vector

$$\begin{array}{}\text{(16)}& {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)=\left[{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{0}}^{J_{\mathrm{2}}}\right(\mathit{\lambda }),{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{1}}^{J_{\mathrm{2}}}(\mathit{\lambda }),\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{K_{J_{\mathrm{2}}}-\mathrm{1}}^{J_{\mathrm{2}}}(\mathit{\lambda }){]}^T\end{array}$$

and the $K_{J_{\mathrm{1}}}\times K_{J_{\mathrm{2}}}$ coefficient matrix

$$\begin{array}{}\text{(17)}& {\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}=\left[\begin{array}{llll}{d}_{\mathrm{0},\mathrm{0}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& d_{\mathrm{0},\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& \mathrm{\dots }& d_{\mathrm{0},K_{J_{\mathrm{2}}}-\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\\ d_{\mathrm{1},\mathrm{0}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& d_{\mathrm{1},\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& \mathrm{\dots }& d_{\mathrm{1},K_{J_{\mathrm{2}}}-\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{⋮}\\ d_{K_{J_{\mathrm{1}}}-\mathrm{1},\mathrm{0}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& d_{K_{J_{\mathrm{1}}}-\mathrm{1},\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}& \mathrm{\dots }& d_{K_{J_{\mathrm{1}}}-\mathrm{1},K_{J_{\mathrm{2}}}-\mathrm{1}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\end{array}\right].\end{array}$$

Considering the computation rules for the Kronecker product “⊗” (cf. Koch, 1999), Eq. (13) can be written as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}f(\mathit{\phi },\mathit{\lambda },t)& =\left({\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\right(\mathit{\lambda })\otimes {\mathit{\varphi }}_{J_{\mathrm{1}}}(\mathit{\phi }\left){)}^T\text{vec}{\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}\right(t)\\ & ={\mathit{\varphi }}_{J_{\mathrm{1}}}^T\left(\mathit{\phi }\right){\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t\right){\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right),\end{array}\end{array}$$

where “vec” refers to the vec operator.

2.2.1 Polynomial B-splines

In the following, we apply polynomial quadratic B-splines

$$\begin{array}{}\text{(19)}& {\mathit{\varphi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}\left(\mathit{\phi }\right):=N_{J_{\mathrm{1}},k_{\mathrm{1}}}^{\mathrm{2}}\left(\mathit{\phi }\right)\end{array}$$

of resolution level J₁∈ℕ₀ and shift k₁=0, 1, …, $K_{J_{\mathrm{1}}}-\mathrm{1}$ to represent the latitude-dependent variations of VTEC. To be more specific, a total of $K_{J_{\mathrm{1}}}={\mathrm{2}}^{J_{\mathrm{1}}}+\mathrm{2}$ B-splines are located along a meridian depending on the latitude $\mathit{\phi }\in [-\mathrm{90}{}^{\circ },\mathrm{90}{}^{\circ }]$. To construct the $K_{J_{\mathrm{1}}}$ B-spline functions, the sequence

$$\begin{array}{}\text{(20)}& \begin{array}{rl}-\mathrm{90}{}^{\circ }={\mathit{\phi }}_{\mathrm{0}}^{J_{\mathrm{1}}}={\mathit{\phi }}_{\mathrm{1}}^{J_{\mathrm{1}}}& ={\mathit{\phi }}_{\mathrm{2}}^{J_{\mathrm{1}}}<{\mathit{\phi }}_{\mathrm{3}}^{J_{\mathrm{1}}}<\mathrm{\dots }<{\mathit{\phi }}_{K_{J_{\mathrm{1}}}}^{J_{\mathrm{1}}}=\\ & ={\mathit{\phi }}_{K_{J_{\mathrm{1}}}+\mathrm{1}}^{J_{\mathrm{1}}}={\mathit{\phi }}_{K_{J_{\mathrm{1}}}+\mathrm{2}}^{J_{\mathrm{1}}}=\mathrm{90}{}^{\circ }\end{array}\end{array}$$

of knot points ${\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}$ is established; the consideration of multiple knot points at the poles is called “endpoint-interpolating” and ensures the closing of the modeling interval. The constant distance between two consecutive knots ${\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}$ and ${\mathit{\phi }}_{k_{\mathrm{1}}+\mathrm{1}}^{J_{\mathrm{1}}}$ for k₁=2, …, $K_{J_{\mathrm{1}}}-\mathrm{1}$ amounts to $\mathrm{180}{}^{\circ }/{\mathrm{2}}^{J_{\mathrm{1}}}$. Following Schumaker and Traas (1991) and Stollnitz et al. (1995b) the normalized quadratic polynomial B-splines can be calculated via the recursive relation

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{N}_{J_{\mathrm{1}},k_{\mathrm{1}}}^n\left(\mathit{\phi }\right)& ={\displaystyle \frac{\mathit{\phi }-{\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}}{{\mathit{\phi }}_{k_{\mathrm{1}}+n}^{J_{\mathrm{1}}}-{\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}}}{N}_{J_{\mathrm{1}},k_{\mathrm{1}}}^{n-\mathrm{1}}\left(\mathit{\phi }\right)\\ & +{\displaystyle \frac{{\mathit{\phi }}_{k_{\mathrm{1}}+n+\mathrm{1}}^{J_{\mathrm{1}}}-\mathit{\phi }}{{\mathit{\phi }}_{k_{\mathrm{1}}+n+\mathrm{1}}^{J_{\mathrm{1}}}-{\mathit{\phi }}_{k_{\mathrm{1}}+\mathrm{1}}^{J_{\mathrm{1}}}}}{N}_{J_{\mathrm{1}},k_{\mathrm{1}}+\mathrm{1}}^{n-\mathrm{1}}\left(\mathit{\phi }\right),\end{array}\end{array}$$

with $n=\mathrm{1},\mathrm{2}$ from the initial values

$$N_{J_{\mathrm{1}},k_{\mathrm{1}}}^{\mathrm{0}}\left(\mathit{\phi }\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if}{\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}\le \mathit{\phi }<{\mathit{\phi }}_{k_{\mathrm{1}}+\mathrm{1}}^{J_{\mathrm{1}}}\text{and}{\mathit{\phi }}_{k_{\mathrm{1}}}^{J_{\mathrm{1}}}<{\mathit{\phi }}_{k_{\mathrm{1}}+\mathrm{1}}^{J_{\mathrm{1}}}\\ \mathrm{0}& \text{otherwise}.\end{array}\right.$$

Note, in Eq. (21) a factor is set to zero if the denominator is equal to zero.

As can be seen from Fig. 2, B-splines are characterized by their compact support or – in other words – they are different from zero only within a small subinterval of length ${\mathrm{\Delta }}_{J_{\mathrm{1}}}\approx \mathrm{3}\cdot h_{J_{\mathrm{1}}}$, where

$$\begin{array}{}\text{(22)}& h_{J_{\mathrm{1}}}={\displaystyle \frac{\mathrm{180}{}^{\circ }}{{\mathrm{2}}^{J_{\mathrm{1}}}+\mathrm{1}}}\end{array}$$

refers to the approximate distance between two consecutive B-splines along the meridian. As the total number $K_{J_{\mathrm{1}}}$ of B-splines depends on the level J₁, finer structures can be modeled by increasing J₁. The numerical value for the level J₁ depends on the global average value Δφ for the input data sampling interval in the latitudinal direction according to

$$\begin{array}{}\text{(23)}& \mathrm{\Delta }\mathit{\phi }<h_{J_{\mathrm{1}}}\end{array}$$

(Schmidt et al., 2011). Solving Eq. (23), considering the level value J₁ from Eq. (22), the inequality

$$\begin{array}{}\text{(24)}& J_{\mathrm{1}}\le {\log }_{\mathrm{2}}\left({\displaystyle \frac{\mathrm{180}{}^{\circ }}{\mathrm{\Delta }\mathit{\phi }}}-\mathrm{1}\right)\end{array}$$

results.

Figure 2 Polynomial B-splines of level J₁=3 with a total number $K_{J_{\mathrm{1}}}={\mathrm{2}}^{\mathrm{3}}+\mathrm{2}=\mathrm{10}$ of B-splines along the meridian. The blue spline function $N_{\mathrm{3},\mathrm{4}}^{\mathrm{2}}\left(\mathit{\phi }\right)$ corresponds to the shift value k₁=4 and covers a subinterval of length ${\mathrm{\Delta }}_{\mathrm{3}}\approx \mathrm{3}\cdot \mathrm{180}{}^{\circ }/\mathrm{9}=\mathrm{60}{}^{\circ }$. The red spline functions $N_{\mathrm{3},\mathrm{0}}^{\mathrm{2}}\left(\mathit{\phi }\right)$ with shift value k₁=0 and $N_{\mathrm{3},\mathrm{9}}^{\mathrm{2}}\left(\mathit{\phi }\right)$ with shift value k₁=9 close the modeling interval at the poles.

2.2.2 Trigonometric B-splines

For modeling the longitudinal variations of VTEC trigonometric B-splines $T_{J_{\mathrm{2}},k_{\mathrm{2}}}^{\mathrm{3}}\left(\mathit{\lambda }\right)$ of order 3 and, depending on the resolution level, J₂∈ℕ₀ and shift k₂=0, 1, …, $K_{J_{\mathrm{2}}}-\mathrm{1}$ are applied. As can be seen from Fig. 3, the total number $K_{J_{\mathrm{2}}}=\mathrm{3}\cdot {\mathrm{2}}^{J_{\mathrm{2}}}$ of trigonometric B-splines are located along the parallels of the chosen spherical coordinate system within the interval $\mathit{\lambda }\in (\mathrm{0}{}^{\circ },\mathrm{360}{}^{\circ })$. Consequently, the first and the last two B-spline functions within the interval $(\mathrm{0}{}^{\circ },\mathrm{360}{}^{\circ })$ have to be completed by the so-called “wrapping around” effect. This constraint allows trigonometric B-splines to be defined in two different ways:

1. Following Schumaker and Traas (1991), Jekeli (2005) and Limberger (2015) periodic trigonometric B-splines can be calculated by means of a recurrence relation similar to Eq. (21). Thereby, additional constraints have to be introduced to force the periodicity of the series coefficients.

2. The second option was introduced by Lyche and Schumaker (2001) and used by Schmidt et al. (2011) and Schmidt (2015b). It will be described in the following in more detail.

To be more specific, the sequence of nondecreasing knot points

$$\begin{array}{}\text{(25)}& \mathrm{0}{}^{\circ }={\mathit{\lambda }}_{\mathrm{0}}^{J_{\mathrm{2}}}<{\mathit{\lambda }}_{\mathrm{1}}^{J_{\mathrm{2}}}<\mathrm{\dots }<{\mathit{\lambda }}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}<\mathrm{\dots }<{\mathit{\lambda }}_{K_{J_{\mathrm{2}}}-\mathrm{1}}^{J_{\mathrm{2}}}<\mathrm{360}{}^{\circ },\end{array}$$

with additional knots

$$\begin{array}{}\text{(26)}& {\mathit{\lambda }}_{K_{J_{\mathrm{2}}}+i}^{J_{\mathrm{2}}}={\mathit{\lambda }}_i^{J_{\mathrm{2}}}+\mathrm{360}{}^{\circ }\text{for}i=\mathrm{0},\mathrm{1},\mathrm{2}\end{array}$$

for considering the periodicity is introduced. Similar to the polynomial B-splines, the distance between consecutive knots ${\mathit{\lambda }}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}$ and ${\mathit{\lambda }}_{k_{\mathrm{2}}+\mathrm{1}}^{J_{\mathrm{2}}}$ for k₂=0, 1, …, $K_{J_{\mathrm{2}}}+\mathrm{1}$ is given as

$$\begin{array}{}\text{(27)}& h_{J_{\mathrm{2}}}={\displaystyle \frac{\mathrm{360}{}^{\circ }}{K_{J_{\mathrm{2}}}}}={\displaystyle \frac{\mathrm{120}{}^{\circ }}{{\mathrm{2}}^{J_{\mathrm{2}}}}};\end{array}$$

thus, the length of the nonzero subinterval of a trigonometric B-spline function $T_{J_{\mathrm{2}},k_{\mathrm{2}}}^{\mathrm{3}}\left(\mathit{\lambda }\right)$ reads ${\mathrm{\Delta }}_{J_{\mathrm{2}}}=\mathrm{3}\cdot h_{J_{\mathrm{2}}}=\mathrm{360}{}^{\circ }/{\mathrm{2}}^{J_{\mathrm{2}}}$. Following Lyche and Schumaker (2001) we define the functions

$$\begin{array}{}\text{(28)}& M_{J_{\mathrm{2}},k_{\mathrm{2}}}\left(\mathit{\lambda }\right)=T_{J_{\mathrm{2}},k_{\mathrm{2}}}^{\mathrm{3}}\left(\mathit{\lambda }\right)=T_{h_{J_{\mathrm{2}}}}^{\mathrm{3}}(\mathit{\lambda }-{\mathit{\lambda }}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}).\end{array}$$

Setting $h_{J_{\mathrm{2}}}=h$ and $\mathit{\lambda }-{\mathit{\lambda }}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}=\mathrm{\Theta }$ for the sake of simplification, the functions $T_{h_{J_{\mathrm{2}}}}^{\mathrm{3}}(\mathit{\lambda }-{\mathit{\lambda }}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}})=T_h^{\mathrm{3}}\left(\mathrm{\Theta }\right)$ can be calculated via

$$\begin{array}{}\text{(29)}& T_h^{\mathrm{3}}\left(\mathrm{\Theta }\right)=\left\{\begin{array}{l}{\displaystyle \frac{{\sin }^{\mathrm{2}}(\mathrm{\Theta }/\mathrm{2})}{\sin (h/\mathrm{2})\sin \left(h\right)}}\text{for}\mathrm{0}\le \mathrm{\Theta }<h\\ {\displaystyle \frac{\mathrm{1}}{\cos (h/\mathrm{2})}}-{\displaystyle \frac{{\sin }^{\mathrm{2}}\left(\right(\mathrm{\Theta }-h)/\mathrm{2})+{\sin }^{\mathrm{2}}\left(\right(\mathrm{2}h-\mathrm{\Theta })/\mathrm{2})}{\sin (h/\mathrm{2})\sin \left(h\right)}}\\ \text{for}h\le \mathrm{\Theta }<\mathrm{2}h\\ {\displaystyle \frac{{\sin }^{\mathrm{2}}\left(\right(\mathrm{3}h-\mathrm{\Theta })/\mathrm{2})}{\sin (h/\mathrm{2})\sin \left(h\right)}}\text{for}\mathrm{2}h\le \mathrm{\Theta }<\mathrm{3}h\\ \mathrm{0}\text{otherwise}.\end{array}\right.\end{array}$$

Finally, we define the basis functions

$$\begin{array}{}\text{(30)}& {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)=\left\{\begin{array}{l}{M}_{J_{\mathrm{2}},k_{\mathrm{2}}}\left(\mathit{\lambda }\right)\text{for}{k}_{\mathrm{2}}=\mathrm{0},\mathrm{\dots },K_{J_{\mathrm{2}}}-\mathrm{3}\\ M_{J_{\mathrm{2}},k_{\mathrm{2}}}\left(\mathit{\lambda }\right)+M_{J_{\mathrm{2}},k_{\mathrm{2}}}(\mathit{\lambda }-\mathrm{360}{}^{\circ })\\ \text{for}{k}_{\mathrm{2}}=K_{J_{\mathrm{2}}}-\mathrm{2},K_{J_{\mathrm{2}}}-\mathrm{1}\end{array}\right.\end{array}$$

introduced in Eq. (14). Figure 3 shows trigonometric B-splines of level J₂=2. With larger values for level J₂ the splines become more narrow and finer structures can be modeled. The choice of the level value J₂ again depends on the input data sampling interval. Analog to Eq. (23), the inequality

$$\begin{array}{}\text{(31)}& \mathrm{\Delta }\mathit{\lambda }<h_{J_{\mathrm{2}}}\end{array}$$

has to be fulfilled, where Δλ denotes the global average value of the data sampling interval in the longitudinal direction. Finally, taking Eq. (27) into consideration, the inequality

$$\begin{array}{}\text{(32)}& J_{\mathrm{2}}\le {\log }_{\mathrm{2}}\left({\displaystyle \frac{\mathrm{120}{}^{\circ }}{\mathrm{\Delta }\mathit{\lambda }}}\right)\end{array}$$

for the level value J₂ is obtained.

Figure 3 $K_{J_{\mathrm{2}}}=\mathrm{3}\cdot {\mathrm{2}}^{\mathrm{2}}=\mathrm{12}$ trigonometric B-splines ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{k_{\mathrm{2}}}^{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)$, according to Eq. (30) for level J₂=2. The blue spline function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{5}}^{\mathrm{2}}\left(\mathit{\lambda }\right)$ with shift value k₂=5 is different from zero only in the subinterval of length ${\mathrm{\Delta }}_{J_{\mathrm{2}}}=\mathrm{360}{}^{\circ }/\mathrm{4}\approx \mathrm{90}{}^{\circ }$. The red basis function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{11}}^{\mathrm{2}}\left(\mathit{\lambda }\right)$ shows the “wrapping-around” effect.

2.3 Spectral resolution of global VTEC models

In Sect. 2.1 we derived the relations between the maximum degree n_max of a SH expansion and the sampling intervals Δφ and Δλ of the input data. In Sect. 2.2 the corresponding relations between the level values J₁ and J₂ of a B-spline expansion and the data sampling intervals have been deduced. The substitution of the expression $\mathrm{180}{}^{\circ }/n_{\max }$ from the inequalities Eq. (12) into Eqs. (24) and (32) yields a total of six inequalities

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{J}_{\mathrm{1}}& \le {\log }_{\mathrm{2}}({\displaystyle \frac{\mathrm{180}{}^{\circ }}{\mathrm{\Delta }\mathit{\phi }}}-\mathrm{1})\le {\log }_{\mathrm{2}}\left(n_{\max }-\mathrm{1}\right),\\ J_{\mathrm{2}}& \le {\log }_{\mathrm{2}}\left({\displaystyle \frac{\mathrm{120}{}^{\circ }}{\mathrm{\Delta }\mathit{\lambda }}}\right)\le {\log }_{\mathrm{2}}\left({\displaystyle \frac{\mathrm{2}\cdot n_{\max }}{\mathrm{3}}}\right).\end{array}\end{array}$$

Given the numerical values 1 to 6 for the B-spline levels J₁ and J₂ Table 1 presents the corresponding largest numerical values for each, the SH degree n_max and the sampling intervals Δφ and Δλ by evaluating the inequalities using Eq. (33).

Table 1 Numerical values for the B-spline levels J₁ and J₂, the maximum SH degree n_max and the input data sampling intervals Δφ and Δλ by evaluating the inequalities from Eq. (33); the left part of the table presents the numbers along a meridian (upper inequalities in Eq. 33), and the right part represents the corresponding numbers along the Equator and its parallels according to the lower inequalities in Eq. (33).

From the spectral point of view the six inequalities from Eq. (33) comprise the following three scenarios:

1. If the global sampling intervals Δφ and Δλ are known, the mid parts of the inequalities from Eq. (33) are given. The maximum degree n_max is calculable from the right-hand side inequalities and may be inserted into the SH expansion in Eq. (10). The left-hand side inequalities yield the two level values J₁ and J₂ which can be inserted into the B-spline expansion from Eq. (13).

2. With a specified numerical value for n_max the right-hand parts of the inequalities from Eq. (33) are given. The data input sampling intervals Δφ and Δλ can be determined from the mid parts of the inequalities. Next, the two numerical values for the level values J₁ and J₂ can be calculated from the left-hand side inequalities and can be inserted into the B-spline expansion in Eq. (13).

3. If the processing time of VTEC maps has to be considered, the level values J₁ and J₂ are subject to certain restrictions; this is due to the fact that the number of numerical operations increases exponentially with the chosen numerical values for the levels. In this case, from the given left-hand side inequalities, the data sampling intervals Δφ and Δλ can be determined from the mid parts. Finally, the right-hand side inequalities yield numerical values for the maximum SH degree n_max.

As already mentioned in the introduction, most of the GIMs produced by the IAACs are based on series expansions in SHs up to a maximum degree of n_max=15. Following the abovementioned second strategy and Table 1, we obtain the approximations J₁=4 (for n_max=17) and J₂=3 (for n_max=12) for the two B-spline levels J₁ and J₂ for this example. Inserting these numbers into the B-spline expansion Eq. (13) yields the spectrally closest representation to the current IGS solutions. A numerical verification of this choice will be presented in Sect. 4.3.

2.4 VTEC output grids

The VTEC GIMs of the IAACs are usually provided with a spatial resolution of $\mathrm{\Delta }\mathrm{\Phi }=\mathrm{2.5}{}^{\circ }$ in the latitudinal direction and $\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{5}{}^{\circ }$ in the longitudinal direction as well as a temporal sampling of ΔT=2 h. Note, the resolution intervals ΔΦ, ΔΛ and ΔT are usually distinct from the sampling intervals Δφ, Δλ and Δt of the observations introduced in Sect. 2.1.

In order to calculate a VTEC value VTEC$(\mathit{\phi },\mathit{\lambda },t)$ at an arbitrary location $P(\mathit{\phi }={\mathit{\phi }}_i+q\cdot \mathrm{\Delta }\mathrm{\Phi },\mathit{\lambda }={\mathit{\lambda }}_k+p\cdot \mathrm{\Delta }\mathrm{\Lambda })$ with $\mathrm{0}\le q\le \mathrm{1}$ and $\mathrm{0}\le p\le \mathrm{1}$ at an arbitrary time moment t, a simple bilinear spatial interpolation from the VTEC values of the four given corner points P(φ_i,λ_k), $P({\mathit{\phi }}_i,{\mathit{\lambda }}_k+\mathrm{\Delta }\mathrm{\Lambda })$, $P({\mathit{\phi }}_i+\mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k)$ and $P({\mathit{\phi }}_i+\mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k+\mathrm{\Delta }\mathrm{\Lambda })$ is performed according to

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\mathrm{VTEC}({\mathit{\phi }}_i& +q\cdot \mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k+p\cdot \mathrm{\Delta }\mathrm{\Lambda },t)\\ & =(\mathrm{1}-q)\cdot (\mathrm{1}-p)\cdot \mathrm{VTEC}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t)\\ & +q\cdot (\mathrm{1}-p)\cdot \mathrm{VTEC}({\mathit{\phi }}_i+\mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k,t)\\ & +p\cdot (\mathrm{1}-q)\cdot \mathrm{VTEC}({\mathit{\phi }}_i,{\mathit{\lambda }}_k+\mathrm{\Delta }\mathrm{\Lambda },t)\\ & +q\cdot p\cdot \mathrm{VTEC}({\mathit{\phi }}_i+\mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k+\mathrm{\Delta }\mathrm{\Lambda },t)\end{array}\end{array}$$

(see, e.g., Schaer et al., 1998, Fig. 4 in this paper).

Figure 4 Schematic representation of the four-point spatial interpolation to calculate the VTEC value at $P({\mathit{\phi }}_i+q\cdot \mathrm{\Delta }\mathrm{\Phi },{\mathit{\lambda }}_k+p\cdot \mathrm{\Delta }\mathrm{\Lambda })$ from the four corner points of the grid cell of interest.

Note, by applying the interpolation formula (34), the quality of the calculated VTEC value decreases with increasing spatial resolution intervals ΔΦ and ΔΛ and depends on the position within the grid cell. In order to improve the quality of the VTEC computation two methods can be used:

1. the chosen model approach, e.g., the SH or the B-spline expansion can be used directly to calculate VTEC values at any arbitrary point P(φ,λ);

2. the resolution intervals ΔΦ and ΔΛ of the output grid can be set to smaller values, e.g., to 1^∘, as was proposed at the IGS workshop 2012.

For the calculation of a VTEC value VTEC$(\mathit{\phi },\mathit{\lambda },t)$ at an arbitrary time moment $t=t_s+r\cdot \mathrm{\Delta }T$ with $\mathrm{0}\le r\le \mathrm{1}$ at a given spatial location P(φ,λ), an interpolation with respect to time can be applied. Commonly, the linear interpolation

$$\begin{array}{}\text{(35)}& \begin{array}{rl}\mathrm{VTEC}(\mathit{\phi },\mathit{\lambda },t)& =(\mathrm{1}-r)\cdot \mathrm{VTEC}(\mathit{\phi },\mathit{\lambda },t_s)\\ & +r\cdot \mathrm{VTEC}(\mathit{\phi },\mathit{\lambda },t_s+\mathrm{\Delta }T)\end{array}\end{array}$$

between the two consecutive maps at epochs t_s and t_s+ΔT is performed (see, e.g., Schaer et al., 1998).

The previously described interpolation methods allow for the calculation of VTEC values VTEC$(\mathit{\phi },\mathit{\lambda },t)$ at any spatial location P(φ,λ) and at any time t. However, for a more accurate calculation of VTEC an increase in the resolution is necessary for both domains. In the following, it is shown that the usage of a MSR based on the B-spline approach in combination with a KF estimation procedure provides the possibility to create VTEC maps of higher spatial and temporal resolution. Consequently, according to Table 1 the calculated VTEC maps cover a wider spectral band, i.e., the numerical value of n_max becomes larger.

3 Multi-scale representation

The B-spline functions as introduced in Sect. 2.2.1 and 2.2.2 allow for the generation of a MSR. To be more specific, B-spline tensor product wavelet functions will be constructed which are intrinsically connected to the resolution levels of the MSR. Usually the MSR is interpreted as viewing a signal under different resolutions, as a microscope does (see, e.g., Schmidt, 2012; Schmidt et al., 2015a; Schmidt, 2015b; Liang, 2017). In all of the aforementioned studies, the MSR is based on a regional 2-D representation of VTEC in terms of tensor products of polynomial B-spline functions only. Within this study, however, we apply the MSR for a global 2-D representation of VTEC in terms of tensor products of polynomial and trigonometric B-spline functions, as described by Lyche and Schumaker (2001) and Schumaker and Traas (1991).

3.1 Pyramid algorithm

Neglecting the time dependency, the B-spline approach Eq. (18) reads

$$\begin{array}{}\text{(36)}& f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })={\mathit{\varphi }}_{J_{\mathrm{1}}}^T\left(\mathit{\phi }\right){\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right).\end{array}$$

In the context of the MSR the vectors ${\mathit{\varphi }}_{J_{\mathrm{1}}}\left(\mathit{\phi }\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)$ are called scaling vectors, and the elements $d_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}$ of the matrix ${\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}$ are denoted as scaling coefficients.

With $J_{\mathrm{1}}^{\prime }=J_{\mathrm{1}}-J,J_{\mathrm{2}}^{\prime }=J_{\mathrm{2}}-J$ and $\mathrm{0}<J\le \min (J_{\mathrm{1}},J_{\mathrm{2}})$ we obtain the 2-D MSR of the target function f(x) introduced in Eq. (7) as

$$\begin{array}{}\text{(37)}& f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })=f_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}^{\prime }}(\mathit{\phi },\mathit{\lambda })+\sum _{j=\mathrm{1}}^J\sum _{\mathit{\vartheta }=\mathrm{1}}^{\mathrm{3}}{g}_{J_{\mathrm{1}}-j,J_{\mathrm{2}}-j}^{\mathit{\vartheta }}(\mathit{\phi },\mathit{\lambda }).\end{array}$$

Following the argumentation of Schmidt et al. (2015a) but considering the polynomial and the trigonometric B-spline functions the low-passed-filtered level $(J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}^{\prime })$ signal $f_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}^{\prime }}(\mathit{\phi },\mathit{\lambda })$ and the band-pass-filtered level $(J_{\mathrm{1}}-j,J_{\mathrm{2}}-j)$ detail signals $g_{J_{\mathrm{1}}-j,J_{\mathrm{2}}-j}^{\mathit{\vartheta }}(\mathit{\phi },\mathit{\lambda })$ can be computed via the following relations:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{f}_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}^{\prime }}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\varphi }}_{J_{\mathrm{1}}^{\prime }}^T\left(\mathit{\phi }\right){\mathbf{D}}_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}^{\prime }}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}^{\prime }}\left(\mathit{\lambda }\right),\\ g_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{1}}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\varphi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right){\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{j_{\mathrm{2}}-\mathrm{1}}\left(\mathit{\lambda }\right),\\ g_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{2}}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\psi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right){\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{j_{\mathrm{2}}-\mathrm{1}}\left(\mathit{\lambda }\right),\\ g_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{3}}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\psi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right){\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{j_{\mathrm{2}}-\mathrm{1}}\left(\mathit{\lambda }\right),\end{array}\end{array}$$

where we introduced the definitions $j_{\mathrm{1}}=J_{\mathrm{1}}-j+\mathrm{1}$ and $j_{\mathrm{2}}=J_{\mathrm{2}}-j+\mathrm{1}$ for $j=\mathrm{1},\mathrm{\dots },J$. Herein, the $K_{j_{\mathrm{1}}-\mathrm{1}}\times \mathrm{1}$ and $K_{j_{\mathrm{2}}-\mathrm{1}}\times \mathrm{1}$ scaling vectors ${\mathit{\varphi }}_{j_{\mathrm{1}}-\mathrm{1}}\left(\mathit{\phi }\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{j_{\mathrm{2}}-\mathrm{1}}\left(\mathit{\phi }\right)$ as well as the $L_{j_{\mathrm{1}}-\mathrm{1}}\times \mathrm{1}$ and $L_{j_{\mathrm{2}}-\mathrm{1}}\times \mathrm{1}$ wavelet vectors ${\mathit{\psi }}_{j_{\mathrm{1}}-\mathrm{1}}\left(\mathit{\phi }\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{j_{\mathrm{2}}-\mathrm{1}}\left(\mathit{\lambda }\right)$ can be calculated by means of the two-scale relations

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\mathit{\varphi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right)& ={\mathit{\varphi }}_{j_{\mathrm{1}}}^T\left(\mathit{\phi }\right){\mathit{P}}_{j_{\mathrm{1}}},\\ {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{j_{\mathrm{2}}-\mathrm{1}}^T\left(\mathit{\lambda }\right)& ={\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{j_{\mathrm{2}}}^T\left(\mathit{\lambda }\right){\stackrel{\mathrm{̃}}{\mathit{P}}}_{j_{\mathrm{2}}},\\ {\mathit{\psi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right)& ={\mathit{\varphi }}_{j_{\mathrm{1}}}^T\left(\mathit{\phi }\right){\mathit{Q}}_{j_{\mathrm{1}}},\\ {\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{j_{\mathrm{2}}-\mathrm{1}}^T\left(\mathit{\lambda }\right)& ={\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{j_{\mathrm{2}}}^T\left(\mathit{\lambda }\right){\stackrel{\mathrm{̃}}{\mathit{Q}}}_{j_{\mathrm{2}}}\end{array}\end{array}$$

with $L_{j_{\mathrm{1}}-\mathrm{1}}=K_{j_{\mathrm{1}}}-K_{j_{\mathrm{1}}-\mathrm{1}}$ and $L_{j_{\mathrm{2}}}=K_{j_{\mathrm{2}}}-K_{j_{\mathrm{2}}-\mathrm{1}}$.

The numerical entries of the $K_{j_{\mathrm{1}}}\times K_{j_{\mathrm{1}}-\mathrm{1}}$ matrix ${\mathit{P}}_{j_{\mathrm{1}}}$ and the $K_{j_{\mathrm{1}}}\times L_{j_{\mathrm{1}}-\mathrm{1}}$ matrix ${\mathit{Q}}_{j_{\mathrm{1}}}$ can be taken from Stollnitz et al. (1995b) or Zeilhofer (2008); the corresponding entries of the $K_{j_{\mathrm{2}}}\times K_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\stackrel{\mathrm{̃}}{\mathit{P}}}_{j_{\mathrm{2}}}$ and the $K_{j_{\mathrm{2}}}\times L_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\stackrel{\mathrm{̃}}{\mathit{Q}}}_{j_{\mathrm{2}}}$ are provided by Lyche and Schumaker (2001).

In Eq. (38) we introduced the $K_{j_{\mathrm{1}}-\mathrm{1}}\times K_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\mathbf{D}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}$ of scaling coefficients $d_{k_{\mathrm{1}},k_{\mathrm{2}}}^{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}$ as well as the $K_{j_{\mathrm{1}}-\mathrm{1}}\times L_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{1}}$, the $L_{j_{\mathrm{1}}-\mathrm{1}}\times K_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{2}}$ and the $L_{j_{\mathrm{1}}-\mathrm{1}}\times L_{j_{\mathrm{2}}-\mathrm{1}}$ matrix ${\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{3}}$ of wavelet coefficients. These four matrices can be calculated via the 2-D downsampling equation

$$\begin{array}{}\text{(40)}& {\displaystyle }\begin{array}{l}\left[\begin{array}{ll}{\mathbf{D}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}& {\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{1}}\\ {\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{2}}& {\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},j_{\mathrm{2}}-\mathrm{1}}^{\mathrm{3}}\end{array}\right]\\ =\left[\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathit{P}}}_{j_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathit{Q}}}_{j_{\mathrm{1}}}\end{array}\right]{\mathbf{D}}_{j_{\mathrm{1}},j_{\mathrm{2}}}\left[\begin{array}{ll}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{P}}}}_{j_{\mathrm{2}}}^T& {\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{Q}}}}_{j_{\mathrm{2}}}^T\end{array}\right],\end{array}\end{array}$$

also known as the 2-D pyramid algorithm. The $K_{j_{\mathrm{1}}-\mathrm{1}}\times K_{j_{\mathrm{1}}}$ matrix ${\stackrel{\mathrm{‾}}{\mathit{P}}}_{j_{\mathrm{1}}}$, the $K_{j_{\mathrm{2}}-\mathrm{1}}\times K_{j_{\mathrm{2}}}$ matrix ${\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{P}}}}_{j_{\mathrm{2}}}$, the $L_{j_{\mathrm{1}}-\mathrm{1}}\times K_{j_{\mathrm{1}}}$ matrix ${\stackrel{\mathrm{‾}}{\mathit{Q}}}_{j_{\mathrm{1}}}$ and the $L_{j_{\mathrm{2}}-\mathrm{1}}\times K_{j_{\mathrm{2}}}$ matrix ${\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{Q}}}}_{j_{\mathrm{2}}}$ can be computed via the identities

$$\begin{array}{}\text{(41)}& {\displaystyle }\left[\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathit{P}}}_{j_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathit{Q}}}_{j_{\mathrm{1}}}\end{array}\right]& {\displaystyle }={\left[\begin{array}{ll}{\mathit{P}}_{j_{\mathrm{1}}}& {\mathit{Q}}_{j_{\mathrm{1}}}\end{array}\right]}^{-\mathrm{1}},\text{(42)}& {\displaystyle }\left[\begin{array}{l}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{P}}}}_{j_{\mathrm{2}}}\\ {\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{Q}}}}_{j_{\mathrm{2}}}\end{array}\right]& {\displaystyle }={\left[\begin{array}{ll}{\stackrel{\mathrm{̃}}{\mathit{P}}}_{j_{\mathrm{2}}}& {\stackrel{\mathrm{̃}}{\mathit{Q}}}_{j_{\mathrm{2}}}\end{array}\right]}^{-\mathrm{1}}\end{array}$$

(see, e.g., Schmidt, 2007). The 2-D pyramid algorithm based on the decomposition Eq. (37) is visualized in Fig. 5. The “zeroth” step transforms the observations $y({\mathit{x}}_{i_j},t_j)$ according to Eqs. (13) and (18) into the elements of the scaling matrix ${\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_j\right)$ as introduced in Eq. (17). The procedure applied will be explained in Sect. 3.2.

Figure 5 A 2-D MSR of the signal $f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })$.

The previously described MSR refers to a successive low-pass filtering of the target function $f(\mathit{\phi },\mathit{\lambda },t)$ into two directions, namely latitude φ and longitude λ, in the same manner. If a signal $f(\mathit{\phi },\mathit{\lambda },t)$ is represented according to Eq. (18) up to the level values J₁ with respect to latitude and J₂ with respect to longitude, i.e., $f(\mathit{\phi },\mathit{\lambda },t)\approx f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda },t)$, the application of the MSR Eq. (37) allows for the computation of low-pass-filtered signal approximations up to the level pairs $(J_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}-\mathrm{1})$, $(J_{\mathrm{1}}-\mathrm{2},J_{\mathrm{2}}-\mathrm{2})$,… . The principal structures of the ionospheric key parameters such as VTEC, however, are usually parallel to the geomagnetic Equator. Consequently, we will additionally deal with a 1-D MSR of the signal $f(\mathit{\phi },\mathit{\lambda },t)$ with respect to the latitude. In this case Eq. (37) reduces to

$$\begin{array}{}\text{(43)}& f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })=f_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })+\sum _{j=\mathrm{1}}^J{g}_{J_{\mathrm{1}}-j,J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda }).\end{array}$$

Thus, signal approximations up to the level pairs $(J_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}})$, $(J_{\mathrm{1}}-\mathrm{2},J_{\mathrm{2}})$, … are obtained. From the four relations in Eq. (38) only the first and the third one have to be considered within the 1-D MSR Eq. (43), namely

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{f}_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\varphi }}_{J_{\mathrm{1}}^{\prime }}^T\left(\mathit{\phi }\right){\mathbf{D}}_{J_{\mathrm{1}}^{\prime },J_{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right),\\ g_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })& ={\mathit{\psi }}_{j_{\mathrm{1}}-\mathrm{1}}^T\left(\mathit{\phi }\right){\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left(\mathit{\lambda }\right)\end{array}\end{array}$$

with $j_{\mathrm{1}}=J_{\mathrm{1}}-j+\mathrm{1}$ for $j=\mathrm{1},\mathrm{\dots },J$, $\mathrm{0}<J\le J_{\mathrm{1}}$ and $J_{\mathrm{1}}^{\prime }=J_{\mathrm{1}}-J$. The $K_{j_{\mathrm{1}}-\mathrm{1}}\times K_{J_{\mathrm{2}}}$ matrix ${\mathbf{D}}_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}$ of scaling coefficients and the $L_{j_{\mathrm{1}}-\mathrm{1}}\times K_{J_{\mathrm{2}}}$ matrix ${\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}^{\mathrm{2}}$ of wavelet coefficients can be calculated from the 1-D downsampling equation

$$\begin{array}{}\text{(45)}& \left[\begin{array}{l}{\mathbf{D}}_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}\\ {\mathbf{C}}_{j_{\mathrm{1}}-\mathrm{1},J_{\mathrm{2}}}^{\mathrm{2}}\end{array}\right]=\left[\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathbf{P}}}_{j_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathbf{Q}}}_{j_{\mathrm{1}}}\end{array}\right]{\mathbf{D}}_{j_{\mathrm{1}},J_{\mathrm{2}}},\end{array}$$

where the matrices ${\stackrel{\mathrm{‾}}{\mathbf{P}}}_{j_{\mathrm{1}}}$ and ${\stackrel{\mathrm{‾}}{\mathbf{Q}}}_{j_{\mathrm{1}}}$ can be computed via Eq. (41). The 1-D pyramid algorithm based on the decomposition Eq. (43) is visualized in Fig. 6.

Figure 6 A 1-D MSR of the signal $f_{J_{\mathrm{1}},J_{\mathrm{2}}}(\mathit{\phi },\mathit{\lambda })$ with respect to the latitude φ. The “zeroth” step on the left-hand side conforms with the one in Fig. 5.

Besides the representation of a signal, e.g., VTEC, by means of approximations on different resolution levels with respect to latitude and longitude, the MSR also allows for the utilization of a powerful “data compression procedure”, as the numerical value of a large number of wavelet coefficients is generally close to zero depending on the signal structure (see, e.g., Zeilhofer, 2008).

3.2 Estimation of unknown model parameters

To estimate the elements of the unknown $K_{J_{\mathrm{1}}}\times K_{J_{\mathrm{2}}}$ matrix ${\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right)$ from VTEC observations $y({\mathit{x}}_{i_s},t_s)$ (cf. Eq. 8) within the “zeroth” step of the MSR we apply Kalman filtering according to Erdogan et al. (2017).

In the linear formulation the Kalman filter consists (1) of the state equation

$$\begin{array}{}\text{(46)}& {\mathit{\beta }}_s={\mathbf{F}}_s{\mathit{\beta }}_{s-\mathrm{1}}+{\mathit{w}}_{s-\mathrm{1}}\end{array}$$

and (2) of the system

$$\begin{array}{}\text{(47)}& {\mathit{y}}_s+{\mathit{e}}_s={\mathbf{A}}_s{\mathit{\beta }}_s\end{array}$$

of observation equations. In Eq. (46 ) the u×1 vector ${\mathit{\beta }}_s=\text{vec}{\mathbf{D}}_{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right)$ – known as the state vector – of the $u=K_{J_{\mathrm{1}}}\cdot K_{J_{\mathrm{2}}}$ unknown scaling coefficients at time moment t_s is predicted from the state vector β_s−1 of the previous time moment t_s−1 by means of the u×u transition matrix F_s and the u×1 vector w_s−1 of the process noise. In Eq. (47) ${\mathit{y}}_s=\left(y\right({\mathit{x}}_{i_s},t_s\left)\right)$ and ${\mathit{e}}_s=\left(e\right({\mathit{x}}_{i_s},t_s\left)\right)$ are the I_s×1 vectors of the observations and the measurement errors, respectively; the (i_s)th row vector ${\mathit{a}}_{i_s}^T$ of the I_s×u coefficient matrix A_s is given by the expression

$$\begin{array}{}\text{(48)}& {\mathit{a}}_{i_s}={\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{J_{\mathrm{2}}}\left({\mathit{\lambda }}_{i_s}\right)\otimes {\mathit{\varphi }}_{J_{\mathrm{1}}}\left({\mathit{\phi }}_{i_s}\right),\end{array}$$

as introduced in Eq. (18). The measurement error vector e_s and the vector w_s of the process noise are assumed to be white noise vectors with expectation values E(e_s)=0 and E(w_s)=0, and fulfill the requirements

$$\begin{array}{}\text{(49)}& \begin{array}{rl}E\left({\mathit{w}}_s{\mathit{w}}_{\mathrm{l}}^T\right)& ={\mathbf{\Sigma }}_{\mathit{w}}{\mathit{\delta }}_{s,l},E\left({\mathit{e}}_s{\mathit{e}}_{\mathrm{l}}^T\right)={\mathbf{\Sigma }}_y{\mathit{\delta }}_{s,l},\\ E\left({\mathit{w}}_s{\mathit{e}}_{\mathrm{l}}^T\right)& =\mathbf{0},\end{array}\end{array}$$

where δ_s,l is the delta symbol that equals 1 for s=l and 0 for s≠l. In Eq. (49) Σ_y and Σ_w are given covariance matrices of the observations and the process noise, respectively.

The solution of the estimation problem as defined in Eqs. (46) and (47) generally consists of the sequential application of a prediction step (time update) and a correction step (measurement update). In the prediction step, the estimated state vector ${\widehat{\mathit{\beta }}}_{s-\mathrm{1}}$ and its covariance matrix $\widehat{D}\left({\widehat{\mathit{\beta }}}_{s-\mathrm{1}}\right)={\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s-\mathrm{1}}$ are propagated from the time moment t_s−1 to the next time moment t_s by means of

$$\begin{array}{}\text{(50)}& {\displaystyle }{\mathit{\beta }}_s^{-}& {\displaystyle }={\mathbf{F}}_s{\widehat{\mathit{\beta }}}_{s-\mathrm{1}},\text{(51)}& {\displaystyle }{\mathbf{\Sigma }}_{\mathit{\beta },s}^{-}& {\displaystyle }={\mathbf{F}}_s{\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s-\mathrm{1}}{\mathbf{F}}_s^T+{\mathbf{\Sigma }}_w,\end{array}$$

where the symbol “–” indicates the predicted quantities. The prediction step is followed by the measurement update

$$\begin{array}{}\text{(52)}& {\displaystyle }{\widehat{\mathit{\beta }}}_s& {\displaystyle }={\mathit{\beta }}_s^{-}+{\mathit{K}}_s({\mathit{y}}_s-{\mathbf{A}}_s{\mathit{\beta }}_s^{-}),\text{(53)}& {\displaystyle }{\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s}& {\displaystyle }=(\mathbf{I}-{\mathit{K}}_s{\mathbf{A}}_s){\mathbf{\Sigma }}_{\mathit{\beta },s}^{-},\end{array}$$

where ${\widehat{\mathit{\beta }}}_s$ and ${\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s}$ are the updated state vector and its covariance matrix, respectively. In Eqs. (52) and (53) the u×I_s Kalman gain matrix

$$\begin{array}{}\text{(54)}& {\mathit{K}}_s={\mathbf{\Sigma }}_{\mathit{\beta },s}^{-}{\mathbf{A}}_s^T({\mathbf{A}}_s{\mathbf{\Sigma }}_{\mathit{\beta },s}^{-}{\mathbf{A}}_s^T+{\mathbf{\Sigma }}_y{)}^{-\mathrm{1}}\end{array}$$

behaves as a weighting factor between the new measurements and the predicted state vector. The chosen step size $t_s-t_{s-\mathrm{1}}$ within the KF determines the maximum temporal resolution of the output.

Using the estimations ${\widehat{\mathit{\beta }}}_s$ and ${\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s}$ from Eqs. (52) and (53), a V×1 vector f_s of function values $f({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)$ at arbitrary locations P(φ_i,λ_k) with i=1, …, I, k=1, …, K and $V=I\cdot K$ can be estimated by

$$\begin{array}{}\text{(55)}& {\displaystyle }{\widehat{\mathit{f}}}_s& {\displaystyle }={\stackrel{\mathrm{‾}}{\mathbf{A}}}_s{\widehat{\mathit{\beta }}}_s,\text{(56)}& {\displaystyle }{\widehat{\mathbf{\Sigma }}}_{f,s}& {\displaystyle }={\stackrel{\mathrm{‾}}{\mathbf{A}}}_s^T{\widehat{\mathbf{\Sigma }}}_{\mathit{\beta },s}{\stackrel{\mathrm{‾}}{\mathbf{A}}}_s,\end{array}$$

where ${\widehat{\mathbf{\Sigma }}}_{f,s}$ is the estimated V×V covariance matrix of the estimation ${\widehat{\mathit{f}}}_s$. The V×u matrix ${\stackrel{\mathrm{‾}}{\mathbf{A}}}_s$ is set up in a similar way to matrix A_s in Eq. (47) with Eq. (48). In the following, we will interpret the function values $f({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)=\mathrm{VTEC}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)$ as VTEC values.

3.3 B-Spline model output

The previously explained procedure allows for the dissemination of two products to the users:

• Product 1: a set of estimated scaling coefficients

  $$\begin{array}{}\text{(57)}& {\widehat{d}}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right){|}_{k_{\mathrm{1}}=\mathrm{0},\mathrm{\dots },K_{J_{\mathrm{1}}}-\mathrm{1},k_{\mathrm{2}}=\mathrm{0},\mathrm{\dots },K_{J_{\mathrm{2}}}-\mathrm{1}}\end{array}$$

  from Eq. (52) at time moments t_s for level values J₁ and J₂ at the spatial positions k₁ in the latitudinal direction and k₂ in the longitudinal direction, respectively, as well as their estimated standard deviations

  $$\begin{array}{}\text{(58)}& {\widehat{\mathit{\sigma }}}_{d;k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right){|}_{k_{\mathrm{1}}=\mathrm{0},\mathrm{\dots },K_{J_{\mathrm{1}}}-\mathrm{1},k_{\mathrm{2}}=\mathrm{0},\mathrm{\dots },K_{J_{\mathrm{2}}}-\mathrm{1}}\end{array}$$

  extracted from the covariance matrix Eq. (53).

• Product 2: estimated VTEC values given as

  $$\begin{array}{}\text{(59)}& {\widehat{\mathrm{VTEC}}}_{J_{\mathrm{1}},J_{\mathrm{2}}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s){|}_{i=\mathrm{1},\mathrm{\dots },I,k=\mathrm{1},\mathrm{\dots },K}\end{array}$$

  according to Eq. (55) at time moments t_s for level values J₁ and J₂ in the latitudinal and longitudinal direction, respectively, calculated at grid points P(φ_i,λ_k) as well as their estimated standard deviations

  $$\begin{array}{}\text{(60)}& {\widehat{\mathit{\sigma }}}_{\mathrm{VTEC}}^{J_{\mathrm{1}},J_{\mathrm{2}}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s){|}_{i=\mathrm{1},\mathrm{\dots },I,k=\mathrm{1},\mathrm{\dots },K},\end{array}$$

  extracted from the covariance matrix Eq. (56).

According to Eq. (35), the time interval ΔT between two consecutive maps of the coefficients Eq. (57) and their standard deviations Eq. (58) or the VTEC grid values Eq. (59) and their standard deviations Eq. (60) at times t_s and t_s+ΔT can be chosen arbitrarily, e.g., as 10 or 15 min, 1 h or 2 h. Following Eq. (34), the coordinates φ_i and λ_k of all of the V grid points P(φ_i,λ_k) are defined as ${\mathit{\phi }}_i=-\mathrm{90}{}^{\circ }+(i-\mathrm{1})\cdot \mathrm{\Delta }\mathrm{\Phi }$ with $\mathrm{\Delta }\mathrm{\Phi }=\mathrm{180}{}^{\circ }/(I-\mathrm{1})$ and ${\mathit{\lambda }}_k=\mathrm{0}{}^{\circ }+(k-\mathrm{1})\cdot \mathrm{\Delta }\mathrm{\Lambda }$ with $\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{360}{}^{\circ }/K$. As previously mentioned, the spatial resolution intervals ΔΦ and ΔΛ are usually chosen as 1^∘, 2.5^∘ or 5^∘, i.e., I=181, 73, 37 and K=360, 144, 72.

The two products, i.e., the set of coefficients or the VTEC grid values reflect the two strategies of dissemination. In case of a SH expansion for RT applications as introduced in Sect. 2.1 the corresponding SH coefficients c_n,m from Eq. (10) can be transferred to the user by means of a RTCM (Radio Technical Commission for Maritime services) standard 1264 message. This message allows for the consideration of SH coefficients, but only up to degree n=16. In the case of the B-spline expansion Eq. (13), however, an encoder procedure for the B-spline coefficients Eq. (57) is necessary, because the user has to evaluate the B-spline model just as in the SH case by substituting the B-spline tensor product Eq. (14) for the SHs Eq. (11). Due to the two restrictions, namely the sole use SH expansions and only up to a maximum degree n_max=16, the RTCM message format for data dissemination has to be urgently discussed and must be set up in a more flexible way (refer to the comments in Sect. 5). To apply the RTCM format in its current form, the VTEC grid values Eq. (59) can alternatively be used as observations $y({\mathit{x}}_{i_s},t_s)$ in Eq. (10) to calculate SH coefficients c_n,m(t_s) by means of a least-squares estimation. This way each GIM can be sent at a high update rate to the user for RT applications.

For Product 2, the VTEC grid values Eq. (59) as well as there standard deviations Eq. (60) are disseminated as VTEC and standard deviation maps, i.e., GIMs, with given spatial resolutions ΔΦ and ΔΛ in the latitudinal and longitudinal direction, respectively, in IONEX format to the user.

4 Numerical investigations

In the following, the described modeling approach developed at DGFI-TUM is applied to real data. To be more specific, we use GPS and GLONASS NRT data in hourly blocks and apply ultra-rapid orbits. A detailed explanation of the data preprocessing and the setup of the full observation equations is presented by Erdogan et al. (2017). The IGS IAACs provide final products based on post-processed GNSS observations and orbits with a latency of more than 1 week. Several IAACs additionally provide rapid products with a latency of 1 d using rapid orbits. An overview on the products used in the this paper is given in Table 2.

Table 2 List of GIM products used in this paper. Information on names, types and latencies are taken from the following references: (1) Roma-Dollase et al. (2017), (2) Orus et al. (2005) and (3) this paper.

For the evaluation of the data we have to define an appropriate coordinate system. Here we follow the standard procedure and use a sun-fixed geomagnetic coordinate system. To be more specific, we identify the coordinate system Σ_E introduced in the context of Eq. (3) with the Geocentric Solar Magnetic (GSM) coordinate system (see, e.g., Laundal and Richmond, 2017). Consequently, the SH and B-spline theory as presented in the previous sections is applied in the orthogonal GSM system. As diurnal variations of the ionosphere are mitigated in this coordinate system, the transition matrix F_s introduced in the state Eq. (46) of the KF can be set to the identity matrix I, i.e., F_s=I. In other words, the dynamic system of the KF is set to a random walk process. Furthermore, for the time update in Eq. (46) we fix the step size $t_s-t_{s-\mathrm{1}}$ to 5 min.

While the scaling coefficients Eq. (57) and their standard deviations Eq. (58) of Product 1 are located within the GSM system, the VTEC values Eq. (59) and their standard deviations Eq. (60) of Product 2 are provided in the aforementioned IONEX format on a regular grid defined in a geographical geocentric Earth-fixed coordinate system. Thus, a coordinate system transformation has to be interposed.

4.1 Validation procedure

For validation purposes we rely on the dSTEC analysis which is currently regarded as the standard method for the quality assessment of VTEC models (see, e.g., Orus et al., 2007; Rovira-Garcia et al., 2015).

This analysis method is based on the calculation of the difference between STEC observations STEC$({\mathit{x}}^S,{\mathit{x}}_R,t_s)$ at discrete time moments t_s according to Eq. (2) and a reference observation STEC$({\mathit{x}}^S,{\mathit{x}}_R,t_{\mathrm{ref}})$ along a specified satellite arc as

$$\begin{array}{}\text{(61)}& \begin{array}{rl}{\mathrm{dSTEC}}_{\mathrm{obs}}\left(t_s\right)& =\mathrm{STEC}({\mathit{x}}^S,{\mathit{x}}_R,t_s)\\ & -\mathrm{STEC}({\mathit{x}}^S,{\mathit{x}}_R,t_{\mathrm{ref}}).\end{array}\end{array}$$

The reference time moment t=t_ref is usually referred to the observation with the smallest zenith angle z=z_ref. In the same manner, the differences

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{\mathrm{dSTEC}}_{\mathrm{map}}\left(t_s\right)& =M\left(z_s\right)\cdot \mathrm{VTEC}({\mathit{x}}_{\mathrm{IPP}},t_s)\\ & -M\left(z_{\mathrm{ref}}\right)\cdot \mathrm{VTEC}\left({\mathit{x}}_{\mathrm{IPP},t_{\mathrm{ref}}}\right)\end{array}\end{array}$$

are calculated by means of Eq. (5) from the VTEC map to be validated. The quality assessment is performed by studying the differences

$$\begin{array}{}\text{(63)}& {\displaystyle }\mathrm{dSTEC}\left(t_s\right){\displaystyle }={\mathrm{dSTEC}}_{\mathrm{obs}}\left(t_s\right)-{\mathrm{dSTEC}}_{\mathrm{map}}\left(t_s\right)\end{array}$$

with expectation value E(dSTEC(t_s))=0.

4.2 Estimation of B-spline multi-scale products

Figure 7 shows the global distribution of the IPPs related to GNSS VTEC observations $y({\mathit{x}}_{\mathrm{IPP}},t_s)=\mathrm{VTEC}({\mathit{x}}_{\mathrm{IPP}},t_s)$ as introduced in Eq. (13) for 6 September 2017 at 13:00 UT. As the B-spline model is set up in the GSM coordinate system and the scaling coefficients are restricted to the state equation

$$\begin{array}{}\text{(64)}& d_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_s\right)=d_{k_{\mathrm{1}},k_{\mathrm{2}}}^{J_{\mathrm{1}},J_{\mathrm{2}}}\left(t_{s-\mathrm{1}}\right)+w\left(t_{s-\mathrm{1}}\right)\end{array}$$

according to Eq. (46), we select $\mathrm{\Delta }\mathit{\phi }=\mathrm{5}{}^{\circ }$ and $\mathrm{\Delta }\mathit{\lambda }=\mathrm{10}{}^{\circ }$ as appropriate values for the global average sampling interval of the input data as introduced at the end of Sect. 2.1. Consequently, the B-spline levels to J₁=5 and J₂=3 are taken from Table 1.

Figure 7 Global distribution of the IPPs from GPS (red dots) and GLONASS (blue stars) measurements for 6 September 2017, at 13:00 UT.

The covariance matrices Σ_y and Σ_w of the observations and the process noise, respectively, as defined in the formulae of Eq.(49), are set up according to Erdogan et al. (2017). In more detail, Σ_y consists of two diagonal block matrices related to GPS and GLONASS VTEC observations. The relative weighting between the blocks, i.e., between GPS and GLONASS, is performed by manually defined variance factors.

Figure 8a shows with $J_{\mathrm{1}}=\mathrm{5},J_{\mathrm{2}}=\mathrm{3}$, $K_{J_{\mathrm{1}}}={\mathrm{2}}^{J_{\mathrm{1}}}+\mathrm{2}=\mathrm{34}$ and $K_{J_{\mathrm{2}}}=\mathrm{3}\cdot {\mathrm{2}}^{J_{\mathrm{2}}}=\mathrm{24}$, the numerical values of the total $\mathrm{816}=\mathrm{34}\cdot \mathrm{24}$ scaling coefficients are

$$\begin{array}{}\text{(65)}& {\widehat{d}}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{\mathrm{5},\mathrm{3}}\left(t_s\right){|}_{k_{\mathrm{1}}=\mathrm{0},\mathrm{\dots },\mathrm{33},k_{\mathrm{2}}=\mathrm{0},\mathrm{\dots },\mathrm{23}}\end{array}$$

according to Eq. (57), estimated by means of Eq. (52). As the shift values k₁ and k₂ determine the location of the scaling coefficients, they can be plotted. Figure 8b shows the corresponding standard deviations as defined in Eq. (58). A test of significance is performed for each of the scaling coefficients according to Koch (1999).

While Fig. 8a and b show the results of Product 1 in the GSM system, Fig. 8c and d depict the corresponding results of Product 2 in a geographical geocentric coordinate system. With the choices $\mathrm{\Delta }\mathrm{\Phi }=\mathrm{2.5}{}^{\circ }$ and $\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{5.0}{}^{\circ }$ for the grid spacing in the latitudinal and longitudinal directions, respectively, Product 2 provides the VTEC grid values

$$\begin{array}{}\text{(66)}& {\widehat{\mathrm{VTEC}}}_{\mathrm{5},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s){|}_{i=\mathrm{1},\mathrm{\dots },\mathrm{73},k=\mathrm{1},\mathrm{\dots },\mathrm{72}}\end{array}$$

and the corresponding standard deviations ${\widehat{\mathit{\sigma }}}_{\mathrm{VTEC}}^{\mathrm{5},\mathrm{3}}$ from Eqs. (59) and (60).

Figure 8 Estimated scaling coefficients (a) and their standard deviations (b) for level values J₁=5 and J₂=3 within the GSM coordinate system. Estimated VTEC values (c) and their standard deviations (d) as GIMs within a geographical coordinate system; all sets calculated for 6 September 2017 at 13:00 UT.

Note that for the visualization of VTEC and their standard deviations in Fig. 8c and d, we computed function values on a much denser grid using the interpolation formula (34).

From the comparison of Fig. 8a and c, it can be stated that the numerical values of the scaling coefficients directly reflect the signal structure, i.e., the signal energy. This fact is the consequence of the localizing character of the B-spline functions. Figure 8b and d reveal that the standard deviations are generally larger where no or only a few GNSS observations from IGS stations are available, namely over the oceans, e.g., the Southern Atlantic, but also over specific continental regions such as the Sahara and the Amazon region.

Figure 8a, i.e., the plot of the set of scaling coefficients ${\widehat{d}}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{\mathrm{5},\mathrm{3}}$ in Eq. (65) can be interpreted as a visualization of the 34×24 matrix D_5,3 defined in Eq. (17) and displayed in the top-left boxes of Figs. 5 and 6 for the 2-D and the 1-D MSR, respectively. Consequently, Fig. 8a and b are the results of the zeroth step within the pyramid algorithm, as explained in Sect. 3. Applying the first step of the 1-D pyramid algorithm, the downsampling Eq. (44) provides both the 18×24 matrix D_4,3 of estimated scaling coefficients ${\widehat{d}}_{k_{\mathrm{1}},k_{\mathrm{2}}}^{\mathrm{4},\mathrm{3}}$ for the level values J₁=4 and J₂=3 as well as the 16×24 matrix ${\mathbf{C}}_{\mathrm{4},\mathrm{3}}^{\mathrm{2}}$ of estimated wavelet coefficients. Consequently, the definition of Product 2 in Sect. 3.3 can be extended to “Multi-Scale Products 2”:

• ophg: estimations with levels J₁=5, J₂=3

  $$\begin{array}{}\text{(67)}& \begin{array}{r}{\widehat{\mathrm{VTEC}}}_{\mathrm{5},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s),{\widehat{\mathit{\sigma }}}_{\mathrm{VTEC}}^{\mathrm{5},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)\\ \mathrm{\Delta }\mathrm{\Phi }=\mathrm{2.5}{}^{\circ },\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{5.0}{}^{\circ };\end{array}\end{array}$$

• oplg: estimations with levels J₁=4, J₂=3

  $$\begin{array}{}\text{(68)}& \begin{array}{r}{\widehat{\mathrm{VTEC}}}_{\mathrm{4},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s),{\widehat{\mathit{\sigma }}}_{\mathrm{VTEC}}^{\mathrm{4},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)\\ {\widehat{g}}_{\mathrm{4},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s),{\widehat{\mathit{\sigma }}}_{\mathrm{g}}^{\mathrm{4},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)\\ \mathrm{\Delta }\mathrm{\Phi }=\mathrm{2.5}{}^{\circ },\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{5.0}{}^{\circ }.\end{array}\end{array}$$

We denote the two Multi-Scale Products 2 as “ophg” and “oplg”', where the first symbols refer to the OPTIMAP processing software, which was developed within a third-party funded project (see Acknowledgements). The “p” is chosen according to the temporal output sampling ΔT of maps, with “t” for ΔT=10 min, “1” for ΔT=1 h and “2” for ΔT=2 h. The third symbol describes the spectral resolution and is selected as “l” for “low” and “h” for “high”. Finally, the last symbols indicates the model domain and is set to “g” for “global”. Furthermore, we want to mention again, that the “ophg” and “oplg” products are all presented in geographical coordinates.

4.3 Comparison of VTEC maps from B-spline and spherical harmonic expansions

As mentioned in the context of Table 1, the B-spline levels J₁=4 for latitude and J₂=3 for longitude fit best to the highest degree n_max=15 of a SH expansion Eq. (10). To be more specific, we compare the multi-scale product “o1lg” with the product “codg” provided by CODE. “codg” is characterized by a SH expansion up to degree n_max=15 and a time interval ΔT=1 h of two consecutive maps (Schaer, 1999).

Figure 9 shows the VTEC and standard deviation maps for 6 September 2017 at 13:00 UT as well as the difference map between “o1lg” and “codg”. Although the structures of the two VTEC maps are rather similar, the difference map shows deviations of up to ±6 TECU. To judge this amount, a comparison of VTEC GIMs from different IAACs was performed (not shown here). This investigation stated that deviations between individual IAAC products can reach ±10 % or even more. Studying the structures within the difference map no larger systematic patterns are visible and, thus, justify our assumption that the quality of “o1lg” is comparable with the quality of the IAAC products. The standard deviation maps in Fig. 9b and d show different structures that are mainly caused by the application of the different estimation strategies, namely KF (“o1lg”) and least-squares estimation (“codg”).

Figure 9 VTEC maps “codg” (a) and “o1lg” (c) as well as their standard deviation maps (b, d); difference map of the two VTEC maps (e); all data for 6 September 2017 at 13:00 UT.

To numerically assess the comparability we apply the dSTEC analysis described in Sect. 4.1. First we define a network of receiver stations which are used in Eq. (61).

Figure 10 Distribution of the 10 IGS receiver stations used for the dSTEC analysis.

The chosen set should not be used within the computation of the VTEC maps. Fulfilling both requirements at the same time is difficult and, thus, the set of stations shown in Fig. 10 contains both independent stations and stations used simultaneously in all VTEC models. As GNSS measurements are taken along the satellite arcs, the corresponding IPPs are located spatially within a grid cell and temporally between the discrete time moments of the “o1lg” and “codg” products. In order to calculate the VTEC values in Eq. (62), the spatial and temporal interpolation formulae (34) and (35) have to be applied. Figure 11 shows the RMS values of the differences Eq. (63) during the time span between 1 and 30 September 2017, at the chosen receiver stations.

Figure 11 RMS values for the “codg” (green) and “o1lg” (blue) products computed at the 10 receiver stations shown in Fig. 10. The values in parentheses in the legend are the average RMS values over all 10 receiver stations for the entire test period between 1 and 30 September 2017.

As can be seen, the RMS values vary between 0.3 and 1.6 TECU. By comparing the RMS values of “o1lg” with a mean RMS value of 0.80 TECU and “codg” with a mean RMS of 0.77 TECU we can state that the quality of these two products is very close to each other.

The results indicate that the overall quality of the NRT product “o1lg” is comparable with that of the final product “codg” including the developed and implemented preprocessing strategies and steps of the GNSS data (see Table 2).

4.4 Assessment of the multi-scale VTEC products

The two multi-scale VTEC products “ophg” and “oplg” have been introduced in the two equation blocks Eqs. (67) and (68). In what follows, we study them during a solar storm of medium intensity on 8 September 2017 and during the strongest storm of the last 10 years, the prominent St. Patrick storm, which occurred on 17 March 2015. Figure 12a, c, and d show the results of the 8 September 2017 event at 19:00 UT and Fig. 12b, d, and e display the corresponding maps for the St. Patrick storm event on 17 March 2015 at 19:00 UT.

Figure 12 Multi-scale VTEC products for solar storm events: high-resolution VTEC map “ophg” for 8 September 2017 (a) and for 17 March 2015 (b); low-resolution VTEC map “oplg” for 8 September 2017 (c) and for 17 March 2015 (d); panels (e) and (f) show the detail signals introduced in equation block Eq. (68) and computed by means of Eq. (44) for the two solar events.

As already mentioned in the context of Eq. (43), it is expected that the detail signal Eq. (44) is dominated by structures parallel to the geomagnetic Equator. The detail signal $g_{\mathrm{4},\mathrm{3}}({\mathit{\phi }}_i,{\mathit{\lambda }}_k,t_s)$ shown in Fig. 12e and f meets these expectations. Especially during the St. Patrick storm event, the detail signal shows strong signatures. It should be mentioned that a large number of estimated wavelet coefficients collected in the matrix $C_{\mathrm{4},\mathrm{3}}^{\mathrm{2}}$ are characterized by absolute values smaller than a given threshold. The neglect of these coefficients allows for a high data compression rate. Consequently, the number of significant coefficients as the outcome of a MSR would go drastically below the number of scaling coefficients within the set Eq. (57) of Product 1; the reader can get an impression of the number of neglected coefficients by paying attention to the light green and light blue colors in Fig. 12e and f. This advantageous feature of the MSR was not studied within this work but will be applied and published in the future.

Next, we focus on the solar storm during September 2017 and study the temporal sampling intervals of different GIMs. In summary, we distinguish between six products of different spectral content and different sampling intervals.

Figure 13 RMS values for the “o2hg”, “o1hg”, “othg”, “o2lg”, “o1lg” and “otlg” products computed at the 10 receiver stations shown in Fig. 10 during September 2017. The values in parentheses in the legend are the average RMS values over all 10 receiver stations for the entire test period between 1 and 30 September 2017.

Figure 13 depicts the RMS values computed by the dSTEC analysis at the stations shown in Fig. 10. It is assumed that a product with a larger sampling interval ΔT is less accurate than a product with a smaller sampling interval. Consequently, the average RMS values of “o2hg” and “o2lg” are larger than the corresponding values for a shorter sampling interval. Furthermore, it is assumed that RMS values for a product with higher B-spline levels, e.g., “othg”, are smaller than for the corresponding product with lower B-spline level values such as “otlg”. By comparing the corresponding color bars in Fig. 13, i.e., orange (“o2lg”) vs. red (“o2hg”), light blue (“o1lg”) vs. blue (“o1hg”) and green (“otlg”) vs. yellow (“othg”), the aforementioned assumptions are confirmed.

The differences in the RMS values of the first three products, “o2lg”, “o1lg” and “otlg”, are caused by their different sampling intervals. Comparing the mean RMS values of 0.92 and 0.80 TECU for “o2lg” and “o1lg”, respectively, we find a relative improvement of approximately 13.0 %. By decreasing the sampling from ΔT=2 h to ΔT=10 min, a further improvement of 16.3 % can be achieved.

Table 3 Relative improvements (in percentage) for a downsizing of the sampling interval of the “o2lg”, “o1lg”, “otlg”, “o2hg”, “o1hg” and “otlg” products.

Comparing the RMS values 0.90, 0.72 and 0.68 TECU of the “o2hg”, “o1hg” and “othg” products, respectively, we find relative improvements of 20 % and 24.4 % by downsizing the sampling interval from 2 to 1 h and finally to 10 min. A summary of the relative improvements is given in Table 3.

In the next step, we compare the quality of the multi-scale products “ophg” and “oplg” directly. First, we compare “o2lg” with “o2hg” and obtain an improvement of approximately 2.2 % . In the same manner, we compare “o1lg” with “o1hg” and “otlg” with “othg” and find that improvements of 10.0 % and 11.7 % can be achieved, respectively. Table 4 shows the results for the comparison of each pair of products; an improvement is indicated by numbers in bold font, and a deterioration is indicated by numbers in italic font. As a consequence, an increase of the numerical value for level J₁, i.e., the enhancement of the spectral resolution with respect to the latitude yields a significant improvement in the RMS values as long as the temporal sampling ΔT is less than 2 h.

Table 4 Results (in percentage) of the comparisons of the high-resolution products “ophg” with the low-resolution products “oplg”. Positive (bold) numbers mean an improvement, and negative (italic) values represent a reduction in the quality.

From the investigations in Sect. 4.3, it could be concluded that the quality of the “o2lg” product is comparable to the quality of the IAAC products. It can be seen from Table 4 that there is a strong improvement of more than 26 % when using the “othg” product instead of the “o2lg” product. It is worth mentioning that both products are based on the same input data and are spatially related to each other by means of the MSR.

4.5 Assessment of high-resolution VTEC models

As the “othg” product outperforms all other products used in the previous sections we now compare it with UPC's “uqrg” product (Roma-Dollase et al., 2017) which provides smaller values in the relative standard deviation of their dSTEC analysis in comparison with the products of other IAACs. “uqrg” is a rapid product and is provided with a sampling interval ΔT=15 min. As the “NKLG” station is not used in the calculation of “uqrg”, it is excluded from the calculation of the overall RMS value shown in the legend.

Figure 14 RMS values for the “uqrg” and “othg” products computed at 9 IGS receiver stations during September 2017. The values in parentheses in the legend are the average RMS values over all 9 receiver stations for the entire test period between 1 and 30 September 2017.

As can be seen, the RMS values vary between 0.5 and 1.8 TECU but are mostly below 1.0 TECU. The dominant RMS value of “uqrg” at the “CHPI” station reduces its quality significantly. If “CHPI” is neglected, the mean RMS value of “uqrg” decrease to 0.59 TECU. Summarizing these investigations, we can state that the overall quality of the two products is very similar. Considering the fact that “othg” is a NRT product with a latency of less than 3 h, it also outperforms “uqrg” which is a rapid product with a latency of around 1 d (see Table 2). However, for a final assessment, further validation studies have to be performed between the different products.

5 Conclusions and outlook

This paper presents an approach to model VTEC solely from NRT GNSS observations by generating a MSR based on B-splines; the unknown model parameters are estimated by means of an KF. Based on this approach, a number of products have been created which differ both in their spectral content and in their temporal resolution. From our investigations we state that the MSR provides B-spline models comparable to the standard GIMs of the IAACs, mostly based on SH expansions up to degree n_max=15. As the core of the numerical study we compare our results with the most prominent VTEC maps of the IAACs to rate the quality. As the dSTEC analysis is the most frequently used validation method, we abandon a comparison with satellite altimetry products here. To summarize the validation studies, it can be stated that the high-resolution “othg” product outperforms all products used within the selected time span of investigation.

Besides the facts, that our models can handle data gaps due to the utilization of localizing basis functions, the application of a KF to include a dynamic prediction procedure and the use of the MSR to create products of different spectral content at the same time, it should be mentioned that DGFI-TUM's products

• are based on NRT GNSS observations only, i.e., are using input data with a latency of less then 3 h (in contrast, “codg” relies on post-processed data with a latency of larger than 3 weeks, and “uqrg” relies on rapid data with a latency of at least 1 d; cf. Table 2);

• rely on specially developed software modules (cf. Fig. 15), e.g., the preprocessing module using ultra-rapid orbits;

• and can be disseminated to users with a delay of 2–3 h.

Figure 15 DGFI-TUM's processing modules, including (blue boxes) the download and preprocessing module for GNSS observations, the modeling module by means of B-splines, MSR and Kalman filtering (orange boxes) with possible output as Product 1 and Product 2 (yellow boxes) and the validation module.

In general, the dissemination of these products to users can be undertaken in two different ways: based on estimated scaling coefficients (Product 1) or by calculated VTEC grid values (Product 2). For RT applications, however, the dissemination in terms of Product 1 is preferred, in particular the usage of the RTCM format. In the scope of the developments in the recent years, RT applications have become more important, e.g., in unmanned or autonomous vehicle development; thus, the restriction of the RTCM message to allow only for SH coefficients needs to urgently be discussed. Particularly from the point of view that there are also other modeling methods, a modification of the RTCM format would be appropriate. The MSR allows for significant data compression to be obtained due the step-wise downsampling of the scaling coefficients according to the pyramid algorithm. Details represented in the signal $f_{J_{\mathrm{1}},J_{\mathrm{2}}}$ of the zeroth step are stored in wavelet coefficients for the following steps (see Fig. 5). A large number of estimated wavelet coefficients are characterized by absolute values smaller than a given threshold and, thus, most of them can be neglected for the reconstruction of the original signal. Hence, the overall number of scaling and wavelet coefficients can be reduced drastically. Considering this powerful feature of data compression, we propose replacing the scaling coefficients of the highest levels with the significant wavelet coefficients of the lower levels for a definition of an alternative and more appropriate format for data dissemination in terms of Product 1.

The results presented encourage the further development of high accuracy VTEC maps. By extending the models by a fourth dimension, i.e., modeling of the electron density directly, inaccuracies due to the mapping function can be avoided. To model the vertical structure of the electron density, additional observations have to be incorporated, e.g., from DORIS, satellite altimetry and ionospheric radio occultations. This would mitigate the inhomogeneity of the data distribution and, in turn, even higher levels of the B-spline expansion can be chosen.