Estimation of Height Changes of Continuous GNSS Stations in the Eastern Anatolia Region During the Seasonal Variation

1. Introduction

Positioning accuracy with GNSS data analysis in geodesy is 1–2 mm in horizontal coordinates and 5–10 mm in vertical coordinates [1,2,3]. There are two main reasons for the weakness in the vertical component. The first is due to the geometric distribution of satellites in the sky. The latter is due to tropospheric path delay, particularly water vapor [4,5,6,7]. In order to increase the usability of height information obtained in today's highly sensitive geodetic studies, it has become necessary to investigate the disruptive effects on GNSS heights since many error sources first manifest themselves in this component. Increasing the sensitivity of GNSS ellipsoidal heights is possible by developing the most appropriate measurement and evaluation strategies according to the experiment and purpose, taking into account the error sources affecting the heights determined by the GNSS technique. GNSS velocities have improved significantly over the last 20 years due to studies on noise characteristics [8,9], ionospheric effects [10], or multipath and geometry effects [11].

Many GNSS time series analyses can be found in the literature review [12,13,14,15,16,17,18]. The first investigation of the permanent GNSS station time series was done in the 1990s [19,20]. Santamaria-Gomez et al. [17,21] developed different methodologies to estimate vertical velocities in the study. The effect of different antenna calibration models, the effect of different tropospheric delay models, and the geometric distribution of stations in the studied network geometry affect vertical velocities. Also, velocity domains were evaluated by analyzing the type and amplitude of noise content in time series.

A more accurate estimation of the seasonal signal is required to distinguish signals produced only by tectonic movement from other signals. For this reason, Ming et al. [22] performed seasonal signal analysis and long-term trend analysis in height time series at IGS stations. Strong seasonal variations in amplitudes of 4–9 mm were observed in the height component at the IGS stations. In previous studies of GNSS time series, the seasonal signal was thought to have annual or semi-annual periods, which can be described by a harmonic pattern of constant amplitude and phase [12,18,20,23].

Time series analysis employs a variety of techniques to separate tectonic displacement signals from seasonal signals, which are a result of the signal's nature and other factors. Davis et al. [15] proposed a linear Kalman filter model to explain the behavior of a seasonal signal, assuming that the stochastic model includes random noise. However, the hypothesis is implausible since many studies have proven that GNSS time series have flicker noise (FN). In particular, GNSS station velocities are often estimated by assuming the presence of white, flicker, and random walk noise [12,20,24].

Significant semi-annual and annual variations are visible in the height time series created using continuous GNSS station data over several years [25]. Seasonal movement from continuous GNSS observations can be a source of error in the GNSS coordinate time series and needs to be modeled to improve the accuracy of GNSS positioning. One of the important parameters affecting the accuracy of the velocity components at each station is the observation (session) time. The purpose of using different observation times is to increase station location accuracy. In order to minimize seasonal effects, it is recommended to estimate velocities with the help of at least 3 years of observations [26,27].

Wang et al. [28] examined the relationship between climate change and the height component of GNSS stations with the data for July and December from GNSS stations that measure continuously from different networks in Taiwan. As a result of the research, GNSS height values showed that the monthly averages of July values are higher than the interval values, and the difference reaches 2 cm, and it reaches 6 cm in the daily averages.

The main objectives and results of our study are to: first, present a velocity field in the Eastern Anatolia region of Turkey from a continuous GNSS network that has operated for almost 6 years. We obtained realistic uncertainties as well by taking the noise characteristics into account. Also, to improve the accuracy of the height component of the GNSS stations, the effects of seasonal variation are investigated. Since it will not be enough to investigate the effects of seasonal changes alone, it is also intended to increase the accuracy of the height component of GNSS stations with meteorological parameters. For this reason, we used two different approaches to see the effects of meteorological parameters on the vertical component. The first is simple linear regression analysis, and another is the ARMA model.

2. Materials and Methods

2.1. Study Area

The Eastern Anatolia Region is a tectonically active region where many geodetic and geophysical scientific studies have been carried out for many years. A part of these scientific studies is to reveal the current velocity field of the region. In this study, time series analyses were made with the MATLAB program after evaluating the data collected at GNSS stations continuously for about 6 years in the 2014–2019 period with a scientific software, GAMIT/GLOBK. It is aimed to reveal the periodic changes in the GNSS observations with the analyses made and to investigate the effects of these changes in the GNSS studies that have been made before and will be made in the Eastern Anatolia Region. Data from 37 continuous GNSS stations were used (Figure 1). The main subject of the study is to determine the periodic effects as a result of the analysis of the time series of these 37 GNSS stations from which data was collected and to develop strategies to improve the height component of GNSS stations.

