# Next Generation Network Real-Time Kinematic Interpolation Segment to Improve the User Accuracy

Authors: Mohammed Ouassou; Anna B. O. Jensen; Jon G. O. Gjevestad; Oddgeir Kristiansen
Journal: International Journal of Navigation and Observation
 (2015)

---

## Abstract

This paper demonstrates that automatic selection of the right interpolation/smoothing method in a GNSS-based network real-time kinematic (NRTK) interpolation segment can improve the accuracy of the rover position estimates and also the processing time in the NRTK processing center. The methods discussed and investigated are inverse distance weighting (IDW); bilinear and bicubic spline interpolation; kriging interpolation; thin-plate splines; and numerical approximation methods for spatial processes. The methods are implemented and tested using GNSS data from reference stations in the Norwegian network RTK service called CPOS. Data sets with an average baseline between reference stations of 60–70 km were selected. 12 prediction locations were used to analyze the performance of the interpolation methods by computing and comparing different measures of the goodness of fit such as the root mean square error (RMSE), mean square error, and mean absolute error, and also the computation time was compared. Results of the tests show that ordinary kriging with the Matérn covariance function clearly provides the best results. The thin-plate spline provides the second best results of the methods selected and with the test data used.


---

## Body

## 1. Introduction

The use of GNSS and network real-time kinematic positioning to achieve GNSS positions with accuracy at the cm-level is increasing rapidly these years. This is partly due to the development and modernization of the GNSS systems themselves (GPS, GLONASS, Galileo, and Beidou), but it is also caused by a general quest for better position accuracy in many user communities.High-accuracy GNSS positioning is based on the carrier phase being observable. Using the notation from [1], the basic observation equation that summarizes the relation between observations and error sources is given as follows:
(1)
Φ
=
λ
-
1
(
r
-
I
+
T
)
+
f
(
δ
t
u
-
δ
t
s
)
+
N
+
ε
ϕ
,
where Φ is the phase observation in cycles, λ is the wavelength in meters/cycle, r is the geometric distance between the receiver and satellite in meters, I is the ionospheric signal delay in meters, T is the tropospheric signal delay in meters, f is the frequency in Hertz, δ
t
u and δ
t
s are the clock errors of, respectively, the receiver and the satellite, N is the initial number of cycles at the first observation epoch (theambiguity), and ε
ϕ is a noise term given in cycles that mainly accounts for multipath (reflected signals) and receiver noise.When using the NRTK technique, a network of reference stations is used to estimate the errors in the positioning process, that is, the effects of the ionosphere and troposphere as well as inaccuracies in the satellite position as provided with the broadcast ephemerids from the satellites.The accuracy of NRTK positioning systems depends on the ability to identify and mitigate the error sources in the system as well as the residual biases. The biases include residual effects from the space segment, signal propagation, environment effects, and receiver noise in the reference network. The mitigation process can be carried out by modeling, estimation, or combinations of observables.The NRTK processing chain can be summarized as follows: the first step is to collect raw measurements from the network of reference stations, solve for the ambiguities within the reference network, and generate error estimates. The next step is to apply the interpolation/smoothing scheme to generate the RTK corrections for the user location. The RTK corrections are then transmitted to users who can then perform real-time positioning with accuracy at the cm-level.Figure1 shows all segments involved in the NRTK processing chain. The figure illustrates the so-called virtual reference station (VRS) concept, which was developed by Landau et al. [2]. Other NRTK standards such as for instance the master auxiliary concept (MAC) also exist [3], but we limit the discussion in this paper to the VRS concept.Figure 1
Network real-time kinematic segments.As the GNSS systems and users become more numerous, the amount of data that needs processing increases as well, which poses some interesting challenges for the NRTK system developers and service providers. This paper focuses on processing large data sets and high quality interpolators/smoothers that can be used to aid the data processing. Let us consider how the RTK processing is carried out. First the user sends his/her position to the control center, and then the network engine chooses a suitable subnetwork which is used to generate corrections, and these corrections are then transmitted back to the user. The first challenge to this model is the number of users, since each user has to be processed independently, and the number of users has increased dramatically in recent years. The solution to this is to construct new models and algorithms. These should be able to process data from large geographical areas, as well as computing the necessary corrections and quality indicators ready for use, so that any RTK user that connects will be served immediately.In other branches of science and engineering, new analysis tools that satisfy these requirements have already been developed: neural networks, machine learning, classification and regression trees, hierarchical models, and so forth. In this paper, some existing interpolation/smoothing methods are applied to real datasets, and the strengths and weaknesses of each method are identified. The results are then used to combine these methods and construct models that describe the observed variations in the data as well as possible.Interpolation methods can be divided into two categories: local methods and global methods. The local methods only use a subset of the data for interpolation, which implies that the required processing time is reduced. Conversely, the global techniques use all the data available to generate predictions. In this paper, both these approaches are considered. Referring to Figure1, the main focus of this paper is directed at the correction interpolation segment and more specifically at the automatic selection of the right interpolation algorithm based on appropriate tests, such that the rover position estimation will be improved.The rest of the paper is organized as follows: Section2 gives a full description of the test data using the Norwegian GNSS network data, known as CPOS, and introduces the variational problem in general. Section 3 covers local interpolation algorithms, specifically the inverse-distance weighted and bilinear/bicubic methods by Akima. Section 4 and the following sections deal with global interpolation methods. First, thin-plate splines and the Bayesian model behind the smoothing are reviewed in this section. Section 5 introduces numerical approximation schemes for Gaussian random fields. Section 6 covers spatial interpolation algorithms and specifically the ordinary kriging method. Section 7: the performance parameters are defined in this section. Section 8: the results from Sections 3–6 are generalized in this section. Section 9 is the conclusion and discussion and covers applications of the results developed in Sections 3–6.


## 2. Test Data and Variational Reconstruction

### 2.1. Test Data

The main success of network real-time kinematic positioning has been the reduction of correlated errors in the network (e.g., ionospheric, tropospheric, and satellite position errors). This type of errors is collectively referred to asdistance-dependent errors and can be subdivided into thedispersive errors which depend on frequency and thenondispersive errors which do not.The size of the network varies with time, as the individual reference stations and satellites may not deliver data for a while, and observations are typically correlated to each other. Modeling the spatial and temporal variations of such a process is too complex to capture the covariance structure of the data, so often we end up imposing stationarity. In this paper, we apply techniques for handling spatial processes in order to capture the covariance structure in the data, such that high quality synthetic data can be provided. The main clue is to employ the right tool from epoch to epoch, based on some appropriate criteria.We prefer to work with real data, and since the real network error estimates were not made available, we decided to analyze the ionospheric path delays for CPOS RTK network, given by absolute TEC values. If the ionosphere data is replaced with the full network corrections, the same algorithms should still function very well. Ionospheric path delay is considered the single largest source of inaccuracy for positioning and navigation, so the quality of the NRTK corrections is strongly affected in the case of moderate to high ionosphere activity. To test the algorithms against each other, a large ionospheric data set from the Norwegian CPOS network is investigated. The data is generated by a first-order geometry-free approach (Section8.1).At the time of writing, the CPOS RTK network contains approximately 180 stations on the Norwegian mainland (for a listing of stations of the CPOS RTK network, see the appendix.). The algorithms were tested with different station configurations (50, 75, and 110 stations), equipped with geodetic dual frequency GNSS receivers, which track both GPS and GLONASS satellites. In this investigation, however, only the GPS portion of the data was used. The distribution of the CPOS RTK network reference stations is given in Figure2.Figure 2
Norwegian CPOS RTK network. Filled triangles mark station locations. Height difference between sites up to 1000 meters.


### 2.2.  Variational Reconstruction

#### 2.2.1. Problem Formulation

Let us assume that the observed noisy measurements{
z
i
} at the locations {
s
i
:
s
i
∈
R
2
} is a random function, with the mean m
(
s
) and variance σ
2. Our goal is then to predict the value at other locations {
s
k
} where we have no observations, under the assumption that the predicted values should resemble its neighbors. To achieve this, we can either interpolate or construct a smooth function g
(
s
) that represents the variation in the data and is robust against outliers.The data that will be modeled is a pure spatiotemporal process, namely, the absolute total electron count (TEC). Assuming weak stationarity of the process under study, the mean and variance are not functions of the spatial locations. The model used to describe the variation in data in this paper, however, is assumed to have the form
(2)
z
s
=
m
s
+
ε
s
.
The mean function m
(
s
), often referred to as thetrend ordeterministic part, determines the large-scale variation in the data. The function ε
(
s
) is called therandom part and determines the small-scale variation. This process model will be assumed in the subsequent discussion of all the different interpolation/smoothing techniques presented in this paper.Some data analysts prefer theCressie decomposition [4, ch. 3] of the observed random field Z
(
s
), which takes the form
(3)
Z
s
=
m
s
+
T
s
+
ε
s
,
where s
∈
R
2 is the spatial location; Z
(
s
) is the observation; m
(
s
)
=
X
β is the trend (the mean component of the model); T
(
s
) is a stationary Gaussian process with variance σ
2 (partial sill), and a correlation function parameterized in its simplest form by ϕ (the range parameter); and finally ε is an error term, with a variance parameter τ
2 (nugget variance).


#### 2.2.2.  Model Parameter Estimations

Once the model is defined, the next step is to estimate the model parameters. In general, this is done numerically by minimizing the negative log-likelihood function. The most used optimization methods are, respectively, the conjugate gradient method, the quasi-Newtonian method, and the Nedler-Mead method. The details of these methods will not be treated in this paper, but the interested reader is referred to references [5, 6].The algorithm may not converge to correct parameter values when called with the default options. The user should therefore try different initial values, and if the parameters have different orders of magnitude, a scaling of the parameters may be necessary. If such problems arise, some possible workarounds include(i)
rescaling data values by dividing by a constant,
(ii)
rescaling coordinates by subtracting values and/or dividing by constants,
(iii)
bootstrapping to accelerate the convergence. This method is used in our implementation of the kriging algorithm in Section6.


#### 2.2.3. Model Validation

In the field of statistical analysis, an appropriate way of analyzing data is to divide it into three distinct subsets. The training dataset is used to construct the model, the validation data is used to check the model, and the last data set is used to challenge the model. The main purpose is to determine whether or not our model is an accurate representation of the real world data. This process is called themodel validation assessment. The most famous methods are the family of cross-validation, generalized maximum likelihood (GML) methods, Akaike information criterion (AIC), Bayesian information criterion (BIC), and so forth.In our implementation, the generalized cross-validation is used to determine the optimal smoothing parameterλ (see Section 4). The computation AIC and BIC are computed in Section 6, when maximum likelihood estimation is used instead of weighted least squares in the kriging algorithm. The GML methods will be used in future work.


