*Article* **Analysis and Discussion on the Optimal Noise Model of Global GNSS Long-Term Coordinate Series Considering Hydrological Loading**

**Yuefan He 1, Guigen Nie 1,2,\*, Shuguang Wu <sup>1</sup> and Haiyang Li <sup>1</sup>**


**Abstract:** The displacement of Global Navigation Satellite System (GNSS) station contains the information of surface elastic deformation caused by the variation of land water reserves. This paper selects the long-term coordinate series data of 671 International GNSS Service (IGS) reference stations distributed globally under the framework of World Geodetic System 1984 (WGS84) from 2000 to 2021. Different noise model combinations are used for noise analysis, and the optimal noise model for each station before and after hydrologic loading correction is calculated. The results show that the noise models of global IGS reference stations are diverse, and each component has different optimal noise model characteristics, mainly white noise + flicker noise (WN+FN), generalized Gauss–Markov noise (GGM) and white noise + power law noise (WN+PL). Through specific analysis between the optimal noise model and the time series velocity of the station, it is found that the maximum influence value of the vertical velocity can reach 1.8 mm when hydrological loading is considered. Different complex noise models also have a certain influence on the linear velocity and velocity uncertainty of the station. Among them, the influence of white noise + random walking noise is relatively obvious, and its maximum influence value in the elevation direction can reach over 2 mm/year. When studying the impact of hydrological loading correction on the periodicity of the coordinate series, it is concluded whether the hydrological loading is calculated or not, and the GNSS long-term coordinate series has obvious annual and semi-annual amplitude changes, which are most obvious in the vertical direction, up to 16.48 mm.

**Keywords:** hydrological loading; long-term coordinate series; noise model; velocity; amplitude

#### **1. Introduction**

The global IGS base station location time series accumulated in the past 20 years have provided valuable basic data for geodesy and geodynamic research. Many scholars have also launched a lot of research on it, including periodic signal analysis, stochastic model research, spatial filtering method and surface loading model research, etc., in GNSS coordinate time series [1–16]. Analyzing it can get the precise movement trend of the station, so as to further study the internal driving mechanism of the station movement. It has far-reaching theoretical significance and broad application prospects to explain plate tectonic movement, establish and maintain a dynamic earth reference frame, study post-ice rebound and sea level changes, and invert the changes in ice and snow quality.

The noise of IGS station mainly comes from some random factors. As the stations are located all over the world, the geographical environment is very different and the sources of noise influence are also different. For example, different types of observation piers, missing time series data, time series of different spans, different geographical environments, etc. may generate some noise [17]. Currently, a large number of researchers generally believe that the characteristics of GNSS coordinate series noise model are WN+FN [18–20].

**Citation:** 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. https:// doi.org/10.3390/rs13030431

Academic Editor: Alex Hay-Man Ng Received: 14 December 2020 Accepted: 20 January 2021 Published: 26 January 2021

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

However, it turns out that this statement is not entirely true. The noise model characteristics of the reference station are actually more complicated. Langbein [3] estimated the noise model of 236 continuous GPS operation stations in California and southern Nevada in 2008; the results showed that the best noise model of FN or random walk noise (RW) station accounted for 50% to 60%, 25%~30% of the station performance are FN+RW or noninteger spectral exponential powerlaw noise (PL). The remaining stations are expressed as band pass+powerlaw noise (BPPL) and first-order Gauss–Markov+random walk noise (FOGMRW). In 2019, Wang et al. [21] analyzed the coordinate time series of 260 continuous GPS stations in the Crustal Movement Observation Network of China (CMONOC). In the noise evaluation, the influence of periodic signals is considered, and the maximum likelihood estimation method is used to discuss the noise characteristics of the residual time series that remove seasonal signals and the coordinate time series that further remove other periodic signals. It is concluded that FN is the main noise related to time. If it is assumed to be pure white noise, its speed uncertainty will be underestimated by 8-10 times. In 2008, Yuan et al. [22] analyzed the noise characteristics of 12 GPS reference stations in Hong Kong. It was concluded that the noise model after the common mode error (CME) could be removed by the principal component, and spatial filtering can be described by the variable white noise (VW) plus flicker noise. In 2010, Jiang et al. analyzed the coordinate series data of the China Continuously Operating Reference Stations (CORS) network in the past 10 years. The results point out that under the framework of the China Geodetic Coordinate System 2000 (CGCS2000), the noise model of coordinate time series appears as WN, FN and RW after principal component filtering. Additionally, if only white noise is considered, the velocity error will be 2~6 times smaller than when colored noise is considered [23].

