1 Introduction

The velocities of the Eurasian plate are described in the area of Europe by means of the permanent station velocities. These sites constitute a so-called velocity model. Interpolation of these discrete points to the remaining areas not covered by the measurement points for the continuous velocity model seems to be an important issue. Such a model will be created in this study for the area of Poland using the velocities obtained from the permanent stations located in Europe, which record the observations from the Global Positioning System (GPS). Worth pointing out is the fact that several of the softwares tested in this research are the black boxes. Therefore the authors compared only the results of a few methods of velocity interpolations rather than concentrating on the algorithms, in an attempt to find the one most appropriate to the present geological situation. The geological situation of the area of Poland, that is, how each of the structures behaves, was a kind of verification of an applied method of interpolation.

The velocities of Central Europe have been determined since the 1990s by means of the epoch campaigns. The GPS campaigns CEGRN (Central European GPS Geodynamic Reference Network—1994, 1995, 1996, 1997, 1999, 2001) were designed to investigate the intraplate velocities of the region between the Baltic Sea and the Mediterranean Sea (Becker et al. 2002; Rogowski and Hefty 1998). In the period of 11 years, eight CERGOP GPS observation campaigns were executed. The velocity vectors were determined for stations which participated in GPS observation campaigns in 1994–2005 (Klek and Rogowski 2006). Several Local Analysis Centres (LAC) took part in the processing aimed at the estimation of the velocities of the permanent and other occupied sites based on the time series of coordinates (Bogusz et al. 2002). From that time, several permanent sites were established. The velocities determined from the long-span time series of coordinates are becoming more and more reliable, but there is still a limited number of factors which can diminish this reliability. Some of them could be real, geophysical effects (Bogusz et al. 2011, 2012a), however some may be the artefacts of the satellite system itself. Among the possible mechanisms Penna and Stewart (2003) showed how non-modelled analysis errors in tidal frequencies in the diurnal and semi-diurnal bands can be efficiently aliased to longer periods. Dong et al. (2002) conducted a breakdown by the different origins of the oscillations in the GPS observations into three groups related to real and artificial effects. A similar analysis concerning daily solutions was presented by Ray et al. (2008). They focused on the non-linear residuals generated in the ITRF2005 combination from weekly combined global frames produced by IGS (International GNSS Service). The Polish Ground-Based Augmentation System (GBAS) ASG-EUPOS has been operating since 2008. As a part of the European project EUPOS, it conducts permanent analyses for stability in position changes for geodynamical studies as well as for the realisation of the ETRS89 reference system (Figurski et al. 2010). Additionally, we should remember that the vertical changes in Central Europe are on the order of single millimetres per year, which corresponds to the accuracy of their determination. The known models of recent vertical movements of the Earth's crust surface in the area of Poland were processed to make the verification on the basis of independent (satellite) measurements. The model was created using the satellite data and it was confronted with earlier (levelling) models. Some differences in the analysed models have been noticed, which indicates the defect of reference distinctly resulting from sea level variability (Kontny and Bogusz 2012).

Interpolation of unknown values on the area among the measured points may be realised by means of at least several methods. The most commonly known techniques include: the Moving Average, Thiessen’s polygons, Triangulation (TIN), Kriging and Spline. Kreemer et al. (2003) made a global model of both velocities and strain rates. They performed a least-squares fitting between the model and observed geodetic velocities as well as modelled and observed geological strain rates. The modelled data were interpolated over a spherical Earth using bi-cubic Bessel splines whereas the observed data contained 3,000 velocity values. Kreemer et al. (1998) used the earthquake strain rates derived from 74 events to estimate a long-term velocity field for the zone between the North American and Pacific plates. The velocity field was estimated for two proposed kinematic models: the first implied the presence of a micro-plate and the other a transform deformation zone (pseudo-plate). The authors obtained the velocity field by the inversion of the earthquake strain rates. Beavan and Haines (2001) derived continuous horizontal velocity and strain fields from 362 GPS stations distributed across New Zealand. The velocity field model was interpolated with the bi-cubic spline functions defined within a curvilinear grid that covered the New Zealand area and was extended into the Australian and Pacific plates. Holt et al. (2000) obtained a velocity field by joint inversion of strain rates of over 230 GPS stations’ velocities. They performed a Kostrov summation of fault slip rates and in this way inferred strain rates within the area of Asia. They then fitted the strain rates with a continuous velocity field using bilinear spline interpolation. Kierulf et al. (2012) obtained the continuous velocity field for Norway. They interpolated the velocities from 140 GNSS stations with the means of Kriging interpolation and the use of the GIA forward model. The authors stated that the statistical interpolation method gave better results inside the analysed network while the geophysical model estimated the velocities better on the edge of the network. Majdanski (2012), using data from POLONAISE’97, CELEBRATION 2000 and SUDETES 2003 experiments, made a 3D model of the velocity distribution in the crust and upper mantle by the interpolation of 2D models and over 30 profiles. He compared the different interpolation techniques and concluded that Kriging is precise enough to model the tectonic events on the area of Poland.

2 Plate Kinematic Models for Calculating the Intraplate Velocities

