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

: Estimating the height component of Global Navigation Satellite System (GNSS) stations is widely known to be more challenging than estimating the horizontal position. In this study, we utilized height time series data from 37 continuous GNSS stations that were part of the Turkish RTK CORS Network called TUSAGA-Active (Turkish National Permanent GNSS Network Active). The data covered the period from 2014 to 2019, and the selection of stations focused on the Eastern Anatolia region of Turkey due to its topographic characteristics and the pronounced inﬂuence of seasonal changes, which facilitated the interpretation of the effects on the height component. The daily coordinates of the GNSS stations were derived using the GAMIT/GLOBK software solution. We identiﬁed statistically signiﬁcant trends, periodic variations, and stochastic components associated with the stations by applying time series analysis to these daily coordinate values. As a result, the vertical velocities of the GNSS stations were determined, along with their corresponding standard deviations. Furthermore, examining the height components of the continuous GNSS stations revealed seasonal effects. We aimed to investigate the potential relationship between these height components and meteorological parameters. The study provides evidence of the interconnectedness between the height components of continuous GNSS stations and various meteorological parameters. Simple linear regression analysis and ARMA time series modeling were utilized to establish this relationship.


Introduction
When using GNSS data analysis in geodesy, the positioning accuracy is 1-2 mm in horizontal coordinates and 5-10 mm in vertical coordinates [1][2][3]. It is common knowledge that estimating the height component of GNSS stations is much more problematic than estimating the horizontal position. There are two main reasons for the weakness in the vertical component. The first reason is the geometric distribution of satellites in the sky. The second reason is tropospheric path delays, particularly caused by water vapor [4][5][6][7]. In order to increase the usability of the height information obtained in today's susceptible geodetic studies, it has become necessary to investigate these disruptive effects on GNSS height. Increasing the sensitivity of GNSS heights is possible by developing the most appropriate measurement and evaluation strategies according to the experiment and purpose, and considering the error sources affecting the heights determined. 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 a review of the literature [12][13][14][15][16][17][18]. The permanent GNSS station time series was first investigated in the 1990s [19,20]. Santamaria-Gomez et al. [17,21] developed different methodologies to estimate vertical velocities. Different antenna calibration models, the effects 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 noise content's time series type and amplitude.
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 various techniques to distinguish between tectonic displacement signals and seasonal signals, which arise from the inherent nature of the signal and other factors. Davis et al. [15] proposed a linear Kalman filter model to explain the behavior of seasonal signals, assuming that the stochastic model includes random noise. However, this hypothesis is implausible, as numerous studies have demonstrated that GNSS time series demonstrate flicker noise. Specifically, GNSS station velocities are often estimated by considering 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 observed in continuous GNSS time series can be a source of error in coordinate time series and thus needs to be modeled to improve positioning accuracy. 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 position 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. The GNSS height values showed that the monthly averages of July values were higher than the interval values, and the difference reached 2 cm. The difference reached 6 cm for the daily average. Dogan et al. [29] conducted a study using 3-day datasets for each month of 2009 from the GPS stations of the Marmara Continuous GPS Network. Their objective was to examine how seasonal changes affect GPS location accuracy. They evaluated 24-h GPS measurements in the dataset with various base lengths, including 4-, 6-, 8-, and 12-h datasets. The study revealed that the accuracy of the observed stations' positions could vary each month in the north-south, east-west, and height directions. Notably, the GPS location accuracy was lowest during the summer (July) and highest in the winter (January and December). The author of [30] analyzed one-year datasets from a mid-latitude GPS station to examine the impact of seasons. The RMS values of the coordinate components were calculated, and then the researcher explored their relationship with meteorological parameters, including temperature, pressure, and humidity. The findings revealed significant correlations between the northern and eastern components of the coordinates and the temperature data. Moreover, it was observed that the RMS values were consistently higher during the spring and summer months, indicating a potential seasonal influence on the GPS station's accuracy.
The authors of [31] used time series analysis to examine the relationship between GPS coordinate changes and seasonal variation. The researchers analyzed the changes in the coordinate components of multiple GPS stations between 1996 and 2016 and studied their amplitude and RMS values. The findings revealed that the amplitude of the time series of a station exhibited a vertical increase towards the end of winter and a vertical collapse towards the end of summer. This observation was attributed to groundwater or uplift following winter snowmelt. Furthermore, the study determined that the horizontal coordinate amplitudes were 2-3 times smaller than the vertical ones. The researchers explained this difference by stating that GPS vertical position accuracy was three times lower than horizontal position accuracy and that vertical component crustal deformations exhibited larger amplitudes. The authors of [32] analyzed crustal deformations by comparing data from eight GPS stations spanning a range of 2.5 to 13 years. The analysis included time series, gravity recovery, and climate test data. The results revealed that no significant seasonal effect was observed in the horizontal component of the GPS stations. However, in the height component, it was determined that the water storage cycle in the climate zone played a significant role. The study found that when more water is stored in the spring, the crust experiences a greater load, resulting in minimal heights. Conversely, when less water is stored in late summer or early autumn, the load on the crust decreases, leading to maximum heights. This seasonal variation in height reached up to 10 mm. Therefore, seasonal variation was the primary cause of fluctuations in the height component.
The main purpose of our study was to present an analysis of the vertical velocity field within the Eastern Anatolia region of Turkey, focusing on its correlation with seasonal variations. This investigation utilizes data from a continuous GNSS network, which has been operating for approximately six years. In order to ensure reliable uncertainties, we took the noise characteristics into account. We recognize that we cannot achieve this objective by only examining seasonal changes. Our research aimed to enhance the accuracy of the vertical component by incorporating meteorological parameters. Therefore, we employed two approaches to evaluate the influence of meteorological parameters on the vertical component. The first approach is simple linear regression analysis, while the second is the ARMA model.

