The Noise Properties and Velocities from a Time-Series of Estonian Permanent GNSS Stations

Abstract

The aim of this study was to estimate the noise properties, velocities, and their uncertainties from a time-series of selected (~9 years long) Estonian continuously operating Global Navigation Satellite System (GNSS) stations. Two software packages based on different processing methods, Gipsy–Oasis and Bernese, were used for daily coordinate calculations. Different methods and software (Tsview, Hector, and MIDAS) were used for coordinate time-series analysis. Outliers were removed using three different strategies. Six different stochastic noise models were used for trend estimation altogether with the analysis of the noise properties of the residual time-series with Hector. Obtained velocities were compared with different land uplift and glacial isostatic adjustment models (e.g., ICE-6G (VM5a), NKG2016LU, etc.). All compared solutions showed similar fit to the compared models. It was confirmed that the best fit to the time-series residuals were with the flicker noise plus white noise model (for the North and East component) and generalized Gauss–Markov model (for Up). Velocities from MIDAS, Tsview, and Hector solutions within the same time-series (Gipsy–Oasis or Bernese) agreed well but velocity uncertainties differed up to four times. The smallest uncertainties were obtained from Tsview; the MIDAS solution produced the most conservative values. Although the East and Up component velocities between Gipsy and Bernese solutions agreed well, the North component velocities were systematically shifted.

1. Introduction

Glacial isostatic adjustment (GIA) is a dominant cause of the long-term intra-continental deformations and motion of the geodetic points in Estonia. Glacial isostatic adjustment-related land uplift signal in precise levelling data and tide gauge observations are quite well studied in Estonia [1,2,3,4,5]. However, no detailed analysis about the Estonian continuously operating or permanent GNSS (cGNSS) stations, their velocities, uncertainties, and noise properties have been published yet, except for a brief paper in Estonian [6]. On the contrary, GIA-related studies using permanent cGNSS stations time-series have a quite long history in Fennoscandia [7,8,9]. Recently, the semi-empirical post-glacial land uplift (LU) model NKG2016LU (Figure 1) was released by the Nordic Geodetic Commission (NKG) [10]. It combines the geodetic observations (repeated levelings and velocities of cGNSS stations) with predicted land motion from the GIA model. Also, a joint NKG GNSS Analysis Centre has been launched by the initiative of NKG. It aims to routinely produce high-quality GNSS solutions (coordinates, velocities, and their uncertainties) for the common needs of the NKG and the Nordic and Baltic countries [11]. Like GNSS, other satellite-based geodetic methods (SLR, DORIS, InSAR) have also proved to provide accurate station position [12,13,14], so their time-series have also been used for various purposes. Some of the examples include determination of plate tectonics, reference frames, geocenter motion, and geohazards [15,16,17,18].

In order to reliably estimate cGNSS stations’ velocities and their uncertainties for geophysical interpretation, the noise in coordinate values needs to be modelled [19,20]. Generally, the noise in GNSS time-series has been described as the combination of colored noise and white noise [21]. It was found that the combination of flicker noise plus white noise models approximates well the noise properties of all three coordinate components [22,23,24,25]. However, other models like power law noise plus white noise or flicker noise plus random walk noise also provide precise solutions for the noise properties of the data [20,21,26]. Langbein [26] and Goudarzi et al. [25] have shown that white noise is associated with the GNSS hardware noise and measurement errors, flicker noise with the GNSS system (i.e., unmodeled tropospheric delay and satellite ephemeris), and random walk noise is generally associated with the random motions of the station’s monument [25,26].

Various software packages for high-precision coordinate calculations (Bernese, Gipsy–Oasis, GipsyX, Gamit/Globk) and GNSS time-series analysis (CATS, Hector, est_noise, MIDAS, Tsview, etc.) are available [27,28,29,30,31,32,33,34]. Some papers that compare Gipsy–Oasis, Bernese, and Gamit/Globk solutions have shown that coordinate differences are in the 2 to 3 mm level in general [35,36,37,38]. Bos et al. [32] and He et al. [39] compared Hector with CATS and found that there was no significant difference between them. Langbein [40] showed that the Hector and est_noise velocity estimates differed less than 0.1 mm/yr.

The aim of this study was to estimate the GNSS velocities, their uncertainties, and noise properties of the coordinate time-series of the selected cGNSS stations in Estonia (Figure 1) by using software packages utilizing different methods (Gipsy–Oasis and Bernese for coordinate calculations, Hector, Tsview, and MIDAS for velocity estimation), and a selection of stochastic noise models. The results were validated using LU and GIA models, including NKG2005LU [41], NKG2016LU, and NKG2016GIA [10], D1 [42], and ICE-6G(VM5a) [43].

The main motivation to conduct this study is related to the ongoing research project PUT1553 of the Estonian Research Council. In that project, the geocentric land uplift together with the sea level rise along the Estonian coast will be estimated in order to establish risks on coastal areas resulting from the sea level rise. For that, the stability of the cGNSS stations (based on noise properties of time-series) was also analyzed in more detail, i.e., how suitable they are for crustal deformation studies.

The paper is organized as follows: in Section 2 an overview of Estonian cGNSS stations is presented; programs and strategies used in coordinate calculations, time-series analysis, and noise properties study are introduced. In Section 3, results of the coordinate calculations, time-series, and noise analysis are presented and discussed. Outcomes from different computation strategies and programs are compared with each other and with velocities from LU and GIA models. Section 4 summarizes findings.

2. Materials and Methods

2.1. Calculation of cGNSS Stations Coordinates

The state network of cGNSS stations in Estonia, named ESTPOS (Figure 2), has gradually evolved since the mid-2000s (the first station was established in 1996). Today it forms a dense network of 29 stations. Most of them are mounted on buildings at least 30 years old, so the long-term stability of these stations is expected. All stations are equipped with the Leica GNSS receivers GR25 and AR25 GNSS choke ring antennas with radome LEIS. Four stations (KURE, SUR4, TOIL, TOR2) also belong to the EUREF Permanent GNSS Network (EPN) [44].