## 2.1. Test Data

The main success of network real-time kinematic positioning has been the reduction of correlated errors in the network (e.g., ionospheric, tropospheric, and satellite position errors). This type of errors is collectively referred to asdistance-dependent errors and can be subdivided into thedispersive errors which depend on frequency and thenondispersive errors which do not.The size of the network varies with time, as the individual reference stations and satellites may not deliver data for a while, and observations are typically correlated to each other. Modeling the spatial and temporal variations of such a process is too complex to capture the covariance structure of the data, so often we end up imposing stationarity. In this paper, we apply techniques for handling spatial processes in order to capture the covariance structure in the data, such that high quality synthetic data can be provided. The main clue is to employ the right tool from epoch to epoch, based on some appropriate criteria.We prefer to work with real data, and since the real network error estimates were not made available, we decided to analyze the ionospheric path delays for CPOS RTK network, given by absolute TEC values. If the ionosphere data is replaced with the full network corrections, the same algorithms should still function very well. Ionospheric path delay is considered the single largest source of inaccuracy for positioning and navigation, so the quality of the NRTK corrections is strongly affected in the case of moderate to high ionosphere activity. To test the algorithms against each other, a large ionospheric data set from the Norwegian CPOS network is investigated. The data is generated by a first-order geometry-free approach (Section8.1).At the time of writing, the CPOS RTK network contains approximately 180 stations on the Norwegian mainland (for a listing of stations of the CPOS RTK network, see the appendix.). The algorithms were tested with different station configurations (50, 75, and 110 stations), equipped with geodetic dual frequency GNSS receivers, which track both GPS and GLONASS satellites. In this investigation, however, only the GPS portion of the data was used. The distribution of the CPOS RTK network reference stations is given in Figure2.Figure 2
Norwegian CPOS RTK network. Filled triangles mark station locations. Height difference between sites up to 1000 meters.


## 2.2.  Variational Reconstruction

### 2.2.1. Problem Formulation

Let us assume that the observed noisy measurements{
z
i
} at the locations {
s
i
:
s
i
∈
R
2
} is a random function, with the mean m
(
s
) and variance σ
2. Our goal is then to predict the value at other locations {
s
k
} where we have no observations, under the assumption that the predicted values should resemble its neighbors. To achieve this, we can either interpolate or construct a smooth function g
(
s
) that represents the variation in the data and is robust against outliers.The data that will be modeled is a pure spatiotemporal process, namely, the absolute total electron count (TEC). Assuming weak stationarity of the process under study, the mean and variance are not functions of the spatial locations. The model used to describe the variation in data in this paper, however, is assumed to have the form
(2)
z
s
=
m
s
+
ε
s
.
The mean function m
(
s
), often referred to as thetrend ordeterministic part, determines the large-scale variation in the data. The function ε
(
s
) is called therandom part and determines the small-scale variation. This process model will be assumed in the subsequent discussion of all the different interpolation/smoothing techniques presented in this paper.Some data analysts prefer theCressie decomposition [4, ch. 3] of the observed random field Z
(
s
), which takes the form
(3)
Z
s
=
m
s
+
T
s
+
ε
s
,
where s
∈
R
2 is the spatial location; Z
(
s
) is the observation; m
(
s
)
=
X
β is the trend (the mean component of the model); T
(
s
) is a stationary Gaussian process with variance σ
2 (partial sill), and a correlation function parameterized in its simplest form by ϕ (the range parameter); and finally ε is an error term, with a variance parameter τ
2 (nugget variance).


### 2.2.2.  Model Parameter Estimations

Once the model is defined, the next step is to estimate the model parameters. In general, this is done numerically by minimizing the negative log-likelihood function. The most used optimization methods are, respectively, the conjugate gradient method, the quasi-Newtonian method, and the Nedler-Mead method. The details of these methods will not be treated in this paper, but the interested reader is referred to references [5, 6].The algorithm may not converge to correct parameter values when called with the default options. The user should therefore try different initial values, and if the parameters have different orders of magnitude, a scaling of the parameters may be necessary. If such problems arise, some possible workarounds include(i)
rescaling data values by dividing by a constant,
(ii)
rescaling coordinates by subtracting values and/or dividing by constants,
(iii)
bootstrapping to accelerate the convergence. This method is used in our implementation of the kriging algorithm in Section6.


### 2.2.3. Model Validation

In the field of statistical analysis, an appropriate way of analyzing data is to divide it into three distinct subsets. The training dataset is used to construct the model, the validation data is used to check the model, and the last data set is used to challenge the model. The main purpose is to determine whether or not our model is an accurate representation of the real world data. This process is called themodel validation assessment. The most famous methods are the family of cross-validation, generalized maximum likelihood (GML) methods, Akaike information criterion (AIC), Bayesian information criterion (BIC), and so forth.In our implementation, the generalized cross-validation is used to determine the optimal smoothing parameterλ (see Section 4). The computation AIC and BIC are computed in Section 6, when maximum likelihood estimation is used instead of weighted least squares in the kriging algorithm. The GML methods will be used in future work.


## 2.2.1. Problem Formulation

Let us assume that the observed noisy measurements{
z
i
} at the locations {
s
i
:
s
i
∈
R
2
} is a random function, with the mean m
(
s
) and variance σ
2. Our goal is then to predict the value at other locations {
s
k
} where we have no observations, under the assumption that the predicted values should resemble its neighbors. To achieve this, we can either interpolate or construct a smooth function g
(
s
) that represents the variation in the data and is robust against outliers.The data that will be modeled is a pure spatiotemporal process, namely, the absolute total electron count (TEC). Assuming weak stationarity of the process under study, the mean and variance are not functions of the spatial locations. The model used to describe the variation in data in this paper, however, is assumed to have the form
(2)
z
s
=
m
s
+
ε
s
.
The mean function m
(
s
), often referred to as thetrend ordeterministic part, determines the large-scale variation in the data. The function ε
(
s
) is called therandom part and determines the small-scale variation. This process model will be assumed in the subsequent discussion of all the different interpolation/smoothing techniques presented in this paper.Some data analysts prefer theCressie decomposition [4, ch. 3] of the observed random field Z
(
s
), which takes the form
(3)
Z
s
=
m
s
+
T
s
+
ε
s
,
where s
∈
R
2 is the spatial location; Z
(
s
) is the observation; m
(
s
)
=
X
β is the trend (the mean component of the model); T
(
s
) is a stationary Gaussian process with variance σ
2 (partial sill), and a correlation function parameterized in its simplest form by ϕ (the range parameter); and finally ε is an error term, with a variance parameter τ
2 (nugget variance).


## 2.2.2.  Model Parameter Estimations

Once the model is defined, the next step is to estimate the model parameters. In general, this is done numerically by minimizing the negative log-likelihood function. The most used optimization methods are, respectively, the conjugate gradient method, the quasi-Newtonian method, and the Nedler-Mead method. The details of these methods will not be treated in this paper, but the interested reader is referred to references [5, 6].The algorithm may not converge to correct parameter values when called with the default options. The user should therefore try different initial values, and if the parameters have different orders of magnitude, a scaling of the parameters may be necessary. If such problems arise, some possible workarounds include(i)
rescaling data values by dividing by a constant,
(ii)
rescaling coordinates by subtracting values and/or dividing by constants,
(iii)
bootstrapping to accelerate the convergence. This method is used in our implementation of the kriging algorithm in Section6.


## 2.2.3. Model Validation

In the field of statistical analysis, an appropriate way of analyzing data is to divide it into three distinct subsets. The training dataset is used to construct the model, the validation data is used to check the model, and the last data set is used to challenge the model. The main purpose is to determine whether or not our model is an accurate representation of the real world data. This process is called themodel validation assessment. The most famous methods are the family of cross-validation, generalized maximum likelihood (GML) methods, Akaike information criterion (AIC), Bayesian information criterion (BIC), and so forth.In our implementation, the generalized cross-validation is used to determine the optimal smoothing parameterλ (see Section 4). The computation AIC and BIC are computed in Section 6, when maximum likelihood estimation is used instead of weighted least squares in the kriging algorithm. The GML methods will be used in future work.


## 3. IDW and Akima Interpolation

### 3.1. IDW Interpolation

The inverse-distance weighted (IDW) scheme is an exact interpolator. It honors the data by assigning weights to all neighboring points according to their distance from the prediction locations
0. Locations that are closer to s
0 receive higher weights, and locations that are far from s
0 are given lower weights; this mechanism is administered by the parameter r in the IDW predictor formula. The user can freely choose the number of observations used to perform the interpolation. This is done by defining a radius R around the prediction location s
0.The IDW predictor scheme is defined as(4)
z
^
(
s
0
)
=
∑
i
=
1
N
z
(
s
i
)
d
0
i
-
r
∑
i
=
1
N
d
0
i
-
r
.
Here, s
0 is the prediction location, N is the number of observations, z
(
s
i
) are the neighboring observations, r is the weight decision parameter, and d is the distance (either spherical or Euclidean).The IDW method is originally due to Shepard [7], which described a global method. All derived IDW methods are either generalizations or variations of this method. The basic Shepard’s method can be expressed as
(5)
F
(
x
,
y
)
=
∑
k
=
1
N
w
k
(
x
,
y
)
f
k
∑
k
=
1
N
w
k
(
x
,
y
)
,
where typically the weight w
k is the inverse Euclidean distance w
k
=
1
/
d
k
2
=
(
x
-
x
k
)
2
+
(
y
-
y
k
)
2
-
1
/
2. We will however define a disk with center (
x
k
,
y
k
) and a radius R and set the weight w
k to zero outside of this disk. A natural scheme suggested by many authors, like, for example, Renka and Brown [8], is given by the expression
(6)
w
k
x
,
y
=
(
R
w
-
d
k
)
+
R
w
d
k
2
,
where
(7)
R
w
-
d
k
+
=
R
w
-
d
k
if
d
k
<
R
w
0
if
d
k
≥
R
w
.Impose the constraints such that(i)
the sum of all weightsw
k inside the disk R
w should be normalized to unity, that is, ∑
‍
w
k
=
1,
(ii)
the predictor is a linear combination of the observations. If the variance of the predictor is then controlled such that it is at a minimum, the IDW behaves almost like the local kriging interpolator (Section 6); however the covariance structure is not preserved.For the implementation, the packagegstat from Edzer Pebesma is used to carry out IDW (see Table 3 for more information).


### 3.2.  Akima Algorithms

