Effects of Spatiotemporal Filtering on the Periodic Signals and Noise in the GPS Position Time Series of the Crustal Movement Observation Network of China

Analysis of Global Positioning System (GPS) position time series and its common mode components (CMC) is very important for the investigation of GPS technique error, the evaluation of environmental loading effects, and the estimation of a realistic and unbiased GPS velocity field for geodynamic applications. In this paper, we homogeneously processed the daily observations of 231 Crustal Movement Observation Network of China (CMONOC) Continuous GPS stations to obtain their position time series. Then, we filtered out the CMC and evaluated its effects on the periodic signals and noise for the CMONOC time series. Results show that, with CMC filtering, peaks in the stacked power spectra can be reduced at draconitic harmonics up to the 14th, supporting the point that the draconitic signal is spatially correlated. With the colored noise suppressed by CMC filtering, the velocity uncertainty estimates for both of the two subnetworks, CMONOC-I (≈16.5 years) and CMONOC-II (≈4.6 years), are reduced significantly. However, the CMONOC-II stations obtain greater reduction ratios in velocity uncertainty estimates with average values of 33%, 38%, and 54% for the north, east, and up components. These results indicate that CMC filtering can suppress the colored noise amplitudes and improve the precision of velocity estimates. Therefore, a unified, realistic, and three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is given. Furthermore, contributions of environmental loading to the vertical CMC are also investigated and discussed. We find that the vertical CMC are reduced at 224 of the 231 CMONOC stations and 170 of them are with a root mean square (RMS) reduction ratio of CMC larger than 10%, confirming that environmental loading is one of the sources of CMC for the CMONOC height time series.


