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

Abstract

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.

1. 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 the North China Plain [4], environmental loading effects [5,6,7] and hydrological cycle [8]. Analyzing the GPS position time series of CMONOC is vital for unbiased velocity estimates, as well as realistic velocity uncertainty estimates, which is a necessity for appropriate interpretations of the above mentioned geophysical phenomena. In addition, analyzing the GPS position time series of CMONOC, a long-running, continental-scale and high-density GPS network, can provide a wealth of information on the periodic signals, noise characteristics, and sources of the common mode component (CMC) for GPS position time series, which can be a reference for similar studies over other regions.

Previous studies have reported that periodic signals are present in GPS position time series, such as the seasonal signals [9], draconitic signal and its harmonics [10], spurious long-period signals caused by unmodeled short-period signals [11], and other periodic signals due to unknown reasons [12]. The draconitic signal was found by Ray et al. [10], with significant peaks at 1.04 cpy and its higher harmonics in the stacked power spectra of 167 International GNSS Service (IGS) stations. It is considered to be mainly related to the repeat times of the GPS constellation in inertial space with respect to the Sun [10], although there are controversial opinions on its sources. For instance, Tregoning and Watson [13] showed that unmodeled diurnal and semidiurnal atmospheric tides can result in spurious draconitic signals. Rodriguez-Solano et al. [14] found that the peak of the 6th draconitic harmonic of the north component can be reduced with the Earth radiation pressure taken into account in GPS data processing. Li et al. [15] found that the draconitic signal could be reduced when considering minor ocean tides. Additionally, King and Watson [16] suggested another possibility that the draconitic signal might be partly due to multipath-related errors. Abraha et al. [17] reported that the draconitic signal is mostly from satellite-related effects with a non-zero contribution of site-specific errors. As the orbital errors and unmodeled geophysical tide loading effects are spatially correlated while the multipath effect is a site-specific error source, we will investigate the spatial correlation of the draconitic signal in this study so that its source can be more reliably speculated.

It is also well known that GPS position time series are temporally correlated [18,19]. The temporal correlation indicates a colored noise which can be described as a power-law noise (PLN) with a spectral index $\kappa$: $P\left(f\right)\propto f^{\kappa }$. White noise (WN, $\kappa =0$), flicker noise (FN, $\kappa =-1$), and random walk noise (RWN, $\kappa =-2$) are three special cases of power-law noise. If the temporal correlation is not taken into consideration and only WN is assumed, the velocity uncertainties might be underestimated with factors of 5–10 [20]. Therefore, a proper description of the noise characteristics is essential to the estimation of realistic velocity uncertainties and to the proper interpretations of geophysical phenomena. To describe the noise characteristics of the GPS position time series, the maximum likelihood estimation (MLE) is commonly used [20,21,22] and a combination of power-law noise plus white noise (PLN+WN) is assumed as the noise model [21].

Apart from the temporal correlation, existing investigations have acknowledged that GPS position time series are spatially correlated in regional GPS networks [23,24]. The common mode position variations, termed as common mode components (CMC), are likely caused by mismodeled and unmodeled environmental loading, GPS errors, and reference frame definition [25]. Filtering out the CMC can improve the signal-to-noise ratio of GPS position time series. Several approaches have been proposed to remove the CMC. One of the most commonly used approaches is stacking filtering, which was first introduced by Wdowinski et al. [23] with an assumption of an equivalent CMC for all the stations and then it was applied to a lot of regions around the world, such as Southern California [26], Nevada-California [27], and the Dead Sea Fault [28]. However, the stacking filtering may not be appropriate to continental-scale networks in which the CMC becomes inhomogeneous. To solve this problem, this approach was modified with the consideration of interstation distances and correlations [29,30]. Another widely used approach to filter out CMC is principle component analysis (PCA) [25]. The PCA approach allows for unequal CMC in the network and it has been applied to continental-scale networks with good performances [31,32,33]. In addition, Teferle et al. [34] mentioned that the stacking filtering is equivalent to the first principle component time series.

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.

2. Materials and Methods

2.1. 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]:

