A Sub-Regional Extraction Method of Common Mode Components from IGS and CMONOC Stations in China

: There is always a need to extract more accurate regional common mode component (CMC) series from coordinate time series of Global Positioning System (GPS) stations, which would be of great beneﬁt to describe the deformation features of the Earth’s surface with more reliability. For this purpose, this paper combines all 11 International Global Navigation Satellite System (GNSS) Service (IGS) stations in China with over 70 stations selected from the Crustal Movement Observation Network of China (CMONOC) to compute CMC series of IGS stations by using a principal component analysis (PCA) method under cases of one whole region and eight sub-regions. The comparison results show that the percentage of ﬁrst-order principal component (PC1) in North, East and Up components increase by 10.8%, 16.1% and 25.1%, respectively, after dividing the whole China region into eight sub-regions. Meanwhile, Root Mean Square (RMS) reduction rates of residual series that have removed CMC also improve obviously after partitioning. In addition, we compute displacements of these IGS stations caused by environmental loadings (including atmospheric pressure loading, non-tidal oceanic loading and hydrological loading) to analyze their contributions to the non-linear variation in GPS coordinate time series. The comparison result shows that the method we raise, PCA ﬁltering in sub-regions, performs better than the environmental loading corrections (ELCs) in improving the signal-to-noise ratio (SNR) of GPS coordinate time series. This paper raises new criteria for selecting appropriate CMONOC stations around IGS stations when computing sub-regional CMC, involving three criteria of interstation distance, geology and self-condition of stations themselves. According to experiments, these criteria are implemental and e ﬀ ective in selecting suitable stations, by which to extract sub-regional CMC with higher accuracy.


Introduction
The continuously operating GNSS reference stations established by the International GNSS Service (IGS) have accumulated a huge amount of GNSS data for over 20 years. The results from these data, coordinate time series of IGS stations, are constantly getting more accurate, accompanied by the gradual improvement of GNSS technology, data-processing methods and error correction models. Topics, such as the global or regional observation of crust deformation and regional fault slips have significant scientific value and practical meaning [1][2][3][4][5][6]. However, the coordinate time series of GPS stations have various influential factors, including one related to time and space is called common mode error (CME) [7,8]. CME series have obvious periodic terms that are probably caused by some geophysical signals [9,10]. Therefore, we use the term "common mode component" (CMC) in this paper instead.
Effective extraction of CMC is quite beneficial in improving the accuracy of coordinate time series, quantifying the station displacement caused by geophysical signals and accurately describing the features of regional earth surface deformation. As a result, recent years have seen CMC become a heated research topic [9][10][11][12][13][14][15][16][17][18][19]. Jiang et al. [20] analyzed 29 stations in the Chinese CORS network and extracted the CMC series by using a PCA method. After filtering CMC from the original coordinate time series, they estimated the velocity and its accuracy of each station. Tian and Shen [15] put forward a correlation-weighted stacking filtering method to remove CMC. This method takes the correlation coefficients of the coordinate components (North, East and Up) among the stations as the weighs in filtering, then the root mean square (RMS) value reduced effectively. Besides, it is available in both large-scale and small-scale GPS networks. Ma et al. [16] employed three kinds of spatial filtering approaches: stacking, PCA and Karhunen-Loeve Expansion (KLE) to extract CMC from stations in the Antarctic Peninsula. The results showed that PCA is an efficient filtering method with the highest accuracy. Zhu et al. [18] calculated the CMC of CMONOC stations through correlation weighted stacking filtering method. Then they quantified the geophysical signals that induced CMC, including environmental loading and thermal expansion. The results announced that environmental loading is one of the main sources of CMC in the vertical component of stations, while the impact of thermal expansion is not as obvious as the former.
Through continuous exploration in recent years, theoretical CMC research has matured, whereas there are still questions worth to discussing further. For example, when we compute the CMC of all the stations in a whole region, the CMC of the sub-regions would be easily regarded as the local individual feature and hence be neglected. This would then lead to decreasing CMC extraction accuracy and inaccurate cognition about the motion features of GPS stations. In addition, according to the existing research results [11,18,21,22], the scale of the research region has direct influence to extracting CMC, so what are the specific characteristics of such influence in China region? Does it have some useful and feasible instructions for the regional CMC extraction mode that this paper proposes? We will explore these topics next.