The reason why WN+FN is considered to be the noise model representing the optimal random characteristics of the station may be analyzed in two ways: (1) When analyzing the time series, a more complex noise model that sufficiently represents the noise characteristics of the reference station is not used; (2) The time scale of the accumulated coordinate series is not long enough to explain the long-period component of the noise model. With the passage of time, the GNSS reference station coordinate time series continue to grow, and the long-period components of noise (such as RW noise with spectral index=2) will become more prominent, which provides favorable conditions for detecting the existence of lowfrequency noise. Therefore, it is meaningful to conduct a more comprehensive noise analysis on GNSS time series and obtain a noise model that can more accurately represent the random characteristics of the reference station.

The GNSS coordinate series contains various ground loading changes related to the elastic deformation of Earth, and is closely related to changes in the atmosphere, seabed pressure, snow, ice, surface water and groundwater [24,25]. Some scholars have deeply analyzed the spatial and temporal distribution characteristics of the vertical displacements of IGS stations caused by the above environmental loadings. Their influence on the nonlinear variation and common mode error (CME) of GPS coordinate time series are estimated. The results show that the average root mean square (RMS) of vertical displacement caused by environmental loadings is 4.0 mm. The impact on one station can reach centimeter level [26]. Collilieux et al. [27] conducted the study based on 748 stations around the world, with many stations concentrated in North America and Europe. They found that when the time series were corrected for loading, the in-phase and out-of-phase coefficients of the annual signal in GPS height coordinates were reduced by 84% and 83 %, respectively. These factors may have a profound influence on the characteristics of the generated noise [28,29]. Among them, the land hydrological loading change is an important factor that causes the periodic vertical deformation of GNSS stations [30–32]. The increase of the hydrological loading on the land will cause the ground to sink, causing the vertical position of the GNSS station to move downward. Decreasing the hydrological loading on the land will cause the surface to rebound, causing the vertical position of the GNSS station to move upward.

The main innovation of this article is to use longer GNSS coordinate time series, adopt more types of noise model combinations and more IGS reference stations, and analyze the optimal noise model characteristics of global regional coordinate time series. At the same time, the displacement of the station caused by the hydrological loading is calculated, and the change analysis of the optimal noise model before and after the hydrological loading correction is given. On this basis, the relationship between the optimal noise model and the velocity, velocity uncertainty and amplitude of the stations are established. The organization structure of this article is as follows. The second part contains GNSS, hydrological data and theoretical methods to be analyzed in this article. The third part analyzes the optimal noise model of the selected station, including the comparison of the optimal noise model before and after hydrological loading correction. In the fourth part, the velocity and velocity uncertainty of the stations and the amplitude of the long-term coordinate series are discussed with different noise models before and after the hydrological loading correction. Finally, the fifth part summarizes the main results.

#### **2. Data and Methods**

#### *2.1. GNSS Data*

The GNSS time series data used in this article are released by the Scripps Orbit and Permanent Array Center (SOPAC), a data analysis center organized by the IGS [33]. This organization not only publishes some raw GNSS time series results, but also publishes its processed (Clean) time series products. In this article, we mainly use the linearized mean GNSS time series (CleanDetrend) after removing the step (coseismic and non-coseismic) and singular values. Seasonal period term is retained. In order to analyze the impact of various noises on it and explain the results of the global optimal noise model after adding hydrological loading.