Unfortunately, the length of observations of most of these 29 stations was too short (less than three years) to calculate reliable velocities [45]. Therefore, 10 cGNSS stations (Figure 1) with the lowest number of jumps and irregularities in time-series and with nine-year-long time-series (2008–2016) were selected for the current study. Two of these stations (KARG and MISS) belong to the private network Trimble VRS Now. Stations of Trimble VRS Now are equipped with the Trimble NetR5 and NetR9 receivers together with the TRM55971.00 antennas.

Two software packages, Gipsy–Oasis II ver. 6.4 (GOA) and Bernese ver. 5.2 (BERN) with different approaches for coordinate calculations were used to derive daily coordinates of the cGNSS stations (Figure 1). The GOA is based on the precise point positioning method (PPP, cf. [27]) whilst BERN uses the classical double difference (DD) approach where a network of cGNSS stations is computed in a single adjustment [33]. The dense network of 56 cGNSS stations in Estonia (ESTPOS stations plus stations of private cGNSS networks) used in BERN calculations are presented in Figure 2.

The precise orbit and clock parameters in ITRF2008 reference frame by Center for Orbit Determination in Europe (CODE) [46] and Jet Propulsion Laboratory (JPL) were used in BERN and GOA, respectively. Vienna Mapping Function (VMF1) was used as an a priori troposphere model to calculate the hydrostatic and wet tropospheric delay. The default GOA and BERN settings were used to estimate troposphere parameters (zenith delay, troposphere gradients). In addition, site specific troposphere parameters were calculated to correct the a priori model. The igs08.atx and epnc_08.atx antenna calibration files were employed to model the phase center variations and offsets in GOA and BERN, respectively. Ocean loading coefficients corrected for solid Earth’s mass center variations were obtained from the FES2004 model [47]. Observations were corrected for second-order ionosphere correction based on the International Reference Ionosphere model IRI-2012b and shell height 600 km. Observations with an elevation angle lower than 10° were not used. The loading effects due to the atmospheric and hydrological mass variations were not modelled in this study. The same observed GNSS data (RINEX files) were processed separately by two institutions: the Estonian Land Board with BERN and the Estonian University of Life Sciences with GOA. Processing parameters are summarized in

Table 1. Parameters used in coordinate calculation by Gipsy-Oasis II 6.4 (GOA) and Bernese 5.2 (BERN) software.

PARAMETER	BERN	GOA
Satellite system	GPS	GPS
Method	double difference (DD)	precise point positioning (PPP)
Elevation mask	10°	10°
Precise orbits and clocks	CODE final products	JPL final products
Antenna calibration file	epnc_08.atx	igs08.atx
Troposphere mapping function	VMF1	VMF1
Ocean loading tide model	FES2004	FES2004
Reference frame	ITRF2008	ITRF2008
Second order ionosphere correction	Yes	Yes

. Further time-series analysis (detection of outliers, noise properties, velocities, and their uncertainties) was carried out by the Estonian University of Life Sciences.

2.2. Time-Series Analysis

Daily coordinates in ITRF2008 reference frame obtained from GOA and BERN solutions were analyzed using the model function [21,48]:

$$y\left(t\right)=y_0+bt+{\displaystyle \sum _{j=1}^2\left\{A_j\sin {\omega }_{j}t+B_j\cos {\omega }_{j}t\right\}}+{\displaystyle \sum _{j=1}^n{C}_{j}H\left(t-T_j\right)+{\epsilon }_y\left(t\right)},$$

(1)

where $y\left(t\right)$ is the position at the epoch $t$ (in years), $y_0$ and $b$ are the position ($t=0$) and velocity parameters, ${\omega }_j$ $\left(j=1,2\right)$ is the angular frequency of an annual and semiannual seasonal term, respectively, $A_j$ and $B_j$ their respective amplitudes of the sine and cosine parts and ${\epsilon }_y$ represents the noise. The offset terms ${\displaystyle \sum _{j=1}^n{C}_{j}H\left(t-T_j\right)}$ are used to model the equipment change (antenna or receiver replacement), repositioning of the antenna or station, earthquakes, and other rapid events in which the $T_j$ is the time of the offset, $n$ is the number of offsets, and coefficient $C_j$ is the magnitude of the offset. The Heaviside step function $H\left(t\right)$ equals 1 for $t\ge 0$ and 0 otherwise. The unknown values of parameters $y_0$, $b$, $A_j$, $B_j$ and $C_j$ are found by fitting the Model (1) to the position time-series of every cGNSS stations.

Firstly, the Matlab tool for GNSS time-series analysis Tsview [28] based on least squares fitting was used. Strengths of this program are: (i) graphical interface, (ii) the possibility to use the database with offsets (antenna and place changes, earthquakes, etc.), (iii) removing outliers with different options (points with large coordinate uncertainties, outliers based on time-series residuals), and (iv) the possibility to perform spectral analysis. One shortcoming, the possibility to use only one correlated noise model (called RealSigma) for velocity uncertainty estimation, can be mentioned. The exact dates of antenna position changes were described for offset calculation. Points with the coordinate uncertainties greater than ±10 mm and outliers detected in time-series by using Tsview’s 3-sigma criterion were removed. The time-series were detrended using annual and semiannual periodic signals. The velocity uncertainties were estimated using RealSigma option.

Secondly, software Hector [32] based on the maximum likelihood estimation (MLE) was employed to allow simultaneous estimation of noise structure together with the velocities and their uncertainties. Hector accepts only stationary noise (with constant noise properties) for matrix operations. Therefore, it operates faster than other similar programs used for time-series analysis, e.g., CATS [30]. Hector allows users to estimate annual and semiannual signals, performing spectral analysis, and to calculate velocities and their uncertainties using white noise and selected correlated noise models as well as a combination of different noise models. Hector also allows the possibility to remove outliers from a time-series based on the interquartile range (IQR) of the residuals from the least squares linear fit. Residuals outside the chosen threshold (IQ-factor multiplied by IQR) are considered as outliers. The recommended value for the IQ-factor is three [49,50].

Two IQ-factors, 3 and 2.2, were tested for removing outliers. The IQ-factor 2.2 was used to obtain comparable results with the Tsview’s 3-sigma level results which is based on the relation: sigma = IQR/1.349 [51]. After removing the outliers, the time-series were detrended using annual and semiannual seasonal signals and analyzed with five different time-correlated noise models (see Section 2.3). In addition, trend estimation, from which outliers were removed with Tsview 3-sigma criterion, was also performed.

