Impacts of Local Effects and Surface Loads on the Common Mode Error Filtering in Continuous GPS Measurements in the Northwest of Yunnan Province, China

While seasonal hydrological mass loading, derived from Gravity Recovery and Climate Experiment (GRACE) measurements, shows coherent spatial patterns and is an important source for the common mode error (CME) in continuous global positioning system (cGPS) measurements in Yunnan, it is a challenge to quantify local effects and detailed changes in daily GPS measurements by using GRACE data due to its low time and spatial resolutions. In this study, we computed and compared two groups of CMEs for nine cGPS sites in the northwest Yunnan province; rCMEs were computed with the residual cGPS time series having high inter-station correlations, while oCMEs were computed with all the GPS time series. The rCMEs-filtered time series had smaller variances and larger root mean square (RMS) reductions than those that were oCMEs-filtered, and when the stations local effects were not removed, spurious transient-like signals occurred. Compared with hydrological mass loading (HYDL), its combination with non-tidal atmosphere pressure and ocean mass reached a better agreement with the CME in the vertical component, with the Nash–Sutcliffe efficiency (NSE) increasing from 0.28 to 0.55 and the RMS reduction increasing from 15.19% to 33.4%, respectively. Our results suggest that it is necessary to evaluate the inter-station correlation and remove the possible noisy stations before conducting CME filtering, and that one should carefully choose surface loading models to correct the raw cGPS time series if CME filtering is not conducted.


Introduction
Continuous GPS (cGPS) measurements not only provide precision position time series of tectonic signals (e.g., plate motion, crustal strain accumulation) [1][2][3][4] but record temporal responses to the variations of surface loading from various sources such as tidal loading, atmospheric pressure, continental water storage [5][6][7][8][9][10][11][12], and unmodeled systematic errors including the mismodelling of satellite orbits, receiver antenna phase center corrections, and multipath effects [6,13]. Both environmental responses and systematic errors will result in spatiotemporal coherent features of cGPS measurements in a rather large area that may extend to more than 2000 km [14]. If the magnitudes of such coherent features are large enough or larger than the tectonic responses, it will be a challenge to detect useful information due to the high uncertainties. As a result, such coherent features are usually called common mode errors (CMEs), which are the main spatially-correlated signals in the cGPS measurements. When CME is extracted and removed from the measurements, the remaining more evidence has suggested that the annual signals in cGPS time series are a response to variations in the Earth's surface mass loading [6,24,29], including changes of non-tidal atmospheric pressure, ocean loading, and hydrological water storage. Furthermore, in these previous studies, the time series of the vertical component were occupied to obtain the CME without ignoring any weakly correlated time series. Taking XIAG (Figure 1a) as an example, XIAG station is located on the south bank of Erhai Lake, which is the second largest reservoir in Yunnan province; the water level fluctuates annually with an amplitude about 3 m and has resulted in significant anomalies in gravity measurements [30]. Such loading effects associated with large local hydrological water storage may further result in significant responses in cGPS measurements [27]. Nevertheless, both the monthly time interval and 400 km spatial resolution of GRACE data are too low to discern these local effects or the detailed changes in daily GPS measurements [24,27,28]. Therefore, it is necessary to quantitatively evaluate the impacts of local effects and loading deformation on the CME filtering for daily cGPS measurements.

