Potential Contributors to CME and Optimal Noise Model Analysis in the Chinese Region Based on Different HYDL Models

Abstract

Optimizing the noise model for global navigation satellite system (GNSS) vertical time series is vital to obtain reliable uplift (or subsidence) deformation velocity fields and assess the associated uncertainties. In this study, by thoroughly considering the effects of hydrological loading (HYDL) that dominates the seasonal fluctuations and common mode error (CME), we analyzed the optimal noise characteristics of GNSS vertical time series at 39 stations spanning from January 2011 to August 2019 in the Chuandian region, southeast of the Qinghai–Tibet Plateau. Our results showed that the optimal noise models without HYDL correction were white noise plus flicker noise (WN + FN), white noise plus power law noise (WN + PL), and white noise plus Gauss–Markov noise (WN + GGM), which accounted for 87%, 10%, and 3% of GNSS stations, respectively. By contrast, the optimal noise models at all stations were WN + FN and WN + PL after correction by different HYDLs. The correlation between CME and HYDL provided by the School and Observatory of Earth Sciences (EOST), namely EOST_HYDL, was 0.63~0.8 and the value of RMS reduction was 18.9~40.3% after removing EOST_HYDL time series from the CME, with a mean value of 31.8%, there is a good correlation and consistency between CME and EOST_HYDL. The absolute value of vertical velocity and its uncertainty with and without EOST_HYDL correction varied from 0.11 to 0.55 mm/a and 0 to 0.23 mm/a, respectively, implying that the effect of HYDL should not be neglected when performing optimal noise model analysis for GNSS vertical time series in the Chuandian region.

1. Introduction

Over the past decades, global navigation satellite system (GNSS) time series have been widely used for crust deformation monitoring [1,2,3,4], regional reference frame establishment and maintenance [5,6], seismic source imaging [7,8,9,10,11], hydrological drought monitoring [12,13,14,15], etc. It is well known that GNSS vertical time series is not only driven by secular tectonic movement, including basin subsidence and mountain uplift, but also has seasonal variations caused by environmental loadings, such as hydrological loading (HYDL), atmospheric loading (ATML), nontidal ocean loading (NTOL), and systematic errors, including inaccurate orbit and troposphere bias [16,17]. To obtain precise vertical GNSS time series, much research effort has focused on spatial filtering [18,19,20] and environmental loading models [21,22,23]. Additionally, an optimal noise model of the GNSS time series is also important for GNSS velocity and its uncertainty estimation, which will provide reliable data to support reasonable geodynamic and geodesy interpretation [24].

It has been demonstrated that the noise associated with GNSS time series can be characterized as white noise (WN) and colored noise including power-law noise (PL), band-pass power-law noise (BPPL), flicker noise (FN), random walk noise (RW), Gauss–Markov model noise (GGM), and first-order Gauss–Markov noise (FOGM) [25,26]. The GNSS velocity and its uncertainty can be inaccurate when using an improper noise model estimated from the GNSS time series. For simplicity, without considering the background environmental effects the noise behaves as WN + FN [24,27,28]. The noise characteristic is more complicated in cases where the diverse environmental loading effects are taken into account. For instance, Li et al. [29] compared the noise characteristics of 38 GNSS stations in the Eurasia region before and after ATML plus HYDL correction; the results showed that the WN percentage increased and the FN percentage decreased. He et al. [30] used various noise model combinations to analyze the noise characteristics of 671 global GNSS stations from 2000 to 2021 before and after HYDL correction; the results showed that the optimal noise models for all stations in three components were diverse, which can be mainly represented by WN + FN, GGM, and WN + PL. Li et al. [31] analyzed the noise characteristics of 11 IGS stations in China’s mainland spanning from 1995 to 2010 before and after environmental loading (HYDL, ATML, and NTOL) effects; the results showed that the optimal noise models were composed of WN + FN, WN + BPPL, WN + FN + RW, WN + PL, and WN + RW + FOGM before the environmental loading correction, which accounted for 55%, 24%, 9%, 9%, and 3%, respectively. After the environmental loading correction, the optimal noise models changed mainly to those represented as WN + BPPL and WN + RW + FOGM.

Clearly, to obtain a reliable and accurate velocity together with its uncertainty, the optimal noise model of GNSS time series considering environmental loading effects in different regions should be established. Located in the southeastern part of the Qinghai–Tibet Plateau, the Chuandian region has many active fault zones and complex geological structures [32], including the Yushu–Ganzi–Xianshuihe fault zone, Xiaojiang fault zone, and Honghe fault zone. Investigating the optimal noise model of GNSS time series affected by environmental loading in the Chuandian region will improve the accuracy of the GNSS velocity and its uncertainty estimation, which will help to reveal the tectonic evolution and crustal deformation of the Qinghai–Tibet Plateau. In fact, numerous studies have focused on the factors influencing GNSS time series in the Chuandian region. For instance, Zhang et al. [33] researched the influence of environmental loading on CME filtering of GNSS time series in the Chuandian region; the results showed that it is necessary to consider suitable environmental loading models to correct the GNSS time series when CME filtering is not conducted. Tan et al. [34] analyzed the potential environmental loading interpretation of CME in the Chuandian region and showed that soil moisture mass loading can interpretation part of CME. Hu et al. [35] found that the uncertainty in GNSS velocity could reach −0.5 mm/a before and after HYDL correction, indicating that the HYDL effect must be considered in the Chuandian region when velocity and its uncertainty estimated from the GNSS time series based on optimal noise model. Furthermore, numerous studies have demonstrated that HYDL is the main factor that induces seasonal variations in GNSS vertical time series in the Chuandian region [36,37,38]. However, significant variance exists among the HYDL from different institutions due to different computational strategies and loading models [39]. For instance, Wu et al. [40] used the HYDL time series provided by various institutions to research its effects on GNSS vertical time series for 633 global stations. The RMS reduction rates were positive at 82.6 and 87.4%. Andrei et al. [41] found that the difference in GNSS vertical time series at a GNSS station from various institutions could reach 3 mm. Previous studies have analyzed the optimal noise models of GNSS time series considering only HYDL effects [30,42,43]; however, the optimal noise models of GNSS time series are diverse due to different HYDL correction, which results in an inaccurate estimation of velocity and its uncertainty. Therefore, it is necessary to investigate the optimal noise model of GNSS time series considering different HYDL effects in the Chuandian region.