Bilinear or bicubic spline interpolation is applied using different versions of algorithms by Akima [9, 10]. Given a set of data points in a plane, our aim is to fit a smooth curve that passes through the given points considered as reference points. The method is local and is based on a piecewise function composed of a set of polynomials and applicable up to the third degree on each interval. The method produces remarkable results with minimum processing time. For a detailed mathematical formulation, please refer to references [9, 10].


#### 3.2.1. Basics of Spline Interpolation/Smoothing

An understanding of the basic elementary building blocks of a 1D spline facilitates the understanding of 2D and 3D splines, for instance, the TPS (Section4).Given real numbers{
t
1
,
t
2
,
…
,
t
n
}
∈
[
a
,
b
], a function g defined on interval [
a
,
b
] is a cubic spline if and only if the following conditions are satisfied:(i)
the functiong is a cubic polynomial on each interval (
a
,
t
1
)
,
(
t
1
,
t
2
)
,
…
,
(
t
n
,
b
);
(ii)
the functiong and its first and second derivatives are continuous at each of the points t
i. Condition (ii) implies that the cubic polynomials from condition (i) fit together on each t
i, where the t
i are calledknots. Together these two conditions imply that g
(
x
) is a function with continuous first and second derivatives on the whole interval [
a
,
b
].For some given real constantsa
i
,
b
i
,
c
i
,
d
i, the cubic spline function g can be expressed as
(8)
g
(
t
)
=
a
i
+
b
i
(
t
-
t
i
)
+
c
i
t
-
t
i
2
+
d
i
t
-
t
i
3
,
where the index i
=
0,1
,
2
,
…
,
n. The end-point knots correspond to the boundaries of the function domain; that is, t
0
=
a and t
n
+
1
=
b.Finding a smoothing spline is not an easy task. Reinsch (1967) proposed an algorithm and showed that the solution of the minimum principle is actually cubic splines. The basic idea is to construct a nonsingular system of linear equations of the second derivativeλ
i of g
^. The resulting equations are computationally efficient because of their banded structure. For an excellent exposition of the material, see also [11].


#### 3.2.2. Output Result from Akima

Figure3 shows the output from Akima with bilinear interpolation.Figure 3
L1-VTEC variation surface generated by AKIMA (bilinear) interpolation. Color scale is expressed in meters. High ionospheric activity (year: 2013, doy: 152), 110 reference stations.


## 3.1. IDW Interpolation

The inverse-distance weighted (IDW) scheme is an exact interpolator. It honors the data by assigning weights to all neighboring points according to their distance from the prediction locations
0. Locations that are closer to s
0 receive higher weights, and locations that are far from s
0 are given lower weights; this mechanism is administered by the parameter r in the IDW predictor formula. The user can freely choose the number of observations used to perform the interpolation. This is done by defining a radius R around the prediction location s
0.The IDW predictor scheme is defined as(4)
z
^
(
s
0
)
=
∑
i
=
1
N
z
(
s
i
)
d
0
i
-
r
∑
i
=
1
N
d
0
i
-
r
.
Here, s
0 is the prediction location, N is the number of observations, z
(
s
i
) are the neighboring observations, r is the weight decision parameter, and d is the distance (either spherical or Euclidean).The IDW method is originally due to Shepard [7], which described a global method. All derived IDW methods are either generalizations or variations of this method. The basic Shepard’s method can be expressed as
(5)
F
(
x
,
y
)
=
∑
k
=
1
N
w
k
(
x
,
y
)
f
k
∑
k
=
1
N
w
k
(
x
,
y
)
,
where typically the weight w
k is the inverse Euclidean distance w
k
=
1
/
d
k
2
=
(
x
-
x
k
)
2
+
(
y
-
y
k
)
2
-
1
/
2. We will however define a disk with center (
x
k
,
y
k
) and a radius R and set the weight w
k to zero outside of this disk. A natural scheme suggested by many authors, like, for example, Renka and Brown [8], is given by the expression
(6)
w
k
x
,
y
=
(
R
w
-
d
k
)
+
R
w
d
k
2
,
where
(7)
R
w
-
d
k
+
=
R
w
-
d
k
if
d
k
<
R
w
0
if
d
k
≥
R
w
.Impose the constraints such that(i)
the sum of all weightsw
k inside the disk R
w should be normalized to unity, that is, ∑
‍
w
k
=
1,
(ii)
the predictor is a linear combination of the observations. If the variance of the predictor is then controlled such that it is at a minimum, the IDW behaves almost like the local kriging interpolator (Section 6); however the covariance structure is not preserved.For the implementation, the packagegstat from Edzer Pebesma is used to carry out IDW (see Table 3 for more information).


## 3.2.  Akima Algorithms

Bilinear or bicubic spline interpolation is applied using different versions of algorithms by Akima [9, 10]. Given a set of data points in a plane, our aim is to fit a smooth curve that passes through the given points considered as reference points. The method is local and is based on a piecewise function composed of a set of polynomials and applicable up to the third degree on each interval. The method produces remarkable results with minimum processing time. For a detailed mathematical formulation, please refer to references [9, 10].


### 3.2.1. Basics of Spline Interpolation/Smoothing

An understanding of the basic elementary building blocks of a 1D spline facilitates the understanding of 2D and 3D splines, for instance, the TPS (Section4).Given real numbers{
t
1
,
t
2
,
…
,
t
n
}
∈
[
a
,
b
], a function g defined on interval [
a
,
b
] is a cubic spline if and only if the following conditions are satisfied:(i)
the functiong is a cubic polynomial on each interval (
a
,
t
1
)
,
(
t
1
,
t
2
)
,
…
,
(
t
n
,
b
);
(ii)
the functiong and its first and second derivatives are continuous at each of the points t
i. Condition (ii) implies that the cubic polynomials from condition (i) fit together on each t
i, where the t
i are calledknots. Together these two conditions imply that g
(
x
) is a function with continuous first and second derivatives on the whole interval [
a
,
b
].For some given real constantsa
i
,
b
i
,
c
i
,
d
i, the cubic spline function g can be expressed as
(8)
g
(
t
)
=
a
i
+
b
i
(
t
-
t
i
)
+
c
i
t
-
t
i
2
+
d
i
t
-
t
i
3
,
where the index i
=
0,1
,
2
,
…
,
n. The end-point knots correspond to the boundaries of the function domain; that is, t
0
=
a and t
n
+
1
=
b.Finding a smoothing spline is not an easy task. Reinsch (1967) proposed an algorithm and showed that the solution of the minimum principle is actually cubic splines. The basic idea is to construct a nonsingular system of linear equations of the second derivativeλ
i of g
^. The resulting equations are computationally efficient because of their banded structure. For an excellent exposition of the material, see also [11].


### 3.2.2. Output Result from Akima

Figure3 shows the output from Akima with bilinear interpolation.Figure 3
L1-VTEC variation surface generated by AKIMA (bilinear) interpolation. Color scale is expressed in meters. High ionospheric activity (year: 2013, doy: 152), 110 reference stations.


## 3.2.1. Basics of Spline Interpolation/Smoothing

An understanding of the basic elementary building blocks of a 1D spline facilitates the understanding of 2D and 3D splines, for instance, the TPS (Section4).Given real numbers{
t
1
,
t
2
,
…
,
t
n
}
∈
[
a
,
b
], a function g defined on interval [
a
,
b
] is a cubic spline if and only if the following conditions are satisfied:(i)
the functiong is a cubic polynomial on each interval (
a
,
t
1
)
,
(
t
1
,
t
2
)
,
…
,
(
t
n
,
b
);
(ii)
the functiong and its first and second derivatives are continuous at each of the points t
i. Condition (ii) implies that the cubic polynomials from condition (i) fit together on each t
i, where the t
i are calledknots. Together these two conditions imply that g
(
x
) is a function with continuous first and second derivatives on the whole interval [
a
,
b
].For some given real constantsa
i
,
b
i
,
c
i
,
d
i, the cubic spline function g can be expressed as
(8)
g
(
t
)
=
a
i
+
b
i
(
t
-
t
i
)
+
c
i
t
-
t
i
2
+
d
i
t
-
t
i
3
,
where the index i
=
0,1
,
2
,
…
,
n. The end-point knots correspond to the boundaries of the function domain; that is, t
0
=
a and t
n
+
1
=
b.Finding a smoothing spline is not an easy task. Reinsch (1967) proposed an algorithm and showed that the solution of the minimum principle is actually cubic splines. The basic idea is to construct a nonsingular system of linear equations of the second derivativeλ
i of g
^. The resulting equations are computationally efficient because of their banded structure. For an excellent exposition of the material, see also [11].


## 3.2.2. Output Result from Akima

Figure3 shows the output from Akima with bilinear interpolation.Figure 3
L1-VTEC variation surface generated by AKIMA (bilinear) interpolation. Color scale is expressed in meters. High ionospheric activity (year: 2013, doy: 152), 110 reference stations.


## 4. Thin Plate Spline Method

### 4.1. Mathematical Preliminaries