The last tested program was MIDAS, which uses the Theil–Sen median trend estimator [34]. According to their study, this approach makes the MIDAS solution robust to effects due to the seasonal signals, outliers, offsets, and heteroscedasticity in GNSS time-series. The step file with the antenna or place changes was created and the time-series were processed using MIDAS’ default options.

2.3. Noise Analysis

Understanding the noise properties of GNSS time-series is crucial for realistic uncertainty estimation of unknown parameters. Several studies have shown that the assumption of the time-uncorrelated white noise leads to underestimation of the velocity uncertainties [20,21,23,25]. Accordingly, it is best to use a combination of the time-correlated colored noise plus white noise to describe noise in GNSS time-series.

The common approach for describing time-correlated colored noise in GNSS time-series is a power law process with power spectrum function:

$$P_x\left(f\right)=P_0\cdot f^k,$$

(2)

where $P_0$ is a normalizing constant, $f$ is temporal frequency, and $k$ is the spectral index to describe the frequency dependence of the process. At $k=0$ the process is a white noise (WN), at $k=-1$ the process is a flicker noise (FN), and at $k=-2$ the process is a random walk (RW). When $k$ is an unknown parameter to be estimated, the process is a general power law noise (PLN) [20].

Another process, called a Gauss–Markov, follows the power spectrum function:

$$P_x\left(f\right)=P_0{\left({\beta }^2+4{\pi }^2{f}^2\right)}^{\frac{k}{2}},$$

(3)

where $\beta$ is the crossover frequency symbolizing the point where the low- and high-spectrum frequencies cross each other. This process is frequency independent for lower frequencies ($k=0$) and approximates the power law process for the higher frequencies. When $k$ is the unknown parameter to be estimated, the general Gauss–Markov (GGM) process is assumed. The special case of the Gauss–Markov noise model is the first order Gauss–Markov (GM), where $k=-2$ [20], which is also an assumption in Tsview’s RealSigma estimation.

For estimation of realistic uncertainties, time-correlated noise properties of GOA and BERN time-series were analyzed using power spectral density (PSD) plots and different noise models. Analysis was made using the program Hector and applied noise models were:

• white noise (WN);
• generalized Gauss–Markov noise (GGM);
• generalized Gauss–Markov noise + white noise (GGM + WN);
• power law noise + white noise (PLN + WN);
• flicker noise + white noise (FN + WN);
• flicker noise + random walk (FN + RW).

As was mentioned above, Hector can only deal with a stationary noise to allow fast inversion techniques. Due to that limitation, flicker noise and random walk noise models can be employed only by using the GGM noise model with a specific parameter $1-\varphi$ [49] or $1-\alpha \delta T$ in terms of Langbein [29], where $\delta T$ is the measurement’s sampling interval and $\alpha$ is the inverse of the time constant of the process. In the present study, the value for the $1-\varphi$ parameter was estimated first from the GGM noise model, and then used as input for noise models involving flicker noise and/or random walk noise.

The relative goodness of the models can be estimated by the Akaike and Bayesian information criterions (AIC and BIC), respectively [49,52,53]. The model with the lowest AIC/BIC value indicates the best fit with the observations [21]. However, sometimes the estimation by the AIC and BIC values favors different models. In that case, a decision about the best noise model can be based on the inspection of the PSD and residual plots.

3. Results and Discussion

3.1. Outliers Removal

The presence of outliers among GOA and BERN time-series was tested using different software and criteria. The tested solutions were named as follows:

• GOA, outliers removed with Tsview and 3-sigma level (solution GOA_1);
• GOA, outliers removed with Hector and IQ-factor 2.2 (solution GOA_2);
• GOA, outliers removed with Hector and IQ-factor 3 (solution GOA_3);
• BERN, outliers removed with Tsview and 3-sigma level (solution BERN_1);
• BERN, outliers removed with Hector and IQ-factor 2.2 (solution BERN_2);
• BERN, outliers removed with Hector and IQ-factor 3 (solution BERN_3).

Results are presented in

Table 2. Amount on the detected outliers (North, East, and Up (NEU) component together) in GOA and BERN time-series using different criteria.

Station Name	GOA Time-Series				BERN Time-Series
	Initial no of Points	Tsview 3-Sigma	Hector IQ = 2.2	Hector IQ = 3	Initial no of Points	Tsview 3-Sigma	Hector IQ = 2.2	Hector IQ = 3
		(GOA_1)	(GOA_2)	(GOA_3)		(BERN_1)	(BERN_2)	(BERN_3)
AUDR	3144	5.6%	3.2%	0.9%	3071	19.5%	8.2%	2.9%
KARD	2822	10.0%	6.0%	2.4%	2765	20.6%	13.1%	7.5%
KARG	2562	7.2%	6.0%	3.4%	2499	9.1%	5.2%	2.5%
KURE	3277	4.5%	2.6%	0.4%	3205	9.0%	5.3%	2.0%
MISS	2626	4.3%	3.2%	0.6%	2539	4.5%	2.6%	0.8%
MVEE	2874	9.8%	6.7%	3.0%	2838	15.8%	12.0%	5.5%
SUR4	3273	4.0%	2.2%	0.3%	3188	10.2%	6.4%	2.9%
TOIL	3272	5.3%	2.7%	0.4%	3168	7.5%	4.2%	1.0%
TOR2	3204	3.9%	2.2%	0.4%	3057	13.4%	4.7%	0.8%
VOR2	2795	22.5%	14.7%	6.4%	2447	29.3%	21.9%	10.9%
Average		7.7%	4.9%	1.8%		13.9%	8.4%	3.7%

.

The Tsview three-sigma rule detected the largest number of outliers from GOA and BERN time-series (Table 2). Although the Hector’s IQ-factor 2.2 should theoretically produce compatible results with the Tsview three-sigma rule, the number of the outliers detected by the former method is almost two times smaller. Hector’s IQ-factor three removes up to three times fewer outliers compared to the IQ-factor 2.2. Note that the number of outliers in the BERN time-series was approximately two times larger, compared to GOA time-series, no matter what software or rejection criterion was used. The reason could be the noise propagation and correlations within dense cGNSS stations network (Figure 2) of the DD solution, which included cGNSS stations with shorter time-series, different monumentations, etc., compared to the PPP solution (Figure 1).

