Global Navigation Satellite System Time Series Research Articles

GCMEbox: A MATLAB toolbox for extracting and analyzing common-mode errors from GNSS time series

High-precision global navigation satellite system (GNSS) time series data contain rich tectonic and non-tectonic motion information. However, such data also have common-mode errors (CME) related to spatiotemporal characteristics. CME will produce bias in accurate and reliable GNSS station positions and velocities estimations, ultimately deciphering erroneous geophysical signals. Currently, CME extraction methods can be divided into three categories: regional stack filtering, reference frame filtering, and statistical signal decomposition. The third method is widely used in geodetic data processing, including GNSS time series analysis, global time-varying gravity signal extraction, interferometric synthetic aperture radar (InSAR) time series analysis, etc. The demand for GNSS time series analysis software with interactive, visualization, and multifunctional features has also increased with the increase in GNSS stations. Therefore, we developed gCMEbox based on MATLAB. The gCMEbox’s primary functions include GNSS time series preprocessing, multiple statistical signal decomposition algorithms for extracting CME, GNSS time series analysis, visualization, velocity field construction, and generation of multiple time series post-processing products. Based on GNSS time series data, this study utilized the gCMEbox for data preprocessing, including format conversion, outlier removal, missing data recovery, etc. By applying statistical signal decomposition methods, gCMEbox successfully extracted CME. The primary objective of this study is to share this toolbox with the scientific community, providing a comprehensive tool for GNSS time series analysis and applications.

Evaluation of recent tectonic movement in northeast Japan by using long-term GNSS and tide gauge measurements

The study investigates tectonic movements in northeast Japan by using long-term (2000–2022) global navigation satellite system (GNSS) and tide gauge measurements. The effect of the 2011 Tohoku-Oki earthquake including the other eight seismic events that occurred within this period is also discussed using GNSS time-series. The result showed the break in GNSS-time series because of occurred earthquakes and pointed out tectonic movements significantly. The GNSS sites located in the central and southern parts of northeast Japan showed that the velocity vectors have strong internal variation and suggest the existence source of alternative deformation because of geological terranes within the region. A least square approach was used, and the trend of sea-level measurements was fitted with the straight line. The obtained results from tide gauge measurements showed a rising trend at almost every site and indicated lithospheric uprising movement because of tectonic activities. This is possible because of the ongoing subduction of the Pacific and Philippine Sea plates beneath the Eurasian and North American plates.

An improved Gaussian process for filling the missing data in GNSS position time series considering the influence of adjacent stations

Due to various unavoidable reasons or gross error elimination, missing data inevitably exist in global navigation satellite system (GNSS) position time series, which may result in many analysis methods not being applicable. Typically, interpolating the missing data is a crucial preprocessing step before analyzing the time series. The conventional methods for filling missing data do not consider the influence of adjacent stations. In this work, an improved Gaussian process (GP) approach is developed to fill the missing data of GNSS time series, in which the time series of adjacent stations are applied to construct impact factors, together with a comparison of the conventional GP and the commonly used cubic spline methods. For the simulation experiments, the root mean square error (RMSE), mean absolute error (MAE) and correlation coefficient (R) are adopted to evaluate the performance of the improved GP. The results show that the filled missing data of the improved GP are closer to the true values than those of the conventional GP and cubic spline methods, regardless of the missing percentages ranging from 5 to 30%, with an interval of 5%. Specifically, the mean relative RMSE and MAE improvements for the improved GP with respect to the conventional GP are 21.2%, 21.3% and 8.3% and 12.7%, 16.2% and 11.01% for the North (N), East (E) and Up (U) components, respectively. In the real experiment, eight GNSS stations are analyzed using improved GP, together with conventional GP and a cubic spline. The results indicate that the first three principal components (PCs) of the improved GP can perverse 98.3%, 99.8% and 77.0% of the total variance for the N, E and U components, respectively. This value is obviously higher than those of the conventional GP and cubic spline. Therefore, we can conclude that the improved GP can better fill in the missing data in GNSS position time series than the conventional GP and cubic spline because of the impacts of adjacent stations.

