Introduction

Today, Global Navigation Satellite System (GNSS) measurements, in particular those from the Global Positioning System (GPS), are fundamental to many geodetic and geophysical investigations (Kreemer et al. 2014; Métivier et al. 2014) and are frequently used during the construction of kinematic reference frames, such as, for example, the International Terrestrial Reference Frame 2014 (ITRF2014) (Altamimi et al. 2016). From the processing or re-processing of the observables of permanently installed GNSS stations, daily or weekly geocentric coordinate solutions are obtained, from which position time series are formed. The primary product from the analysis of these time series is often the linear rate of change, or velocity, and the associated uncertainty (Zhang et al. 1997).

The velocities are assumed to represent the linear movement of the earth’s crust due to tectonic plate motions (Larson et al. 1997; Drewes 2009; Altamimi et al. 2012) or the viscoelastic relaxation associated with glacial isostatic adjustment (Johansson et al. 2002; Bradley et al. 2009). In addition, almost all sites within the global network of GPS stations also show nonlinear periodic motions, which have been associated primarily with the seasonal effects on earth. However, some stations may also exhibit nonlinear and non-periodic character motions (Shih et al. 2008; Bogusz 2015). This may be due to the elastic response of the earth’s crust from the rapid ice mass loss at the polar ice sheets or mountain glaciers (Wahr et al. 2013), being located in the deforming zones near plate boundaries (Wdowinski et al. 2004) or areas of oil, gas, or groundwater extraction (Munekane et al. 2004). In this study, we will not deal with cases showing such nonlinear and non-periodic behavior.

Blewitt and Lavallée (2002) were one of the first to investigate the effect of periodic signals on GPS time series, and today it is widely acknowledged that such seasonal terms affect GPS and other geodetic time series (Bos et al. 2010; Davis et al. 2012; Bogusz and Figurski 2014). The causes are almost completely understood and can be grouped into categories suggested by Dong et al. (2002): real geophysical effects of atmospheric (Tregoning and van Dam 2005), hydrological (van Dam et al. 2001) or ocean loadings (van Dam et al. 2012) with thermal expansion (Romagnoli et al. 2003; Xu et al. 2017) and numerical artifacts of navigation satellite systems. Penna and Stewart (2003) described aliased periodic signals in the coordinate time series due to under-sampling of residual diurnal and semidiurnal tidal signatures. Griffiths and Ray (2013) extended this analysis for the largest waves in the International Earth Rotation and Reference Frame Service (IERS) Conventions 2010 diurnal and semi-diurnal tidal polar motion model, assuming 24-h sampling. The second technique-related error is associated with satellite orbits (draconitics). Agnew and Larson (2007) found that for daily sampling rates in GPS-derived coordinates this period will alias to a frequency of 1.04333 cpy (cycles per year). Ray et al. (2008) compared harmonics obtained using techniques such as GPS, Very Long Baseline Interferometry, and Satellite Laser Ranging when they discovered an anomalous peak in the GPS-derived time series, which was not present in other series. They explained it to be related to the interval required for the constellation to repeat its inertial orientation with respect to the sun (GPS year). Amiri-Simkooei (2013), analyzing Jet Propulsion Laboratory (JPL) data, obtained a period of 351.6 ± 0.2 days. Abraha et al. (2017) proved that most of the power in draconitic period is satellite induced. Finally, multipath (King et al. 2012), insufficient modeling of antennas (Sidorov 2016), errors in the network, the adjustment that transfers from fiducial stations, or the inclusion of the scale (Tregoning and van Dam 2005) contributes to constellation-specific periodic signals in the time series.

From the above, it is clear that the deterministic model of a GNSS position time series needs to include parameters for both the linear and periodic motions at a given station (Bevis and Brown 2014). Furthermore, in a true geodetic approach, the model must provide means to obtain the most realistic uncertainties associated with the parameter estimates in order to provide confidence limits at a given significance level. In this respect, we model the time series, which have been pre-processed for outliers, offsets, and gaps, with least squares estimation (LSE) as:

$$ x\left( t \right) = x_{0} + v \cdot t + \sum\limits_{i = 1}^{n} {\left[ {{\text{AC}}^{i} \cdot \cos \left( {\omega_{i}^{T} \cdot t} \right) + {\text{AS}}^{i} \cdot \sin \left( {\omega_{i}^{T} \cdot t} \right)} \right]} + \varepsilon \left( t \right), $$
(1)