The plate motion model can be described as the rotation around the Euler pole with the spherical coordinates \( \varphi_{E} \) and \( \lambda_{E} \). The angular velocity ω for each plate p is estimated from (Altamimi et al. 2007; Drewes 2009):

$$ X_{\text{l}} =\omega \times X_{\text{i}} $$
(1)

where X l is a point’s velocity vector and X i is a position vector for each point.

Velocity components v n and v e of the velocity vector can be described by (Hefty 2007):

$$ v_{\text{n}} = \omega \times \sin (\lambda_{\text{P}} - \lambda_{\text{E}} ) \times \cos \varphi_{\text{E}} $$
(2)
$$ v_{\text{e}} = \omega \times [\cos \varphi_{\text{E}} \times \sin \varphi_{\text{P}} \times \cos (\lambda_{\text{P}} - \lambda_{\text{E}} ) - \sin \varphi_{\text{E}} \times \cos \varphi_{\text{P}} ] $$
(3)

The comparison of Euler parameters between different plate motion models for the Eurasian plate (Fig. 1) is presented in Table 1.

Fig. 1
figure 1

The Euler parameters

Full size image
Table 1 The Euler parameters for different plate motion models
Full size table

The discrepancies between velocities’ reductions using the various plate motion models on the example of the South-Western part of Poland are shown in Fig. 2.

Fig. 2
figure 2

Intraplate velocities obtained using various plate motion models

Full size image

3 Data

The Polish Active Geodetic Network (ASG-EUPOS) is the multifunctional precise satellite positioning system established by the Head Office of Geodesy and Cartography in 2008. The adjusted network consists of over 130 Polish foreign incorporated sites.

A project concerning the construction of modules of support and development of the ASG-EUPOS system was started in 2010 (Figurski et al. 2011). One of the tasks is to construct a geodynamic module for interpreting the results of processing satellite observations for contemporary research of the Earth’s dynamics in the area of Central Europe. A key issue of this module is to develop a continuous model of velocities based on measurements from permanent stations. The studies made so far concerned horizontal (Bogusz et al. 2012b) and vertical velocity fields (Kontny and Bogusz 2012) in a discrete form.

The determination of the time series was made using the Bernese version 5.0 GPS Software (Dach et al. 2007) and normal equations from the EPN (EUREF Permanent Network) and ASG-EUPOS (more than 300 GNSS sites presented in Fig. 3) sites. The detailed description was given in the paper by Bogusz et al. (2012b).

Fig. 3
figure 3

Location of permanent stations recording GNSS observations in Europe

Full size image

The correct determination of the intraplate velocities is doubtful when based upon using various models of plate motion. In our studies we have used the NUVEL-1A-NNR geological model (De Mets et al. 1990) for calculating the residua of velocities expressed in the ITRF2005 reference frame (International Terrestrial Reference Frame) (Altamimi et al. 2007).

Figure 4 presents histograms of the permanent stations’ velocity residua of the ASG-EUPOS system for the North–South component (N) and the East–West (E) component for the NUVEL (Northwest University Velocity) tectonic plate motion model, used for interpolation (the horizontal axis denotes the stations’ velocity residua in mm/year while the vertical axis shows the number of observations with the same value). Most of the observations have values within the range of 0–3 mm/year for the N component and within the range of −2 to 0 mm/year for the E component. The value distributions comply with the normal distribution (red line). This may indicate the stability of the lithosphere within the Eurasian plate. The statistics in the form of the maximum and minimum values (Table 2) were prepared for the observations involving the area of Poland in order to make a comparison between the data and the results of interpolation made by the means of the specific methods.

Fig. 4
figure 4

Histograms of the ASG-EUPOS permanent stations’ velocity residua for the NUVEL model (left: North–South component, right: East–West component)

Full size image
Table 2 The statistics of the ASG-EUPOS permanent stations’ velocity residua for the NUVEL model
Full size table

Figure 5 presents the velocity residua of the ASG-EUPOS system stations with respect to the NUVEL model. Most of the intraplate velocities of the presented permanent stations are directed to the northwest. Anomalous directions of the velocity residua are indicated by green ellipses. An anomalous small value of the velocity residuum for the DRWP (Drawsko Pomorskie) station is indicated by the navy blue ellipse while orange ellipses indicate stations with velocity residua both smaller than others and differently directed.

Fig. 5
figure 5

Intraplate velocities of the ASG-EUPOS stations

Full size image

4 The Tested Software

Three packages of software were tested within the research: GMT (The Generic Mapping Tools), Golden Software’s Surfer package and ArcGIS. GMT (Wessel and Smith 1998) is free software for computation results visualising in a wide range of projections and for various kinds of results processing, e.g. filtering or interpolating. GMT is also equipped with a vast free GIS database, e.g. coastlines, rivers or state borders. The software saves the produced maps in the PostScript (PS) or the Encapsulated PostScript (EPS) formats according to the user’s command input from the command line level. It also enables the user to write scripts and automate the tasks run by the software (Wessel and Smith 2011). Golden Software’s Surfer package (license number: Golden Software’s Surfer package 10.0, WIG WAT, nr WS-113545-ta4t) is a programme developed for visualising data; therefore it is often used for creating maps and terrain surface modelling. A wide set of embedded interpolation algorithms enables us to select an optimum one for the input data set (Golden Software, Version 10.0). ArcGIS (license number: ESRI_SENTINEL_KEY = 3712999727004) is a software designed for creating, editing, browsing and analysing geographic information in the form of maps and data organised in tables. The programme enables the user to analyse data from the point of view of e.g. their spatial distribution or relations between them. GIS software is now not only a tool for the visualisation of the calculations, but also allows for advanced statistical analysis of phenomena in relation to their location on the Earth’s surface (Burrough 2001; Burrough and Mc Donnell 1998; Johnston et al., 2001).