3.2. Noise Properties

Noise properties of six different time-series obtained from different outlier removing methods (see Section 3.1) were analyzed using six different noise models mentioned in Section 2.3. The relative goodness of the time-series was determined using AIC and BIC values provided by Hector in order to select the best time-series within GOA and BERN solutions for further analysis. It turned out that AIC/BIC values of the solution GOA_1 were ~2% and 5% smaller on average, compared to the solutions GOA_2 and GOA_3, respectively. The AIC/BIC values of the solution GOA_2 were ~3% smaller compared to the solution GOA_3.

The differences of AIC/BIC values were larger for the BERN time-series: the values of the solution BERN_1 were ~4% and ~11% smaller on average, compared to the solutions BERN_2 and BERN_3, respectively. The AIC/BIC values of the solution BERN_2 were ~6% smaller compared to the solution BERN_3. Further analysis of the noise properties was conducted with the time-series where outliers were removed using software Hector and the IQ-factor 2.2 (solutions GOA_2 and BERN_2). This was a compromise between the AIC/BIC values and a percentage of the removed outliers (Table 2). Apparently, removing many data points flagged as outliers improves the fitting results (and lowers the AIC/BIC values), but the real signals searched by trend modelling could then be lost.

The PSD plots of the time-series residuals are presented in Figure 3 and Figure 4. The comparison of the plots shows that the North, East, and Up component of the BERN time-series tend to follow power law noise (linear slope in logarithmic scale) in general. However, more white noise (flat spectrum) seems to appear in the BERN time-series, compared to GOA. It can also be seen, that the GOA’s Up component noise was flatter on low frequencies than the BERN’s Up component. The PSD plots of all cGNSS stations are presented in Table S1 of the Supplementary Materials.

From spectral analysis of residual time-series, the strongest power was found in the one-year period. In the North component of GOA_2 solution, power peaked at roughly 26–27 cpy, corresponding to the period of the ~14-day cycle. The 14-day cycle was observed in cGNSS time-series previously [54]. It approximately follows the aliased period resulting from the unmodeled diurnal O1 and semidiurnal M2 tides [55].

For detailed analysis of the noise properties, separate PSD plots of the residual time-series for each station were estimated and fitted with theoretical noise models. Examples of such plots are presented in Figure 5 and Figure 6 for stations TOIL and MISS. Accordingly, the best suited noise models for GOA time-series are GGM or GGM + WN (since the fraction of the WN is very small) for the North component, GGM, GGM + WN, or PLN + WN for the East component, and GGM or FN + RW for the Up component. For the BERN time-series, the best suited noise models were GGM for the North component, GGM, GGM + WN, or PLN + WN for the East component, and GGM or GGM + WN for the Up component.

Since the visual selection of the best fitting noise model can be ambiguous, relative goodness of the noise models were estimated based on the AIC and BIC values provided by Hector

Table 3. Number of times that each noise model had the lowest Akaike or Bayesian information criterion (AIC or BIC) values for the North, East, and Up component based on the GOA time-series solution GOA_2.

Noise Model	North			East			Up			Total SUM	% (NEU Component)
	AIC	BIC	SUM	AIC	BIC	SUM	AIC	BIC	SUM		AIC	BIC
FN + RW	0	0	0	1	1	2	5	7	12	14	20.0%	26.7%
FN + WN	3	8	11	8	9	17	0	0	0	28	36.7%	56.7%
GGM	7	2	9	1	0	1	5	3	8	18	43.3%	16.7%
GGM + WN	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%
PLN + WN	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%
WN	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%

FN + RW—combination of flicker noise and random walk, FN + WN—combination of flicker noise and white noise, GGM—generalized Gauss-Markov noise, GGM + WN—combination of generalized Gauss-Markov noise and white noise, PLN + WN—combination of power law noise and white noise, WN—white noise.

and

Table 4. Number of times that each noise model had the lowest AIC or BIC values for the North, East, and Up component based on the BERN time-series solution BERN_2.

Noise Model	North			East			Up			Total SUM	% (NEU Component)
	AIC	BIC	SUM	AIC	BIC	SUM	AIC	BIC	SUM		AIC	BIC
FN + RW	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%
FN + WN	4	9	13	5	9	14	1	5	6	33	33.3%	76.7%
GGM	0	0	0	3	1	4	7	5	12	16	33.3%	20.0%
GGM + WN	6	1	7	2	0	2	2	0	2	11	33.3%	3.3%
PLN + WN	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%
WN	0	0	0	0	0	0	0	0	0	0	0.0%	0.0%

show how many times the corresponding noise model had the lowest AIC or BIC value for each component, and total supremacy for North, East, and Up (NEU) components in percent (sum of the NEU’s AIC or BIC divided by number of coordinate components × number of stations, i.e., 30, for more details see Reference [21]). From these we can see that the AIC and BIC values were not consistent with each other, especially for the North component. There is some dispute on model selection when AIC and BIC results are controversial [56,57]. Therefore, an attempt was made to select the best noise model by looking at both criteria together. The sum of the numbers shows the number of times the corresponding noise model had the lowest AIC or BIC value (Table 3 and Table 4), indicating that the best noise models for the GOA time-series are FN + WN for North and East, and FN + RW or GGM for the Up component. However, He et al. [21] has shown that the noise models which include RW can provide a better fit with the long time-series. The overall best noise model for the GOA_2 solution was FN + WN following by GGM, based on the sum of the AIC and BIC values of all components.

The same analysis for the BERN time-series showed that the best noise models were FN + WN for the North and East components, and GGM for the Up component. The overall best noise model was FN + WN followed by GGM, also for the GOA time-series. This confirms previous findings from numerous studies [20,21,22,23] that the FN + WN is the most suitable noise model to describe stochastic processes in long GNSS time-series.