On the stochastic significance of peaks in the least-squares wavelet spectrogram and an application in GNSS time series analysis

In this paper, the mathematical derivation of the underlying probability distribution function for the normalized least-squares wavelet spectrogram is presented. The impact of empirical and statistical weights on the estimation of the spectral peaks and their significance are demonstrated from the statistical point of view both theoretically and practically. The simulation results show an improvement of approximately 0.02mm (RMSE) for annual signal estimation when statistical weights are considered in the least-squares wavelet analysis (LSWA). The weighted LSWA estimates the signals more accurately than the ordinary LSWA for different percentage amount of missing data. As a real-world application, Global Navigation Satellite Systems (GNSS) time series for a station in Rome, Italy are analyzed. The analyses of the GNSS time series provided by different agencies for the same station reveal statistically significant annual peaks, more significant in 2010 but less significant between 2018 and 2020, while the higher frequency components show different spectral patterns over time. A declining trend of approximately −0.42 mm/year since 2004 is estimated for the GNSS height time series, likely due to gradual land subsidence. The results not only highlight the advantages of LSWA but can also help to better understand the uncertainties involved in signal estimation.

Deep Learning CNN-GRU Method for GNSS Deformation Monitoring Prediction

Hydraulic structures are the key national infrastructures, whose safety and stability are crucial for socio-economic development. Global Navigation Satellite System (GNSS) technology, as a high-precision deformation monitoring method, is of great significance for the safety and stability of hydraulic structures. However, the GNSS time series exhibits characteristics such as high nonlinearity, spatiotemporal correlation, and noise interference, making it difficult to model for prediction. The Neural Networks (CNN) model has strong feature extraction capabilities and translation invariance. However, it remains sensitive to changes in the scale and position of the target and requires large amounts of data. The Gated Recurrent Units (GRU) model could improve the training effectiveness by introducing gate mechanisms, but its ability to model long-term dependencies is limited. This study proposes a combined model, using CNN to extract spatial features and GRU to capture temporal information, to achieve an accurate prediction. The experiment shows that the proposed CNN-GRU model has a better performance, with an improvement of approximately 45%, demonstrating higher accuracy and reliability in predictions for GNSS deformation monitoring. This provides a new feasible solution for the safety monitoring and early warning of hydraulic structures.

Selection of noise models for GNSS coordinate time series based on model averaging algorithm

In the field of global navigation satellite system (GNSS) time series noise analysis, appropriately modeling the noise components plays an important role in determining the velocity of GNSS sites and quantifying the uncertainty associated with the velocity estimation. Over the years, researchers have focused on only one optimal noise model, while other noise models that show similar performance to the optimal model have been ignored. We investigated whether these ignored noise models can be made use of to describe the noise in the GNSS time series after applying a model averaging algorithm. The experimental data were derived from 28 International GNSS Service (IGS) sites in the California region of the United States and 110 IGS sites worldwide. The results showed that for the GNSS time series of 28 IGS sites in the California, 79%, 68%, and 75% of the site components can be applied the model averaging algorithm in the east/north/up (E/N/U) directions, respectively. Based on it, the east direction showed the best performance, with 50% of the site components obtaining more conservative velocity uncertainty after applying the model averaging algorithm compared to the optimal noise model. For GNSS time series of 110 IGS stations worldwide, the model averaging algorithm demonstrates excellent performance in all the E/N/U directions. In the E/N/U directions, 86%, 94%, and 57% of the site components can apply the model averaging algorithm. Building upon this, 77%, 65%, and 62% of the site components achieve more conservative velocity uncertainty in the E/N/U directions compared to the optimal noise model. To fully validate the feasibility of the model averaging algorithm, we also tested GNSS time series of varying lengths and different thresholds of the model averaging algorithm. In summary, the model averaging algorithm performs exceptionally well in the noise analysis of GNSS time series. It helps prevent overly optimistic estimation results.