Before data analysis, this paper needs to filter the long-term coordinate series data of 671 IGS reference stations distributed globally under the framework of WGS84 from 2000 to 2021. The screening criteria are as follows: first, within the 21 years, the length of the GNSS coordinate series shall not be less than 3 years to ensure the integrity of the long time series and the reliability of the data processing results. Then, delete sites with obvious abnormal nonlinear motion, including obvious post-earthquake deformation caused by a major earthquake, local surface deformation or abnormal motion caused by other unknown causes. Through the elimination of the above steps, 671 sites around the world are finally selected for data processing. The time length distribution of the selected GNSS stations in the past 21 years are shown in Figure 1. The longest time series can reach 20.54 years, the shortest time length is 3.47 years, the average time length can reach 13.25 years, and the average data integrity rate can be up to 63.1%. In this article, the average data integrity rate equals average time series length of selected sites divided by the total length of GNSS coordinate time series (21 years).

#### *2.2. Hydrological Loading Data*

Environmental loadings including hydrological loading, atmospheric loading, nontidal ocean loading and other environmental factors are closely related to the elastic deformation of the earth, and these factors will also have a significant impact on the noise characteristics generated in the GNSS time series [25,28]. The impact of this loadings are not considered in the GNSS coordinate time series data processing under the WGS84 framework provided by SOPAC. Among them, the change of hydrological loading is an important factor that causes periodic vertical deformation of GNSS stations, as well as an important cause of surface subsidence and rebound. Therefore, it is particularly important to study the influence of hydrological loading on the GNSS coordinates time series [30,31].

**Figure 1.** Distribution of 671 giobal GNSS reference stations and their coordinate time series length in 21 years (2000–2021).

At present, many organizations provide hydrological loading products, for example, German Research Centre for Geosciences (GFZ) (http://rz-vm115.gfz-potsdam.de:8080 /repository), International Mass Loading Service (IMLS) (http://massloading.net) products released by National Aeronautics and Space Administration (NASA) and School and Observatory of Earth Sciences (EOST) (http://loading.u-strasbg.fr/) loading products from university of Strasbourg. This article first selects the hydrological loading (HYDL) data products provided by GFZ [25,34] and estimate the displacement of global GNSS stations caused by hydrological loading, and subtract the influence of hydrological loading from the GNSS coordinate series. Then, analyze the optimal noise model of the global GNSS coordinate series taking into account the hydrological loading. Finally, based on the optimal noise model combination, we will further discuss its influence on the long-term coordinate series velocity, velocity uncertainty and amplitude of the IGS stations.

The hydrological loading data selected in this paper are based on the center of figure (CF) frame, with a time resolution of 24 h and a spatial resolution of 0.5◦ × 0.5◦. The hydrological loading data model calculated by GFZ is the Land Surface Discharge Model (LSDM) [35]. The physics and parameterization of the LSDM is based on the research of Hagemann and Düumenil (1998). LSDM includes soil moisture, snow cover, shallow groundwater, and surface water stored in rivers and lakes [36].

#### *2.3. Theoretical Method*

The GNSS coordinate sequence is represented by two parts: a deterministic model and a stochastic model [37]. The deterministic model is composed of trend, periodic terms (including anniversaries, semi-anniversaries, etc.), offset terms, etc. The stochastic model is the noise that this article will analyze. The GNSS coordinate time series is represented by the following functional model.

$$y(t\_i) = y\_0 + vt\_i + \sum\_{k=1}^{q} a\_k \sin(2\pi f\_k t\_i) + b\_k \cos(2\pi f\_k t\_i) + \sum\_{j=1}^{n\_g} g\_j H(t\_i - T\_{\xi j}) + r(t\_i) \tag{1}$$

In the above formula, *y*<sup>0</sup> is the intercept; *v* is the linear trend in unit per year, where a year is defined as 365.25 days; *ti* is the epoch of the GNSS time series; *ak*, and *bk* are the amplitude of periodic signals and *fk* is the corresponding frequency; *gj* and *Tgj* is the offset and corresponding epoch, respectively; *r* refers to the residual time series, which can be various combinations of noise models. *H* is the Heaviside step function:

$$\begin{cases} \quad H(t) = 0, t < 0 \\ \quad H(t) = 0.5, t = 0 \\ \quad H(t) = 1, t > 0 \end{cases} \tag{2}$$

In this article, we use the Hector software to analyze the optimal noise model of the long-term coordinate series of the IGS reference station [38]. Hector software can estimate linear trend terms, high-order polynomials, seasonal terms, periodic terms and a variety of noise model combinations in coordinate time series. It uses the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) information standards to select the optimal noise model [39,40]. Both use the maximum likelihood estimate as a starting point, but in order to avoid overfitting, measures are taken to add parameters. The AIC and BIC information standards are expressed by the following formula.

