Regional Ionosphere Delay Models Based on CORS Data and Machine Learning

2.1 Selection of Study Region and Data

Figure 2 shows the locations of stations that belong to the CORS (blue dots) and the EPN (red dots) networks whose dual-frequency GPS observations were used in this study to generate VTEC models. Two VTEC models were developed based on the regions covered, namely, RIM IONO_BH and RIM IONO_WB.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

Map documenting the locations of the dual-frequency stations of the CORS (blue dots) and the EPN (red dots) networks whose observations were used for VTEC modeling. Stations used for the RIM validation are highlighted with an inner black dot with a white rim; their names are indicated on the map.

To estimate the RIM IONO_BH, the station network was selected to include nine Bosnia and Herzegovina Positioning Service (BIHPOS) stations each located approximately 80 km from the central EPN station in Sarajevo. An orange ellipse on the map marks this area. The RIM IONO_WB was derived from GPS observations obtained from the following CORS networks: the Albanian GNSS Permanent Stations (AlbGNSS), BIHPOS, the Croatian Positioning System (CROPOS), the Macedonian Positioning System (MAKPOS), and the Slovenia-Geodesy-Navigation-Location (SIGNAL). We also used observations from eight EPN stations within this region. Observations from the CORS networks located in Serbia and Kosovo were not available from the network providers and are not freely available to the public; thus this information was not used in this study. Also, historical data from the CORS network in Montenegro were not preserved for these study periods and thus not available for use in this study. Consequently, CORS observations from these countries were not included in the RIM IONO_WB. The network selected for the RIM IONO_WB included approximately 80 CORS and EPN stations between approximately 40° N to 47° N and 13° E to 23° E. RIM IONO_WB was developed as an independent model that would be applicable for each of the participating countries mentioned above. Therefore, the RIM IONO_WB includes files with the extension “ION” that contain separate sets of ionospheric corrections for each country. The reference point of the Taylor series expansion was set approximately in the mid-range of each country included in the study. Permanent stations that were used to validate the newly-developed RIMs are highlighted on the map with inner black and white circles (Figure 2). These stations were not used for estimating the regional models. The following stations shown on the map were used to validate the RIM IONO_BH: EPN SRJV and BIHPOS SARA (which are very close to one another), BIHPOS SEKO, and, additionally, two stations outside the modeled region, EPN POZE and EPN DUB2. Stations used for the RIM IONO_WB and the RIM IONO_WB_AI validation included EPN GSR1, EPN POZE, EPN SRJV, IGEWE TIRA, and EPN ORID.

2.2 Regional Ionosphere Modeling with Bernese GNSS Software

Bernese GNSS Software version v.5.2 (Dach et al., 2015) was used to process GNSS data, estimate ionosphere models and perform positioning. Before generating ionosphere models, several pre-processing steps must be carried out. Table 1 summarizes routines and the processing steps according to the order of their execution. Additional details on each routine and processing step in the Bernese GNSS software are described by Dach et al. (2015). The ionosphere models were estimated within the IONEST routine using the Bernese-recommended values (Dach et al., 2015) listed in Table 2. This routine works only on the assumption of GPS zero-difference observations.

Ionosphere mapping was performed on the undifferenced (zero-difference) level by analyzing the so-called geometry-free linear combination (L₄) of phase observations which are formed by subtracting observables at different frequencies. This means that all frequency-independent effects such as the satellite-receiver geometrical range, clock errors, and tropospheric delay, among others, were canceled out, while the ionospheric effect remained (Ciraolo et al., 2007). Thus, the geometry-free linear combination includes the ionospheric delay and can be used to estimate the ionosphere model. With this approach, the carrier phases are not fitted to code (pseudo-range) observations; this means that no code-phase leveling was applied.

The observation equation for zero-difference phase observations can be outlined as shown in Equation (4):

4