Extraction of periodic signals in Global Navigation Satellite System (GNSS) vertical coordinate time series using the adaptive ensemble empirical modal decomposition method

Abstract. Empirical modal decomposition (EMD) is an efficient tool for extracting a signal from stationary or non-stationary time series and is enhanced in stability and robustness by ensemble empirical mode decomposition (EEMD). Adaptive EEMD further improves computational efficiency through adaptability in the white noise amplitude and set average number. However, its effectiveness in the periodic signal extraction in Global Navigation Satellite System (GNSS) coordinate time series regarding the inevitable missing data and offset issues has not been comprehensively validated. In order to thoroughly investigate their impacts, we simulated 5 years of daily time series data with different missing data percentages or a different number of offsets and conducted them 300 times for each simulation. The results show that high accuracy could reach the overall random missing rate below 15 % and avoid consecutive misses exceeding 30 d. Meanwhile, offsets should be corrected in advance regardless of their magnitudes. The analysis of the vertical components of 13 stations within the Australian Global Sea Level Observing System (GLOSS) monitoring network demonstrates the advantage of adaptive EEMD in revealing the time-varying characteristics of periodic signals. From the perspectives of correlation coefficients (CCs), root mean square error (RMSE), power spectral density indices (κ) and signal-to-noise ratio (SNR), the means for adaptive EEMD are 0.36, 0.81, −0.18 and 0.48, respectively, while for least squares (LS), they are 0.27, 0.86, −0.50 and 0.23. Meanwhile, a significance test of the residuals further substantiates the effectiveness in periodic signal extraction, which shows that there is no annual signal remaining. Also, the longer the series, the higher the accuracy of the reasonable extracted periodic signal concluded via the significance test. Moreover, driving factors are more effectively facilitated by the time-varying periodic characteristics compared with the constant periodic signal derived by LS. Overall, the application of adaptive EEMD could achieve high accuracy in analyzing GNSS time series, but it should be based on properly dealing with missing data and offsets.

A machine learning approach for slow slip event detection using GNSS time-series

Extracting tectonic transient displacements on the Earth’s surface from Global Navigation Satellite System (GNSS) time series remains a challenge, because GNSS station displacements depend on multiple processes occurring simultaneously, along with noise that obscures low-magnitude transient signals. We present a novel method for automatic detection of slow slip events (SSEs) in time series of a GNSS network by training a supervised machine learning (ML) model for classification. The proposed methodology detects both temporally and spatially the signatures of SSEs or regional transients within a GNSS network. The time series of a GNSS network were transformed into grayscale images, from which descriptors, including Bag of Visual Words (BoW) and Extended Local Binary Patterns (ELBP), were extracted. These descriptors served as input features for two distinct ML models: Support Vector Machines (SVM) and Artificial Neural Networks (NN). To train and test the ML classification model, two 3-year synthetic datasets were generated, one with GNSS networks featuring slow slip events (SSEs) of varying location, duration, onset time, and magnitude, and the other without SSEs, resulting in positive and negative sets, respectively. For each GNSS network, an image was created by combining the east and north components of the time series, which have been previously detrended and common mode error filtered. Each image is further divided into sub-images corresponding to 60 days time windows, in order to temporarily detect the existence of a transient. For training and testing, the datasets were separated into 75% for training and 25% for testing, each with 50% positive and 50% negative cases. In the final step, we analyze the positively classified images, representing the time windows in which the classifier detected transients. Within each of these windows, we identify the network’s time series with the highest velocity, indicating the stations and geographic area where the detected transients occurred. The test results demonstrate that both ML models achieved high performance using both ELBP and BoW descriptors as features. Finally, our ML models were validated on a real dataset with a transient signal recorded before the 2014 Iquique earthquake in Chile, and they effectively detected this anomalous signal. The proposed method can effectively detect transient signals related to SSEs with high accuracy, sensitivity, and specificity in both the test and instrumentally recorded datasets.

