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

Abstract

Mass redistribution of the atmosphere, oceans, and terrestrial water storage generates crustal displacements which can be predicted by environmental loading models and observed by the Global Positioning System (GPS). In this paper, daily height time series of 235 GPS stations derived from a homogeneously reprocessed Crustal Movement Observation Network of China (CMONOC) and corresponding loading displacements predicted by the Deutsche GeoForschungsZentrum (GFZ) are compared to assess the effects of loading corrections on the nonlinear variations of GPS time series. Results show that the average root mean square (RMS) of vertical displacements due to atmospheric, nontidal oceanic, hydrological, and their combined effects are 3.2, 0.6, 2.7, and 4.0 mm, respectively. Vertical annual signals of loading and GPS are consistent in amplitude but different in phase systematically. The average correlation coefficient between loading and GPS height time series is 0.6. RMS of the GPS height time series are reduced by 20% on average. Moreover, an investigation of 208 CMONOC stations with observing time spans of ~4.6 years shows that environmental loading corrections lead to an overestimation of the GPS velocity uncertainty by about 1.4 times on average. Nevertheless, by using a common mode component filter through principal component analysis, the dilution of velocity precision due to environmental loading corrections can be compensated.

1. Introduction

Global Positioning System (GPS) height time series record complex nonlinear variations, which have been demonstrated to be caused by GPS errors, environmental loading, and other effects [1]. Therefore, a comprehensive understanding of the nature and characteristics of the nonlinear variations of GPS height time series is important to the investigations of geophysical phenomena [2], the analysis of GPS error sources [3], and the modeling for GPS height time series [4]. Moreover, an accurate removal of the nonlinear variations from the GPS height time series is vital for the accuracy of secular velocity estimates [5,6,7,8].

GPS errors lead to artificial displacements of GPS stations. The GPS errors are caused by a deficiency of the models and strategies used in GPS data processing, such as the tropospheric delay model [9,10], high-order ionospheric delay model [11,12], unmodeled daily and subdaily tides [13,14], GPS draconitic signal [3], multipath and geometry effects [15], earth radiation pressure model [16] reference frame definition [17], etc.

On the contrary, environmental loading leads to real displacements of the crust. The main sources of environmental loading that displace the earth’s surface are the mass redistribution of the atmosphere, oceans, and terrestrial water storage. A number of studies have quantified the contributions of environmental loading displacements based on surface loading models (SLMs) and evaluated their influences on the nonlinear displacements of the GPS stations. With regard to atmospheric loading (ATML), van Dam et al. [18] estimated its influences on GPS height time series. In that paper, the authors found that the variance of GPS height time series could be reduced with a maximum ratio of 24% by accounting for ATML. Petrov and Boy [19] as well as Tregoning and van Dam [20] found that displacements of the earth’s surface caused by ATML could reach 20 mm and 18 mm, respectively. Later, van Dam et al. [21] reported that the nonlinear vertical displacements of GPS stations could be better explained by ATML displacements after consideration of topography. As for nontidal oceanic loading (NTOL), van Dam et al. [22] found that it could cause vertical displacements with a maximum root mean square (RMS) of 2 mm. Williams and Penna [23] investigated NTOL effects on 17 GPS stations close to the southern North Sea and found that NTOL effects are comparable in size with ATML effects. Later, van Dam et al. [24] showed that the scatters of height time series are mitigated for 70% of stations when considering NTOL. With regard to hydrological loading (HYDL), van Dam et al. [25] investigated the height time series of 147 GPS stations and found significant reductions of scatter at 92 stations after HYDL correction. Fritsche et al. [26] performed a quantitative analysis of 208 stations worldwide. They found that displacements due to HYDL predicted with the WaterGAP Global Hydrology Model (WGHM) are slightly overestimated on average. Dill and Dobslaw [27] found that the scatter of GPS height time series with ATML and NTOL effects corrected can be explained by HYDL up to 54%. Moreover, the nonlinear variations of GPS height time series are usually characterized by obvious seasonal (annual and semiannual) signals. Dong et al. [1] found the vertical annual amplitudes of GPS height time series can be partially explained by environmental loading for 90 out of 128 globally distributed GPS stations. Recently, the impact of environmental loading on GPS time series and on the reference frame have been reviewed [28]. In addition to the environmental loading effects, GPS stations are also sensitive to other geophysical effects, such as the thermal expansion of their monuments and nearby bedrock. Yan et al. [29] estimated the thermal expansion of 86 globally distributed GPS stations and found that thermal expansion can cause GPS height changes with a few millimeters.