As the four used noise models were combinations of white and colored noise, the fraction of WN as well as RW noise in combined models were estimated also. The fraction of the WN was quite small in the GOA time-series, ranging from 0.4 to 7.6%. On the other hand, the percentage of the WN was quite high in the BERN time-series, ranging from 20 to 42%. The fraction of the RW noise was almost zero in the BERN time-series. However, the GOA time-series showed a very large proportion of the RW noise in the Up component (90.4%). The reasons for these controversial findings concerning the RW noise are currently unknown. One possible explanation could be that some station-dependent effects common to all stations in our Bernese adjustment were reduced during DD processing, but not in PPP.

An important characteristic of the noise properties is spectral index $k$ (2). The spectral indexes were estimated from PLN + WN and GGM noise models during the velocity estimation process in Hector. There was a good agreement between GOA_2 and BERN_2 as well as between PLN + WN and GGM solutions concerning North and East component spectral indexes. The average spectral index for North and East components were −0.8 and −0.9, respectively, which are close to the spectral index of a flicker noise. The Up component spectral indexes of BERN time-series agreed between the PLN + WN and GGM noise models, being −0.7 on average. However, spectral indexes of the Up component of the GOA time-series disagreed between the PLN + WN and GGM solutions. Since the GOA time-series’ PSD of the Up component was quite flat on low frequencies, the PLN model did not fit with the observations and the estimated spectral index was equal to the spectral index of the FN (−1.0). The spectral index estimated from the GGM model was almost equal to the spectral index of the RW noise (−2.0). As we saw from Table 3, FN + RW noise was the preferred noise model and RW was the dominant noise for the Up component of the GOA time-series. Calculated spectral indexes are presented in Table S2 of the Supplementary Materials.

3.3. Time-Series Analysis

All time-series (GOA and BERN solutions) were visually inspected using Tsview. The visual inspection of the residual time-series of some stations revealed large residual periodic signals. As a pre-eminent example, the residual time-series of station TOR2 is discussed here. Peak-to-peak amplitude of the residual signal of this station was up to 10 mm, especially in the East component. It is interesting that the phases of the residual signal were changed ~90° after the applied break (due to the station re-installation to a new building a few hundred meters away from the old site), even though amplitudes of the residual signal have remained the same. Several reasons for such apparent seasonal residual signal, like mismodelling of the propagating medium effects, bedrock (or building) thermal expansion, antenna phase center modelling error and local multipath, can be found in the literature [58].

The amplitudes and phases of annual and semiannual periodic signals estimated for all cGNSS stations (except for TOR2) from the GOA and BERN time-series can be found in Tables S3 and S4 of the Supplementary Materials. The GOA and BERN time-series (N, E components) in TOR2 present a clear phase shift of the annual signal in 2011 for an unknown reason. Because of such a phase shift, only the linear trend was used to evaluate the velocity components of TOR2. The time-series plots of all cGNSS station are presented in Table S5 of the Supplementary Materials.

After visual inspection of the time-series, trend and uncertainties were estimated (using RealSigma option) with the Tsview. Results are presented in Table S6 of the Supplementary Materials. Next, trend and uncertainties were estimated with Hector using six noise models (listed in Section 2.3) and with the MIDAS. Results are discussed in the next section.

3.4. Comparison and Validation of the Velocity Estimates

The velocities and their uncertainties by using different noise models and outlier removal strategies were compared to inspect the consistency among results. In addition, all velocity solutions were compared with the velocities based on the LU and GIA models: NGK2005LU [41], NKG2016LU and NKG2016GIA [10], ICE-6G(VM5a) [43], and D1 [42], and with the GNSS velocity solution for five stations (AUDR, KURE, SUR4, TOIL, TOR2) used in NKG2016LU (hereafter referred as BIFROST2016 [10]). For the models NGK2005LU, NKG2016LU, and D1 only Up component velocities were available for comparison. The NKG2016GIA, ICE-6G(VM5a), and BIFROST2016 include horizontal velocities also. The velocities of cGNSS stations interpolated from the abovementioned models are presented in Tables S7–S12 of the Supplementary Materials. Due to the different reference frame issues, only Up component velocities were compared with the LU and GIA models (reference frame for NKG2005LU is ITRF2005, for NKG2016LU and BIFROST2016 is ITRF2008. The latter coincides with the reference frame used in our GNSS solution. The GIA model uses its own geocentric (Earth’s center of mass) reference frame). Considering the small size of Estonian territory, constant vertical shift with the cGNSS solution due to the difference of reference frames could be expected. The statistics of the cGNSS stations’ Up component velocity residuals from the models are presented in Section 3.4.1, Section 3.4.2, Section 3.4.3, Section 3.4.4.

3.4.1. Gipsy–Oasis Time-Series and Hector Solution

Velocities and their uncertainties from different solutions (GOA_1, GOA_2, and GOA_3) showed only small differences. The difference between the solutions was within 0.02 ± 0.07 (the mean and standard deviation) and 0.04 ± 0.03 mm/yr for velocities and uncertainties, respectively. This indicates that different strategies for outlier removal had insignificant impact on velocity estimates. This could be related to the length of time-series used in the current study. Klos et al. [59] found that different criterions of outlier removal lead to different results; however, their time-series were shorter.

Larger influence on the estimated velocities and their uncertainties was found by using different noise models. As demonstrated in the previous section, the models describing the stochastic processes of selected GNSS time-series should be either FN + WN, FN + RW or GGM (Table 3). The differences of velocities between FN + WN, FN + RW, GGM, and PLN + WN solutions are presented in Figure 7. There are only slight differences between the FN + WN and FN + RW as well as of FN + WN and GGM solutions. This applies also to the differences in the North and East components obtained by FN + WN and PLN + WN models. Differences between FN + WN and PLN + WN results are larger for the Up component, especially for some stations. This is related to the fact that the power spectrum of the Up component of the GOA time-series was relatively flat on low frequencies, and therefore, the PLN model is less suitable to describe its noise (Figure 5).

The differences between the velocity uncertainties of selected noise models of the GOA_2 solution are presented in Figure 8. Differences are relatively small between the pairs of FN + WN and FN + RW models, and FN + WN and GGM models. However, differences were larger between the FN + WN and PLN + WN model Up component. The reason for the uncertainty differences is the same as for the velocity differences: PLN models do not fit with the observations for the Up component of the GOA time-series (Figure 3).