Study Area
The Eastern Anatolia Region has active tectonic activity that has been the focus of numerous geodetic and geophysical scientific studies. One aspect of these studies was to determine the current velocity field of the region. We evaluated the data collected from GNSS stations for approximately 6 years, from 2014 to 2019, using a scientific software called GAMIT/GLOBK [33]. We utilized data from 37 continuous GNSS stations, as shown in Figure 1. This study performed time series analyses with MATLAB software (version R2015b, Natick, MA, USA: The MathWorks, Inc.). Through these analyses, we aimed to identify periodic changes in GNSS height component observations and investigate their implications for previous and future GNSS studies in the Eastern Anatolia Region.
This study used 37 continuous GNSS stations from the Turkey National Permanent GNSS Network-Active (TUSAGA-Active). The selected stations mainly focused on the Eastern Anatolia region of Turkey, chosen for its favorable topographic characteristics that facilitate a clearer interpretation of the height component effects. Additionally, the region's pronounced seasonal variations offer greater sensitivity than other regions within the country. By comparing GNSS measurements taken in different months and years, we investigated the impact of seasonal variation on station location accuracy. Figure 1 shows the locations of the GNSS stations, with eight stations placed on the terrace, twenty-one on the roof, and three on concrete poles. The stations on the ground consist of 2 m long concrete poles, while the ones on the terraces are 3 m, and the steel poles on the roofs are 4 m. Although the stations on the building were expected to be stable, the examination of the time series of the stations revealed 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 selecting these stations, attention was paid to ensuring that the data was continuous and in working condition. Appl  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 selecting these stations, attention was paid to ensuring that the data was continuous and in working condition.