GPS Coordinate Time Series
The data that this paper uses to calculate the CMC of IGS stations in China region is mainly from Scripps Orbit and Permanent Array Center (SOPAC), one of IGS data analysis centers, and from China Earthquake Administration (CEA) that provides CMONOC data. Table 1 shows the related information about their products that we use in this paper. SOPAC products are from the address "ftp://garner.ucsd.edu/pub/timeseries/measures" which contains "raw" series that include outliers, "cleaned" series that has excluded outliers and "detrend" series that has removed linear trend. CMONOC products are acquired from "ftp://ftp.cgps.ac.cn/products/position/gamit/" which includes coordinate time series of 260 CMONOC stations that spread evenly throughout China. They are also divided into "raw" and "detrend" series. The daily loosely constraint solutions of these CMONOC stations are processed by the software GAMIT10.4. The detailed processing strategy emerges in the manual whose address is "ftp://ftp.cgps. ac.cn/doc/processing_manual.pdf". Figure 1 shows all the IGS and CMONOC stations in China. The data that this paper uses is the "cleaned" series of 11 IGS stations and "detrend" series of over 70 selected CMONOC stations.

Environmental Loading Data
The mass redistribution of atmosphere, ocean and terrestrial water reserves can cause changes in the Earth's surface mass loading, thus leading to surface deformation [23][24][25]. This is called environmental loading deformation. Surface loading models (SLM) can be established based on geophysical observations and Earth models to estimate surface displacements caused by environmental loadings.
This paper estimates the displacement of IGS stations caused by environmental loading based on the surface loading products provided by GFZ [26]. GFZ provides global surface elastic deformation data in grids form, including atmospheric pressure loading (NTAL), non-tidal oceanic loading (NTOL) and hydrological loading (HYDL), with a spatial resolution of 0.5° × 0.5°. The temporal resolution is 3 hours for NTAL and NTOL data, and 24 hours for HYDL data. The input data GFZ calculates NTAL, NTOL and HYDL are surface pressure data from European Centre for Medium-Range Weather Forecasts (ECMWF), Max Planck Institute Ocean Model (MPIOM, [27]) and land surface discharge model (LSDM, [28]). Based on the environmental loading products (netCDF files) of GFZ, we adopt a bicubic interpolation to interpolate certain station locations and then to compute the environmental loading series of each CMONOC station. In order to unify the temporal

Environmental Loading Data
The mass redistribution of atmosphere, ocean and terrestrial water reserves can cause changes in the Earth's surface mass loading, thus leading to surface deformation [23][24][25]. This is called environmental loading deformation. Surface loading models (SLM) can be established based on geophysical observations and Earth models to estimate surface displacements caused by environmental loadings.
This paper estimates the displacement of IGS stations caused by environmental loading based on the surface loading products provided by GFZ [26]. GFZ provides global surface elastic deformation data in grids form, including atmospheric pressure loading (NTAL), non-tidal oceanic loading (NTOL) and hydrological loading (HYDL), with a spatial resolution of 0.5 • × 0.5 • . The temporal resolution is 3 h for NTAL and NTOL data, and 24 h for HYDL data. The input data GFZ calculates NTAL, NTOL and HYDL are surface pressure data from European Centre for Medium-Range Weather Forecasts  [27]) and land surface discharge model (LSDM, [28]). Based on the environmental loading products (netCDF files) of GFZ, we adopt a bicubic interpolation to interpolate certain station locations and then to compute the environmental loading series of each CMONOC station. In order to unify the temporal resolution of these environmental loading displacement series, NTAL and NTOL displacement series take the average every six epochs to get the daily time series. Then we add these three loadings together to acquire the sum environmental loading (SUML) of 11 IGS stations in China. The surface elastic deformation data used takes the Center of Earth's Figure (CF) as the reference frame [29].

Methodology
At present, the popular method for CMC extraction is spatiotemporal filtering, which includes regional stacking [7,21,30], PCA [10,31,32] and KLE [9]. The CMC series is contained in the residual time series of stations, so when we get the original coordinate time series, the first thing is always to acquire the residual series by removing the linear and periodic terms.