In this section, our main interest is not to construct a functionm
(
s
i
) that exactly interpolates the data z
i at n distinct points but to find an attractive way to smooth noisy data. The method of thin-plate splines (TPS) will be used for this purpose.Duchon [12] was the first to build the theoretical foundation for the TPS method. The name TPS comes from the physical situation of bending a thin surface, where the method minimizes the bending energy of a thin plate fixed at the data sites. For our application, the TPS method is used to minimize the cost function:
(9)
J
m
(
f
)
=
∫
R
d
D
m
f
(
s
)
2
d
s
,
where D
m is a vector of partial differentiation operators of order m. In the two-dimensional case, that is, when d
=
2, m
=
2, and f
=
f
(
x
,
y
), the TPS penalty function can be written as
(10)
J
2
f
=
∬
R
2
∂
2
f
∂
x
2
2
+
2
∂
2
f
∂
x
∂
y
2
+
∂
2
f
∂
y
2
2
d
x
d
y
.LetΔ
4
(
x
,
y
) denote the differential operator in the integrand of (10). The thin-plate spline f
(
x
,
y
), which is the solution to the variational problem of minimizing the penalty J
2
(
f
), can then be found by solving the biharmonic equation Δ
4
f
(
x
,
y
)
=
0.The goal is to find the functionf
λ in Sobolev space [13, p. 250] that minimizes the following expression:
(11)
S
(
f
)
=
n
-
1
∑
i
=
1
n
‍
z
(
s
i
)
-
f
(
s
i
)
2
+
λ
J
2
(
f
)
,
where n is the total number of observations, J
2 is a smoothness penalty (the cost function), and λ is the smoothing parameter. The smoothing parameter is a trade-off factor between the rate of change of the residual error and local variation. Optimal minimization of S
(
f
) results in a good compromise between smoothness and goodness of fit.Once the curve approximation of the data has been constructed, generating values at any location, where no observations are available, is accomplished by simply indexing the variablesx and y and fetching the corresponding value. This is a major advantage of smoothing methods over interpolation methods; no extra interpolations are required after the curve has been constructed for a given epoch.Validation is carried out by the GCV (Section2.2.3). Let the cross-validation function g
^
n
,
λ
[
k
], with k
=
1,2
,
…
,
n, be defined as
(12)
g
^
n
,
λ
[
k
]
=
arg
min
⁡
g
^
∈
W
2
(
m
)
∑
i
=
1
n
‍
z
i
-
g
^
x
i
2
+
λ
∫
0
1
‍
g
^
m
x
i
2
d
x
,
where λ
>
0. The OCV (ordinary cross-validation) and OCV MSE (ordinary cross-validation mean square error) V
0
(
λ
), respectively, are defined as
(13)
OCV
λ
=
arg
min
⁡
λ
∈
R
+
V
0
λ
,
V
0
λ
=
1
n
∑
i
=
1
n
‍
g
^
n
,
λ
k
x
i
-
z
i
2
.The determination of the GCV (general cross-validation) goes as follows. First, the expression forV
0 has to be rewritten. There exists an n
×
n matrix A
(
λ
), the smoothing/influence/sensitivity matrix with the property. Consider
(14)
g
^
n
,
λ
(
x
1
)
⋮
g
^
n
,
λ
(
x
n
)
=
A
λ
Z
,
such that V
0
(
λ
) can be written as
(15)
V
0
λ
=
n
-
1
∑
j
=
1
n
‍
a
k
j
z
j
-
z
k
2
1
-
a
k
k
2
,
where k
,
j
∈
{
1,2
,
…
,
n
} and a
k
j is element {
k
,
j
} of A
(
λ
).Definition 1 (generalized cross-validation (GCV)).
LetA
(
λ
) be the smoothing matrix defined in (14); then the GCV function is given by the expression
(16)
GCV
λ
=
n
-
1
I
-
A
λ
z
2
n
-
1
tr
⁡
I
-
A
λ
2
.


### 4.2. Estimation of the Smoothing Parameterλ

The smoothing parameterλ plays a central role in the TPS method. By adjusting the value of λ, one can get the desired level of smoothness at the cost of accuracy at the data sites. When we set this parameter to zero, the problem is reduced to an interpolation with no smoothing. On the other hand, when the smoothing parameter tends to infinity, the method yields a plane which is least-square fitted to the data. The smoothness penalty method can be chosen by any criteria, but the most popular criterion is GCV (generalized cross-validation), also known as the “left-out one” method. The GCV criterion selects the smoothing parameter λ that minimizes the GCV function, equation (16), that is λ
^
=
arg
min
⁡
λ
∈
R
+
⁡
GCV
(
λ
).The GCV functionV
(
λ
) is the predicted mean square error and can be viewed as a weighted version of the OCV
(
λ
)
=
V
0
(
λ
):
(17)
V
λ
=
1
n
∑
i
=
1
n
‍
g
^
n
,
λ
i
x
i
-
z
i
2
W
k
λ
,
W
k
λ
=
I
-
a
k
k
λ
n
-
1
tr
⁡
I
-
A
λ
2
.In geodesy, it is often interesting to estimate the accuracyf
λ. Two loss functions are considered: the mean square prediction error T
(
λ
), and the stricter Sobolev error is defined as W
(
λ
)
=
f
λ
-
f
W
2,
(18)
W
λ
=
∫
a
b
‍
f
λ
-
f
2
d
G
+
∫
a
b
‍
f
λ
m
-
f
m
2
d
G
.
The performance of an estimator is often well characterized by the risk function, defined as the expectation value of the loss function:
(19)
R
T
λ
=
E
T
λ
,
R
W
λ
=
E
W
λ
.In this analysis, the GCV is used to estimate the smoothing parameterλ. Figure 12 shows the smoothed surface generated by the TPS with GCV.For implementation, the CRAN packagergcvpack is used to implement the TPS algorithm (see Table 3 for more information).


## 4.1. Mathematical Preliminaries

In this section, our main interest is not to construct a functionm
(
s
i
) that exactly interpolates the data z
i at n distinct points but to find an attractive way to smooth noisy data. The method of thin-plate splines (TPS) will be used for this purpose.Duchon [12] was the first to build the theoretical foundation for the TPS method. The name TPS comes from the physical situation of bending a thin surface, where the method minimizes the bending energy of a thin plate fixed at the data sites. For our application, the TPS method is used to minimize the cost function:
(9)
J
m
(
f
)
=
∫
R
d
D
m
f
(
s
)
2
d
s
,
where D
m is a vector of partial differentiation operators of order m. In the two-dimensional case, that is, when d
=
2, m
=
2, and f
=
f
(
x
,
y
), the TPS penalty function can be written as
(10)
J
2
f
=
∬
R
2
∂
2
f
∂
x
2
2
+
2
∂
2
f
∂
x
∂
y
2
+
∂
2
f
∂
y
2
2
d
x
d
y
.LetΔ
4
(
x
,
y
) denote the differential operator in the integrand of (10). The thin-plate spline f
(
x
,
y
), which is the solution to the variational problem of minimizing the penalty J
2
(
f
), can then be found by solving the biharmonic equation Δ
4
f
(
x
,
y
)
=
0.The goal is to find the functionf
λ in Sobolev space [13, p. 250] that minimizes the following expression:
(11)
S
(
f
)
=
n
-
1
∑
i
=
1
n
‍
z
(
s
i
)
-
f
(
s
i
)
2
+
λ
J
2
(
f
)
,
where n is the total number of observations, J
2 is a smoothness penalty (the cost function), and λ is the smoothing parameter. The smoothing parameter is a trade-off factor between the rate of change of the residual error and local variation. Optimal minimization of S
(
f
) results in a good compromise between smoothness and goodness of fit.Once the curve approximation of the data has been constructed, generating values at any location, where no observations are available, is accomplished by simply indexing the variablesx and y and fetching the corresponding value. This is a major advantage of smoothing methods over interpolation methods; no extra interpolations are required after the curve has been constructed for a given epoch.Validation is carried out by the GCV (Section2.2.3). Let the cross-validation function g
^
n
,
λ
[
k
], with k
=
1,2
,
…
,
n, be defined as
(12)
g
^
n
,
λ
[
k
]
=
arg
min
⁡
g
^
∈
W
2
(
m
)
∑
i
=
1
n
‍
z
i
-
g
^
x
i
2
+
λ
∫
0
1
‍
g
^
m
x
i
2
d
x
,
where λ
>
0. The OCV (ordinary cross-validation) and OCV MSE (ordinary cross-validation mean square error) V
0
(
λ
), respectively, are defined as
(13)
OCV
λ
=
arg
min
⁡
λ
∈
R
+
V
0
λ
,
V
0
λ
=
1
n
∑
i
=
1
n
‍
g
^
n
,
λ
k
x
i
-
z
i
2
.The determination of the GCV (general cross-validation) goes as follows. First, the expression forV
0 has to be rewritten. There exists an n
×
n matrix A
(
λ
), the smoothing/influence/sensitivity matrix with the property. Consider
(14)
g
^
n
,
λ
(
x
1
)
⋮
g
^
n
,
λ
(
x
n
)
=
A
λ
Z
,
such that V
0
(
λ
) can be written as
(15)
V
0
λ
=
n
-
1
∑
j
=
1
n
‍
a
k
j
z
j
-
z
k
2
1
-
a
k
k
2
,
where k
,
j
∈
{
1,2
,
…
,
n
} and a
k
j is element {
k
,
j
} of A
(
λ
).Definition 1 (generalized cross-validation (GCV)).
LetA
(
λ
) be the smoothing matrix defined in (14); then the GCV function is given by the expression
(16)
GCV
λ
=
n
-
1
I
-
A
λ
z
2
n
-
1
tr
⁡
I
-
A
λ
2
.


## 4.2. Estimation of the Smoothing Parameterλ

The smoothing parameterλ plays a central role in the TPS method. By adjusting the value of λ, one can get the desired level of smoothness at the cost of accuracy at the data sites. When we set this parameter to zero, the problem is reduced to an interpolation with no smoothing. On the other hand, when the smoothing parameter tends to infinity, the method yields a plane which is least-square fitted to the data. The smoothness penalty method can be chosen by any criteria, but the most popular criterion is GCV (generalized cross-validation), also known as the “left-out one” method. The GCV criterion selects the smoothing parameter λ that minimizes the GCV function, equation (16), that is λ
^
=
arg
min
⁡
λ
∈
R
+
⁡
GCV
(
λ
).The GCV functionV
(
λ
) is the predicted mean square error and can be viewed as a weighted version of the OCV
(
λ
)
=
V
0
(
λ
):
(17)
V
λ
=
1
n
∑
i
=
1
n
‍
g
^
n
,
λ
i
x
i
-
z
i
2
W
k
λ
,
W
k
λ
=
I
-
a
k
k
λ
n
-
1
tr
⁡
I
-
A
λ
2
.In geodesy, it is often interesting to estimate the accuracyf
λ. Two loss functions are considered: the mean square prediction error T
(
λ
), and the stricter Sobolev error is defined as W
(
λ
)
=
f
λ
-
f
W
2,
(18)
W
λ
=
∫
a
b
‍
f
λ
-
f
2
d
G
+
∫
a
b
‍
f
λ
m
-
f
m
2
d
G
.
The performance of an estimator is often well characterized by the risk function, defined as the expectation value of the loss function:
(19)
R
T
λ
=
E
T
λ
,
R
W
λ
=
E
W
λ
.In this analysis, the GCV is used to estimate the smoothing parameterλ. Figure 12 shows the smoothed surface generated by the TPS with GCV.For implementation, the CRAN packagergcvpack is used to implement the TPS algorithm (see Table 3 for more information).


## 5.  Numerical Approximation Methods