Velocity uncertainties compared to the WN model were 6.2 (from 1.5 to 15 depending on the noise model) times larger on average using colored noise models. The largest differences were with the PLN + WN model, which is expected, since the PLN model did not fit with the observations (Figure 5).

All compared solutions showed also similar differences from LU and GIA models. The Hector solution with the IQ-factor 2.2 (solution GOA_2) demonstrated slightly better overall performance. As this solution revealed also small AIC and BIC values compared to the other solutions, it was selected for further investigation to select the best performing noise model within this solution. Results of the comparison with the different LU and GIA models are shown in

Table 5. The mean and standard deviation in mm/yr of the Up velocity residuals from different LU and GIA models based on the solution GOA_2.

Noise Model	NKG2005-LU [41]	NKG2016-LU [10]	NKG2016-GIA [10]	D1 [42]	ICE-6G-(VM5a) [43]	BIFROST-2016 [10]
WN	0.08 ± 0.36	0.16 ± 0.40	0.21 ± 0.63	−0.30 ± 0.46	−1.27 ± 0.86	0.27 ± 0.22
PLN + WN	0.08 ± 0.49	0.17 ± 0.56	0.21 ± 0.76	−0.30 ± 0.58	−1.26 ± 0.85	0.16 ± 0.27
GGM	0.08 ± 0.38	0.17 ± 0.43	0.21 ± 0.66	−0.30 ± 0.48	−1.26 ± 0.85	0.26 ± 0.23
GGM + WN	0.08 ± 0.38	0.17 ± 0.43	0.21 ± 0.65	−0.30 ± 0.48	−1.26 ± 0.85	0.26 ± 0.23
FN + WN	0.08 ± 0.37	0.16 ± 0.42	0.21 ± 0.64	−0.30 ± 0.47	−1.26 ± 0.85	0.27 ± 0.22
FN + RW	0.08 ± 0.38	0.16 ± 0.43	0.21 ± 0.66	−0.30 ± 0.48	−1.26 ± 0.85	0.26 ± 0.23

. Systematic shift (mean of the residuals) can be seen with almost all compared models, which is mostly related to the different reference frames between cGNSS solution and tested models, but also with the used noise models. For example, smaller systematic shift with BIFROST2016 solution was obtained using PLN + WN solution. This is presumably related to the fact that the same noise model was also used in the BIFROST2016 solution.

Based on the AIC and BIC values (Table 3), as well as on comparison with the LU and GIA models (Table 5), the most accurate velocity solution was obtained using the FN + WN noise model. The best noise model, however, was different for the N, E, and U components. For the N and E components, it was the FN + WN, and for the U component, it was the FN + RW or GGM noise model, respectively. The velocities and uncertainties from solutions FN + WN (N, E) and GGM (U) are given in Table S6 of the Supplementary Materials.

3.4.2. Bernese Time-Series and Hector Solution

The same comparison was made with the BERN time-series. Similar to the GOA case, velocities and their uncertainties from different solutions (BERN_1, BERN_2, and BERN_3) showed only small differences. The difference between solutions was within 0.01 ± 0.03 and 0.01 ± 0.01 mm/yr for velocities and uncertainties, respectively, which shows that the different strategies for outliers’ removal had practically no impact on the BERN time-series velocity or uncertainty estimates.

Again, larger influence on the estimated velocities and their uncertainties was noticed from the used noise model. For the BERN time-series, the models describing the stochastic processes of the selected GNSS time-series should be either FN + WN, GGM or GGM + WN (Table 4). The differences of velocities between FN + WN, FN + RW, GGM, and PLN + WN noise models are presented in Figure 9. There are only slight differences between the pairs of selected noise models. This is related to the fact that contrary to the GOA time-series, the power spectrum of the Up components of the BERN time-series are not as flat on low frequencies, and therefore, the PSD of PLN and FN models are more similar (Figure 6).

The same applies to the differences between the velocity uncertainties of the selected noise models (Figure 10) Differences are relatively small between FN + WN and FN + RW, FN + WN and GGM, and FN + WN and PLN + WN models. Differences are larger for only some stations, especially in the Up component (stations MVEE, TOR2, VOR2, cf. Figure 1). The anomalies in time-series of these stations can be observed also (Table S5 of the Supplementary Materials). The reason for the small uncertainty differences is the same as for the velocity differences: PSD of PLN, FN and GGM noise models of the BERN time-series are closer to each other, compared to the GOA time-series (Figure 6).

Velocity uncertainties compared to the WN model were 6.5 (from 3.2 to 10.9 depending on the noise model) times larger on average using colored noise models. The largest differences were found with the FN + RW model, which was expected, since this noise model fit less with the observations (Figure 6).

The results of comparing velocities with the LU and GIA models (

Table 6. The mean and standard deviation in mm/yr of the Up velocity residuals from different LU and GIA models based on the solution BERN_2.

Noise Model	NKG2005-LU [41]	NKG2016-LU [10]	NKG2016-GIA [10]	D1 [42]	ICE-6G-(VM5a) [43]	BIFROST-2016 [10]
WN	0.17 ± 0.40	0.26 ± 0.45	0.30 ± 0.67	−0.21 ± 0.51	−1.17 ± 0.86	0.28 ± 0.20
PLN + WN	0.10 ± 0.41	0.19 ± 0.47	0.24 ± 0.69	−0.27 ± 0.52	−1.24 ± 0.85	0.24 ± 0.17
GGM	0.14 ± 0.42	0.23 ± 0.47	0.27 ± 0.69	−0.23 ± 0.53	−1.20 ± 0.86	0.26 ± 0.19
GGM + WN	0.15 ± 0.42	0.23 ± 0.48	0.28 ± 0.70	−0.23 ± 0.53	−1.19 ± 0.85	0.27 ± 0.19
FN + WN	0.14 ± 0.42	0.22 ± 0.48	0.27 ± 0.69	−0.24 ± 0.53	−1.20 ± 0.86	0.26 ± 0.18
FN + RW	0.14 ± 0.42	0.22 ± 0.48	0.27 ± 0.70	−0.24 ± 0.53	−1.20 ± 0.85	0.26 ± 0.18

) were quite similar to the GOA results. Although the solution with outliers removed using Tsview with three-sigma level, solution BERN_1 showed better performance based on the AIC and BIC values, comparison with the LU and GIA models revealed quite similar residuals for all solutions. Therefore, the Hector solution with the IQ-factor 2.2 (solution BERN_2) was chosen for further selection of the best performing noise model within this solution.