Introduction
The Crustal Movement Observation Network of China (CMONOC) is an important scientific infrastructure in China, which initially contained 27 Continuous Global Positioning System (GPS) stations in 1999 and then was extended to 260 GPS stations since 2011.CMONOC has provided useful datasets to facilitate the investigations of various geophysical phenomena, such as global plate motion [1], regional tectonic deformation [2], uplift of the Tibetan Plateau [3], land subsidence of Furthermore, periodic signals contribute to the temporal noise in GPS position time series [35,36] and both of them can be affected by CMC filtering.Therefore, CMC filtering has a further effect on velocity estimates and their uncertainties.Teferle et al. [37] argued that the main periodic terms (annual and semi-annual) need to be taken into account when computing the CMC so that the velocity bias due to CMC filtering can be mitigated.
GPS is one of the most important measurement techniques for crustal movement.Since the early 1990s, GPS has been applied to measure the crustal deformation in China.Wang et al. [38] firstly reported the present-day crustal deformation in China constrained by GPS measurements.Then, the crustal deformation in China was further investigated [2,39,40].However, these works are mainly based on campaign GPS measurements with which the transient signals are difficult to observe and the precision and accuracy of the GPS velocity estimates is limited.Owing to the extension of the CMONOC network to 260 continuous GPS stations and ≈2000 campaign GPS sites in 2011, the crustal deformation in China can be described with higher precision, as well as with more detailed information [3,[41][42][43].However, the noise characteristics of only a portion of the CMONOC stations were investigated with MLE in previous studies [42].Thus, a unified, realistic, and three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is still lacking.
In this paper, we homogeneously process the daily GPS observations of 231 CMONOC stations and filter out the CMC by using the PCA method in Section 2. By comparing the time series before and after CMC filtering, we evaluate the effects of CMC filtering on the stations' noise characteristics and velocity uncertainty estimates in Section 3. As a result, a unified three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is given in Section 4. Additionally, contributions of environmental loading to the CMC of the CMONOC network are quantified.Conclusions are drawn in Section 5.

GPS Data Processing
We collected the daily GPS observations of 260 CMONOC stations from February 1999 to August 2015.The observations were then processed with GAMIT/GLOBK software (Ver.10.5, http://geoweb.mit.edu/gg/)developed by Massachusetts Institute of Technology, Scripps Institution of Oceanography and Harvard University, the United States [44].Station coordinates, orbital parameters, tropospheric delay parameters, and Earth orientation parameters were estimated.Vienna Mapping Function 1 (VMF1) [45] and a priori zenith hydrostatic delay from the European Centre for Medium-Range Weather Forecasts (ECMWF) [46] were adopted for the correction of tropospheric delay.International Geomagnetic Reference Field 11 (IGRF11) [47] and ionospheric data from the Centre for Orbit Determination in Europe (CODE) [48] were used for the correction of high order ionospheric delay.The Finite Element Solutions 2004 (FES2004) model [49] was adopted for modeling the Ocean tide loading.The observations were weighted according to their elevation angles and post-fit phase residuals.The cutoff elevation angle was set as 10 • .Meanwhile, we also processed 73 globally evenly-distributed IGS stations as shown in Figure A1 with the same strategy.To speed up the calculation, the 260 CMONOC stations and 73 global IGS stations were divided into eight and two interwoven subnetworks, respectively.Afterwards, all of the daily loosely constrained solutions produced by GAMIT were combined and aligned to the IGb08 reference frame [50] by GLOBK with a six parameter (three translation and three rotation parameters) transformation which minimizes the coordinate differences of the globally distributed IGS stations.
We excluded 29 stations with observation time spans less than two years or with poor quality.Thus, only 231 CMONOC stations were analyzed.Their locations and time spans are shown in Figure 1.Among the 231 CMONOC stations, 27 of them were generally observed for ≈16.5 years and termed as CMONOC-I, and the other 204 stations were observed since 2011.0 and termed as CMONOC-II in this paper.The resulting position time series for the north, east and up (NEU) components were then handled separately and modeled with the following function [53]: = − = 0 for < 0 = 0.5 for = 0 = 1 for > 0 (3) in which R and R are the reference coordinate and epoch, respectively; is the velocity; is the epoch of GPS observations; is the number of offset events; is the ordinal number of the offset events, = 1, …, ; is the time difference between the epoch of GPS observations and the epoch of the th offset event as defined in Equation ( 2); H( ) is the Heaviside step function used The resulting position time series for the north, east and up (NEU) components were then handled separately and modeled with the following function [53]: ∆t l = 0 for t < t EQl ∆t l = t − t EQl for t > t EQl (4) in which x R and t R are the reference coordinate and epoch, respectively; v is the velocity; t is the epoch of GPS observations; n j is the number of offset events; j is the ordinal number of the offset events, j = 1, . . ., n j ; τ j is the time difference between the epoch of GPS observations t and the epoch of the jth offset event t j as defined in Equation ( 2); H(τ j ) is the Heaviside step function used to model the offsets as defined in Equation ( 3); b j is the jth offset; s k and c k are the coefficients of seasonal signals; ω k = 2π/p k ; p 1 = 1 year and p 2 = 1/2 year; n l is the number of logarithmic transients caused by earthquakes; l is the ordinal number of the post-seismic deformation events, l = 1, . . ., n l ; t EQl is the epoch of the lth earthquake event; ∆t l is the time since the lth earthquake occurred as defined in Equation ( 4); and T l is the transient timescale parameter set as T l = 1 year as suggested [53].
The offsets due to earthquakes or device changes were identified by checking the earthquake catalogues and station log files and then inspecting the time series.The earthquakes which lead to the offsets and/or post-seismic deformation of these CMONOC position time series are listed in Table A1.Some other offsets caused by unknown reasons were also identified by visual inspection.All identified offsets were then estimated with Equation (1).Moreover, time series of three north and nine east components show obvious post-seismic deformation caused by the 2011 Mw 9.0 Tohoku-Oki earthquake [54] and the 2001 Mw 7.8 Kokoxili earthquake [55].An example of stations affected by co-and post-seismic deformation is shown in Figure A2.Outliers in the residual time series were identified and cleaned according to the three interquartile range (IQR) thresholds [56].

Regional Filtering of CMC
Residual position time series in regional GPS network include CMC, local effects and random errors, in which CMC is one of the major spatiotemporal correlated error sources [25].To extract and remove the CMC, the PCA technique is introduced with an assumption that the spatiotemporal patterns of local effects and random errors are separable for those of the CMC [25].Later on, PCA has been widely used to extract CMC for regional GPS networks with even continental-scales [31][32][33].
Here we generally followed the PCA technique adopted by Dong et al. [25] and Serpelloni et al. [31] to extract and filter out the CMC from the CMONOC residual positions time series for NEU components.CMONOC-I and CMONOC-II were handled separately in this procedure due to significant differences between their time spans.Their NEU components were also handled separately.The main procedures are as follows: 1.
Obtain the residual position time series of all stations according to Equation (1) with trend, offsets, seasonal, and post-seismic terms removed and construct a m × n residual data matrix X, in which m and n are the number of epochs and stations, respectively.

3.
Decompose the symmetric matrix B as B = VΛV T and sort the eigenvectors to rank the eigenvalues in descending order.

4.
Consider a linear transformation A = XV, thus we have X = AV T .The columns of A and rows of V T are termed as principle components (PC) and spatial responses (SR), respectively.The ith PC and SR are termed as "mode i" together.

5.
Normalize each SR and scaling the corresponding PC by SR = SR/a and PC = a•PC, in which a is the component with the maximum absolute value in this SR.

6.
Select the significant PCs to calculate the CMC.CMC = A j V T j .The columns of A j are the first j PCs that are considered to be significant and rows of V T j are the corresponding SR.
To avoid the CMC from being contaminated, we considered the time series with abnormally large root mean square (RMS) (>3, >3, and >8 mm in north, east, and up directions, respectively) as outliers and we excluded them from the calculation of CMC.Moreover, if a station showed an anomalously large SR while its neighboring stations were with SRs close to zero, this station component was also excluded.Thus, the PCA may be iterated several times to remove outliers before the final operation.The excluded stations were filtered with the average CMC of the network [31].Furthermore, epochs with observations less than one third of the total number of stations were also excluded.
Figure 2 shows the cumulative percentage of the first 15 CMONOC-I PC eigenvalues and the first 25 CMONOC-II PC eigenvalues.For CMONOC-I, the first PCs represent 41%, 33%, and 35% of the total variances for the NEU components, respectively.For CMONOC-II, the first PCs represent 42%, 36%, and 37% of the total variances for the three components, respectively.We also find that, for both CMONOC-I and CMONOC-II, their second and even higher order PCs have significant eigenvalues, indicating that the power has spread into several PCs.
Figure 3 shows the first two modes for the NEU components of CMONOC-I.For all three components of CMONOC-I, their first PCs show significant responses with spatially uniform distribution whereas their second PCs show indistinctive responses at most of the stations.Figure 4 shows the first four modes for the up component of CMONOC-II.The first principle component (PC1) shows nearly uniform responses in space with an average of 71%, representing the mean part of CMC.This result is consistent to Serpelloni et al. [31], but completely different from Ming et al. [33].Further investigation is still needed to interpret the difference.The responses of PC2 and PC3 oscillate interactively in space, with a NE-SW gradual transition for the responses of PC2 and a SE-NW gradual transition for those of PC3.The PC2 and PC3 are likely to represent some parts of the CMC to fit the sub-regional signal.As for PC4, only a small subset of the stations shows significant response, indicating it just represents local effects.The first four modes for the north and east components of CMONOC-II show similar patterns to those of the up component and they are displayed in Figures A3 and A4, respectively.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 32 42%, 36%, and 37% of the total variances for the three components, respectively.We also find that, for both CMONOC-I and CMONOC-II, their second and even higher order PCs have significant eigenvalues, indicating that the power has spread into several PCs. Figure 3 shows the first two modes for the NEU components of CMONOC-I.For all three components of CMONOC-I, their first PCs show significant responses with spatially uniform distribution whereas their second PCs show indistinctive responses at most of the stations.Figure 4 shows the first four modes for the up component of CMONOC-II.The first principle component (PC1) shows nearly uniform responses in space with an average of 71%, representing the mean part of CMC.This result is consistent to Serpelloni et al. [31], but completely different from Ming et al. [33].Further investigation is still needed to interpret the difference.The responses of PC2 and PC3 oscillate interactively in space, with a NE-SW gradual transition for the responses of PC2 and a SE-NW gradual transition for those of PC3.The PC2 and PC3 are likely to represent some parts of the CMC to fit the sub-regional signal.As for PC4, only a small subset of the stations shows significant response, indicating it just represents local effects.The first four modes for the north and east components of CMONOC-II show similar patterns to those of the up component and they are displayed in Figures A3 and A4, respectively.Until now, the CMC has not been rigorously defined, as the nature and spatial distribution of the CMC are very complex and might be network-dependent.Thus, there is no generally accepted consensus on the selection of significant PCs to calculate CMC [25].For a small GPS network with a spatial scale of 5° × 5°, Dong et al. [25] considered a PC to be significant if more than 50% stations' responses are larger than 25% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.Whereas for a continental-scale GPS network of the Euro-Mediterranean region (between longitude 10°W-35°E and latitude 30°N-60°N), Serpelloni et al. [31] proposed a criterion that a significant PC should have uniform or gradual varied responses in space as the network is larger than the wavelength of the CMC and the CMC might partly leak into several PCs.The extents of CMONOC-I and CMONOC-II (between longitude 70°E-140°E and latitude 15°N-55°N) are close to that of the Euro-Mediterranean GPS network.Therefore, the criterion proposed by Serpelloni et al. [31] is followed in our analysis.
Here we slightly modify the quantitative criteria proposed by Dong et al. [25] so that the selection also meets the criterion proposed by Serpelloni et al. [31].For CMONOC-I, we consider a PC to be significant if more than 50% stations have SRs with absolute values larger than 50% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.For CMONOC-II, we consider a PC is significant if more than 30% stations have SRs with absolute values larger than 30% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.Consequently, only the first PCs are treated as CMC for the NEU components of CMONOC-I, while the first three PCs are selected for the three components of CMONOC-II.We set stricter selection criteria for CMONOC-I because its number of stations is relatively small (only 22 stations used for the NEU components, respectively) and its interstation distances are very large (usually >500 km), so it is difficult to distinguish the CMC and local effects if only a fraction of stations show significant responses.To illustrate the effects of CMC filtering, unfiltered and filtered height time series of three CMONOC-II stations and their CMC are shown in Figure Until now, the CMC has not been rigorously defined, as the nature and spatial distribution of the CMC are very complex and might be network-dependent.Thus, there is no generally accepted consensus on the selection of significant PCs to calculate CMC [25].For a small GPS network with a spatial scale of 5 • × 5 • , Dong et al. [25] considered a PC to be significant if more than 50% stations' responses are larger than 25% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.Whereas for a continental-scale GPS network of the Euro-Mediterranean region (between longitude 10 • W-35 • E and latitude 30 • N-60 • N), Serpelloni et al. [31] proposed a criterion that a significant PC should have uniform or gradual varied responses in space as the network is larger than the wavelength of the CMC and the CMC might partly leak into several PCs.The extents of CMONOC-I and CMONOC-II (between longitude 70 • E-140 • E and latitude 15 • N-55 • N) are close to that of the Euro-Mediterranean GPS network.Therefore, the criterion proposed by Serpelloni et al. [31] is followed in our analysis.
Here we slightly modify the quantitative criteria proposed by Dong et al. [25] so that the selection also meets the criterion proposed by Serpelloni et al. [31].For CMONOC-I, we consider a PC to be significant if more than 50% stations have SRs with absolute values larger than 50% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.For CMONOC-II, we consider a PC is significant if more than 30% stations have SRs with absolute values larger than 30% and the corresponding eigenvalues exceed 1% of the summation of all eigenvalues.Consequently, only the first PCs are treated as CMC for the NEU components of CMONOC-I, while the first three PCs are selected for the three components of CMONOC-II.We set stricter selection criteria for CMONOC-I because its number of stations is relatively small (only 22 stations used for the NEU components, respectively) and its interstation distances are very large (usually >500 km), so it is difficult to distinguish the CMC and local effects if only a fraction of stations show significant responses.To illustrate the effects of CMC filtering, unfiltered and filtered height time series of three CMONOC-II stations and their CMC are shown in Figure A5.These stations, Mulei (XJML, 43.8To illustrate the effects of the number of PCs used to estimate the CMC and to test the criteria proposed in this paper, Figure 5 compares the interstation correlation coefficients for the residuals of the CMONOC-II stations' up component time series without CMC filtering or filtered with CMCs derived from the first, the first three and the first five PCs, respectively.Compared with the unfiltered time series, the interstation correlation coefficients are reduced with CMC filtering using only the first PC.However, the average interstation correlation coefficients still range from −0.15 to 0.31, indicating that the first PC is not enough to filter out the CMC of the CMONOC-II up component time series.As for the residuals of the time series filtered with the first three and the first five PCs, their average interstation correlation coefficients are rather small and their differences are insignificant, indicating that the first three PCs are enough to filter out the CMC of the CMONOC-II up component time series.This result suggests that the criteria used to extract CMC is proper here.To illustrate the effects of the number of PCs used to estimate the CMC and to test the criteria proposed in this paper, Figure 5 compares the interstation correlation coefficients for the residuals of the CMONOC-II stations' up component time series without CMC filtering or filtered with CMCs derived from the first, the first three and the first five PCs, respectively.Compared with the unfiltered time series, the interstation correlation coefficients are reduced with CMC filtering using only the first PC.However, the average interstation correlation coefficients still range from −0.15 to 0.31, indicating that the first PC is not enough to filter out the CMC of the CMONOC-II up component time series.As for the residuals of the time series filtered with the first three and the first five PCs, their average interstation correlation coefficients are rather small and their differences are insignificant, indicating that the first three PCs are enough to filter out the CMC of the CMONOC-II up component time series.This result suggests that the criteria used to extract CMC is proper here.

Spatial Correlation
Figure 6 illustrates the interstation correlation coefficients between the residual time series of different CMONOC stations for the unfiltered and filtered time series, respectively.The interstation correlation coefficients for the NEU components of the unfiltered time series exhibit a similar behavior: the correlation coefficients decrease gradually with increasing interstation distances and with average correlation coefficients of 0.37, 0.30 and 0.32 for the three components, respectively.Moreover, the correlation coefficients between stations can be as large as 0.8 (average value around 0.5) with interstation distance of 500 km, and then decrease to zero with an interstation distance of about 4000 km for all the three components.These results validate the existence of spatial correlation in the unfiltered residual time series for all three component time series, which is consistent with other regional GPS networks, such as in Mexico [29] and also in China [57].
After CMC filtering with PCA, the average correlation coefficients are only 0.1 for short interstation distances (tens of kilometers), and reduce to zero rapidly with interstation distances of about 700 km for all the three components.That is to say, some local scale errors still remain in the filtered time series as we conservatively adopted only the first several PCs to represent the regional CMC.However, the significant reduction of the interstation correlation coefficients indicates that PCA is an effective filtering technique for suppressing the CMC in such a large regional GPS network and the criteria for choosing significant PCs used in this paper are appropriate.

Spatial Correlation
Figure 6 illustrates the interstation correlation coefficients between the residual time series of different CMONOC stations for the unfiltered and filtered time series, respectively.The interstation correlation coefficients for the NEU components of the unfiltered time series exhibit a similar behavior: the correlation coefficients decrease gradually with increasing interstation distances and with average correlation coefficients of 0.37, 0.30 and 0.32 for the three components, respectively.Moreover, the correlation coefficients between stations can be as large as 0.8 (average value around 0.5) with interstation distance of 500 km, and then decrease to zero with an interstation distance of about 4000 km for all the three components.These results validate the existence of spatial correlation in the unfiltered residual time series for all three component time series, which is consistent with other regional GPS networks, such as in Mexico [29] and also in China [57].
After CMC filtering with PCA, the average correlation coefficients are only 0.1 for short interstation distances (tens of kilometers), and reduce to zero rapidly with interstation distances of about 700 km for all the three components.That is to say, some local scale errors still remain in the filtered time series as we conservatively adopted only the first several PCs to represent the regional CMC.However, the significant reduction of the interstation correlation coefficients indicates that PCA is an effective filtering technique for suppressing the CMC in such a large regional GPS network and the criteria for choosing significant PCs used in this paper are appropriate.

Periodic Signal
By investigating and comparing the periodic signals in the unfiltered and filtered time series, we address several questions in this section.For instance, which periodic signals (e.g., the draconitic signal) are significant in the CMONOC time series?Do these signals contribute to the CMC?Is the draconitic signal spatially correlated in the CMONOC network?As the significant difference of observing time span between CMONOC-I (≈16.5 years) and CMONOC-II (≈4.6 years) may lead to different resolutions in detecting periodic signals, these two subnetworks are treated separately.Firstly, the power spectra of each time series was generated by Lomb-Scargle periodogram technique [58,59] with an oversampling factor of 4 [60].Then, the spectra for the NEU components of the unfiltered and filtered CMONOC-I stations were stacked respectively.Similarly, the spectra for the NEU components of the unfiltered and filtered CMONOC-II stations were stacked, respectively.These so-called stacked power spectra were also adopted by [10,13,16] to detect anomalous periodic signals in GPS position time series.

Periodic Signal
By investigating and comparing the periodic signals in the unfiltered and filtered time series, we address several questions in this section.For instance, which periodic signals (e.g., the draconitic signal) are significant in the CMONOC time series?Do these signals contribute to the CMC?Is the draconitic signal spatially correlated in the CMONOC network?As the significant difference of observing time span between CMONOC-I (≈16.5 years) and CMONOC-II (≈4.6 years) may lead to different resolutions in detecting periodic signals, these two subnetworks are treated separately.Firstly, the power spectra of each time series was generated by Lomb-Scargle periodogram technique [58,59] with an oversampling factor of 4 [60].Then, the spectra for the NEU components of the unfiltered and filtered CMONOC-I stations were stacked respectively.Similarly, the spectra for the NEU components of the unfiltered and filtered CMONOC-II stations were stacked, respectively.These so-called stacked power spectra were also adopted by [10,13,16] to detect anomalous periodic signals in GPS position time series.
In theory, if two periodic signals (with periods of T1 and T2, respectively, and T1 > T2) are to be distinguishable, the time span of the time series must be no shorter than T1T2/(T1−T2), according to the Rayleigh criterion [61].In reality, longer time series may be needed due to the noise in the observations [12].Thus, at least 25.4 years are needed to distinguish the annual and draconitic signals, which exceeds the time span of our GPS observations.Nevertheless, the CMONOC-I time series can distinguish the annual and draconitic signals from their second and higher harmonics (at least 12.7 In theory, if two periodic signals (with periods of T 1 and T 2 , respectively, and T 1 > T 2 ) are to be distinguishable, the time span of the time series must be no shorter than T 1 T 2 /(T 1 −T 2 ), according to the Rayleigh criterion [61].In reality, longer time series may be needed due to the noise in the observations [12].Thus, at least 25.4 years are needed to distinguish the annual and draconitic signals, which exceeds the time span of our GPS observations.Nevertheless, the CMONOC-I time series can distinguish the annual and draconitic signals from their second and higher harmonics (at least 12.7 years is needed) The CMONOC-II time series can distinguish from their 6th and higher harmonics (at least 4.2 years are needed).
Figure 7 shows the stacked power spectra of the unfiltered and filtered CMONOC-I residual time series with annual and semi-annual signals removed.The slopes of the stacked spectra in the log-log plot (equal to the spectral index κ) are close to −1 (FN, κ = −1) at low frequencies and close to 0 (WN, κ = 0) at high frequencies indicating a FN + WN character on average.To focus on the draconitic harmonics and some other periodic signals, a zoom-in of Figure 7 is shown in Figure 8. Peaks at the frequencies of the annual harmonics can be seen at the third annual harmonics for all the three components, but not visible for its higher harmonics.However, significant peaks at the frequencies of the first eight draconitic harmonics are clearly visible.Moreover, peaks can also be found at the frequencies of some other draconitic harmonics, such as the 10th for the up component and the 11th, 12th, and 14th for the north component.In addition, peaks are obvious at the frequencies of 24.75, 25.74, and 26.74 cpy, equal to periods of 14.76, 14.19, and 13.66 days, for all the three components (except the 25.74 cpy for the north component).These signals are likely to be the spurious long-period signals caused by the aliasing effects of the unmodeled short-period tidal signals (diurnal O1 tide and semidiurnal M2 tide) [11,62].The 13.66-day period spectra might also be related to the semi-monthly Mf tide.As for the filtered solution, several peaks at the frequencies of the draconitic harmonics are reduced in their stacked spectra, especially at the frequencies of the 6th draconitic harmonic for the north and east components and at the frequencies of the 8th draconitic harmonic for the up component.Moreover, the peaks at the frequencies of the 7th and 10th draconitic harmonics for the up component even become insignificant.In addition, the tri-annual signal and the periodic signals with frequencies of 24.75, 25.74, and 26.74 cpy are also reduced.Furthermore, some unknown signals are also reduced significantly, such as 3.32 cpy for the up component and 3.26 and 4.42 cpy for the north component.These unknown signals will be investigated in our future work.
Figure 9 shows the stacked power spectra of the CMONOC-II residual time series.Similar to CMONOC-I, the behavior of all the three components of CMONOC-II's spectra follows like FN at low frequencies and like WN at high frequencies.Several peaks at the frequencies of the draconitic harmonics can also be seen but with broader bands.A zoom-in of Figure 9 is also shown in Figure 10.As the time span of the CMONOC-II is ≈4.6 years, the annual and draconitic signals can only be distinguished from their 6th harmonic onwards in theory.Thus, the broad peaks for the third to 5th draconitic harmonics are surely leaked by the corresponding annual harmonics.Moreover, only peaks at the 6th, 8th, 12th, and 14th draconitic harmonics for the north component, 6th to 9th draconitic harmonics for the east component, and 6th draconitic harmonics for the up component are identified to be caused by the draconitic signal.Moreover, the stacked power spectra of the filtered CMONOC-II time series show that the peaks are reduced significantly around the frequencies of the annual and draconitic harmonics, and only the peaks at the frequencies of 6th to 9th draconitic harmonics are still distinguishable.Moreover, the signals with frequencies of 24.75 cpy for the north and up components, and frequency of 25.74 cpy for the east component are also reduced significantly.
frequencies of the first eight draconitic harmonics are clearly visible.Moreover, peaks can also be found at the frequencies of some other draconitic harmonics, such as the 10th for the up component and the 11th, 12th, and 14th for the north component.In addition, peaks are obvious at the frequencies of 24.75, 25.74, and 26.74 cpy, equal to periods of 14.76, 14.19, and 13.66 days, for all the three components (except the 25.74 cpy for the north component).These signals are likely to be the spurious long-period signals caused by the aliasing effects of the unmodeled short-period tidal signals (diurnal O1 tide and semidiurnal M2 tide) [11,62].The 13.66-day period spectra might also be related to the semi-monthly Mf tide.As for the filtered solution, several peaks at the frequencies of the draconitic harmonics are reduced in their stacked spectra, especially at the frequencies of the 6th draconitic harmonic for the north and east components and at the frequencies of the 8th draconitic harmonic for the up component.Moreover, the peaks at the frequencies of the 7th and 10th draconitic harmonics for the up component even become insignificant.In addition, the tri-annual signal and the periodic signals with frequencies of 24.75, 25.74, and 26.74 cpy are also reduced.Furthermore, some unknown signals are also reduced significantly, such as 3.32 cpy for the up component and 3.26 and 4.42 cpy for the north component.These unknown signals will be investigated in our future work.
Figure 9 shows the stacked power spectra of the CMONOC-II residual time series.Similar to CMONOC-I, the behavior of all the three components of CMONOC-II's spectra follows like FN at low frequencies and like WN at high frequencies.Several peaks at the frequencies of the draconitic harmonics can also be seen but with broader bands.A zoom-in of Figure 9 is also shown in Figure 10.As the time span of the CMONOC-II is ≈4.6 years, the annual and draconitic signals can only be distinguished from their 6th harmonic onwards in theory.Thus, the broad peaks for the third to 5th draconitic harmonics are surely leaked by the corresponding annual harmonics.Moreover, only peaks at the 6th, 8th, 12th, and 14th draconitic harmonics for the north component, 6th to 9th draconitic harmonics for the east component, and 6th draconitic harmonics for the up component are identified to be caused by the draconitic signal.Moreover, the stacked power spectra of the filtered CMONOC-II time series show that the peaks are reduced significantly around the frequencies of the annual and draconitic harmonics, and only the peaks at the frequencies of 6th to 9th draconitic harmonics are still distinguishable.Moreover, the signals with frequencies of 24.75 cpy for the north and up components, and frequency of 25.74 cpy for the east component are also reduced significantly.

Noise Characteristics
To investigate the effects of CMC filtering on the noise characteristics of CMONOC time series, the Hector software (Ver.1.7, http://segal.ubi.pt/hector/) with MLE was employed.Hector is a fast GPS position time series analysis software developed by Bos et al. [22] to meet the urgent demands of noise characterization for the rapidly increasing numbers of GPS stations and the accumulated GPS observations in recent years.A combination noise model of the PLN + WN is used as it is able to offer detailed description of the noise characteristics with the spectral index being estimated.Moreover, since CMONOC-I and CMONOC-II have significantly different observing time spans, they are investigated separately in this section.
Figure 11 illustrates the spatial distribution characteristics of the WN amplitudes, PLN amplitudes, spectral indices and velocity uncertainties of the unfiltered time series estimated with the PLN + WN model.The WN amplitudes of all the three components increase from north to south and reach their maximums at the southernmost China, indicating a general latitude dependence, which is consistent with [19,20].However, the latitude dependence is not obvious for the PLN amplitudes, spectral indices and velocity uncertainties, though there are several scattered southern stations with extraordinary large values.In addition, the WN amplitudes generally vary up to 2 mm and up to 6 mm for the horizontal and vertical components, respectively.The PLN amplitudes generally vary from 2 to 6 mm/year −κ/4 and from 6 to 18 mm/year −κ/4 for the horizontal and vertical components, respectively.The spectral indices are generally from −0.6 to −1.6 for the horizontal and from −0.6 to −1.2 for vertical components, respectively.The velocity uncertainties generally vary up to 0.8 mm/year and up to 2.4 mm/year for the horizontal and vertical components, respectively.
As a comparison, Figure 12 illustrates the spatial distribution characteristics of the WN amplitudes, PLN amplitudes, spectral indices and velocity uncertainties of the filtered time series estimated with the PLN + WN model.After CMC filtering, the PLN amplitudes and velocity uncertainties are generally reduced significantly.However, some stations in Yunnan province (between longitude 97°E-107°E and latitude 20°N-30°N) still show significant PLN amplitudes and velocity uncertainties.Most of these stations are excluded from the calculation of CMC because of

Noise Characteristics
To investigate the effects of CMC filtering on the noise characteristics of CMONOC time series, the Hector software (Ver.1.7, http://segal.ubi.pt/hector/) with MLE was employed.Hector is a fast GPS position time series analysis software developed by Bos et al. [22] to meet the urgent demands of noise characterization for the rapidly increasing numbers of GPS stations and the accumulated GPS observations in recent years.A combination noise model of the PLN + WN is used as it is able to offer detailed description of the noise characteristics with the spectral index being estimated.Moreover, since CMONOC-I and CMONOC-II have significantly different observing time spans, they are investigated separately in this section.
Figure 11 illustrates the spatial distribution characteristics of the WN amplitudes, PLN amplitudes, spectral indices and velocity uncertainties of the unfiltered time series estimated with the PLN + WN model.The WN amplitudes of all the three components increase from north to south and reach their maximums at the southernmost China, indicating a general latitude dependence, which is consistent with [19,20].However, the latitude dependence is not obvious for the PLN amplitudes, spectral indices and velocity uncertainties, though there are several scattered southern stations with extraordinary large values.In addition, the WN amplitudes generally vary up to 2 mm and up to 6 mm for the horizontal and vertical components, respectively.The PLN amplitudes generally vary from 2 to 6 mm/year −κ/4 and from 6 to 18 mm/year −κ/4 for the horizontal and vertical components, respectively.The spectral indices are generally from −0.6 to −1.6 for the horizontal and from −0.6 to −1.2 for vertical components, respectively.The velocity uncertainties generally vary up to 0.8 mm/year and up to 2.4 mm/year for the horizontal and vertical components, respectively.
As a comparison, Figure 12 illustrates the spatial distribution characteristics of the WN amplitudes, PLN amplitudes, spectral indices and velocity uncertainties of the filtered time series estimated with the PLN + WN model.After CMC filtering, the PLN amplitudes and velocity uncertainties are generally reduced significantly.However, some stations in Yunnan province (between longitude 97 • E-107 • E and latitude 20 • N-30 • N) still show significant PLN amplitudes and velocity uncertainties.Most of these stations are excluded from the calculation of CMC because of their abnormal non-linear variations and just filtered with the average CMC of the network.Their abnormal non-linear variations might be related to the droughts that occurred in Yunnan [8].Table 1 lists the mean and standard deviations of the WN amplitudes, PLN amplitudes, spectral indices, velocity uncertainties, and RMS of the CMONOC-I time series before and after CMC filtering.After the filtering, the average WN amplitudes are slightly reduced in horizontal but enlarged in the vertical direction, indicating that the WN component is probably correlated to site-specific errors that cannot be removed by CMC filtering.Nevertheless, their average PLN amplitudes show significant reductions for all the three components of the filtered time series, with average reduction ratios of 27%, 20%, and 27% for the NEU components, respectively.Moreover, the average spectral indices are slightly reduced for the north and up components.As for the velocity uncertainty estimates, they are reduced by 4%, 18%, and 16% on average for the NEU components, respectively, Furthermore, obvious RMS reductions can also be noticed from Table 1, indicating an improved signal-to-noise ratio.Table 1 lists the mean and standard deviations of the WN amplitudes, PLN amplitudes, spectral indices, velocity uncertainties, and RMS of the CMONOC-I time series before and after CMC filtering.After the filtering, the average WN amplitudes are slightly reduced in horizontal but enlarged in the vertical direction, indicating that the WN component is probably correlated to site-specific errors that cannot be removed by CMC filtering.Nevertheless, their average PLN amplitudes show significant reductions for all the three components of the filtered time series, with average reduction ratios of 27%, 20%, and 27% for the NEU components, respectively.Moreover, the average spectral indices are slightly reduced for the north and up components.As for the velocity uncertainty estimates, they are reduced by 4%, 18%, and 16% on average for the NEU components, respectively, Furthermore, obvious RMS reductions can also be noticed from Table 1, indicating an improved signal-to-noise ratio.Table 2 lists the mean and standard deviations of the WN amplitudes, PLN amplitudes, spectral indices, velocity uncertainties and RMS of the CMONOC-II time series before and after CMC filtering.The average WN amplitudes are also slightly reduced in horizontal but enlarged in vertical after the filtering, which are similar to the changes of WN amplitudes for the CMONOC-I time series.However, compared with the reduction ratios of CMONOC-I stations' PLN amplitude, the CMONOC-II stations' PLN amplitude show much larger reductions for all the three components after CMC filtering, with average reduction ratios of 47%, 46%, and 52% for the NEU components, respectively.Moreover, the CMONOC-II stations' velocity uncertainty estimates also obtain greater reduction ratios after CMC filtering, with average reduction ratios of 33%, 38%, and 54% for the NEU components, respectively.Table 2 lists the mean and standard deviations of the WN amplitudes, PLN amplitudes, spectral indices, velocity uncertainties and RMS of the CMONOC-II time series before and after CMC filtering.The average WN amplitudes are also slightly reduced in horizontal but enlarged in vertical after the filtering, which are similar to the changes of WN amplitudes for the CMONOC-I time series.However, compared with the reduction ratios of CMONOC-I stations' PLN amplitude, the CMONOC-II stations' PLN amplitude show much larger reductions for all the three components after CMC filtering, with average reduction ratios of 47%, 46%, and 52% for the NEU components, respectively.Moreover, the CMONOC-II stations' velocity uncertainty estimates also obtain greater reduction ratios after CMC filtering, with average reduction ratios of 33%, 38%, and 54% for the NEU components, respectively.To quantitatively estimate the effects of the filtering on geophysical parameters, we compare the station velocities between the filtered and unfiltered time series estimated with the PLN + WN model.Their differences are displayed in Figure 13.The majority of the differences are within ±0.3 mm/year and within ±0.6 mm/year for the horizontal and vertical components, respectively.To quantitatively estimate the effects of the filtering on geophysical parameters, we compare the station velocities between the filtered and unfiltered time series estimated with the PLN + WN model.Their differences are displayed in Figure 13.The majority of the differences are within ±0.3 mm/year and within ±0.6 mm/year for the horizontal and vertical components, respectively.

Characteristics of the Draconitic Harmonics
Figures 7-10 have demonstrated that CMC filtering with PCA is capable to decrease the draconitic harmonics for both CMONOC-I and CMONOC-II.Compared with CMONOC-I, CMONOC-II shows a more significant reduction.We attribute this difference to the fact that the spatial correlations between the CMONOC-II stations are more significant as their interstation distances are much smaller and, thus, more CMC components are reduced by the filtering.In addition, a considerable part of the CMONOC-I stations are not included in the calculation of CMC and just corrected with the average responses.In general, these results indicate that the CMC in the CMONOC time series contains the draconitic harmonic and some other periodic signals, thus supporting the conclusion drawn by Amiri-Simkooei [63] and Amiri-Simkooei et al. [12] that the draconitic signal is correlated in space.However, these results cannot exclude the possibility that the draconitic signal is partly due to site-specific effects as suggested by King and Watson [16] because there are still several peaks visible at the frequencies of the draconitic harmonics and the estimation of the power spectrum has uncertainty.Moreover, these results are consistent to the conclusions of Abraha et al. [17].

Characteristics of the Draconitic Harmonics
Figures 7-10 have demonstrated that CMC filtering with PCA is capable to decrease the draconitic harmonics for both CMONOC-I and CMONOC-II.Compared with CMONOC-I, CMONOC-II shows a more significant reduction.We attribute this difference to the fact that the spatial correlations between the CMONOC-II stations are more significant as their interstation distances are much smaller and, thus, more CMC components are reduced by the filtering.In addition, a considerable part of the CMONOC-I stations are not included in the calculation of CMC and just corrected with the average responses.In general, these results indicate that the CMC in the CMONOC time series contains the draconitic harmonic and some other periodic signals, thus supporting the conclusion drawn by Amiri-Simkooei [63] and Amiri-Simkooei et al. [12] that the draconitic signal is correlated in space.However, these results cannot exclude the possibility that the draconitic signal is partly due to site-specific effects as suggested by King and Watson [16] because there are still several peaks visible at the frequencies of the draconitic harmonics and the estimation of the power spectrum has uncertainty.Moreover, these results are consistent to the conclusions of Abraha et al. [17].

Changes of the Noise Characteristics Due to CMC Filtering
In Section 3.3, the CMONOC-I and CMONOC-II time series have been analyzed with the PLN + WN, respectively.For both of the CMONOC-I and CMONOC-II, their PLN amplitudes and velocity uncertainties are reduced after CMC filtering, consistent with prior studies [31,57,64].However, as the two subnetworks have significant differences in the time span of observations, the number of stations, and interstation distances, they may show some differences in their noise characteristics before and after CMC filtering.Before CMC filtering, the average WN amplitudes, PLN amplitudes and RMS of the CMONOC-II time series are lower than those of the CMONOC-I time series, which may be due to the fact that the quality of the latest observations are better than the early observations.While the average velocity uncertainties of the CMONOC-II time series are larger than those of the CMONOC-I time series as their time span of observations are much shorter.Comparing the effects of CMC filtering on CMONOC-I and CMONOC-II time series, one can find that the CMONOC-II time series show much larger reductions in PLN amplitude and velocity uncertainty estimates for all three components, indicating that the CMONOC-II time series benefit more from CMC filtering.This may be because that the CMONOC-II subnetwork has shorter interstation distances and the first three PCs are used for the calculation of its CMC, whereas only the first PC is used to calculate of the CMONOC-I subnetwork.
Table A2 lists the velocity estimates and their uncertainties for the unfiltered and filtered NEU component time series for the CMONOC stations estimated with the PLN + WN model.Moreover, the horizontal and vertical velocity field of the CMC filtered CMONOC time series estimated with the PLN + WN model are shown in Figures 14 and 15, respectively.The velocity estimates of CMONOC stations provide abundant information on the tectonic deformation over the mainland China and we will focus on it in our future work.

Changes of the Noise Characteristics Due to CMC Filtering
In Section 3.3, the CMONOC-I and CMONOC-II time series have been analyzed with the PLN + WN, respectively.For both of the CMONOC-I and CMONOC-II, their PLN amplitudes and velocity uncertainties are reduced after CMC filtering, consistent with prior studies [31,57,64].However, as the two subnetworks have significant differences in the time span of observations, the number of stations, and interstation distances, they may show some differences in their noise characteristics before and after CMC filtering.Before CMC filtering, the average WN amplitudes, PLN amplitudes and RMS of the CMONOC-II time series are lower than those of the CMONOC-I time series, which may be due to the fact that the quality of the latest observations are better than the early observations.While the average velocity uncertainties of the CMONOC-II time series are larger than those of the CMONOC-I time series as their time span of observations are much shorter.Comparing the effects of CMC filtering on CMONOC-I and CMONOC-II time series, one can find that the CMONOC-II time series show much larger reductions in PLN amplitude and velocity uncertainty estimates for all three components, indicating that the CMONOC-II time series benefit more from CMC filtering.This may be because that the CMONOC-II subnetwork has shorter interstation distances and the first three PCs are used for the calculation of its CMC, whereas only the first PC is used to calculate CMC of the CMONOC-I subnetwork.
Table A2 lists the velocity estimates and their uncertainties for the unfiltered and filtered NEU component time series for the CMONOC stations estimated with the PLN + WN model.Moreover, the horizontal and vertical velocity field of the CMC filtered CMONOC time series estimated with the PLN + WN model are shown in Figures 14 and 15, respectively.The velocity estimates of CMONOC stations provide abundant information on the tectonic deformation over the mainland China and we will focus on it in our future work.

Environmental Loading Effects on the CMC
Environmental loading is one of the potential candidates for the CMC (Dong et al., 2006).In this section, the contributions of environmental loading to the CMC of the CMONOC network is quantified.As environmental loading effects are usually most significant in vertical, only the CMONOC height time series are analyzed.We firstly obtained the vertical environmental loading displacements for the CMONOC stations predicted with loading products from German Research Centre for Geosciences (GFZ) [66], in which atmospheric loading, non-tidal oceanic loading and hydrological loading effects are taken into account.Environmental loading displacements at a station are shown in Figure A6 as an example.
We then calculated and compared the RMS of CMC for the CMONOC height time series before and after environmental loading correction as shown in Figure 16a,b.RMS reduction ratios of the corresponding CMC are shown in Figure 16c.For the CMONOC height time series before and after environmental loading correction, the RMS of CMC range from 1.4 mm to 5.0 mm and from 0.3 mm to 4.7 mm with medians of 3.4 mm and 2.7 mm, respectively.In addition, the according RMS reduction ratios of CMC range between −9.1% and 77.8%, with a median of 18.4% at station Wuqia (XJWU, 39.7°N, 75.2°E).To illustrate and compare the time series before and after environmental loading correction, station XJWU is selected as an example as shown in Figure 17.Before loading correction, the RMS of station XJWU's unfiltered and filtered height time series are 4.1 mm and 2.8 mm, respectively.After loading correction, the RMS of station XJWU's unfiltered and filtered height time series are 3.6 mm and 2.6 mm, respectively.Accordingly, the RMS of CMC for station XJWU's height time series decrease from 3.0 mm to 2.4 mm after the environmental loading correction.
In general, 224 of the 231 stations show a decrease in the RMS of CMC after the environmental loading correction and 170 of them are with a RMS reduction ratio of CMC larger than 10%, confirming that environmental loading is one of the sources of CMC for the CMONOC height time series, which is consistent with Zhu et al. [6].

Environmental Loading Effects on the CMC
Environmental loading is one of the potential candidates for the CMC (Dong et al., 2006).In this section, the contributions of environmental loading to the CMC of the CMONOC network is quantified.As environmental loading effects are usually most significant in vertical, only the CMONOC height time series are analyzed.We firstly obtained the vertical environmental loading displacements for the CMONOC stations predicted with loading products from German Research Centre for Geosciences (GFZ) [66], in which atmospheric loading, non-tidal oceanic loading and hydrological loading effects are taken into account.Environmental loading displacements at a station are shown in Figure A6 as an example.
We then calculated and compared the RMS of CMC for the CMONOC height time series before and after environmental loading correction as shown in Figure 16a,b.RMS reduction ratios of the corresponding CMC are shown in Figure 16c.For the CMONOC height time series before and after environmental loading correction, the RMS of CMC range from 1.4 mm to 5.0 mm and from 0.3 mm to 4.7 mm with medians of 3.4 mm and 2.7 mm, respectively.In addition, the according RMS reduction ratios of CMC range between −9.1% and 77.8%, with a median of 18.4% at station Wuqia (XJWU, 39.7 • N, 75.2 • E).To illustrate and compare the time series before and after environmental loading correction, station XJWU is selected as an example as shown in Figure 17.Before loading correction, the RMS of station XJWU's unfiltered and filtered height time series are 4.1 mm and 2.8 mm, respectively.After loading correction, the RMS of station XJWU's unfiltered and filtered height time series are 3.6 mm and 2.6 mm, respectively.Accordingly, the RMS of CMC for station XJWU's height time series decrease from 3.0 mm to 2.4 mm after the environmental loading correction.
In general, 224 of the 231 stations show a decrease in the RMS of CMC after the environmental loading correction and 170 of them are with a RMS reduction ratio of CMC larger than 10%, confirming that environmental loading is one of the sources of CMC for the CMONOC height time series, which is consistent with Zhu et al. [6].

Conclusions
Daily position residual time series of 231 CMONOC GPS stations have been analyzed to investigate the effects of CMC filtering with PCA on the periodic signals and noise characteristics.The source of the draconitic signal in the CMONOC GPS position time series was discussed.The noise characteristics of the CMONOC-I (27 stations, ≈16.5 years) and CMONOC-II (204 stations, ≈4.6 years) stations before and after CMC filtering have been compared, and possible reasons for their differences have been discussed, such as the differences in time span of observations, number of stations, and interstation distances.Therefore, a unified, realistic, and three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is given.Furthermore, contributions of environmental loading to the CMC of CMONOC were quantified as well.The following are the main findings: 1.The stacked power spectra of the CMC-filtered CMONOC residual time series show that peaks are reduced significantly at frequencies of tri-annual, draconitic harmonics up to 14th, and 24.75, 25.74, and 26.74 cpy, indicating that the CMC of the CMONOC network contains draconitic

Conclusions
Daily position residual time series of 231 CMONOC GPS stations have been analyzed to investigate the effects of CMC filtering with PCA on the periodic signals and noise characteristics.The source of the draconitic signal in the CMONOC GPS position time series was discussed.The noise characteristics of the CMONOC-I (27 stations, ≈16.5 years) and CMONOC-II (204 stations, ≈4.6 years) stations before and after CMC filtering have been compared, and possible reasons for their differences have been discussed, such as the differences in time span of observations, number of stations, and interstation distances.Therefore, a unified, realistic, and three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is given.Furthermore, contributions of environmental loading to the CMC of CMONOC were quantified as well.The following are the main findings: 1.The stacked power spectra of the CMC-filtered CMONOC residual time series show that peaks are reduced significantly at frequencies of tri-annual, draconitic harmonics up to 14th, and 24.75, 25.74, and 26.74 cpy, indicating that the CMC of the CMONOC network contains draconitic harmonics and some other periodic signals.These results support the view that the draconitic

Conclusions
Daily position residual time series of 231 CMONOC GPS stations have been analyzed to investigate the effects of CMC filtering with PCA on the periodic signals and noise characteristics.The source of the draconitic signal in the CMONOC GPS position time series was discussed.The noise characteristics of the CMONOC-I (27 stations, ≈16.5 years) and CMONOC-II (204 stations, ≈4.6 years) stations before and after CMC filtering have been compared, and possible reasons for their differences have been discussed, such as the differences in time span of observations, number of stations, and interstation distances.Therefore, a unified, realistic, and three-dimensional CMONOC GPS velocity field estimated with the consideration of colored noise is given.Furthermore, contributions of environmental loading to the CMC of CMONOC were quantified as well.The following are the main findings: 1.
The stacked power spectra of the CMC-filtered CMONOC residual time series show that peaks are reduced significantly at frequencies of tri-annual, draconitic harmonics up to 14th, and 24.75, 25.74, and 26.74 cpy, indicating that the CMC of the CMONOC network contains draconitic harmonics and some other periodic signals.These results support the view that the draconitic signal is spatially correlated.However, the possibility that the draconitic signal is caused by site-specific effects cannot be ruled out because of the weakened, but still visible, peaks at the frequencies of the draconitic harmonics.

2.
For the unfiltered time series, the velocity uncertainties of the CMONOC stations estimated with an assumption of the PLN + WN model generally vary up to 0.8 mm/year and up to 2.4 mm/year for the horizontal and vertical components, respectively.After CMC filtering, the average white noise amplitudes are slightly reduced in horizontal but enlarged in vertical for both the CMONOC-I (≈16.Table A1.The earthquakes which caused co-seismic and/or post-seismic deformation of the CMONOC stations (http://www.globalcmt.org/CMTsearch.html[25,26]).

32 Figure 1 .
Figure 1.Location map of the selected 231 continuous GPS stations investigated in this study.Squares and circles indicate CMONOC-I and CMONOC-II stations, respectively.Plus and multiplication signs indicate stations with their east and north components affected by post-seismic deformation (PSD), respectively.Several earthquakes in China and surrounding regions are plotted with information from the Global CMT Catalog Search (http://www.globalcmt.org/CMTsearch.html[51,52]).

Figure 1 .
Figure 1.Location map of the selected 231 continuous GPS stations investigated in this study.Squares and circles indicate CMONOC-I and CMONOC-II stations, respectively.Plus and multiplication signs indicate stations with their east and north components affected by post-seismic deformation (PSD), respectively.Several earthquakes in China and surrounding regions are plotted with information from the Global CMT Catalog Search (http://www.globalcmt.org/CMTsearch.html[51,52]).

Figure 2 .
Figure 2. Cumulative percentage of the first 15 CMONOC-IPC eigenvalues (left panel) and the first 25 CMONOC-II PC eigenvalues (right panel).

Figure 2 .
Figure 2. Cumulative percentage of the first 15 CMONOC-IPC eigenvalues (left panel) and the first 25 CMONOC-II PC eigenvalues (right panel).

Figure 3 .
Figure 3. Temporal and spatial responses of the first two PCs for the north (upper panels), east (middle panels), and up (lower panels) components of CMONOC-I.Black dots indicate the excluded stations.

Figure 3 .
Figure 3. Temporal and spatial responses of the first two PCs for the north (upper panels), east (middle panels), and up (lower panels) components of CMONOC-I.Black dots indicate the excluded stations.

Figure 4 .
Figure 4. Temporal and spatial responses of the first four PCs for the up component of the CMONOC-II.Black dots indicate the excluded stations.

Figure 4 .
Figure 4. Temporal and spatial responses of the first four PCs for the up component of the CMONOC-II.Black dots indicate the excluded stations.

32 Figure 5 .
Figure 5. Correlation coefficients between the residual time series of different CMONOC-II stations' up component time series as a function of interstation distances in kilometers (left panels) and histograms of the interstation correlation coefficients (right panels).The interstation correlation coefficients of the unfiltered residual time series or residual time series filtered with CMCs derived from the first, the first three, and the first five PCs are shown as blue, green, orange, and purple dots.Moreover, they are smoothed by a boxcar smoother with a full width of 50 km and shown as the dark blue, dark green, dark orange, and dark purple lines, respectively.

Figure 5 .
Figure 5. Correlation coefficients between the residual time series of different CMONOC-II stations' up component time series as a function of interstation distances in kilometers (left panels) and histograms of the interstation correlation coefficients (right panels).The interstation correlation coefficients of the unfiltered residual time series or residual time series filtered with CMCs derived from the first, the first three, and the first five PCs are shown as blue, green, orange, and purple dots.Moreover, they are smoothed by a boxcar smoother with a full width of 50 km and shown as the dark blue, dark green, dark orange, and dark purple lines, respectively.

32 Figure 6 .
Figure 6.Correlation coefficients between the residual time series of different CMONOC stations as a function of interstation distances in kilometers (left panels) and histograms of the interstation correlation coefficients (right panels).The interstation correlation coefficients of the unfiltered and filtered time series are shown as blue and orange dot marks, respectively.Moreover, they are smoothed by a boxcar smoother with a full width of 50 km and shown as the blue and red lines, respectively.

Figure 6 .
Figure 6.Correlation coefficients between the residual time series of different CMONOC stations as a function of interstation distances in kilometers (left panels) and histograms of the interstation correlation coefficients (right panels).The interstation correlation coefficients of the unfiltered and filtered time series are shown as blue and orange dot marks, respectively.Moreover, they are smoothed by a boxcar smoother with a full width of 50 km and shown as the blue and red lines, respectively.

Figure 7 .
Figure 7. Stacked power spectra of the unfiltered (blue) and filtered (red) CMONOC-I residual time series with annual and semi-annual signals removed.The stacked power spectra for the up and east components have been shifted upward for clarity (by factors of 5000 and 100, respectively).Vertical black and gray lines indicate the annual (1 cpy) and draconitic (1.04 cpy) oscillations and their harmonics up to the 4th.Vertical green lines indicate the frequencies of 24.75, 25.74, and 26.74 cpy.

Figure 7 .
Figure 7. Stacked power spectra of the unfiltered (blue) and filtered (red) CMONOC-I residual time series with annual and semi-annual signals removed.The stacked power spectra for the up and east components have been shifted upward for clarity (by factors of 5000 and 100, respectively).Vertical black and gray lines indicate the annual (1 cpy) and draconitic (1.04 cpy) oscillations and their harmonics up to the 4th.Vertical green lines indicate the frequencies of 24.75, 25.74, and 26.74 cpy.

Figure 8 .
Figure 8.A zoom-in of Figure 7. Vertical dashed gray lines indicate harmonics of the draconitic oscillations (1.04 cpy) up to the 14th.

Figure 9 .
Figure 9. Same as Figure 7, but for the CMONOC-II time series.

Figure 8 . 32 Figure 8 .
Figure 8.A zoom-in of Figure 7. Vertical dashed gray lines indicate harmonics of the draconitic oscillations (1.04 cpy) up to the 14th.

Figure 9 .
Figure 9. Same as Figure 7, but for the CMONOC-II time series.

Figure 9 .
Figure 9. Same as Figure 7, but for the CMONOC-II time series.

Figure 10 .
Figure 10.Same as Figure 8, but for the CMONOC-II time series.

Figure 10 .
Figure 10.Same as Figure 8, but for the CMONOC-II time series.

32 Figure 11 .
Figure 11.The WN amplitudes for the NEU components (a-c); PLN amplitudes for the NEU components (d-f); spectral indices for the NEU components (g-i) and velocity uncertainties for the NEU components (j-l) of the unfiltered time series estimated with the PLN + WN model.

Figure 11 .
Figure 11.The WN amplitudes for the NEU components (a-c); PLN amplitudes for the NEU components (d-f); spectral indices for the NEU components (g-i) and velocity uncertainties for the NEU components (j-l) of the unfiltered time series estimated with the PLN + WN model.

Figure 12 .
Figure 12.Same as Figure 11, but for the CMC filtered time series.

Figure 12 .
Figure 12.Same as Figure 11, but for the CMC filtered time series.

Figure 13 .
Figure 13.Velocity differences between the filtered and unfiltered time series.

Figure 13 .
Figure 13.Velocity differences between the filtered and unfiltered time series.

Figure 14 .
Figure 14.Horizontal velocity field of the CMC filtered CMONOC stations estimated with the PLN + WN model.The reference frame is IGb08.Error ellipses reflect a 95% confidence interval.The tectonic blocks are generated according to Zhang et al. [65].

Figure 14 .
Figure 14.Horizontal velocity field of the CMC filtered CMONOC stations estimated with the PLN + WN model.The reference frame is IGb08.Error ellipses reflect a 95% confidence interval.The tectonic blocks are generated according to Zhang et al. [65].

Figure 15 .
Figure 15.Same as Figure 14, but for the vertical velocity field.

Figure 15 .
Figure 15.Same as Figure 14, but for the vertical velocity field.

32 Figure 16 .
Figure 16.RMS of the CMC before (a) and after (b) environmental loading correction and the corresponding RMS reduction ratios (c) for the vertical displacement time series of the 231 CMONOC stations.

Figure 17 .
Figure 17.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the vertical displacement time series of station Wuqia (XJWU, 39.7°N, 75.2°E) before (a) and after (b) environmental loading correction.

Figure 16 . 32 Figure 16 .
Figure 16.RMS of the CMC before (a) and after (b) environmental loading correction and the corresponding RMS reduction ratios (c) for the vertical displacement time series of the 231 CMONOC stations.

Figure 17 .
Figure 17.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the vertical displacement time series of station Wuqia (XJWU, 39.7°N, 75.2°E) before (a) and after (b) environmental loading correction.

Figure 17 .
Figure 17.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the vertical displacement time series of station Wuqia (XJWU, 39.7 • N, 75.2 • E) before (a) and after (b) environmental loading correction.

Figure A1 .
Figure A1.Location map of the globally distributed 73 IGS stations which are generally selected from the core network of IGb08 reference frame (ftp://igs.org/pub/station/coord/IGb08_core.txt).The black and blue dots indicate two subnetwork stations, respectively.The red pentagrams indicate the common stations of these two subnetworks.

Figure A2 .
Figure A2.An example of stations affected by co-and post-seismic deformation.Station Hegang (HLHG) is located in Northeast China and about 1500 km away from the epicenter (red solid circle) of the Mw 9.0 Tohoku-Oki Earthquake.The gray curve indicate the fitting curve of east component modeled with Equation (1).The black curves indicate the fitting curves of the NEU components modeled with Equation (1), but without the logarithmic term.The post-seismic deformation is not significant for north or up components of HLHG and, thus, there is no need to be taken into consideration in the function model.As for the east component of HLHG, obvious post-seismic deformation can be noticed and, thus, a logarithmic term needs to be included in the function model.

Figure A2 .
Figure A2.An example of stations affected by co-and post-seismic deformation.Station Hegang (HLHG) is located in Northeast China and about 1500 km away from the epicenter (red solid circle) of the Mw 9.0 Tohoku-Oki Earthquake.The gray curve indicate the fitting curve of east component modeled with Equation (1).The black curves indicate the fitting curves of the NEU components modeled with Equation (1), but without the logarithmic term.The post-seismic deformation is not significant for north or up components of HLHG and, thus, there is no need to be taken into consideration in the function model.As for the east component of HLHG, obvious post-seismic deformation can be noticed and, thus, a logarithmic term needs to be included in the function model.

Figure A3 .
Figure A3.Same as Figure 4, but for the north component.Figure A3.Same as Figure 4, but for the north component.

Figure A3 .
Figure A3.Same as Figure 4, but for the north component.Figure A3.Same as Figure 4, but for the north component.

Figure A4 .
Figure A4.Same as Figure 4, but for the east component.

Figure A5 .
Figure A5.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the selected three vertical position time series of the CMONOC-II stations, in which station XJML, HAHB and HECX are characterized by the maximum, medium and minimum RMS reduction in CMONOC-II after the filtering, respectively.

Figure A4 . 32 Figure A4 .
Figure A4.Same as Figure 4, but for the east component.

Figure A5 .
Figure A5.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the selected three vertical position time series of the CMONOC-II stations, in which station XJML, HAHB and HECX are characterized by the maximum, medium and minimum RMS reduction in CMONOC-II after the filtering, respectively.

Figure A5 .
Figure A5.Unfiltered time series (red, shifted upward by 150 mm) and their residual time series (orange, shifted upward by 110 mm) obtained with Equation (1), common mode error (green, shifted upward by 70 mm), and filtered time series (blue, shifted upward by 30 mm), and their residual time series (purple) obtained with Equation (1) for the selected three vertical position time series of the CMONOC-II stations, in which station XJML, HAHB and HECX are characterized by the maximum, medium and minimum RMS reduction in CMONOC-II after the filtering, respectively.

Table 1 .
Mean and standard deviations of the WN amplitude (σ W N ), PLN amplitude (σ PL ), spectral index, and velocity uncertainty (σ vel. ) and the RMS of the unfiltered and filtered CMONOC-I time series estimated with the PLN + WN model.

Table 2 .
Same as Table 1, but for the CMONOC-II time series.

Table 1 .
Mean and standard deviations of the WN amplitude (σ ), PLN amplitude (σ ), spectral index, and velocity uncertainty (σ .) and the RMS of the unfiltered and filtered CMONOC-I time series estimated with the PLN + WN model.

Table 2 .
Same as Table 1, but for the CMONOC-II time series.
5 years) and CMONOC-II (≈4.6 years) stations.Nevertheless, the average power-law noise amplitudes are significantly suppressed by CMC filtering.Therefore, the velocity uncertainty estimates of north, east and up components for both the CMONOC-I and CMONOC-II stations are reduced.Compared with CMONOC-I, the CMONOC-II stations obtain greater reduction ratios in velocity uncertainty estimates with average values of 33%, 38%, and 54% for the north, east, and up components, respectively.These results indicate that CMC filtering can suppress the colored noise amplitudes and improve the precision of velocity estimates.3.After environmental loading correction, vertical CMC are reduced at 224 of the 231 CMONOC stations.In addition, 170 of them are with an RMS reduction ratio of CMC larger than 10%, confirming that environmental loading is one of the sources of CMC for the CMONOC height time series.

Table A2 .
Three-dimensional velocity estimates and their uncertainties of the CMC filtered 231 CMONOC stations estimated with the PLN + WN model.The reference frame is IGb08.

Table A2 .
Three-dimensional velocity estimates and their uncertainties of the CMC filtered 231 CMONOC stations estimated with the PLN + WN model.The reference frame is IGb08.