In this study, 37 continuous GNSS stations covering the 2014–2019 date range from the Turkey National Permanent GNSS Network-Active (TUSAGA-Active) were used. The stations were chosen from the Eastern Anatolia region of Turkey because it is easier to interpret the effects of the height component due to their topographic features and also because seasonal changes are more effective there than in other parts of the country. By comparing GNSS measurements taken in various months and years, it was possible to investigate the impact of seasonal variation on station location accuracy. The locations of the GNSS stations are given in Figure 1. Eight of the stations used in the application were placed on the terrace; twenty-one of them are on the roof; and three of them are concrete poles. The stations installed on the ground are 2 meters-long concrete pole; the poles installed on the terraces are 3 meters; and the steel poles installed on the roofs are 4 meters. Normally, the stations on the building are expected to be more stationary, but when the time series of the stations are examined, it is seen that this is not the case.

Furthermore, average temperature, pressure, relative humidity, wind speed, and precipitation data from 37 Meteorology General Directorate stations in the study area were used to show the relationship between meteorological parameters and the GNSS height component. In the selection of these stations, attention was paid to ensuring that the data was continuous and in working condition.

2.2. GNSS Time Series in Height Component

Time series analyses are performed for two main purposes: First of all, it is to make predictions for the future by explaining the past with the help of observation values from past periods. The other purpose is to determine the effect of any factor that is effective in the series. This effect can be a momentary event or a set of events spread over a certain period of time. In this way, the degree of impact is determined, and measures are taken or regulations are made for the future. Some changes and irregularities can be seen over the course of the series. These changes and irregularities are caused by the existence of four factors (components) and the difference in the direction and intensity of their effects. These are the trends: seasonal fluctuations, cyclical fluctuations, erratic movements: noise [29,30].

The time series of GNSS stations are obtained as a result of daily processes. Linear, periodic, and irregular movements of reference stations can be determined by time series analysis. The time series of these reference stations can be written for horizontal and vertical components as follows [31,32];

$$\stackrel{\_}{x}(t_i)=\underset{\mathrm{Trend} \mathrm{Component}}{\underbrace{{\displaystyle \sum }_{k=1}^m{a}_t{t}_i^{k-1}}}+\underset{\mathrm{Periodic} \mathrm{Component}}{\underbrace{{\displaystyle \sum }_{s=1}^r\left[b_s\cos (2\pi f_s{t}_i)+c_s\sin (2\pi f_s{t}_i)\right]}}+\underset{\mathrm{Stochastic} \mathrm{Component}}{\underbrace{{\displaystyle \sum }_{j=1}^p{\alpha }_j\stackrel{\_}{x}(t_{i j})+{\displaystyle \sum }_{jj=1}^q{\beta }_{jj}\epsilon (t_{i\hspace{0.17em}jj})}}+\underset{\mathrm{Residual}}{\underbrace{\epsilon (t_i)}}$$

(1)

where $a_k$ is the trend component parameter, $b_s,\hspace{0.17em}{c}_s$ are amplitude of periodic signals and $f_s$ is the corresponding frequency, ${\alpha }_j$ is an autoregressive (AR) model parameter, ${\beta }_{jj}$ is moving-average (MA) model parameter and $\epsilon \left(t_i\right)$ is residual (noise) to which we pay more attention.

The trend component in the time series represents the time-dependent long-period changes of the series [33]. The Least Squares Method (LSM) separates the trend component from the series in the time series. However, LSM cannot predict the parameters of the time series components in the equation with outliers. The main disadvantage of LSM is its sensitivity to outliers [31].

An outlier is an observation that significantly deviates from other observations in the sample [34]. Prior to the application, the coordinate time series were cleaned for outliers. Outliers should be removed from the time series before starting the analysis. Outliers in the time series were investigated using a robust method. According to the predicted model, outliers were determined and removed from the time series, then gaps in the data were calculated and the time series were reconstructed [35].