Joint Bayesian Modeling of Velocity Break Points, Noise Characteristics, and Their Uncertainties in GNSS Time Series: Far‐Field Velocity Anomalies Concurrent With Magmatic Activity in Iceland

AbstractEstimating the timing of velocity changes (break points) in Global Navigation Satellite System (GNSS) coordinate time series is required for understanding various Earth processes and how they may couple with each other. We jointly estimate break points, noise parameters and their uncertainties in GNSS time series using Bayesian interference. Synthetic data experiments demonstrate that time‐correlated noise can cause spurious estimates of break points for small velocity change and suggest the lower limits of velocity change for reliable estimates. A case study at the Krafla volcanic system shows increased rift‐perpendicular movementof 7.6–9.8 mm/yr, preceding by 23–77 days the 2014 Bárðarbunga rifting episode that occurred ∼120 km to the south (maximum likelihood estimates). These far‐field velocity anomalies occur in the same period as enhanced near‐field seismicity and deformation before the dike intrusion at Bárðarbunga, possibly indicating the coupling between the two volcanic systems through a deep partial melt zone.

Periodic signal extraction of GNSS height time series based on adaptive singular spectrum analysis

Singular spectrum analysis is widely used in geodetic time series analysis. However, when extracting time-varying periodic signals from a large number of Global Navigation Satellite System (GNSS) time series, the selection of appropriate embedding window size and principal components makes this method cumbersome and inefficient. To improve the efficiency and accuracy of singular spectrum analysis, this paper proposes an adaptive singular spectrum analysis method by combining spectrum analysis with a new trace matrix. The running time and correlation analysis indicate that the proposed method can adaptively set the embedding window size to extract the time-varying periodic signals from GNSS time series, and the extraction efficiency of a single time series is six times that of singular spectrum analysis. The method is also accurate and more suitable for time-varying periodic signal analysis of global GNSS sites.

Function Model Based on Nonlinear Transient Rheology of Rocks: An Analysis of Decadal GNSS Time Series After the 2011 Tohoku‐oki Earthquake

AbstractThe postseismic Global Navigation Satellite System (GNSS) time series after the 2011 Tohoku‐oki earthquake is often empirically explained by one or more logarithmic and exponential decay functions. However, a function model with these decay functions struggles to illuminate various deformation mechanisms. Here, we propose a new function model incorporating laboratory‐derived constitutive laws (power‐law flow and velocity‐strengthening friction) utilized in mechanical postseismic models of the last decade. We demonstrated that our function model accurately predicts the decade‐long time series of inland GNSS stations, including coastal areas where significant postseismic uplift continues now. Moreover, it decomposes them into displacements due to viscoelastic relaxation and afterslip, similar to results produced by previous stress‐dependent postseismic models of the 2011 Tohoku‐oki earthquake. This physics‐based function model is an effective tool to ease the process of forecasting long‐term GNSS time series by dividing them into two dominant mechanisms at the plate interface and the surrounding viscoelastic mantle.

Wavelet-like denoising of GNSS data through machine learning. Application to the time series of the Campi Flegrei volcanic area (Southern Italy)

The great potential of the Global Navigation Satellite System (GNSS) in monitoring ground deformation is widely recognized. As with other geophysical data, GNSS time series can be significantly noisy, hiding elusive ground deformation signals. Several denoising techniques have been proposed to improve the signal-to-noise ratio over the years. One of the most effective denoising techniques has been proved to be multi-resolution decomposition through the discrete wavelet transform. However, wavelet analysis requires long data sets to be effective, as well as long computation times, that hinder its use as a real or near real-time monitoring tool. We propose training by a Convolutional Neural Network (CNN) to perform the equivalent of wavelet analysis to overcome these limitations. Once trained, the CNN model provides answers within seconds, making it feasible as a real-time data analysis tool. Our Machine Learning algorithm is tested on daily GNSS time series collected in the Campi Flegrei caldera (Southern Italy), which is a highly volcanic risk area. Without significant gaps, the retrieved RMSE and R2 values vary in the ranges 0.65–0.98 and 0.06–0.52 cm, respectively. These results are encouraging, as they hint at the possibility of applying this methodology in more effective real-time monitoring solutions for active volcanoes.

