High-Resolution Ionosphere Corrections for Single-Frequency Positioning

The ionosphere is one of the main error sources in positioning and navigation; thus, information about the ionosphere is mandatory for precise modern Global Navigation Satellite System (GNSS) applications. The International GNSS Service (IGS) and its Ionosphere Associated Analysis Centers (IAAC) routinely provide ionospheric information in terms of global ionosphere maps (final GIM). Typically, these products are modeled using series expansion in terms of spherical harmonics (SHs) with a maximum degree of n = 15 and are based on post processed observations from Global Navigation Satellite Systems (GNSS), as well as final satellite orbits. However, precise applications such as autonomous driving or precision agriculture require real-time (RT) information about the ionospheric electron content with high spectral and spatial resolution. Ionospheric RT-GIMs are disseminated via Ntrip protocol using the SSR VTEC message of the RTCM. This message can be streamed in RT, but it is limited for the dissemination of coefficients of SHs of lower degrees only. It allows the dissemination of SH coefficients up to a degree of n = 16. This suits to most the SH models of the IAACs, but higher spectral degrees or models in terms of B-spline basis functions, voxels, splines and many more cannot be considered. In addition to the SHs, several alternative approaches, e.g., B-splines or Voxels, have proven to be appropriate basis functions for modeling the ionosphere with an enhanced resolution. Providing them using the SSR VTEC message requires a transfer to SHs. In this context, the following questions are discussed based on data of a B-spline model with high spectral resolution; (1) How can the B-spline model be transformed to SHs in order to fit to the RTCM requirements and (2) what is the loss of detail when the B-spline model is converted to SHs of degree of n = 16? Furthermore, we discuss (3) what is the maximum necessary SH degree n to convert the given B-spline model and (4) how can the transformation be performed to make it applicable for real-time applications? For a final assessment, we perform both, the dSTEC analysis and a single-frequency positioning in kinematic mode, using the transformed GIMs for correcting the ionospheric delay. The assessment shows that the converted GIMs with degrees n ≥ 30 coincide with the original B-spline model and improve the positioning accuracy significantly.