The diversity of climate and geography in China leads to complex environmental loading which deflects the crust. Thus, investigation of the characteristics of the environmental loading effects in China is important to the interpretation of related geophysical phenomena and to the understanding of the nonlinear variations of GPS height time series over this region. Wang et al. [30] reported that ATML and HYDL effects are remarkable in the height time series of 25 investigated GPS stations in China. Jiang et al. [31] investigated 11 GPS stations in China and found that environmental loading corrections could not only reduce annual amplitudes but also change the optimal noise model of the GPS height time series. Recently, Gu et al. [32] evaluated contributions of environmental loading effects on the nonlinear variations of 224 Crustal Movement Observation Network of China (CMONOC) GPS stations. In that paper, the authors reported poor consistency between the vertical annual amplitudes of their GPS and environmental loading time series. They attributed the poor consistency to the uncertainties in the hydrological model (GLDAS Noah025) they used, in which the river flow is ignored. Indeed, the estimations of environmental loading effects are with considerable uncertainty due to the differences of input data, models, and methods used [33]. Currently, GeoForschungsZentrum (GFZ) provides an environmental loading product with the consideration of river flow. Thus, in this paper, we will investigate whether the nonlinear variations of GPS stations in China can be better explained by the GFZ loading products.

Previous studies have found that corrections of environmental loading affect not only the scatter of GPS position time series but also the uncertainty of their secular velocity estimates. Santamaría-Gómez and Mémin [34] reported that environmental loading displacements are characterized by a power-law process which results in an overestimation of GPS velocity uncertainties if loading effects are corrected directly. Recently, Klos et al. [35] examined 376 GPS stations and found that the overestimation of GPS velocity uncertainties can be avoided if environmental loading displacements are modeled with Improved Singular Spectrum Analysis (ISSA). However, the ISSA needs a relatively long data span, such as the threshold of 10 years adopted in that paper. As many GPS stations in the world have only existed for less than 10 years, one purpose of this paper is to find an approach of avoiding the overestimation of velocity uncertainties for the GPS stations with a relatively shorter data span. The efficiency of this approach will be verified with the CMONOC GPS stations.

In this paper, we start by analyzing the spatiotemporal characteristics of GFZ-predicted ATML, NOTL, and HYDL effects on 235 CMONOC GPS stations. We then evaluate the consistency between the environmental loading results and the reprocessed daily CMONOC GPS height time series [36,37,38,39,40,41,42]. Afterwards, we confirm and quantify the loss of velocity precision caused by environmental loading correction. Finally, we try to address the problem of losing velocity precision by employing principal component analysis (PCA) to the GPS height time series with loading correction.

2. Data and Methods

2.1. GPS Height Time Series

We obtained the daily height time series of the 235 CMONOC GPS stations with an observation time span of 2.0–16.5 years from a homogeneous reprocessing of their GPS observations by using GAMIT/GLOBK (Ver. 10.5) software [43]. The strategies of GPS data processing have been described in detail in Jiang et al. [44], in which a set of latest models were adopted. For instance, tropospheric delay was modeled with the state-of-the art Vienna Mapping Function 1 (VMF1) [45] and the a priori zenith hydrostatic delay from European Centre for Medium-Range Weather Forecasts (ECMWF) [46]. Higher-order ionospheric delay was modeled with International Geomagnetic Reference Field 11 (IGRF11) [47] and ionospheric data from the Centre for Orbit Determination in Europe (CODE) [48]. Solid Earth tides, ocean tides, and pole tides were conventionally modeled and corrected [49], but S1 and S2 atmospheric tides were left uncorrected. The locations and time spans of the selected stations are displayed in Figure 1. Among these stations, 27 of them are with time spans longer than 8 years and were termed as CMONOC-I in this paper. The rest of the 208 stations were termed as CMONOC-II.

2.2. Environmental Loading Time Series