$$\ln(L) = -\frac{1}{2} \left[ Nln(2\pi) + ln det(\mathbb{C}) + r^T \mathbb{C}^{-1} r \right] \tag{3}$$

In the above formula, *N* represents the actual number of observations, and *C* refers to the covariance matrix, which can be decomposed into the following formula.

$$\mathcal{C} = \sigma^2 \overline{\mathcal{C}} \tag{4}$$

*C* represents the sum of multiple noise models, and σ represents the standard deviation, which can be estimated by the residual term in the following formula.

$$
\sigma = \sqrt{\frac{r^T \overline{C}^{-1} r}{N}} \tag{5}
$$

Due to det*cA* = *cNdetA*, Where *A* refers to any matrix, *detA* refers to the determinant of matrix *A*. Therefore, you can get

$$\ln(L) = -\frac{1}{2}\left[\mathrm{Nl}\ln(2\pi) + \mathrm{lndet}(\overline{\mathsf{C}}) + 2\mathrm{Nl}\ln(\sigma) + \mathrm{N}\right] \tag{6}$$

The parameter k refers to the sum of a design matrix and the noise model parameters and the variance of the driving white noise. For example, if the white noise + power law noise model is used to estimate the linear trend term, 5 parameters are estimated, including the normal offset, the linear trend term, the power spectrum index, the variance of driving white noise and the difference between power law noise and white noise (k = 2 + 2 + 1 = 5).

Then,

$$AIC = 2k + 2ln(L)\tag{7}$$

$$BIC = kln(N) + 2ln(L)\tag{8}$$

AIC and BIC are commonly used optimal noise model evaluation indicators in the maximum likelihood estimation analysis of Hector software [41]. From the two calculation formulae, it can be found that the better noise model should have a smaller AIC or BIC value. Since BIC considers the number of parameters to be estimated in the model, it is more suitable for comparing noise models with different numbers of parameters to be estimated. In this article, the optimal noise model is obtained by using BIC evaluation criteria. There are many noise models to choose from in the Hector software, including PL, WN, FN, RW, GGM, first-order autogressive noise (AR1) and a combination model of these noises.

#### **3. Results**

#### *3.1. Before Hydrological Loading Correction*

As far as the global GNSS single-day solution time series are concerned, if no measures are taken to reduce the spatial correlation of the station, it can be concluded that WN + FN is the most suitable stochastic noise model reflecting global GNSS coordinates [2]. Due to the different environments of each measuring station, the sources of noise may also be different, and their noise characteristics may not be exactly the same. This paper selects six noise models, namely WN, GGM, WN+FN, WN+PL, WN+RW and WN+GGM. According to the evaluation criteria of the optimal noise model in Section 2.3 above, use Hector software to obtain the minimal BIC value corresponding to the N, E, and U directions of 671 sites around the world, then the noise model corresponding to this minimal BIC value is the optimal noise model of the selected site. Finally, statistics are made on the optimal noise models corresponding to 671 sites around the world, and the proportion of the global optimal noise models before the hydrological loading correction is obtained, as shown in Figure 2 below.

**Figure 2.** The optimal noise model distribution in the North (N), East (E), Up (U) directions and total components of global International Global Navigation Satellite System (GNSS) Service (IGS) stations before hydrological loading correction. (**a**), (**b**) and (**c**) respectively show the optimal noise model percentage in the N, E and U directions. The optimal noise model percentage of total components shown in (**d**).

It can be seen from Figure 2 that the noise characteristics of global IGS reference station components are diversified, mainly represented by WN+FN and GGM models. Especially in the N and E directions, the WN+FN model accounts for 51.3% and 53.4% of the world's optimal noise models, respectively. In the U direction, the GGM model has the largest proportion of the optimal noise model, which is as high as 51.6%. Next is the WN+PL noise model, and other combined noise models also occupy a certain proportion. The N, E, and U components of most stations exhibit different noise characteristics. As shown in Figure 2d above, the total components refer to the sum of the three directions of N, E, and U. Figure 2d shows the accumulation of the stations in these three directions to obtain