GNSS Time Series in Height Component
The time series analyses were performed for two main purposes: first, to make predictions for the future by explaining the past with the help of observation values from past periods. The second purpose is to determine the effect of any effective factor in the series. This effect can be a momentary event or a set of events spread over a certain time period. This way, the degree of impact is determined, measures are taken, or regulations are made for the future. Some changes and irregularities can be seen in the series. These changes and irregularities were caused by four factors (components) and there were differences in the direction and intensity of their effects. These are trend, seasonality, cyclicality, and noise [34,35].
The time series of the GNSS stations are obtained from 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 the horizontal and vertical components as follows [36,37]; where is the trend component parameter, and are the amplitude of periodic signals, is the corresponding frequency, is the autoregressive model parameter, is the moving average model parameter, and is the residual (noise) to which we pay more attention.

GNSS Time Series in Height Component
The time series analyses were performed for two main purposes: first, to make predictions for the future by explaining the past with the help of observation values from past periods. The second purpose is to determine the effect of any effective factor in the series. This effect can be a momentary event or a set of events spread over a certain time period. This way, the degree of impact is determined, measures are taken, or regulations are made for the future. Some changes and irregularities can be seen in the series. These changes and irregularities were caused by four factors (components) and there were differences in the direction and intensity of their effects. These are trend, seasonality, cyclicality, and noise [34,35].
The time series of the GNSS stations are obtained from 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 the horizontal and vertical components as follows [36,37]; where a k is the trend component parameter, b s and c s are the amplitude of periodic signals, f s is the corresponding frequency, α j is the autoregressive AR(p) model parameter, β jj is the moving average MA(q) model parameter, and ε(t) i is the residual (noise) to which we pay more attention.
The trend component in the time series represents the time-dependent long-period changes [38]. 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 [36].
An outlier is an observation that significantly deviates from other observations in the sample [39]. 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, the gaps in the data were calculated, and the time series were reconstructed [40].
Outliers in time series were eliminated from the series by using Bi-square weighted robust predictor methods. New values were 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 weights, whereas points far from the line get reduced weights. Robust fitting with Bi-square weights uses an iteratively re-weighted LS algorithm and follows the weight function [40]. The weight functions of robust estimators in Bi-square are: where r is a tuning constant equal to 4.685.
Outliers from the observations were removed, and the gaps were determined using the model. Subsequently, the series was reconstructed. Figure 2 shows that the new dataset for the RZE1 height is compatible with the raw dataset. In Figure 2, the black lines represent the raw data, the red stars denote the outliers, and the blue lines illustrate the new dataset with the outliers removed. Then, the Fast Fourier Transform (FFT) was applied to the height time series to calculate the frequencies and corresponding period values where the trend component was removed [41]. Since the peak values in the graphs vary, the frequency range was adjusted to ensure graph suitability. Figure 3 displays the power-frequency relationship of GNSS stations with different ellipsoid heights.
changes [38]. 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 [36].
An outlier is an observation that significantly deviates from other observations in the sample [39]. 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, the gaps in the data were calculated, and the time series were reconstructed [40].
Outliers in time series were eliminated from the series by using Bi-square weighted robust predictor methods. New values were 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 weights, whereas points far from the line get reduced weights. Robust fitting with Bi-square weights uses an iteratively re-weighted LS algorithm and follows the weight function [40]. The weight functions of robust estimators in Bi-square are: where is a tuning constant equal to 4.685.
Outliers from the observations were removed, and the gaps were determined using the model. Subsequently, the series was reconstructed. Figure 2 shows that the new dataset for the RZE1 height is compatible with the raw dataset. In Figure 2, the black lines represent the raw data, the red stars denote the outliers, and the blue lines illustrate the new dataset with the outliers removed. Then, the Fast Fourier Transform (FFT) was applied to the height time series to calculate the frequencies and corresponding period values where the trend component was removed [41]. Since the peak values in the graphs vary, the frequency range was adjusted to ensure graph suitability. Figure 3 displays the power-frequency relationship of GNSS stations with different ellipsoid heights.   In order to better understand and interpret the results obtained, the periods included in each time series were calculated using the equation from the frequencies with high power values in the power-frequency graphs of the stations. We observed annual periodic movement in the vertical component only at the SIVS1 and TOK1 stations; semiannual periodic movement at all stations except the ADIY1, ERZ2, and TUF1 stations; seasonal periodic movement at all stations except the ARTV, BING, DIV1, FASA, GEME, GUMU, GURU, and TRBN GNSS stations; and most daily periodic motion was observed at all GNSS stations [14,42].