where x 0, v, ACi, ASi, n are the intercept, velocity, cosine, and sine terms of ith harmonic and the number of harmonics of the angular velocity ω of period T, respectively. The term ε(t) contains the residuals and all variations not explicitly modeled and, therefore, disregarded by this station motion model. It is now commonly known that the residuals do not follow a strict random behavior, i.e., a white noise with spectral index κ = 0, but follow that of a combination of colored and white noise, with the latter often being modeled as a power-law process (Agnew 1992). This temporally correlated noise is of a form of flicker noise with spectral index κ = − 1 or of random-walk process with spectral index κ = − 2. Flicker noise is present in the GPS position time series due to mismodeling in GNSS satellites orbits, Earth Orientation Parameters, large-scale atmospheric or hydrospheric effects. It has been already described by Williams et al. (2004) or Amiri-Simkooei et al. (2007) who stated that a combination of flicker and white noise is the preferred noise model for GPS data. Teferle et al. (2008) emphasized that their time series showed evidence of power-law noise close to the flicker. Caporali (2003) stated that the power spectral densities of time series prove that flicker noise is preferred for most stations to represent the data at frequencies below 6 cpy. For frequencies higher than 6 cpy, the spectrum tends to become a white noise. Santamaria-Gomez et al. (2011) concluded that a combination of power-law and white noise is the preferred one for a global data set. Klos et al. (2015) stated that this combination is valid also for weekly data. A random-walk process might be present in GPS position time series due to local environment and monumentation. Johnson and Agnew (1995) were the first to state that the type of monument may have an impact on the long-term correlation present in GPS data. Beavan (2005) compared two monuments: a concrete pillar and a deep drilled braced monument and stated that there is no significant difference in the stochastic part of GPS position time series collected. However, the analysis was only performed for 4.5 years of data. Even when a perfect monument is considered with the environment being transparent to the signal, the noise would be flicker and arise from low-frequency fluctuations of the satellite clocks (Dutta and Horn 1981). Beyond these features, there are also some seasonal features that add more correlated noise to GNSS data. As stated by Johnson and Agnew (1995), if time series were too short, it is more difficult to detect any change due to correlation than it is for long data. As is now widely acknowledged, the stochastic properties of the residuals significantly influence the magnitude of the uncertainties associated with the parameters (Zhang et al. 1997; Langbein and Johnson 1997; Mao et al. 1999; Williams 2003), and the whiter the noise, the faster the velocity uncertainty decays with increasing time series length. Therefore, the noise character has the greatest impact on velocity uncertainty and errors of the deterministic model estimated at the same time.

We provide a general dilution of precision (GDP) of the velocity uncertainties being the ratio of uncertainties of velocities arising from two different assumptions of the deterministic model. The first of them assumes a linear velocity only, while the second is a combination of linear velocity and periodic components. Each of the assumptions is presented with certain noise models starting from white, moving on to flicker and end up with random-walk, the most extreme case for GNSS position time series. The determined errors of the velocities are being discussed along the two previous papers of Blewitt and Lavallée (2002) and Bos et al. (2010) which provided the first overview of this topic. Using simulated data we show that our approach has additional advantages to GNSS time series analysis in comparison to the ones mentioned before. We deliver the estimates of a combined effect of periodic signals and noise model assumed in residuals, to show how much one may underestimate the uncertainties of velocities when inappropriate assumptions were made prior to the analysis. Finally, we apply the newly developed formulae to real GNSS time series of selected ITRF2014 stations. Beyond the GNSS position time series, our conclusions are also valid for any other type of position time series, as those collected by Very Long Baseline Interferometry (VLBI), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), Satellite Laser Ranging (SLR), or tide gauges.

Dilution of precision

Blewitt and Lavallée (2002) developed a model to calculate the bias level when one does not account for periodic signals of annual frequency. The velocity bias expressed in mm/year introduced by them was a zero-crossing oscillatory function of data span, tending to zero for infinitely long time series. According to this function they discovered the “zero-bias theorem” for unbiased velocity near integer-plus-half years and introduced the ratio of uncertainties of velocities v 1 and v 2. They called it the dilution of precision (DP) and estimated its value as (Blewitt and Lavallée 2002):