the total optimal noise model percentage. For example, WN+FN, which accounts for the largest proportion of the optimal noise model, can reach 46.4%. It defines to the sum of the number of stations where the optimal noise model in the three directions of N, E, and U is WN+FN divided by the total number of stations in three directions (671\*3). Therefore, when looking for a physical explanation for the change characteristics of time series, each component needs to be treated differently.

#### *3.2. Hydrological Loading Displacement*

This paper uses the hydrological loading products of GFZ to obtain the hydrological loading displacement values of the long-term coordinate time series of 671 IGS stations worldwide from 2000 to 2021. First, we selected three IGS stations at different latitudes, namely ACP1 (longitude 280.0501, latitude 9.3714), ALGO (longitude 281.9286, latitude 45.9558), ALRT (longitude 297.6595, latitude 82.4943). Then, we drew the displacement series diagram of their hydrological loading from 2000 to 2021. As shown in Figure 3, in the horizontal direction, the hydrological loading displacement values of these three stations are all small, and the maximum value is less than 2 mm. Especially for the station of ALRT, the maximum horizontal hydrological loading displacement is only 0.93 mm. However, it can be seen in Figure 3 that the hydrological loading influence value of these three stations is larger in the vertical direction. For the low-latitude station ACP1, the maximum vertical hydrological loading can reach 6.46 mm, and there is a periodic fluctuation in the long-term series. For the mid-latitude site ALGO, the maximum vertical hydrological loading can reach 9.48 mm, and the long-term series shows that its hydrological loading influence value has a downward trend. For the high-latitude site ALRT, the maximum vertical hydrological loading can reach 6.29 mm, but it has an upward trend in the long-term series. Regardless of whether the vertical hydrological loading series of these three stations rises, falls, or changes periodically, it may be related to the surrounding environment, including rainfall, soil moisture, etc.

**Figure 3.** Distribution diagram of hydrological loading displacement series at the three IGS stations ACP1, ALGO, and ALRT from 2000 to 2021.

Since the time span is as long as 21 years, this article counts the RMS value in the N, E, and U directions of the selected 671 IGS reference stations from 2000 to 2021 due to hydrological loading. Then, the maximum RMS value of the three directions of N, E, and U is represented by a histogram. As shown in Figure 4, the magnitude of the hydrological loading in the N and E directions are relatively small. Among them, the RMS value of more than 90% of the stations are less than 1 mm, the maximum RMS value in the E direction can reach 2.07 mm, and the maximum RMS value in the N direction can reach 2.15 mm. In the vertical direction, the magnitude of the influence of hydrological loading is relatively large, and the RMS value can reach 13.48 mm. Especially in the Amazon River Basin in South America, the Mississippi River Basin in North America, the Congo River Basin in Africa, the Yangtze River Basin in China and the Indochina Peninsula, the magnitude of the hydrological load is relatively large. Therefore, when considering the influence factors of the displacement of the reference frame, the hydrological loading must be considered.

**Figure 4.** A statistical diagram of the RMS value of hydrological loading displacement in the N, E, and U directions and their maximum RMS values. (**a**), (**b**) and (**c**) respectively show the RMS value of hydrological loading displacement in the N, E, and U directions and their maximum RMS values are shown in (**d**).

#### *3.3. After Hydrological Loading Correction*

Research shows that environmental loading, including hydrological loading, atmospheric loading and non-tidal ocean loading are the main factors that cause non-linear changes in stations [42]. Among them, the hydrological loading change is an important reason for the periodic vertical changes of GNSS stations. The GNSS time series itself contains the influence of loading. This article deducts it on the basis of it, and then analyzes the influence of hydrological loading on GNSS time series. Therefore, based on the above calculation results, this paper corrects the hydrological loading of these 671 IGS reference stations, re-analyzes the noise of the corrected coordinate time series, and obtains the optimal noise model of the global IGS reference station N, E, U components after deducting the influence of hydrological loading. Figure 5 shows the proportional distribution of the optimal noise model for each station component.