Simple Linear Regression
Expressing the relationship between a single independent variable and the dependent variable with a linear function is defined as simple linear regression analysis [43]. (3) The value of represents the dependent variable when 0, which can also be understood as the point where the vertical axis intersects with . On the other hand, denotes the regression coefficient in the model. The error is the difference between the dependent variable's actual value and its predicted value based on the model. 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. We observed annual periodic movement in the vertical component only at the SIVS1 and TOK1 stations; semi-annual periodic movement at all stations except the ADIY1, ERZ2, and TUF1 stations; seasonal periodic movement at all stations except the ARTV, BING, DIV1, FASA, GEME, GUMU, GURU, and TRBN GNSS stations; and most daily periodic motion was observed at all GNSS stations [14,42].

Simple Linear Regression
Expressing the relationship between a single independent variable x and the dependent variable Y with a linear function is defined as simple linear regression analysis [43].
The value of a represents the dependent variable Y when x = 0, which can also be understood as the point where the vertical axis intersects with Y. On the other hand, b denotes the regression coefficient in the model. The error (ε) is the difference between the dependent variable's actual value Y and its predicted value Y based on the model. In regression analysis, we aim to estimate the independent variable Y by considering the dependent x variables. However, the results obtained from these estimations often require further refinement or correction. The key question we need to address is the extent of inaccuracy in our estimates. In other words, we must determine the gap between actual and predicted values. The following concepts are used to evaluate the regression model [44].

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.
The mean squared error represents the average difference between the actual and predicted values in the dataset. It measures the variance of the residuals.
(best value = 0; worst value = +∞). The root mean squared error is the square root of the mean squared error. It measures the standard deviation of the residuals.
(best value = 0; worst value = +∞). The coefficient of determination or R-squared represents the proportion of the variance in the dependent variable that the linear regression model explains. It is a scale-free score; irrespective of the values being small or large, the value of R square will be less than one.
(worst value = −∞; best value = +1) where ∼ Y is the mean value of Y. RMSE and R-squared quantify how well a linear regression model fits a dataset. The RMSE indicates how well a regression model can predict the value of a response variable in absolute terms. In contrast, R-squared indicates how well the predictor variables can explain the variation in the response variable [45].

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 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(p) and MA(q) models, it contains (p + q) terms and is expressed as ARMA(p, q) [46].
θ j a t−j + a t t = 1, 2, . . . , T Here, a t is known as the normal white noise process. It has a zero mean and variance σ 2 . T is the amount of data in the time series [47]. AR parameters should satisfy the condition for stationarity, and MA parameters should satisfy the condition for inevitability [48].
The 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. The 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. Here, is known as the normal white noise process. It has a zero mean and variance 2 .
is the amount of data in the time series [47]. parameters should satisfy the condition for stationarity, and parameters should satisfy the condition for inevitability [48].
The 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. The 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.    The 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 the ARTV, ERZ2, GIRS, MUUS, and SIRT1 stations. The RMS values of the vertical component were computed with all meteorological data. Figure 5 compares the height components obtained from the ARMA model and GPS measurements. Histograms are particularly useful for understanding the shape of the data distribution. In a histogram, the type of data being measured is represented on the horizontal axis, and the vertical axis represents how many observations are in each bin. The Appl. Sci. 2023, 13, x FOR PEER REVIEW 9 of 18 Figure 5 compares the height components obtained from the ARMA model and GPS measurements. Histograms are particularly useful for understanding the shape of the data distribution. In a histogram, the type of data being measured is represented on the horizontal axis, and the vertical axis represents how many observations are in each bin. The histograms of all stations had a bell-shaped distribution. The model is quite compatible and makes predictions close to the true values.