Vertical displacement time series due to environmental loading were predicted by the GFZ loading service [25] (http://www.gfz-potsdam.de/en/esmdata/loading/). GFZ provides elastic deformations due to ATML, NTOL, and HYDL on a 0.5° × 0.5° regular global grid with temporal resolutions of 3 h (for ATML, NTOL) and 24 h (for HYDL). The loading time series were referred to the Center of Figure (CF) by using the patched Green’s function approach with load Love numbers from the “ak135” elastic Earth model [50]. The input geophysical data of GFZ-predicted ATML, NTOL, and HYDL were derived from ECMWF operational surface pressure data, Max Planck Institute Ocean Model (MPIOM) [51], and Hydrological Land Surface Discharge Model (LSDM) [52], respectively. Vertical loading time series at a specific station were obtained by bicubic interpolation within the global grid. To be consistent with the time resolution of the GPS time series, the ATML and NTOL time series with 3 h intervals were averaged into daily time series.

2.3. Time Series Analysis and Comparison Methods

To obtain the nonlinear variations of the GPS height time series, the following functional model was used [53]:

$$x\left(t\right)=x_R+v\left(t-t_R\right)+{\displaystyle \sum }_{j=1}^{n_j}{b}_{j}H\left(t-t_j\right)+\epsilon$$

(1)

$$\begin{gathered} H\left(\tau \right)=0 \mathrm{for} \tau <0 \\ H\left(\tau \right)=0.5 \mathrm{for} \tau =0 \\ H\left(\tau \right)=1 \mathrm{for} \tau >0 \end{gathered}$$

(2)

in which $x_R$ and $t_R$ are the reference height and epoch, respectively; $v$ is the secular velocity; $t$ is the epoch; $H\left(\tau \right)$ is the Heaviside step function as defined in Equation (2); $b_j$ and $t_j$ are the offset and its corresponding epoch, respectively; and$\epsilon$ is the residual. Offsets in the GPS height time series caused by earthquake, device change, and unknown reasons were identified by checking earthquake catalogues, station logs, and visual inspection, respectively.

RMS was used to evaluate the power of loading time series and the nonlinear variations of the GPS height time series:

$$\mathrm{RMS}=\sqrt{\frac{{\epsilon }^T\epsilon }{n}}$$

(3)

in which $\epsilon$ is the residuals obtained from Equation (1) and $n$is the number of data points.

To assess the effects of environmental loading corrections on mitigating the nonlinear variations of the GPS height time series, the RMS reduction ratio was introduced:

$${\mathrm{RMS}}_{reduction}=\frac{\mathrm{RMS}\left[{\mathrm{H}}_{\mathrm{before}}\left]-\mathrm{RMS}\right[{\mathrm{H}}_{\mathrm{after}}\right]}{\mathrm{RMS}\left[{\mathrm{H}}_{\mathrm{before}}\right]}$$

(4)

in which ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ are the GPS time series before and after loading correction, respectively, and $\mathrm{RMS}\left[{\mathrm{H}}_{\mathrm{before}}\right]$ and $\mathrm{RMS}\left[{\mathrm{H}}_{\mathrm{after}}\right]$are the RMSs of the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series, respectively.

Additionally, to investigate seasonal signals of the GPS height time series, the following function was used:

$$x\left(t\right)=x_R+v\left(t-t_R\right)+{\displaystyle \sum }_{j=1}^{n_j}{b}_{j}H\left(t-t_j\right)+{\displaystyle \sum }_{k=1}^2{A}_{k}sin\left({\omega }_{k}t+{\phi }_k\right)+\epsilon$$

(5)

in which $A_k$ and ${\phi }_k$ are the amplitude and initial phase of seasonal signals, ${\omega }_k=2\pi /{\tau }_k$, ${\tau }_1=1$ year and ${\tau }_2=1/2$ year. Other quantities are the same as those shown in Equation (1). ${\phi }_k$ was converted into the day of year (doy) when its seasonal signal reached the maximum which was considered as the phase of the seasonal signal in this paper. The seasonal signals of environmental loading time series were also obtained by using Equation (4) but without the estimation of offsets.

Similar to the RMS reduction ratio, the reduction ratio of annual amplitude was also introduced to assess the effects of environmental loading corrections on mitigating the annual variations of the GPS height time series. The reduction ratio of annual amplitude was calculated as follows:

$${\mathrm{A}}_1\_reduction=\frac{{\mathrm{A}}_1\left[{\mathrm{H}}_{\mathrm{before}}\right]-{\mathrm{A}}_1\left[{\mathrm{H}}_{\mathrm{after}}\right]}{{\mathrm{A}}_1\left[{\mathrm{H}}_{\mathrm{before}}\right]}$$

(6)

in which ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ indicate the GPS height time series before and after loading corrections, respectively, and ${\mathrm{A}}_1\left[{\mathrm{H}}_{\mathrm{before}}\right]$and ${\mathrm{A}}_1\left[{\mathrm{H}}_{\mathrm{after}}\right]$are the annual amplitudes of the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series, respectively.

In this study, the Hector software (Ver. 1.6, http://segal.ubi.pt/hector/) was used to estimate the station velocities with maximum likelihood estimation [54]. The noise model was assumed to be power-law plus white noise (PL + WN) in the estimation, which has been widely used in noise analysis of GPS time series by many studies [55,56].

Spatiotemporal correlated components existing in the residual height time series of regional GPS networks are termed as common mode component (CMC). PCA is an approach which has been widely used to filter out CMC, so as to improve the signal-to-noise ratio of GPS height time series [2,57]. The theory and procedures of PCA used to extract the CMC in GPS height time series are described in detail in Dong et al. [57] and generally followed in this paper to calculate the CMC for the residual height series of CMONOC-II.

3. Results

3.1. Quantitative Assessment of Environmental Loading Effects

Figure 2 illustrates the RMS of the vertical displacements due to ATML, NTOL, HYDL, and the summation of their contributions (SumL). The average RMS of ATML, NTOL, HYDL, and SumL for the 235 stations are 3.2, 0.6, 2.7, and 4.0 mm, respectively, indicating that ATML contributes the largest share of the environmental loading at these stations on average. However, the magnitude of the RMS of ATML displacements are very different between these stations and a clear geographic dependence can be noticed. In northwest, northeast, and mid-east China, the RMS of ATML range between 3.0 mm and 5.0 mm. In contrast, in the Qinghai–Tibet Plateau, the RMS of ATML are less than 3.0 mm due to the fact that the atmosphere is much thinner and less active in the Qinghai–Tibet Plateau, which has a mean altitude higher than 4000 m. Moreover, for 34 coastal stations located within 100 km distance of the coast (see Figure 1), most of their RMS of ATML are lower than 3.5 mm, apparently lower than those of their neighboring inland regions due to the inverted barometer (IB) effect. However, stations located close to the west coast of the Bohai Bay still show considerably large RMS of ATML, indicating a more complicated atmosphere–ocean interaction in this region which needs further investigation. Compared with ATML, NTOL has a different spatial pattern characterized by relatively large RMS in coastal regions and very small RMS in inland regions. The RMS of NTOL for the 34 coastal stations range between 0.7 mm and 3.5 mm, significantly larger than those for the 201 inland stations (beyond 100 km distance of the coast), which has an average RMS of 0.4 mm. In particular, the largest RMS of NTOL can be found in the Jiaodong and Liaodong peninsulas. With respect to HYDL, its RMS are less than 3.0 mm in North China and increase gradually to more than 5.0 mm in southwest China, where they are strongly influenced by the Southwest Monsoon. Moreover, the RMS of HYDL can be as large as 6.7 mm in the middle and lower Yangtze River, which are significantly larger than in the neighboring regions. Furthermore, the RMS of SumL range between 2.0 mm to 5.0 mm for most of the stations, aside from stations located in southwest China and in the middle and lower Yangtze River, where they are dominated by large HYDL.

Figure 3 illustrates the vertical displacement time series due to ATML, NTOL, HYDL, and SumL predicted with GFZ loading products and GPS height time series (black) for the selected six CMONOC stations where the RMS of ATML, NTOL, and HYDL reach the maxima or minima. Among the 235 GPS stations that we investigated, stations XJTC (Tacheng, 46.8°N, 82.9°E), XJJJ (Jijitai, 42.8°N, 94.3°E), and XJML (Mulei, 43.8°N, 90.3°E) are characterized by the largest RMS of ATML (4.6 mm), the least RMS of NTOL (0.3 mm), and the least RMS of HYDL (0.5 mm), respectively. These three stations are located in the northeast of the Xinjiang region, which is the hinterland of the Eurasian continent. The northeast Xinjiang is dominated by the strong Mongolia and Siberian high pressure in winter which leads to the largest RMS of ATML in China. Moreover, the RMS of NTOL are minimal in Xinjiang, as it is very far away from the ocean. In addition, Xinjiang is characterized by a typical continental semiarid climate [58], with minimal precipitation for the whole year because of its far distance from the ocean and the prevention of water vapor by the Qinghai–Tibet Plateau and several east–west mountain ranges, such as Kunlun Mountains, Tianshan Mountains, and Altai Mountains. This might explain why HYDL shows the lowest RMS in Xinjiang.

LNJZ (Dalian, 39.0°N, 121.7°E) is a representative station with the largest RMS of NTOL (3.5 mm). This station is located in the Liaodong peninsula, which lies between the Bohai Sea and Yellow Sea. At Dalian, NTOL variations are very strong in winter but relatively weak in summer. Its ATML shows an obvious annual variation with an RMS of 3.0 mm. However, the ATML time series reaches its peak in summer and thus balances out some part of the NTOL variations. HYDL is rather weak at this coastal area with an RMS of 1.6 mm. Consequently, the RMS of SumL for the station LNJZ is only 4.1 mm.

Station YNRL (Ruili, 24.0°N, 97.8°E) is located in Southwestern China and is characterized by the largest RMS of HYDL, with a value of 8.5 mm. In Southwest China, precipitation shows significant annual variations due to the Southwest Monsoon, which transports a great amount of water vapor from the Indian Ocean and the Bay of Bengal [44]. In particular, the annual precipitation for the YNRL station is more than 1400 mm on average during the period of 2011.0–2015.0 [obtained from National Meteorological Information Center of China Meteorological Administration, http://data/cma.cn/]. Compared with HYDL, the ATML and NTOL variations at the YNRL station are much weaker, with RMS of 2.4 and 0.4 mm, respectively. Thus, the SumL time series for this station mainly reflects the contributions from HYDL.

Station YONG (Yongxing, 16.8°N, 112.3°E) is situated at the Yongxing Island of the South China Sea. Due to the strong IB effect at this island, the RMS of ATML is only 0.8 mm, which is the lowest among the 235 stations investigated. The NTOL and HYDL variations are also very small, with RMS of 1.0 and 1.7 mm, respectively. Moreover, the RMS of SumL is 2.0 mm at station YONG, whereas the RMS of its GPS height time series is 10.6 mm. This inconsistency needs further investigation.

3.2. Annual Signals Analysis and Comparison

The left panels of Figure 4 illustrate the vertical annual amplitudes of the ATML, NTOL, HYDL, and SumL time series predicted with GFZ loading products for the CMONOC stations. Comparing Figure 4 with Figure 2, one can find that the geographic distribution patterns of the annual amplitudes of ATML, NTOL, HYDL, and SumL are similar to those of their own RMS. Their mean annual amplitudes are 3.6, 0.4, 3.3, and 4.5 mm, respectively. In addition, the right panels of Figure 4 show the vertical annual phases of the ATML, NTOL, HYDL, and SumL time series for the CMONOC stations. We can see that annual phases of all these quantities also show evident geographic dependences. The annual phases of ATML range between May and June and are earlier in the Qinghai–Tibet Plateau and later in Southwest China. The peaks of the annual signals of NTOL are in January for Northeast China and gradually delayed towards the south, with the latest peaks in July for South China. Annual phases of HYDL are in autumn and winter for Northwest China, whereas they are in early spring for most of other regions. Annual phases of SumL are delayed from southwest to northwest, with the earliest phase in February for Southwest China and latest peaks in July for Northwest China.

Figure 5a displays the vertical annual amplitudes of the ${\mathrm{H}}_{\mathrm{before}}$ time series. The mean annual amplitude of the ${\mathrm{H}}_{\mathrm{before}}$ time series is 4.6 mm, which is very close to the SumL time series predicted by GFZ with a value of 4.5 mm. The geographic distribution patterns of the vertical annual amplitudes of the ${\mathrm{H}}_{\mathrm{before}}$ and SumL time series are also almost the same. To further investigate the consistency between the vertical annual amplitudes of the ${\mathrm{H}}_{\mathrm{before}}$ and SumL time series, a linear fitting line is determined with a weighted total least squares algorithm [59] and shown in Figure 6. The slope of the fitting line is 0.86, which is close to the ideal case of “Slope = 1”.

In addition, Figure 5b displays the vertical annual phases of the ${\mathrm{H}}_{\mathrm{before}}$ time series. The annual phases of the ${\mathrm{H}}_{\mathrm{before}}$ time series are similar to those of the SumL time series in geographic distribution pattern, which is delayed from southwest to northwest. However, it can be also observed that the annual phases of the ${\mathrm{H}}_{\mathrm{before}}$ time series are later than those of the SumL on the whole. To illustrate the phase differences between the vertical annual signals of the SumL and ${\mathrm{H}}_{\mathrm{before}}$ time series, the phase differences are displayed and statistically analyzed in Figure 7. The mean and standard deviations of their annual phase differences are −20.3 and 22.0 days. Moreover, 86% of the 235 stations show negative differences, especially for stations located in Southwest, Northwest, and North China. The phase differences may be caused by mismodeled or unmodeled environmental loading signals in the GFZ products, other unmodeled geophysical signals (e.g., thermal expansions of GPS monuments and nearby bedrock), and GPS errors [1].

The differences between the vertical annual signals of the SumL and ${\mathrm{H}}_{\mathrm{before}}$ time series result in annual signals remained in the ${\mathrm{H}}_{\mathrm{after}}$ time series as displayed in Figure 5c and Figure 8. With environmental loading effects having been taken into consideration, a significant reduction of annual amplitude can be noticed at the majority of the stations. Nevertheless, apparent annual amplitudes are still remaining for stations located in Southwest China and in the middle and lower Yangtze River. These results suggest that the GFZ-predicted environmental loading product is unable to remove the annual signals of GPS time series completely, even though the SumL and ${\mathrm{H}}_{\mathrm{before}}$ time series have a good consistency in annual amplitude. Therefore, the annual signals still need to be estimated when modelling the ${\mathrm{H}}_{\mathrm{after}}$ time series.

3.3. Consistency Assessment with Correlation Coefficients and RMS Reduction

Figure 9 illustrates the correlation coefficients between the SumL and ${\mathrm{H}}_{\mathrm{before}}$ time series for the vertical component of the CMONOC stations. The correlation coefficients generally vary between 0.1 and 0.8, with larger correlation coefficients in North China and the southern Qinghai–Tibet Plateau where the vertical displacement is dominated by ATML and HYDL, respectively. However, the correlation coefficients are relatively poor in Southeast China, where the magnitudes of ATML and HYDL are comparable. The average correlation coefficient is 0.6.

Figure 10 illustrates the RMS reduction ratios for the vertical component of the CMONOC stations with loading corrections from GFZ, with a maximum of about 40% and a mean reduction rate of 20%. Among all the stations, only stations CQCS (Changshou, 29.9°N, 107.2°E) and GDZH (Zhuhai, 22.3°N, 113.6°E) show increased RMSs after the loading correction. CQCS is close to the Yangtze River and is about 400 km upstream from the Three Gorges Dam. GDZH is a coastal station located in the Pearl River Delta. The increase of their RMSs after the loading correction suggest a poor performance of the GFZ loading products at these two stations which needs further investigation. In Southwest China, the RMS reduction ratios generally range between 8% and 24%, although the annual amplitudes predicted by GFZ loading products and observed by GPS are comparable. The poor RMS reductions for these stations might be related to the remaining annual signals in the ${\mathrm{H}}_{\mathrm{after}}$ time series as shown in Figure 5c. Further investigations are still needed to find out the reason.

3.4. Stacked Power Spectra Analysis

We analyzed the spectrum characteristics of the ATML, NOTL, HYDL, and SumL time series by using stacked power spectra [3] to investigate the influences of environmental loading corrections on the spectra property and velocity uncertainties of relatively short GPS time series (<5 years). Only the selected 208 CMONOC-II stations were investigated as their GPS observation data spans are relatively short, with values of ~4.6 years.

We firstly obtained the power spectra of each loading time series by using the Lomb–Scargle periodogram technique [60] with an oversampling factor of 4 [61]. Then, we stacked the power spectra of the ATML, NOTL, HYDL, and SumL time series, respectively. The stacked power spectra of each of the abovementioned loading quantities were obtained and are displayed in Figure 11. Regarding the stacked power spectra of the ATML, a sudden drop at frequencies higher than 40 cpy can be seen, indicating an autoregressive behavior. This result confirms the findings of earlier studies [19,26,35]. Likewise, the frequency characteristics of NTOL’s stacked power spectra also follow an autoregressive process. However, the stacked power spectra of the HYDL illustrate a stochastic behavior close to random walk noise. Furthermore, the stacked power spectra of the SumL are characterized by autoregressive processes, which are similar to those of ATML and NTOL.

In the next step, we compared the residuals of the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series fitted with Equation (5). Their stacked power spectra are displayed in Figure 12. The stacked power spectrum of the ${\mathrm{H}}_{\mathrm{before}}$ time series is characterized by a power-law noise with the slope close to −1 (flicker noise) at low frequencies and close to 0 (white noise) at high frequencies. With the environmental loading signals subtracted, the stacked power spectrum of the ${\mathrm{H}}_{\mathrm{after}}$ time series is closer to horizontal at high frequencies and an evident power reduction can be noticed at the frequency band between 10 cpy and 50 cpy. A similar result is also found by Klos et al. [35], in which the phenomenon of power loss is considered as an artifact. We also notice that the phenomenon of power loss changes the noise property of the GPS height time series, which could influence the estimates of velocity uncertainty.

Table 1. Mean and standard deviations of the spectral index (κ) and velocity uncertainty (${\mathsf{\sigma }}_{vel.}$) of the ${\mathrm{H}}_{\mathrm{before}}$, ${\mathrm{H}}_{\mathrm{after}}$, CMC-filtered ${\mathrm{H}}_{\mathrm{before}}$, and CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series for the selected 208 CMONOC-II stations estimated with PL + WN noise model. VD and DP indicate the velocity difference and dilution of precision with respect to the ${\mathrm{H}}_{\mathrm{before}}$ time series.

	${\mathbf{H}}_{\mathbf{before}}$	${\mathbf{H}}_{\mathbf{after}}$	CMC-filtered ${\mathbf{H}}_{\mathbf{before}}$	CMC-filtered ${\mathbf{H}}_{\mathbf{after}}$
κ	−0.8 ± 0.1	−1.2 ± 0.2	−0.7 ± 0.3	−0.8 ± 0.3
${\mathsf{\sigma }}_{vel.}$ (mm/year)	0.9 ± 0.3	1.2 ± 0.5	0.4 ± 0.3	0.4 ± 0.4
VD (mm/year)	-	0.1 ± 0.4	−0.1 ± 0.3	0.1 ± 0.4
DP	-	1.4 ± 0.5	0.5 ± 0.3	0.4 ± 0.3

lists the mean and standard deviations of the spectral index and velocity uncertainty for the CMONOC-II ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series estimated with PL + WN noise model. We can see that after the loading correction, the average spectral index drops from −0.8 to −1.2. Moreover, the average velocity uncertainty is 0.9 and 1.2 mm/year for the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series, respectively, indicating a significant overestimation of the velocity uncertainty caused by the loading correction. A detailed discussion will be given in the next section.

4. Discussion

4.1. The Effects of GFZ Loading Products on the Nonlinear Variations

In this section, we discuss the effects of GFZ loading products on the nonlinear variations of the CMONOC stations with a comparison to a previous work, Gu et al. [32], which used other loading products. As shown in Figure 6, the slope of the fitting line is 0.86, which is close to the ideal case of “Slope = 1” and significantly superior to the slope of 0.4 obtained from Gu et al. [32]. In that paper, the authors attributed the poor consistency between the vertical annual amplitudes of their environmental loading and GPS time series to the uncertainties in the hydrological model (GLDAS Noah025) they used, especially for Southwest China, where the vertical annual amplitudes of environmental loading signals they derived are significantly smaller than those of GPS time series. Southwest China is located within the drainage basins of the Yangtze (Jinsha), Mekong (Lancang), and Salween (Nujiang), with a monsoon climate, and is thus very sensitive to the HYDL variations due to the river flow. The influence of river flow has been taken into consideration in the GFZ-predicted HYDL but ignored in GLDAS Noah025 model. Thus, the HYDL and SumL signals may be better predicted by products from GFZ than those used in that paper. However, the average correlation coefficients between the SumL and ${\mathrm{H}}_{\mathrm{before}}$ time series for the vertical component of the CMONOC stations is 0.6, worse than the value of 0.76 obtained from Gu et al. [32]. This relatively poor correlation coefficient might be due to the significant annual phase differences illustrated in Figure 7, as the correlation coefficient between two signals is sensitive to their phase difference.

4.2. The Impacts of Environmental Loading Corrections on Velocity Uncertainties

In this section, we discuss the impacts of environmental loading corrections on the velocity uncertainties of the CMONOC-II stations. To quantify the changes of velocity uncertainties, we adopted the dilution of precision (DP) as an index [5,35], which is defined as:

$$\mathrm{DP}=\frac{{\sigma }_{vel.}\left[\mathrm{H}\right]}{{\sigma }_{vel.}\left[{\mathrm{H}}_{\mathrm{before}}\right]}$$

(7)

in which ${\mathsf{\sigma }}_{vel.}\left[\mathrm{H}\right]$ indicates the velocity uncertainty of one solution and ${\mathsf{\sigma }}_{vel.}\left[{\mathrm{H}}_{\mathrm{before}}\right]$ indicate the velocity uncertainty of the ${\mathrm{H}}_{\mathrm{before}}$ solution as a reference. If the DP of one solution is larger than 1, it means the velocity uncertainly is overestimated with respect to that of the ${\mathrm{H}}_{\mathrm{before}}$ solution.

Table 1 lists the mean and standard deviations of the spectral index and velocity uncertainty for the CMONOC-II ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series estimated with PL + WN noise model. We can see that after the loading correction, the average DP is 1.4, indicating a significant overestimation of velocity uncertainty caused by the loading correction. This overestimation of velocity uncertainty might arise from the power loss and noise property changes of the GPS time series after loading correction as shown in Figure 12.

Figure 13b illustrates the DP of the ${\mathrm{H}}_{\mathrm{after}}$ solution. Among the 208 CMONOC-II stations investigated, 165 of them are characterized with DP values larger than 1.0, which indicates that the loading correction leads to an increase of velocity uncertainty for 79% of the stations. Moreover, 28 out of the 208 stations are characterized by DP values larger than 2.0. The DP values of station GSLX (Longxi, 35.0°N, 104.6°E) and QHGC (Gangcha, 37.3°N, 100.1°E) are even larger than 3.0. Whereas for stations located in Southwest China dominated by HYDL, their DP values are generally less than the average value of 1.4. Therefore, we can infer that HYDL is not the major cause of the overestimation of velocity uncertainty. As for North and Northwest China, the DP values are generally larger than 1.4. ATML signals are stronger than the NTOL and HYDL signals in these regions. Thus, the overestimation of velocity uncertainty may be caused by the correction of ATML effects [35]. However, this argument does not apply to Northeast China, which is also dominated by ATML rather than HYDL, but is characterized with DP values generally less than 1.0. Therefore, we cannot simply judge the overestimation of velocity uncertainty by the domination of ATML. The overestimation of the velocity uncertainty should be estimated and confirmed case by case.

4.3. Impact of CMC Filtering on the Velocity Uncertainty

As known from the existing literature, a kind of spatiotemporal correlated component, termed as CMC, exists in the residual height time series of regional GPS networks. The CMC might arise from mismodeled and unmodeled geophysical effects, GPS errors, and reference frame definition [62]. Additionally, the ${\mathrm{H}}_{\mathrm{after}}$ time series might also suffer from the CMC since the geophysical effects are corrected imperfectly, and the errors due to GPS data processing methods, error models, and reference frame definition are not accurate enough. Therefore, we tried to solve the problem of overestimation of velocity uncertainty caused by the environmental loading correction by a CMC filter with PCA approach.

We first confirmed the CMC in the CMONOC network by checking the interstation correlations as shown in Figure 14. The mean interstation correlation coefficient of the ${\mathrm{H}}_{\mathrm{after}}$ time series is 0.29, which is slightly lower than that of the ${\mathrm{H}}_{\mathrm{before}}$ time series with a value of 0.33. These results confirm the existence of CMC in both the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series and indicates that environmental loading correction can only reduce the interstation correlation coefficient for the CMONOC stations moderately. We now discuss two questions. One is whether the CMC filter for the ${\mathrm{H}}_{\mathrm{after}}$ time series with PCA approach can compensate the precision loss of the velocity estimates that arises from loading correction. The other is whether a loading correction is necessary for the time series with their CMC to be filtered.

CMC filters for the ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series are performed with the PCA approach that has been described in the Section 2. Figure 14 shows that both the mean interstation correlation coefficients of the CMC-filtered ${\mathrm{H}}_{\mathrm{before}}$ and CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series are reduced to zeroes, indicating that PCA is an effective CMC filter approach for both ${\mathrm{H}}_{\mathrm{before}}$ and ${\mathrm{H}}_{\mathrm{after}}$ time series. From Table 1 we can find that the mean spectral index of the ${\mathrm{H}}_{\mathrm{after}}$ time series turns from −1.2 to −0.8 after applying the CMC filter, which is equal to that of the ${\mathrm{H}}_{\mathrm{before}}$ time series. Moreover, the mean velocity uncertainty of the ${\mathrm{H}}_{\mathrm{before}}$, ${\mathrm{H}}_{\mathrm{after}}$, and CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series are 0.9, 1.2, and 0.4 mm/year, respectively. Figure 13f shows the DP of the CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ solution with respect to the ${\mathrm{H}}_{\mathrm{before}}$ time series. The DP values of the CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series are larger than 1.0 for only 5 stations and less than 1.0 for the other 203 stations, and with a mean value of 0.4, whereas the mean DP of the ${\mathrm{H}}_{\mathrm{after}}$ time series with respect to the ${\mathrm{H}}_{\mathrm{before}}$ solution is 1.4. These results indicate that the CMC filter is capable of not only compensating the precision loss of velocity estimates due to the loading correction but also improving the precision of velocity estimates.

Table 1 also shows that the average velocity uncertainties of the CMC-filtered ${\mathrm{H}}_{\mathrm{before}}$ and CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series are the same, with a value of 0.4 mm/year. Nevertheless, Figure 13c,e show that the numbers of stations with absolute values of their velocity differences larger than 0.5 mm/year are 19 and 32 for the CMC-filtered ${\mathrm{H}}_{\mathrm{before}}$ and CMC-filtered ${\mathrm{H}}_{\mathrm{after}}$ time series, respectively. These results indicate that, compared with the solely CMC-filtered time series, the time series that applies loading correction firstly and CMC filtering secondly obtains no further improvement in velocity uncertainty but a larger opportunity to yield extreme velocity differences with respect to the ${\mathrm{H}}_{\mathrm{before}}$ time series.

5. Conclusions

Environmental loading effects contribute much to the nonlinear variations of GPS height time series. Thus, an accurate model and removal of the environmental loading effects is important to the improvement of the signal-to-noise ratio of the GPS time series. By using GFZ-predicted loading data and daily height time series of 235 GPS stations derived from a homogeneously reprocessed CMONOC network, this paper assessed the effects of environmental loading correction on the nonlinear variations of the GPS height time series. Moreover, the problem of overestimation of GPS velocity uncertainties caused by environmental loading correction is confirmed and solved by the proposed approach of CMC filter. This approach breaks through the restriction of the ISSA approach on the requirement of the length of time span (>10 years) and has been proven to be efficient for GPS time series with relatively short time span (~4.6 years).

Results show that the average RMS of vertical crustal displacements due to the integrated environmental loading effects is 4.0 mm, in which nontidal atmospheric loading comprises the most significant effect with an average RMS of 3.2 mm. Hydrological loading also plays an important role but its effect is unbalanced in China, with an RMS ranging between 0.5 mm and 8.5 mm. Effects of nontidal oceanic loading is only significant for coastal and island stations, especially for the stations located in the Jiaodong and Liaodong peninsulas, with RMS reaching up to 3.5 mm. Regarding the annual signals, the GFZ-predicted integrated vertical environmental loading signals fit well with the GPS height time series in annual amplitude, with a slope of 0.86 for the fitting line. However, a systematic difference can be noticed between their annual phases. With respect to the consistency between the loading and GPS height time series, their average correlation coefficient and RMS reduction ratio is 0.6 and 20%, respectively.

In addition, spectral analysis results show that hydrological loading time series are dominated by random walk noise, while atmospheric and nontidal oceanic loading time series are characterized by an autoregressive behavior. Therefore, environmental loading corrections would change the statistical properties of GPS height time series and the derived velocity uncertainty. An investigation of 208 CMONOC-II stations with observing time spans of ~4.6 years shows that the velocity uncertainty for 79% of the stations are enlarged after the loading correction. Then, PCA is applied to mitigate the CMC in the residual time series. Results show that the average dilution of velocity precision of the loading corrected time series is reduced from 1.4 to 0.4 with PCA applied, indicating that PCA is an effective approach to compensate the dilution of velocity precision due to loading correction. However, compared with the solely CMC-filtered time series, a time series with loading correction firstly and CMC filter secondly obtains no further improvement in velocity uncertainty but at a larger opportunity to yield extreme velocity differences with respect to the time series with neither CME filter nor environmental loading correction. These results indicate that environmental loading correction is not necessary for the CMONOC GPS height time series if their CMC will be filtered out.