A Temporal and Spatial Denoising Method for Intelligent Settlement Sensing System

Global Navigation Satellite System (GNSS) has the ability to provide timely structural settlement information and been widely installed on the critical infrastructures. However, the inevitable noise negatively affects the accurate assessment, measurement and evaluation of the status. Although some methods have been developed in previous studies, fully mining the temporal and spatial correlations is still necessary to further analyze. Therefore, we propose a new denoising method based on truncated high-order singular value decomposition to reduce the noise in multivariate GNSS signals. Using synthetic signals and real-world GNSS signals, the proposed method is evaluated by comparison with some benchmark models. The results show that the proposed method can effectively reduce both white noise (WN) and flicker noise (FN). Moreover, experiments on the signal-noise ratio (SNR) from 1 to 10 demonstrate that the proposed method can achieve stable performances for GNSS time series denoising. The organization of GNSS signals as the high-dimensional tensor has been proved to be an effective analytical tool to mine the complicated spatial and temporal correlations.

Geodetic evidence of land subsidence in Cirebon, Indonesia

The subsidence phenomena occur in numerous cities in Indonesia (Jakarta, Bandung, and Semarang) and have been observed over the last few decades using geodetic techniques such as interferometric synthetic aperture radar (InSAR) and the Global Navigation Satellite System (GNSS). Three variables are considered the leading causes of substantial subsidence in those cities: urban development, excessive groundwater removal, and geological context. Considering these considerations, we assess surface displacement in the surrounding area of Cirebon, Indonesia, which has similar background circumstances. We used 43 Sentinel-1A images covering from 2014 to 2021 and processed them using the small baseline subset (SBAS) method obtaining 78 interferograms for time series analysis. A weighted correction from two weather models was used to correct tropospheric errors. The results from the InSAR time series were compared with the continuous GNSS time series (2010–2021) processed using GAMIT/GLOBK GNSS processing softwares. Cirebon experienced a slow land subsidence rate with an average magnitude of -2.6 mm/yr. Furthermore, we underline that more significant subsidence rates were found in several locations. Conic-shaped subsidence with a diameter of around 1 km was observed in a highly inhabited area, with a maximum rate of up to -17 mm/yr near the cone’s center. We also found more prominent land subsidence along the coastline with magnitudes up to -25 mm/yr. Surface lowering with a magnitude of up to -32 mm/yr also affects power plant regions close to Cirebon city. Meanwhile, the fastest subsidence rate was observed over the salt evaporation field area of about -50 mm/yr. We attribute the natural and anthropogenic factors to causing the subsidence in this area.

GNSS_Vel_95CI.py: A Python Module for Calculating the Uncertainty of GNSS-Derived Site Velocity

GNSS_Vel_95CI.py is an open-source Python-3 module for calculating the 95% confidence interval (95% CI) for the site velocity derived from global navigation satellite systems (GNSS) daily positions, which are often affected by time-correlated noises. The detailed methodology for calculating the 95% CI is documented in a recent article (Wang 2022) and initially programmed in Fortran. However, few young researchers (e.g., graduate students) are familiar with Fortran now. In an effort to support a broader user community, we have realized the method in the Python programming language, which has been commonly taught in current college curriculums. Through the use of this module, researchers and engineers can focus on the applications of GNSS time series, rather than on coding and data processing. In particular, the module has the option to output the autocorrelation function (ACF) and the GNSS time-series decomposition results: the linear, nonlinear, seasonal, and residual components, which allow students and researchers having little coding experience to conduct advanced GNSS time-series analysis. The module is versatile, easy to install through the pip installer and GitHub, and simple to use. An example Python program is provided to illustrate the use of this module.