The analysis of various interpolation methods and programmes is carried out to determine the range of the possible solutions. The analyses of vertical velocities of the Earth’s crust on the area of Poland conducted by Kowalczyk (2009) indicated the significance of this problem.

5 The Investigated Interpolation Methods

The following methods of data interpolation were tested within the research project:

  1. (a)

    in the GMT programme: the Nearest Neighbor method, the Triangulation method (TIN) and the Surface method;

  2. (b)

    in the Surfer programme: the Minimum Curvature method, Inverse Distance to a Power and Kriging;

  3. (c)

    in the ArcGIS programme: Spline Interpolation method and the Kriging method and the most appropriate semi-variogram model was selected.

Selection of the optimum interpolation method should depend on: the explicitness of the interpolation algorithms, complexity and time needed for completing interpolation within the specific programme, and the capability to compute the errors of the interpolated values. Worth pointing out is the fact that in the described study the computed errors are not based on the errors of velocity vectors determined by the GNSS observations, but are rather the evidence of the efficiency of the selected interpolation method.

Some of the tested methods have an inexplicit interpolation function. When a selected method is implicit we cannot say for sure how it really works. It is very important to know how the algorithm calculates the interpolated values or their errors to compare it with other tested methods. Some of the described methods are popular ones, so we can suppose that the algorithm looks a certain way. The other ones are not well known and it is hard to say if their way of calculation is appropriate to our goal (how they treat the velocities and the distances between given points and the interpolation nodes), thus their efficiency can only be estimated by comparing the results.

5.1 The Minimum Curvature Method

Figure 6 presents the results of interpolating the examined data by means of the Minimum Curvature method, which is analogous to a thin, linearly elastic plate passing through each of the values with a minimum amount of bending. This method generates the smoothest possible surface while attempting to honour the data as closely as possible (Golden Software). Anomalous values of the interpolated residua of the velocity appeared in the North-East (red ellipse). They may be caused by the lack of data beyond the visualised area—the method adjusts to all of the existing data but it fails in areas where the number of data significantly decreases. The green ellipses indicate the impact of the anomalously directed velocity residua of the permanent stations on the interpolated nodes, while the navy blue and orange ellipses show the impact of the small values of the velocity residua of the permanent stations on the obtained results. The Minimum Curvature method was tested in the Surfer programme. It has an inexplicit form of the interpolation function, and the programme does not enable us to compute the errors of the interpolated values.

Fig. 6
figure 6

Interpolation by means of the Minimum Curvature method

Full size image

5.2 The Kriging Method

The Kriging method is a linear estimation algorithm (Webster and Olivier 2007). The estimated value V in an interpolated point is a sum of the products of the weight λi and the values vi in the data points:

$$ V{ = }\sum\limits_{\text{i}}^{N} {\lambda_{\text{i}} \times v_{\text{i}} } $$
(4)

The estimation is unbiased when the weights are made to sum to 1:

$$ \sum\limits_{\text{i}}^{N} {\lambda_{\text{i}} } { = 1} $$
(5)

The weights are calculated from the Kriging equations:

$$ \lambda \,{ =\, }A^{ - 1} \times {\text{b}} $$
(6)

where:

\( {\text{A = }}\left[ {\begin{array}{*{20}c} {\gamma (d_{ 1 , 1} )} & {\gamma ( {\text{d}}_{ 1 , 2} )} & \ldots & {\gamma (d_{ 1 ,N} )} & 1\\ {\gamma (d_{ 2 , 1} )} & {\gamma ( {\text{d}}_{ 2 , 2} )} & \ldots & {\gamma (d_{ 2 ,N} )} & 1\\ \ldots & \ldots & \ldots & \ldots & \ldots \\ {\gamma (d_{N , 1} )} & {\gamma ( {\text{d}}_{N , 2} )} & \ldots & {\gamma (d_{N,N} )} & 1\\ 1& 1& \ldots & 1& 0\\ \end{array} } \right] \) is a matrix of semivariances depending on the distance between the data points,

\( \lambda\, {=}\left[ {\begin{array}{*{20}c} {\lambda_{ 1} } \\ {\lambda_{ 2} } \\ \ldots \\ {\lambda_{N} } \\ \psi \\ \end{array} } \right] \) is a vector of weights and ψ is a Lagrange multiplier,

\( {\text{b = }}\left[ {\begin{array}{*{20}c} {\gamma (d_{\text{i,1}} )} \\ {\gamma (d_{\text{i,2}} )} \\ \ldots \\ {\gamma (d_{{{\text{i,}}N}} )} \\ 1\\ \end{array} } \right] \) is a vector of the semivariances depending on the distance between the interpolated point and the measurements. In this work, the stationarity of the input data was assumed. The semivariance function γ(di) is adjusted based on empirical semivariance values calculated from the input data.