Outliers in time series are eliminated from the series by using Bi-square weighted robust predictor methods. New values are determined for the eliminated values in the series. Also, the loss of data can be determined by this method. Bi-square weights minimize a weighted sum of squares, where the weight given to each observation point depends on how far the point is from the fitted line. Points near the line get full weight, whereas points far from the line get reduced weight. Robust fitting with bi-square weights uses an iteratively re-weighted LS algorithm and follows the weight function [35].

$$w_i= \{\begin{array}{l}for \left|v_i\right|\le r {\left[1-{\left(\frac{v_i}{r}\right)}^2\right]}^2\\ for \left|v_i\right|>r 0\end{array} r=4.685$$

(2)

Figure 2. Outlier test for RZE1 station.

Figure 2. Outlier test for RZE1 station.

Outliers of the stations were removed, the gaps were calculated according to the model, and then the frequencies contained in the series and the period values corresponding to these frequencies were calculated by applying Fast Fourier Transform (FFT) to the height time series where the trend component was removed [36]. Since the peaking values in the graphs differ, the frequency range has been changed according to the suitability of the graphs. Figure 3 shows the power-frequency relationship of GNSS stations of different ellipsoid heights.

In order to better understand and interpret the results obtained, the periods included in each time series were calculated using the $f=1/T$ equation from the frequencies with high power values in the power-frequency graphs of the stations. Annual periodic movement in the vertical only at SIVS1 and TOK1 stations; semi-annual periodic movement at all stations except ADIY1, ERZ2, and TUF1 stations; seasonal periodic movement at all stations except ARTV, BING, DIV1, FASA, GEME, GUMU, GURU, and TRBN GNSS stations; and much daily periodic motion was observed at all GNSS stations [14,37].

2.3. Simple Linear Regression

The relationship between a single independent variable (x) and the dependent variable (Y) expressed with a linear function is defined as simple linear regression analysis [38].

$$Y_i=a+b{x}_i+{\epsilon }_i$$

(3)

where the value of $a$ is the value of the dependent variable $Y$ when $x=0$ or in other words, the site where the vertical axis intersects with $Y$, while the value of $b$ is expressed as the regression coefficient. The difference between the true value $Y$ of the dependent variable and its predicted value $\stackrel{\wedge }{Y}$ from the model is the error $\left(\epsilon \right)$.

$$\epsilon =Y_i-(a+b{x}_i)=Y_i-\stackrel{\wedge }{Y_i}$$

(4)

In the regression, we try to estimate the independent variable Y by taking the dependent X variables. The results we get as a result of estimation are usually either incomplete or wrong. The real question we have to ask here is how wrong it is. In other words, what we really need to do is find the distance between the actual values and the predicted values. The following are the concepts used to evaluate the regression model [39].

2.3.1. Model Performance

The mean absolute error represents the average of the absolute difference between the actual and predicted values in the dataset. It measures the average of the residuals in the dataset.

$$MAE=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left|y_i-\stackrel{\wedge }{y}\right|}$$

(5)

(best value = 0; worst value = +∞)

where; $\stackrel{\wedge }{y}$ predicted value of $y$ and $\overline{y}$ means value of $y$.

The mean squared error represents the average of the squared difference between the original and predicted values in the data set. It measures the variance of the residuals.

$$MSE=\frac{1}{N}{{\displaystyle \sum _{i=1}^N\left(y_{i-}\stackrel{\wedge }{y}\right)}}^2$$

(6)

(best value = 0; worst value = +∞)

The root mean squared error is the square root of mean squared error. It measures the standard deviation of residuals.

$$RMSE=\sqrt{MSE}=\sqrt{\frac{1}{N}{{\displaystyle \sum _{i=1}^N\left(y_{i-}\stackrel{\wedge }{y}\right)}}^2}$$

(7)

(best value = 0; worst value = +∞)

The coefficient of determination or R-squared, represents the proportion of the variance in the dependent variable that is explained by the linear regression model. It is a scale-free score, irrespective of the values being small or large, the value of R square will be less than one.

$$R^2=1-\frac{{\displaystyle \sum {\left(y_i-\stackrel{\wedge }{y}\right)}^2}}{{\displaystyle \sum {\left(y_i-\stackrel{\wedge }{y}\right)}^2}}$$

(8)

(worst value = −∞; best value = +1)

Both RMSE and R-squared quantify how well a linear regression model fits a dataset. The RMSE tells how well a regression model can predict the value of a response variable in absolute terms, while R-squared tells how well the predictor variables can explain the variation in the response variable [40].