where L₄ is the geometry-free phase observable, α = 4.03 · 10¹⁷ ms⁻²TECU⁻¹ is a constant, f₁, f₂ are the frequencies associated with the carriers L₁ and L₂, β · ϕ represents the wind-up term associated with the right-handed polarized GPS signal (typically a centimeter-level term), B₄ = λ₁ · B₁ – λ₂ · B₂ is a constant phase bias caused by the initial phase ambiguities (i.e, unknown integer number of cycles) B₁ and B₂ with their corresponding wavelengths λ₁ and λ₂. The initial phase ambiguities B₁ and B₂ are estimated as real-value parameters within the initial least-squares adjustment using the phase observations L₁ and L₂ on both frequencies. The initial phase ambiguity has the same value provided that no loss of signal lock occurs. At least one parameter (B₄) has to be solved for each receiver and satellite pass since it contains phase ambiguities. For cycle slip detection, differences between two satellites for the same epoch are formed from the measurements in the observation files. If cycle slips are detected, the geometry–free linear combination L₄ is checked in order to determine their size in both frequencies. A new ambiguity is set up at the corresponding epoch of the detected cycle slip. Pre-processing options were set in the IONEST routine to define a new ambiguity parameter B₄ for each detected cycle slip. The resulting geometry-free linear combination in Equation (4) allows us to estimate VTEC as a function of geographic latitude ϕ_IPP and the Sun-fixed longitude s_IPP at the intersection point of the line-of-sight between the satellite and the receiver with the ionospheric layer. This intersection point is known as the ionospheric pierce point (IPP).

The regional VTEC model is based on the two-dimensional (2-D) Taylor series expansion shown in Equation (5):

5

where n, m are the degree values of the 2-D Taylor series expansion in geographic latitude and Sun-fixed longitude, n_max, m_max are the maximum degree values of the Taylor series expansion in geographic latitude and Sun-fixed longitude (Table 2), ϕ₀, s₀ denote the coordinates of the origin of the Taylor expansion and c_nm stands for the unknown coefficients of the Taylor series expansion, i.e., the regional ionosphere model parameters to be estimated (Dach et al., 2015; Wild, 1994).

The Sun-fixed longitude s is related to the local solar time (LT) as shown in Equation (6):

6

where UT is an abbreviation for the Universal Time, LT is local time, and λ_IPP denotes the geographical longitude of the IPP; all values are in radians.

The latitude ϕ_IPP and the longitude λ_IPP as shown in Equation (7) and Equation (8):

7

8

where ϕ_k, λ_k are the latitude and longitude of the receiver k, A_z stands for the azimuth from the receiver to the satellite, and ψ is the angle between the lines joining the center of the Earth with the IPP and receiver location. It can be calculated as shown in Equation (9):

9

where E_l stands for the elevation angle.

Ionosphere models were derived with a temporal sampling of one hour with the Bernese default extension ION and later combined into a 24-hour file for each day, namely the RIM IONO_BH and RIM IONO_WB, as described in Section 2.1.

When applying the ionosphere model presented in Equation (5), the ionospheric range corrections (in meters) for the zero-difference GPS observations of the i-th frequency can be computed as shown in Equation (10):

10

where the negative sign is used for phase observations and the positive sign for code observations.

2.3 Regional Ionosphere Modeling with Machine Learning

Machine learning represents a branch of artificial intelligence (AI) that can approximate the function between inputs and outputs based on rules/relationships that an AI system has “learned” from the data during the learning/training phase. This characteristic distinguishes machine learning from traditional modeling/programming approaches which require an extensive list of rules describing relationships between inputs and outputs to be explicitly specified (Natras & Schmidt, 2021). More specifically, machine learning addresses the problem of finding an approximation function that maps inputs (called predictors, features, or independent variables) to one or more outputs (called responses or dependent variables). This study focuses on the estimation of a single output (VTEC). Therefore, the output is presented below as a single-column vector. To start, a dataset was prepared that contained the measurements of the input feature and the output y for a set of observations i(i = 1, 2, …, n) as indicated in Equation (11):

11

Artificial neural networks (ANNs) represent a state-of-the-art technique with widespread applications in many fields. The ANN usually consists of an input layer with input variables, an output layer with output variables, and one or more hidden layers between them. Hidden layers allow the model to create a complex nonlinear mapping function between the network inputs and outputs by introducing a nonlinear activation function (Sharma et al., 2020). To estimate the values of the neurons in a current layer, parameters (weights) need to be added to the values of the neurons from the previous layers. The weight matrices can be defined as shown in Equation (12):

12

Where the matrix W contains all weights in the neural network. W^(k–1) are weight matrices assigned to p neurons from the previous layer (k–1) in order to estimate values of L neuron(s) in the current hidden or output layer (k) with k = (2, 3, …, K), where k = 1 corresponds to the input layer. w_{l, o} is the bias term (l = 1, 2…, L). In the matrix X we added an element x_i,o = 1 to each row. It is convenient to include this element, because of the bias term w_o as this will permit us to present Equation (13) in vector form as a scalar product. Thus, the linear model will predict the output y_i given a vector of input features x_i and weights w_l with l = 1, i.e., the model consists of an input and an output layer, which can be expressed as shown in Equation (13):