Introduction
The delay on electromagnetic signals traveling through the atmosphere is mostly caused by free electrons which are available within the ionosphere between approximately 50 km to 1000 km. In fact, the so called ionospheric delay d ion , which can be approximated better than 99.9%, affects the propagation of GNSS signals between a satellite S and a receiver R, and is one of the largest error sources in positioning and navigation [1]. It can cause errors of several tens of meters in single frequency positioning for a frequency of f = 1575.42 MHz (L 1 GPS carrier frequency) [2][3][4]. The magnitude of the signal delay depends on the frequency f and on the number of free electrons N e , which disturb the propagation of the signal. Using the ionospheric linear-combination (LC) [5], the user of a dual-frequency receiver can determine the so-called Slant Total Electron Content (STEC) as the integral of the total number of electrons acting on the signal at points P(x) along the ray-path, with position vector x = r cos ϕ cos λ, cos ϕ sin λ, sin ϕ T . Note that the coordinate triple (ϕ, λ, r) comprises the latitude ϕ, the longitude λ and the radial distance r within a geocentric coordinate system Σ E . Single-frequency receivers, however, require external information about the state of the ionosphere to increase the accuracy in positioning. The International GNSS Service (IGS) routinely provides ionospheric delay corrections in terms of global ionosphere maps (GIM) representing the Vertical Total Electron Content (VTEC) as the integrated electron density along the height h between the altitude boundaries h 1 and h 2 of the ionosphere. The GIMs are typically based only on GNSS observations and on the assumption of the Single Layer Model (SLM), which allows the transformation assuming a mapping function M(z), which solely depends on the zenith angle z [6,7]. In the SLM, it is assumed that all electrons in the ionosphere are concentrated in an infinitesimal thin spherical layer of radius R = R e + H, where R e is the mean radius of the Earth and H is the layer height above the Earth's surface. In order to represent VTEC as a continuous function as in Equation (2), the mapped VTEC-given at the position x IPP of the so called Ionosphere Pierce Points (IPP), as the intersection point of the ray-path and the single layer-is used as input to different modeling approaches for GIMs. The official IGS GIM is a combination of independent GIMs [3,8], provided by the Ionosphere Associated Analysis Centers (IAAC) [9]. The IAAC, the Center for Orbit Determination in Europe (CODE), the European Space Operations Center of the European Space Agency (ESOC/ESA), the Chinese Academy of Sciences (CAS), the Canadian Geodetic Survey of Natural Resources Canada (EMRG) and the Wuhan University (WHU) generate GIMs using spherical harmonic (SH) series expansions [8,10]. The GIM of the Universitat Politècnica de Catalunya (UPC) is based on a discretization technique in terms of voxels [11][12][13], whereas the GIM of the Jet Propulsion Laboratory (JPL) is based on a spline approach. The OPTIMAP group (see Acknowledgments) uses basis functions in terms of B-splines for the modeling of the GIMs [14][15][16][17].
It is necessary to distinguish between "final", "rapid", "ultra-rapid" or "real-time" GIMs. The classification is based on the underlying input data, see Table 1. Final GIMs, for instance, are usually based on post-processed observations, satellite orbits and receiver and satellite clocks, which are generally available after 2-3 weeks and the final GIMs can thus be provided with the corresponding product latency. In this regard, we classify rapid GIMs with a latency of one day, ultra-rapid GIM with a latency of 2-3 h, near real-time (NRT) with a latency of about 15 min and finally a real-time (RT) GIM-due to processing times-with a latency of some seconds. As it can be seen from Table 1, NRT and RT products are based on the same input data, but the distinction is based on product latency, which differs due to the computational burden. The IAAC's GIMs are disseminated using the IONosphere map EXchange (IONEX) format. The ASCII-based IONEX format was developed and modified by Schaer et al. (1998) [6] and supports the dissemination of VTEC grids as epoch files or daily files. The spatial sampling of the grid points can be chosen arbitrarily but is typically fixed with sampling intervals of ∆Φ = 2.5 • and ∆Λ = 5 • in latitude ϕ and longitude λ [18], respectively. The socalled epoch IONEX file contains one single VTEC map for an arbitrary epoch, while the daily IONEX file provides several VTEC maps at consecutive epochs with a temporal sampling of ∆T = 2 h, 1 h, 15 min or 10 min. The temporal sampling of consecutive maps within one IONEX file is typically fixed with 1-2 h for final products, up to 15 min for rapid-and 10 min for ultra-rapid products [15]. Since IONEX is a grid-based format, it is particularly flexible, allowing the maps to be provided without considering the underlying modeling approach, e.g., SHs, voxels or B-splines. However, with a spatial sampling of 2.5 • × 5 • , there are 5184 grid points to be calculated for each snapshot map and thus, with decreased temporal sampling intervals, the size of the IONEX file increases. For the dissemination of GIMs in RT, it is obvious that the calculation of such VTEC grids and thus, the creation of IONEX, even epoch IONEX is not adequate. The dissemination of GIMs in RT requires datastream-based formats, see for instance Caissy et al. [19].
A more convenient data format is provided by the Radio Technical Commission for Maritime Services (RTCM) and their Real-Time GNSS Data Transmission Standard RTCM 3.0. The format consists of individual State Space Representation (SSR) messages providing corrections of biases, orbits and clocks for each GNSS. The IM201 submessage type valid for all GNSS includes the SSR Ionosphere VTEC Spherical Harmonics corrections [20]. Centre national d'études spatiales (CNES), CAS and UPC developed an RT combination which is based on this SSR VTEC message for an experimental IGS RT-GIM [21]. The SSR VTEC message allows the transmission of SH coefficients up to a maximum degree of n = 16, by means of the Networked Transport of RTCM via Internet Protocol (Ntrip) to the user.
The SSR VTEC message has two drawbacks, namely, the restriction to coefficients of a series expansion in terms of SHs and their limitation to the maximum degree n max = 16. Consequently, GIMs generated with alternative modeling approaches, such as B-splines or voxels, cannot be considered. However, B-splines and voxels have proven to be appropriate candidates for ionosphere modeling as well [11,22,23]. Roma-Dollase et al. (2017) [24] compared the performance of the GIMs of the IAACs and showed that the voxels-based GIM, the 'uqrg', provided by UPC has been performing with higher accuracy. Goss et al. (2019) [15] derived the relation between the B-spline approach and the SHs in the frequency domain and generated GIMs with different spectral resolutions based on a Multi Resolution Representation (MRR) in terms of B-splines. To be more specific, the maximum degree n max of the SH series expansion defines its cutoff frequency and thus, the minimum wavelength which can be represented. In the B-spline case, the so-called levels define the corresponding minimum wavelength. A B-spline-based model was developed, generating high resolution GIMs with a cutoff frequency comparable to n max = 33.
These high-resolution B-spline models do not coincide with the RTCM standards and has not yet been used for single-frequency positioning. Hence, it remains the question, how can high-resolution GIMs improve the positioning and how can they be made available to users. This paper presents a study based on a ultra-rapid and high-resolution GIM, generated by a B-spline model. This GIM is transferred to an SH expansion considering the (1) spectral, spatial and temporal resolution as well as (2) the computational time, which is needed for the transformation. An accuracy assessment based on the dSTEC analysis [24] and an assessment in terms of single-frequency positioning using the open source software RTKLIB [25] is performed.
The paper is outlined as follows; in Section 2, we introduce the different modeling approaches for GIMs. Section 2.4 gives an overview about the currently available data formats for the dissemination of GIMs. In Section 3, a methodology for the transfer is described considering the above mentioned two requirements about the resolutions and computational time. In this regard, we discuss in Sections 3.1 and 3.2 the usage of the so-called Reuter grid [26,27] to generate pseudo-observations for the estimation of SH coefficients in Section 3.3.
The developed approach is numerically tested and assessed in Section 4. A final summary about the applicability of the developed approach and the conclusions found are given in Section 5.