Seasonal Variation
In this study, we investigated the impact of seasonal variation on the accuracy of the vertical positioning of GNSS through the analysis of continuous GNSS stations' time series data. It is widely recognized that weather changes and the presence of water vapor in the atmosphere can influence the position accuracy of GNSS, leading to fluctuations in GNSS height measurements. Additionally, the height component is particularly susceptible to variations in air passage. To assess these effects, we conducted time series analysis on the daily coordinate values of the stations. This analysis allowed us to identify statistically significant trends, periodic patterns, and stochastic components in the station data. Our analysis yielded valuable insights, including determining annual vertical velocities and calculating standard deviations associated with these velocities.
For the selected stations based on ellipsoid heights, we calculated the velocity and standard deviation values for the height component monthly, seasonally, and yearly. Notably, as the ellipsoid height increased, we observed a decrease in the velocity and its corresponding standard deviation. Specifically, the station with the lowest ellipsoidal height exhibited the minimum vertical velocity values during winter, while the station with the highest ellipsoidal height demonstrated this minimum during autumn. Furthermore, the station with the lowest ellipsoidal height showcased the minimum standard deviation values during winter, whereas the station with the highest ellipsoidal height displayed this minimum during summer. Our findings show that velocity changes resulting from

Seasonal Variation
In this study, we investigated the impact of seasonal variation on the accuracy of the vertical positioning of GNSS through the analysis of continuous GNSS stations' time series data. It is widely recognized that weather changes and the presence of water vapor in the atmosphere can influence the position accuracy of GNSS, leading to fluctuations in GNSS height measurements. Additionally, the height component is particularly susceptible to variations in air passage. To assess these effects, we conducted time series analysis on the daily coordinate values of the stations. This analysis allowed us to identify statistically significant trends, periodic patterns, and stochastic components in the station data. Our analysis yielded valuable insights, including determining annual vertical velocities and calculating standard deviations associated with these velocities.
For the selected stations based on ellipsoid heights, we calculated the velocity and standard deviation values for the height component monthly, seasonally, and yearly. Notably, as the ellipsoid height increased, we observed a decrease in the velocity and its corresponding standard deviation. Specifically, the station with the lowest ellipsoidal height exhibited the minimum vertical velocity values during winter, while the station with the highest ellipsoidal height demonstrated this minimum during autumn. Furthermore, the station with the lowest ellipsoidal height showcased the minimum standard deviation values during winter, whereas the station with the highest ellipsoidal height displayed this minimum during summer. Our findings show that velocity changes resulting from seasonal variations can be significant, particularly in high-precision geodetic surveys. Therefore, it is crucial to account for these effects when conducting such surveys. Figure 6 shows the seasonal, monthly, and annual calculated velocities at the stations. The minimum velocity for RZE1 and DIYB was observed in December during the winter. The minimum velocity for GEME was observed in August in the summer, while it was observed in November for HINI in the autumn. On the other hand, the maximum velocities were observed at stations in different seasons. RZE1, with the lowest ellipsoid height, showed the highest velocity in autumn in September, while HINI, with the highest ellipsoid height, showed its maximum velocity in December in winter. Figure 7 shows the stations' seasonal, monthly, and annual calculated standard deviations. The minimum standard deviation was observed at the RZE1 and DIYB stations in the winter. The minimum standard deviation for GEME was observed in the spring and for HINI in the summer. The maximum standard deviation for RZE1, DIYB, and GEME was observed in the summer, while for HINI, this value was observed in the autumn. seasonal variations can be significant, particularly in high-precision geodetic surveys. Therefore, it is crucial to account for these effects when conducting such surveys. Figure 6 shows the seasonal, monthly, and annual calculated velocities at the stations. The minimum velocity for RZE1 and DIYB was observed in December during the winter. The minimum velocity for GEME was observed in August in the summer, while it was observed in November for HINI in the autumn. On the other hand, the maximum velocities were observed at stations in different seasons. RZE1, with the lowest ellipsoid height, showed the highest velocity in autumn in September, while HINI, with the highest ellipsoid height, showed its maximum velocity in December in winter. Figure 7 shows the stations' seasonal, monthly, and annual calculated standard deviations. The minimum standard deviation was observed at the RZE1 and DIYB stations in the winter. The minimum standard deviation for GEME was observed in the spring and for HINI in the summer. The maximum standard deviation for RZE1, DIYB, and GEME was observed in the summer, while for HINI, this value was observed in the autumn.   Seasonal effects and various meteorological factors such as daily average temperature, pressure, relative humidity, wind speed, and precipitation can affect periodic movements in the altitude component of GNSS stations. Considering the relationship between these factors and the stations' height component, 2014 and 2019 meteorological station data were obtained from the General Directorate of Meteorology. Table 2 presents the performance results of the simple linear regression model for each meteorological parameter. The table includes the employed model and the selected criteria for performance evaluation. Upon reviewing the results, it is evident that the model's performance improved for all parameters as the ellipsoid height increased. When  Table 2 presents the performance results of the simple linear regression model for each  meteorological parameter. The table includes the employed model and the selected criteria for performance evaluation. Upon reviewing the results, it is evident that the model's performance improved for all parameters as the ellipsoid height increased. When the R-squared values approached one and the MAE and RMSE values approached zero, the model performed successfully in all parameters. The estimation of the parameters was carried out with the least squares method. The hypothesis tests determined the regression model's suitability and the regression coefficients' significance.