Machine Learning Models Applied to a GNSS Sensor Network for Automated Bridge Anomaly Detection

Structural health monitoring (SHM) based on global navigation satellite systems (GNSS) is an interesting solution to provide absolute positions at different locations of a structure in a global reference frame. In particular, low-cost GNSS stations for large-scale bridge monitoring have gained increasing attention these last years because recent experiments showed the ability to achieve a subcentimeter accuracy for continuous monitoring with adequate combinations of antennas and receivers. Technical solutions now allow displacement monitoring of long bridges with a cost-effective deployment of GNSS sensing networks. In particular, the redundancy of observations within the GNSS network with various levels of correlations between the GNSS time series makes such monitoring solution a good candidate for anomaly detection based on machine learning models, using several predictive models for each sensor (based on environmental conditions, or other sensors as input data). This strategy is investigated in this paper based on GNSS time series, and an anomaly indicator is proposed to detect and locate anomalous structural behavior. The proposed concepts are applied to a cable-stayed bridge for illustration, and the comparison between multiple tools highlights recurrent neural networks (RNN) as an effective regression tool. Coupling this tool with the proposed anomaly detection strategy enables one to identify and localize both real and simulated anomalies in the considered data set.

Water-Vapour Monitoring from Ground-Based GNSS Observations in Northwestern Argentina

The Central Andes in northwestern Argentina are characterized by steep topographic and climatic gradients. The humid foreland areas at 1 km asl elevation rapidly rise to over 5 km in the eastern Cordillera, and they form an orographic rainfall barrier on the eastern windward side. This topographic setting combined with seasonal moisture transport through the South American monsoon system leads to intense rainstorms with cascading effects such as landsliding and flooding. In order to better quantify the dynamics of water vapour transport, we use high-temporal-resolution global navigation satellite system (GNSS) remote sensing techniques. We are particularly interested in better understanding the dynamics of high-magnitude storms with high water vapour amounts that have destructive effects on human infrastructure. We used an existing GNSS station network with 12 years of time series data, and we installed two new ground stations along the climatic gradient and collected GNSS time series data for three years. For several stations we calculated the GNSS signal delay gradient to determine water vapour transport direction. Our statistical analysis combines in situ rainfall measurements and ERA5 reanalysis data to reveal the water vapour transport mechanism for the study area. The results show a strong relationship between altitude and the water vapour content, as well as between the transportation pathways and the topography.

Iteration empirical mode decomposition method for filling the missing data of GNSS position time series

The Global Navigation Satellite System (GNSS) can provide the daily position time series for the geodesy and geophysical studies. However, due to various unpredictable factors, such as receiver failure or bad observation conditions, missing data inevitably exist in GNSS position time series. Most traditional time series analysis methods require the time series should be completed. Therefore, filling the missing data is a valuable step before analyzing the GNSS time series. In this study, a new method named Iteration Empirical Mode Decomposition (Iteration EMD) is proposed to fill the missing data in GNSS position time series. The simulation experiments are performed by randomly removing different missing percentages of the synthetic time series, with the added different types noise. The results show that Iteration EMD approach performs well regardless of high or low missing percentage. When the missing percentage increases from 5 % to 30 % with a step of 5 %, all the Root Mean Square Errors (RMSE) and Mean Absolute Errors (MAE) of Iteration EMD are smaller than Interpolation EMD. The relative improvements at different percentages of Iteration EMD relative to Interpolation EMD are significant, especially for the high missing percentage. The real GNSS position time series of eight stations were selected to further evaluate the performance of Iteration EMD with an average missing percentage 8.15 %. Principal Component Analysis (PCA) was performed on the filled time series, which is used to assess the interpolation performance of Iteration EMD and Interpolation EMD. The results show that Iteration EMD can preserve variance 75.9 % with the first three Principal Components (PC), more than 66.5% of interpolation EMD. Therefore, we can conclude that Iteration EMD is an efficient interpolation method for GNSS position time series, which can make full use of available information in existing time series to fill the missing data.