Global Ionosphere Map Generation
Representing the continuous 3-D function VTEC(ϕ, λ, t) as introduced in Equation (2) by means of a series expansion reads with the space-dependent basis functions φ k (x) and the unknown time-dependent series coefficients c k (t). Given the i s = 1, 2, . . . , I s observations y(x i s , t s ) with observation errors e(x i s , t s ) at discrete time moments t s = t 0 + s · ∆t with s ∈ N 0 and temporal sampling ∆t, Equation (4) can be rewritten as Herein, the total number N + 1 of terms depends, according to the sampling theorem on a sphere, on the spatial sampling intervals ∆ϕ in latitude and ∆λ in longitude for all observations at time t s . The truncation error describes the higher order signal parts of VTEC(x i s , t s ) which cannot be modeled by the series expansion (5) and is therefore neglected in the following.

Modeling Coordinate System
VTEC exhibits a strongly time-varying phenomenon. The ionosphere shows seasonal and annual, as well as daily and subdaily variations. These regular variations follow spatially the geomagnetic equator, therefore the modeling of the ionosphere is typically applied in the geocentric solar magnetic (GSM) coordinate system, which results in much slower variations of VTEC and allows for a more precise representation; detailed information about the GSM coordinate system can be found in [28]. Note, in the sequel of this paper, the models are set up in the GSM coordinate system, while plots and figures are represented in the Earth-fixed geographic coordinate system.

Spherical Harmonics Model
Using the SH approach, the observation Equation (5) can be rewritten as where Y n,m (x i s ) = P n,|m| (sin ϕ) · cos mλ if m ≥ 0 sin |m|λ if m < 0 (8) are the SHs of degree n = 0, . . . , n max and order m = −n, . . . , n, with P n,|m| (sin ϕ) being the normalized associated Legendre functions. The maximum degree n max represents the cutoff frequency of the model representation and is usually defined as a measure for the spectral content of a signal on a sphere. According to the sampling theorem on a sphere, the maximum degree n max is determined from and defines the total number N = (n max + 1) 2 of terms in Equation (7) and, thus, the total number of unknown series coefficients c n,m (t). The SHs are basis functions which are globally defined, as can be seen from Equation (8), this means that they are different from zero almost everywhere on the sphere and require a homogeneous distribution of observations [22]. However, the observations y(x i s , t s ) at the IPP positions are usually given with a globally inhomogeneous distribution with small sampling intervals over the continents and larger sampling intervals in oceanic regions, thus, ∆ϕ and ∆λ in (9) have to be interpreted as global average values. In contrast to SHs, B-splines are localizing basis functions, i.e., they are different from zero only in small intervals, see the red colored B-spline functions in Figure 1. This means that for the estimation of unknown series coefficients, each tensor product (11) must be supported by observations to avoid singularities, [14][15][16][17]. In this regard, the B-spline levels J 1 and J 2 have to be chosen carefully. Similar to the calculation of the maximum degree (9), the levels depend on the mean sampling intervals ∆ϕ and ∆λ and the inequalities have to be fulfilled, see [15,22,23]. Note, the higher the levels are chosen, the more spline functions are distributed along the latitude ϕ and the longitude λ and the finer are the structures in terms of wavelengths L ϕ = 2 · ∆ϕ in latitude and L λ = 2 · ∆λ in longitude that can be represented [15,16].