**Figure 5.** The optimal noise model distribution in the N, E, U directions and total components of global IGS stations after hydrological loading correction. (**a**), (**b**) and (**c**) respectively show the optimal noise model percentage in the N, E and U directions and the optimal noise model percentage of total components are shown in (**d**).

From Figure 5a–d, it can be seen that, compared with Figure 2a–d, after calculating the hydrological loading correction, the percentage of the WN+FN model in the N, E, and U directions obviously increases, and the percentage of GGM model somewhat decreases. However, in the N and E directions, the optimal noise model, WN+FN, still accounts for the largest proportion with 55.4% and 56.8%, respectively. Compared with the case without employing hydrological loading correction, the percentage increased by 4.1% and 3.4%, respectively. The second largest proportion is WN+PL. In the U direction, the GGM model still has the largest proportion of the optimal noise model, accounting for 47.6%, followed by the WN+PL model. Regardless of whether the impact of hydrological loading is considered, in the long-term coordinate series analysis of global IGS reference stations, the smaller proportion of the optimal noise model are the combination of WN+GGM and WN+RW models.

#### **4. Discussion**

#### *4.1. Velocity and Velocity Uncertainty Analysis*

Studies have shown that the estimated values of velocity and velocity uncertainty in GNSS time series will vary with the noise model used [43–45]. Therefore, in the velocity error estimation, the influence of different noise models should be considered. According to the above analysis, in the global GNSS long-term coordinate series, WN+FN is the optimal noise model. The results obtained in this article are the same as those obtained by many previous scholars [46,47]. In this article, we used the CleanDetrend data provided by SOPAC, and the trend item information has been deduced from the GNSS time series. We added the optimal noise model combination WN+FN, using Hector software to extract the velocity of the selected 671 IGS long-term coordinate series. The uncertainty of the velocity is also generated immediately, which actually represents the fluctuation of the

velocity value. When considering the optimal noise model combination WN+FN, compare the difference of the velocity and its uncertainty relative to SOPAC. That is, when only WN + FN is considered, the velocity and velocity uncertainty difference distribution of each station components are shown in Figure 6 below:

**Figure 6.** The difference between the velocity and its uncertainty estimation under the optimal noise model of global IGS stations. (**a**,**c**,**e**) show the difference between the velocity under the optimal noise model of global IGS stations. (**b**,**d**,**f**) show the difference between the velocity uncertainty estimation under the optimal noise model of global IGS stations.

In order to more clearly and intuitively see the specific size of the velocity difference and the velocity uncertainty difference of each direction in each interval, this article lists the percentage of the interval of velocity and velocity uncertainty of each station. As shown in Table 1, in the direction of N, E and U, the absolute value of the velocity difference is more than 85% within 1 mm, and only less than 10% of the stations have the absolute value of the velocity difference between 1 and 2 mm, and some of them are larger than 2 mm. As shown in Figure 6a,c,e, these sites over 2 mm are marked with black and white dots. We can see that most of the black and white dots are near the coastal area. For example, ANTC (longitude 288.4680◦, latitude −37.3387◦), MIZU (longitude 141.1328◦, latitude 39.1352◦) stations, etc. The reason for analysis may be related to the imperfect simulation of noise models in these coastal areas. For the velocity uncertainty difference, in the N and E directions, 86% and 87% of the stations velocity uncertainty differences are less than 0.4 mm/year, and the maximum value is less than 1 mm. In the U direction, the influence value is slightly larger than that in the N and E directions. A total of 87% of the stations velocity uncertainty differences are less than 0.8 mm, and the maximum value is less than 2 mm. If you consider the effects of rounding errors, the calculated results can be considered reasonable.

**Table 1.** Distribution of station velocity and velocity uncertainty difference.


In order to further study the influence of hydrological loading on GNSS long-time series velocity and its uncertainty, considering only the WN+FN noise, this paper analyzes the differences in the velocities and its uncertainties of the 671 IGS stations before and after the hydrological loading correction, and draws the percentage diagrams in the N, E, and U directions. The result shows that it has basically no effect on the velocity uncertainty, but it will cause varying degrees of velocity changes at the station, especially in the vertical direction. The statistical results are shown in Figure 7. It can be seen that in the N direction, after the hydrological loading correction, the velocity influence value can reach up to 0.9 mm, and 97.4% of the values are between 0 and 0.2 mm. In the E direction, the maximum influence value can reach 0.5 mm, and 89.8% of the influence values are between 0 and 0.1 mm. In the U direction, we can see that the influence value is more obvious relative to the horizontal direction. Its maximum value can reach 1.8 mm, most of the sites are between 0 and 0.4 mm, and the percentage can account for 82.6%. Therefore, when estimating the vertical velocity, the influence of hydrological loading must be considered.