13

where e_i is an error. The weights in W are determined during the training phase using the method of least squares to minimize the sum of the squares of the error e_i and as shown in Equation (14):

14

In a fully connected network, all neurons in one layer are connected to all neurons in the next layer; this is also known as a multilayer perceptron network (Ramchoun et al., 2016). A feed-forward network indicates that the connections between the neurons are all in one direction (from input to output) and hence the information is fed-forward from one layer to the next. Afterward, the backpropagation is used to adjust the weights of the ANN using a stochastic gradient descent method (Bottou, 1991) to minimize R(W), as shown in Equation (14). In this study, we tested different setups of an ANN architecture and inputs. These specifications are presented in Table 3. In total, six ANN and two Random Forest (RF) models have been developed. One ANN model does not contain a hidden layer (ANN1); all other ANN models have hidden layers. The ANN1 model is an example of the linear model, where the outputs are estimated as described by Equation (13). The architecture of the ANN5 model is shown in Figure 3. ANN5 is a fully-connected feed-forward network with an input layer, three hidden layers, and an output layer. The number of neurons in the first layer corresponds to the length of the vector x_j. Neurons in each layer can be referred to as activation neuron vector . The input layer can be expressed as . Hidden layers compute the derived features from to . The last layer provides the output . Activation neurons in the hidden layer can be calculated as described in Equation (15):

15

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Representation of the architecture of the ANN5 model with one input, three hidden, and one output layer (right). An additional bias unit was added to the input and the hidden layers. ANN architectures were drawn using the web-based tool NN-SVG (LeNail, 2019).

where q is the activation function (Hastie et al., 2009). The activation function used in hidden layers is known as a ReLu (Rectified Linear unit) function (Sharma et al., 2020). ReLu is a non-linear activation function that is widely used in deep learning models because of its simplicity and effectiveness. It is defined in Equation (16):

16

meaning that the function will generate a value of zero for any negative value of z; for any positive value, it will return that value. The network output can be defined as shown in Equation (17):

17

The optimal hyper-parameters (number of neurons, number of hidden layers, learning rate for the stochastic gradient descent, and so on) were selected by testing different values within a certain range and observing the changes in R(W). Hyperparameters were searched within the following intervals: number of hidden units = [5, 10, 20, 30, 40, 60], number of hidden layers: [1, 2, 3, 4], epochs = [100, 200, 300, 400, 500], batch size = [50, 100, 200, 300, 400, 500]. Batch size refers to the number of training examples propagated through the network in one iteration (forward/backward pass). The number of epochs represents the number of complete passes through all the training examples. Hyperparameters that minimize R(W) were selected. The optimal number of hidden units was identified as 10 with three hidden layers. This architecture was kept constant in all ANN models from Table 3 tested; additional fine-tuning resulted in no significant differences. The addition of more neurons and more layers to the ANN increases the complexity of the interactions between layers and neurons. While this can be perceived as a positive, too much complexity can lead to overfitting the training data as the machine learns from its noise and is thus not capable of generalizing information from new examples. Choosing the appropriate complexity of an ANN requires careful attention when optimizing the model. Additionally, the Random Forest (RF) algorithm (Breiman, 2001), which builds an ensemble of decision trees, can be used to compare results and determine which input features are relevant. The RF algorithm provides the possibility of easily estimating the relative importance or contribution of each input feature (Breiman, 2001). A detailed explanation of how a decision tree and RF are built to address the VTEC problem can be found in Natras et al. (2022a). The RF model includes 300 decision trees; the quality of splits within the trees is measured with the mean squared error.

Most of the input data for these models consists of the coefficients c_nm, i.e. regional ionosphere parameters from Equation (5), which is estimated for each hour in each country within the study region (Section 2.1) by processing the GNSS observations in the Bernese GNSS Software. During the training phase, the geographical coordinates belonging to the origin of the coefficient expansion were provided with the output data of the VTEC at these positions. Information on latitude and longitude was then added to extract spatial interactions. Additional input data were introduced to each new model. The information regarding time was added to extract temporal dependencies and features. For the ANN models, the sine and cosine components are calculated to preserve their cyclic significance, as shown in Equation (18):