$$ {\text{DP}}\left( \tau \right) = \left[ {1 - \frac{6}{{\left( {\pi \cdot f \cdot \tau } \right)^{2} }} \cdot \frac{{\left( {\cos \left( {\pi \cdot f \cdot \tau } \right) - \frac{{\sin \left( {\pi \cdot f \cdot \tau } \right)}}{\pi \cdot f \cdot \tau }} \right)^{2} }}{{1 - \frac{{\cos \left( {\pi \cdot f \cdot \tau } \right) \cdot \sin \left( {\pi \cdot f \cdot \tau } \right)}}{\pi \cdot f \cdot \tau }}}} \right]^{{ - \frac{1}{2}}} , $$
(2)

where τ is the time span and f denotes the frequency of the periodic terms. Figure 1 shows the DP value which increases toward infinity when the time span is shorter than 1 year. A number of maxima of oscillations in DP can be noticed for integer years. These come from periodic terms in (1) defined by f. The time span of observations is given by τ. For certain epochs starting from 1.5 years with a step of 1-year the DP = 1. It results from the fact that the velocity uncertainties for a model with linear velocity and a model with seasonal terms added are equal to each other. Also, for a time span longer than 3.5 years, the DP < 1.05, which means that the difference between both variances is below 5%. However, Blewitt and Lavallée (2002) inputted only white noise, and as a consequence assumed that the character of a stochastic part ε has a little impact on the estimated uncertainties.

Fig. 1
figure 1

(recomputed from Blewitt and Lavallée 2002)

Dilution of precision for models: with linear velocity and linear velocity plus annual term.

Full size image

It is expected that when we add a power-law noise, the difference between both variances will need more time to fall below 5%, while the power-law noise process adds temporal correlation to the time series. That is why Bos et al. (2010) discussed the results of Blewitt and Lavallée (2002) by empirically analyzing the effect of periodic signals using six stations, but assumed a combination of colored and white noise. They noticed that the way the stochastic part is described might be of greater importance than that of the seasonal part in the deterministic model. This is due to the fact that velocity uncertainty depends mostly on the spectral index and amplitude of the power-law process. They concluded that the knowledge of the noise characteristics of GNSS time series is crucial when velocities and their uncertainties are determined. Furthermore, they observed a shift in the minimum of the DP from the integer-plus-half years toward integer-plus-a-quarter years position. We confirm the results found by Bos et al. (2010) and also deliver the mathematic formulae of DP for time series with power-law noise combined with white one.

In order to estimate the DP, it is first necessary to obtain an expression for the covariance of the estimates in terms of the noise covariance. For this, let us consider the vector of residuals ε defined as the difference between the model parameters and the data fitted into (1):

$$ {\varvec{\upvarepsilon}} = {\mathbf{A\theta }} - {\mathbf{x}} $$
(3)

where A is the model or design matrix, \( {\varvec{\uptheta}} \) is a vector with parameters of the model, x stands for the observational data or measurements, and ε is the vector of residuals.

When the parameters of the model are determined by means of maximum likelihood estimation (MLE; Langbein 2012) and assuming a Gauss distribution for such estimates, they are computed as:

$$ {\hat{\mathbf{\theta }}} = \left( {{\mathbf{A}}^{T} \cdot {\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1} \cdot {\mathbf{A}}} \right)^{ - 1} \cdot {\mathbf{A}}^{T} \cdot {\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1} \cdot {\mathbf{x}} $$
(4)

with C εε being the covariance matrix for the residuals ε.

Then, using (3) to replace x in (4), and rearranging terms, it yields:

$$ {\hat{\mathbf{\theta }}} - {\varvec{\uptheta}} = \left( {{\mathbf{A}}^{T} \cdot {\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1} \cdot {\mathbf{A}}} \right)^{ - 1} \cdot \left( {{\mathbf{A}}^{T} \cdot {\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1} } \right){\kern 1pt} \,{\varvec{\upvarepsilon}} $$
(5)

which can be used to compute the second moment or covariance of the residuals in terms of the covariance noise as follows:

$$ {\mathbf{C}}_{{\hat{\theta }\hat{\theta }}} = \left( {{\mathbf{A}}^{T} \cdot {\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1} \cdot {\mathbf{A}}} \right)^{ - 1} . $$
(6)