**Figure 7.** The percentage diagram of velocity influence value before and after hydrological loading correction.

Compared with the WN+FN model, different complex noise models will also have a certain impact on the station velocity and its uncertainty, as shown in Figure 8 below, the vertical direction is larger and the horizontal direction is smaller. Among them, regarding the influence of velocity, in the directions of N and E, the remaining four noise models are basically the same, and more than 85% of the station velocity differences are within 1 mm in absolute value. However, in the U direction, the ratio of the WN+RW velocity difference between (−1, 1) is significantly smaller than other noise models, accounting for only 68.7%, while in (−2, 1) and (1, 2), the proportion between them is slightly larger. The difference between these four noise models is mainly reflected in the velocity uncertainty. For the velocity uncertainty difference, in the directions of N, E, and U, the GGM noise model accounts for the largest proportion between (0, 0.2), accounted for 76.5%, 79.1%, 64.5%, respectively, followed by WN+GGM, WN+PL noise model combination. Among them, WN+RW has the largest impact level, and the velocity uncertainty difference of more than 85% of the stations are greater than 1 mm/year. Especially in the U direction, the velocity uncertainty difference is greater, and more than 95% of the stations are greater than 2 mm/year. It can be concluded that for the establishment of mm-level high-precision reference frames and plate motion analysis, we need to take into account the differences brought about by different noise models.

#### *4.2. Amplitude Analysis*

To study the relationship between the hydrological loading and the amplitude of the selected global IGS reference stations. In this paper, the annual and half-annual amplitudes of the stations considering the influence of hydrological loading under different noise models are estimated. The results under the optimal noise model are shown in Figure 9. It can be seen that no matter whether the hydrological loading is corrected or not, there are significant annual and semi-annual amplitudes at each station. Without considering the influence of hydrological loading, the U direction is the largest, and its maximum annual amplitude can reach 16.48 mm, and 95.7% of the stations have an annual amplitude between 0 and 8 mm. The N and E directions are next. In the N direction, the maximum annual amplitude can reach 9.87 mm, and in the E direction, the maximum annual amplitude can reach 10.07 mm. A total of 92.1% and 91.8% of the station sizes are between 0 and 2 mm. In addition to the influence of the annual amplitude, the half-annual amplitude also has a significant change, which is still more affected in the U direction, and its maximum value can reach 5.11 mm. Compared with the vertical direction, the horizontal half-annual amplitude is smaller. In the N and E directions, 93.4% and 90.6% of the stations are less than 1 mm.

After the hydrological loading is corrected, when only the WN+FN noise model is considered, it will have a certain impact on the annual and half-annual amplitude of the stations. Figure 10 plots the annual and half-annual amplitude differences before and after the hydrological loading correction in each direction. It can be seen that after the hydrological loading correction is added, it still has a relatively large impact on the U direction. The maximum annual amplitude difference can reach 8.08 mm, and the annual amplitude difference of most stations are between −3 and 3 mm, which can account for 97.5%. Compared with the annual amplitude, the half-annual amplitude difference is smaller. In the U direction, the maximum value can reach 1.71 mm, and the absolute value of the half-annual amplitude difference of 98.6% of the stations are less than 1 mm. In the horizontal direction, whether it is the annual amplitude or the half-annual amplitude, its influence value is relatively small. In the N direction, the maximum annual amplitude difference is 1.04 mm, and the absolute value of the annual amplitude difference in 99.7% of the stations are within 1 mm. The maximum half-annual amplitude difference is 0.21 mm, and the absolute value of the half-annual amplitude difference in 98.3% of the stations are less than 0.2 mm. In the E direction, the maximum annual amplitude difference is 1.01 mm, and the absolute value of the annual amplitude difference in 99.9% of the stations are within 1 mm. The maximum half-annual amplitude difference is 0.32 mm, and the