Numerical approximation techniques will assist us in processing huge data sets with convergence. The main idea is based on the pioneering work of Besag [14].Let us assume that our observations at different locations follow a multivariate Gaussian distribution with meanμ
z and variance-covariance Σ
z. Then the continuously Gaussian fields have the distribution
(20)
f
(
z
)
=
2
π
-
n
/
2
Σ
z
-
1
exp
⁡
-
1
2
z
-
μ
z
T
Σ
z
-
1
(
z
-
μ
z
)
.Approximating the continuous Gaussian random field by the discrete Gauss-Markov random field is accomplished by introducing the Markov property. This is done as follows: we say that two locationss
i and s
j are conditionally independent if and only if
(21)
x
i
⊥
x
j
∣
x
-
(
i
,
j
)
.
This property is very important when constructing the precision matrix Q of the GMRF. That is, if we know what happens nearby, we can ignore everything that lies further away. Consider
(22)
x
i
⊥
x
j
∣
x
-
(
i
,
j
)
⟹
Q
i
,
j
=
0
.That is, element(
i
,
j
) of Q is zero if the process at location i is conditionally independent of a process at j given the process at all locations except {
i
,
j
}. Figure 4 illustrates the concept of the GMRF.Figure 4
Mesh construction on top of which the GMRF is built. Red circles represent the configured reference stations. The mesh in this figure is used to construct the neighborhoods which are key to specifying GMRF conditional independence structure.The sparse precision matrix makes the GMRF computationally effective, but it is difficult to construct reasonable precision matrices. As a conclusion, the GMRF is a Gaussian field with a sparse precision matrixQ
=
Σ
-
1. For an excellent description of the theory and applications of GMRF, the reader is referred to, for example, Rue and Held [15].The integrated nested Laplace approximation (INLA) method developed by Håvard Rue is used to implement the GMRF (see Table3 for more information).


## 6. Kriging Interpolator

The kriging interpolator is a linear spatial interpolation algorithm and is primarily used in geostatistics. In recent years, however, the interpolator has been applied in many new areas, such as geophysics and climate data analysis.Given the observationsZ
(
s
i
)
i
=
1
N, we want to predict the value of Z
(
s
0
) where no observations have been made. Our goal is to find an estimator Z
^
0
=
Z
^
(
s
0
)
=
∑
i
=
1
N
w
i
Z
(
s
i
) such that the following requirements are met.(i) Unbiasedness. This means that E
Z
(
s
0
)
=
E
Z
0 and is accomplished if ∑
i
=
1
N
‍
w
i
=
1 and the mean is stationary.(ii) Minimum Prediction Variance. We make some assumptions about the mean value of the random field Z
(
s
). If the mean is unknown but constant across the entire region of interest, we haveordinary kriging. Otherwise, the method is known assimple kriging.Any estimator that meets the conditions of unbiasedness and minimum prediction variance is said to be a BLUP (best linear unbiasedness predictor). Let us examine the components of the MSPE (mean square prediction error). Consider(23)
σ
S
0
2
=
var
Z
0
-
Z
^
0
=
var
Z
0
+
var
Z
^
0
-
2
cov
⁡
Z
0
,
Z
^
0
=
σ
2
+
var
∑
i
=
1
N
‍
w
i
Z
i
-
2
cov
⁡
Z
0
,
∑
i
=
1
N
‍
w
i
Z
i
=
σ
2
+
∑
i
=
1
N
‍
∑
j
=
1
N
‍
w
i
w
j
cov
⁡
Z
i
,
Z
j
-
2
∑
i
=
1
N
‍
w
i
cov
⁡
Z
i
,
Z
0
=
σ
2
+
∑
i
=
1
N
‍
∑
j
=
1
N
‍
w
i
w
j
C
i
j
-
2
∑
i
=
1
N
‍
w
i
C
i
0
.We want to minimizevar
Z
(
s
0
)
-
Z
^
(
s
0
) subject to the constraint ∑
i
=
1
N
‍
w
i
=
1.The procedure is well defined by the method of Lagrange multipliers. Form the LagrangianL,
(24)
L
=
var
Z
s
0
-
Z
^
s
0
2
-
2
λ
∑
i
=
1
N
‍
w
i
-
1
.We then take the partial derivatives ofL with respect to the weights w
i and to λ, set the equations to be equal to zero, and solve them; we get
(25)
w
1
⋮
w
N
λ
=
C
11
⋯
C
1
N
1
⋮
⋱
⋮
⋮
C
1
N
⋯
C
N
N
1
1
⋯
1
0
-
1
C
10
⋮
C
N
0
1
.Equation (25), which is the kriging equation, is used to compute the weights. The computation of weights is based on the covariances among locations in the sample (region of interest) and the covariances between sample locations and the location to be predicted. To be specific.(1)
Covariances among the locations in the sample:(26)
C
i
j
=
C
s
i
,
s
j
=
Cov
Z
i
,
Z
j
∀
i
,
j
.
The covariance matrix of theN sample values read
(27)
C
=
C
11
⋯
C
1
N
⋮
⋱
⋮
C
1
N
⋯
C
N
N
.
(2)
Covariances between the sample locations and the prediction point:(28)
C
i
0
=
C
s
i
,
s
0
=
Cov
Z
s
i
,
Z
s
0
∀
i
.
The vector of covariances between the sample locations and the prediction point read(29)
c
=
C
10
C
20
⋯
C
N
0
T
.Equation (25) becomes
(30)
w
λ
=
C
1
1
T
0
-
1
c
1
,
where w is 1
×
N vector of weights and 1
=
[
1
⋯
1
] is a vector of the same dimensions.


### 6.1. Directional Effects

Another form of nonstationarity lies in the covariance structure. One specific way to relax the stationarity assumption is to allow directional effects. For instance, the correlation decay rate at increasing distances may be allowed to depend on the relative orientation between pairs of locations. The simplest form of directional effects in the covariance structure is calledgeometrical anisotropy. This arises when a stationary covariance structure is transformed by a differential stretching and rotation of the coordinate axes. Hence, geometrical anisotropy is defined by two additional parameters. Algebraically, a model with geometrical anisotropy in spatial coordinates x
=
(
x
1
,
x
2
) can be converted to a stationary model in coordinates y
=
(
y
1
,
y
2
) by the transformation
(31)
y
1
y
2
=
x
1
x
2
cos
⁡
(
ψ
A
)
-
sin
(
ψ
A
)
sin
(
ψ
A
)
cos
⁡
(
ψ
A
)
1
0
0
ψ
R
-
1
.
ψ
A is called theanisotropy angle and ψ
R
>
1 theanisotropy ratio. The direction with the slowest correlation decay is called theprincipal axis.


### 6.2.  Choice of Covariance Function

The spatial correlation between measurements at different locations is described by the semivariogram functions:(32)
γ
h
=
1
2
var
Z
s
i
-
Z
s
j
=
1
2
var
Z
s
i
+
var
Z
s
j
-
2
cov
⁡
Z
s
i
,
Z
s
j
=
C
0
+
C
h
,
where C
(
0
) is the variance and C
(
h
) is the covariance. The variogram and the covariance contain the same information and can be used interchangeably.In this study, the spatial correlation functionC
(
h
) is defined by the Matérn function and is given by
(33)
C
h
=
σ
2
2
v
-
1
Γ
v
κ
h
v
K
v
κ
h
.h
=
s
i
-
s
j
∈
R
+ is the Euclidean spatial distance between locations. K
v is the modified Bessel function of the second kind [16], the order v
>
0 measures the degree of the smoothness of the process, and κ is the scaling parameter related to the distance of decorrelation (dependency becomes almost 0).C
(
h
) is obtained from spectral densities [17, p. 31] of the form
(34)
f
(
ω
)
=
ϕ
α
2
+
ω
2
-
v
+
1
/
2
v
>
0
,
ϕ
>
0
,
α
>
0
.Figure5 shows the empirical semivariogram γ
(
h
) with the Matérn covariance function, which fits the L1-VTEC data well. It also works in a wide range of circumstances; including low, moderate, and high ionospheric activities, tested with several different reference station configurations, more specifically 75, 100, and 115 stations.Figure 5
An example of the Matérn covariance function, using data with moderate ionosphere activity. Empirical variogram (small circles) and Matérn model fitted by maximum likelihood (solid line). The parameter values used were: smoothnessκ
=
0.5, partial sill σ
2
=
0.185, and range parameter ϕ
=
4.05.


### 6.3. Computation of the Inverse Matrix

The kriging equation (25) requires the inverse of the covariance matrix to be computed, and this is detrimental to the performance of the algorithm for large data sets. The operation may occasionally even fail to invert the matrix. Numerical methods with optimization algorithms will help us avoid this, for instance, factorization methods, ill-conditioned test, and other suitable methods.


## 6.1. Directional Effects

Another form of nonstationarity lies in the covariance structure. One specific way to relax the stationarity assumption is to allow directional effects. For instance, the correlation decay rate at increasing distances may be allowed to depend on the relative orientation between pairs of locations. The simplest form of directional effects in the covariance structure is calledgeometrical anisotropy. This arises when a stationary covariance structure is transformed by a differential stretching and rotation of the coordinate axes. Hence, geometrical anisotropy is defined by two additional parameters. Algebraically, a model with geometrical anisotropy in spatial coordinates x
=
(
x
1
,
x
2
) can be converted to a stationary model in coordinates y
=
(
y
1
,
y
2
) by the transformation
(31)
y
1
y
2
=
x
1
x
2
cos
⁡
(
ψ
A
)
-
sin
(
ψ
A
)
sin
(
ψ
A
)
cos
⁡
(
ψ
A
)
1
0
0
ψ
R
-
1
.
ψ
A is called theanisotropy angle and ψ
R
>
1 theanisotropy ratio. The direction with the slowest correlation decay is called theprincipal axis.


## 6.2.  Choice of Covariance Function

The spatial correlation between measurements at different locations is described by the semivariogram functions:(32)
γ
h
=
1
2
var
Z
s
i
-
Z
s
j
=
1
2
var
Z
s
i
+
var
Z
s
j
-
2
cov
⁡
Z
s
i
,
Z
s
j
=
C
0
+
C
h
,
where C
(
0
) is the variance and C
(
h
) is the covariance. The variogram and the covariance contain the same information and can be used interchangeably.In this study, the spatial correlation functionC
(
h
) is defined by the Matérn function and is given by
(33)
C
h
=
σ
2
2
v
-
1
Γ
v
κ
h
v
K
v
κ
h
.h
=
s
i
-
s
j
∈
R
+ is the Euclidean spatial distance between locations. K
v is the modified Bessel function of the second kind [16], the order v
>
0 measures the degree of the smoothness of the process, and κ is the scaling parameter related to the distance of decorrelation (dependency becomes almost 0).C
(
h
) is obtained from spectral densities [17, p. 31] of the form
(34)
f
(
ω
)
=
ϕ
α
2
+
ω
2
-
v
+
1
/
2
v
>
0
,
ϕ
>
0
,
α
>
0
.Figure5 shows the empirical semivariogram γ
(
h
) with the Matérn covariance function, which fits the L1-VTEC data well. It also works in a wide range of circumstances; including low, moderate, and high ionospheric activities, tested with several different reference station configurations, more specifically 75, 100, and 115 stations.Figure 5
An example of the Matérn covariance function, using data with moderate ionosphere activity. Empirical variogram (small circles) and Matérn model fitted by maximum likelihood (solid line). The parameter values used were: smoothnessκ
=
0.5, partial sill σ
2
=
0.185, and range parameter ϕ
=
4.05.