2.4. The Autoregressive Moving-Average (ARMA)

$ARMA$ models are used for modeling stationary data and are a combination of $AR$ and $MA$ models. In these models, the observation value of any period of a time series is expressed as a linear combination of a certain number of previous observation values and the error term. If the $ARMA$ model is a combination of $AR\left(p\right)$ and $MA\left(q\right)$ model, it contains $p+q$ terms and is expressed as $ARMA\hspace{0.17em}\left(p,q\right)$ [41].

$$x_t={\displaystyle {\sum }_{i=1}^p{\varphi }_i}{x}_{t-i}-{\displaystyle {\sum }_{j=1}^q{\theta }_j}{a}_{t-j}+a_t;t=1,\mathrm{......},T$$

(9)

Here, $a_t$ is known as normal white noise process. It has a zero mean and variance ${\sigma }^2$. $T$ is the amount of data in the time series [42]. $AR$ parameters should satisfy the condition for stationarity and $MA$ parameters should satisfy the conditions for inevitability [43].

ARMA model was applied to the height component of GNSS stations in the same region by using the daily average temperature, pressure, relative humidity, wind speed, and precipitation values measured at 37 meteorology stations in the study area. Model results for four GNSS stations are given in Figure 4 as an example, and the model results for all stations subject to the study are given in Table 1.

ARMA model results are generally consistent with the time series of all stations. The RMS values are generally less than 5 mm in the vertical directions, except for ARTV, ERZ2, GIRS, MUUS, and SIRT1 stations. The RMS values of the vertical component were computed with all meteorological data.

3. Results and Discussion

3.1. Seasonal Variation

In this study, the effects of seasonal variation on the vertical position accuracy of GNSS calculated by time series analysis of continuous GNSS stations were investigated. Weather changes and water vapor in the atmosphere affect the position accuracy of GNSS and cause fluctuations in GNSS height values. It is also known that the height component has more air passage changes. By applying time series analysis to the daily coordinate values of the stations, statistically significant trends, periodic, and stochastic components of the stations were determined. As a result of the analysis, the vertical annual velocities of the stations and the standard deviations of the velocities were determined.

For the stations determined according to the ellipsoid heights, the velocity and standard deviation values of the height component were calculated for each month, season, and year. As the ellipsoid height increases, the velocity and its standard deviation values decrease. While the minimum vertical velocity values are observed for the station with the lowest ellipsoidal height in winter, for the station with the highest ellipsoidal height in autumn, the minimum standard deviation values are determined in winter for the station with the lowest ellipsoidal height, and in summer for the station with the highest ellipsoidal height. According to the results obtained, the coordinate displacements caused by seasonal variation may be important, and their effects should be considered, especially in high-precision geodetic surveys [44].

Figure 6. Vertical velocity variations of GNSS stations (The minimum velocity for RZE1 station in winter (December), for DIYB station in winter (December), for GEME station in summer (August), and for HINI station in autumn (November). The maximum velocity for RZE1 station in autumn (September), for DIYB station in autumn (October), for GEME station in spring (April), and for HINI station in winter (December)).

Figure 6. Vertical velocity variations of GNSS stations (The minimum velocity for RZE1 station in winter (December), for DIYB station in winter (December), for GEME station in summer (August), and for HINI station in autumn (November). The maximum velocity for RZE1 station in autumn (September), for DIYB station in autumn (October), for GEME station in spring (April), and for HINI station in winter (December)).

Figure 7. Standard deviation variations for vertical velocities of GNSS stations (The minimum standard deviation for RZE1 GNSS station in winter (December), for DIYB GNSS station in winter (February), for GEME GNSS station in spring (March), and for HINI GNSS station in summer (June). The maximum standard deviation for RZE1 GNSS station in summer (July), for DIYB GNSS station in summer (August), for GEME GNSS station in summer (August), and for HINI GNSS station in autumn (September)).

Figure 7. Standard deviation variations for vertical velocities of GNSS stations (The minimum standard deviation for RZE1 GNSS station in winter (December), for DIYB GNSS station in winter (February), for GEME GNSS station in spring (March), and for HINI GNSS station in summer (June). The maximum standard deviation for RZE1 GNSS station in summer (July), for DIYB GNSS station in summer (August), for GEME GNSS station in summer (August), and for HINI GNSS station in autumn (September)).