Acquisition of the Residual Time Series
The original coordinate time series of GPS stations consists of a linear trend, periodic terms (annual and semi-annual), offsets and noises, so the mathematical model that reflects the displacement at coordinate components (North, East and Up) can be expressed as follows [33,34]: where x R and t R are the reference coordinate and epoch; t are epochs of GPS time series; v means linear velocity; H(t) is the Heaviside step function and used to model offsets with amplitudes b j ; n b is the number of offsets and t j are epochs of these offsets; A k , w k , ϕ k are the amplitudes, angular velocities and initial phases of periodic signals, respectively; n A is the number of periodic signals, k = 1 for annual signal and k = 2 for semi-annual signal; n L is the number of logarithmic functions used to model post-seismic deformation; a k and T k are two parameters of the logarithmic function; t k are the epochs of earthquakes that caused logarithmic displacement; v(t) is the residual series that represents the difference between the observations and the fitting values. It is also used for extracting CMC series. For offsets contained in coordinate time series, firstly we remove some of them (including earthquakes and antenna/receiver changing) based on the epochs defined by CEA. Then we manually inspect the remaining offsets in stations one by one, as it is a better method to distinguish them [35]. If some large earthquakes near a station happened right in the same epochs, we assume that these offsets are probably caused by earthquakes and remove them in the GPS time series. Besides, a few station components have obvious post-seismic relaxations, and we use logarithmic functions to fit them. We use least squares fitting (LSF) method to calculate the linear trend, periodic terms and post-seismic relaxations in Equation (1), then we obtain the residual series of each station component. For outliers remaining in residual series, three-time interquartile range (3IQR) method is used for culling. After acquiring the residual series, then we extract CMC from them.

Principal Component Analysis
PCA is a popular approach to reduce the dimensionality of high-dimensional data. The principle is to transform the original data into a set of linearly independent components and then obtain the principal components (PCs) that represent the biggest contributors to the original data. In the application of GPS coordinate time series, PCA decomposes the residual series of GPS stations into PCs in the temporal domain and the corresponding eigenvectors in the spatial domain. Then the spatiotemporal variation of the GPS network can be expressed to the largest extent through the first several PCs. For a regional GPS network with n stations that spanning m days, their residual coordinate time series compose a matrix: B is the corresponding covariance matrix with elements that are defined as: The symmetric matrix B can be decomposed as: where V T (n × n), formed by eigenvectors, is a matrix with orthonormal rows; Λ is a diagonal matrix composed by non-zero eigenvalues of matrix B; generally in real GPS data, matrix B is full in rank, which means r = n. The eigenvalues and eigenvectors of matrix B can be written as , v n is a set of orthonormal bases that help to expand X(t i , x j ): where a k (t) is the k-th PC and can be derived by: where v k (x) is the eigenvector of matrix B. PCs represent the temporal variations and the eigenvectors represent the spatial responses to the PCs. The decomposition result of PCA is arranged in descending order, which means the first-order principal component (PC1) always retains the largest amount of information about the regional GPS network. PC1 reflects the common temporal characteristics of the whole network, while the high-order PCs are usually related to local or individual site effects. The eigenvectors of these PCs, however, reflect the spatial distribution of the corresponding temporal change. As a result, we calculate CMC by: where p refers to the number of PCs that remains to be ascertained in real data processing. PCA has the advantage of simplifying questions and making them intuitive because it converts a question from high dimensions to lower dimensions. Besides, it obtains the features of the CMC spatiotemporal distribution completely through the original data itself, unlike the regional stacking filtering method that firstly supposes that the spatial distribution of CMC is even. Due to these merits, we employ PCA to compute the CMC in the research region.