VTEC Products
Given the (N + 1) × 1 vector c s = [ c 0 , c 1 , . . . , c N ] T of estimated series coefficients c k (t s ) for a time moment t s , the estimated values VTEC(ϕ v , λ v , t s ) at arbitrary points with latitudes ϕ v and longitudes λ v , with v = 1, . . . , V, can be calculated by Consequently, two strategies for disseminating ionospheric corrections can be set up. For each one, an appropriate product type is defined, which provides the information about the state of the ionosphere, namely • Product type 1: estimated series coefficients c k (t s ) with k = 0, . . . , N. It is assumed that a user has an appropriate converter for evaluating the respective series expansion (4).
In that way, the corrections can directly be calculated by the user. • Product type 2: VTEC grid values VTEC(ϕ l , λ r , t s ) for l = 1, . . . , L latitudes ϕ l and r = 1, . . . , R longitudes λ r , where L is the number of points along the meridians and R is the number of grid points along the circles of latitudes with arbitrary sampling ∆Φ = ϕ l+1 − ϕ l and ∆Λ = λ r+1 − λ r , respectively. This implies that a user has to interpolate between the provided grid points in order to obtain the required correction.
In the following, we distinguish between grid-based Product type 2 formats (Section 2.5.1) and coefficient-based Product type 1 formats (Section 2.5.2).

Dissemination Formats for Positioning
To improve positioning by GNSS, ionospheric effects on the signals have to be corrected. Ionospheric corrections are therefore typically provided in terms of VTEC maps. The type of application (real-time, post-processing) determines the Product type to be used. For a proper application of the corrections, the preliminary position of the receiver x R and the position of the GNSS satellites x S are required. These allow the determination of the position vector x IPP of the IPP and subsequently the determination of the ionospheric correction either from Product type 1 or from Product type 2.

IONEX Format
The IONEX format was introduced by Schaer et al. in 1998 [6] and is up to now commonly used for the dissemination of GIMs of any type of latency. The format is particularly advantageous, because it allows for provision of additional information such as used receiver stations as well as DCBs. Usually, the VTEC values are given in an Earth-fixed geographic coordinate system on grids with a spatial sampling of ∆Φ = 2.5 • in latitude and ∆λ = 5 • in longitude for snapshot maps with a temporal sampling of ∆T = 2 h.
When providing corrections using grid-based formats, the user must first perform an interpolation between the grid points and the temporal snapshot maps [6,15]. In order to calculate the values VTEC(ϕ IPP = ϕ l + q · ∆Φ, λ IPP = λ r + p · ∆Λ) for the position of the IPPs with 0 ≤ q ≤ 1 and 0 ≤ p ≤ 1 for time t s , the simple bilinear spatial interpolation from VTEC values of the surrounding four points, as well as the temporal interpolation for any arbitrary time t = t s + τ · ∆T with 0 ≤ τ ≤ 1 according to Schaer et al. (1998) [6] can be applied. By applying the interpolation, the quality of the calculated VTEC values depends on the position of VTEC(ϕ l + q · ∆Φ, λ r + p · ∆Λ, t) within the grid cell, on the spatial sampling intervals ∆Φ and ∆Λ, as well as the temporal sampling ∆T. The interpolation is applied in Sun-fixed coordinate system in order to take into account the Earth's rotation, see [6].

Coefficient Based Format
The currently preferred method for dissemination of ionospheric corrections for RT positioning is based on SH coefficients, follows Product type 1 and is based on RTCM streams. The SSR VTEC message allows the dissemination of estimated coefficients c n,m (t s ) of a SH model with additional information about degree n and order m of the underlying SH model, as well as the time moments t s and the height H of the single layer model for which the coefficients are valid. According to the message specifications the coefficients must be provided in Sun-fixed longitude with a phase shift of 2 h to the approximate VTEC maximum at 14:00 local time. This means that a coordinate transformation must be carried out for the values VTEC(ϕ IPP , λ IPP , t s ), which first provides a transformation from GSM to the Sun-fixed geographic coordinate system and then a rotation around the spin-axis of the Earth by 2 h. The SSR VTEC message is designed for the dissemination of SH coefficients and allows them only up to degree and order n = 16. This means that higher resolution models cannot be distributed via the SSR VTEC message. Another restriction of the message is that models not based on SHs require always a transformation to SHs.

Requirements for GIM Formats
The strategy for the dissemination and thus the choice of Product type 1 or Product type 2 of the estimated VTEC information is subject to certain restrictions, depending on their further use, namely: . . .
For precise positioning in post processing mode: to be more specific, as long as the data are not immediately used by a user, the size of the data is less important and the transformation may take longer. Which means, the VTEC information can be prepared with high resolution and quality and be provided to the user. It does not matter which dissemination strategy is chosen and through which platform (FTP, streaming) the information reaches the user. . . .
For precise navigation and positioning in RT: in this case, the selection of possible dissemination formats is limited. Hence, the user must receive the ionospheric information within seconds and with high precision, in order to correct the GNSS measurements used for positioning. Currently, only the SSR VTEC message is used for this purpose. Consequently, the dissemination strategy with Product type 1 has to be chosen.