18

where HoD is the hour of the day. Decision tree-based algorithms, such as RF, do not require time information to be split into sine and cosine components (Boussard et al., 2017). Also, RF considers just one input feature at a time as the splitting variable when building a tree (Breiman, 2001). Therefore, it will fail to process sine and cosine components simultaneously, although they are expected to be considered as one system. Data on solar and geomagnetic activity such as the solar flux F10.7 and the geomagnetic indices of Kp and Dst were also introduced gradually. Eventually, the models (ANN6 and RF2) were trained solely on information regarding time, location, and solar and geomagnetic activity. The data were shuffled and randomly divided into training (80%) and validation datasets (20%). Using the training data, we built models that can approximate the function between inputs and outputs and estimate the outcome for new inputs. Using validation data, values for the hyperparameters were selected. All input features were scaled (i.e. standardized) to obtain data with a mean value of zero and a standard deviation of one. This is standard procedure in many learning algorithms that are sensitive to the scale of the input features (Zheng & Casari, 2018). This ensures that all inputs are treated equally, even if the variables have different scales. Also, gradient descent can converge faster, meaning that the optimal parameters for each neuron can be located more quickly. Feature standardization is defined in Equation (19):

19

where x_i,j is the data point i from j-th input feature, is the mean of the j-th input feature (over all data points from test and validation datasets together), σ_j is the standard deviation of the j-th input feature, is the resulting standardized data point i from j-th input feature. Scaling is performed on each input feature independently. Mean and standard deviation values are then stored for later data to perform the consistent transformation.

A correlation heatmap is created from a 2D correlation matrix of the input features and the output (Figure 4). Coefficients c_nm of the regional ionosphere are presented as Cnm in the heatmap. Five coefficients were estimated for each country, as explained in Section 2.1. Mean correlations for coefficients are presented. The coefficient c₀₀ has the strongest positive correlation to VTEC, followed by c₁₀ and c₀₂, both of which have a moderate negative relationship to VTEC. The F10.7 index has a positive moderate relationship, while the Dst and Kp indices have a weaker relationship with VTEC. The lowest correlations involve HoD, and latitude and longitude, which were added to extract temporal and spatial features.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Correlation heatmap of the input features and the output of the machine learning model.

Figure 5 and Table 4 show that the temporal, solar, and geomagnetic activity information increase the ANN model accuracy, respectively, in terms of the RMSE and the correlation coefficients (CCs) between the VTEC model output and VTEC from Equation (5). These input data are useful for deriving temporal, solar, and geomagnetic features and relationships. Moreover, removing the regional ionosphere coefficients from the AI models (ANN6 and RF2) increases the RMSE approximately three times and decreases the CCs by approximately 0.1. The ANN6 and RF2 models are not as accurate as the linear model with the regional ionosphere coefficients (ANN1). These results indicate that estimating the VTEC using only the spatiotemporal (latitude, longitude, and HoD), geomagnetic, and solar parameters (F10.7, Kp, and Dst) will not be sufficient. These data alone cannot provide an accurate description of VTEC. However, in conjunction with more influential, descriptive parameters of VTEC variations, such as regional ionosphere coefficients, one can significantly improve the model. Based on these results, the ANN5 model was selected as optimal and henceforth referred to as the RIM IONOWB_AI model.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

The RMSE (left) and correlation coefficients (CCs, right) for training and validation of ANN and RF models.

Knowing the underlying relationships between the inputs and the output of the AI model is useful for interpreting what the AI model has learned and which input features have been selected as relevant. The relative importance of the input features was calculated as a root square of the sum of the squared improvements over the nodes in a tree generated by RF when a specific input feature x_j was selected (Breiman, 2001). The importance of the input features was estimated using the RF1 and RF2 models (Figure 6). The relative importance of each coefficient is calculated as the mean importance of the regional ionosphere coefficients estimated for five countries. For the RF1 model, the first coefficient c₀₀ contributes most to the VTEC output, followed by HoD and F10.7 index. Contributions from other coefficients were at similar levels. Smaller contributions were provided by the Dst index, longitude, latitude, and Kp index. These parameters extract spatial dependence and provide additional geomagnetic-dependent VTEC variations and relationships, but do not have a dominant influence in defining the VTEC variations with respect to the other inputs. Because the coefficients were calculated from the CORS and EPN observations, they already contain regional ionosphere information related to the origin of the Taylor series expansion and therefore, they have a large impact on the model output. When training a model with no regional ionosphere parameters, the highest contribution (nearly 70%) comes from HoD. The other most influential inputs are the F10.7, Dst, and Kp indices. The contributions of these inputs are higher in the RF2 model than in the RF1 model. This is because the RF2 model does not include input parameters that describe the regional ionosphere. Regional ionosphere coefficients correlate much more closely with the VTEC than any of the other input features. The least critical contributions come from latitude and longitude. This may be because the F10.7, Kp, and Dst indices are not spatially dependent. These results provide us with some intuition of how the AI models learn from these datasets.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Relative importance of input variables when training AI models with (left) and without coefficients (right) estimated using Random Forest (RF).