Estimation variance (specifically Kriging variance) may be obtained as:

$$ {{\upsigma}}^{ 2} \,{=}\,\sum\limits_{\text{i = 1}}^{\rm N} {\lambda_{\text{i}} \times \gamma (d_{\text{i}} ) { + }\,\psi } $$
(7)

In the Surfer programme the Kriging method is realised by defining a model of the semivariogram that best adjusts to the location of the processed data. A semivariogram is a statistically-based, quantitative description of a surface’s roughness. It is a function of a separation vector: this includes both distance and direction, or Δx and Δy. The variogram function yields the mean dissimilarity between points separated by the specified vector (dissimilarity is measured by the squared difference in the Z values) (Golden Software).

All available semivariogram models in the Surfer software were tested, and the linear one defined as (Golden Software):

$$ \gamma ( {\text{h) = C(h)}} $$
(8)

was empirically recognised as the best. C here is the vertical scale for the specific components of the variogram, while h is the relative separation distance used in an isotropic setting, computed by the following equation:

$$ {\text{h = }}\frac{{\sqrt {\Updelta {\text{x}}^{ 2} \,{ + }\,\Updelta {\text{y}}^{ 2} } }}{\text{A}} $$
(9)

The Length (A) parameter defines how rapidly the variogram components change with the increasing separation distance. The Length (A) parameter for a variogram component is used to scale the physical separation distance. [Δx Δy] is the separation vector (in map coordinates).

Figures 7 and 8 present a comparison of interpolation results obtained by means of the Kriging method using various semivariogram models. The linear model was empirically recognised to be the best. The remaining models, i.e. the rational square, gave unexpected values of the interpolated velocity residua. The navy blue ellipse indicates the impact that the small velocity residuum for the DRWP station has on the obtained interpolation results. The impact of the anomalously directed intraplate velocities of the system stations has been recorded for only one case (green ellipse). This means that the Kriging method is immune to local disturbances of the velocity field, perfectly adjusts to the spatial distribution of the values of the velocity residua and displays the trends in the data.

Fig. 7
figure 7

Interpolation by means of the Kriging method with the linear model of the semivariogram

Full size image
Fig. 8
figure 8

Interpolation by means of the Kriging method with the rational square model of the semivariogram

Full size image

An advantage of the Kriging method in the Surfer programme is the capability to compute the errors of the interpolated values. The errors are computed as the difference between the interpolated and the measured values (Golden Software):

$$ \Updelta_{\text{int}}\,{ = }\,v_{\text{int}} - v_{\text{meas}} $$
(10)

where:

Δint :

an interpolation error,

v int :

residuum of the interpolated velocity,

v meas :

measurement result

The following parameters have influence on the errors of the interpolated values: the number of sectors (in our study, 1 sector was used) into which the programme divides the area into the vicinity of the interpolation node, the maximum number of points from each and all sectors (in both cases the entire number of data) which the programme uses for interpolation, the minimum number of measurements in all sectors (1 datum), and the number of empty sectors around the interpolation node above which the node is attributed to zero (1 sector). Moreover, the various radiuses of the area from which the data are used for interpolation, called the vicinity of the interpolated node, may be considered.

Fig. 9
figure 9

Kriging with the linear variogram model, the interpolated values and the interpolation errors

Full size image

The Kriging method has an inexplicit interpolation function but the Surfer programme computes the errors of the interpolated values, which is a significant advantage of the method (Fig. 9).

5.2.1 The Development of the Kriging Model

The Kriging method has the capability to generalise and predict trends in velocities, and thereby remove the local effects that are not related to the tectonic movements. The sensitivity of the method tested in ArcGIS depends on its parameters that were input at several steps of the modelling. In the development phase of empirically selecting the appropriate semivariogram, the optimum interval (lag size) and the maximum distance (number of lags) for which the semivariance is computed have to be selected.

The semivariance indicates the degree of similarity between point values at a given distance h and has been calculated from a set of values based on the equation (Burrough 2001; Zhang and Yao 2008):

$$ \gamma (h ) = \frac{ 1}{ 2N} \times \sum\limits_{\text{i = 1}}^{N} {\left( {\rm V (x_{\text{i}} )- V (x_{\text{i + h}} )} \right)} $$
(11)

where N is the number of samples, V(xi) and V(xi+h) are the measured values at a distance of h. The theoretical model of the semivariogram that fits the empirical data the best has been selected on the basis of a cross-validation procedure, which allows one to evaluate how accurately the model is able to determine the approximated values (Webster and Olivier 2007). Several models of semivariogram have been tested and the adjustment was assessed by three statistics:

  • the mean error ME:

    $$ {\text{ME = }}\frac{ 1}{N} \times \sum\limits_{\text{i = 1}}^{N} {\left( {V (x_{\text{i}} )- \widehat{V} (x_{\text{i}} )} \right)} $$
    (12)

  • the mean squared error MSE:

    $$ {\text{MSE = }}\frac{ 1}{N} \times \sum\limits_{\text{i = 1}}^{N} {\left( {V (x_{\text{i}} )- \widehat{V} (x_{\text{i}} )} \right)^{ 2} } $$
    (13)

  • the mean squared deviation ratio MSDR:

    $$ {\text{MSDR = }}\frac{ 1}{N} \times \sum\limits_{\text{i = 1}}^{N} {\frac{{\left( {V ( {\text{x}}_{\text{i}} )- \widehat{V} (x_{\text{i}} )} \right)^{ 2} }}{{\widehat{\sigma }^{ 2} (x_{\text{i}} )}}} $$
    (14)