In addition, the velocity values of the stations were calculated for different years, and a decrease was observed in the height component depending on the observation duration. As the observation duration for the height component increases, both the velocity values and their standard deviation values decrease. In order to avoid velocity estimation errors completely, the data length should be more than 4.5 years.

It has been determined that the height components of the GNSS stations used in the application make periodic movements. Considering that this may be due to seasonal effects, daily average temperature, pressure, relative humidity, wind speed, and precipitation values between 2014 and 2019 of meteorological stations in the cities where the GNSS stations are located were obtained from the General Directorate of Meteorology. By using these parameters, the relationship of the stations with the height component was investigated.

Table 2. Regression results for GNSS stations of different ellipsoidal heights (m).

temperature	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.781	0.012	0.0880		0.01	0.1135	70.69
DIYB	0.864	0.009	0.0943		0.01	0.1173	773.67
GEME	0.850	0.008	0.1072		0.02	0.1371	1214.79
HINI	0.839	0.008	0.1006		0.02	0.1302	1742.62
pressure	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.845	0.019	0.1106		0.02	0.1311	70.69
DIYB	0.801	0.014	0.0987		0.02	0.1250	773.67
GEME	0.894	0.009	0.1122		0.02	0.1458	1214.79
HINI	0.864	0.008	0.1004		0.02	0.1288	1742.62
relative humidity	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.843	0.019	0.1087		0.02	0.1136	70.69
DIYB	0.797	0.014	0.0946		0.01	0.1184	773.67
GEME	0.894	0.009	0.1021		0.02	0.1338	1214.79
HINI	0.861	0.008	0.0986		0.02	0.1248	1742.62
wind speed	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.845	0.019	0.1029		0.02	0.1267	70.69
DIYB	0.794	0.014	0.0968		0.01	0.1223	773.67
GEME	0.893	0.009	0.1096		0.02	0.1404	1214.79
HINI	0.861	0.008	0.0967		0.02	0.1225	1742.62
precipitation	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.847	0.019	0.1122		0.02	0.1371	70.69
DIYB	0.796	0.014	0.0962		0.01	0.1212	773.67
GEME	0.895	0.009	0.1095		0.02	0.1396	1214.79
HINI	0.858	0.008	0.0952		0.01	0.1201	1742.62

Table 2. Regression results for GNSS stations of different ellipsoidal heights (m).

temperature	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.781	0.012	0.0880		0.01	0.1135	70.69
DIYB	0.864	0.009	0.0943		0.01	0.1173	773.67
GEME	0.850	0.008	0.1072		0.02	0.1371	1214.79
HINI	0.839	0.008	0.1006		0.02	0.1302	1742.62
pressure	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.845	0.019	0.1106		0.02	0.1311	70.69
DIYB	0.801	0.014	0.0987		0.02	0.1250	773.67
GEME	0.894	0.009	0.1122		0.02	0.1458	1214.79
HINI	0.864	0.008	0.1004		0.02	0.1288	1742.62
relative humidity	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.843	0.019	0.1087		0.02	0.1136	70.69
DIYB	0.797	0.014	0.0946		0.01	0.1184	773.67
GEME	0.894	0.009	0.1021		0.02	0.1338	1214.79
HINI	0.861	0.008	0.0986		0.02	0.1248	1742.62
wind speed	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.845	0.019	0.1029		0.02	0.1267	70.69
DIYB	0.794	0.014	0.0968		0.01	0.1223	773.67
GEME	0.893	0.009	0.1096		0.02	0.1404	1214.79
HINI	0.861	0.008	0.0967		0.02	0.1225	1742.62
precipitation	R-squared	Standard error	MAE		MSE	RMSE	Height
RZE1	0.847	0.019	0.1122		0.02	0.1371	70.69
DIYB	0.796	0.014	0.0962		0.01	0.1212	773.67
GEME	0.895	0.009	0.1095		0.02	0.1396	1214.79
HINI	0.858	0.008	0.0952		0.01	0.1201	1742.62

The estimation of the parameters was carried out with the least squares method. Hypothesis tests determined the suitability of the regression model and the significance of the regression coefficients.

3.2. Determination of Noise Patterns

Noise analysis is needed to reveal the characteristics of the stations and to investigate the time series more deeply. [12,13,20,24,29] emphasized that GNSS time series contain time-independent white noise (WN), time-dependent flicker noise (FN), and random walk noise (RWN). In all these analyses, it was observed that the vertical component had more noise than the horizontal component, and the periodic motion of the vertical component was more pronounced.