## 6.3. Computation of the Inverse Matrix

The kriging equation (25) requires the inverse of the covariance matrix to be computed, and this is detrimental to the performance of the algorithm for large data sets. The operation may occasionally even fail to invert the matrix. Numerical methods with optimization algorithms will help us avoid this, for instance, factorization methods, ill-conditioned test, and other suitable methods.


## 7. Performance Parameters

In order to carry out the performance analysis of each individual algorithm, an averaging weighted reference signalX was constructed. It is defined as a linear combination of values generated by algorithms with different weights, that is, X
=
∑
i
=
1
5
‍
w
i
Al
g
i, under the normalization constraint ∑
i
=
1
5
‍
w
i
=
1. Five algorithms are involved to construct the reference signal X.The weights are chosen according to algorithm performance measured in terms of minimum value and stability of variance, functionality, correctness, and processing time. Figure6 shows the variance of two algorithms. We see that ordinary kriging has a minimum and stable variance; therefore its weight is higher than for the Akima bicubic spline. Table 1 summarizes the weight assignment for the algorithms.Table 1
Definition of algorithm weights.
Algorithm type
Weights
Comments
OK
w
1
=
.
25
Produces trustworthy results with MCF
TPS
w
2
=
.
25
A real competitor of kriging, always delivers, even when the covariance structure is poor
GMRF
w
3
=
.
25
Trustworthy results, like OK and TPS
Akima
w
4
=
.
15
Handles small variations perfectly. The weight is reduced compared to OK, GMRF, and TPS
IDW
w
5
=
.
10
No covariance structure is preserved. The assigned weight is reduced compared to OK, GMRF, TPS, and AkimaFigure 6
The ordinary kriging algorithm (BLUP, pink curve) has a stable and minimum variance compared to the AKIMA algorithm (green curve). Normal ionospheric activity (year 2013, doy: 153), 110 reference stations.


### 7.1. Quality of Service Parameters Definitions

For each one of the quality of service (QoS) parameters whose values are negotiable, the worst case performance must be specified. In some cases, the minimum or the maximum values are preferable, in other cases the averaged value.The criteria chosen for performance evaluation in this paper are based on comparing the reference signalX to the output Y
k from algorithm k. Analysis is based on statistical monitoring and detecting the changes in spatial location, scale, and level. The full list is given in Table 2.Table 2
Quality of service (QoS) parameter definitions.
Parameters
Explanations
MAE
Mean absolute error
MSE
Mean square error
RMSE
Root mean square error
NSE
Nash-Sutcliffe efficiency
KGE
Kling-Gupta efficiency
R
Pearson correlation coefficient
β
Changes in locations (mean)
α
Changes in scale (variance)
CT
Computation time is set to 5 seconds. In NRTK processing the values most used are:
ionosphere: 10 seconds
geometrical: 15–20 secondsTable 3
Implementation packages.
Algorithm
Package name
Download
Comments
Kriging
geoR
CRAN
Prediction by OK
(i) AIC computation
(ii) BIC computation
(iii) MLE estimation
(iv) WLS estimation
(v) Anisotropy
GMRF
INLA
CRAN
Prediction by GMRF
TPS
rgcvpack
CRAN
Prediction by TPS
(i) GCV to estimateλ
AKIMA
akima
CRAN
Bilinear interpolation
Bicubic interpolation
IDW
gstat
CRAN
IDW interpolation
Normality
Moments
CRAN
Jarque-Bera normal test
Smoothing
KernSmooth
CRAN
Wand Kernel smoothingThe required length of time series before we can carry out the goodness of fit is a critical parameter. With the data sets used for testing, values in the range of 100–200 epochs were acceptable.All algorithms compete about the QoS, the one with highest score is selected as the winner, and the corrections from this algorithm are used. 12 locations (can be regarded as VRS) are chosen inside the CPOS RTK network for testing, and one location is chosen randomly for each run to compute the QoS. The mathematical definitions of the QoS parameters are given in Table2.(i) Mean Absolute Error (MAE). MAE measures the average absolute error and is defined below. Ideally, this value should be as small as possible. Consider
(35)
MAE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
.(ii) Mean Square Error (MSE). This measures the average squared error and is defined below. This value should also be as close to zero as possible. Consider
(36)
MSE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
2
.(iii) Root Mean Square Error (RMSE). RMSE between reference signals X
i and Y
i gives the standard deviation of the algorithm prediction error and minimum value is preferable. Consider
(37)
RMSE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
2
1
/
2
.(iv) Nash-Sutcliffe Efficiency (NSE) [18]. NSE determines the relative magnitude of the noise variance compared to the observed data variance. Unity means a perfect match, and zero means that algorithm predictions are as accurate as the mean of the observed information, while negative values imply that the observed mean is better than the predicted one. Consider
(38)
NSE
=
1
-
∑
i
=
1
N
X
i
-
Y
i
2
∑
i
=
1
N
X
i
-
X
-
i
2
.(v) Kling-Gupta Efficiency (KGE) [19]. KGE was developed by Gupta et al. as a goodness of fit that decomposes the NSE to facilitate the analysis of correlation, bias, and variability. Consider
(39)
KGE
=
1
-
R
-
1
2
α
-
1
2
β
-
1
2
1
/
2
.Three components are involved in computation of this index.(a)
R is the Pearson product moment correlation coefficient, which ideally should tend to unity. This quantity is defined by the expression
(40)
R
=
N
∑
‍
X
i
Y
i
-
∑
‍
X
i
∑
‍
Y
i
N
∑
X
i
2
-
∑
‍
X
i
2
N
∑
‍
Y
i
2
-
∑
‍
Y
i
2
1
/
2
.
(b)
β represents the change in locations. This index is defined as the ratio between distribution locations (means), and the ideal value is unity. Consider
(41)
β
=
μ
O
μ
C
.
(c)
Variability ratio (VR) represents changes in scale (variances). This index is defined as the ratio between distribution standard deviations, and the ideal value is again unity. Consider(42)
VR
=
σ
C
/
μ
C
σ
O
/
μ
O
.(vi) Computation Time (CT). Algorithms with high quality data and minimum CT are preferable.(vii) Coefficient of Determination (R2). 0
⩽
R
2
⩽
1 and gives the portion of the variance of one variable that is predictable from the other variable.(viii) Spearman Correlation Coefficient
(
ρ
). -
1
⩽
ρ
⩽
1 is a nonparametric test used to measure the degree of associations between two variables.


## 7.1. Quality of Service Parameters Definitions

For each one of the quality of service (QoS) parameters whose values are negotiable, the worst case performance must be specified. In some cases, the minimum or the maximum values are preferable, in other cases the averaged value.The criteria chosen for performance evaluation in this paper are based on comparing the reference signalX to the output Y
k from algorithm k. Analysis is based on statistical monitoring and detecting the changes in spatial location, scale, and level. The full list is given in Table 2.Table 2
Quality of service (QoS) parameter definitions.
Parameters
Explanations
MAE
Mean absolute error
MSE
Mean square error
RMSE
Root mean square error
NSE
Nash-Sutcliffe efficiency
KGE
Kling-Gupta efficiency
R
Pearson correlation coefficient
β
Changes in locations (mean)
α
Changes in scale (variance)
CT
Computation time is set to 5 seconds. In NRTK processing the values most used are:
ionosphere: 10 seconds
geometrical: 15–20 secondsTable 3
Implementation packages.
Algorithm
Package name
Download
Comments
Kriging
geoR
CRAN
Prediction by OK
(i) AIC computation
(ii) BIC computation
(iii) MLE estimation
(iv) WLS estimation
(v) Anisotropy
GMRF
INLA
CRAN
Prediction by GMRF
TPS
rgcvpack
CRAN
Prediction by TPS
(i) GCV to estimateλ
AKIMA
akima
CRAN
Bilinear interpolation
Bicubic interpolation
IDW
gstat
CRAN
IDW interpolation
Normality
Moments
CRAN
Jarque-Bera normal test
Smoothing
KernSmooth
CRAN
Wand Kernel smoothingThe required length of time series before we can carry out the goodness of fit is a critical parameter. With the data sets used for testing, values in the range of 100–200 epochs were acceptable.All algorithms compete about the QoS, the one with highest score is selected as the winner, and the corrections from this algorithm are used. 12 locations (can be regarded as VRS) are chosen inside the CPOS RTK network for testing, and one location is chosen randomly for each run to compute the QoS. The mathematical definitions of the QoS parameters are given in Table2.(i) Mean Absolute Error (MAE). MAE measures the average absolute error and is defined below. Ideally, this value should be as small as possible. Consider
(35)
MAE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
.(ii) Mean Square Error (MSE). This measures the average squared error and is defined below. This value should also be as close to zero as possible. Consider
(36)
MSE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
2
.(iii) Root Mean Square Error (RMSE). RMSE between reference signals X
i and Y
i gives the standard deviation of the algorithm prediction error and minimum value is preferable. Consider
(37)
RMSE
=
1
N
∑
i
=
1
N
‍
X
i
-
Y
i
2
1
/
2
.(iv) Nash-Sutcliffe Efficiency (NSE) [18]. NSE determines the relative magnitude of the noise variance compared to the observed data variance. Unity means a perfect match, and zero means that algorithm predictions are as accurate as the mean of the observed information, while negative values imply that the observed mean is better than the predicted one. Consider
(38)
NSE
=
1
-
∑
i
=
1
N
X
i
-
Y
i
2
∑
i
=
1
N
X
i
-
X
-
i
2
.(v) Kling-Gupta Efficiency (KGE) [19]. KGE was developed by Gupta et al. as a goodness of fit that decomposes the NSE to facilitate the analysis of correlation, bias, and variability. Consider
(39)
KGE
=
1
-
R
-
1
2
α
-
1
2
β
-
1
2
1
/
2
.Three components are involved in computation of this index.(a)
R is the Pearson product moment correlation coefficient, which ideally should tend to unity. This quantity is defined by the expression
(40)
R
=
N
∑
‍
X
i
Y
i
-
∑
‍
X
i
∑
‍
Y
i
N
∑
X
i
2
-
∑
‍
X
i
2
N
∑
‍
Y
i
2
-
∑
‍
Y
i
2
1
/
2
.
(b)
β represents the change in locations. This index is defined as the ratio between distribution locations (means), and the ideal value is unity. Consider
(41)
β
=
μ
O
μ
C
.
(c)
Variability ratio (VR) represents changes in scale (variances). This index is defined as the ratio between distribution standard deviations, and the ideal value is again unity. Consider(42)
VR
=
σ
C
/
μ
C
σ
O
/
μ
O
.(vi) Computation Time (CT). Algorithms with high quality data and minimum CT are preferable.(vii) Coefficient of Determination (R2). 0
⩽
R
2
⩽
1 and gives the portion of the variance of one variable that is predictable from the other variable.(viii) Spearman Correlation Coefficient
(
ρ
). -
1
⩽
ρ
⩽
1 is a nonparametric test used to measure the degree of associations between two variables.