Determination of Noise Patterns
Noise analysis is needed to reveal the stations' characteristics and investigate the time series more deeply. Refs. [12,13,20,24,34] 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 investigated using the CATS (GPS Coordinate Time Series Analysis Software) software [49]. CATS software employs Maximum Likelihood Estimation (MLE) to fit a multi-parameter model to time series data, such as fixed GNSS station coordinates. Through MLE, it becomes feasible to simultaneously estimate the amplitudes of various noise sources and unknown factors (e.g., velocity, periodic signals, and coordinate pulses) [12,41].
Three noise types and combinations of these noise types were used in the study 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, the most appropriate noise model for the height component is given in Table 3. When the MLE values were examined, it was seen that the WN + FN model was the most dominant model among the three noise models, as suggested by [12,20,24]. Again, the WN + FN noise model was more suitable for all stations than the WN + RWN noise model [50]. In addition, when the results were examined, the values of WN + FN and WN + RWN noise combinations for all stations were quite close. 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 can be seen that the lowest noise amplitudes were for WN. The reason for the WN amplitude differences could be the use of different GNSS receivers at the stations. A significant amount of WN amplitude and seasonal atmospheric mass distributions were thought to reflect more general physical-based processes, such as atmospheric noise. A small amount of RWN amplitude results from more regional-scale processes, such as atmospheric effects [12,25,51]. Appl  As a result of the analyses, we obtained velocities for the height component of the stations. When examining the annual velocity of the height coordinate components, we observed a positive vertical velocity at the ADIY, ADIY1, EKIZ, ERGN, ERZR, and HINI GNSS stations. The station with the highest vertical velocity was ADIY1, measuring 4.92 mm/year, while the smallest vertical velocity was observed at station ERZR, measuring 0.03 mm/year. Conversely, the remaining stations exhibited a negative vertical velocity. The smallest negative vertical velocity was −0.03 mm/year at the SIRT GNSS station, while the UDE1 GNSS station had the highest vertical velocity of −7.10 mm/year. The vertical velocities of all GNSS stations are presented in Figure 9 and Table A1.  As a result of the analyses, we obtained velocities for the height component of the stations. When examining the annual velocity of the height coordinate components, we observed a positive vertical velocity at the ADIY, ADIY1, EKIZ, ERGN, ERZR, and HINI GNSS stations. The station with the highest vertical velocity was ADIY1, measuring 4.92 mm/year, while the smallest vertical velocity was observed at station ERZR, measuring 0.03 mm/year. Conversely, the remaining stations exhibited a negative vertical velocity. The smallest negative vertical velocity was −0.03 mm/year at the SIRT GNSS station, while the UDE1 GNSS station had the highest vertical velocity of −7.10 mm/year. The vertical velocities of all GNSS stations are presented in Figure 9 and Table A1. As a result of the analyses, we obtained velocities for the height component of the stations. When examining the annual velocity of the height coordinate components, we observed a positive vertical velocity at the ADIY, ADIY1, EKIZ, ERGN, ERZR, and HINI GNSS stations. The station with the highest vertical velocity was ADIY1, measuring 4.92 mm/year, while the smallest vertical velocity was observed at station ERZR, measuring 0.03 mm/year. Conversely, the remaining stations exhibited a negative vertical velocity. The smallest negative vertical velocity was −0.03 mm/year at the SIRT GNSS station, while the UDE1 GNSS station had the highest vertical velocity of −7.10 mm/year. The vertical velocities of all GNSS stations are presented in Figure 9 and Table A1.