In this study, we analyzed the noise characteristics of GNSS vertical time series in the Chuandian region considering different hydrological loading from various institutions. In addition, we also investigated the relationship between CME and HYDL time series and the effects of HYDL on GNSS time series. The structure of this paper is as follows: Section 2 describes the GNSS vertical time series and HYDL data, noise model, and principal component analysis methods. Section 3 reports the results of CME extraction and optimal noise analysis before and after three HYDL corrections. Section 4 discusses the potential factor of CME about HYDL in the Chuandian region and the changes in velocity and its uncertainty affected by HYDL. Section 5 presents the conclusions of this paper.

2. Data and Methods

2.1. GNSS Vertical Time Series

We selected the daily GNSS vertical time series of 39 stations in the Chuandian region, which is derived from the crustal movement observation network of China (CMONOC) spanning from January 2011 to August 2019, the GNSS vertical time series of all stations in the Chuandian region was provided by ftp://ftp.cgps.ac.cn/products/position/ (accessed on 13 November 2022). Figure 1 displays the detailed spatial distribution of all stations in the Chuandian region. The daily GNSS vertical time series was calculated using the GAMIT/GLOBK software [44]. First, the unconstrained daily GNSS solutions for all stations were calculated using GAMIT software, which contains station positions, baseline, and satellite orbits. Main strategies for GAMIT solution used are as follows: (1) The elevation cutoff angle was set at 10° for all stations. (2) The effects of ocean tides and solid Earth tides were corrected using the FES2004 [45] and IERS10 [46] models, respectively. (3) Antenna offsets of receivers were corrected using absolute antenna phase calibration models. Then, unconstrained daily GNSS solutions were transformed into the International Reference Frame (ITRF) 2014 by seven parameter transformations of translation, rotation, and scale through GLOBK software. More detailed information about the processing of GNSS data from CMONOC can be obtained at ftp://ftp.cgps.ac.cn/doc/ (accessed on 13 November 2022).

There were some outliers, position offsets, and missing data in the GNSS vertical time series. Therefore, we had to preprocess the GNSS vertical time series before noise characteristic analysis and CME extraction. The main strategies for preprocessing used in this study are as follows: (1) Outliers were detected and removed using the interquartile range (IQR) rule. (2) Position offsets in GNSS vertical time series, which are caused by earthquakes or instrument changes, were corrected by the least squares fitting (LSF) method. (3) The regularized expectation maximization method (RegEM) [47], which consider objects and correlation among all GNSS stations, was used to interpolate the missing data. This method is used extensively for gap-filling in missing GNSS time series [35,48,49]. Figure 2 shows the position offset correction results for the YNRL station and the interpolation results for the SCDF and SCNN stations. It can be seen from Figure 2c,d that the interpolation data showed good consistency with the measured GNSS time series, which displayed obvious seasonal variations.

2.2. Hydrological Loading Time Series

The HYDL time series used in this paper was obtained from the School and Observatory of Earth Sciences (EOST, http://loading.u-strasbg.fr/ (accessed on 13 November 2022)) [50], German Research Center for Geosciences (GFZ, http://rz-vm115.gfz-potsdam.de:8080/repository (accessed on 13 November 2022)) [51], and the International Mass Loading Service (IMLS, http://massloading.net (accessed on 13 November 2022)) [52]. EOST models the global grids using spherical harmonics formalism, the HYDL time series provided by EOST, namely EOST_HYDL, was calculated using the ERA interim model [53], which was estimated from the datasets of ECMWF reanalysis, with a temporal–spatial resolution of ${0.5}^{\circ }\times {0.5}^{\circ }\times 6\mathrm{h}$; this hydrological model contains the water in soil moisture and snow. GFZ used the Green’s function approach to model global grids; the HYDL time series provided by GFZ, namely GFZ_HYDL, was calculated using the land surface discharge model (LSDM) [54], with a temporal–spatial resolution of ${0.5}^{\circ }\times {0.5}^{\circ }\times 24\mathrm{h}$; this hydrological model contains the water in soil moisture, ice, snow, and runoff and drainage. IMLS used spherical harmonics and load love numbers (LLN) to model the global grids; the HYDL time series provided by IMLS, namely IMLS_HYDL, was calculated using the global Earth-observing system forward processing instrumental team (GEOSFPIT) [55], with a temporal–spatial resolution of ${0.5}^{’}\times {0.5}^{’}\times 3\mathrm{h}$; this hydrological model included soil moisture, water content in the canopy, and water equivalent of snow cover. We used a bicubic interpolation method to calculate GFZ_HYDL time series at each station based on latitude and longitude in the Chuandian region. Based on gridded files, we calculated the EOST_HYDL time series using the interpolation method at each station. IMLS provides precomputed series of loading time series for on-demand computation service; we provided the XYZ coordinates at each station to the server for calculating the IMLS_HYDL time series in the Chuandian region. To unite the daily resolution with the GNSS vertical time series, the HYDL time series obtained from EOST and IMLS needed to be averaged into daily time series. In addition, the EOST_HYDL, GFZ_HYDL, and IMLS_HYDL time series were linked to the center-of-figure (CF) reference frame for consistency with the GNSS vertical time series.

2.3. Noise Model

GNSS time series can be composed of two parts: deterministic and stochastic models [56]. The former mainly consist of a linear trend, seasonal items (including annual and semi-annual), and position offsets. The stochastic model includes the noise, which was used to analyze the noise characteristics of the GNSS time series in this study. The function model of the GNSS time series can be expressed using the following equation:

(1)$\begin{array}{l}y\left(t_i\right)=a+b{t}_i+c\sin (2\pi t_{\mathrm{i}})+d\cos (2\pi t_i)+e\sin (4\pi t_i)+f\cos (4\pi t_i)\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{j=1}^{n_g}{g}_{j}H(t_i}-T_{gj})+v_i\end{array}$

where $t_i$ represents the unit of the year; the coefficients $a$ and $b$ represent the initial GNSS position and linear movement rate, respectively; $c,d,e$, and $f$ describe the annual and semiannual amplitudes in the periodic term; $g_j$ represents the position offsets caused by earthquakes or instrument changes; $H(t_i-T_{gj})$ is the Heaviside step function; and $v_i$ represents the GNSS residual time series, which is composed of different noise models.

Hector software [57], which can calculate the linear trend, seasonal items (mainly including annual and semiannual), and various noise model combination in GNSS time series, was used to investigate the noise characteristics of all stations in the Chuandian region. The Akaike information criterion (AIC) and Bayesian information criterion (BIC) information standards were used to choose the optimal noise model [58,59]. Furthermore, the maximum likelihood was used as a starting point for estimation, to avoid overfitting; it also added penalties for the parameters. The optimal noise model usually has a smaller BIC or AIC value when using the Hector software to analyze the noise characteristic of GNSS vertical time series [30]. Compared to the AIC, the BIC accounts for more parameters. Therefore, we used the BIC evaluation standards to choose the optimal noise model. Hector software can use many types of noise models, mainly including WN, FN, GGM, PL, RW, and their combinations. More detailed information concerning the noise models from the Hector software can be found at http://segal.ubi.pt/hector/manual_1.7.2.pdf (accessed on 13 November 2022).