Methodology
For the dissemination of a B-spline-based RT ionospheric correction to a user by means of Product type 1 and the currently intended data format, i.e., the SSR VTEC message, a transformation of the B-spline model into SH coefficients is required. For the transformation of the B-spline model, it must be considered that: 1. The transformation should be done with high precision, i.e., since B-splines and SHs are characterized by different features, a high degree n max of the SH expansion needs to be considered. 2. The time which is needed for the transformation is limited. Hence, the processing time including the generation of pseudo-observations-the evaluation of Equation (  18)-and the estimation of the SH coefficients (cf. Equation (20)) should be completed within seconds.

Point Distribution on a Sphere
The transformation of the B-spline model to SHs requires input data generated by the B-spline model which reflect the variations in the GIM globally, i.e., VTEC values given as pseudo-observations on a global grid. As a matter of fact, a reliable computation of SH series coefficients requires homogeneously distributed input data. Since a regular grid as described in Section 2.5 does not fulfill this criterion, a Reuter grid is used instead [26]. The Reuter grid provides equi-distributed points P(ϕ l , λ l,r ) on the sphere with spherical coordinates ϕ l = −90 • + l ∆ϕ l and λ l, for r = 0, . . . , γ l − 1 and l = 1, . . . , γ − 1 with γ ∈ N and the sampling intervals ∆ϕ l = 180 • γ in latitude direction. Furthermore, means the total number of grid points along the circle of latitude ϕ l . Reforming of Equation (17) results in the sampling interval ∆λ l of the grid points in longitude direction. For ϕ l = 0, it applies ∆ϕ l = ∆λ l . More information about the Reuter grid can be found in [27,30]. The left part of the Table 2 shows for different values of γ the values of the spatial sampling intervals ∆ϕ l along the meridian and ∆λ l at the equator, as well as the total number V of Reuter grid points. The right part of Table 2 provides, for each degree n max , the necessary mean sampling intervals ∆ϕ and ∆λ from Equation (9) and the corresponding total number of coefficients N.
The Reuter grid is particularly advantageous for the RT conversion of the B-spline model to SHs, because it provides a homogeneous distribution of the pseudo-observations and fulfills at the same time the conditions of the inequalities (9). Furthermore, a small number V of pseudo observations allows a transformation with short processing time. The values ∆ϕ l and ∆λ l for Reuter grids of γ = 16,21,25,31,35 in the left part of Table  2 fulfill the conditions (9) for SHs with maximum degrees n max = 15, 20, 24, 30, 34 in the right part.

Pseudo Observations from B-spline Model Output
The pseudo observations are generated by means of the given of levels J 1 and J 2 with u = K J 1 · K J 2 for the time moment t s . We rewrite Equation (13) for the B-spline case as with the V × u design matrix A s consisting B-spline tensor products according Equation (11). f s is a V × 1 vector comprising the values VTEC R (ϕ l , λ r , t s ), with the index R indicating the VTEC values on a Reuter grid.

Estimation of SH Coefficients from Reuter Grid
The series coefficients c n,m are estimated by parameter estimation from V > N pseudo observationsf s . By rewriting the observation Equation (7), the linear equation system can be established, with the V × 1 consistency vector e s , the (N + 1) × 1 vector c = [c 0,0 , c 1,0 , . . . , c N,N ] T of series coefficients, as well as the V × (N + 1) design matrix X s comprising the functions Y n,m (ϕ l , λ r ) according to (8). In case of a matrix X s with full column rank, the problem can be solved by Substituting Equation (18) for the observations vectorf s the transformation between the B-spline and the SH coefficients is established and reads with the (N + 1) × u transformation matrix T s . With the estimated coefficients c s , the Product type 1 and Product type 2 can be distributed to users (cf. Section 2.4). Figure 3 shows a flowchart on the methodology of the previously derived transformation method. Thereby the generation of pseudo-observations VTEC R (ϕ l , λ r , t s ) on a Reuter grid by means of given B-spline coefficients β s is followed by the estimation of the SH coefficients. The different sets β s and c s of coefficients can be considered for the dissemination in terms of Product type 1, however, only the set c s can be disseminated by means of the SSR VTEC message. Furthermore, both sets allow the generation of Product type 2. The two versions VTEC(ϕ l , λ r , t s ) generated by B-splines and VTEC SH (ϕ l , λ r , t s ) generated by means of the transformation method are depicted on the right hand side of the flowchart in Figure 3. Both grids can be disseminated by means of the IONEX file format.