CMC Extraction
In this section, we use PCA to extract CMC of IGS stations in two modes. The first extraction mode is one whole extraction from all IGS stations in China, a direct and most popular method. Another extraction mode that this paper proposes is eight sub-regional extractions, taking near CMONOC stations and tectonic units into consideration. Our purpose is to gain a CMC series with more accuracy.  Figure 1 shows all eleven IGS stations in China. The longest baseline is TWTF to GUAO, with 3654 km while the shortest is TNML to TCMS, only 8 m. The average length of these baselines is 2016 km, so the GPS network is medium-scale [22]. Next, we use PCA to extract CMC from all IGS stations in China, in order to analyze the amount of CMC in this medium-scale network. Figure 2 illustrates the effect of different time spans of the coordinate time series on the percentage of PC1 when uniformly extracting the CMC of the IGS station in China. In general, the length of the coordinate time series is inversely proportional to the percentage of PC1. The longer the time span is, the smaller the percentage of PC1. This can be understood as the longer the coordinate time series accumulated by each station, the less significant the common components in the residual series. Specifically, when the time series span is less than 10 years, the percentage of PC1 varies greatly with different time spans. When the time span is longer than 10 years, it gradually becomes stable. is 2016 km, so the GPS network is medium-scale [22]. Next, we use PCA to extract CMC from all IGS stations in China, in order to analyze the amount of CMC in this medium-scale network. Figure 2 illustrates the effect of different time spans of the coordinate time series on the percentage of PC1 when uniformly extracting the CMC of the IGS station in China. In general, the length of the coordinate time series is inversely proportional to the percentage of PC1. The longer the time span is, the smaller the percentage of PC1. This can be understood as the longer the coordinate time series accumulated by each station, the less significant the common components in the residual series. Specifically, when the time series span is less than 10 years, the percentage of PC1 varies greatly with different time spans. When the time span is longer than 10 years, it gradually becomes stable. The percentages of PC1 in North and East components are about 34%, and about 26% in Up component. For the coordinate time series of a nine years' span (2010-2019), Figure 3 shows the percentage of all the 11 PCs in North, East and Up component, which is also the contribution rate of each PC to the residual time series. Table 2 gives information about normalized eigenvectors of top three PCs. They are unitless values calculated through the eigenvalue of one station dividing the largest eigenvalue among these stations.  For identifying the number of PC used for computing CMC, we can choose it both from eigenvalues and eigenvectors. According to the treatment of Dong et al. [9], we regard the component as PC if most sites (more than 50%) have significant normalized responses (larger than 0.25). From Figure 3, we know that the contribution percentage of PC1 to residual series in North, East and Up is 43.8%, 38.8% and 29.6% respectively. Table 2 shows that the normalized eigenvectors of PC1 are over For the coordinate time series of a nine years' span (2010-2019), Figure 3 shows the percentage of all the 11 PCs in North, East and Up component, which is also the contribution rate of each PC to the residual time series. Table 2 gives information about normalized eigenvectors of top three PCs. They are unitless values calculated through the eigenvalue of one station dividing the largest eigenvalue among these stations. is 2016 km, so the GPS network is medium-scale [22]. Next, we use PCA to extract CMC from all IGS stations in China, in order to analyze the amount of CMC in this medium-scale network. Figure 2 illustrates the effect of different time spans of the coordinate time series on the percentage of PC1 when uniformly extracting the CMC of the IGS station in China. In general, the length of the coordinate time series is inversely proportional to the percentage of PC1. The longer the time span is, the smaller the percentage of PC1. This can be understood as the longer the coordinate time series accumulated by each station, the less significant the common components in the residual series. Specifically, when the time series span is less than 10 years, the percentage of PC1 varies greatly with different time spans. When the time span is longer than 10 years, it gradually becomes stable. The percentages of PC1 in North and East components are about 34%, and about 26% in Up component. For the coordinate time series of a nine years' span (2010-2019), Figure 3 shows the percentage of all the 11 PCs in North, East and Up component, which is also the contribution rate of each PC to the residual time series. Table 2 gives information about normalized eigenvectors of top three PCs. They are unitless values calculated through the eigenvalue of one station dividing the largest eigenvalue among these stations.  For identifying the number of PC used for computing CMC, we can choose it both from eigenvalues and eigenvectors. According to the treatment of Dong et al. [9], we regard the component as PC if most sites (more than 50%) have significant normalized responses (larger than 0.25). From Figure 3, we know that the contribution percentage of PC1 to residual series in North, East and Up is 43.8%, 38.8% and 29.6% respectively. Table 2 shows that the normalized eigenvectors of PC1 are over  For identifying the number of PC used for computing CMC, we can choose it both from eigenvalues and eigenvectors. According to the treatment of Dong et al. [9], we regard the component as PC if most sites (more than 50%) have significant normalized responses (larger than 0.25). From Figure 3, we know that the contribution percentage of PC1 to residual series in North, East and Up is 43.8%, 38.8% and 29.6% respectively. Table 2 shows that the normalized eigenvectors of PC1 are over 0.25 on each component of more than 50% sites, while the normalized eigenvectors of PC2 meet the above demand in North component but fail in East and Up. Consequently, the CMC on North is PC1 and PC2, while the CMC on East and Up components is PC1 only.

One Extraction from all IGS Stations
Considering the fact that the contribution percentage of PC1 that represents regional common deformation is small, the CMC of all IGS stations in China region probably fails to be extracted accurately. From the above equations, CMC extraction is influenced by the number of GPS stations and distances among them. Therefore, it is highly necessary to combine some CMONOC stations to extract CMC of these IGS stations in higher accuracy. When deciding the available CMONOC stations, this paper takes two factors into account, the distance away from an IGS station and division of tectonic units. As the region size has a direct influence on CMC extraction, we should ascertain the appropriate region size and study the features of such influence in China.