2.4. Principal Component Analysis

Principal component analysis (PCA) is an orthogonal decomposition algorithm that decomposes a group time series into a set of linear spatial and temporal orthogonal eigenvectors and corresponding PCs. PCA can reduce the dimensionality of a group time series through some PCs, which can represent the most informative internal structure in the original time series [60].

For a regional GNSS network of n stations, the daily residual time series with a span of m days, in which the linear and seasonal items were removed from the original GNSS time series using the LSF method, can be defined as $X\left(t_i,x_j\right)\left(i=1,2,\dots ,m;j=1,2,\dots ,n\right)$. The covariance matrix B can be expressed using the following equation:

(2)$b_{i,j}=\frac{1}{m-1}{\displaystyle \sum _{k=1}^{m}X\left(t_k,x_i\right)}X\left(t_k,x_j\right)$

The symmetric matrix B can be decomposed as:

(3)$B=V{\Lambda }{V}^T$

where the eigenvector matrix $V^T$ is an $n\times n$ matrix and ${\Lambda }$ is a matrix with $k$ nonzero diagonal eigenvalues. For most of the GNSS time series data, the symmetric matrix B is a full-rank matrix. From linear algebra, the daily residual time series $X\left(t_i,x_j\right)$ can be defined using the orthogonal function basis $v$ as follows:

(4)$X\left(t_i,x_j\right)={\displaystyle \sum _{k=1}^n{a}_k}\left(t_i\right)v_k\left(x_j\right)$

where $a_k\left(t_i\right)$ is the kth PC of the matrix $X$, which represents a temporal signature of the mode, and can also be expressed using the following equation:

(5)$a_k\left(t_i\right)={\displaystyle \sum _{j=1}^{n}X\left(t_i,x_j\right)}{v}_k\left(x_j\right)$

where $v_k\left(x_j\right)$ is the $a_k\left(t_i\right)$ corresponding eigenvector, which represents the spatial response of the mode. We determined the $p$ PCs according to the descending eigenvalues, and the CME was expressed as:

(6)$CME\left(t_i,s_j\right)={\displaystyle \sum _{k=1}^p{a}_k\left(t_i\right)}{v}_k\left(x_j\right)$

3. Results

3.1. Comparison of Different HYDL Time Series with GNSS

To better investigate the relationship between HYDL time series and GNSS time series, we used the LSF method to remove the linear trend from the GNSS vertical time series. In addition, the time series induced by atmospheric loading (ATML) and nontidal oceanic loading (NTOL) was also removed from the GNSS vertical time series. GNSS vertical time series and corresponding HYDL time series provided by EOST, GFZ, and IMLS at some GNSS stations (SCML, SCXC, YNLJ, and YNYS) are shown in Figure 3. It can be seen from the figure that the HYDL time series and GNSS vertical time series showed an obvious seasonal variation. Moreover, the range of GNSS vertical time series at all GNSS stations in the Chuandian region was −20~20 mm, while the range of HYDL time series was −10~10 mm, implying that HYDL cannot fully explain the seasonal variations in GNSS vertical time series.

The root mean squares (RMS) reduction and correlation [48] were also used to quantitatively investigate the relationship between the GNSS vertical time series and the corresponding HYDL time series from the Chuandian region. The correlations between the GNSS vertical time series and HYDL time series (EOST_HYDL, GFZ_HYDL, and IMLS_HYDL) are shown in Figure 4. The correlations between the HYDL time series (EOST_HYDL, GFZ_HYDL, and IMLS_HYDL) and the GNSS vertical time series was 0.27~0.74, 0.27~0.75, and 0.25~0.73 with a mean value of 0.6, 0.6, and 0.54, respectively, indicating that the HYDL time series (EOST_HYDL, GFZ_HYDL, and IMLS_HYDL) showed good correlation with the GNSS vertical time series for most GNSS stations in the Chuandian region. Some stations have low correlation, which may have been affected by unknown local effects, GNSS system error, and noise [42,61].

From Figure 5, the value of RMS reduction between the HYDL time series (EOST_HYDL, GFZ_HYDL, and IMLS_HYDL) and GNSS vertical time series was 2~31%, −10~32%, and 1.5~27% with a mean value of 18.9%, 16.7%, and 16%, respectively. The RMS reduction between EOST_HYDL, IMLS_HYDL, and GNSS vertical time series was positive for all GNSS stations in the Chuandian region. According to correlation and RMS reduction results, the seasonal signals in GNSS vertical time series at most GNSS stations in the Chuandian region were affected by HYDL. Furthermore, EOST_HYDL had the most positive effect on seasonal variations in GNSS vertical time series when compared with GFZ_HYDL and IMLS_HYDL.

3.2. CME Extraction

Some spatial correlation noise, called CME, exists in regional GNSS networks, which can be derived from environmental loading effects, GNSS symmetric error, and GNSS observing environmental effects [62]. It has been proven that PCA is a reliable and useful method for extracting CME from regional GNSS networks [63,64]. Thus, the PCA method was used to extract the CME from 39 GNSS stations in the Chuandian region. Figure 6 displays the normalized spatial eigenvector of the first three principal components and the contribution rate for all principal component eigenvalues. Figure 6a indicates that the spatial response of PC1 showed little fluctuation, and the values of the spatial eigenvectors in PC1 were more than 0.25 for all GNSS stations. The spatial response of PC2 and PC3 fluctuated sharply, and thus could not satisfactorily reflect the common variation in the Chuandian region. In addition, the values of the spatial eigenvectors in PC2 and PC3 were less than 0.25 for most GNSS stations. From Figure 6b, the contribution rate of the first three PCs was 47%, 5%, and 4%, respectively.

The criteria for extracting the CME from the regional GNSS network can be found in [60]; the GNSS stations (>50%) had obvious normalized responses (>0.25) and the contribution rate of PC eigenvalues was greater than 1% of the sum of all contribution rates. Thus, the first PC model can be defined as the CME in the Chuandian region. The value of CME at most GNSS stations in the Chuandian region was −10~10 mm, which showed significant seasonal variations. The changes in GNSS vertical time series after filtering by PCA at some GNSS stations (SCML, SCXC, YNLJ, and YNYS) are shown in Figure 7. After filtering by PCA, the GNSS vertical time series were smoother and more convergent, implying that PCA methods can effectively improve the signal-to-noise ratio of GNSS vertical time series.