The cross-validation numerical results are presented in Tables 3 and 4 where parameters such as mean error, root mean squared error and root mean squared deviation ratio are given. The most frequently used models of semivariograms (Webster and Olivier 2007) like circular, pentaspherical, spherical, exponential and Gaussian were tested. The three models that fit the data the best are presented in the Tables 3 and 4.

Table 3 Cross-validation numerical results (North component—VN)
Full size table
Table 4 Cross-validation numerical results (East component—VE)
Full size table

Based on the cross-validation numerical results, the semivariogram model (Gaussian model) has been selected, with the following values:

  • mean error is close to zero (−0.015 mm/year for the N component and −0.002 mm/year for the E component),

  • root mean squared error is the smallest (0.708 mm/year for the N component and 0.696 mm/year for the E component),

  • root mean squared deviation ratio is close to one (1.095 for the N component and 1.118 for E component).

The function of the semivariances in the Gaussian model (Burrough 1986; Mc Bratney and Webster 1986) has the form of:

$$ \gamma ( {\text{h)}} = C_{ 0} \,{ + }\,C ( 1- e^{{ - (\frac{\text{h}}{\text{a}} )^{ 2} }} ) $$
(15)

where C 0 is a nugget value, C is a sill, h is the lag and a is the range.

In Figs. 10, 11, 12, 13 the semivariograms for V N (Figs. 10, 11) and V E (Figs. 12, 13) components are shown in direction of the minimum and maximum ranges. In the corners of the figures the pie charts of the empirical semivariance in square cells that are one lag wide have been placed. The pie charts show changes of empirical semivariance values with regard to the direction, which justifies the use of the anisotropic model. The applied model of anisotropy is a zonal variant, for which semivariograms for the specific directions have different ranges. The ellipses show the range in the model in a given direction.

Fig. 10
figure 10

The values of the semivariances for the North component in the NW–SE direction

Full size image
Fig. 11
figure 11

The values of the semivariances for the North component in the perpendicular direction

Full size image
Fig. 12
figure 12

The values of the semivariances for the East component in the NE-SW direction

Full size image
Fig. 13
figure 13

The values of the semivariances for the East component in the perpendicular direction

Full size image

Semivariogram parameters (C 0—nugget value, C—partial sill, h—lag, and a—range (minimum, maximum)) for the V N component are C 0 = 0.362 mm/year, C = 0.256 mm/year, h = 55 km, a is from 324 to 440 km respectively. For the V E component parameters take the following values: C 0 = 0.316 mm/year, C = 0.334 mm/year, h = 55 km and a ranges from 324 to 440 km.

5.2.1.1 The North Component V N

For the V N component of the velocity, the semivariances change in a smaller range in the direction of NW–SE (Fig. 10) while there is a significant increase in the perpendicular direction (Fig. 11). The values of the semivariances are calculated using the points located inside the ellipse of the semi-axes of 440 km and 324 km divided into 8 sectors. The azimuth of the longer axis of the ellipse is 164°.

5.2.1.2 The East Component V E

For the V E component of the velocity, the semivariances change in a smaller range in the direction of NE-SW (Fig. 12) while there is a significant increase in the perpendicular direction (Fig. 13). The values of the semivariances are calculated using the points inside the ellipse of the semi-axes of 440 km and 324 km divided into 8 sectors. The azimuth of the longer axis of the ellipse is 52°.

Application of the anisotropy model enables the better fitting of the model to the data. It takes into consideration the differentiation of the velocity with respect to both the distance and direction.

The directions of the ellipses coincide with the directions of the velocity gradients. For the V N component, the longer axis of the ellipse (A = 164°) coincides with the direction of the smaller gradient, and the shorter one coincides with the direction of the greater gradient. It is similar for the V E component while the longer axis’ azimuth is 52°.

It is worth mentioning that the azimuths of the ellipses agree with the location of the faults in the T–T (Teisseyre-Tornquist) Zone (approx. azimuth is 132°) and the respective perpendicular directions. It is probable that they will coincide with the main deformations that will be determined in the area of Poland.

Figures 14 and 15 present the velocities’ components in the ITRF2005 reference frame with the assumption of the isotropic and anisotropic models. In the mentioned figures the variation of velocities on the area of Poland is shown. From the analysis of the isolines it appeared that the values of north velocities’ components VN decrease in the North-East direction (the South-West part of Poland moves slower than the North-East one) whereas the VE components increase in the South-East direction (the North-West part of Poland moves slower than the South-East one).

Fig. 14
figure 14