Regional ionosphere parameters estimated from the CORS and EPN observations for each country included in the development of RIM IONO_WB were combined with spatiotemporal (latitude, longitude, and HoD), solar (F10.7 index), and geomagnetic (Kp and Dst indices) parameters via an AI technique using the ANN algorithm; this has resulted in the RIM IONOWB_AI model (Figure 7). This model can estimate VTEC for any location in the western Balkans.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

Flowchart leading to the development of RIMs. First, GPS data from local CORS and EPN observations were processed in the Bernese GNSS software to form geometry-free linear combinations of phase observations which were then used to estimate coefficients representing the regional ionosphere parameters that form the basis of the RIM IONOBH and the RIM IONOWB. Regional ionosphere parameters were then fed into the ANN along with the spatial and temporal parameters as well as solar and geomagnetic indices, resulting in the RIM IONOWB_AI model.

2.4 Selection of the Study Period

Three time periods were selected for the study:

• March 20–26, 2014,

• March 15–20, 2015, and

• March 20–26, 2018.

These periods were chosen because they included different phases of the solar cycle and varying ionosphere activity. Specifically, the solar maximum occurred in 2014, the declining phase of the solar cycle was in 2018, and a severe geomagnetic storm occurred in March 2015. There were large ionosphere disturbances over the study region that resulted in VTEC changes of 60% to 150% on the day of the storm in March 2015 and between 50% and 80% on the following days compared to its regular variability (Natras et al., 2023). Furthermore, all periods evaluated were around the spring equinox, when one can expect to find the highest number of electrons within the ionosphere over the study region (Natras et al., 2023).

The RIM IONOBH model was estimated for the first and third study periods, while the RIM IONOWB was generated for the first and second study periods. Therefore, while both models were established for the period of the solar maximum in 2014, different time frames were chosen for the second study period, including the solar minimum period (March 2018) for the RIM IONOBH and the geomagnetic storm period (March 2015) for the RIM IONOWB. This was because the models were developed independently of each other as part of two different projects implemented at different research institutes with different objectives, i.e., to study RIMs in the period of solar maximum (2014) and minimum (2018), and to analyze RIMs during a severe geomagnetic storm (2015). Furthermore, data were not available from all CORS networks for all periods examined; this resulted in different second study periods.

2.4.1 Overview of Solar and Geomagnetic Activity