Results and Discussion
Subsequently, the developed approach will be applied for the period from 2 September 2017 to 12 September 2017. This period was chosen because of the varying ionospheric activity. Hence, the period is in the decreasing phase of the last solar cycle with a moderate number of sunspots. In addition, solar flares occurred between September 4 to September 8 with the consequence of a geomagnetic storm and increased values for the Kp index up to the value 8.
Given are the B-spline coefficients in the Sun-fixed coordinate system of the B-spline series expansion with the levels J 1 = 5 and J 2 = 3 for the mentioned period. At this stage, it should be mentioned that the present B-spline model according to Goss et al. [15] and Erdogan et al. [17] was generated by means of hourly GNSS observations and ultra-rapid GNSS orbits and that it is classified as an ultra-rapid GIM according to the definition from the introduction.
According to [15] and to Table 3, an expansion in terms of polynomial B-splines with the level J 1 = 5 can be identified by a cutoff frequency of n max = 33 in SHs. An expansion in terms of trigonometric B-splines with the level J 2 = 3 corresponds to a cutoff frequency n max = 12. The ionosphere usually shows structures that follow the geomagnetic equator and thus, stronger variations occur in latitude direction. Therefore, a higher spectral resolution is chosen in the latitude direction and a lower spectral resolution in the longitude direction. Since the spectral representation in latitude and longitude directions are very different for the given B-spline model, different cases are considered in the following to test and validate the determined approach. Table 3. Different values for B-spline levels J 1 and J 2 with their corresponding values for the cutoff frequency n max and the minimum wavelengths L ϕ and L λ (see Goss et al., 2019 and 2020 [15,16]).   Thereafter, five different transformation cases are considered and the SH coefficients are estimated by the parameter estimation as described in Section 3.3. Thereby, the requirements for the transformation, which were defined in the beginning of Section 3 were taken into account. The values for γ, as well as for the highest degree n max of the SH series expansion were used according to Section 3.1 and Table 2. Table 4 presents the five cases, which are analyzed in the following for their quality and feasibility. For quality estimation we follow the flowchart shown in Figure 3 and generate both the original B-spline GIM for VTEC(ϕ l , λ r , t s ) and the transformed version VTEC SH (ϕ l , λ r , t s ) on a regular grid. The quality characteristics, RMS (%, TECU), mean value δ mean , maximum value δ max and minimum value δ min are determined based on the deviations, i.e., the differences δ(ϕ l , λ r , t s ) = VTEC(ϕ l , λ r , t s ) − VTEC SH (ϕ l , λ r , t s ) between the original GIM and the SH GIM. The applicability of the approach is executed on the basis of the necessary processing time ∆t per epoch, which is needed for the transformation. For comparison reasons, the transformations have been computed on a workstation with 64 GB RAM, and a 8 core processor of 3.2 GHz clock rate. Thereby ∆t comprises the evaluation time for the pseudo-observations as well as the estimation of the SH coefficients. All values given in Table 4 are averages obtained from transformed version computed with a ∆T = 10 min temporal sampling within the designated time span.

Polynomial B-Splines Trigonometric B-Splines
In the following, we discuss the given cases in more detail. In the first case, with V = 317 pseudo-observations given on a Reuter grid with γ = 16, an SH series expansion with degree n max = 15 and N = 256 coefficients has to be estimated. The average processing time for the transformation took ∆t = 1.43 s. However, the transformed version with degree n max = 15 produces systematic errors compared to the original GIM. The RMS value of the deviations shows with 1.31 TECU significant differences. The relative RMS value provides the RMS in percentage with 9.23% for the first case. The first column in Figure 5 shows the SH GIM in the top panel and the deviation map in the bottom panel. The deviation map represents mainly stripes in east-west direction, which indicate the difference in the spectral resolution for latitude direction between the original GIM corresponding to n max = 33 and representing the minimum wavelengths of L ϕ = 10.9 • and the SH GIM of n max = 15, representing minimum of wavelengths L ϕ = 24 • (cf. Table 3). Visible are deviations with maximum values of δ max = 8.22 TECU and minimum δ min = −8.06 TECU, whereas the average δ mean is close to 0 TECU. In the second case, both the grid resolution and the degree of the SH expansion are increased to γ = 21 and n max = 20, respectively. This causes an insignificant increase in computation time to ∆t = 1.85 s, but a significant improvement in the transformation accuracy to a relative RMS of 5.83% (cf. Table 4). For the third case, the processing time doubles compared to the first case with ∆t = 3.12 s, but the error of the transformation is reduced by a half and a relative RMS of 4.19% can be achieved. This improvement is also visible in the second column of Figure 5. Due to the increased degree of the SHs, finer structures in latitude direction can be represented better than in the first case. Therefore, the stripes in east-west direction in the deviation map show a smaller extension in latitude direction and occur with decreased magnitude of δ max = 5.23 TECU and δ min = −4.21 TECU.   Table 4 for September 8, at 12:00 UT.
Further improvements can be achieved by increasing the degree of the SH expansion for the transformation, but with a further extension of the processing time per epoch. The transformation using n max = 34 needs an average transformation time of ∆t = 8.23 s per epoch. The corresponding SH GIM is shown in the top right panel in Figure 5. The eastwest stripes in the deviation map below are mainly visible in the area of the equatorial anomaly with maximum values of δ max = 1.79 TECU and minimum values of δ min = −1.9 TECU. The SH GIM describes a relative RMS of 1.83% compared to the original GIM. For these five cases, it can be concluded that the quality of the transformation improves with increasing degree n max for the SH series expansion.
Accordingly, Figure 6 shows on the left the relative RMS values and on the right the RMS values of the deviations for the period from 2 September 2017 and 12 September 2017 with decreasing magnitude.