The North velocity component with the assumption of the isotropic (left map) and the anisotropic (right map) model

Full size image
Fig. 15
figure 15

The East velocity component with the assumption of the isotropic (left map) and the anisotropic (right map) model

Full size image

The differences among the results in modelling were found by analysing the distances between the isolines. In the anisotropic model the distances decrease. The greater gradient appeared in the main directions, and the interpolated points were calculated in these directions with the use of the weights (customised to the specifics of the phenomenon) that differed from the ones that were used in the isotropic model.

5.3 The Spline Interpolation Method

The Spline Interpolation method is a minimum curvature function, which passes through all data points with the accuracy of their mean errors. In the data points the function’s curvature has the minimum value, while the course of the function between the points is close to linear. All data points have their influence on the value of the interpolated point. Each station is involved in determining the vectors “a” of the function’s coefficients described by Eq. (18) separately for each component. The number of the coefficients is equal to the number of stations. The closest stations have the greatest impact on the value of the interpolated point and the effect decreases logarithmically. In this model, the isotropy of the input data was assumed.

The Spline used for interpolation is described by the following equation (Osada 1995):

$$ {\text{f}(n,e,a) = a}_{{L + 1}} + a_{{L{ + 2}}} \times n + a_{{L{ + 3}}} \times {{e + }}\frac{ 1}{ 2} \times \sum\limits_{\text{i = 1}}^{\text{L}} { [ (N_{\text{i}} - n )^{ 2} { + (}E_{\text{i}} - e )^{ 2} ]} \times {\text{ln[(}}N_{\text{i}} - n )^{ 2} { + (}E_{\text{i}} - {{e)}}^{ 2} ]\times a_{\text{i}} $$
(16)

where:

n, e :

horizontal coordinates of the interpolated point;

Ni, Ei :

horizontal coordinates of the i-th permanent station.

Function parameters ai were determined from solutions of the system equations with the following conditions:

$$ \sum\limits_{\text{i = 1}}^{\text{L}} {a_{\text{i}} { = 0}} ,\sum\limits_{\text{i = 1}}^{\text{L}} {a_{\text{i}} \times n_{\text{i}} { = 0}} ,\sum\limits_{\text{i = 1}}^{\text{L}} {a_{\text{i}} \times e_{\text{i}} { = 0}} $$
(17)

Parameters a i were obtained separately for V N and V E from the following equation:

$$ a\,{ = }\,A^{ - 1} \times f $$
(18)

where:

a :

vector of the function parameters,

A :

matrix of the coefficients,

f :

vector of the permanent stations’ velocities.

The variance–covariance matrix of the a i parameters is calculated using the following equation:

$$ C_{\text{a}} \,{ = }\,A^{ - 1} \times C_{\text{V}} \times (A^{ - 1} )^{\text{T}} $$
(19)

To assess the accuracy of the interpolated velocities values V, the error propagation law is used:

$$ {{\upsigma}}_{\text{V}}^{ 2} \,{ = }\,F^{\text{T}} \times C_{{{\text{a}}_{\text{N}} }} \times F $$
(20)

The velocity field obtained by means of the Spline Interpolation shows well-matched trends between the nearest ASG-EUPOS sites (Fig. 16), which is necessary to interpret the local effects of the ground surface deformations. To obtain regional trends of horizontal movements, it is necessary to exclude points with clearly different values and velocity directions or use a method of data filtering. Examples of the influence of local movements on the value of the velocity field may be observed particularly: in the southern part near the WODZ (Wodzisław Śląski) and KATO (Katowice) points (orange ellipses), in the northwestern part near the DRWP point (navy blue ellipse) and in the eastern part near the CHEL (Chełm) point. The green ellipses show the great influence of the unexpected sites’ velocities on the interpolated values (the proof of fitting the method through all data points).

Fig. 16
figure 16

The Spline method and its errors

Full size image

6 Comparison of Results

The comparison of interpolation results obtained by means of all the above discussed methods enabled the authors to select the optimum one for creating a continuous model of the velocity field. Table 5 presents the interpolation results of the tested methods.

Table 5 Comparison of interpolation results
Full size table

The differences between all the results obtained by means of the tested interpolation methods and the data are at the level of a few tenths of mm/year. The Nearest Neighbor method was rejected by the authors despite the fact that the interpolation results are very similar to the measurement data. Its disadvantage is locality—it does not take into consideration the spatial variability of the discussed phenomenon—and discontinuity—the nearest measurement data have the greatest impact on the interpolated nodes without taking the farther stations into consideration. The disadvantage of the Triangulation method is the fact that the interpolation results depend only on three measurement data. Therefore, this method will not be used for interpolation in the creation of the continuous model of the velocity field. The Inverse Distance to a Power method does not take into consideration, as the Nearest Neighbor method does, the spatial variability of the discussed phenomenon. Interpolation by means of the Kriging method, due to connecting points of similar values, perfectly displays the spatial variability of the investigated phenomenon. Additionally, it enables the computation of the errors of the interpolated values, which makes it the optimum method for the interpolation of investigated observations for the continuous model of the velocity field. Interpolation by means of the Surface method yields results similar to the measurement data. Despite this, it has not been selected by the authors as the optimum method for interpolating the velocity residua of permanent stations in the area of Poland because changing the values of its parameters does not result in changes in the interpolated values. The Spline method interpolates the velocity residua for the continuous model perfectly adjusting to the measurement data. Additionally, the ArcGIS software enables one to compute the errors of the interpolated values which, due to this feature of the method, are minimum. However, the method has a disadvantage—it does not take into consideration the spatial variability of the investigated phenomenon—the measurement datum closest to the interpolated node has the greatest impact on the interpolated value.