Based on the AIC and BIC values (Table 4) as well as on comparison with the LU and GIA models (Table 6), it was then decided that the overall best velocity solution was obtained using the FN + WN noise model. Again, as with the GOA case, the best noise model was slightly different for the N, E, and U components, based on the AIC and BIC values. For the N and E components, it was the FN + WN, and for the U component the GGM noise model, respectively. The velocities and uncertainties from solutions FN + WN (N, E) and GGM (U) are given in Table S6 of the Supplementary Materials.

3.4.3. MIDAS Solution

Similar comparison with the LU and GIA models was made also with the velocities obtained from the MIDAS solution (

Table 7. The mean and standard deviation in mm/yr of the Up velocity residuals from different LU and GIA models based on the MIDAS and Tsview solutions.

Noise Model	NKG2005-LU [41]	NKG2016-LU [10]	NKG2016-GIA [10]	D1 [42]	ICE-6G-(VM5a) [43]	BIFROST-2016 [10]
GOA (MIDAS)	0.20 ± 0.31	0.29 ± 0.34	0.33 ± 0.55	−0.18 ± 0.42	−1.14 ± 0.88	0.34 ± 0.21
GOA (Tsview)	0.07 ± 0.36	0.15 ± 0.40	0.20 ± 0.62	−0.31 ± 0.46	−1.28 ± 0.86	0.28 ± 0.22
BERN (MIDAS)	0.05 ± 0.43	0.13 ± 0.46	0.18 ± 0.64	−0.33 ± 0.51	−1.30 ± 0.91	0.11 ± 0.27
BERN (Tsview)	0.16 ± 0.39	0.25 ± 0.45	0.29 ± 0.67	−0.22 ± 0.51	−1.18 ± 0.86	0.28 ± 0.19

). It could be seen that the fit to the models (based on the standard deviations) was slightly better using the GOA time-series, although systematic biases (based on the mean) in the majority were smaller with the BERN time-series. The results were also similar with the Hector solutions (Table 5 and Table 6), which shows that MIDAS was a considerable alternative to the Hector and CATS. The velocities and uncertainties from the MIDAS solution are given in Table S6 of the Supplementary Materials.

The differences between the MIDAS GOA and GOA_2 solutions were smallest using the FN + WN solution (within 0.05 ± 0.16 mm/yr on average) for the velocities and the PLN + WN solution (0.06 ± 0.03 mm/yr on average) for uncertainties. The differences between MIDAS BERN and BERN_2 solutions were smallest with the GGM + WN solution (within 0.05 ± 0.11 mm/yr on average) for the velocities and the PLN + WN solution (0.12 ± 0.04 mm/yr on average) for uncertainties. Considering the fact that the PLN + WN noise model did not fit with the Up component of GOA_2 and BERN_2 solutions, (Figure 5 and Figure 6), it can be concluded that the MIDAS solution provides too large estimates for Up component velocity uncertainties (Table S6 of Supplementary Materials).

3.4.4. Tsview Solution

As a by-product from the solutions GOA_1 and BERN_1, where outliers were removed from the time-series using Tsview and three-sigma criterion and then processed using Hector, velocities and their uncertainties (using RealSigma option) were obtained also from the Tsview solution. Obtained velocities were compared also with LU and GIA models (Table 7). As with all previous solutions, fit to the GOA time-series was slightly better.

Compared to the previous solutions, differences between Tsview solutions were smallest with the WN solution for the velocities (within 0.01 ± 0.02 mm/yr) and with the GGM + WN solution for uncertainties (also within 0.01 ± 0.02 mm/yr). The latter can be explained by the fact that both, Tsview’s RealSigma algorithm and Hector’s GGM noise models assume that the noise process is the first order Gauss–Markov process [28,49].

3.5. Comparison of Different Softwares

Velocities and their uncertainties from best solutions of GOA and BERN time-series (GOA_2 and BERN_2) were compared with each other as well as with the MIDAS and Tsview results (velocities and uncertainties are presented in Table S6 of the Supplementary Materials). Differences between different solutions (Tsview, Hector or MIDAS) within the same time-series (GOA or BERN) showed good consistency, being within −0.01 ± 0.14, −0.02 ± 0.15, 0.00 ± 0.18 mm/yr on average for North, East, and Up components, respectively. The largest velocity differences were found with the MIDAS solution. Velocity uncertainty differences were within 0.00 ± 0.11, −0.01 ± 0.10, and 0.01 ± 0.45 mm/yr on average for North, East, and Up components, respectively. Again, the largest velocity uncertainties were obtained with the MIDAS solution.

However, significant systematic difference −0.47 ± 0.08 mm/yr in North component velocity between GOA_2 and BERN_2 solutions was found (GOA_2 minus BERN_2). Note that this finding was not consistent with the previous studies, where differences between GOA and BERN solutions remained within ±0.1 mm/yr [21,37]. It seems that the North velocity component of either GOA, BERN or both was biased. However, the BIFROST2016 solution showed smaller residuals of North with the GOA time-series. This bias needs further verifying. Differences in the East and Up components were more or less random: average differences were 0.17 ± 0.13 and 0.06 ± 0.12 mm/yr, respectively. Differences between GOA_2 and BERN_2 velocities using best noise models FN + WN for the N and E components, and GGM for the U component, are presented in Figure 11.

Differences between velocity uncertainties of GOA_2 and BERN_2 solutions depended on the noise model used. The largest differences were between the PLN + WN model’s Up component uncertainties, 0.60 ± 0.08 mm/yr on average. The reason for that was explained before: the GOA time-series Up component’s noise did not fit as well to the PLN model as did the BERN time-series Up component. The uncertainty differences between GOA_2 and BERN_2 solutions based on best noise models are presented in Figure 11.