Validation
Subsequently, the original B-spline GIM and its transformed versions are validated using the dSTEC analysis and tested for their ability to correct the ionospheric disturbances in single frequency positioning. Additionally, the GIMs of the IAACs CODE from Berne, Switzerland and UPC from Barcelona, Spain are used for comparison.

dSTEC Analysis
The dSTEC analysis is one of the most commonly used validation method for GIMs. It is based on the comparison of differenced STEC observations dSTEC obs of a receiver and satellite pair along its phase-continuous arc, with the differenced STEC values dSTEC map computed from the GIM to be validated as with expectation value E(dSTEC(t s )) = 0. Note, the STEC values from the GIMs are obtained by applying the mapping from Equations (3) and by the applying the interpolations from Section 2.5.1 according to Schaer et al. (1998) [6] for the position and time of the IPP of the observations along the satellite arc. More detailed information on the dSTEC analysis can be found in [15,18,24]. For the calculation of dSTEC obs , the receiver stations shown in Figure 7 are used. These are either independent of the GIMs to be validated or have contributed to all of their calculations. Table 5 provides a summary of the results of the dSTEC analysis. The original B-spline GIM of DGFI-TUM is called 'othg'. The naming follows the definition in [15], where 'o' stands for the OPTIMAP processing software, which was developed in a thirdparty project (see Acknowledgments). The second digit describes the temporal output sampling with 't' for ten minutes. The 'h' describes the high spectral resolution and the 'g' the global expansion. The transformed versions of the 'othg', coinciding with the different cases of the previous section are named accordingly, but with the respective degree of SHs used for the transformation in the second and third digit. The last two columns in Table 5 show the results for the GIMs 'codg' and 'uqrg' provided by CODE and UPC, respectively. All GIMs used in the dSTEC analysis are given in IONEX format with spatial sampling of ∆ϕ = 2.5 • in latitude and ∆λ = 5 • in longitude. All other specifications for the different GIMs are depicted in Table 5.  The last row shows the RMS values of dSTEC(t s ), Equation (23), for each GIM given as the average over all stations and all observations during the period in September 2017. They indicate the performance of the individual GIMs during the designated period and for the selected stations. As usual the 'uqrg' performs best in the selection of GIMs with an RMS of 0.85 TECU, followed by 'othg' with an RMS of 0.91 TECU. The transformed versions show a trend that was already shown in Figure 6; the RMS values of transformed versions with larger n max approach the value of the original 'othg'. With an RMS of 0.92 TECU, the 'o34g' has a negligible difference in the accuracy to 'othg'. This confirms their approximate agreement in the spectral resolution with n max = 34 for 'o34g' and the cutoff frequency of n max = 33 (cf. Table 3) in the latitude direction of 'othg'. A comparison of the GIMs 'o15g' and 'codg', both based on an SH series expansion of n max = 15, shows that a higher degree for the transformation is necessary, since the original GIM 'othg' provides a higher and 'o15g' a less quality than the 'codg' during the period of investigation. The estimated coordinate components X est (t s ), Y est (t s ) and Z est (t s ) are subtracted from the actual coordinates X(t s ), Y(t s ) and Z(t s ) provided by the IGS [33] and the deviations of the 3D position are determined as (24) Figure 8 shows the time series of the differences S(t s ) for 2 September 2017 (left) and 8 September 2017 (right). A significant difference can be seen in S(t s ) between the two days. On 8 September 2017, due to the high geomagnetic activity the ionospheric corrections are more difficult to determine and thus, the positioning accuracy decreases and the deviations increase. This trend is especially noticeable for the stations BOGT and APSA close to the equator. It can also be observed that their deviations increase during local noon, i.e., between 15:00 and 20:00 UT for BOGT and between 20:00 and 24:00 UT for ASPA. The large deviations of the position of BOGT and ASPA between 00:00 and 05:00 UT on September 8 are due to the high geomagnetic variations with increased geomagnetic index of Kp = 8. The variations are not as significant at the station WTZR, which was located on the night side of the Earth at that time. The daily variations of S in Figure  8 show that 'othg', 'o30g' and 'o34g', in light green, blue and red, respectively, exhibit a similar variation. Stronger deviations and mostly larger values S(t s ) can be recognized for the GIMs 'o15g', 'o20g' and 'o24g' with the dashed lines. The GIM 'codg' shows a similar behavior during the days. However, the 'uqrg' performs better in single frequency positioning and sometimes shows lower values than the 'othg', the 'o30g' and the 'o34g', except at local noon at station WTZR, where it shows large deviations.
A detailed evaluation of the performance of the different GIMs is performed using the respective RMS and the averageS values of the time series S(t s ), which are depicted in Table 6.   The color scheme in Table 6 implies the lowest and highest values of the RMS and S in green and red, for the respective station and day. It can be seen that for all stations there are differences in the positioning accuracy between the two days examined. Hence, for DOY 251, the RMS and average values are increased. There is an additional trend which shows, that for the selected stations and days the high resolution GIMs, 'othg' and 'uqrg' allow a correction of the ionospheric disturbances in a way that leads to a positioning with increased accuracy (see the green colors). Furthermore, the poor performance of the 'o15g'-which has mostly highlighted values in red-confirms the result from Section 4.1.1, that a transformed version a with higher degree of SH series expansion is necessary to achieve the quality of the original GIM. As the degree of SHs for the transformation increases, the accuracy of positioning increases when using their ionosphere corrections. The values written in bold in Table 6 allow the conclusion that a transformation with at least the maximum degree n max = 30 is necessary to achieve the quality of 'othg'.
It should be pointed out, that the pure SH model, the 'codg', can correct the ionospheric disturbances for single frequency positioning better than the transformed version 'o15g'. However, in the example shown here, 'codg' cannot achieve accuracy of the 'othg', the 'uqrg' and the transformed versions 'o30g' and 'o34g'.