absolute value of the half-annual amplitude difference in 99.4% of the stations are less than 0.2 mm. From this, it can be concluded that the calculated hydrological loading will indeed cause the annual and semi-annual movements of IGS reference stations in the global region. The impact varies from station to station and is related to the geographical environment around the stations. The annual amplitude is larger than the half-annual amplitude, but it cannot completely reduce the annual and half-annual amplitudes of the GNSS long time coordinate series, especially the half-annual amplitude.

**Figure 8.** The distribution of velocity and velocity uncertainty differences under different noise models. (**a**) shows the distribution of velocity differences under different noise models and (**b**) shows the distribution of velocity uncertainty differences under different noise models.

**Figure 9.** Distribution of the annual amplitudes before hydrological loading correction are shown (**a**,**c**,**e**). (**b**,**d**,**f**) show the half-annual amplitude before hydrological loading correction. (unit: mm).

In order to more clearly explain the influence of the hydrological loading correction on the annual and half-annual amplitudes, Table 2 lists the percentages of the annual and half-annual amplitude differences in each interval in the N, E, and U directions in detail.

Figure 11 shows the statistics of the annual and half-annual amplitude difference before and after the hydrological loading correction in various directions under different noise models. It can be seen, no matter which noise model is selected before and after the hydrological loading correction, that the influence of the annual amplitude and half-annual amplitude is larger in the vertical direction and smaller in the horizontal direction. Therefore, it can

be concluded when estimating the influence of the hydrological loading on the annual and semi-annual amplitude difference of the IGS station that the noise models including WN+FN, WN+PL, WN+RW, GGM and other noise models have negligible influence.

**Figure 10.** Statistics of annual and half-annual amplitude differences before and after hydrological loading correction. (**a**,**c**,**e**) show the statistics of annual amplitude differences before and after hydrological loading correction. The half-annual amplitude differences before and after hydrological loading correction are shown in (**b**,**d**,**f**). (unit: mm).


**Table 2.** Annual and half-annual amplitude difference percentage table before and after hydrological loading correction.

**Figure 11.** Statistics of annual and half-annual amplitude difference under different noise models before and after hydrological loading correction. (**a**) shows the statistics of annual amplitude difference under different noise models before and after hydrological loading correction and (**b**) shows the half-annual amplitude difference under different noise models before and after hydrological loading correction. (unit: mm).

#### **5. Conclusions**

In this paper, we used the GNSS coordinate time series of 671 IGS stations around the world from 2000 to 2021 to analyze the world's optimal noise model. More importantly, the influence of hydrological loading on the GNSS coordinates time series is discussed, and the following conclusions are obtained:


However, there are currently many research institutions that have calculated the hydrological loading. In this article, only the hydrological loading products of GFZ are discussed, and other hydrological loading products are not explained. Therefore, in future research, multiple hydrological loading models should be compared and analyzed, and the influence of hydrological loading on the GNSS coordinates time series should be studied at a deeper level.

**Author Contributions:** Conceptualization, Y.H. and G.N.; methodology, Y.H. and S.W.; software, Y.H.; validation, G.N. and H.L.; formal analysis, H.L.; resources, G.N.; writing—original draft preparation, Y.H.; writing—review and editing, S.W.; visualization, H.L.; supervision, G.N.; funding acquisition, G.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the National Natural Science Foundation of China (Grant No. 41074023). National key Research and Development Scheme Strategic International Coopearation in Science and Technology Innovtion Program, and the funding number is No.2018YFE0206500.

**Acknowledgments:** The 671 long-time coordinate series calculated in this article from 2000 to 2021 were provided by SOPAC (ftp://garner.ucsd.edu/archive/garner/timeseries/measures/ats/). Hydrological loading data were provided by the GFZ agency from Germany (http://rz-vm115 .gfz-potsdam.de:8080/ repository). The Hector software (http://segal.ubi.pt/hector/) was used to calculate and analyze GNSS long-term coordinate series. We used GMT software and Origin software to plot the calculation results. We express our heartfelt thanks to these organizations and software providers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**