The GOA_2 and BERN_2 FN + WN and GGM model solutions were also compared with the corresponding MIDAS results. Firstly, comparison with the GOA time-series (Figure 12) shows small systematic differences in velocities and larger systematic differences in the Up component of velocity uncertainties. As was explained above, it seems that MIDAS produces too large velocity uncertainties for Up component of GOA time-series.

Secondly, comparison of the BERN time-series (Figure 13) displays more random differences in velocities, compared to the GOA time-series (Figure 12), and systematic differences in velocity uncertainties, especially in the Up component, although these differences are smaller compared to the GOA time-series.

The GOA time-series’ MIDAS solution was compared also with the BERN time-series’ MIDAS solution (Figure 14). The same systematic difference in the North and East components’ velocity is obvious, as in Figure 11. Also, systematic differences in velocity uncertainties are obvious, especially in the Up component.

4. Conclusions

This paper gave an overview of Estonian cGNSS station time-series, their noise properties, velocities, and their uncertainties. Daily coordinates of ten cGNSS stations were calculated using different processing methods (e.g., PPP and DD using GOA and BERN software). The obtained nine-year-long time-series (2008–2016) were analyzed using Tsview, Hector, and MIDAS software. Different outlier removal strategies were used to clean the time-series, after which selection of the colored noise models (e.g., white, power law, etc.) were tested for estimation of the stations’ velocities and their uncertainties. In addition, noise properties of the time-series were obtained.

Velocity and its estimated uncertainty based on cGNSS time-series was affected by the choice of the computational strategy (e.g., PPP or DD), used software, and stochastic noise model. The results show that velocities and their uncertainties depend on the used method (PPP or DD) to produce daily position time-series. Differences reached up to 0.5 mm/yr in the North component. Differences were systematic and suggest that the North component velocity of GOA or BERN or both solutions was biased. This bias needs further verifying. For example, the influence of selected fiducial stations and use of the whole network of 56 cGNSS stations in the BERN adjustment on the coordinates of analyzed 10 stations can be tested. Differences in the East and Up components were much smaller and did not exceed 0.2 mm/yr.

In addition, various outlier removal strategies to clean the time-series were tested. It was found that different outlier removal strategies had very low or no impact on velocities and their uncertainties. However, it had an impact on the AIC and BIC values, which were used to select the best solution amongst solutions from different noise models.

Results also showed that the used noise model had a rather significant impact on velocity and their uncertainty estimates. The impact was larger for the GOA time-series, which can be explained with the “more colored” noise of the GOA time-series, i.e., the absolute values of spectral indexes were larger compared to the BERN time-series. The largest uncertainties were obtained for the Up component using the PLN + WN model. However, these uncertainties are rather unrealistic due to the fact that the PLN model was not fitting well with the residual power spectral densities of the GOA time-series, which was flatter on low frequencies than the PLN model expected. The largest velocity uncertainties in BERN_2 solutions were obtained from the FN + RW model. Again, the main reason for that was the misfit between the power spectral densities of the time-series residuals and the curve of the theoretical noise model (Figure 6). Uncertainties from the colored noise models showed six times larger uncertainties on average, compared to uncorrelated white noise. The power spectral density plots of the residuals (Figure 3, Figure 4, Figure 5 and Figure 6) indicate large spectral index values for power-law type of noise. As a conclusion, the pure white noise model probably gives unrealistic results in time-series analysis of GNSS data.

Results of the noise properties analysis showed that a significantly larger proportion of WN was estimated in the BERN time-series, compared to the GOA time-series. On the other hand, a significant fraction of the RW noise was found in the GOA time-series Up component. In general, RW noise can be related to local station-dependent effects and monument instability [60,61,62]. However, RW noise was not present in the BERN time-series. One explanation could be that some station-dependent effects, common to all stations, were reduced during DD processing, but not in PPP.

Within the same time-series (GOA or BERN), velocity estimates between Tsview, Hector, and MIDAS were small, being within ±0.2 mm/yr on average. However, differences between velocity uncertainties of Tsview, Hector, and MIDAS solutions were larger, ranging up to 1 mm/yr for GOA_2 solution. The MIDAS solution produced the largest uncertainties, followed by the PLN + WN model uncertainties. This indicates that the uncertainty estimated with MIDAS was too pessimistic for the GOA time-series than the one estimated with a FN + WN or GGM noise model, in particular for the Up component.

Results of the velocity estimation were verified using LU models NKG2005LU, NKG2016LU, and GIA models NKG2016GIA, D1, and ICE-6G(VM5a). Also, velocities from BIFROST2016 solution were used for comparison. Amongst compared LU models, the best fit to the velocities was obtained with the NKG2005LU model; the fit to the newest model NGK2016LU was slightly worse.

Among compared GIA models, the North, East, and Up components’ velocities were available only for the NKG2016GIA and ICE-6G(VM5a) models. Up component velocities were available for analysis for the D1 model. However, due to the different reference frame issues, only Up component velocity comparison was carried out with GIA models. The best fit to the cGNSS time-series was obtained with the D1 model. The systematic shift and standard deviation of the residuals from GOA_2, BERN_2, MIDAS, and Tsview solutions with the majority of compared GIA and LU models were similar. The standard deviation and systematic shift were largest with the ICE-6G(VM5a) model. Systematic shift with the cGNSS velocities is related to the different reference frame issues.

Fit to the BIFROST2016 solution was better with the BERN_2 solution. Fit of the MIDAS solution to the LU and GIA models was similar with the GOA_2 and BERN_2 solutions. The best fit was obtained with the NKG2005LU model and the BIFROST2016 solution. Tsview and MIDAS solutions based on the GOA time-series showed better fit to the models compared to the BERN time-series.

It was found that, generally, the best noise model describing the data was a combination of FN and WN. Although our sample of GNSS stations was small, we found that the different noise models for estimation of the horizontal and vertical components and their uncertainties of the GNSS time-series are preferred. Based on the fit between theoretical and empirical noise models, FN plus WN model was preferred for estimation of the North and East components’ velocities and their uncertainties. GGM was a favored noise model for estimation of the Up component.