Modelling and prediction of GNSS time series using GBDT, LSTM and SVM machine learning approaches

Global navigation satellite system (GNSS) site coordinate time series provides essential data for geodynamic and geophysical studies, realisation of a regional or global geodetic reference frames, and crustal deformation research. The coordinate time series has been conventionally modelled by least squares (LS) fitting with harmonic functions, alongside many other analysis methods. As a key limitation, the traditional modelling approaches simply use the functions of time variable, despite good knowledge of various underlying physical mechanisms responsible for the site displacements. This paper examines the use of machine learning (ML) models to reflect the effects or residential effects of physical variables related to Sun and the Moon ephemerides, polar motion, temperature, atmospheric pressure, and hydrology on the site displacements. To form the ML problem, these variables are constructed as the input vector of each ML training sample, while the vertical displacement of a GNSS site is regarded as the output value. In the evaluation experiments, three ML approaches, namely the gradient boosting decision tree (GBDT) approach, long short-term memory (LSTM) approach, and support vector machine (SVM) approach, are introduced and evaluated with the time series datasets collected from 9 GNSS sites over the period of 13 years. The results indicate that all three approaches achieve similar fitting precision in the range of 3–5 mm in the vertical displacement component, which is an improvement in over 30% with respect to the traditional LS fitting precision in the range of 4–7 mm. The prediction of the vertical time series with the three ML approaches shows the precision in the range of 4–7 mm over the future 24- month period. The results also indicate the relative importance of different physical features causing the displacements of each site. Overall, ML approaches demonstrate better performance and effectiveness in modelling and prediction of GNSS time series, thus impacting maintenance of geodetic reference frames, geodynamics, geophysics, and crustal deformation analysis.

Quantitative Evaluation of Environmental Loading Products and Thermal Expansion Effect for Correcting GNSS Vertical Coordinate Time Series in Taiwan

We explored the driving factors of nonlinear signals in vertical coordinate sequences of stations in a Taiwan global navigation satellite system (GNSS) network, including atmospheric loading (ATML), hydrological loading (HYDL), and non-tidal ocean loading (NTOL) effects. At the same time, we used the finite element analysis software MIDAS to quantify the vertical displacements of different types of monuments due to the thermal expansion effect, including deep drilled braced (DDB) and short drilled braced (SDB). By quantitatively comparing the correction results of GNSS time series with different single mass loading models, we found that there was little difference in the correction of different environmental loading products. We compared different combinations of each loading product to correct the GNSS time series, and finally selected the best combination suitable for Taiwan GNSS network, that is, ATML (GFZ_ECMWF IB) + HYDL (IMLS_MERRA2) + NTOL (IMLS_MPIOM06). We found that the spatial and temporal models of ATML and NTOL are very similar, with non-tidal atmospheric loading and non-tidal ocean loading working together, a pattern that may be related to tropical cyclones. Both models also showed good correction effect on GNSS stations in the western plain of Taiwan, but with limited correction effect in the eastern part of Taiwan. This may be due to the influence of the subtropical monsoon climate in Taiwan and the barrier of the central mountain range, resulting in obvious differences between eastern and western Taiwan. The hydrological loading was found to act in the opposite way to the thermal expansion effect in the temporal domain, indicating that some displacements in hydrological loading may cancel out displacements caused by the thermal expansion effect. This aspect of displacement is not included in the hydrological loading model but should be considered when accurately estimating the temporal and spatial variation of water storage capacity in Taiwan using GNSS observed displacements.