Loading Deformation Due to Atmosphere, Non-Tidal Ocean, and Hydrological Water Mass
For each GPS station, we obtained the daily elastic displacement time series resulting from surface loads of non-tidal atmospheric and oceanic pressure and hydrological mass variations with gridded loading displacements stored on a regular 0.5° × 0.5° global grid with 24 h hydrological water storage (HYDL) and 3 h non-tidal atmospheric pressure (NTAL) and ocean loading (NTOL) provided by the GeoForschungsZentrum (GFZ) loading service (https://isdc.gfz-potsdam.de/esmdata/). The elastic surface deformations are calculated in the spatial domain by convolving the loading of Green's function in the center of the Earth's figure frame (CF) with the loaded Love numbers using modelled mass distributions from the models ECMWF (European Center for Medium-Range Weather Forecasts), MPIOM (Max-Planck-Institute for Meteorology Ocean Model), and Land Surface Discharge Model [32]. Please refer to Dill and Dobslaw [32] and the Geo Forschungs Zentrum (GFZ) loading service (https://isdc.gfz-potsdam.de/esmdata/) for more detailed information about the loading deformation datasets and mass distribution models.

Evaluation of Similarity among GPS Time Series
The Pearson's correlation coefficients (R) is used to measure the statistical relationship between two measurements. For this purpose, as a case study, we took nine cGPS stations in the northwest of Yunnan Province into account to stress the importance of considering the local effects before CME. We first evaluated the correlations among the position time series of these cGPS measurements. With the highly correlated seven stations as a supervised group, we obtained the CMEs time series by weighted stacking for three components of these stations, and then the supervised results were compared with the unsupervised one. Finally, we compared the CMEs with the modeled daily loading displacements, induced by non-tidal atmospheric and oceanic pressures, and surface hydrological mass variations. Our results suggest that it is necessary to evaluate the inter-station correlation and remove the possible noisy Sensors 2020, 20, 5408 4 of 15 stations before conducting CME filtering, and that one should carefully choose surface loading models to correct the raw cGPS time series instead of using CME filtering.

Processing of GPS Datasets and Post-Processing of the GPS Residual Daily Time Series
Daily vertical displacements of nine cGPS stations (Figure 1) from the Crustal Movement Observation Network of China (CMONOC) were processed with GAMIT/GLOBK software [31] in an ITRF2008 reference frame [2]. In the processing, we used the VMF1 tropospheric mapping function to estimate the hydrostatic and wet zenith delays. The effects from polar motion, solid Earth tides, and ocean tides were removed in the data processing. The XIAG site has been continuously observed since the late 1990s, while the other stations have been continuously measured since mid-2010 and/or mid-2011. The locations of all the stations are shown in Figure 1 and Table 1. In the tectonic background, the long-term rates of three components of a continuous GPS station remain constant. For simplification, a linear rate v 0 and initial position x 0 are estimated with least squares regression, while the seasonal terms are retained as we regard the common periodicity as a CME. The residuals' time series of each component r(t) can be obtained by removing the linear trend from the measurements.
where t is the time and t 0 the origin of time, x 0 is the initial coordinate at time t = t 0 , v 0 is a constant, the linear velocity of the point, and ε is the noise term. The noise term is observational white noise ε ∼ N 0, σ 2 , where σ 2 is the error variance of daily solution.
We applied the median-interquartile range (IQR) algorithm [16] to detect and remove the outliers in the postfit residuals. In this cleaning process, the IQR of a dataset is defined as the difference between its 75th and 25th percentiles, the values are outliers when the absolute values of difference between the dataset and its median are larger than three times the IQR. The postfit residual position time series are cleaned separately in the east, north, and vertical (up) directions.

Loading Deformation Due to Atmosphere, Non-Tidal Ocean, and Hydrological Water Mass
For each GPS station, we obtained the daily elastic displacement time series resulting from surface loads of non-tidal atmospheric and oceanic pressure and hydrological mass variations with gridded loading displacements stored on a regular 0.5 • × 0.5 • global grid with 24 h hydrological water storage (HYDL) and 3 h non-tidal atmospheric pressure (NTAL) and ocean loading (NTOL) provided by the GeoForschungsZentrum (GFZ) loading service (https://isdc.gfz-potsdam.de/esmdata/). The elastic surface deformations are calculated in the spatial domain by convolving the loading of Green's function in the center of the Earth's figure frame (CF) with the loaded Love numbers using modelled mass distributions from the models ECMWF (European Center for Medium-Range  [32]. Please refer to Dill and Dobslaw [32] and the Geo Forschungs Zentrum (GFZ) loading service (https://isdc.gfz-potsdam.de/esmdata/) for more detailed information about the loading deformation datasets and mass distribution models.

Evaluation of Similarity among GPS Time Series
The Pearson's correlation coefficients (R) is used to measure the statistical relationship between two measurements.
where y p (i) and y q (i) are the values of p-th and q-th measurements at the i-th epoch, with mean values of y p and y q , respectively. N is the time length of observations. The R values of 1 and −1 denote a perfect positive and negative correlation between two measurements, respectively. Because the correlation coefficient value provides no information about the magnitude discrepancy between the two datasets, we further computed the Nash-Sutcliffe Efficiency (NSE) [33] to evaluate the consistency between CMEs and loading deformation. The NSE is a function to quantitatively evaluate the amplitude discrepancy, which indicates how well the observed data fits the simulated data. NSE ranges from infinity to 1. When NSE = 1, it means that these two variables have perfect consistency; when NSE = 0, it indicates that simulated data are equivalent with the mean of the observation data. Whereas if the NSE value is less than 0, the discrepancy of amplitudes between two variables becomes unacceptable.
where y obs (i) and y pred (i) are the values of the i-th observation and prediction, respectively; y obs and y pred are the mean values of observation and prediction, respectively, and N is the number of stations.

Estimate of the CMEs in cGPS Position Time Series
When the correlation coefficients are larger than 0.5, the postfit residuals are then stacked to compute the CME of each component by weighting with the inverse of the error variance of the daily solutions.
where r i,j and σ 2 i,j are the values of residuals and error variance of daily solutions of i-th station at j-th epoch, respectively.
In order to test whether the low correlated station will affect the CMEs, we set two groups of stations to compute the CMEs for all three components. In the first group, all nine stations were involved in computation of CMEs. While in the second group, only the seven stations that showed high inter-correlation were selected. Hereinafter, the CMEs derived from the first and the second groups are referred to as oCMEs and rCMEs, respectively.

Evaluation of CMEs
For the purpose of evaluating the reliability of these two CME corrections, we calculated the variance of the corrected time series for the east, north, and up components. The lower value of variance of the corrected time series indicates a better correction. The F-test was applied to evaluate which variance was smaller for the postfit residuals after correction of oCMEs and rCMEs. First, we set the null hypothesis H0 = 0 by assuming these two variances are statistically equal to each other. The F-statistic is the ratio of two variances, i.e., F = σ 2 rCME /σ 2 oCME , where σ 2 rCME and σ 2 oCME are variances for the residual time series after rCME and oCME corrections. If the F-statistic is larger than the upper critical value or less than the lower critical value, the null hypothesis is false; and meanwhile if the p-value is also less than the significant level, e.g., 5%, then we will reject the null hypothesis. In this case, if the F-statistic is less than the upper critical value, the variance of the numerator is less than that of the denominator; whereas if the F-statistic is larger than the lower critical value, the variance of the numerator is larger than that of denominator. Accordingly, we can make a decision to select the better result to correct the original GPS time series.
Furthermore, the RMS reduction is also used to evaluate how well the CME fits the residual position time series.
where RMS GPS and RMS GPS−CME are the RMSs of the residual position time series and the CMEs for three components, respectively. RMS reduction can also be used to test what percentage the loading deformation contributes to the CME, in which case the subscript GPS and CME can be replaced with CME and Loading, respectively.

Residual cGPS Time Series
The residual time series of three components, after removal of a best fitting linear trend, are shown in Figure 1 for all nine cGPS stations. For the up component, the residual time series show strong annual and semi-annual variations at all stations. Table 1 shows the amplitudes of annual and semi-annual variations of the three components. It is notable that, for XIAG and YNJD, the residual time series of the east and north components are much more scattered than those of all the other stations. Furthermore, there are abnormal transient-like signals in both the east and north components of XIAG during the period from 2015 to 2016. As such, these two stations were suspected to be significantly affected by local effects, and the reliability in the CME corrections may be also vulnerable. We further computed two types of CMEs by evaluating the inter-station correlation and show the difference between these two CME corrections.  Table 1 show the inter-station correlation coefficients of three (east, north, and up) components for 36 GPS station pairs. The correlation coefficients are relatively low for the north and east components. For the east component, only five station pairs show correlation coefficients larger than 0.5. For the north component, the correlation coefficients are greater than 0.5 for 19 station pairs. In comparison, the up component shows larger correlation coefficients, with 34 station pairs showing correlation coefficients larger than 0.5, and there are more than 26 station pairs with correlation coefficients larger than 0.65. Figure 3 shows that the correlation coefficients decrease with the increasing inter-station distances. Although Figure 3a shows an overall negative linear trend, the variance is too large to be expressed as a linear trend. While after ignoring the low coefficients, both north and up components show linear trends with correlation coefficients decreasing with increasing distances. This confirms that the correlation coefficient becomes smaller when the sub-network is extended over larger areas. coefficients larger than 0.65. Figure 3 shows that the correlation coefficients decrease with the increasing inter-station distances. Although Figure 3a shows an overall negative linear trend, the variance is too large to be expressed as a linear trend. While after ignoring the low coefficients, both north and up components show linear trends with correlation coefficients decreasing with increasing distances. This confirms that the correlation coefficient becomes smaller when the sub-network is extended over larger areas.   Table 1.

Comparison between the oCMEs and rCMEs
Figure 4a-c shows the oCMEs and rCMEs for three components, respectively. In order to display and distinguish between oCMEs and rCMEs more easily, we add an offset to each component; the offsets were ±3.5, ±4.5, and±10.5 for east, north, and up components, respectively; also shown are the modeled annual and semi-annual variations in blue curves, which show that the vertical CMEs have the most remarkable seasonal variations, the north CMEs the second highest variation, while the east CMEs have minor seasonality. Such discrepancies further reflect the different inter-station correlations for the three-component displacement time series. Because of the similarities in both amplitude and phase for each component of oCMEs and rCMEs, it is not easy to distinguish between them by visual inspection. We computed the variances of postfit residuals after rCMEs and oCMEs filtering, i.e., (rCME) and (oCME) after removing a combination of linear trend and annual and semi-annual variations. Our calculations show that the east and up components are significantly

Comparison between the oCMEs and rCMEs
Figure 4a-c shows the oCMEs and rCMEs for three components, respectively. In order to display and distinguish between oCMEs and rCMEs more easily, we add an offset to each component; the offsets were ±3.5, ±4.5, and ±10.5 for east, north, and up components, respectively; also shown are the modeled annual and semi-annual variations in blue curves, which show that the vertical CMEs have the most Sensors 2020, 20, 5408 8 of 15 remarkable seasonal variations, the north CMEs the second highest variation, while the east CMEs have minor seasonality. Such discrepancies further reflect the different inter-station correlations for the three-component displacement time series. Because of the similarities in both amplitude and phase for each component of oCMEs and rCMEs, it is not easy to distinguish between them by visual inspection. We computed the variances of postfit residuals after rCMEs and oCMEs filtering, i.e., σ 2 (rCME) and σ 2 (oCME) after removing a combination of linear trend and annual and semi-annual variations. Our calculations show that the east and up components are significantly different based on the F-test and p-values at the significant level of 5% (Figure 4a,c). Nevertheless, the postfit variances of rCME and oCME for the north component are statistically the same based on both F-test and p-value, with a value larger than 0.05 (Figure 4b).

CME Corrections of GPS Position Time Series
We corrected the GPS time series for three components with oCMEs and rCMEs. We took YNYA as an example to show the differences between oCME correction (yellow circles) and rCME correction (gray circles) for east, north, and up components ( Figure 5). On the whole, the oCME-corrected time series are more scattered than those of the rCME correction. More specifically, after rCMEs correction, the variances of filtered time series of east, north, and up components decreased from 3.9 mm 2 , 2.9 mm 2 , and 36.1 mm 2 to 2.2 mm 2 , 1.3 mm 2 , and 17.1 mm 2 , respectively, which are statistically smaller than those with oCME correction (Table 2). Table 2 and Figure 6 show the variances of postfit residuals of east, north, and up components for all nine cGPS stations. In Figure 6, the bars in blue, green, and red denote variances of the postfit residuals for the original time series and the corrected time series with oCME and rCME, respectively. The F-test was also applied to evaluate whether the variance after rCMEs correction was less than that of the oCME correction. Except for XIAG and YNYL, nearly all the other seven stations showed obvious reduction of variances of postfit residuals for all three components when rCMEs were used instead of

CME Corrections of GPS Position Time Series
We corrected the GPS time series for three components with oCMEs and rCMEs. We took YNYA as an example to show the differences between oCME correction (yellow circles) and rCME correction (gray circles) for east, north, and up components ( Figure 5). On the whole, the oCME-corrected time series are more scattered than those of the rCME correction. More specifically, after rCMEs correction, the variances of filtered time series of east, north, and up components decreased from 3.9 mm 2 , 2.9 mm 2 , and 36.1 mm 2 to 2.2 mm 2 , 1.3 mm 2 , and 17.1 mm 2 , respectively, which are statistically smaller than those with oCME correction (Table 2).  The bars in blue, green, and red denote variances of the postfit residuals for the original time series and the time series after oCME and rCME filtering, respectively. Except for XIAG and YNYL, most of the variances after rCME correction are less than those corrected with oCME for all three components. The numbers labeling the horizontal axis denote the stations listed in Table 1. Comparison between rCME correction (gray circles) and oCME correction (yellow circles) for (a) up, (b) north, and (c) east components of YNYA; also shown are the weighted residuals (blue circles) of three components of XIAG, which are weakly correlated with all the other stations. The curves in blue are non-linear trends. The R oCME and R rCME are correlation coefficients between residual time series of XIAG and that after oCME and rCME corrections of YNYA. The dashed red curves in oCME-corrected north and east components denote the spurious transient-like signals. σ 2 dGPS denotes variance of detrended GPS time series. The superscripts with #0 denote that variances are statistically equal to each other after rCME and oCME corrections, while those with #1 denote that variances with oCME correction are less than those with rCME correction. For the remainder, the variances with rCME correction are less than those with oCME correction. Table 2 and Figure 6 show the variances of postfit residuals of east, north, and up components for all nine cGPS stations. In Figure 6, the bars in blue, green, and red denote variances of the postfit residuals for the original time series and the corrected time series with oCME and rCME, respectively.
The F-test was also applied to evaluate whether the variance after rCMEs correction was less than that of the oCME correction. Except for XIAG and YNYL, nearly all the other seven stations showed obvious reduction of variances of postfit residuals for all three components when rCMEs were used instead of oCMEs. Four exceptions of statistically equal variances were for the east and up components of YNYL, the north component of YNZD, and the up component of YNJD (Table 2, Figure 6).

Figure 5.
Comparison between rCME correction (gray circles) and oCME correction (yellow circles) for (a) up, (b) north, and (c) east components of YNYA; also shown are the weighted residuals (blue circles) of three components of XIAG, which are weakly correlated with all the other stations. The curves in blue are non-linear trends. The R and R are correlation coefficients between residual time series of XIAG and that after oCME and rCME corrections of YNYA. The dashed red curves in oCME-corrected north and east components denote the spurious transient-like signals. Figure 6. The variances of postfit residuals for (a) east, (b) north, and (c) up components. The bars in blue, green, and red denote variances of the postfit residuals for the original time series and the time series after oCME and rCME filtering, respectively. Except for XIAG and YNYL, most of the variances after rCME correction are less than those corrected with oCME for all three components. The numbers labeling the horizontal axis denote the stations listed in Table 1. The bars in blue, green, and red denote variances of the postfit residuals for the original time series and the time series after oCME and rCME filtering, respectively. Except for XIAG and YNYL, most of the variances after rCME correction are less than those corrected with oCME for all three components. The numbers labeling the horizontal axis denote the stations listed in Table 1.
Furthermore, the RMS reductions of oCMEs and rCMEs for all three components are shown in Table 3. It shows that when rCMEs are used instead of oCMEs, the RMS reductions are increased for the three components of almost all the highly inter-correlated seven stations except for YNZD, whose RMS reduction of oCME is slightly larger than that of rCME in the north component. Table 3. Root mean square (RMS) reductions of CME corrections for GPS time series.

Station Name
East RMS Reduction (%) North RMS Reduction (%) Up RMS Reduction (%) rCME oCME rCME oCME rCME oCME  Figure 7 shows rCMEs (gray open circles) stacked from the highly correlated residual position time series; the annual amplitudes are 0.45 mm, 1.05 mm, and 8.70 mm for the east, north, and up components, respectively. The stacked displacements are loaded from HYDL, NTAL, and NTOL. In these three loading terms, the HYDL loadings have the largest annual amplitudes of 0.57 mm, 1.15 mm, and 7.28 mm for the east, north, and up components, respectively. In comparison, the effects due to the non-tidal atmosphere and ocean are minor, with the annual amplitudes of 2.61 mm and 0.23 mm in vertical displacements, respectively, while the annual amplitudes are less than 0.1 mm for both east and north components.  dGPS denotes variance of detrended GPS time series. The superscripts with #0 denote that variances are statistically equal to each other after rCME and oCME corrections, while those with #1 denote that variances with oCME correction are less than those with rCME correction. For the remainder, the variances with rCME correction are less than those with oCME correction. Furthermore, the RMS reductions of oCMEs and rCMEs for all three components are shown in Table 3. It shows that when rCMEs are used instead of oCMEs, the RMS reductions are increased for the three components of almost all the highly inter-correlated seven stations except for YNZD, whose RMS reduction of oCME is slightly larger than that of rCME in the north component. For each loading term, the NSE with the rCMEs and the corresponding RMS reductions are listed in Table 4. It is notable that there are apparent phase lags between the loading responses and the rCMEs of the three components. As a result, neither the RMS reduction nor the NSE show significant agreement between the loading effect and rCMEs. Although the HYDL loading can contribute a 15.19% RMS reduction to the CME of the up component, it is relatively low compared with the consistency of their annual amplitudes. We thus compared the rCMEs and combinations of these three loading terms. Figure 8 shows both rCMEs (gray open circles) and the combined loading deformation (green open circles) and the modeled annual and semi-annual variations for three components. As expected, the combined loading deformation in the up component has the most remarkable seasonal variations and agrees much better with the CME than before combination; the NSE increases from 0.28 to 0.55, and RMS reduction increases from 15.19% to 33.4% (Table 4). In contrast, the phase lag in the east component is too large to reach a good agreement between them, resulting in both a negative NSE and RMS reduction. For the north component, the increases in both NSE and RMS reduction are also minor.

Discussions
In this study, we chose seven sites that show high inter-station correlations to conduct CME filtering for a local network and discarded the other two sites to avoid any possible bias due to misleading weighting factors. Figure 6 and Table 2 show that the rCME filtering is much better than that of the oCMEs in which the low-correlated sites are included. Since XIAG is weakly correlated with all the other stations, we further examined if there are any abnormal signals when the residual time series of XIAG are involved. For each component, we estimated the non-linear trends (curves in blue) for the two corrected time series and the weighted residuals (blue circles) of XIAG, and then computed the correlation coefficients between the weighted time series of XIAG and these two corrections. Figure 6a shows a weak correlation in the up component. In comparison, Figures 6a,b show that the weighted residual time series of XIAG are negatively correlated with the oCME correction in both the east and north components with coefficients of −0.732 and −0.626, respectively, indicating a large percent of oCME correction is contributed from the weighted residual time series of these two components of XIAG. Furthermore, compared with the rCME correction, a transient-like signal is present in the oCME corrections during the period from summer 2014 to the beginning of 2016 in the north component (Figure 5b). During the time from 2015 to 2016, a similar spurious transient is present in the oCME-corrected east component (Figure 5c). Although transient deformation may be present due to the intense drought [10,24], spurious transients will lead to misunderstanding of the physical processes. east components, respectively. The loading deformation is a combination of elastic response to NTAL, NTOL, and HYDL provided by Dill and Dobslaw [32]. The RMS% denotes the RMS reduction for each component. The curves in blue and red denote the modeled annual and semi-annual variations for rCMEs and loading deformation, respectively.

Discussions
In this study, we chose seven sites that show high inter-station correlations to conduct CME filtering for a local network and discarded the other two sites to avoid any possible bias due to misleading weighting factors. Figure 6 and Table 2 show that the rCME filtering is much better than that of the oCMEs in which the low-correlated sites are included. Since XIAG is weakly correlated with all the other stations, we further examined if there are any abnormal signals when the residual time series of XIAG are involved. For each component, we estimated the non-linear trends (curves in blue) for the two corrected time series and the weighted residuals (blue circles) of XIAG, and then computed the correlation coefficients between the weighted time series of XIAG and these two corrections. Figure 6a shows a weak correlation in the up component. In comparison, Figure 6a,b show that the weighted residual time series of XIAG are negatively correlated with the oCME correction in both the east and north components with coefficients of −0.732 and −0.626, respectively, indicating a large percent of oCME correction is contributed from the weighted residual time series of these two components of XIAG. Furthermore, compared with the rCME correction, a transient-like signal is present in the oCME corrections during the period from summer 2014 to the beginning of 2016 in the north component (Figure 5b). During the time from 2015 to 2016, a similar spurious transient is present in the oCME-corrected east component (Figure 5c). Although transient deformation may be present due to the intense drought [10,24], spurious transients will lead to misunderstanding of the physical processes.
Our results suggest that if any individual station affected by local effects is involved in the CME correction, spurious transient deformation may be present in the filtered time series ( Figure 5). It is notable that even for such a regional network, whose largest inter-station distance is less than 400 km, the supervised model, by simply ignoring sites with low inter-station correlation, it is possible to avoid bias from local effects and also make a significant reduction of error in the corrected residual time series.
Regardless of the bias from local effects of low-correlated sites, our results of rCMEs and oCMEs show remarkable annual and semi-annual variations (Figures 4, 7 and 8), particularly in the up and north components. Several previous studies have suggested that the seasonal variations in the up component of cGPS measurements in Yunnan are mainly caused by the hydrological mass redistribution, which has also been detected from GRACE gravity measurements [22,24,[27][28][29]. Such inference is physically reliable because the Yunnan province is located in the margin of the South Asian monsoon; the large amount of precipitation during the summer monsoon season will cause a significant deflection of the Earth's surface. Every year, the monsoon season occurs from June through September, nevertheless, the vertical displacements reach the lowest trough in October, forming a phase lag as shown in Figure 8. In comparison, Figure 8a shows that the atmospheric pressure applied on the surface has a phase lag later than that of the surface deformation. As a result, although the hydrological mass loading has the largest contribution to the magnitude, only when combined with NTAL and NTOL loadings does the combined loading deformation reach a good agreement with the CME (NSE = 0.55). Such agreement suggests that the seasonal variations in CME result from entire surface loading rather than the hydrological mass itself.
Such large-scale seasonal hydrological mass variations associated with the South Asian monsoon have caused seasonal deformations with larger annual amplitudes recorded in cGPS measurements in the Himalayas and Southeast Asia [8,9,11,34]. However, there are several sites showing mismatches with the loading deformation because of local effects such as groundwater irrigation [9] as well as flood [34]. These studies suggest that the cGPS measurements in horizontal directions are much more sensitive to the shorter spatial wavelengths of regional loads and local effects [27]. As such, a loading deformation model based on in situ measurements or surface mass redistribution of finer resolution is required to study any individual local effects.

Conclusions
In this study, we computed and compared two sets of CMEs for a local cGPS network in Norhtwest Yunnan province, China; the rCMEs were computed with the residual time series having high inter-station correlations, while a second CME, oCMEs, used all the residual time series. Our results show that after rCMEs filtering, the corrected residual time series have smaller variances and larger RMS reductions with respect to the original residual time series than those obtained using oCME filtering, indicating that rCMEs filtering is more efficient than oCMEs filtering. Since the local effects may be present in the low correlated stations, when the oCME filtering is applied, the filtered results show transient-like signals in both the east and north components. These results suggest that if any individual stations affected by local effects are involved in the CME correction, spurious transient deformation may be present in the filtered time series. As a result, it is necessary to evaluate the inter-station correlation and remove the possible noisy stations before conducting CME filtering.
We further compared the rCMEs with loading deformation induced by hydrological mass redistribution and non-tidal atmosphere pressure and ocean mass. The results showed that the combined loading deformation in the up component had the most remarkable seasonal variations, and only when these three loading terms were combined could the loading deformation in the up component reach a good agreement with the CME. Compared with HYDL itself, the NSE and RMS reduction with combined loading increased from 0.28 to 0.55 and from 15.19% to 33.4%, respectively, suggesting that the CME of the up component is mainly the result of the combined surface loading rather than hydrological water storage itself. Nevertheless, the loading deformation does not agree well with rCMEs of the east and north components, mostly because of the spatial resolution difference between these two types of measurements. In this case, one should carefully choose surface loading models to correct the raw cGPS time series if CME filtering is not conducted.