## 8. Implementation and Analyses

Packages used in the implementation are downloaded from the Comprehensive R Archive Network (CRAN). Table3 gives a full description of each package.


### 8.1. Data Preprocessing

The data used in this investigation is ionospheric data obtained from the Norwegian CPOS RTK network. The TEC values are generated using the GFA (geometry-free approach) algorithm. The algorithm takes two steps to process.(a)
From IGS [20], the global ionospheric model (GIM) is available. The ionospheric path delay and differential code biases (DCBs) for each satellite are retrieved from the IONEX (IONosphere map EXchange) format. This information is used to estimate the hardware biases of the reference stations, by using the code observable.
(b)
From the previous step, we then use the biases in phase measurement and compute the ionospheric delay. The procedure is described in more detail in [21].


### 8.2. Interpretation of Results

The test results show that ordinary kriging with the Matérn covariance function is the most appropriate choice under normal circumstances and produces a smooth solution with acceptable accuracy. The Matérn covariance function is well-behaved even for nonstationary fields and is governed by three parameters: location, scale, and shape.Stein [17, p. 12] recommended the use of the Matérn model due to its flexibility (ability to model the smoothness of physical processes and possibility to handle nonstationarity).The processing time of this algorithm increases as the number of observations increase. Another approach is to exclude observations that are far away from the interpolation points and use only a subset of the data for the interpolation. This approach is calledlocal kriging.In order to increase the convergence of OK, we have incorporated the bootstrapping algorithm on historical data to get a very good guess for initial values. Figure9 illustrates the concept.One major obstacle of this algorithm is the computation of the inverse matrix in the kriging equation (25). Using the numerical approximation of the inverse matrix, the computation time will improve considerably, as mentioned previously in Section 6.3.The WLSE (weighted least square estimation) algorithm is preferable to maximum likelihood or restricted maximum likelihood and works in most cases, regardless of the distribution of the observations. If the observations have a Gaussian distribution, the WLS and ML/REML yield the same results.We are often faced with a nonstationary process where we are interested in estimating the spatial covariance for the entire random field. Guttorp and Sampson [22] proposed a two-step approach for solving this problem, a nonparametric algorithm to estimate the spatial covariance structure for the entire random field without assuming stationarity. The interested reader is referred to [23, pp. 93–95].When the covariance structure preserves sparsity, numerical approximation methods are preferable to all other methods, as they require less memory and computation time.TPS algorithm is preferred when performing smoothing rather than interpolating data.


### 8.3. Algorithms Delay Comparison

In this subsection, the delays caused by different algorithms are investigated. The test results are shown in Figure7. The number of observations varied between 700 and 3000 ionospheric piercing points (IPPs). The data used for the plot are the first 255 epochs with a resolution of 30 seconds from high ionospheric activities, from the year 2013 and day of year 152.The left panel shows the CT from global methods, the GMRF (blue curve), TPS (red), and OK (brown). The right panel shows the computation time (CT) caused by the local interpolation methods, the AKIMA (brown) and IDW (blue).
(a)
(b)The delay caused by local methods is shown on the right of Figure7 and is much lower compared to global methods.The GMRF has the highest delay over the OK and the TPS. The only challenge of TPS is to select a good smoothing parameter,λ. The modified cross-validation, the generalized cross-validation, and robust GCV all work well.The IDW methods are local interpolation techniques and use only a subset of the data set to perform the interpolation. The benefit of these methods is the reduced computation time.


### 8.4. QoS Results

Statisticians are usually more interested in smoothing data than interpolating it. When the data is noisy, the TPS smoothing scheme works best. One major advantage of this approach is that once the curve that represents the variation in the data is constructed, we can retrieve the value at any other location without reinterpolating the data set.Figure8 shows the result for an arbitrary prediction point with coordinates (
lon
,
lat
)
=
(
5.0,57.0
). The reference signal X is compared to the predicted values generated by the ordinary kriging algorithm with MCF. The computed quality of service parameters (QoS) are presented below the plot.Figure 8
QoS computation. Time series of predicted values by ordinary kriging (black curve) and the reference signal (blue curve). Goodness of fit parameters are given below the plot.Figure 9
Bootstrapping to accelerate the convergence of OK. Partial sill (σ
2, blue curve) and range parameter (ϕ, orange curve).The results are summarized in Table4 where the QoS parameters are provided for each of the interpolation algorithms tested. An arbitrary epoch has been picked for the test. High scores are highlighted in bold font. The result shows that the ordinary kriging has the best performance. The TPS comes in second place and is the only real competitor to the ordinary kriging for this case. As kriging has the best performance, the corrections from this algorithm will be used to generate synthetic observations of the user in the field. This comparison to determine the best interpolation algorithm is performed for each epoch.Table 4
Computed QoS by different algorithms.
QoS Parameters
Akima  bilinear
Akima  bicubic
TPS
GMRF
Kriging
ME
−0.050
0.030
0.010
−0.030
0.020
MAE
0.090
0.300
0.070
0.090
0.080
MSE
0.020
0.410
0.020
0.020
0.020
RMSE
0.150
0.640
0.130
0.140
0.140
NSE
0.660
0.440
0.790
0.480
0.740
R
2
0.780
0.660
0.820
0.840
0.780
ρ
.920
.930
.940
.950
.950
KGE
0.925
0.589
0.911
0.955
0.976
R
0.958
0.829
0.953
0.959
0.980
β
0.997
1.002
0.994
1.005
1.001
α
0.937
0.626
1.076
1.016
1.013


### 8.5. Practical Considerations

#### 8.5.1. Ordinary Kriging

The library geoR from CRAN (Comprehensive R Archive Network) is used to implement the spatial interpolation. In order to produce correct results, the data analyst must handle many small exceptions that may otherwise result in unexpected output.(i) Performance Analysis. The variance of the estimated grids are analyzed. If the variance is very small, this ensures stability of the algorithm.(ii) Parameter Estimations. In order to increase/accelerate the convergence of ordinary kriging, we have incorporated the bootstrapping algorithm on historical data to get very good estimates of the initial values. Figure 9 shows the first 255 estimates for the parameters partial sill (
σ
2
) and range parameter (
ϕ
), for moderate ionospheric activity, and the network configuration with 75 reference receivers.(iii) Anisotropy. The optimal estimated kriging weights (negative weights) and variances are very sensitive to anisotropy. Our aim is to ensure that the spatial process does not depend on direction. The geometry anisotropy correction is applied by transforming a set of coordinates according to the geometric anisotropy parameters. The package gstat does not provide automatic fitting of anisotropy parameters, while thegeoR package transforms/backtransforms a set of coordinates according to the estimated geometric anisotropy parameters.(iv) Normality Test. The MLE (maximum likelihood estimation) procedure requires that the observations are Gaussian distributed, this assumption is violated in most cases. Therefore the Jarque-Bera test is used as a test of the normality and is based on the third and fourth moments of a distribution, called skewness and kurtosis coefficients; the interested is referred to [24]. If the test fails, the weighted least square method is used to estimate the parameters. Figures 10 and 11 from a configuration with 100 sites and high ionospheric activity confirm that the L1-VTEC distribution is not normally distributed.Figure 10
The quantile-quantile (q-q) plot for visualizing normality test. Jarque-Bera test failed and maximum likelihood is not used to estimate parameters.Figure 11
Nonparametric smoothing with an Epanechnikov kernel is used to determine the L1-VTEC distribution. The distribution is not Gaussian. The weighted least square is used in this case to estimate parameters.Figure 12
L1-VTEC variation surface generated by the TPS algorithm. The modes are well represented in this case.Based on the tests and checks mentioned above, the ordinary kriging is assigned a weight of0
.
25 when computing the QoS values.


#### 8.5.2. Test Results

Figure12 shows the smoothed curve generated by TPS software with GCV. Once the curve is determined, we can easily retrieve any value inside the coverage area without extra computation compared to other interpolation methods. In addition, this describes the variation very well.Once the smoothing parameterλ is determined by GCV, TPS is the real competitor of the kriging algorithm and the weight assigned to it is 0.25.


## 8.1. Data Preprocessing

The data used in this investigation is ionospheric data obtained from the Norwegian CPOS RTK network. The TEC values are generated using the GFA (geometry-free approach) algorithm. The algorithm takes two steps to process.(a)
From IGS [20], the global ionospheric model (GIM) is available. The ionospheric path delay and differential code biases (DCBs) for each satellite are retrieved from the IONEX (IONosphere map EXchange) format. This information is used to estimate the hardware biases of the reference stations, by using the code observable.
(b)
From the previous step, we then use the biases in phase measurement and compute the ionospheric delay. The procedure is described in more detail in [21].


## 8.2. Interpretation of Results

The test results show that ordinary kriging with the Matérn covariance function is the most appropriate choice under normal circumstances and produces a smooth solution with acceptable accuracy. The Matérn covariance function is well-behaved even for nonstationary fields and is governed by three parameters: location, scale, and shape.Stein [17, p. 12] recommended the use of the Matérn model due to its flexibility (ability to model the smoothness of physical processes and possibility to handle nonstationarity).The processing time of this algorithm increases as the number of observations increase. Another approach is to exclude observations that are far away from the interpolation points and use only a subset of the data for the interpolation. This approach is calledlocal kriging.In order to increase the convergence of OK, we have incorporated the bootstrapping algorithm on historical data to get a very good guess for initial values. Figure9 illustrates the concept.One major obstacle of this algorithm is the computation of the inverse matrix in the kriging equation (25). Using the numerical approximation of the inverse matrix, the computation time will improve considerably, as mentioned previously in Section 6.3.The WLSE (weighted least square estimation) algorithm is preferable to maximum likelihood or restricted maximum likelihood and works in most cases, regardless of the distribution of the observations. If the observations have a Gaussian distribution, the WLS and ML/REML yield the same results.We are often faced with a nonstationary process where we are interested in estimating the spatial covariance for the entire random field. Guttorp and Sampson [22] proposed a two-step approach for solving this problem, a nonparametric algorithm to estimate the spatial covariance structure for the entire random field without assuming stationarity. The interested reader is referred to [23, pp. 93–95].When the covariance structure preserves sparsity, numerical approximation methods are preferable to all other methods, as they require less memory and computation time.TPS algorithm is preferred when performing smoothing rather than interpolating data.