3.3. Optimal Noise Model Analysis before and after Different HYDL Correction

The five noise models: white noise WN; white noise plus random walk noise (WN + RW); white noise plus flicker noise plus random walk noise (WN + FN + RW); white noise plus flicker noise (WN + FN); and white noise plus power law noise (WN + PN), were used to analyze the noise characteristics of the 39 GNSS stations in the Chuandian region. The optimal noise models for all stations are shown in Figure 8. The optimal noise models were diverse but mainly WN + FN, which accounted for 87% of all GNSS stations. The optimal noise models of WN + PL and WN + GGM accounted for 10% and 3%, respectively.

Numerous studies have proved that velocity uncertainty is underestimated when only using the WN model to analyze the noise characteristics of GNSS vertical time series [65]. Thus, in this research, we estimated the velocity uncertainty from GNSS vertical time series based on the optimal noise model and WN model. The results of velocity uncertainty are shown in Figure 9. The value of velocity uncertainty estimated from the WN model was 0.02~0.06 mm/a with a mean of 0.04 mm/a, while the optimal noise model provided a value of 0.15~1.12 mm/a with a mean of 0.47 mm/a, indicating that the velocity uncertainty associated with the WN model was seriously underestimated in comparison with the optimal noise model at all stations in the Chuandian region.

Environmental loading, mainly including HYDL, ATML, and NTOL, can affect the stochastic models of GNSS vertical time series [66]. Among these loadings, HYDL is the main reason for the induction of seasonal variations in GNSS vertical time series in the Chuandian region. Referring to Section 3.1, the seasonal variations in GNSS vertical time series affected by HYDL, which were derived from various institutions, were diverse. Therefore, we reanalyzed the noise characteristics of the 39 GNSS stations after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction.

Figure 10 and Table 1 display the results of the optimal noise models after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction. The results showed that the optimal noise models for all stations were WN + FN and WN + PL after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction. Comparing with Figure 8, it can be seen that the optimal noise models of WN + FN obviously decreased, in which the percentage decreased by 23%, 13%, and 25%, respectively. The WN + PL also obviously increased, in which the percentage increased by 26%, 16%, and 28%, respectively. The optimal noise models before and after HYDL correction were mainly WN + FN and WN + PL. Our results are similar to those of other studies [30,35]. Hu et al. [35] investigated the optimal noise models before and after HYDL correction at 27 GNSS stations in the Yunnan region; the results showed the optimal noise models are diverse, mainly WN + FN and WN + PL. He et al. [30] analyzed the optimal noise models of global GNSS stations considering HYDL correction; the results showed the optimal noise models were WN + FN and WN + PL at most stations.

4. Discussion

4.1. Potential Hydrological Interpretation of the CME

Referring to Section 3.1, the seasonal variations in GNSS vertical time series caused by EOST_HYDL were more obvious in comparison with GFZ_HYDL and IMLS_HYDL. Therefore, the EOST_HYDL time series was used to further investigate the potential factor of CME. Figure 11 shows the CME and corresponding EOST_HYDL time series at some GNSS stations (SCML, SCXC, YNLJ, and YNYS), which displayed significant seasonal variations with values of −10~10 mm. There was good agreement between the CME and EOST_HYDL time series. In addition, we also used RMS reduction and correlation to quantitatively investigate the relationship between CME and the corresponding EOST_HYDL time series.

The results of correlation and RMS reduction between the CME and EOST_HYDL time series are shown in Figure 12 and Table 2. We observed good correlation and consistency between CME and EOST_HYDL. From Figure 12a, the value of correlation between CME and EOST_HYDL was 0.63~0.8 with a mean value of 0.74. From Figure 12b, the value of RMS reduction was 18.9~40.3% with a mean value of 31.8% at all GNSS stations after removal of the EOST_HYDL time series from the CME in the Chuandian region, implying that HYDL is the main factor causing the CME in the Chuandian region.

4.2. Changes in Velocity and Its Uncertainty due to Different HYDL Effects

To investigate the effects of HYDL on the estimation of velocity and its uncertainty from GNSS vertical time series in the Chuandian region, we estimated the velocity and its uncertainty from the GNSS vertical time series based on optimal noise models before and after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction. The results for the estimation of velocity and its uncertainty difference are shown in Figure 13. The absolute value of velocity difference was 0.11~0.55, 0.02~0.35, and 0.01~0.29 mm/a with an absolute mean value of 0.15, 0.16, and 0.13 mm/a, respectively. The absolute value of velocity uncertainty difference was 0~0.23, 0~0.28, and 0~0.26 mm/a, respectively, implying that the effect of HYDL cannot be ignored when using optimal noise models to estimate the velocity and its uncertainty from GNSS vertical time series in the Chuandian region.

5. Conclusions

This study was carried out based on the GNSS vertical time series of 39 stations spanning from January 2011 to August 2019 and the corresponding HYDL time series provided by EOST, GFZ, and IMLS in the Chuandian region. Firstly, we compared the GNSS vertical time series with the EOST_HYDL, GFZ_HYDL, and IMLS_EOST time series. Secondly, we used the PCA method to extract the CME from the GNSS vertical time series and discussed the potential source of CME. Finally, we analyzed the noise characteristics of the GNSS vertical time series from all stations before and after EOST_HYDL, GFZ_HYDL, and IMLS_EOST correction and discussed their effects on estimation of the velocity and its uncertainty. Some detailed conclusions are as follows:

• Both HYDL time series and GNSS vertical time series show an obvious seasonal variation. The value of correlation between the GNSS vertical time series and EOST_HYDL, GFZ_HYDL, and IMLS_HYDL was 0.27~0.74, 0.27~0.75, and 0.25~0.73 with a mean value of 0.6, 0.6, and 0.54; the value of RMS reduction was 2~31%, −10~32%, and 1.5~27% with a mean value of 18.9%, 16.7%, and 16% after removing EOST_HYDL, GFZ_HYDL, and IMLS_HYDL time series from the GNSS, respectively. Compared with GFZ_HYDL and IMLS_HYDL, EOST_HYDL had the most positive effect on seasonal variations in GNSS vertical time series in the Chuandian region.

