A Comprehensive Analysis of Environmental Loading Effects on Vertical GPS Time Series in Yunnan, Southwest China

Abstract

Seasonal variations in the vertical Global Positioning System (GPS) time series are mainly caused by environmental loading, e.g., hydrological loading (HYDL), atmospheric loading (ATML), and nontidal oceanic loading (NTOL), which can be synthesized based on models developed by various institutions. A comprehensive comparison among these models is essential to extract reliable vertical deformation data, especially on a regional scale. In this study, we selected 4 HYDL, 5 ATML, 2 NTOL, and their 40 combined products to investigate their effects on seasonal variations in vertical GPS time series at 27 GPS stations in Yunnan, southwest China. These products were provided by the German Research Center for Geosciences (GFZ), School and Observatory of Earth Sciences (EOST), and International Mass Loading Service (IMLS). Furthermore, we used the Cross Wavelet Transform (XWT) method to analyze the relative phase relationship between the GPS and the environmental loading time series. Our result showed that the largest average Root-Mean-Square (RMS) reduction value was 1.32 mm after removing the deformation associated with 4 HYDL from the vertical GPS time series, whereas the RMS reductions after 5 ATML and 2 NTOL model corrections were negative at most stations in Yunnan. The average RMS reduction value of the optimal combination of environmental loading products was 1.24 mm, which was worse than the HYDL (IMLS_GEOSFPIT)-only correction, indicating that HYDL was the main factor responding for seasonal variations at most stations in Yunnan. The XWT result showed that HYDL also explained the annual variations reasonably. Our finding implies that HYDL (IMLS_GEOSFPIT) contributes the most to the environmental loading in Yunnan, and that the ATML and NTOL models used in this paper cannot be effective to correct seasonal variations.

1. Introduction

Vertical GPS time series not only contain tectonic signals (e.g., volcanic motion, plate motion) [1,2,3], which present long-term linear motion, but they also record nonlinear variation signals, mainly annual and semi-annual seasonal variations, which have been demonstrated to be affected by environmental loading [4,5,6], e.g., hydrological loading (HYDL), atmospheric loading (ATML), and nontidal oceanic loading (NTOL). At present, numerous HYDL, ATML, and NTOL products are provided by several institutes, including the German Research Center for Geosciences (GFZ), School and Observatory of Earth Sciences (EOST), and International Mass Loading Service (IMLS), which have been widely used to investigate their effects on seasonal variations in vertical GPS time series [7,8,9,10]. Seasonal variation signals are often regarded as noise in the geodetic data used for tectonic research. Thus, an analysis of the HYDL, ATML, and NTOL effects on seasonal variations in GPS data is important for obtaining precise vertical GPS velocity.

Seasonal variations caused by HYDL, ATML, and NTOL may differ within the same region. Li et al. [11] investigated HYDL, ATML, and NTOL effects on the seasonal variations at Crustal Movement Observation Network of China (CMONOC) stations, and the median RMS reduction rates corrected by the optimal HYDL, ATML, and NTOL models in GPS time series were 6.56%, 17.17%, and 1.44%, respectively. Yuan et al. [12] used 235 GPS stations within 2~16.5 years from the CMONOC to study the environmental loading effects on the seasonal variations in vertical GPS time series. The results showed that the average RMS caused by HYDL, ATML, NTOL, and their combined effect was 2.7, 0.6, 3.2, and 4 mm, respectively. In addition, the HYDL, ATML, and NTOL time series in the same region provided by various institutions may differ due to different geophysical models and methods [13], such as the Green function and spherical harmonics. For example, Wu et al. [14] used the HYDL products provided by GFZ and EOST to study their seasonal variation effects at 633 global stations. The RMS reduction rates were positive at 82.6 and 87.4%, respectively. Andrei et al. [15] compared the time series at ABOA station with environmental loading time series and found that the difference in NTOL provided by GFZ and EOST reached 3 mm. Li et al. [11] used five HYDL, six ATML, and five NTOL products provided by GFZ, EOST, and IMLS, respectively, to quantitatively evaluate their effects on seasonal variations at CMONOC stations. The results showed that the difference in RMS reduction reached 20% based on an optimal combination of HYDL, ATML, and NTOL products. Therefore, it is very important to choose suitable environmental loading products provided by various institutions for correcting seasonal variations in vertical GPS time series in the same region.

Yunnan is situated in the southeast of the Tibet Plateau, an area that has tectonic activity and serious earthquakes [16]. Therefore, investigating the environmental loading effects on GPS is of great importance for obtaining a reliable velocity field that can contribute to research into the tectonic evolution of the Tibet Plateau. At present, many scholars are using environmental loading to study their effects on vertical GPS time series in Yunnan [17,18,19], and they concluded that HYDL is the main factor causing seasonal variations in this region; however, the seasonal variations corrected by various HYDL products may differ. In addition, the difference between the various HYDL, ATML, and NTOL products for correcting the seasonal variations have not been comprehensively analyzed, but the effect of ATML is an indisputable contributor. Tan et al. [20] used independent component analysis to separate the different sources of deformation from the 22 GPS stations in Yunnan, and the results showed that annual variations in the first independent component may have been caused by soil moisture and atmospheric loading. Zhang et al. [21] investigated the environmental loading effects on the Common Mode Error (CME) in GPS data in the northwest of Yunnan. Compared with HYDL, the combination of the HYDL, ATML, and NTOL time series provided by GFZ had better consistency with the CME in vertical GPS time series. However, Liu et al. [22] used GPS data to estimate Terrestrial Water Storage (TWS) variations, in which the ATML and NTOL time series provided by IMLS were removed from the GPS. The TWS variations were anomalies, which may have been caused by inaccurate ATML correction. Zhan et al. [23] investigated the vertical annual motion, and motion obtained by the ATML time series provided by EOST was quite different from the GPS in amplitude and phase; therefore, the ATML is not the main factor causing the annual motion.

Based on the above studies, we comprehensively analyzed different HYDL, ATML, NTOL, and their combined effects on vertical GPS time series in Yunnan. We used the vertical GPS time series span from January 2011 to August 2019 for 27 stations in Yunnan and compared them with 4 HYDL, 5 ATML, 2 NTOL, and their 40 combined products provided by GFZ, EOST, and IMLS. The correlation and RMS reduction analysis were only used in the time domain to evaluate the environmental loading effects on seasonal variations. Furthermore, the XWT method was used to analyze the relative phase relationship between vertical GPS time series and environmental loading deformation in time–frequency space.