Influence of Region Size on CMC Extraction
When exploring the influence of region size on CMC extraction, we choose one satisfactory experimental region. Given the fact that Northeast Block of China is the most stable among all blocks in China, without an obvious interior fault zones [36], it proves to be the most suitable research region. Moreover, most CMONOC stations in Northeast Block have similar monument types and geological features, which helps to approximately regard the distance as the unique variable that influences the accuracy of CMC extraction.
The spatial distribution of IGS and CMONOC stations in Northeast Block of China is shown in Figure 4. When carrying out this experiment, we employ a strategy that takes CHAN, an IGS station, as the region center and the radius is 300, 400, 500, 600, 700, 800, 900 and 1050 km respectively. A series of concentric circles contain 19 CMONOC stations that locate in northeastern three provinces and eastern part of Inner Mongolia. The information of the total 20 stations is given in Table 3   We use PCA to decompose the residual series of stations that locate in circles with different radiuses. Figure 5 demonstrates the relation between region size and the contribution percentage of PC1 in extracting CMC. In this region, the contribution percentage of PC1 in horizontal components (North and East) remains steady as the region size increases gradually. In the East, it has a slight  We use PCA to decompose the residual series of stations that locate in circles with different radiuses. Figure 5 demonstrates the relation between region size and the contribution percentage of PC1 in extracting CMC. In this region, the contribution percentage of PC1 in horizontal components (North and East) remains steady as the region size increases gradually. In the East, it has a slight decrease from 300 km to 1050 km, while in the North it is nearly unchanged and arrives 50% slower. However, the contribution percentage of PC1 in vertical component (Up) has a strong relation with region size. It drops significantly along with the increase of distance, from 78% to only 49%. Furthermore, we notice that the change trend of the Up component falls dramatically between 600 km and 700 km, with the figure going from 66% to 58%. In summary, CMC can be extracted accurately in a region of 600 km, and the contribution percentages of PC1 in North, East and Up components are more than 50%.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 17 decrease from 300 km to 1050 km, while in the North it is nearly unchanged and arrives 50% slower. However, the contribution percentage of PC1 in vertical component (Up) has a strong relation with region size. It drops significantly along with the increase of distance, from 78% to only 49%. Furthermore, we notice that the change trend of the Up component falls dramatically between 600 km and 700 km, with the figure going from 66% to 58%. In summary, CMC can be extracted accurately in a region of 600 km, and the contribution percentages of PC1 in North, East and Up components are more than 50%.

Sub-Regional Extractions Considering Tectonic Units
It is necessary to combine the CMONOC stations around IGS stations in order to acquire more accurate CMC from IGS stations in the China region. Above we discussed the influence of the region size on CMC extraction, and this may provide a criterion to choose suitable CMONOC stations. However, the Northeast Block is the most stable throughout China and the established criterion may not necessarily hold for other blocks. In cases where the region is less stable, we suggest that there is a strong need for another criterion, namely geology, to further restrain selecting stations. In this paper, we raise second-class tectonic units as the geology criterion [37,38] because CMC can be extracted more accurately in a relatively small-sized and stable region according to Section 3.2. What's more, the self-condition of stations should also be combined together with these selection criteria, such as the completeness of the coordinate time series, time span, and whether there is some short-term drastic land subsidence or some unexplained offsets.
Criteria for selecting CMONOC stations are protocoled in this paper as follows. First, interstation distance criterion: The distance between one CMONOC station and IGS station does not exceed 600 km. Second, geology criterion: CMONOC stations and IGS station are in the same secondclass tectonic unit. Third, self-condition criterion: It means that coordinate time series have a time span of over 5 years with 60% completeness, and meanwhile without apparent tectonic movements.
According to these criteria, some stations are discarded. They include XZCY and XZGZ whose distances to LHAZ are more than 600 km in the Lhasa Sub-block. TJBH and TJWQ in Tianjin are also abandoned where there is drastic land subsidence because of groundwater exploitation. There are also some stations with unexplained offsets, such as SXCZ and NMBT that are discarded. After the selecting process, we get the sub-regional division of IGS stations as shown in Figure 6. Table 4 gives information about the eight divided sub-regions. Table 4. Information of eight divided sub-regions when extracting CMC of IGS stations combining CMONOC stations. Note that GUAO and URUM are in the same sub-region, and three IGS stations in Taiwan (TCMS, TNML and TWTF) are in the same sub-region.