The most suitable noise model for the height component of the stations was examined with the CATS software [45]. CATS is software that uses Maximum Likelihood Estimation (MLE) to fit a multi-parameter model to time series such as fixed GNSS station coordinates. Using MLE, it is possible to estimate the amplitudes of multiple noises and the unknown (e.g., velocity, periodic signals, and coordinate pulses) simultaneously [12,36].

Three noise types and combinations of these noise types were used in the study carried out to reveal the correlated noise models. First, the noise was assumed to be only WN, then a combination of WN+FN and WN+RWN was used. The preferred noise pattern was determined to be one of these three combinations. In the second stage, the amplitudes and spectral indices of the noise model were presented simultaneously with the WN.

As a result of the analysis made, the most appropriate noise model for the height component is given in Table 3. When the MLE values are examined, it is seen that the WN+FN model is the most dominant model among the three noise models, as suggested by [12,20,24]. Again, for all stations, the WN+FN noise model was found to be more suitable than the WN+RWN noise model [46]. In addition, when the results are examined, the values of WN+FN and WN+RWN noise combinations for all stations are quite close to each other.

Amplitudes are one of the important parameters for noise analysis. The amplitudes indicate the magnitude of the appropriate noise pattern present in the available data. The calculated amplitudes for the noise models are given in Figure 8. It is seen that the lowest noise amplitudes are for WN. It is considered that the reason for the WN amplitude differences is the use of different GNSS receivers at the stations. A significant amount of WN amplitudes and seasonal atmospheric mass distributions are considered to reflect more general physical-based processes, such as atmospheric noise. A very small amount of RWN amplitude also results from more regional-scale processes such as atmospheric effects [12,25,47].

As a result of the analysis, the velocities for the height component of the stations were obtained. When the annual velocity of the height coordinate components of the stations is examined, a vertical velocity value in the (+) direction was observed at ADIY, ADIY1, EKIZ, ERGN, ERZR, and HINI GNSS stations; that is, the height value increased. While the highest vertical velocity value was seen at station ADIY1 with a value of 4.92 mm/year, the smallest vertical velocity value was obtained at station ERZR with a value of 0.03 mm/year. At other stations, a vertical velocity value in the (-) direction was observed, that is, the height value decreased. The smallest vertical velocity value was -0.03 mm/year at SIRT GNSS station, and the vertical velocity was -7.10 mm/year at UDE1 GNSS station. The vertical velocities of GNSS stations are given in Figure 9.

4. Conclusions

The goal of this paper is to find a way to increase the accuracy of the height component of continuous observation GNSS stations. For the vertical displacement analysis of the region, we selected 37 GNSS stations and processed approximately 6 years of continuous data.

First of all, data arrangements were made, and an outlier test was used in MATLAB ver. R2015b (The MathWorks, Inc., USA). After data editing was completed, time series were plotted, and anomalies were investigated. The seasonal behavior of GNSS stations was examined, and strong seasonal and interannual signals were found at all stations. It is also clearly seen that GNSS signals have a lot of noise. CATS software is used for spectral analysis. As expected, the most appropriate noise model for the station was found to be the WN+FN noise model. Then, the relationship between seasonal variation and vertical time series was examined. As a result of the analysis, the vertical seasonal and annual velocities of the stations and the standard deviations of the velocities were determined. The results are variable for stations with different ellipsoid heights. However, as expected, as the observation time increased, the velocity and standard deviation values improved. Furthermore, when the height components of continuous GNSS stations were examined, it was seen that there were seasonal effects, and it was investigated whether the height components were related to meteorological parameters. Simple linear regression analysis was performed to determine how dependent the height components of the continuous GNSS stations were on meteorological parameters. As a result of the analysis, the height components of the continuous GNSS stations are dependent on meteorological parameters such as temperature, pressure, relative humidity, wind speed, and precipitation. In addition, height component time series analysis of continuous GNSS stations was performed by using Autoregressive Moving Average (ARMA) models from linear time series methods. As a result of the study, the performance of the ARMA modeling results again indicated the dependence of the height component of the continuous GNSS stations on the meteorological parameters.

The consideration of the above-mentioned analysis may improve the quality of the station coordinate estimates, especially for the up component.