Conclusions
This paper aimed to reveal the relationship between the height component of the GNSS stations that make continuous observations and seasonal variation and meteorological parameters. For the vertical displacement analysis of the region, we selected 37 GNSS stations and processed approximately 6 years of continuous data.
First, we organized the data and used MATLAB to conduct an outlier test. This test helped identify and handle data points that significantly deviated from expected patterns or exhibited unusual behavior. After editing the data, we plotted the time series and investigated anomalies. We examined the seasonal behavior of GNSS stations and found strong seasonal signals at all stations. Additionally, we observed that GNSS signals exhibited noise. To perform spectral analysis, we utilized CATS software in this study. The most appropriate noise model for the stations was the WN+FN noise model, as expected. We then investigated the relationship between seasonal variation and the vertical time series. Based on our analyses, we calculated the seasonal and annual vertical velocities for each station, along with their respective standard deviations. The results showed variations among stations with different ellipsoid heights. As anticipated, the velocity and standard deviation values improved with longer observation periods.
Furthermore, we observed seasonal effects in the height components of continuous GNSS stations, prompting us to investigate their relationship with meteorological parameters. A simple linear regression analysis was conducted to assess this dependency, examining the relationship between the height components of continuous GNSS stations and meteorological factors such as temperature, pressure, relative humidity, wind speed, and precipitation. The analysis revealed a significant correlation between the height components and these meteorological parameters. When the results were examined, it can be seen that the model's performance increased for all parameters as the height of the ellipsoid increased. When the R-squared values approached one, and the MAE and RMSE values approached zero, the model performed successfully in all parameters. In addition to simple linear regression, we employed ARMA models. The ARMA modeling results further confirmed the dependence of the height components on meteorological parameters. The ARMA model results generally aligned well with the time series of all stations, with RMS values typically below 5 mm in the vertical direction, except for the ARTV, ERZ2, GIRS, MUUS, and SIRT1 stations.
The results from our study, including the simple linear regression analysis and ARMA modeling, provide strong evidence supporting the association between the height components of continuous GNSS stations and various meteorological parameters. Based on the analyses above, significant improvements in the quality of velocity estimates for continuous GNSS stations, particularly in the vertical component, can be achieved by either working with meteorological parameters or selecting the appropriate season (time) for conducting measurements based on the study's objective.  Data Availability Statement: All the results obtained at the end of our research have been thoroughly presented in the paper. The data presented in this study are available on request from the corresponding author.