This paper is organized as follows: Section 2 describes the data and methods, including the GPS and environmental loading model and the RMS reduction, correlation, and XWT methods. Section 3 reports the results of HYDL, ATML, NTOL, and their combined effects on the GPS and the XWT analysis. Section 4 discusses the influence of the noise model, velocity, and the velocity uncertainty due to environmental loading in vertical GPS time series and its effects on the seasonal variations in Yunnan. Section 5 summarizes the main finding.

2. Data and Methods

2.1. Data

2.1.1. GPS Data

The vertical GPS time series of 27 stations in Yunnan from the CMONOC spanning from January 2011 to August 2019 was used in this study. The location distribution of 27 GPS stations is shown in Figure 1. The vertical GPS time series was calculated using GAMIT/GLOBK software. First, GAMIT solved the daily solution, which estimated station location, satellite orbit parameters, and the variance-covariance matrix [24]. Then, GLOBK adjusted the network to obtain the daily coordinate time series. The vertical GPS time series of 27 GPS stations was obtained by ftp://ftp.cgps.ac.cn/products/position/ (accessed on 22 May 2022).

To obtain a high-precision vertical GPS time series for further analysis, the vertical GPS time series needed to be preprocessed [25] (e.g., outliers removed, offsets correction, missing data interpolation). The outliers were removed using the interquartile range method, the offsets were corrected by Least Squares Fitting (LSF), and the interpolation of missing data was performed via the Regularized Expectation Maximization (RegEM) method [26]. It has been widely used to interpolate missing data in GPS data [27,28]. The RegEM method takes into account the correlations and objects between all stations and relies on data models without requiring prior information. The results of RegEM interpolation at some GPS stations (YNMH, YNLA, YNXP, YNCX) are shown in Figure 2. It can be seen that the result was consistent with the overall motion in the vertical GPS time series. To better analyze HYDL, ATML, NTOL, and their combined effects on seasonal variations, the linear trends were removed from the vertical GPS time series by the LSF method.

2.1.2. Environmental Loading Model

The mass redistribution of HYDL, ATML, and NTOL can result in seasonal deformation in the Earth’s surface, which can be detected by GPS data [29]. Thus, the environmental loading products can be used to evaluate the effects on seasonal variations in vertical GPS time series. GFZ provides the surface deformation on global grid data [30], which were calculated by the Land Surface Discharge Model (LSDM) [31], European Centre for Medium-Range Weather Forecasts (ECMWF), and the Max-Planck Institute Ocean Model (MPIOM) [32], respectively. The HYDL, ATML, and NTOL time series of the station were obtained through bicubic interpolation based on the global grid data with latitude and longitude. The HYDL and ATML products provided by EOST were only used in this work and were calculated by the ERA_interim model [33], which was estimated from the datasets of ECMWF reanalysis, and another ATML product was calculated by the ECMWF(IB) model. The HYDL, ATML, and NTOL products provided by IMLS were calculated by the Global Earth Observing System Forward Processing Instrumental Team (GEOSFPIT), Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA2) [34], and MPIOM06. To unite the time resolution with GPS data, all environmental loading time series under the center of the figure frame were averaged into daily time series. More detailed information about the HYDL, ATML, and NTOL products provided by GFZ, IMLS, and EOST are shown in Table 1.

The 4 HYDL (GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim), 5 ATML (GFZ_ECMWF, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim, EOST_ECMWF(IB)), 2 NTOL (GFZ_EMPIOM and IMLS_MPIOM06) products were used in this paper.

Figure 3 compares the time series induced by environmental loading provided by GFZ, IMLS, and EOST with vertical GPS time series at YNCX station. As shown in Figure 3, the HYDL, ATML, and GPS time series had a significant seasonal variation. The range of time series at YNCX station was −20~20 mm, while the range of HYDL, ATML, and NTOL time series was −10~10, −5~5, and −1~1 mm, respectively, indicating that environmental loading cannot fully explain the seasonal variations at YNCX station.

2.2. Methods

2.2.1. Quantitative Evaluation Metrics

We used the correlation and RMS reduction [35] to quantitatively evaluate HYDL, ATML, and, NTOL, and their combined effects on seasonal variations in vertical GPS time series. The correlation coefficient was calculated:

(1)$r=\frac{{\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}}\sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}$

where $x_i,\stackrel{\_}{x},y_i,\stackrel{\_}{y}$ are the daily and average values of the GPS and environmental loading time series, respectively, and $n$ is the number of observation epochs for each station. The RMS reduction was calculated:

(2)$DRMS=RMS(GPS)-RMS(GPS-Loading)$

where $RMS(GPS),RMS(GPS-Loading)$ are the $RMS$ values of the vertical GPS time series before and after HYDL, ATML, and NTOL or their combined correction.

2.2.2. Cross Wavelet Transform

XWT combines wavelet analysis and the cross spectrum, which has been widely used to analyze the phase relationship between two time series in time–frequency space [36]. The XWT of two time series $X_n$ and $Y_n$ is defined as:

(3)$W^{XY}=W^X{W}^{Y*}$

where $Y^{*}$ refers to the complex conjugation of $W^{Y*}$. The circular mean of the phase can be adopted to measure the phase relationship between $X_n$ and $Y_n$. The circular mean of a set of angles $\left({\alpha }_i,i=1,2,3,\dots n\right)$ was calculated by the following equation:

(4)$a_m=\arg \left(X,Y\right)=\arg \left[\frac{1}{n}{\displaystyle \sum _{i=1}^n\cos \left({\alpha }_i\right),\frac{1}{n}{\displaystyle \sum _{i=1}^n\sin \left({\alpha }_i\right)}}\right]$

The circular standard deviation $s=\sqrt{-2\ln \left(R/n\right)}\left(R=\sqrt{X^2+Y^2}\right)$ can be used as an index to measure the scatter of angles around the mean. Further processing of $a_m$ can show the correlation of $X_n$ and $Y_n$ directly. The XWT-based semblance [37] was calculated:

(5)${\rho }_i=\cos \left({\alpha }_m\right)$

where the value of $\rho$ ranges from −1 (completely inverse correlation) to 1 (completely positive correlation).

3. Results

3.1. Quantitative Assessment of Environmental Loading Effects

3.1.1. HYDL Effects

Figure 4 shows the correlation between HYDL products (GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim) and vertical GPS time series for all GPS stations. The correlation value of YNGM station was lowest. The range of correlation value between GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim, and vertical GPS time series was 0.02~0.7, 0.1~0.74, 0.09~0.73, and −0.07~0.71 with an average value of 0.51, 0.54, 0.55, and 0.49, respectively, indicating that the time series derived from HYDL products correlated with the vertical GPS time series at most stations.

As shown in Figure 5, the RMS reduction at the YNGM station was negative after 4 HYDL products were corrected, while other stations were positive. The range of RMS reduction from the vertical GPS time series corrected by GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, and EOST_ERA_interim was −2.19~2.18, −0.64~2.38, −0.88~2.56, and −1.32~1.95 mm with an average value of 0.79, 1.25, 1.32, and 1.02 mm, respectively. The difference between GFZ_LSDM and IMLS_GEOSFPIT reached 0.53 mm, and the largest average RMS reduction value was 1.32 mm, corresponding to the IMLS_GEOSFPIT product, which implied that HYDL has a significant effect on the seasonal variation in the GPS time series. This result is consistent with those of Sheng et al. [18] and Hao et al. [17]. According to the correlation and RMS reduction among the 4 HYDL products correction, the HYDL (IMLS_GEOSFPIT) product had the best positive effect on correcting the seasonal variations.

From Figure 5, the RMS reduction value of GPS stations (YNHZ, YNDC, KMIN, YNTH, YNML, YNMZ, YNWS, YNJP) near the Xiaojiang fault in eastern Yunnan was smaller than for the other stations after IMLS_GEOSFPIT correction. Rainfall in Yunnan is mainly influenced by the warm and humid air over the South China Sea and the Bay of Bengal from south to north [19], and the terrain gradually increases from south to north and from west to east. When the air reaches the northern and eastern parts of Yunnan, it has less water vapor, so precipitation decreases. Therefore, the TWS variations are not obvious near the Xiaojiang fault, and HYDL cannot adequately explain the seasonal variations in vertical GPS time series. The YNGM station is located in a region with obvious TWS variations, but the RMS reduction value was negative, which means that the anomalous of YNGM station may have been affected by unknown local effects, GPS error, and noise.

3.1.2. ATML Effects

As shown in Figure 6, the correlation range between the 5 ATML products and the vertical GPS time series for all GPS stations had a similar pattern. The time series caused by the ATML showed a lower correction with the vertical GPS time series.

The results of RMS reduction are shown in Figure 7. The RMS reduction value was negative at most GPS stations in Yunnan, the range of RMS reduction values of vertical GPS time series corrected by ATML products (GFZ_ECMWF, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim, EOST_ECMWF(IB)) was −0.68~0.34, −0.56~0.33, −0.61~0.35, −0.73~0.34, and −0.69~0.35 mm with an average value of −0.13, −0.07, −0.1, −0.17, and −0.14 mm, respectively. Above the correction and RMS reduction, the 5 ATML products could not correct the seasonal variations in vertical GPS time series.

3.1.3. NTOL Effects

The 2 NTOL products provided by GFZ and IMLS were only used in this research. As shown in Figure 8 and Figure 9, the value of correction and RMS reduction was relatively small. The range of correlation was −0.14~0.13 and −0.18~0.14, respectively. The range of RMS reduction was −0.07~0.04 and −0.09~0.05 mm with an average value of 0.004 and 0 mm, respectively. The 2 NTOL products could not effectively correct the seasonal variations in vertical GPS time series.

3.1.4. Optimal Combination of HYDL, ATML, and NTOL Effects

We calculated the total environmental loading time series induced by HYDL, ATML, and NTOL, provided by GFZ, EOST, and IMLS, respectively. It produced a combination of 4 HYDL, 5 ATML, and 2 NTOL products. There were 40 different combinations of environmental loading products in all, and the range of mean correlation and RMS reduction was 0.5~0.57 and 0.91~1.24 mm, respectively.

To further analyze the combined loading effects on seasonal variations in vertical GPS time series, the first 10 and last 10 combined results sorted by the value of mean RMS reduction was listed in Table 2. The optimal combination of environmental loading products was HYDL (GFZ_LSDM), ATML (EOST_ERA_interim), and NTOL (IMLS_MPIOM06), of which the mean correlation and RMS reduction value were 0.566 and 1.24 mm, respectively, whereas the worst combination was HYDL (IMLS_MERRA2), ATML (EOST_ERA_interim), and NTOL (EOST_ERA_interim). The difference in mean correlation and RMS reduction value between the optimal and worst combination was 0.07 and 0.34 mm, respectively. Among the 40 combined results, the combed loading products, including HYDL (GFZ_LSDM), caused more obvious seasonal variations in the vertical GPS time series, which was superior to the combined loading products, including HYDL (IMLS_MERRA2). From Section 3.1.1, the mean value of the RMS reduction corrected by HYDL (GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim) products was 0.79, 1.25, 1.32, and 1.02 mm, respectively. The HYDL (IMLS_MERRA2, IMLS_GEOSFPIT) products provided by IMLS were superior to other HYDL products for correcting the seasonal variations. However, the optimal combined loading products did not contain the HYDL product provided by IMLS. One possible explanation is that the combination of individual best products may have resulted in a worse mass conservation problem.

3.2. Cross Wavelet Transform Analysis

The correlation and RMS reduction were only used in the time domain to analyze the environmental loading effects on the seasonal variations in GPS data. To further analyze the relationship between the deformation induced by environmental loading and the vertical GPS time series, the XWT method was used because it can analyze the relative phase relationship between the GPS and HYDL, ATML, NTOL, as well as their combined time series in time–frequency space. From Section 3.1.1, the mean value of RMS reduction corrected by HYDL (IMLS_GEOSFPIT) was 1.32 mm, which was superior to the optimal combination of environmental loading products. Hence, the HYDL (IMLS_GEOSFPIT), ATML (IMLS_GEOSFPIT), NTOL (IMLS_GEOSFPIT), and their combined products were used for XWT analysis.