Figure 17 presents the interpolation results of the tested methods. The maximum values of the components N and E of the intraplate velocities with respect to the NUVEL model were interpolated using the Minimum Curvature method because it adjusts the interpolated surface to the measurement data to the largest extent and the distances between the reference stations outside Poland significantly increase (the interpolated velocity residua adjusted to the nearest data). The minimum value of the interpolated E component of the velocity residua was obtained using the Surface method. The changing of the parameters describing the method in the GMT programme did not result in the change of the interpolation results. Results most similar to the measurement data were obtained by means of the Nearest Neighbor method—the interpolated nodes depend only on the closest measurement datum. The Minimum Curvature, Kriging and Spline Functions methods were recognised as optimum for interpolating the velocity residua of the permanent stations. Kriging is immune to local disturbances of residua and it displays the trend of the investigated phenomenon in the discussed area very well. Spline Functions interpolation adjusts perfectly to all data. These features of the mentioned interpolation methods determine the values of the obtained velocity residua and interpolation errors in the nodes. Empirical analysis of the interpolation results obtained using the Minimum Curvature method proved that except for the North-Eastern and South-Western areas, the method yielded the identical results as the Kriging method (Fig. 18). The differences between these two methods may result from a significant decrease of measurement data density outside Poland (the Minimum Curvature method adjusts to measurement data).

Fig. 17
figure 17

Comparison of the discussed interpolation methods (mm/year)

Full size image
Fig. 18
figure 18

Comparison of the Minimum Curvature and the Kriging methods

Full size image

The measurement data are marked in red, the interpolation results using the Kriging method are marked in green, and the interpolation results using the Minimum Curvature method are marked in blue. The indicated errors concern the interpolation using the Kriging method.

To compare the two interpolation methods, Kriging and Spline, the differences between the interpolated velocity residua were determined:

$$ \Updelta v_{\rm N} = v_{\rm N}^{\text{kriging}} - v_{\rm N}^{\text{spline}} $$
(21)
$$ \Updelta v_{\text{E}}\, { = }\,v_{\text{E}}^{\text{kriging}} - v_{\text{E}}^{\text{spline}} $$
(22)

Additionally, errors in the interpolated nodes were compared using the following relations:

$$ m_{{{{\Updelta}}_{\text{int}} \rm N}} = \pm \sqrt {\left( {\frac{{\partial \Updelta v_{\rm N} }}{{\partial v_{\rm N}^{\text{kriging}} }}m_{{v_{\rm N}^{\text{kriging}} }} } \right)^{ 2} { + }\left( {\frac{{\partial \Updelta v_{\rm N} }}{{\partial v_{N}^{\text{spline}} }}m_{{v_{\rm N}^{\text{spline}} }} } \right)^{ 2} } { = } \pm \sqrt {m_{{v_{\rm N}^{\text{kriging}} }}^{ 2} + m_{{v_{\rm N}^{\text{spline}} }}^{ 2} } $$
(23)
$$ m_{{{{\Updelta}}_{\text{int}} {\text{E}}}} \,{ = }\, \pm \sqrt {\left( {\frac{{\partial \Updelta v_{\text{E}} }}{{\partial v_{\text{E}}^{\text{kriging}} }}m_{{v_{\text{E}}^{\text{kriging}} }} } \right)^{ 2} { + }\left( {\frac{{\partial \Updelta v_{\text{E}} }}{{\partial v_{\text{E}}^{\text{spline}} }}m_{{v_{\text{E}}^{\text{spline}} }} } \right)^{ 2} } { = } \pm \sqrt {m_{{v_{\text{E}}^{\text{kriging}} }}^{ 2} { + }\,m_{{v_{\text{E}}^{\text{spline}} }}^{ 2} } $$
(24)

where \( m_{{{{\Updelta}}_{\text{int}} N}} \) is the value of the ellipse’s semi-axis in the South-North direction, \( m_{{{{\Updelta}}_{\text{int}} {\text{E}}}} \) the value of the ellipse’s semi-axis in the West-East direction.

Figure 19 presents the differences between the interpolation results obtained by means of the Kriging and Spline methods. Except for the Southern area of Poland, the differences between the interpolated values of the velocity residua Δv are slight—about 1 mm/year. The significant change of the differences in Southern Poland may be caused by anomalous values of the velocity in the measurement data in the vicinity of the interpolation nodes and by other characteristic features of the compared methods. The differences in the semi-axes of the error ellipses of the interpolated nodes are larger than 0.5 mm/year.

Fig. 19
figure 19

Differences between the Kriging and spline interpolation methods

Full size image

7 Summary