$$\begin{array}{ll}x\left(t\right)=& x_{\mathrm{R}}+v\left(t-t_{\mathrm{R}}\right)+{\displaystyle \sum _{j=1}^{n_j}}{b}_{j}H\left({\tau }_j\right)\\ & +{\displaystyle \sum _{k=1}^2}\left[s_k\sin \left({\omega }_{k}t\right)+c_k\cos \left({\omega }_{k}t\right)\right]+{\displaystyle \sum _{l=1}^{n_l}}{a}_l\log \left(1+\mathsf{\Delta }{t}_l/T_l\right)+\epsilon \end{array}$$

(1)

$${\tau }_j= t-t_j$$

(2)

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

(3)

$$\begin{array}{c}\mathsf{\Delta }{t}_l=0\mathrm{for}t<t_{EQl}\\ \mathsf{\Delta }{t}_l=t-t_{EQl}\mathrm{for}t>t_{EQl}\end{array}$$

(4)

in which $x_{\mathrm{R}}$ and $t_{\mathrm{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$; ${\tau }_j$ is the time difference between the epoch of GPS observations $t$ and the epoch of the $j$th offset event $t_j$ as defined in Equation (2); H(${\tau }_j$) is the Heaviside step function used to model the offsets as defined in Equation (3); $b_j$ is the $j$th offset; $s_k$ and $c_k$ are the coefficients of seasonal signals; ${\omega }_k=2\pi /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 $l\mathrm{th}$ earthquake event; $\mathsf{\Delta }{t}_l$ is the time since the $l\mathrm{th}$ 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. The earthquakes which caused co-seismic and/or post-seismic deformation of the CMONOC stations (http://www.globalcmt.org/CMTsearch.html [51,52]).

Time (Greenwich Mean Time)	Epicenter			Magnitude (Mw)
	Longitude (°)	Latitude (°)	Region
20 September 1999	120.8	24.2	Nantou, Taiwan, China	7.6
14 November 2001	92.9	35.8	Kokoxili, Qinghai, China	7.8
26 December 2004	94.3	3.1	Indian Ocean	9.0
8 October 2005	73.5	34.4	Pakistan	7.6
12 May 2008	104.1	31.4	Wenchuan, Sichuan, China	7.9
6 October 2008	90.5	29.7	Dangxiong, Xizang, China	6.3
11 March 2011	143.1	37.5	Tohoku-Oki, Japan	9.1
20 April 2013	103.1	30.2	Lushan, Sichuan, China	6.6
21 July 2013	104.3	34.6	Minxian, Gansu, China	6.0
12 February 2014	82.6	36.2	Yutian, Xinjiang, China	6.9
25 April 2015	85.3	27.9	Nepal	7.9
12 May 2015	86.1	27.7	Nepal	7.2

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

2.2. 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:

• 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\times n$ residual data matrix $\mathrm{X}$, in which $m$ and $n$ are the number of epochs and stations, respectively.
• Calculate the covariance matrix $\mathrm{B}={\mathrm{X}}^{\mathrm{T}}\mathrm{X}/\left(m-1\right)$.
• Decompose the symmetric matrix $\mathrm{B}$ as ${\mathrm{B}=\mathrm{V}\mathsf{\Lambda }\mathrm{V}}^{\mathrm{T}}$ and sort the eigenvectors to rank the eigenvalues in descending order.
• Consider a linear transformation $\mathrm{A}=\mathrm{XV}$, thus we have $\mathrm{X}={\mathrm{AV}}^{\mathrm{T}}$. The columns of $\mathrm{A}$ and rows of ${\mathrm{V}}^{\mathrm{T}}$ are termed as principle components (PC) and spatial responses (SR), respectively. The $i$th PC and SR are termed as “mode $i$” together.
• Normalize each SR and scaling the corresponding PC by $\mathrm{SR}=\mathrm{SR}/\mathrm{a}$ and $\mathrm{PC}=\mathrm{a}·\mathrm{PC}$, in which $\mathrm{a}$ is the component with the maximum absolute value in this SR.
• Select the significant PCs to calculate the CMC. $\mathrm{CMC}={\mathrm{A}}_j{\mathrm{V}}_j^{\mathrm{T}}$. The columns of ${\mathrm{A}}_j$ are the first $j$ PCs that are considered to be significant and rows of ${\mathrm{V}}_j^{\mathrm{T}}$ 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 Figure A3 and Figure 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 A5. These stations, Mulei (XJML, 43.8°N, 90.3°E), Hebi (HAHB, 35.7°N, 114.5°E) and Cangxian (HECX, 38.5°N, 116.9°E), are characterized by the maximum, medium, and minimum vertical RMS reductions in CMONOC-II after the filtering, respectively.

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.

3. Results

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

3.2. 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 T₁ and T₂, respectively, and T₁ > T₂) are to be distinguishable, the time span of the time series must be no shorter than T₁T₂/(T₁−T₂), 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, $\kappa =-1$) at low frequencies and close to 0 (WN, $\kappa =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.

3.3. 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. Mean and standard deviations of the WN amplitude (${\mathsf{\sigma }}_{WN}$), PLN amplitude (${\mathsf{\sigma }}_{PL}$), spectral index, and velocity uncertainty (${\mathsf{\sigma }}_{vel.}$) and the RMS of the unfiltered and filtered CMONOC-I time series estimated with the PLN + WN model.

	Unfiltered Time Series			Filtered Time Series
	N	E	U	N	E	U
${\mathsf{\sigma }}_{WN}$ (mm)	1.0 ± 0.3	1.1 ± 0.5	2.0 ± 2.0	0.9 ± 0.3	1.0 ± 0.5	2.7 ± 1.7
${\mathsf{\sigma }}_{PL}$ (mm/year−κ/4)	4.2 ± 0.9	4.1 ± 1.1	14.9 ± 3.1	3.1 ± 1.2	3.4 ± 1.5	11.1 ± 4.2
$\kappa$	−1.2 ± 0.1	−1.3 ± 0.2	−0.9 ± 0.1	−1.3 ± 0.2	−1.3 ± 0.2	−1.0 ± 0.2
${\mathsf{\sigma }}_{vel.}$ (mm/year)	0.2 ± 0.1	0.2 ± 0.1	0.4 ± 0.2	0.2 ± 0.1	0.2 ± 0.2	0.3 ± 0.2
RMS (mm)	2.4 ± 1.2	2.6 ± 1.3	6.6 ± 2.3	2.1 ± 1.3	2.3 ± 1.6	5.7 ± 2.7

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. Same as Table 1, but for the CMONOC-II time series.

	Unfiltered Time Series			Filtered Time Series
	N	E	U	N	E	U
${\mathsf{\sigma }}_{WN}$ (mm)	0.7 ± 0.3	0.9 ± 0.4	1.3 ± 1.4	0.6 ± 0.3	0.7 ± 0.4	1.9 ± 1.2
${\mathsf{\sigma }}_{PLN}$ (mm/year−κ/4)	3.4 ± 0.7	3.3 ± 0.8	12.5 ± 1.6	1.9 ± 0.9	1.9 ± 1.0	6.1 ± 2.0
$\kappa$	−0.9 ± 0.2	−1.0 ± 0.2	−0.9 ± 0.1	−1.0 ± 0.4	−1.0 ± 0.4	−0.7 ± 0.3
${\mathsf{\sigma }}_{vel.}$ (mm/year)	0.3 ± 0.2	0.3 ± 0.2	0.9 ± 0.3	0.2 ± 0.3	0.2 ± 0.2	0.4 ± 0.3
RMS (mm)	1.5 ± 0.4	1.6 ± 0.5	5.0 ± 0.9	1.1 ± 0.6	1.1 ± 0.6	3.6 ± 1.1

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.

4. Discussion

4.1. Characteristics of the Draconitic Harmonics

Figure 7, Figure 8, Figure 9 and Figure 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].

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

Station	Lat. (°)	Lon. (°)	Velocity (mm/year)			Velocity Uncertainty (mm/year)
			N	E	U	N	E	U
AHAQ	117.0	30.6	−11.17	33.03	0.72	0.06	0.07	0.36
AHBB	117.3	32.9	−10.78	32.79	0.18	0.10	0.06	0.40
BJFS	115.9	39.6	−10.69	30.20	1.83	0.24	0.18	0.16
BJGB	117.2	40.7	−10.48	29.78	1.57	0.05	0.05	0.23
BJSH	116.2	40.3	−11.28	29.55	0.82	0.09	0.07	0.15
BJYQ	116.0	40.4	−10.59	30.66	1.14	0.09	0.06	0.19
CHUN	125.4	43.8	−11.77	25.89	0.27	0.12	0.14	0.45
CQCS	107.2	29.9	−9.15	34.68	0.27	0.09	0.17	0.19
DLHA	97.4	37.4	1.03	37.23	0.70	0.15	0.32	0.23
DXIN	100.2	41.0	−4.58	30.72	1.15	0.13	0.10	0.18
FJPT	119.8	25.5	−11.71	31.22	0.40	0.30	0.47	0.28
FJWY	118.0	27.6	−11.50	33.03	0.93	0.16	0.14	0.47
FJXP	120.0	26.9	−11.46	32.96	1.06	0.08	0.15	0.23
GDSG	113.6	24.8	−10.59	33.60	−0.06	0.11	0.44	0.57
GDZH	113.6	22.3	−11.18	33.68	5.75	0.07	0.09	0.40
GDZJ	110.3	21.2	−10.50	33.67	2.06	0.20	0.23	0.40
GSAX	95.8	40.5	−0.24	31.27	0.21	0.12	0.09	0.13
GSDH	94.7	40.1	0.48	30.81	−0.01	0.17	0.26	0.19
GSDX	104.6	35.6	−6.69	36.40	0.68	0.06	0.04	0.19
GSGL	102.9	37.5	−4.78	34.37	1.96	0.45	0.48	1.07
GSGT	99.8	39.4	−3.38	31.61	1.24	0.09	0.08	0.18
GSJN	105.8	35.5	−8.49	34.68	0.67	0.18	0.08	0.44
GSJT	104.1	37.2	−5.29	34.75	1.09	0.08	0.25	0.13
GSJY	98.2	39.8	−1.93	31.51	2.02	0.17	0.16	0.19
GSLX	104.6	35.0	−7.01	35.52	0.91	0.10	0.11	0.26
GSLZ	103.7	36.1	−0.75	40.94	−2.66	1.75	0.75	0.23
GSMA	102.1	34.0	−4.25	41.41	1.31	0.35	0.13	0.27
GSML	100.8	38.4	−2.92	33.02	0.68	0.11	0.12	0.19
GSMQ	103.1	38.6	−5.50	32.48	−0.30	0.37	0.06	0.71
GSMX	104.0	34.4	−5.92	38.08	0.83	0.12	0.14	0.36
GSPL	106.6	35.5	−8.95	35.06	1.04	0.05	0.61	0.48
GSQS	106.2	34.7	−7.71	33.77	0.34	0.22	0.20	0.37
GSWD	104.8	33.4	−8.15	36.17	−0.24	0.39	0.31	0.50
GUAN	113.3	23.2	−11.03	32.64	0.55	0.20	0.32	0.97
GXBH	109.2	21.7	−8.41	33.30	2.20	0.31	0.14	0.45
GXBS	106.7	23.9	−9.07	34.60	−1.07	0.16	0.12	0.65
GXHC	108.0	24.7	−9.30	34.54	1.07	0.08	0.11	0.36
GXNN	108.1	22.6	−9.22	33.85	2.07	0.14	0.09	0.55
GXWZ	111.2	23.5	−9.43	34.24	−1.57	0.16	0.16	0.47
GZFG	107.7	28.0	−9.42	33.32	1.89	0.37	0.60	0.82
GZGY	106.7	26.5	−9.01	34.94	0.37	0.05	0.20	0.70
GZSC	104.9	26.6	−8.73	34.68	1.46	0.08	0.10	0.31
HAHB	114.5	35.7	−10.81	32.43	−0.54	0.06	0.09	0.42
HAJY	112.4	35.2	−10.35	33.14	0.19	0.15	0.09	0.56
HAQS	114.0	32.8	−10.61	34.24	−1.64	0.04	0.08	0.40
HBES	109.5	30.3	−9.99	34.10	1.94	0.19	0.11	0.81
HBJM	112.2	31.1	−10.73	34.77	0.83	0.19	0.34	0.38
HBXF	112.0	32.0	−10.48	33.69	1.20	0.10	0.06	0.58
HBZG	111.0	30.8	−9.26	34.57	0.10	0.47	0.44	0.49
HECC	115.8	40.9	−10.64	29.85	1.28	0.10	0.05	0.21
HECD	117.9	41.0	−9.85	29.81	0.71	0.08	0.12	0.36
HECX	116.9	38.5	−9.74	30.82	−24.12	0.43	0.38	2.34
HELQ	114.3	38.2	−10.48	32.11	4.24	0.08	0.14	0.49
HELY	114.7	37.4	−10.37	33.23	1.17	0.20	0.14	0.25
HETS	118.3	39.7	−6.33	27.70	2.39	0.77	0.89	1.43
HEYY	114.2	40.1	−10.90	32.47	1.01	0.12	0.39	0.31
HEZJ	114.9	40.8	−9.91	30.08	2.05	0.34	0.19	0.29
HIHK	110.2	20.0	−10.43	33.02	−1.25	0.37	0.13	0.57
HISY	109.5	18.2	−8.15	34.07	2.68	0.26	0.63	0.75
HLAR	119.7	49.3	−11.32	25.79	1.44	0.10	0.09	0.30
HLHG	130.2	47.4	−15.90	24.38	0.29	0.78	1.45	0.54
HLWD	126.1	48.7	−12.24	25.21	0.25	0.44	0.41	0.83
HNLY	113.6	28.2	−10.35	33.72	0.03	0.08	0.08	0.59
HNMY	109.8	27.9	−9.71	34.51	0.68	0.06	0.20	0.29
HRBN	126.6	45.7	−11.76	24.98	−0.45	0.16	0.18	0.21
JIXN	117.5	40.1	−10.12	28.17	1.74	0.12	0.06	0.17
JLCL	123.5	44.6	−11.62	28.38	−0.55	0.26	0.23	0.35
JLYJ	129.5	42.9	−16.81	26.47	−0.11	0.32	0.70	0.44
JSLS	119.4	31.3	−11.39	33.56	−1.92	0.27	0.12	1.14
JSLY	119.5	34.7	−12.05	32.06	0.47	0.06	0.08	0.26
JSNT	120.9	32.0	−11.78	32.83	−1.33	0.10	0.13	0.31
JSYC	120.0	33.4	−9.76	32.46	−11.28	0.44	0.27	1.46
JXJA	115.1	26.7	−11.13	33.93	−0.44	0.06	0.31	0.37
KMIN	102.8	25.0	−16.73	33.53	−1.67	0.39	0.97	0.75
LALB	102.2	19.9	−6.25	31.40	2.23	0.40	0.28	0.34
LALX	105.0	18.2	−8.74	31.52	2.68	0.40	0.44	0.47
LHAS	91.1	29.7	15.50	46.66	1.79	0.12	0.15	0.38
LNDD	124.3	40.0	−10.68	26.09	0.13	0.27	1.12	0.37
LNJZ	121.7	39.1	−12.03	29.42	0.63	0.29	0.11	0.42
LNYK	122.6	40.7	−11.64	27.76	1.20	0.12	0.26	0.31
LUZH	105.4	28.9	−9.56	34.65	0.73	0.10	0.19	0.12
NMAG	122.6	43.3	−11.00	28.74	1.67	0.10	0.15	0.33
NMAL	120.1	43.9	−10.92	28.74	1.20	0.12	0.10	0.36
NMAY	101.7	39.2	−5.25	31.82	1.51	0.23	0.05	0.12
NMAZ	105.7	38.8	−6.63	31.84	0.92	0.26	0.21	0.33
NMDW	117.0	45.5	−10.53	28.40	−0.24	0.12	0.08	0.29
NMEJ	101.1	42.0	−4.50	31.41	1.12	0.23	0.03	0.26
NMEL	111.9	43.6	−9.61	29.52	1.70	0.17	0.30	0.17
NMER	123.7	50.6	−11.28	25.12	0.71	1.92	1.37	0.50
NMTK	111.3	40.3	−10.57	31.58	1.61	0.18	0.40	0.41
NMWH	106.8	39.7	−9.62	32.42	3.75	0.12	0.19	0.29
NMWJ	108.1	41.3	−9.43	30.09	0.86	0.12	0.27	0.29
NMWL	122.0	46.0	−10.96	27.57	1.12	0.17	0.17	0.56
NMZL	116.0	42.2	−10.20	29.60	1.33	0.28	0.03	0.20
NXHY	105.6	36.6	−6.55	34.59	−0.93	0.10	0.15	0.18
NXYC	106.3	38.5	−7.34	32.04	−4.25	0.68	0.62	0.49
NXZW	105.2	37.6	−5.62	33.25	1.52	0.05	0.13	0.16
QHBM	100.7	32.9	−4.94	44.90	1.02	0.07	0.05	0.18
QHDL	98.1	36.3	2.53	38.02	0.58	0.08	0.11	0.27
QHGC	100.1	37.3	−2.10	37.58	0.63	0.14	0.07	0.34
QHGE	94.8	36.1	4.53	37.16	0.20	0.14	0.13	0.65
QHLH	93.3	38.7	3.12	34.62	2.86	0.50	0.15	0.38
QHMD	98.2	34.9	0.50	45.52	1.62	0.07	0.09	0.42
QHME	101.4	37.5	−3.12	37.09	1.35	0.10	0.14	0.40
QHMQ	100.2	34.5	−1.12	42.30	0.51	0.17	0.05	0.22
QHMY	90.8	38.5	7.16	35.51	1.06	0.14	0.07	0.22
QHQL	100.2	38.2	−3.13	34.70	1.67	0.15	0.34	0.30
QION	109.8	19.0	−10.93	31.64	0.39	0.20	0.15	0.28
SCBZ	106.7	31.8	−9.13	34.33	1.21	0.04	0.06	0.20
SCDF	101.1	31.0	−12.25	43.38	1.60	0.09	0.48	0.33
SCGY	105.9	32.4	−8.33	34.56	0.58	0.19	0.79	0.38
SCGZ	100.0	31.6	−9.30	47.06	−0.60	0.11	0.42	0.33
SCJL	101.5	29.0	−17.82	39.77	2.03	0.10	0.06	0.52
SCJU	104.5	28.2	−7.61	34.73	0.51	0.17	0.27	0.84
SCLH	100.7	31.4	−11.09	46.19	1.16	0.06	0.19	0.28
SCMB	103.5	28.8	−8.99	35.25	0.22	0.13	0.21	0.32
SCML	101.3	27.9	−19.49	48.87	−2.91	0.07	0.24	0.39
SCMN	102.2	28.3	−16.29	39.10	0.32	0.05	0.09	0.20
SCMX	103.8	31.7	−6.48	45.06	9.28	0.24	0.31	0.34
SCNC	105.9	31.0	−9.58	34.82	1.28	0.10	0.19	0.12
SCNN	102.7	27.1	−16.34	37.93	0.67	0.61	0.36	0.27
SCPZ	101.7	26.5	−17.80	36.01	1.06	0.09	0.08	0.15
SCSM	102.4	29.2	−12.19	38.08	1.86	0.05	0.06	0.34
SCSN	105.6	30.5	−8.83	34.30	1.50	0.04	0.11	0.21
SCSP	103.6	32.6	−10.03	41.08	1.43	0.14	0.15	0.34
SCTQ	102.8	30.1	−9.66	35.94	2.92	0.18	0.29	0.45
SCXC	99.8	28.9	−16.99	41.52	2.82	0.05	0.12	0.34
SCXD	102.4	28.3	−14.15	38.71	0.17	0.04	0.09	0.19
SCXJ	102.4	31.0	−6.05	40.86	0.30	0.35	0.24	0.19
SCYX	102.5	28.7	−13.09	38.12	−0.91	0.13	0.16	0.30
SCYY	101.5	27.4	−17.60	38.62	0.80	0.10	0.10	0.16
SDCY	119.5	36.8	−11.46	31.38	1.84	0.11	0.07	0.47
SDJX	116.4	35.4	−11.06	32.74	1.82	0.13	0.08	0.30
SDLY	118.3	35.0	−12.42	31.71	2.72	0.94	0.28	0.70
SDQD	120.3	36.1	−11.34	31.78	1.24	0.05	0.05	0.19
SDRC	122.4	37.2	−11.38	31.25	−2.27	0.10	0.14	0.65
SDYT	121.4	37.5	−10.96	30.29	0.68	0.17	0.07	0.18
SDZB	118.0	36.8	−11.28	32.24	0.92	0.27	0.17	0.42
SHA2	121.2	31.1	−11.77	31.93	−1.55	0.15	0.17	0.48
SHAO	121.2	31.1	−13.32	32.16	−1.91	0.25	0.24	0.27
SNAK	108.8	32.8	−9.81	34.62	1.45	0.25	0.28	0.40
SNMX	106.7	33.1	−10.57	35.06	−1.39	0.68	0.28	0.82
SNTB	107.3	34.1	−9.60	35.28	1.06	0.06	0.10	0.43
SNXY	108.4	35.2	−9.27	34.18	1.28	0.06	0.06	0.37
SNYA	109.5	36.6	−10.47	32.08	2.16	0.12	0.59	0.27
SUIY	130.9	44.4	−12.49	24.52	−1.14	0.11	0.18	0.40
SXCZ	113.2	36.2	−9.58	30.04	1.34	0.25	0.29	0.23
SXGX	111.9	36.3	−9.69	32.54	3.02	0.36	0.10	0.61
SXKL	111.6	38.8	−8.84	32.73	2.97	0.12	0.06	0.71
SXLF	111.4	36.1	−10.14	33.63	1.77	0.08	0.09	0.39
SXLQ	114.0	39.4	−10.32	32.39	1.87	0.05	0.05	0.25
SXTY	112.4	37.7	−9.21	33.87	−0.28	0.47	0.58	0.51
SXXX	111.2	35.1	−8.69	33.03	1.38	0.44	0.14	0.47
SXYC	112.9	37.6	−10.62	32.55	2.26	0.09	0.09	0.55
TAIN	117.1	36.2	−11.73	30.76	0.92	0.10	0.05	0.27
TASH	75.2	37.8	23.46	25.44	1.34	0.15	0.15	0.35
TJBD	117.4	39.7	−11.65	30.62	1.11	0.20	0.04	0.29
TJBH	117.7	39.1	−11.92	31.96	−16.82	0.49	0.24	0.28
URU2	87.6	43.8	5.17	30.83	1.41	0.18	0.23	0.79
URUM	87.6	43.8	6.89	32.52	1.63	0.34	0.36	0.43
WUHN	114.4	30.5	−10.73	32.28	0.03	0.28	0.20	0.51
WUSH	79.2	41.2	15.17	30.30	1.16	0.12	0.17	0.21
XIAA	109.0	34.2	−6.77	30.06	−5.10	0.30	0.47	0.31
XIAG	100.3	25.6	−16.68	29.96	0.67	0.14	0.25	0.33
XIAM	118.1	24.4	−12.82	32.57	0.61	0.14	0.12	0.26
XJAL	88.1	47.9	3.72	28.83	0.31	0.37	0.30	1.16
XJBC	78.8	39.8	17.37	29.21	0.49	0.11	0.10	0.21
XJBE	86.9	47.7	3.34	29.20	−0.14	0.09	0.12	0.25
XJBL	75.0	38.7	21.97	24.10	−0.41	0.25	0.35	0.65
XJBY	83.7	42.8	7.75	31.59	1.08	0.15	0.26	0.33
XJDS	84.9	44.3	6.68	31.88	0.72	0.08	0.09	0.20
XJFY	89.5	47.0	3.15	29.58	0.20	0.24	0.10	0.35
XJHT	79.0	37.2	18.38	26.28	0.31	0.22	0.12	0.24
XJJJ	94.3	42.8	0.54	32.55	1.45	0.07	0.05	0.17
XJKC	83.0	41.7	10.92	32.23	−1.12	0.09	0.10	0.24
XJKE	86.2	41.8	8.06	31.22	−0.40	0.17	0.28	0.53
XJML	90.3	43.8	3.39	32.19	1.13	0.08	0.14	0.16
XJQH	91.0	46.2	2.54	31.35	0.51	0.11	0.13	0.15
XJQM	85.5	38.1	8.63	28.83	−0.10	0.32	0.19	0.20
XJRQ	88.2	39.0	6.62	29.97	−0.67	0.06	0.04	0.11
XJSH	86.1	44.2	5.00	31.66	1.49	0.07	0.07	0.22
XJSS	90.3	42.9	3.84	32.27	−0.27	0.07	0.06	0.12
XJTC	82.9	46.8	2.38	28.37	0.00	0.15	0.21	0.39
XJTZ	83.7	39.0	11.41	29.25	0.54	0.31	0.04	0.18
XJWL	86.7	42.9	6.84	32.33	1.10	0.17	0.05	0.21
XJWQ	81.0	45.0	4.22	29.32	1.22	0.14	0.06	0.22
XJWU	75.2	39.7	14.86	28.82	−0.20	0.18	0.26	0.37
XJXY	83.3	43.4	7.17	30.65	1.99	0.05	0.05	0.40
XJYC	77.4	37.9	20.92	27.25	−0.80	0.18	0.11	0.29
XJYT	82.0	36.4	13.85	25.58	−0.42	0.30	0.20	0.34
XJZS	80.4	42.7	7.99	31.48	1.70	0.09	0.16	0.14
XNIN	101.8	36.6	−3.76	37.88	−0.36	0.17	0.13	0.23
XZAR	87.2	29.3	21.06	40.03	3.06	0.19	0.39	0.33
XZCD	97.2	31.1	−4.74	50.20	1.93	0.11	0.15	0.71
XZCY	97.5	28.7	−9.34	42.22	4.20	0.17	0.18	0.52
XZDX	91.1	30.5	13.77	47.03	2.31	0.29	0.10	0.59
XZGE	80.1	32.5	17.96	30.82	0.99	0.08	0.10	0.30
XZNQ	92.1	31.5	9.10	50.73	0.86	0.29	0.18	0.70
XZRK	88.9	29.2	21.54	42.39	2.27	0.18	0.16	0.32
XZRT	79.7	33.4	17.72	30.21	1.48	0.28	0.21	0.76
XZYD	88.9	27.4	29.84	39.90	−0.26	0.56	0.34	0.90
XZZB	84.2	29.7	22.35	36.60	2.00	0.60	0.26	0.68
XZZF	86.9	28.4	24.02	40.01	4.74	0.31	0.22	0.98
YANC	107.4	37.8	−9.43	31.99	1.60	0.11	0.09	0.21
YNCX	101.5	25.0	−16.38	31.63	1.40	0.06	0.08	0.16
YNDC	103.2	26.1	−12.33	34.99	1.82	0.12	0.09	0.50
YNHZ	103.3	26.4	−11.45	34.91	3.49	0.29	0.21	0.64
YNJD	100.9	24.4	−16.11	29.06	1.91	0.54	0.49	0.43
YNJP	103.2	22.8	−8.42	34.33	0.87	0.35	0.15	0.46
YNLA	100.0	22.6	−9.92	28.61	−0.78	0.41	0.70	0.83
YNLC	100.1	23.9	−13.21	29.85	0.34	0.17	0.71	0.22
YNLJ	100.0	26.7	−18.19	32.54	1.16	0.07	0.06	0.13
YNMH	100.5	21.9	−9.93	28.26	1.05	0.76	1.01	0.57
YNMJ	101.7	23.4	−14.35	31.01	1.98	0.39	0.25	0.30
YNML	103.4	24.4	−9.85	35.53	1.15	2.33	0.76	0.44
YNMZ	103.4	23.4	−9.13	33.72	−0.25	0.34	0.27	0.32
YNRL	97.8	24.0	−7.87	23.50	1.10	0.26	0.30	0.64
YNSD	99.2	24.7	−13.70	26.88	1.85	0.09	0.08	0.33
YNSM	101.0	22.7	−13.65	29.20	−1.84	2.73	1.00	0.43
YNTC	98.4	25.0	−10.34	24.65	1.95	0.09	0.43	0.31
YNTH	102.8	24.1	−13.11	33.04	2.40	0.13	0.19	0.21
YNWS	104.2	23.4	−10.57	32.14	3.82	0.87	0.89	0.18
YNYA	101.3	25.7	−16.45	33.30	2.15	0.04	0.10	0.15
YNYM	101.9	25.7	−16.97	33.83	1.72	0.27	0.13	0.18
YNYS	100.8	26.7	−16.97	35.98	1.95	0.08	0.17	0.17
YNZD	99.7	27.8	−20.66	34.02	−1.76	0.13	0.18	0.38
YONG	112.3	16.8	−10.80	29.73	0.76	0.40	0.40	0.95
ZHNZ	113.1	34.5	−11.40	32.76	2.42	0.06	0.13	0.33
ZJJD	119.3	29.5	−11.71	33.07	0.56	0.07	0.10	0.75
ZJWZ	120.8	27.9	−11.67	32.37	0.10	0.07	0.13	0.38
ZJZS	122.0	30.1	−11.29	32.32	0.05	0.11	0.11	0.27

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 Figure 14 and Figure 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.

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

5. 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:

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