Solar and geomagnetic indices are presented in Figure 8 for periods March 20–26, 2014 (left), March 15–20, 2015 (middle), and March 20-26, 2018 (right). These indices include the sunspot number R, radio flux F10.7 of the Sun’s emission at the 10.7 cm wavelength (Covington, 1969), disturbance storm time Dst (Sugiura, 1964), and Kp (Chapman & Bartels, 1962) obtained from the NASA/GSFC OMNI data set via OMNIWeb (https://omniweb.gsfc.nasa.gov/form/dxl.html). The first study period represents a solar maximum that reached its peak in April 2014. Therefore, the number of sunspots was increased throughout, ranging from 112 to 151, and F10.7 had values larger than 150 sfu. By contrast, the geomagnetic conditions were mostly quiet with Kp values below 4 with small fluctuations in the Dst index at approximately zero. The second investigation period included the strongest geomagnetic storm of the solar cycle 24, known as the St. Patrick’s Day Storm, that occurred on March 17, 2015. This severe storm was characterized by a main phase in which the Kp index reached a value of 8 and the Dst index was less than −200 nT. This was followed by a recovery phase that began on March 18, 2015, and lasted for several days until the geomagnetic field returned to normal conditions. It is worth mentioning that both the number of sunspots (from 20 to 53) and the solar radio flux (about 110 sfu) were significantly lower compared to those measured during the first period. The third study period featured low levels of solar activity during the decline phase of the solar cycle toward its minimum. This period was characterized by no sunspots and an F10.7 less than 70 sfu. The geomagnetic activity was calm to active, with no geomagnetic storms.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 8

Overview of solar and geomagnetic indices for the three study periods. Left panel: March 20-26, 2014; middle panel: March 15-20, 2015; right panel: March 20-26, 2018. From top to bottom: R sunspot number (SN), solar radio flux F10.7 in sfu (solar flux units), Dst in nT, Kp (Quiet kp · 10 < 30, Moderate 30 ≤ kp · 10 < 40, Active 40 ≤ kp · 10<50, Storm kp · 10 ≥ 50).

2.5 Validation of the Regional Ionosphere Models

The RIMs IONOBH and IONOWB were validated against other ionosphere models, including:

• GIMs final products from CODE and IGS (from https://cddis.nasa.gov/archive/gnss/products/ionex/)

• European RIM from the GIOMO model, based on weighting functions and developed at the Department of Geodesy and Geoinformation, Vienna University of Technology (Magnet, 2019)

• European RIM based on polynomial B-spline functions developed at DGFI-TUM as a two-step VTEC model (TSM-Product 2c) named as OTHR model (Goss et al., 2020).

• NeQuick2 (computed via https://t-ict4d.ictp.it/nequick2/nequick-2-web-model)

• GPS broadcasted model Klobuchar. VTEC values were calculated from broadcast ephemeris data (https://cddis.nasa.gov/archive/gnss/data/daily/).

The climatological empirical model NeQuick2 (Nava et al., 2008) is the latest version of the NeQuick ionosphere electron density model developed for the computation of slant electron density profiles and TEC in the ionosphere at a given height, geocentric latitude, and geocentric longitude. Information on solar activity is provided to the model by the daily solar radio flux F10.7. The NeQuick2 model takes into account the contribution of the plasmasphere to VTEC up to 40,000 km. Although the electron density in the plasmasphere is much lower than in the ionosphere, under some conditions, such as at night and during periods of low solar activity, the contribution of the plasmasphere may represent a larger proportion of the VTEC. The relative global contribution of the plasmasphere to the VTEC depends on latitude and solar activity, with a minimum contribution of about 10% during daytime hours and a maximum of up to 60% at night.

The VTEC data of the RIM GIOMO model were provided for the study periods in March 2014 and March 2018 during the project that led to the development of the IONOBH. By contrast, VTEC data from the RIM OTHR model were available for all study periods.

The ionosphere models were also evaluated for precise point positioning (PPP) by processing single-frequency (L₁) observations for selected GNSS stations using the Bernese GNSS Software. The following cases of L₁ PPP solutions were carried out:

• L₁ positioning solution without ionospheric corrections,

• L₁ positioning solution with ionospheric corrections from the final GIM CODE, and

• L₁ positioning solution with ionospheric corrections from the RIMs IONOBH and IONOWB.

However, there is no interface in the Bernese software that can be used to process broadcasted ionospheric delay models. Therefore, we could not apply corrections from the Klobuchar model. The following processing steps were implemented for the PPP method: data preprocessing (RNXSMT), import of data into the Bernese format (RXOBV3), preparing orbit and Earth’s orientation information (POLUPD, PRETAB, ORBGEN), data preprocessing 2 (CODSPP, GPSEST, RESRMS, SATMRK), performing a solution for epoch parameters, and/or creation of normal equations NEQ (GPSEST) and an NEQ-based final session solution (ADDNEQ2). More details on each step are included in Dach et al. (2015). Positioning errors were estimated as differences between “true” positions and the single-frequency positioning results expressed as north, east, and up components. The weekly combined EPN solutions (https://www.epncb.oma.be/ftp/product/combin/) were used as “true” positions for the EPN stations. For the CORS stations, “true” positions were estimated with the dual-frequency PPP method in the Bernese GNSS software. For the analysis, errors were expressed as 1-D RMS vertical and 2-D RMS horizontal position errors, as well as, 3-D RMS position errors as shown in Equation (20):

20a

20b

20c

where ΔE_i, ΔN_i, ΔU_i are errors in the east, north, and up (vertical) components, respectively, of the i-th position estimate sample and n is the total number of position estimate samples.