Sub-Regional Extractions Considering Tectonic Units
It is necessary to combine the CMONOC stations around IGS stations in order to acquire more accurate CMC from IGS stations in the China region. Above we discussed the influence of the region size on CMC extraction, and this may provide a criterion to choose suitable CMONOC stations. However, the Northeast Block is the most stable throughout China and the established criterion may not necessarily hold for other blocks. In cases where the region is less stable, we suggest that there is a strong need for another criterion, namely geology, to further restrain selecting stations. In this paper, we raise second-class tectonic units as the geology criterion [37,38] because CMC can be extracted more accurately in a relatively small-sized and stable region according to Section 3.2. What's more, the self-condition of stations should also be combined together with these selection criteria, such as the completeness of the coordinate time series, time span, and whether there is some short-term drastic land subsidence or some unexplained offsets.
Criteria for selecting CMONOC stations are protocoled in this paper as follows. First, interstation distance criterion: The distance between one CMONOC station and IGS station does not exceed 600 km. Second, geology criterion: CMONOC stations and IGS station are in the same second-class tectonic unit. Third, self-condition criterion: It means that coordinate time series have a time span of over 5 years with 60% completeness, and meanwhile without apparent tectonic movements.
According to these criteria, some stations are discarded. They include XZCY and XZGZ whose distances to LHAZ are more than 600 km in the Lhasa Sub-block. TJBH and TJWQ in Tianjin are also abandoned where there is drastic land subsidence because of groundwater exploitation. There are also some stations with unexplained offsets, such as SXCZ and NMBT that are discarded. After the selecting process, we get the sub-regional division of IGS stations as shown in Figure 6. Table 4 gives information about the eight divided sub-regions.  We use PCA to decompose the residual series of the stations in each sub-region, and contribution percentage of PC1 is shown in Figure 7. After analyzing the eigenvalues and normalized eigenvectors in all sub-regions, we define PC1 as CMC in each sub-region. From Figure 7, apart from a few coordinate components, the contribution percentages of PC1 are over 40% in most components. Especially in Lhasa Sub-block, the indicators on three coordinate components are all over 70%, representing that most of CMC there can be extracted precisely by PCA filtering. The average values on the N, E and U component are 54.6%, 54.9% and 54.7%, which increases by 10.8%, 16.1% and 25.1% respectively, compared with the results of mode 1. From this, we know that this indicator, contribution percentages of PC1, can be an important assessment criterion of CMC extraction.   We use PCA to decompose the residual series of the stations in each sub-region, and contribution percentage of PC1 is shown in Figure 7. After analyzing the eigenvalues and normalized eigenvectors in all sub-regions, we define PC1 as CMC in each sub-region. From Figure 7, apart from a few coordinate components, the contribution percentages of PC1 are over 40% in most components. Especially in Lhasa Sub-block, the indicators on three coordinate components are all over 70%, representing that most of CMC there can be extracted precisely by PCA filtering. The average values on the N, E and U component are 54.6%, 54.9% and 54.7%, which increases by 10.8%, 16.1% and 25.1% respectively, compared with the results of mode 1. From this, we know that this indicator, contribution percentages of PC1, can be an important assessment criterion of CMC extraction.

Comparison of Two CMC Extraction Modes
In addition to the contribution percentage of PC1 as an indicator, there is another indicator which can be used to evaluate the effect of CMC extraction under the above two modes. This indicator is RMS reduction rate of residual series before and after PCA filtering (shown in Equation 7). Figure  8 shows RMS reduction rate of residual time series in Chinese IGS stations under cases of mode 1, mode 2, and the difference between the two modes:  Figure 8b). The results show that PCA filtering is indeed an effective method to improve the signalto-noise ratio (SNR) of GPS coordinate time series. Meanwhile, CMC extraction mode that proposed in this paper (mode 2) is more robust than overall filtering of an entire region (mode 1). Figure 8c shows the difference in RMS reduction rate for the two modes. There are 26 coordinate components (accounting for 78.8%) have positive values in the difference of RMS reduction rate, indicating that mode 2 is a more effective filtering method on most coordinate components.

Comparison of Two CMC Extraction Modes
In addition to the contribution percentage of PC1 as an indicator, there is another indicator which can be used to evaluate the effect of CMC extraction under the above two modes. This indicator is RMS reduction rate of residual series before and after PCA filtering (shown in Equation (7)). Figure 8 shows RMS reduction rate of residual time series in Chinese IGS stations under cases of mode 1, mode 2, and the difference between the two modes:

Comparison of Two CMC Extraction Modes
In addition to the contribution percentage of PC1 as an indicator, there is another indicator which can be used to evaluate the effect of CMC extraction under the above two modes. This indicator is RMS reduction rate of residual series before and after PCA filtering (shown in Equation 7). Figure  8 shows RMS reduction rate of residual time series in Chinese IGS stations under cases of mode 1, mode 2, and the difference between the two modes:  Figure 8b). The results show that PCA filtering is indeed an effective method to improve the signalto-noise ratio (SNR) of GPS coordinate time series. Meanwhile, CMC extraction mode that proposed in this paper (mode 2) is more robust than overall filtering of an entire region (mode 1). Figure 8c shows the difference in RMS reduction rate for the two modes. There are 26 coordinate components (accounting for 78.8%) have positive values in the difference of RMS reduction rate, indicating that mode 2 is a more effective filtering method on most coordinate components.

Relevant Geophysical Signals
The nonlinear variations of GPS coordinate time series contain a variety of geophysical signals. Current research focuses are on the displacements of GPS stations under the impact of surface environmental loadings such as atmospheric pressure, oceanic tide and hydrology [19,39].  Figure 8b). The results show that PCA filtering is indeed an effective method to improve the signal-to-noise ratio (SNR) of GPS coordinate time series. Meanwhile, CMC extraction mode that proposed in this paper (mode 2) is more robust than overall filtering of an entire region (mode 1). Figure 8c shows the difference in RMS reduction rate for the two modes. There are 26 coordinate components (accounting for 78.8%) have positive values in the difference of RMS reduction rate, indicating that mode 2 is a more effective filtering method on most coordinate components.

Relevant Geophysical Signals
The nonlinear variations of GPS coordinate time series contain a variety of geophysical signals. Current research focuses are on the displacements of GPS stations under the impact of surface environmental loadings such as atmospheric pressure, oceanic tide and hydrology [19,39].
In this section, based on the grid data of loading deformation provided by GFZ, we adopt bicubic interpolation to calculate the displacement series at 11 IGS stations in China. Taking CHAN station as an example, Figure 9 shows the station displacement caused by NTAL, HYDL, NTOL and the sum of these displacements, SUML. In order to further explore the relation between environmental loading and CMC, and the source of CMC, we remove the annual and semi-annual terms in Equation (1) to obtain the residual series, with the seasonal periodic term in them. Thus, it helps to study the characteristics of non-linear variations in GPS coordinate time series [33]. Using the same methods and steps as in the previous section, we get the CMC series in mode 2. These series contain significant seasonal periodic signals.