The cross wavelet power spectrum between the GPS and HYDL, ATML, NTOL, and their combined time series at YNCX station, respectively, are shown in Figure 10. The relative phase relationship between two time series is described by the arrow (with anti-phase pointing left and in-phase pointing right; the downward-pointing arrow represents the vertical GPS time series leading the environmental loading time series by 90°). The phase angle will vary little or be consistent [38] between the GPS and environmental loading time series, which are physically related. From Figure 10, the obvious common power was located in the period range of 256~512 days in the whole time span between the GPS and HYDL, ATML, NTOL, and their combined time series. In addition, the obvious common power was located in the period range of 128~256 days for the whole time span between the GPS and HYDL time series, and the obvious common power was located in the period range of 128~256 days for a smaller time interval between the GPS and ATML, NTOL, and their combined time series, respectively. Moreover, as shown in Figure 10a, the relative phasing angle of the period less than 128 days fluctuated violently, sub-seasonal peaks of the vertical GPS time series may be caused by some technique errors (e.g., unmodeled high-order ionospheric path delays error and tropospheric delay error) [38,39]. The XWT results of other stations were the same as the YNCX station. The fluctuation range of the relative phase angle was small in the period range of 256~512 days corresponding to the annual cycle, and the period range of 256~512 days was more obvious than the other period. Thus, the mean XWT-based semblance close to 1 year was calculated in this work, unveiling more effectively the phase relationship between the GPS and HYDL, ATML, NTOL, and their combined time series. The results are shown in Figure 11.

As shown in Figure 11, the mean values of the XWT-based semblance at YNCX station were 0.98, 0.31, 0.28, and 0.97, respectively, implying that the ATML and NTOL models used in this paper were not the main factor causing the annual variations, while annual variations were caused by HYDL. From Figure 11a–c, the range of mean XWT-based semblance between the GPS and HYDL time series was 0.9~1 at all GPS stations (except for KMIN and YNGM), and the value of KMIN and YNGM stations was 0.79 and 0.37, respectively. The mean XWT-based semblance between GPS and ATML and NTOL time series at YNGM station was 0.95 and 0.98, respectively, while the range of other stations having a similar range pattern was −0.53~0.35 and −0.59~0.39, respectively. From Figure 11d, the number of stations with a mean XWT-based semblance value between 0.9 and 1 was 15, implying that the annual variations changes in the vertical GPS time series were affected by ATML and NTOL at some stations. Above the results, the HYDL was the main factor causing the annual variations at most stations. Meanwhile, ATML and NTOL could only explain the annual variations at YNGM station.

To verify the reliability of the XWT approach, we also used the LSF method to calculate the annual phase between the GPS and HYDL, ATML, NTOL, and their combined time series, respectively. The difference in annual phase results is shown in Figure 12. The range between the GPS and HYDL, ATML, NTOL, and their combined time series was −7.0°~6.0°, −7.9°~4.6°, −7.8°~4.3°, and −6.9°~5.8° with an average value of −0.76°, −1.6°, −2.0°, and −0.81°, respectively, indicating that the XWT method can be adapted to analyze the relative phase relationship between the GPS and the deformation induced by environmental loading.

4. Discussion

4.1. Changes in the Vertical GPS Time Series Due to Environmental Loading

Environmental loading can affect the stochastic part in the vertical GPS time series. According to Section 3.1.1, HYDL can explain the seasonal variations at all stations (except for YNGM), and the mean value of RMS reduction corrected by HYDL (IMLS_GEOSFPIT) was superior to the optimal combined loading products. Therefore, we analyzed the noise characteristics, velocity, and velocity uncertainty estimated from vertical GPS time series (except for YNGM) after being corrected by the HYDL (IMLS_GEOSFPIT). Hector software [40] was used to investigate the noise characteristics before and after HYDL correction. The noise characteristics estimated from the vertical GPS time series may differ after HYDL correction. Therefore, we selected different noise models [41]: White Noise (WN), White Noise and Flicker Noise (WN + FN), White Noise and Power-law Noise (WN + PL), White Noise and Flicker Noise and Random Walk Noise (WN + FN + RW), White Noise and Gauss Markov Model Noise (WN+GGM). The results of noise characteristics before and after HYDL correction are shown in Table 3 and Figure 13.

As Table 3 and Figure 13 show, the optimal noise model at most stations was WN + FN. Among all the stations, the number in WN + FN before HYDL correction was 22, accounting for 85% of all stations, while the number of stations in WN + FN after HYDL correction was 20, accounting for 77% of all stations. The WN + FN noise percentage decreased in 3 of 26 stations (YNSD, YNSM, YNTH). The number of the stations in WN+PL before and after HYDL correction was 3 and 4, respectively, while WN+GGM before and after HYDL correction was 1 and 2, indicating that HYDL correction changed the optimal noise model slightly for most stations.

The results of velocity and velocity uncertainty estimated from the vertical GPS time series using white noise had a certain error at most stations [42,43]. To obtain the velocity and velocity uncertainty precisely, the combination of different noise was used to estimate velocity and velocity uncertainty from the vertical GPS time series for all stations (except for YNGM). As shown in Figure 14, the range of velocity and velocity uncertainty difference was −0.25~0.25 mm/a at most stations, while the maximum value reached ~−0.5 mm/a, corresponding to the YNXP station and indicating that HYDL effects cannot be ignored when the velocity and velocity uncertainty are estimated from the vertical GPS time series.

4.2. The Effects Factors of Vertical GPS Time Series on the Seasonal Variations in Yunnan