## 8.3. Algorithms Delay Comparison

In this subsection, the delays caused by different algorithms are investigated. The test results are shown in Figure7. The number of observations varied between 700 and 3000 ionospheric piercing points (IPPs). The data used for the plot are the first 255 epochs with a resolution of 30 seconds from high ionospheric activities, from the year 2013 and day of year 152.The left panel shows the CT from global methods, the GMRF (blue curve), TPS (red), and OK (brown). The right panel shows the computation time (CT) caused by the local interpolation methods, the AKIMA (brown) and IDW (blue).
(a)
(b)The delay caused by local methods is shown on the right of Figure7 and is much lower compared to global methods.The GMRF has the highest delay over the OK and the TPS. The only challenge of TPS is to select a good smoothing parameter,λ. The modified cross-validation, the generalized cross-validation, and robust GCV all work well.The IDW methods are local interpolation techniques and use only a subset of the data set to perform the interpolation. The benefit of these methods is the reduced computation time.


## 8.4. QoS Results

Statisticians are usually more interested in smoothing data than interpolating it. When the data is noisy, the TPS smoothing scheme works best. One major advantage of this approach is that once the curve that represents the variation in the data is constructed, we can retrieve the value at any other location without reinterpolating the data set.Figure8 shows the result for an arbitrary prediction point with coordinates (
lon
,
lat
)
=
(
5.0,57.0
). The reference signal X is compared to the predicted values generated by the ordinary kriging algorithm with MCF. The computed quality of service parameters (QoS) are presented below the plot.Figure 8
QoS computation. Time series of predicted values by ordinary kriging (black curve) and the reference signal (blue curve). Goodness of fit parameters are given below the plot.Figure 9
Bootstrapping to accelerate the convergence of OK. Partial sill (σ
2, blue curve) and range parameter (ϕ, orange curve).The results are summarized in Table4 where the QoS parameters are provided for each of the interpolation algorithms tested. An arbitrary epoch has been picked for the test. High scores are highlighted in bold font. The result shows that the ordinary kriging has the best performance. The TPS comes in second place and is the only real competitor to the ordinary kriging for this case. As kriging has the best performance, the corrections from this algorithm will be used to generate synthetic observations of the user in the field. This comparison to determine the best interpolation algorithm is performed for each epoch.Table 4
Computed QoS by different algorithms.
QoS Parameters
Akima  bilinear
Akima  bicubic
TPS
GMRF
Kriging
ME
−0.050
0.030
0.010
−0.030
0.020
MAE
0.090
0.300
0.070
0.090
0.080
MSE
0.020
0.410
0.020
0.020
0.020
RMSE
0.150
0.640
0.130
0.140
0.140
NSE
0.660
0.440
0.790
0.480
0.740
R
2
0.780
0.660
0.820
0.840
0.780
ρ
.920
.930
.940
.950
.950
KGE
0.925
0.589
0.911
0.955
0.976
R
0.958
0.829
0.953
0.959
0.980
β
0.997
1.002
0.994
1.005
1.001
α
0.937
0.626
1.076
1.016
1.013


## 8.5. Practical Considerations

### 8.5.1. Ordinary Kriging

The library geoR from CRAN (Comprehensive R Archive Network) is used to implement the spatial interpolation. In order to produce correct results, the data analyst must handle many small exceptions that may otherwise result in unexpected output.(i) Performance Analysis. The variance of the estimated grids are analyzed. If the variance is very small, this ensures stability of the algorithm.(ii) Parameter Estimations. In order to increase/accelerate the convergence of ordinary kriging, we have incorporated the bootstrapping algorithm on historical data to get very good estimates of the initial values. Figure 9 shows the first 255 estimates for the parameters partial sill (
σ
2
) and range parameter (
ϕ
), for moderate ionospheric activity, and the network configuration with 75 reference receivers.(iii) Anisotropy. The optimal estimated kriging weights (negative weights) and variances are very sensitive to anisotropy. Our aim is to ensure that the spatial process does not depend on direction. The geometry anisotropy correction is applied by transforming a set of coordinates according to the geometric anisotropy parameters. The package gstat does not provide automatic fitting of anisotropy parameters, while thegeoR package transforms/backtransforms a set of coordinates according to the estimated geometric anisotropy parameters.(iv) Normality Test. The MLE (maximum likelihood estimation) procedure requires that the observations are Gaussian distributed, this assumption is violated in most cases. Therefore the Jarque-Bera test is used as a test of the normality and is based on the third and fourth moments of a distribution, called skewness and kurtosis coefficients; the interested is referred to [24]. If the test fails, the weighted least square method is used to estimate the parameters. Figures 10 and 11 from a configuration with 100 sites and high ionospheric activity confirm that the L1-VTEC distribution is not normally distributed.Figure 10
The quantile-quantile (q-q) plot for visualizing normality test. Jarque-Bera test failed and maximum likelihood is not used to estimate parameters.Figure 11
Nonparametric smoothing with an Epanechnikov kernel is used to determine the L1-VTEC distribution. The distribution is not Gaussian. The weighted least square is used in this case to estimate parameters.Figure 12
L1-VTEC variation surface generated by the TPS algorithm. The modes are well represented in this case.Based on the tests and checks mentioned above, the ordinary kriging is assigned a weight of0
.
25 when computing the QoS values.


### 8.5.2. Test Results

Figure12 shows the smoothed curve generated by TPS software with GCV. Once the curve is determined, we can easily retrieve any value inside the coverage area without extra computation compared to other interpolation methods. In addition, this describes the variation very well.Once the smoothing parameterλ is determined by GCV, TPS is the real competitor of the kriging algorithm and the weight assigned to it is 0.25.


## 8.5.1. Ordinary Kriging

The library geoR from CRAN (Comprehensive R Archive Network) is used to implement the spatial interpolation. In order to produce correct results, the data analyst must handle many small exceptions that may otherwise result in unexpected output.(i) Performance Analysis. The variance of the estimated grids are analyzed. If the variance is very small, this ensures stability of the algorithm.(ii) Parameter Estimations. In order to increase/accelerate the convergence of ordinary kriging, we have incorporated the bootstrapping algorithm on historical data to get very good estimates of the initial values. Figure 9 shows the first 255 estimates for the parameters partial sill (
σ
2
) and range parameter (
ϕ
), for moderate ionospheric activity, and the network configuration with 75 reference receivers.(iii) Anisotropy. The optimal estimated kriging weights (negative weights) and variances are very sensitive to anisotropy. Our aim is to ensure that the spatial process does not depend on direction. The geometry anisotropy correction is applied by transforming a set of coordinates according to the geometric anisotropy parameters. The package gstat does not provide automatic fitting of anisotropy parameters, while thegeoR package transforms/backtransforms a set of coordinates according to the estimated geometric anisotropy parameters.(iv) Normality Test. The MLE (maximum likelihood estimation) procedure requires that the observations are Gaussian distributed, this assumption is violated in most cases. Therefore the Jarque-Bera test is used as a test of the normality and is based on the third and fourth moments of a distribution, called skewness and kurtosis coefficients; the interested is referred to [24]. If the test fails, the weighted least square method is used to estimate the parameters. Figures 10 and 11 from a configuration with 100 sites and high ionospheric activity confirm that the L1-VTEC distribution is not normally distributed.Figure 10
The quantile-quantile (q-q) plot for visualizing normality test. Jarque-Bera test failed and maximum likelihood is not used to estimate parameters.Figure 11
Nonparametric smoothing with an Epanechnikov kernel is used to determine the L1-VTEC distribution. The distribution is not Gaussian. The weighted least square is used in this case to estimate parameters.Figure 12
L1-VTEC variation surface generated by the TPS algorithm. The modes are well represented in this case.Based on the tests and checks mentioned above, the ordinary kriging is assigned a weight of0
.
25 when computing the QoS values.


## 8.5.2. Test Results

Figure12 shows the smoothed curve generated by TPS software with GCV. Once the curve is determined, we can easily retrieve any value inside the coverage area without extra computation compared to other interpolation methods. In addition, this describes the variation very well.Once the smoothing parameterλ is determined by GCV, TPS is the real competitor of the kriging algorithm and the weight assigned to it is 0.25.


## 9. Conclusion

A significant improvement of the rover position estimation can be achieved by applying the right interpolation/smoothing algorithm at the NRTK interpolation segment. This will reduce the information loss under prediction of the user error level and will provide high quality of virtual reference station data from epoch to epoch.Five methods have been suggested to generate the rover correction. The study shows that the kriging interpolator, the ordinary kriging with the Matérn covariance function, is the most appropriate choice for weighted spatial linear interpolation, while TPS is a strong competitor of OK when the aim is to smooth, not to interpolate, the data. After performing matrix sparsity tests, the GMRF is computationally effective, requires less memory, and produces good results as TPS and OK.For local methods the covariance structure is in general not conserved. For gentle variation in data, the Akima with bicubic method is an appropriate choice because it is the real spline method. While IDW is stable, it is inaccurate and in addition does not conserve the covariance structure of the process under study.One major benefit of these techniques is that there is no need for any prior estimation of the spatial dependency, as in the case of Bayesian analysis (e.g., Kalman filter).


## 10. Discussions

(1)
As we mentioned in the Introduction, processing large data sets is a challenge of the future, and our suggestion for how to handle this is formulated as follows. First of all, we already have enough mathematical tools to do the job, so we do not need to develop new ones. These tools can be considered as elementary building blocks in the hands of the data analyst/modeler. The main challenge is to know that the right tools exist, what they can do for us, what their strengths and weaknesses are, and how to combine them in appropriate ways to describe the observed variations as well as possible. Figure13 shows the number of users connected at the same time and the historical data of the users using the CPOS services in Norway. Both curves increase exponentially in a period of one decade, and if the development continues to follow the same pattern, the existing tools will not be sufficient to process the large data sets.
(2)
Data quality and quantity are important to perform reliable statistical analysis, and elementary checks are necessary before starting analysis.
(3)
In geodesy and geophysical data analysis, the Gauss-Markov model and Kalman Filter are often considered when modeling and when state estimation is necessary. Since new navigation satellite system (e.g., Galileo, Beidou) in addition to the old GPS and GLONASS becomes operation, massive data sets need to be processed in real-time, so we are experiencing a computational paradigm shift.
(4)
To avoid information loss between the correction and interpolation segments, picking the right algorithm for the job is essential for improving the user position errors.Figure 13
Number of Norwegian CPOS RTK users (blue curve) increased from 1 to ~3000 and the users connected to the system at the same time (red curve) increased from 1 to ~500.



---

*Source: 346498-2015-02-19.xml*