Relevant Geophysical Signals
The nonlinear variations of GPS coordinate time series contain a variety of geophysical signals. Current research focuses are on the displacements of GPS stations under the impact of surface environmental loadings such as atmospheric pressure, oceanic tide and hydrology [19,39]. . Among these three loadings, NTAL contributes the most to SUML; the second contributor is HYDL; and NTOL contributes the least to SUML. SUML represents obvious seasonal variations, and its amplitude can be up to 10mm in some years.
In this section, based on the grid data of loading deformation provided by GFZ, we adopt bicubic interpolation to calculate the displacement series at 11 IGS stations in China. Taking CHAN station as an example, Figure 9 shows the station displacement caused by NTAL, HYDL, NTOL and the sum of these displacements, SUML. In order to further explore the relation between environmental loading and CMC, and the source of CMC, we remove the annual and semi-annual terms in Equation (1) to obtain the residual series, with the seasonal periodic term in them. Thus, it helps to study the characteristics of non-linear variations in GPS coordinate time series [33]. Using the same methods and steps as in the previous section, we get the CMC series in mode 2. These series contain significant seasonal periodic signals. . Displacements on CHAN station caused by geophysical signals, including atmospheric loading (NTAL), hydrological loading (HYDL), non-tidal oceanic loading (NTOL) and their sum (SUML). Note that the scale of the vertical axis in (c) is smaller than (a,b,d) The displacements caused by environmental loadings in vertical component (Up) are much bigger than that in horizontal components (North and East). Among these three loadings, NTAL contributes the most to SUML; the second contributor is HYDL; and NTOL contributes the least to SUML. SUML represents obvious seasonal variations, and its amplitude can be up to 10 mm in some years.
PCA filtering and environmental loading corrections (ELCs) are two commonly used methods for analyzing non-linear variations of GPS coordinate time series because both can effectively reduce SNR of GPS coordinate time series. In order to evaluate whether PCA filtering under mode 2 or ELCs is more effective for improving SNR of the GPS coordinate time series (especially in the vertical component), we calculate RMS reduction rate of residual series before and after PCA filtering and ELCs. The results are shown in Figure 10. As shown in Figure 10a, RMS reduction rates of all 33 coordinate components after PCA filtering are positive values, while 23 coordinate components after ELCs are positive, and RMS reduction rate of the remaining 30.3% of the coordinate components corrected by environmental loadings become larger (shown as some columns are below the horizontal axis in Figure 10b). This change in RMS reduction rates is consistent with the results from Jiang et al. [30]. Especially for the three vertical coordinate components shown in Figure 10b, which have negative RMS reduction rates. Given the fact that there is an apparent phase difference between the periodic terms of environmental loading displacement and GPS coordinate time series in vertical component, there is less RMS reduction rate for ELCs than PCA filtering. This reveals the high efficiency of PCA filtering under mode 2 that we demonstrate in this paper, and this PCA filtering under mode 2 is superior to ELCs in improving the SNR of the GPS coordinate time series. loading and CMC, and the source of CMC, we remove the annual and semi-annual terms in Equation (1) to obtain the residual series, with the seasonal periodic term in them. Thus, it helps to study the characteristics of non-linear variations in GPS coordinate time series [33]. Using the same methods and steps as in the previous section, we get the CMC series in mode 2. These series contain significant seasonal periodic signals. PCA filtering and environmental loading corrections (ELCs) are two commonly used methods for analyzing non-linear variations of GPS coordinate time series because both can effectively reduce SNR of GPS coordinate time series. In order to evaluate whether PCA filtering under mode 2 or ELCs is more effective for improving SNR of the GPS coordinate time series (especially in the vertical component), we calculate RMS reduction rate of residual series before and after PCA filtering and ELCs. The results are shown in Figure 10. As shown in Figure 10a, RMS reduction rates of all 33 coordinate components after PCA filtering are positive values, while 23 coordinate components after ELCs are positive, and RMS reduction rate of the remaining 30.3% of the coordinate components corrected by environmental loadings become larger (shown as some columns are below the horizontal axis in Figure 10b). This change in RMS reduction rates is consistent with the results from Jiang et al. [30]. Especially for the three vertical coordinate components shown in Figure 10b, which have negative RMS reduction rates. Given the fact that there is an apparent phase difference between the periodic terms of environmental loading displacement and GPS coordinate time series in vertical component, there is less RMS reduction rate for ELCs than PCA filtering. This reveals the high efficiency of PCA filtering under mode 2 that we demonstrate in this paper, and this PCA filtering under mode 2 is superior to ELCs in improving the SNR of the GPS coordinate time series.

Conclusions
Combining some selected CMONOC stations, we employ a PCA filtering method to compute regional CMC of IGS stations in China under two extraction modes, one whole extraction of all eleven IGS stations and eight sub-regional extractions combined with near CMONOC stations. The following conclusions may be drawn: (1) The average PC1 contribution percentage of eight sub-regions in North, East and Up components are 54.6%, 54.9% and 54.7%, which increase by 10.8%, 16.1% and 25.1%, respectively,

Conclusions
Combining some selected CMONOC stations, we employ a PCA filtering method to compute regional CMC of IGS stations in China under two extraction modes, one whole extraction of all eleven IGS stations and eight sub-regional extractions combined with near CMONOC stations. The following conclusions may be drawn: (1) The average PC1 contribution percentage of eight sub-regions in North, East and Up components are 54.6%, 54.9% and 54.7%, which increase by 10.8%, 16.1% and 25.1%, respectively, compared with the result of one whole region. The regional CMC extracted in sub-regional mode are more significant, and RMS reduction rates on more coordinate components are positive values after CMC extraction. Moreover, the differences of RMS reduction rates are positive on 78.8% coordinate components, indicating that PCA filtering in a sub-regional mode is a more effective filtering method for most coordinate components. This will help to improve SNR in GPS coordinate time series, and thus to reveal the surface deformation characteristics of the IGS station more realistically and accurately. (2) Regional CMC is a spatiotemporal-dependent displacement series. It is affected by factors, such as the region size and time span of the coordinate time series. When the CMC of the IGS station in China is uniformly extracted, the length of the coordinate time series is inversely proportional to the percentage of PC1. The longer the time span is, the smaller the percentage of PC1. When the time series span is greater than 10 years, the percentage of PC1 in North and East components are about 34%, and about 26% in the Up component. In the Northeast Block of China, regional