Summary and Conclusions
The dissemination of RT ionosphere information is currently based on the RTCM message 1264, which can provide SH coefficients up to a maximum degree of n max = 16. The IGS started to provide a combined RT-GIM, collecting the GIMs from CAS, CNES and UPC using Ntrip protocol and the SSR VTEC message [21]. Since the SSR VTEC message is currently the only format type for dissemination the RT-GIMs of the IAACs, only SH models can be considered. Other models which are not based on SHs need to be converted and suffer from degeneration in the accuracy of the GIMs. Especially, highresolution models cannot be converted appropriately to SHs when a low degree n max = 16 is used in the transformation. This means that an adaptation or extension of the existing dissemination formats, as already objected by Goss et al. (2019) and (2020) [15,16], must be carried out.
This paper presents a novel approach on the transformation of GIMs modeled by means of a B-spline series expansion to SH coefficients. It should be noted that the developed method can also be applied to models based on voxels or alternative basis functions.
The developed approach is setup in a way that the processing time during transformations are acceptable for RT applications and the transformed SH GIMs maintain the quality of the original GIM.
However, the numerical investigations are based on ultra-rapid GIMs modeled by a B-spline series expansion according to [15]. The given ultra-rapid GIMs only serve to test the developed approach and to assess the accuracy and quality of the transformations and the latency can further be decreased. The investigation period covers the days 2 September to 12 September 2017 with a geomagnetic storm on September 8. In a first step, a case study considering five cases with different degrees n max for the transformation to SH is shown. It was found that with an increased value n max ≥ 30, the transformations converge sufficiently with the original B-spline model. However, for a transformation of the given B-spline model with level values of J 1 = 5 and J 2 = 3 a degree of n max > 34 is required.
In the assessment comprising a dSTEC analysis and a single-frequency PPP, the models of CODE (SHs of n max = 15-'codg') and UPC (voxel and Kriging model-'uqrg') were used for comparison with the B-spline GIM 'othg' and its transformed versions. For both assessment methods, the high-resolution GIMs 'othg' and 'uqrg' as well as the transformed versions with degrees n max ≥ 30 performed the best. The single-frequency PPP shows a discrepancy between the transformed SH GIMs of lower degree n max and the original GIM 'othg', but also between 'codg' and the high-resolution solutions 'othg' and 'uqrg'.
Finally, it can be concluded that for a transformation of the B-spline model to SHs, the maximum degree n max has to be adapted accordingly. In this way, high quality ionospheric corrections can be provided to single-frequency users for the correction of positioning and navigation. This means that for the dissemination of high-resolution GIMs, an extension of the limited degree of the SH coefficients which is allowed by SSR VTEC message is required and accomplished urgently.

Data Availability Statement:
The data presented in this study are openly available in ftp://cddis. gsfc.nasa.gov/.