In general, the vector of the residuals ε is a linear combination of colored and white noise:

$$ \varepsilon_{i} = \sum\limits_{j = 0}^{N - 1} {h_{j} \cdot w_{i - j} } + v_{i} $$
(7)

where the first term on the right side is a convolution of white noise, \( w_{i} \in {\rm N}\left( {0,\sigma_{pl} } \right) \), and the second term is a white noise process, \( v_{i} \in {\rm N}\left( {0,\sigma_{wn} } \right) \), with σ pl and σ wn standing for the amplitude parameters of the colored and white noise processes, respectively. N is the number of data, while h is defined by the recursive formula (Bos et al. 2008):

$$ \begin{aligned} h_{0} & = 1 \\ h_{i} & = \left( {\frac{ - \kappa }{2} + i - 1} \right)\frac{{h_{i - 1} }}{i},\quad \, \forall i > 0 \\ \end{aligned} $$
(8)

with κ being the spectral index ranging from −3 to 1, which, when appropriately assigned, may characterize the stochastic part of the time series (Mandelbrot and Van Ness 1968).

From (7) the covariance of the residuals is:

$$ {\mathbf{C}}_{\varepsilon \varepsilon ;kl} = \sigma_{pl}^{2} \sum\limits_{j = 0}^{l} {h_{j} \cdot h_{{j + \left| {k - l} \right|}} } + \sigma_{wn}^{2} \cdot \delta_{kl} $$
(9)

or, in matrix notation:

$$ {\mathbf{C}}_{\varepsilon \varepsilon } = \sigma_{pl}^{2} \cdot {\mathbf{L}} \cdot {\mathbf{L}}^{T} + \sigma_{wn}^{2} \cdot {\mathbf{I}} $$
(10)

L is a lower triangular Toeplitz matrix with coefficients defined as follows:

$$ L_{ij} = \left\{ {\begin{array}{*{20}c} {h_{i - j} } & {\forall \left( {i - j} \right) \ge 0} \\ 0 & {\forall \left( {i - j} \right) < 0} \\ \end{array} } \right. $$
(11)

where \( h_{i - j} \) are the same coefficients as in (8) that depend on the spectral index κ. Finally, by inserting (10) in (6), we get the following expression for the covariance matrix of the estimated parameters, \( {\hat{\mathbf{\theta }}} \):

$$ {\mathbf{C}}_{{\hat{\theta }\hat{\theta }}} = \left[ {{\mathbf{A}}^{T} \cdot \left( {\sigma_{pl}^{2} \cdot {\mathbf{E}}\left( \kappa \right) + \sigma_{wn}^{2} \cdot {\mathbf{I}}} \right)^{ - 1} \cdot {\mathbf{A}}} \right]^{ - 1} $$
(12)

where E is the part of the covariance matrix from the power-law process that depends on the spectral index \( \kappa \):

$$ {\mathbf{E}}\left( \kappa \right) \equiv {\mathbf{L}} \cdot {\mathbf{L}}^{T} . $$
(13)

Thus, the covariance of the estimates is expressed in terms of the spectral index \( \kappa \) and the amplitudes of the colored and white noise processes, σ pl and σ wn , respectively.

Finally, note that the process defined by (7) and (8) is a fractionally differenced Gaussian noise process with spectral density (Hosking 1981):

$$ P(f) = \frac{{2^{\kappa } \sigma_{pl} }}{{ ( {\text{sin(}}\pi \cdot f ) )^{ - \kappa } }} $$
(14)

which, as \( f \to 0 \), yields a power-law:

$$ P(f) \propto f^{\kappa } $$
(15)

Henceforth, we will refer to the color noise process as power-law throughout the paper.

Models

If we assume that the deterministic model of GNSS time series follows a linear trend with a vector of parameters built as:

$$ {\varvec{\uptheta}} = \left[ {x_{0} ,v} \right]^{T} , $$
(16)

where v is the slope or velocity and x 0 is the intercept, the model matrix A is created as:

$$ {\mathbf{A}} = \left[ {\begin{array}{*{20}c} 1 & {t_{0} } \\ \vdots & \vdots \\ 1 & {t_{N - 1} } \\ \end{array} } \right]. $$
(17)

Then, the variance of velocity estimate can be computed by inserting (17) into (12).

The proper modeling of the seasonal terms may directly influence on the reliability of determined parameters, such as the velocity. Indeed, as Blewitt and Lavallée (2002) already showed, the seasonal terms may influence the uncertainty of velocity. However, they did not consider the power-law character of residuals but assumed them to follow a white noise process. Bos et al. (2010) noticed that not only do the seasonal terms affect the linear velocity when being improperly removed, but the noise properties are much more important for a reliable estimation of the velocity error. These directly affect the uncertainty of velocity due to different shapes of the covariance matrix, as in (9–10).

Let us assume now a deterministic model with linear velocity and annual term:

$$ {\mathbf{A}} = \left[ {\begin{array}{*{20}c} 1 & {t_{0} } & {\cos \left( {2 \cdot \pi \cdot f \cdot t_{0} } \right)} & {\sin \left( {2 \cdot \pi \cdot f \cdot t_{0} } \right)} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & {t_{N - 1} } & {\cos \left( {2 \cdot \pi \cdot f \cdot t_{N - 1} } \right)} & {\sin \left( {2 \cdot \pi \cdot f \cdot t_{N - 1} } \right)} \\ \end{array} } \right]. $$
(18)

The vector of parameters is built here as:

$$ {\varvec{\uptheta}} = \left[ {x_{0} ,v,AC,AS} \right]^{T} , $$
(19)

where AC and AS are the annual cosine and sine terms, respectively. Similar to the previous model, substitution of (18) in (12) yields the covariance for the estimates in (19).

Simulated series

In order to test the general dilution of precision formulas derived above, we carried out a number of evaluations. For these, we simulated time series of up to a maximum length of 25 years, which is at the time of writing the longest term of available GNSS time series. Synthetic data were created based on two approaches for the deterministic part: The first included annual and semiannual terms, and the second included all significant periods, i.e., all tropical and draconitic terms up to ninth harmonic, plus fortnightly and the Chandlerian period since the latter may be present at some stations at island and coastal sites (Richard Gross, private communication, 2015); see Bogusz and Klos (2016) for more details. Then, we added temporal correlation in the form of a combination of power-law and white noise, built a noise covariance matrix as in (9) and transferred it to estimate the variances of the estimated parameters as in (12). We assumed four different spectral indices equal to 0, − 1, − 1.5, and − 2 that indicate pure white, pure flicker, fractional Brownian motion, and pure random-walk noise, respectively, with \( \sigma_{pl}^{2} = 1 \) and \( \sigma_{wn}^{2} = 1 \). White noise was added to flicker, fractional Brownian motion, and random-walk when assumed.

Figure 2 shows the error of the velocity for white, flicker and random-walk noise. For time series shorter than 70 days, the error of velocity is much higher for white and flicker noise assumptions than for random-walk noise. This situation changes when the data becomes longer than 70 days. The error of velocity is the smallest for the white noise assumption. Random-walk noise delivers the greatest errors of parameters. The longer the data is, the greater is the difference between errors delivered for random-walk noise and flicker or white noise assumptions. Having assumed white, flicker and random-walk noises, for 20 years of data one will obtain an error of velocity equal to 10−3, 10−2, and 10−1 mm/year, respectively.

Fig. 2
figure 2

Error of velocity \( \sigma_{v}^{{}} \) (in mm/year) for different lengths of time series in double logarithmic scale. The integer spectral indices are examined. We assumed that \( \sigma_{pl}^{2} = 1 \) and \( \sigma_{wn}^{2} = 1 \)

Full size image

We estimated the relative differences in velocity variances for two deterministic models: with linear velocity denoted as \( \sigma_{v1}^{2} \) and with linear velocity plus seasonal terms denoted as \( \sigma_{v2}^{2} \), as

$$ \Delta \sigma_{v}^{2} = \frac{{\sigma_{v1}^{2} - \sigma_{v2}^{2} }}{{\sigma_{v1}^{2} }}. $$
(20)

Adding seasonal terms results in oscillations in the estimated velocity error when the length of time series changes. The velocity variances are computed from the diagonal terms of the covariance matrix (12). These oscillations, of course, are not so obvious any more when the spectral index of power-law dependencies increases, as was also noticed in Bos et al. (2010). The largest differences between two models are being observed for the random-walk noise assumption. As the spectral index was shifted from 0 toward − 2, the differences between velocity uncertainties will, of course, enlarge much more. Due to the above, the largest differences between two models are being expected for badly monumented stations (Langbein and Johnson 1997; Beavan 2005; Klos et al. 2016).

The uncertainties estimated will differ when seasonal terms are included or are not included. This difference is computed as a ratio between both values. We called it the general dilution of precision (GDP) in accordance with Blewitt and Lavallée (2002) and as in (2), but taking into consideration a power-law noise in the stochastic part, as was underlined in (15). We have adopted two approaches to estimate the seasonal part. First, the widely used annual and semiannual terms were subtracted from the GPS position time series. Then, the approach, which is consistent with Bogusz and Klos (2016), was assumed: tropical and draconitics up to their ninth harmonics, plus Chandlerian and fortnightly terms were modeled. Figures 3 and 4 show a GDP for white, flicker, random-walk and fractional Brownian motion of spectral indices equal to 0, − 1, − 2, and − 1.5, respectively.

Fig. 3
figure 3

General dilution of precision (GDP) for white noise (blue), flicker noise (red), fractional Brownian motion of spectral index equal to − 1.5 (orange) and random walk (black) plotted for a deterministic model containing a linear velocity plus annual and semiannual oscillations

Full size image
Fig. 4
figure 4

General dilution of precision (GDP) for white noise (blue), flicker noise (red), fractional Brownian motion of spectral index equal to − 1.5 (orange), and random walk (black) plotted for a deterministic model containing a linear velocity plus extended model of periodicities of all tropical and draconitic terms up to ninth harmonic plus fortnightly and the Chandlerian period, according to Bogusz and Klos (2016)

Full size image

We computed the GDP values for harmonics of 1 year. They increased the local maxima in the GDP plot. It can be noticed that in comparison to Blewitt and Lavallée (2002) who considered the annual term, the semiannual signal also increases the local minimum of GDP with a white noise assumption. It means that adding the power-law noise to the annual curve or adding a semiannual term to white noise causes an increase in the velocity uncertainty even at those points where the estimated velocity should not be biased. The more seasonal terms are added to the series, the more biased is the velocity uncertainty, especially for short time scales. Having compared Figs. 3, 4, and 5, it is clear that the type of deterministic model affects the velocity uncertainty and makes GDP to reach the value of 5% (we adopted this number following Blewitt and Lavallée 2002) after 9 years, rather than 4, as was expected for annual plus semiannual terms (a case of white and flicker noise). The value of 5% means an increase in velocity uncertainty of 0.025 mm/year when a typical error of velocity of 0.5 mm/year is considered (Bruyninx et al. 2013). However, with the increasing demand on velocities from reference frame (0.1 mm/year, Plag and Pearlman 2009) and sea level applications, we argue that even lower changes than 5% of GDP could be considered as significant. All the above shows is that periodic terms affect the velocity uncertainty much more at short timescales than they do for long-term data, leaving the values of velocity unbiased at the same time. With the increasing time span of observations, the assumption of seasonal signals becoming less important is validated. Here, the power-law character of the residuals plays a crucial role in determining the velocity uncertainty. In this way, 7, 9, and 17 years are enough for white, flicker, and random-walk noise, respectively, to decrease GDP below 5%, to omit periodic oscillations in the GNSS-derived time series and to take only noise model into consideration. So, providing a long time series, the periodicities we assumed will not influence the velocity uncertainty as much as noise would. This is why we should focus on obtaining the best estimate for the spectral index of each geodetic time series such as position, sea level, or zenith total delay.

Fig. 5
figure 5

Comparison between GDPs estimated for white noise (blue), flicker noise (red), fractional Brownian motion of spectral index equal to − 1.5 (orange) and random walk (black) plotted for deterministic model containing (1) linear velocity plus annual and semiannual oscillations (dashed lines) and (2) linear velocity plus extended model of periodicities from Bogusz and Klos (2016) (solid lines)

Full size image

Real GNSS time series

In order to confirm our theoretical approach with real data, we used position time series of continuous GPS stations. The time series we chose contributed to the newest International Terrestrial Reference Frame (ITRF2014). The GPS measurements were processed in the network solution named “repro2” by the International GNSS Service (IGS) (Rebischung et al. 2016, http://acc.igs.org/reprocess2.html). In this study, we picked 898 stations with time series of different lengths from 3 to 23 years that showed no specific or unusual behavior regarding long-term nonlinearity or earthquakes. This is to ensure that our results are not compromised by stations, which, for example, did not follow the assumption of a linear evolution in the coordinates. The station distribution is indicated in Fig. 6. The daily time series of the North, East, and Up components were pre-processed for outliers, offsets, and gaps if necessary and then analyzed with the reformulated maximum likelihood estimation (MLE) as implemented in the Hector software (Bos et al. 2013). We generated two different models to be fitted to GPS position time series. In the first case, the stochastic model was assumed to be a pure white noise with just linear velocity in the deterministic part. In the second case, a combination of power-law and white process with the full deterministic model was employed. To generate the GDP values, we estimated the ratio between (a) and (b), meaning that we fit (a) the appropriate noise model which has been already found to be preferred for GPS data and the most proper deterministic model taking all periodicities into account against (b) which is a simple linear fit with white noise assumption. The above exercise shows a clear combined effect of periodicities and power-law noise on the uncertainties of linear velocities.

Fig. 6
figure 6

Spectral indices of power-law noise estimated for selected ITRF2014 GPS stations for the North, East, and Up components

Full size image

The median amplitudes of the annual term are at the level of 1.65, 1.78, and 4.22 ± 0.20 mm for the North, East, and Up components, respectively. Almost all median amplitudes of the second, third, and fourth harmonics of the annual term fall below 1 mm with errors varying between 0.10 and 0.25 mm for horizontal and vertical changes.

Figure 6 presents a stochastic character: a power-law spectral index of examined stations. A clear spatial dependence may be noticed for few areas. In general, spectral indices estimated for ITRF2014 stations vary between − 1.6 and − 0.3, meaning we deal with different spectral characters of residuals. Spectral indices higher than − 0.5 were found for the shortest time series examined, confirming the findings of Williams et al. (2004) who emphasized, that white noise may cover a power-law character of residuals when time series are not long enough. For the North and East components, no spatial dependencies were found, expect for the area of Europe, for which the spectral indices we found are much more consistent than for any other part of the world. For the Up component, a clear spatial dependence may be found for Europe, indicating a strong impact of a Baltic Sea on the residuals of GPS position time series. That has been already also noticed by Klos and Bogusz (2017), explaining lower indices found for the coastal areas of the Baltic Sea by the impact of the sea. Also, stations situated within the areas of tundra and tropics are characterized by lower spectral indices than the rest of the stations. However, the detailed description of possible causes of such a character falls outside the scopes of this study.

Figure 7 presents GDP values for a globally distributed set of stations of series of different lengths from 3 to 23 years. From these figures, we can easily notice, that the GDPs for all components remain between 1 and 10 for a majority of stations with medians of 3.1, 3.2, and 2.8, respectively, for the North, East, and Up components. For 20, 17, and 19 North, East, and Up components, respectively, the GDP fall within 5%, as indicated by Blewitt and Lavallée (2002). For 5, 8, and 2 North, East, and Up components, respectively, the GDP values are equal to 1, meaning that there was no difference in velocity uncertainty when only white noise with a linear velocity was assumed and when a full deterministic model was employed along with the preferred noise model. The GDP estimates vary between 1 and 10 for a set of 756, 739 and 696 stations for the North, East, and Up components, respectively, of the ITRF2014 GPS position time series.

Fig. 7
figure 7

General Dilution of precision (GDP) for the North, East, and Up components. GDP means here a ratio between two assumptions of models being fitted into position time series: (1) a linear velocity plus a pure white noise, (2) a linear velocity, all seasonal terms and a combination of power-law and white noise model

Full size image

The results delivered for real GPS position time series are in good agreement with the theoretical formulae derived above. No spatial dependencies were found for the North and East components. The greatest deviations from 1 were observed for stations with large gaps within the data or significant variations in the amplitudes of seasonal signals, which were not modeled and therefore transferred to the stochastic part and also, for those, where the noise character was estimated to vary between flicker and random-walk noise. This means that the character of the noise, which affects the residuals of GPS data, influences the value of GDP for time series of sufficient length. For short series, the seasonal signals which were added to the deterministic model affect the GDP more than the character of the residuals, as was shown for the synthetic dataset. For a few stations, the GDP estimates are higher than 10. This was noticed for stations with large gaps inside the time span of observations. A clear influence of the Baltic Sea was noticed for the GDPs estimated for the Up component. All values are higher than 10, meaning that the uncertainties of estimates are highly affected by the type of noise by which time series are characterized.

Conclusions and discussion

Not accounting for the seasonal signals results in an increase in the autocorrelation or temporally correlated noise within the time series. This, in turn, influences the stochastic model when the periodic signals are not properly modeled. In this way, if we were certain about the presence of seasonal signals in the GNSS time series and did not model them, the residuals would resemble more a flicker noise. This would lead to increased uncertainties for all parameter estimates. Again, when the seasonal signals are properly modeled, another issue that may cause an artificial increase in the velocity uncertainty is an improperly assumed noise model itself. When flicker and random walk are being compared to a fractional Brownian motion, there may be an underestimation of the error bounds by a factor of two.

Some points can be easily noticed and raised for deeper discussion from the presented results. As long as the spectral index increases, the amplitudes of oscillations also increase. This arises from the fact that any power-law process with κ < 0 brings a correlation between amplitudes of seasonal terms and velocity. In this way, the GDP value is much higher for any time series length considered. The strong peaks of oscillations as seen in the GDP are indicated for short time scales, especially for the random-walk case. The applied oscillations play a significant role, even much more important than the a priori assumed noise character. The noise character starts to become important for time series longer than 9 years. The local minima and maxima of GDP are also being enlarged together with a change of the spectral index from 0 to − 2. This shows that the GDP may differ from integer-plus-half years as found by Blewitt and Lavallée (2002), who considered only white noise. This is clearly noticed in the case of random walk and has already been empirically confirmed by Bos et al. (2010). In this research, we provided mathematic formulas for the findings of Bos et al. (2010) and confirmed their correctness with synthetic and real GPS data and focused on power-law plus white noise, which is generally regarded as the best estimate for the stochastic model of GNSS time series. We showed that periodic signals are more important for short time scales, whereas the stochastic noise plays a significant role when the length of the time series increases. Also, with increasing spectral indices, the GDP decreases more slowly. We discussed a previously published approach, which indicated that 3.5 years of data are enough for the GDP to fall below 5%. When more seasonal signals and their harmonics were added to the deterministic model that include periodicities of all tropical and draconitic terms up to ninth harmonic plus fortnightly and the Chandlerian period, the GDP requires 9 years to fall below 5% for a white and flicker noise model. We have also discovered that the noise character starts to become more important than the periodic signals for time series longer than 9 years. And finally, Blewitt and Lavallée (2002) used the value of 5% to calculate the minimum velocity bias. However, this value is disputable. With the increasing demand on velocities, we argue that even smaller change in GDP could be considered as significant. This means that 7 and 9 years of continuous observations could be considered as the threshold for white and flicker noise, while 17 years for random walk to make the GDP to decrease to below 5% and to omit periodic signals in the GNSS-derived time series, taking only the noise model into consideration.

For real GPS data, the GDP values vary between 0 and 20, when a full deterministic model plus a proper noise model are considered against the assumption of linear velocity plus a white noise model. This clearly shows the combined effect which different assumptions for the mathematical model have on the uncertainties of the estimates. When the more realistic model is used, the closer to the truth and proper uncertainties are obtained.

Naturally, these assumptions are also valid for other kinds of data, such as VLBI, SLR, DORIS, or tide gauges. VLBI and SLR position time series are both characterized by the stochastic part being quite close to white noise (Feissel-Vernier et al. 2007; Botai et al. 2011), meaning that the type of noise would not have much impact on the uncertainties of the estimates. In these cases, a number of seasonals being fitted into position time series play a crucial role in the GDP estimates. As was shown by Collilieux et al. (2007) that the VLBI and SLR position time series are characterized by a number of significant oscillations, which have to be taken into account when a realistic model is to be employed. Riddell et al. (2017) found that the uncertainties of X, Y, and Z translations increase of 0.10 mm/year at a minimum, when a power-law noise is employed to describe them. The uncertainties of parameters delivered from DORIS series with the preferred noise model described as white plus flicker noise (Williams and Willis 2006; Khelifa et al. 2012) will be influenced much more by the character of noise than VLBI and SLR position series.