• The value of CME at most stations was −10~10mm, which displayed significant seasonal variations. There was good agreement between the CME and EOST_HYDL time series. The value of correlation between CME and EOST_HYDL was 0.63~0.8, the value of RMS reduction was 18.9~40.3% with a mean value of 31.8% after removing EOST_HYDL time series from the CME, indicating that the HYDL effect is one of the principal factors causing the CME in the Chuandian region.

• The optimal noise models before HYDL correction were WN + FN, WN + PL, and WN + GGM, while the optimal noise models were WN + FN and WN + PL after corrected by HYDL. The absolute value of velocity difference was 0.11~0.55, 0.02~0.35, and 0.01~0.29 mm/a before and after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction, the absolute value of velocity uncertainty difference was 0~0.23, 0~0.28, and 0~0.26 mm/a, respectively. Thus, the influence of HYDL on the estimation of velocity and its uncertainty from GNSS vertical time series in the Chuandian region cannot be ignored.

• In this study, we used five common models to analyze the noise characteristics of 39 GNSS stations in the Chuandian region. However, the stochastic noise model of characteristics in the GNSS time series is very complex. Therefore, in the future we will use additional noise models to investigate the noise characteristics of GNSS vertical time series. Although the HYDL time series provided by EOST, GFZ, and IMLS were used to research the noise characteristics of the GNSS vertical time series, the HYDL time series used in this paper could not detect the effect of groundwater. Therefore, in the future we will use the GRACE model to investigate the effect of groundwater on GNSS vertical time series.

Author Contributions