The seasonal variations in vertical GPS time series are mainly influenced in two different ways [19]: a technical factor about the GPS calculation or observation (e.g., GPS draconitic signal effects or reference frame effects), and a geophysical factor about the redistribution of mass loading (e.g., HYDL, ATML, NTOL, bedrock thermal expansion, poroelastic deformation). Tan et al. [20] used independent component analysis to separate the draconitic signal at 22 GPS stations, and removed them from vertical GPS time series, causing the annual amplitude to be reduced by ~1 mm on average. Zhan et al. [19] investigated the reference frame influence on the annual amplitude, and the results showed a value limited to 1 mm with the 132 frame stations and the nonlinear reference frame. Yan et al. [44] studied the thermal effects on the vertical GPS time series in China, and the results showed that the thermal effects are limited in low-latitude. The annual amplitude of KMIN and XIAG stations in Yunnan was 0.4 mm. Thus, the thermal expansion had a slight effect on the vertical GPS time series in Yunnan. In this paper, we made a comparison between vertical GPS time series and 4 HYDL, 5 ATML, and 2 NTOL products. According to Section 3.1, the mean value of RMS reduction after HYDL, ATML, and NTOL correction in the best products was 1.32, −0.07, and 0.004 mm, respectively, implying that HYDL was the main factor causing the seasonal variations at most stations in Yunnan. The ATML and NTOL had a negative effect on the seasonal variations, but the range of the ATML time series was −5~5 mm. Considering that Yunnan is located in the southeastern part of the Tibet Plateau, the geological structure is very complex. In addition, the establishment of the ATML model was based on the global observation data. Measured data are less used in China. Therefore, the ATML model was often inaccurately estimated, because it does not usually reflect the real atmospheric pressure in Yunnan. Because it is inland, NTOL had almost no effect on the seasonal variations.

5. Conclusions

In this paper, we comprehensively analyzed the environmental loading effects on the vertical GPS time series of 27 stations in Yunnan using different HYDL, ATML, NTOL, and their combined products provided by GFZ, EOST, and IMLS, which related to the correction, RMS reduction, annual timescale, noise characteristic, velocity, and its uncertainty. Some conclusions are as follows:

The mean correlation between time series derived from HYDL (GFZ_LSDM, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim) and the vertical GPS time series was 0.51, 0.54, 0.55, and 0.49, respectively, while the range of RMS reduction after HYDL correction was −2.19~2.18, −0.64~2.38, −0.88~2.56, and −1.32~1.95 mm with an average value of 0.79, 1.25, 1.32, and 1.02 mm, respectively. According to the correlation and RMS reduction result, HYDL was the main factor causing seasonal variations at most GPS stations in Yunnan. Moreover, the IMLS_GEOSFPIT product represented the best choice for correcting seasonal variations at most stations. The mean RMS reduction after ATML (GFZ_ECMWF, IMLS_MERRA2, IMLS_GEOSFPIT, EOST_ERA_interim, EOST_ECMWF(IB)) correction was −0.13, −0.07, −0.1, −0.17, and −0.14 mm, respectively, while the range of RMS reduction after NTOL (GFZ_EMPIOM, IMLS_MPIOM06) correction was −0.07~0.04 and −0.09~0.05 mm with an average value of 0.004 and 0 mm, respectively, indicating that the ATML and NTOL models used in this paper cannot be effective to correct seasonal variations. Among the 40 different combinations of loading products correction, the range of mean correlation and RMS reduction value was 0.5~0.57, and 0.91~1.24 mm, respectively. The mean RMS reduction value of the optimal combined loading products, corresponding to HYDL (GFZ_LSDM), ATML (EOST_ERA_interim), and NTOL (IMLS_MPIOM06), was 1.24 mm, which was worse than that of the HYDL (IMLS_GEOSFPIT)-only correction, implying that the ATML and NTOL models used in this paper may have increased the error when combined with the HYDL model.

Common powers were located in the period of 256~512 days in the whole time span between the GPS and HYDL, ATML, NTOL, and their combined time series, and the range of mean XWT-based semblance between the GPS and HYDL time series was 0.9~1 at most GPS stations, while the range between the GPS and ATML and NTOL time series was very low. The results showed that HYDL explained the annual variations at most stations, while ATML and NTOL did so at YNGM station.

The optimal noise model at most stations was WN + FN before and after HYDL correction. The maximum value of velocity uncertainty reached −0.5 mm/a, implying that the HYDL effect cannot be ignored. Although we used four different HYDL products to investigate the effect on the vertical GPS time series, they could not detect the effect of groundwater. Therefore, we will use GRACE data to investigate their effects on GPS data in the future.

Author Contributions

All the authors participated in editing and reviewing the manuscript. S.H., K.C., H.Z., C.X., T.W., Z.Y. and Q.Z. conceived and designed the experiments; S.H. performed the experiments; S.H., K.C., H.Z., C.X., T.W., Z.Y. and Q.Z. analyzed and interpreted the results. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