All the previous research on obtaining the continuous velocity field applied Kriging or Delaunay triangulation to the different GPS data sets. Grenerczy (2003) created the continuous velocity plate field based on data from regional networks, for the whole period of each station activity. The following networks were considered in the studies: CEGRN (the Central European GPS Geodynamic Reference Network), CRODYN (the Croatian and Slovene Geodynamic Network), HGRN (the Hungarian GPS Geodynamic Reference Network), IGS (The International GNSS Service) and SAGET (Satellite Geodynamic Traverses). The continuous velocity field was created by interpolating data for the whole region of Central Europe with two interpolation methods: Kriging and Delaunay triangulation with linear interpolation. Both methods gave similar effects. The obtained intraplate velocity values ranged from 0.1 up to 5 mm/year depending on the region of interpolation. The produced GPS-derived contour map of crustal velocity magnitudes in Central Europe was in very good agreement with the earthquake activity map. Grenerczy (2004) described the accuracy of GPS network assumption in the Hungarian area and made the new method of analysis independent of the scale coefficient between frames. Linear regression and its errors were checked on a few baselines of the HGRN network. Statistical tests were made and their results were favourable with respect to values, correlation and standard deviations of estimated parameters. A velocity map based on the reference station data was made, without making any unnecessary transformations, just minimising the velocities on the stable part of the European plate, using seven stations. Based on the map of velocities, the map of the continuous velocity field was made by interpolating velocities via the Kriging method with the nugget effect. Hefty (2007) demonstrated a method of combining regional and local velocity fields into a homogenous system consistent with the ITRF2000 and then reduced for the APKIM2000 model. The described horizontal velocities resulted from seven various regional and local analyses of 192 GPS sites by means of epoch campaigns and permanent observations from the region of Central Europe, Adriatic region and the Balkan Peninsula. Several characteristic blocks with significant horizontal velocity trends for the Central and Southeastern Europe were proved. The two remarkable movements oriented northwards and southwards were shown with the magnitude of 3÷5 mm/year for the Adriatic region and 3 mm/year for the Southern Carpathian Mountains and Eastern part of the Balkan Peninsula. The results of the interpolation were presented in 2° × 1° and 1° × 0.5° grids with 2σ error ellipses.

This research concerned testing seven methods for the interpolation of the intraplate velocities by means of three programmes. The Nearest Neighbor, Triangulation and Surface methods were tested by means of the free GMT programme. The GMT programme does not enable one to compute the errors of the interpolated values, and its interpolation algorithms are inexplicit. The Minimum Curvature, Inverse Distance to a Power and Kriging methods were tested by means of the Surfer programme. Surfer’s interpolation algorithms are also inexplicit. The possibility of computing errors of the interpolated values is available only for the Kriging method. Data processing using the discussed software is neither complicated nor a time-consuming task. Selection of the optimum interpolation method should be dependent on the explicitness of the interpolation function and on the capability to determine the errors of the interpolated values. The Spline method was tested using the ArcGIS software. It yielded values of interpolated intraplate velocities similar to the Kriging method. However, the errors in the interpolation nodes were much smaller because this method perfectly adjusts to the measurement data which is its disadvantage—it does not take into consideration the spatial distribution of the investigated phenomenon. Statistics in the form of maximum, minimum and mean values of the velocity residuum components were prepared for the measurement data and for all interpolation methods. The Kriging Function was selected to be used for modelling the continuous velocity field in the area of Poland on the basis of the analyses and conducted tests. In the authors’ opinion, Kriging is the best method of modelling velocity residua to obtain trends of the investigated phenomenon. This method is immune to local disturbances of the velocity field. The difference between the arithmetic mean of the components of velocity residua computed from the measurement data and those interpolated by means of the Kriging method is at the level of hundredth of a millimeter/y. Interpolation by means of the Green i in the GMT programme yields almost identical results as Kriging, and the values of the interpolated intraplate velocities with respect to the NUVEL model using Spline Functions proved compatibility in independent investigations with interpolated values using the Kriging Functions at the level of 1 mm/year. The Minimum Curvature method may alternately be applied for the Kriging method. Empirical analysis of the results proved that both methods yield identical results for appropriate density of measurement data. The impact of the small intraplate velocity at one site (DRWP) is observed in all results of the tested interpolation methods. Anomalous directions of velocity residua at some of the permanent stations cause differences from the North-East trend of the velocity residua vectors. This occurs on a greater or a smaller area depending on the characteristics of the tested interpolation method.

This research on the interpolation method was based on the residua from the NUVEL-1A NNR model, which seems to be suitable for such application. However, on the basis of other work conducted by the authors it could be concluded that for geokinematic applications the APKIM method (The Actual Plate Kinematic and Deformation Model) is more appropriate. This enabled the authors to develop a continuous intraplate velocity field (Figs. 20, 21). In the future, a discontinuous one will be created in which interpolation only within the geologic units of Europe will be performed for interpretation purposes. Moreover, the creation of a continuous velocity field on the area of Europe will lead to obtaining the continuous strain field in future research.

Fig. 20
figure 20

Intraplate velocities using APKIM2005 (IGN) model

Full size image
Fig. 21
figure 21

The continuous intraplate velocity field for Europe with the use of APKIM2005 (IGN) model

Full size image