All authors participated in editing and reviewing the manuscript. S.H., K.C., H.Z., T.W., Z.Y. and Q.Z. conceived and designed the experiments; S.H. performed the experiments; S.H., K.C., H.Z., 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 GNSS vertical time series data (http://data.earthquake.cn (accessed on 13 November 2022)) used in this research were provided by the 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 13 November 2022)), EOST (http://loading.u-strasbg.fr/ (accessed on 13 November 2022)), and IMLS (http://massloading.net (accessed on 13 November 2022)). GMT software was used to plot the figures in this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviation

ATML	atmospheric loading
AIC	Akaike information criterion
BIC	Bayesian information criterion
BPPL	band-pass power-law noise
CME	common mode error
CMONOC	crustal movement observation network of China
EOST	School and Observatory of Earth Sciences
FN	flicker noise
FOGM	first-order Gauss–Markov noise
GGM	Gauss–Markov noise
GFZ	German Research Center for Geosciences
GNSS	global navigation satellite system
GEOSFPIT	global Earth observing system forward processing instrumental team
HYDL	hydrological loading
IMLS	International Mass Loading Service
IQR	interquartile range
LSF	least squares fitting
LSDM	land surface discharge model
LLN	load Love numbers
NTOL	nontidal ocean loading
PL	power law noise
PCA	principal component analysis
RW	random walk noise
RMS	root mean squares
RegEM	regularized expectation maximization methods
WN	white noise

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Spatial distribution and data availability of 39 GNSS stations in the Chuandian region.

Figure 2. GNSS vertical time series of the YNRL station before (a) and after (b) position offset correction, and RegEM interpolation results for the SCDF (c) and SCNN (d) stations. The black line represents the original GNSS vertical time series of the YNRL station, the red line represents the GNSS vertical time series modeled by the LSF method, and the black and red dots represent the original GNSS vertical time series and interpolation data, respectively.

Figure 3. Comparison of the HYDL time series provided by EOST, GFZ, and IMLS with GNSS vertical time series at SCML (a), SCXC (b), YNLJ (c), and YNYS (d) stations.

Figure 4. Correlation between GNSS vertical time series and HYDL time series provided by EOST (a), GFZ (b), and IMLS (c) for all stations in the Chuandian region.

Figure 5. RMS reduction between GNSS vertical time series and HYDL time series provided by EOST (a), GFZ (b), and IMLS (c) for all stations in the Chuandian region.

Figure 6. (a) The normalized spatial eigenvectors of the first three principal components and (b) contribution rate of all principal component eigenvalues.

Figure 7. Unfiltered and PCA-filtered GNSS vertical time series for SCML (a), SCXC (b), YNLJ (c), and YNYS (d) stations.

Figure 8. Distribution of optimal noise models before HYDL correction.

Figure 9. Velocity uncertainty before HYDL correction associated with the WN model and optimal noise model.

Figure 10. Distribution of optimal noise models after HYDL correction provided by EOST (a), GFZ (b), and IMLS (c).

Figure 11. Comparison of CME with HYDL time series provided by EOST at SCML (a), SCXC (b), YNLJ (c), and YNYS (d) stations. The black line represents the CME and the red line represents the EOST_HYDL time series.

Figure 12. Correlation (a) and RMS reduction (b) between CME and HYDL time series provided by EOST at all stations in the Chuandian region.

Figure 13. Absolute value of velocity and its uncertainty difference before and after HYDL correction for all stations in the Chuandian region: (a) represents the absolute value of velocity difference corrected by EOST_HYDL, GFZ_HYDL, and IMLS_HYDL; (b) represents the absolute value of velocity uncertainty difference corrected by EOST_HYDL, GFZ_HYDL, and IMLS_HYDL.

Table 1

The optimal noise model estimated from GNSS vertical time series before and after EOST_HYDL, GFZ_HYDL, and IMLS_HYDL correction.

Stations	Before HYDL Correction	After EOST_HYDL Correction	After GFZ_HYDLCorrection	After IMLS_HYDLCorrection
KMIN	WN + PL	WN + PL	WN + PL	WN + PL
SCDF	WN + FN	WN + PL	WN + PL	WN + PL
SCGZ	WN + PL	WN + PL	WN + PL	WN + PL
SCJL	WN + FN	WN + FN	WN + FN	WN + FN
SCLH	WN + FN	WN + PL	WN + PL	WN + PL
SCMB	WN + FN	WN + FN	WN + FN	WN + FN
SCML	WN + FN	WN + FN	WN + FN	WN + FN
SCMN	WN + FN	WN + PL	WN + PL	WN + PL
SCNN	WN + FN	WN + FN	WN + FN	WN + FN
SCPZ	WN + FN	WN + PL	WN + FN	WN + PL
SCSM	WN + PL	WN + PL	WN + PL	WN + PL
SCXC	WN + FN	WN + FN	WN + FN	WN + FN
SCXD	WN + FN	WN + FN	WN + FN	WN + FN
SCXJ	WN + FN	WN + PL	WN + PL	WN + PL
SCYX	WN + FN	WN + PL	WN + FN	WN + PL
SCYY	WN + FN	WN + FN	WN + FN	WN + FN
XIAG	WN + FN	WN + FN	WN + FN	WN + FN
YNCX	WN + FN	WN + FN	WN + FN	WN + FN
YNDC	WN + FN	WN + FN	WN + FN	WN + FN
YNHZ	WN + FN	WN + FN	WN + FN	WN + FN
YNJD	WN + FN	WN + FN	WN + FN	WN + FN
YNJP	WN + FN	WN + FN	WN + FN	WN + FN
YNLA	WN + GGM	WN + FN	WN + FN	WN + FN
YNLC	WN + FN	WN + FN	WN + FN	WN + FN
YNLJ	WN + FN	WN + FN	WN + FN	WN + FN
YNMJ	WN + FN	WN + FN	WN + FN	WN + FN
YNML	WN + FN	WN + FN	WN + FN	WN + FN
YNMZ	WN + FN	WN + FN	WN + FN	WN + FN
YNRL	WN + FN	WN + FN	WN + FN	WN + FN
YNSD	WN + FN	WN + FN	WN + PL	WN + PL
YNSM	WN + FN	WN + FN	WN + FN	WN + FN
YNTH	WN + FN	WN + PL	WN + FN	WN + PL
YNWS	WN + FN	WN + FN	WN + FN	WN + FN
YNXP	WN + FN	WN + PL	WN + PL	WN + PL
YNYA	WN + FN	WN + FN	WN + FN	WN + FN
YNYL	WN + FN	WN + FN	WN + FN	WN + FN
YNYM	WN + PL	WN + PL	WN + PL	WN + PL
YNYS	WN + FN	WN + PL	WN + FN	WN + PL
YNZD	WN + FN	WN + PL	WN + FN	WN + PL

Table 2

Correlation and RMS reduction between CME and EOST_HYDL time series at all stations in the Chuandian region.

Stations	Correlation	RMS Reduction (%)	Stations	Correlation	RMS Reduction (%)	Stations	Correlation	RMS Reduction (%)
KMIN	0.79	38.05	SCXJ	0.63	19.68	YNML	0.77	33.82
SCDF	0.64	18.90	SCYX	0.71	27.90	YNMZ	0.73	31.96
SCGZ	0.64	20.06	SCYY	0.78	36.92	YNRL	0.72	29.86
SCJL	0.71	27.11	XIAG	0.79	37.94	YNSD	0.74	32.63
SCLH	0.63	20.21	YNCX	0.80	40.26	YNSM	0.77	35.85
SCMB	0.70	28.06	YNDC	0.78	36.93	YNTH	0.78	35.00
SCML	0.76	29.36	YNHZ	0.77	22.28	YNWS	0.71	24.08
SCMN	0.73	31.24	YNJD	0.79	38.58	YNXP	0.80	35.20
SCNN	0.77	32.90	YNJP	0.70	28.96	YNYA	0.80	39.25
SCPZ	0.79	37.71	YNLA	0.76	32.13	YNYL	0.74	30.95
SCSM	0.69	26.80	YNLC	0.77	36.17	YNYM	0.79	39.14
SCXC	0.72	30.99	YNLJ	0.77	36.04	YNYS	0.79	37.86
SCXD	0.73	31.68	YNMJ	0.77	36.61	YNZD	0.74	32.46

References

1. Cheng, P.; Cheng, Y.; Wang, X.; Xu, Y. Update China geodetic coordinate frame considering plate motion. Satell. Navig.; 2021; 2, 2. [DOI: https://dx.doi.org/10.1186/s43020-020-00032-w]

2. Richter, A.; Ivins, E.; Lange, H.; Mendoza, L.; Schröder, L.; Hormaechea, J.L.; Casassa, G.; Marderwald, E.; Fritsche, M.; Perdomo, R. Crustal deformation across the Southern Patagonian Icefield observed by GNSS. Earth Planet. Sci. Lett.; 2016; 452, pp. 206-215. [DOI: https://dx.doi.org/10.1016/j.epsl.2016.07.042]

3. Wang, Q.; Zhang, P.-Z.; Freymueller, J.T.; Bilham, R.; Larson, K.M.; Lai, X.a.; You, X.; Niu, Z.; Wu, J.; Li, Y. Present-day crustal deformation in China constrained by global positioning system measurements. Science; 2001; 294, pp. 574-577. [DOI: https://dx.doi.org/10.1126/science.1063647] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11641493]

4. Wang, W.; Zhao, B.; Wang, Q.; Yang, S. Noise analysis of continuous GPS coordinate time series for CMONOC. Adv. Space Res.; 2012; 49, pp. 943-956. [DOI: https://dx.doi.org/10.1016/j.asr.2011.11.032]

5. Huisman, L.; Teunissen, P.; Hu, C. GNSS precise point positioning in regional reference frames using real-time broadcast corrections. J. Appl. Geod.; 2012; 6, pp. 15-23. [DOI: https://dx.doi.org/10.1515/jag-2011-0006]

6. Legrand, J.; Bergeot, N.; Bruyninx, C.; Wöppelmann, G.; Bouin, M.-N.; Altamimi, Z. Impact of regional reference frame definition on geodynamic interpretations. J. Geodyn.; 2010; 49, pp. 116-122. [DOI: https://dx.doi.org/10.1016/j.jog.2009.10.002]

7. Lyu, M.; Chen, K.; Xue, C.; Zang, N.; Zhang, W.; Wei, G. Overall subshear but locally supershear rupture of the 2021 Mw 7.4 Maduo earthquake from high-rate GNSS waveforms and three-dimensional InSAR deformation. Tectonophysics; 2022; 839, 229542. [DOI: https://dx.doi.org/10.1016/j.tecto.2022.229542]

8. Chen, K.; Avouac, J.-P.; Aati, S.; Milliner, C.; Zheng, F.; Shi, C. Cascading and pulse-like ruptures during the 2019 Ridgecrest earthquakes in the Eastern California Shear Zone. Nat. Commun.; 2020; 11, 22. [DOI: https://dx.doi.org/10.1038/s41467-019-13750-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31911581]

9. Ruhl, C.; Melgar, D.; Grapenthin, R.; Allen, R. The value of real-time GNSS to earthquake early warning. Geophys. Res. Lett.; 2017; 44, pp. 8311-8319. [DOI: https://dx.doi.org/10.1002/2017GL074502]

10. Chen, K.; Ge, M.; Babeyko, A.; Li, X.; Diao, F.; Tu, R. Retrieving real-time co-seismic displacements using GPS/GLONASS: A preliminary report from the September 2015 M w 8.3 Illapel earthquake in Chile. Geophys. J. Int.; 2016; 206, pp. 941-953. [DOI: https://dx.doi.org/10.1093/gji/ggw190]

11. Metivier, L.; Collilieux, X.; Lercier, D.; Altamimi, Z.; Beauducel, F. Global coseismic deformations, GNSS time series analysis, and earthquake scaling laws. J. Geophys. Res. Solid Earth; 2014; 119, pp. 9095-9109. [DOI: https://dx.doi.org/10.1002/2014JB011280]

12. Huang, L.; Mo, Z.; Xie, S.; Liu, L.; Chen, J.; Kang, C.; Wang, S. Spatiotemporal characteristics of GNSS-derived precipitable water vapor during heavy rainfall events in Guilin, China. Satell. Navig.; 2021; 2, 13. [DOI: https://dx.doi.org/10.1186/s43020-021-00046-y]

13. Jiang, Z.; Hsu, Y.-J.; Yuan, L.; Cheng, S.; Feng, W.; Tang, M.; Yang, X. Insights into hydrological drought characteristics using GNSS-inferred large-scale terrestrial water storage deficits. Earth Planet. Sci. Lett.; 2022; 578, 117294. [DOI: https://dx.doi.org/10.1016/j.epsl.2021.117294]

14. White, A.M.; Gardner, W.P.; Borsa, A.A.; Argus, D.F.; Martens, H.R. A review of GNSS/GPS in hydrogeodesy: Hydrologic loading applications and their implications for water resource research. Water Resour. Res.; 2022; 58, e2022WR032078. [DOI: https://dx.doi.org/10.1029/2022WR032078] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36247691]

15. Zhao, Q.; Ma, X.; Yao, W.; Liu, Y.; Du, Z.; Yang, P.; Yao, Y. Improved drought monitoring index using GNSS-derived precipitable water vapor over the loess plateau area. Sensors; 2019; 19, 5566. [DOI: https://dx.doi.org/10.3390/s19245566]

16. Liu, B.; Dai, W.; Liu, N. Extracting seasonal deformations of the Nepal Himalaya region from vertical GPS position time series using Independent Component Analysis. Adv. Space Res.; 2017; 60, pp. 2910-2917. [DOI: https://dx.doi.org/10.1016/j.asr.2017.02.028]

17. Yan, J.; Dong, D.; Bürgmann, R.; Materna, K.; Tan, W.; Peng, Y.; Chen, J. Separation of sources of seasonal uplift in China using independent component analysis of GNSS time series. J. Geophys. Res. Solid Earth; 2019; 124, pp. 11951-11971. [DOI: https://dx.doi.org/10.1029/2019JB018139]

18. Li, W.; Li, Z.; Jiang, W.; Chen, Q.; Zhu, G.; Wang, J. A New Spatial Filtering Algorithm for Noisy and Missing GNSS Position Time Series Using Weighted Expectation Maximization Principal Component Analysis: A Case Study for Regional GNSS Network in Xinjiang Province. Remote Sens.; 2022; 14, 1295. [DOI: https://dx.doi.org/10.3390/rs14051295]

19. Wang, H.; Li, W.; Shu, C.; Shum, C.; Li, F.; Zhang, S.; Zhang, Z. Assessment of spatiotemporal filtering methods towards optimising crustal movement observation network of China (CMONOC) GNSS data processing at different spatial scales. All Earth; 2022; 34, pp. 107-119. [DOI: https://dx.doi.org/10.1080/27669645.2022.2098611]

20. Zhou, W.; Ding, K.; Liu, P.; Lan, G.; Ming, Z. Spatiotemporal Filtering for Continuous GPS Coordinate Time Series in Mainland China by Using Independent Component Analysis. Remote Sens.; 2022; 14, 2904. [DOI: https://dx.doi.org/10.3390/rs14122904]

21. Jiang, W.; Li, Z.; van Dam, T.; Ding, W. Comparative analysis of different environmental loading methods and their impacts on the GPS height time series. J. Geod.; 2013; 87, pp. 687-703. [DOI: https://dx.doi.org/10.1007/s00190-013-0642-3]

22. Liu, B.; Ma, X.; Xing, X.; Tan, J.; Peng, W.; Zhang, L. Quantitative Evaluation of Environmental Loading Products and Thermal Expansion Effect for Correcting GNSS Vertical Coordinate Time Series in Taiwan. Remote Sens.; 2022; 14, 4480. [DOI: https://dx.doi.org/10.3390/rs14184480]

23. Yuan, P.; Li, Z.; Jiang, W.; Ma, Y.; Chen, W.; Sneeuw, N. 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]

24. Mao, A.; Harrison, C.G.; Dixon, T.H. Noise in GPS coordinate time series. J. Geophys. Res. Solid Earth; 1999; 104, pp. 2797-2816. [DOI: https://dx.doi.org/10.1029/1998JB900033]

25. Langbein, J. Noise in GPS displacement measurements from Southern California and Southern Nevada. J. Geophys. Res. Solid Earth; 2008; 113. [DOI: https://dx.doi.org/10.1029/2007JB005247]

26. 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; 10, [DOI: https://dx.doi.org/10.1029/2003JB002741]

27. He, X.; Bos, M.; Montillet, J.; Fernandes, R. 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]

28. Zhang, J.; Bock, Y.; Johnson, H.; Fang, P.; Williams, S.; Genrich, J.; Wdowinski, S.; Behr, J. Southern California Permanent GPS Geodetic Array: Error analysis of daily position estimates and site velocities. J. Geophys. Res. Solid Earth; 1997; 102, pp. 18035-18055. [DOI: https://dx.doi.org/10.1029/97JB01380]

29. Li, Z.; Cao, L.; Jiang, S. Comprehensive analysis of Mass Loading Effects on GPS Station Coordinate Time Series Using Different Hydrological Loading Models. IEEE Access; 2021; [DOI: https://dx.doi.org/10.1109/ACCESS.2021.3067381]

30. He, Y.; Nie, G.; Wu, S.; Li, H. Analysis and discussion on the optimal noise model of global GNSS long-term coordinate series considering hydrological loading. Remote Sens.; 2021; 13, 431. [DOI: https://dx.doi.org/10.3390/rs13030431]

31. Li, Z.; Jiang, W.; Liu, H.; Qu, X. Noise model establishment and analysis of IGS reference station coordinate time series inside China. Acta Geod. Cartogr. Sin.; 2012; 41, pp. 496-503.

32. Liang, S.; Gan, W.; Shen, C.; Xiao, G.; Liu, J.; Chen, W.; Ding, X.; Zhou, D. 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]

33. Zhang, K.; Wang, Y.; Gan, W.; Liang, S. 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]

34. Tan, W.; Chen, J.; Dong, D.; Qu, W.; Xu, X. Analysis of the Potential Contributors to Common Mode Error in Chuandian Region of China. Remote Sens.; 2020; 12, 751. [DOI: https://dx.doi.org/10.3390/rs12050751]

35. Hu, S.; Chen, K.; Zhu, H.; Xue, C.; Wang, T.; Yang, Z.; Zhao, Q. A Comprehensive Analysis of Environmental Loading Effects on Vertical GPS Time Series in Yunnan, Southwest China. Remote Sens.; 2022; 14, 2741. [DOI: https://dx.doi.org/10.3390/rs14122741]

36. 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]

37. 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.

38. 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]

39. Li, C.; Huang, S.; Chen, Q.; Dam, T.v.; Fok, H.S.; Zhao, Q.; Wu, W.; Wang, X. 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]

40. 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]

41. Andrei, C.-O.; Lahtinen, S.; Nordman, M.; Näränen, J.; Koivula, H.; Poutanen, M.; Hyyppä, 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]

42. Hu, S.; Wang, T.; Guan, Y.; Yang, Z. 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.

43. Michel, A.; Santamaría-Gómez, A.; Boy, J.-P.; Perosanz, F.; Loyer, S. Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models. Remote Sens.; 2021; 13, 4523. [DOI: https://dx.doi.org/10.3390/rs13224523]

44. 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]

45. Lyard, F.; Lefevre, F.; Letellier, T.; Francis, O. Modelling the global ocean tides: Modern insights from FES2004. Ocean Dyn.; 2006; 56, pp. 394-415. [DOI: https://dx.doi.org/10.1007/s10236-006-0086-x]

46. McCarthy, D.D.; Petit, G. IERS Conventions (2003); IERS Technical Note No. 32 Publisher of the Federal Agency for Cartography and Geodesy: Frankfurt am Main, Germany, 2004.

47. 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]

48. Li, F.; Li, W.; Zhang, S.; Lei, J.; Zhang, Q.; Yuan, L. Spatiotemporal filtering for regional GNSS network in Antarctic Peninsula using independent component analysis. Chin. J. Geophys.; 2019; 62, pp. 3279-3295.

49. Li, W.; Li, F.; Zhang, S.; Lei, J.; Zhang, Q.; Yuan, L. 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]

50. Gegout, P.; Boy, J.-P.; Hinderer, J.; Ferhat, G. Modeling and observation of loading contribution to time-variable GPS sites positions. Int. Assoc. Geod. Symp.; 2010; 135, pp. 651-659.

51. Dill, R.; Dobslaw, H. 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]

52. Petrov, L. The international mass loading service. Int. Assoc. Geod. Symp.; 2017; 146, pp. 79-83.

53. Berrisford, P.; Dee, D.; Fielding, K.; Fuentes, M.; Kallberg, P.; Kobayashi, S.; Uppala, S. The ERA-interim archive. ERA Rep. Ser.; 2009.

54. Dill, R. Hydrological Model LSDM for Operational Earth Rotation and Gravity Field Variations. Scientific Technical Report STR; GFZ: Potsdam, Germany, 2008.

55. Suarez, M.J.; Rienecker, M.; Todling, R.; Bacmeister, J.; Takacs, L.; Liu, H.; Gu, W.; Sienkiewicz, M.; Koster, R.; Gelaro, R. The GEOS-5 Data Assimilation System-Documentation of Versions 5.0.1, 5.1.0, and 5.2.0; NASA Technical Report Series on Global Modeling and Assimilation, NASA/T-2008-104606, 27; NASA Goddard Space Flight Center: Greenbelt, MA, USA, 2008.

56. Bevis, M.; Brown, A. Trajectory models and reference frames for crustal motion geodesy. J. Geod.; 2014; 88, pp. 283-311. [DOI: https://dx.doi.org/10.1007/s00190-013-0685-5]

57. Bos, M.; Fernandes, R.; Williams, S.; Bastos, L. 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]

58. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control.; 1974; 19, pp. 716-723. [DOI: https://dx.doi.org/10.1109/TAC.1974.1100705]

59. Schwarz, G. Estimating the dimension of a model. Ann. Stat.; 1978; 6, pp. 461-464. [DOI: https://dx.doi.org/10.1214/aos/1176344136]

60. Dong, D.; Fang, P.; Bock, Y.; Webb, F.; Prawirodirdjo, L.; Kedar, S.; Jamason, P. Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis. J. Geophys. Res. Solid Earth; 2006; 111, [DOI: https://dx.doi.org/10.1029/2005JB003806]

61. Ray, J.; Altamimi, Z.; Collilieux, X.; van Dam, T.J.G.s. Anomalous harmonics in the spectra of GPS position estimates. GPS Solut.; 2008; 12, pp. 55-64. [DOI: https://dx.doi.org/10.1007/s10291-007-0067-7]

62. Wdowinski, S.; Bock, Y.; Zhang, J.; Fang, P.; Genrich, J. Southern California permanent GPS geodetic array: Spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake. J. Geophys. Res. Solid Earth; 1997; 102, pp. 18057-18070. [DOI: https://dx.doi.org/10.1029/97JB01378]

63. Li, W.; Shen, Y.; Li, B. Weighted spatiotemporal filtering using principal component analysis for analyzing regional GNSS position time series. Acta Geod. Geophys.; 2015; 50, pp. 419-436. [DOI: https://dx.doi.org/10.1007/s40328-015-0100-1]

64. Shen, Y.; Li, W.; Xu, G.; Li, B. Spatiotemporal filtering of regional GNSS network’s position time series with missing data using principle component analysis. J. Geod.; 2014; 88, pp. 1-12. [DOI: https://dx.doi.org/10.1007/s00190-013-0663-y]

65. He, X.; Yu, K.; Montillet, J.-P.; Xiong, C.; Lu, T.; Zhou, S.; Ma, X.; Cui, H.; Ming, F. GNSS-TS-NRS: An Open-source MATLAB-Based GNSS time series noise reduction software. Remote Sens.; 2020; 12, 3532. [DOI: https://dx.doi.org/10.3390/rs12213532]

66. Wu, S.; Nie, G.; Liu, J.; Xue, C.; Wang, J.; Li, H.; Peng, F. Analysis of deterministic and stochastic models of GPS stations in the crustal movement observation network of China. Adv. Space Res.; 2019; 64, pp. 335-351. [DOI: https://dx.doi.org/10.1016/j.asr.2019.04.032]