The vertical GPS time series data (http://data.earthquake.cn, accessed on 22 May 2022) used in this paper was provided by “China Earthquake Networks Center, National Earthquake Data Center”. The environmental loading data were provided by the GFZ (http://rz-vm115.gfz-potsdam.de:8080/repository, accessed on 22 May 2022), EOST (http://loading.u-strasbg.fr/, accessed on 22 May 2022), and IMLS (http://massloading.net, accessed on 22 May 2022). The GMT and Origin software was used to plot the figures in this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The location distribution of 27 GPS stations in Yunnan represented by red dots. The red line represents the boundary of Yunnan. The black line represents the Xiaojiang Fault.

Figure 2. The RegEM interpolation result of YNMH, YNLA, YNXP, and YNCX stations. The red and green dots represent the vertical GPS time series and interpolation data, respectively.

Figure 3. The environmental loading time series provided by GFZ, IMLS, EOST, and vertical GPS time series at YNCX station: (a) four HYDL products, (b) five ATML products, (c) two NTOL products.

Figure 4. The correlation between four HYDL products and vertical GPS time series at 27 GPS stations.

Figure 5. The RMS reductions at 27 GPS stations after four different HYDL products correction.

Figure 6. The correlation between five ATML products and vertical GPS time series at 27 GPS stations.

Figure 7. The RMS reductions at 27 GPS stations after five ATML products correction.

Figure 8. The correlation between two NTOL products and the vertical GPS time series at 27 GPS stations.

Figure 9. The RMS reductions at 27 GPS stations after two NTOL products correction.

Figure 10. The cross wavelet power spectrum between GPS and HYDL (a), ATML (b), NTOL (c), and their combined (d) time series at YNCX station. The color bar represents the spectrum energy. The thin, black lines represent the COI delimiting the region not influenced by edge effects. Thick, black contours represent the significance test at the 95% level.

Figure 11. The mean XWT-based semblance between GPS and HYDL (a), ATML (b), NTOL (c), and their combined (d) time series for all stations at the period closet to 1 year.

Figure 12. The difference in annual phase between the GPS and HYDL (a), ATML (b), NTOL (c), and their combined (d) time series for all stations by the LSF and XWT methods.

Figure 13. Optimal noise model distribution before (a) and after (b) HYDL correction.

Figure 14. The difference in velocity and velocity uncertainty obtained by optimal noise model before and after HYDL correction.

Table 1

Environmental loading product provided by GFZ, IMLS, and EOST.

Organizations	Type	Model	Spatiotemporal Resolution	Time Span
GFZ	HYDL	LSDM	0.5° × 0.5°/24 h	1976–present
	ATML	ECMWF	0.5° × 0.5°/3 h
	NTOL	EMPIOM	1° × 1°/3 h
IMLS	HYDL	MERRA2	2′ × 2′/3 h	1980–present
	HYDL	GEOSFPIT	2′ × 2′/3 h	2000–present
	ATML	MERRA2	2′ × 2′/6 h	1980–present
	ATML	GEOSFPIT	2′ × 2′/3 h	2000–present
	NTOL	MPIOM06	2′ × 2′/3 h	1980–present
EOST	HYDL	ERA interim	0.5° × 0.5°/6 h	1979–2019
	ATML	ERA interim	0.5° × 0.5°/6 h
	ATML	ECMWF(IB)	0.5° × 0.5°/3 h	2000–present

Table 2

The mean correlation and RMS reduction between the vertical GPS time series with the combination of environmental loading products.

Combination			Correlation	RMS Reduction(mm)
HYDL	ATML	NTOL
GFZ_LSDM	EOST_ERA_interim	IMLS_MPIOM06	0.57	1.24
GFZ_LSDM	EOST_ECMWF(IB)	IMLS_MPIOM06	0.57	1.24
GFZ_LSDM	EOST_ERA_interim	GFZ_EMPIOM	0.57	1.24
GFZ_LSDM	EOST_ECMWF(IB)	GFZ_EMPIOM	0.57	1.23
GFZ_LSDM	GFZ_ECMWF	IMLS_MPIOM06	0.57	1.23
GFZ_LSDM	IMLS_GEOSFPIT	IMLS_MPIOM06	0.57	1.23
GFZ_LSDM	EOST_ERA_interim	GFZ_EMPIOM	0.57	1.23
GFZ_LSDM	EOST_ERA_interim	GFZ_EMPIOM	0.57	1.22
EOST_ERA_interim	EOST_ERA_interim	GFZ_EMPIOM	0.53	1.22
EOST_ERA_interim	IMLS_MERRA2	IMLS_MPIOM06	0.53	1.22
IMLS_MERRA2	IMLS_MERRA2	GFZ_EMPIOM	0.51	1.00
IMLS_MERRA2	IMLS_MERRA2	IMLS_MPIOM06	0.51	0.99
IMLS_MERRA2	IMLS_GEOSFPIT	GFZ_EMPIOM	0.51	0.99
IMLS_MERRA2	IMLS_GEOSFPIT	IMLS_MPIOM06	0.51	0.98
IMLS_MERRA2	EOST_ECMWF(IB)	GFZ_EMPIOM	0.51	0.93
IMLS_MERRA2	GFZ_ECMWF	IMLS_MPIOM06	0.51	0.93
IMLS_MERRA2	EOST_ECMWF(IB)	IMLS_MPIOM06	0.50	0.92
IMLS_MERRA2	GFZ_ECMWF	IMLS_MPIOM06	0.50	0.92
IMLS_MERRA2	EOST_ERA_interim	GFZ_EMPIOM	0.50	0.92
IMLS_MERRA2	EOST_ERA_interim	IMLS_MPIOM06	0.50	0.91

Table 3

The optimal noise model estimated from vertical GPS time series before and after HYDL correction.

Stations	Before Correction	After Correction	Stations	Before Correction	After Correction
KMIN	WN + PL	WN + PL	YNMZ	WN + FN	WN + FN
XIAG	WN + FN	WN + FN	YNRL	WN + FN	WN + FN
YNCX	WN + FN	WN + FN	YNSD	WN + FN	WN + GGM
YNDC	WN + FN	WN + FN	YNSM	WN + FN	WN + PL
YNHZ	WN + FN	WN + FN	YNTC	WN + FN	WN + FN
YNJD	WN + FN	WN + FN	YNTH	WN + FN	WN + PL
YNJP	WN + FN	WN + FN	YNWS	WN + FN	WN + FN
YNLA	WN + PL	WN + FN	YNXP	WN + FN	WN + FN
YNLC	WN + FN	WN + FN	YNYA	WN + FN	WN + FN
YNLJ	WN + FN	WN + FN	YNYL	WN + FN	WN + FN
YNMH	WN + FN	WN + FN	YNYM	WN + FN	WN + FN
YNMJ	WN + FN	WN + FN	YNYS	WN + GGM	WN + GGM
YNML	WN + FN	WN + FN	YNZD	WN + PL	WN + PL

References

1. Le Mével, H.; Feigl, K.L.; Córdova, L.; DeMets, C.; Lundgren, P.J. Evolution of unrest at Laguna del Maule volcanic field (Chile) from InSAR and GPS measurements, 2003 to 2014. Geophys. Res. Lett.; 2015; 42, pp. 6590-6598. [DOI: https://dx.doi.org/10.1002/2015GL064665]

2. Pan, Y.; Hammond, W.C.; Ding, H.; Mallick, R.; Jiang, W.; Xu, X.; Shum, C.; Shen, W.J.J. GPS Imaging of Vertical Bedrock Displacements: Quantification of Two-Dimensional Vertical Crustal Deformation in China. J. Geophys. Res. Solid Earth; 2021; 126, e2020JB020951. [DOI: https://dx.doi.org/10.1029/2020JB020951]

3. Wang, W.; Qiao, X.; Yang, S.; Wang, D.J. Present-day velocity field and block kinematics of Tibetan Plateau from GPS measurements. Geophys. J. Int.; 2017; 208, pp. 1088-1102. [DOI: https://dx.doi.org/10.1093/gji/ggw445]

4. Kreemer, C.; Blewitt, G.J.J. Robust estimation of spatially varying common-mode components in GPS time-series. J. Geod.; 2021; 95, pp. 1-19. [DOI: https://dx.doi.org/10.1007/s00190-020-01466-5]

5. Li, Z.; Chen, W.; van Dam, T.; Rebischung, P.; Altamimi, Z.J.J. Comparative analysis of different atmospheric surface pressure models and their impacts on daily ITRF2014 GNSS residual time series. J. Geod.; 2020; 94, pp. 1-20. [DOI: https://dx.doi.org/10.1007/s00190-020-01370-y]

6. Nicolas, J.; Verdun, J.; Boy, J.-P.; Bonhomme, L.; Asri, A.; Corbeau, A.; Berthier, A.; Durand, F.; Clarke, P.J. Improved Hydrological Loading Models in South America: Analysis of GPS Displacements Using M-SSA. Remote Sens.; 2021; 13, 1605. [DOI: https://dx.doi.org/10.3390/rs13091605]

7. Vandam, T.M.; Blewitt, G.; Heflin, M.B. Atmospheric pressure loading effects on Global Positioning System coordinate determinations. J. Geophys. Res. Solid Earth; 1994; 99, pp. 23939-23950. [DOI: https://dx.doi.org/10.1029/94JB02122]

8. Williams, S.; Penna, N.J. Non-tidal ocean loading effects on geodetic GPS heights. Geophys. Res. Lett.; 2011; 38, [DOI: https://dx.doi.org/10.1029/2011GL046940]

9. Wu, Y.; Zhao, Q.; Zhang, B.; Wu, W.J. Characterizing the seasonal crustal motion in Tianshan area using GPS, GRACE and surface loading models. Remote Sens.; 2017; 9, 1303. [DOI: https://dx.doi.org/10.3390/rs9121303]

10. Xiang, Y.; Yue, J.; Li, Z.J.A. Joint analysis of seasonal oscillations derived from GPS observations and hydrological loading for mainland China. Adv. Space Res.; 2018; 62, pp. 3148-3161. [DOI: https://dx.doi.org/10.1016/j.asr.2018.08.028]

11. Li, C.; Huang, S.; Chen, Q.; Dam, T.V.; Fok, H.S.; Zhao, Q.; Wu, W.; Wang, X.J. Quantitative evaluation of environmental loading induced displacement products for correcting GNSS time series in CMONOC. Remote Sens.; 2020; 12, 594. [DOI: https://dx.doi.org/10.3390/rs12040594]

12. Yuan, P.; Li, Z.; Jiang, W.; Ma, Y.; Chen, W.; Sneeuw, N.J. Influences of environmental loading corrections on the nonlinear variations and velocity uncertainties for the reprocessed global positioning system height time series of the crustal movement observation network of China. Remote Sens.; 2018; 10, 958. [DOI: https://dx.doi.org/10.3390/rs10060958]

13. Farrell, W.E. Deformation of the Earth by surface loads. Rev. Geophys.; 1972; 10, pp. 761-797. [DOI: https://dx.doi.org/10.1029/RG010i003p00761]

14. Wu, S.; Nie, G.; Meng, X.; Liu, J.; He, Y.; Xue, C.; Li, H. Comparative Analysis of the Effect of the Loading Series from GFZ and EOST on Long-Term GPS Height Time Series. Remote Sens.; 2020; 12, 2822. [DOI: https://dx.doi.org/10.3390/rs12172822]

15. Andrei, C.-O.; Lahtinen, S.; Nordman, M.; Näränen, J.; Koivula, H.; Poutanen, M.; Hyyppä, J.J. GPS time series analysis from aboa the finnish antarctic research station. Remote Sens.; 2018; 10, 1937. [DOI: https://dx.doi.org/10.3390/rs10121937]

16. Liang, S.; Gan, W.; Shen, C.; Xiao, G.; Liu, J.; Chen, W.; Ding, X.; Zhou, D.J. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. J. Geophys. Res. Solid Earth; 2013; 118, pp. 5722-5732. [DOI: https://dx.doi.org/10.1002/2013JB010503]

17. Hao, M.; Freymueller, J.T.; Wang, Q.; Cui, D.; Qin, S. Vertical crustal movement around the southeastern Tibetan Plateau constrained by GPS and GRACE data. Earth Planet. Sci. Lett.; 2016; 437, pp. 1-8. [DOI: https://dx.doi.org/10.1016/j.epsl.2015.12.038]

18. Sheng, C.-Z.; Gan, W.-J.; Liang, S.-M.; Chen, W.-T.; Xiao, G.-R. Identification and elimination of non-tectonic crustal deformation caused by land water from GPS time series in the western Yunnan province based on GRACE observations. Chin. J. Geophys.; 2014; 57, pp. 42-52.

19. Zhan, W.; Li, F.; Hao, W.; Yan, J. Regional characteristics and influencing factors of seasonal vertical crustal motions in Yunnan, China. Geophys. J. Int.; 2017; 210, pp. 1295-1304. [DOI: https://dx.doi.org/10.1093/gji/ggx246]

20. Tan, W.; Dong, D.; Chen, J.J. Application of independent component analysis to GPS position time series in Yunnan Province, southwest of China. Adv. Space Res.; 2022; 69, pp. 4111-4122. [DOI: https://dx.doi.org/10.1016/j.asr.2022.03.016]

21. Zhang, K.; Wang, Y.; Gan, W.; Liang, S.J. Impacts of Local Effects and Surface Loads on the Common Mode Error Filtering in Continuous GPS Measurements in the Northwest of Yunnan Province, China. Sensors; 2020; 20, 5408. [DOI: https://dx.doi.org/10.3390/s20185408] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32967242]

22. Liu, B.; Yu, W.; Dai, W.; Xing, X.; Kuang, C.J. Estimation of Terrestrial Water Storage Variations in Sichuan-Yunnan Region from GPS Observations Using Independent Component Analysis. Remote Sens.; 2022; 14, 282. [DOI: https://dx.doi.org/10.3390/rs14020282]

23. Zhan, W. Study on Vertical Crustal Motion in Chinese Mainland and Typical Areas Based on Continuous GPS. Wuhan Wuhan Univ.; 2017; 210, pp. 1295-1304.

24. Li, Y.J. Analysis of GAMIT/GLOBK in high-precision GNSS data processing for crustal deformation. Earthq. Res. Adv.; 2021; 1, 100028. [DOI: https://dx.doi.org/10.1016/j.eqrea.2021.100028]

25. Zhao, B.; Huang, Y.; Zhang, C.; Wang, W.; Tan, K.; Du, R. Crustal deformation on the Chinese mainland during 1998–2014 based on GPS data. Geod. Geodyn.; 2015; 6, pp. 7-15. [DOI: https://dx.doi.org/10.1016/j.geog.2014.12.006]

26. Schneider, T. Analysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values. J. Clim.; 2001; 14, pp. 853-871. [DOI: https://dx.doi.org/10.1175/1520-0442(2001)014<0853:AOICDE>2.0.CO;2]

27. Hu, S.; Wang, T.; Guan, Y.; Yang, Z.J. Analyzing the seasonal fluctuation and vertical deformation in Yunnan province based on GPS measurement and hydrological loading model. Chin. J. Geophys. Chin. Ed.; 2021; 64, pp. 2613-2630.

28. Li, W.; Li, F.; Zhang, S.; Lei, J.; Zhang, Q.; Yuan, L.J. Spatiotemporal filtering and noise analysis for regional GNSS network in Antarctica using independent component analysis. Remote Sens.; 2019; 11, 386. [DOI: https://dx.doi.org/10.3390/rs11040386]

29. Liu, R.; Li, J.; Fok, H.S.; Shum, C.; Li, Z.J.S. Earth surface deformation in the north China plain detected by joint analysis of GRACE and GPS data. Sensors; 2014; 14, pp. 19861-19876. [DOI: https://dx.doi.org/10.3390/s141019861]

30. Dill, R.; Dobslaw, H.J. Numerical simulations of global-scale high-resolution hydrological crustal deformations. J. Geophys. Res. Solid Earth; 2013; 118, pp. 5008-5017. [DOI: https://dx.doi.org/10.1002/jgrb.50353]

31. Dill, R. Hydrological Model LSDM for Operational Earth Rotation and Gravity Field Variations; Scientific Technical Report STR08/09; GFZ German Research Centre For Geosciences: Potsdam, Germany, 2008.

32. Marsland, S.J.; Haak, H.; Jungclaus, J.H.; Latif, M.; Röske, F.J. The Max-Planck-Institute global ocean/sea ice model with orthogonal curvilinear coordinates. Ocean Model.; 2003; 5, pp. 91-127. [DOI: https://dx.doi.org/10.1016/S1463-5003(02)00015-X]

33. Berrisford, P.; Dee, D.; Fielding, K.; Fuentes, M.; Kallberg, P.; Kobayashi, S.; Uppala, S.J. The ERA-interim archive. ERA Rep. Ser.; 2009; 1, pp. 1-16.

34. Reichle, R.H.; Draper, C.S.; Liu, Q.; Girotto, M.; Mahanama, S.P.; Koster, R.D.; De Lannoy, G.J. Assessment of MERRA-2 land surface hydrology estimates. J. Clim.; 2017; 30, pp. 2937-2960. [DOI: https://dx.doi.org/10.1175/JCLI-D-16-0720.1]

35. Gu, Y.; Yuan, L.; Fan, D.; You, W.; Su, Y.J. Seasonal crustal vertical deformation induced by environmental mass loading in mainland China derived from GPS, GRACE and surface loading models. Adv. Space Res.; 2017; 59, pp. 88-102. [DOI: https://dx.doi.org/10.1016/j.asr.2016.09.008]

36. Grinsted, A.; Moore, J.C.; Jevrejeva, S.J. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes Geophys.; 2004; 11, pp. 561-566. [DOI: https://dx.doi.org/10.5194/npg-11-561-2004]

37. Cooper, G.; Cowan, D.J.C. Geosciences. Comparing time series using wavelet-based semblance analysis. Comput. Geosci.; 2008; 34, pp. 95-102. [DOI: https://dx.doi.org/10.1016/j.cageo.2007.03.009]

38. Peng, Y.; Scales, W.A.; Hartinger, M.D.; Xu, Z.; Coyle, S.J.S.N. Characterization of multi-scale ionospheric irregularities using ground-based and space-based GNSS observations. Satellite Navigation; 2021; 2, pp. 1-21. [DOI: https://dx.doi.org/10.1186/s43020-021-00047-x]

39. Yang, F.; Meng, X.; Guo, J.; Yuan, D.; Chen, M.J.S.N. Development and evaluation of the refined zenith tropospheric delay (ZTD) models. Satellite Navigation; 2021; 2, pp. 1-9. [DOI: https://dx.doi.org/10.1186/s43020-021-00052-0]

40. Bos, M.; Fernandes, R.; Williams, S.; Bastos, L.J. Fast error analysis of continuous GNSS observations with missing data. J. Geod.; 2013; 87, pp. 351-360. [DOI: https://dx.doi.org/10.1007/s00190-012-0605-0]

41. He, X.; Bos, M.S.; Montillet, J.-P.; Fernandes, R.; Melbourne, T.; Jiang, W.; Li, W.J. Spatial Variations of Stochastic Noise Properties in GPS Time Series. Remote Sens.; 2021; 13, 4534. [DOI: https://dx.doi.org/10.3390/rs13224534]

42. He, X.; Bos, M.; Montillet, J.; Fernandes, R.J. Investigation of the noise properties at low frequencies in long GNSS time series. J. Geod.; 2019; 93, pp. 1271-1282. [DOI: https://dx.doi.org/10.1007/s00190-019-01244-y]

43. Williams, S.D.; Bock, Y.; Fang, P.; Jamason, P.; Nikolaidis, R.M.; Prawirodirdjo, L.; Miller, M.; Johnson, D.J. Error analysis of continuous GPS position time series. J. Geophys. Res. Solid Earth; 2004; 109, [DOI: https://dx.doi.org/10.1029/2003JB002741]

44. Yan, H.M.; Chen, W.; ZHU, Y.Z.; ZHANG, W.M.; Zhong, M.; Liu, G.Y. Thermal effects on vertical displacement of GPS stations in China. Chin. J. Geophys.; 2010; 53, pp. 252-260. [DOI: https://dx.doi.org/10.1002/cjg2.1492]