*Article* **Potential Contributors to Common Mode Error in Array GPS Displacement Fields in Taiwan Island**

**Xiaojun Ma 1,2, Bin Liu 1,2,\*, Wujiao Dai 3, Cuilin Kuang <sup>3</sup> and Xuemin Xing 1,2**


**Abstract:** The existence of the common mode error (CME) in the continuous global navigation satellite system (GNSS) coordinate time series affects geophysical studies that use GNSS observations. To understand the potential contributors of CME in GNSS networks in Taiwan and their effect on velocity estimations, we used the principal component analysis (PCA) and independent component analysis (ICA) to filter the vertical coordinate time series from 44 high-quality GNSS stations in Taiwan island in China, with a span of 10 years. The filtering effects have been evaluated and the potential causes of the CME are analyzed. The root-mean-square values decreased by approximately 14% and 17% after spatio-temporal filtering using PCA and ICA, respectively. We then discuss the relationship between the CME sources obtained by ICA and the environmental loads. The results reveal that the independent displacements extracted by ICA correlate with the atmospheric mass loading (ATML) and land water storage mass loading (LWS) of Taiwan in terms of both its amplitude and phase. We then use the white noise plus power law noise model to quantitatively estimate the noise characteristics of the pre- and post-filtered coordinate time series based on the maximum likelihood estimation criterion. The results indicate that spatio-temporal filtering reduces the amplitude of the PL and the periodic terms in the GPS time series.

**Keywords:** GPS time-series analysis; common mode error; independent component analysis; seasonal signals; surface mass loading

#### **1. Introduction**

Continuous global navigation satellite system (GNSS) networks have been widely used in the fields of geodesy and geophysics. The permanently operating ground reference stations provide not only a millimeter-level 3-D coordinate time series under a global reference frame, but also provide important observational data and quantitative constraints for studies involving crustal deformation and dynamic processes on a global and regional scale [1,2]. However, the GNSS data processing procedure is affected by many factors, and numerous errors occur in the GNSS terminal products [3–5]. Geophysical signals can be difficult to separate from noise and unknown error sources because they are hidden in a detrended residual coordinate time series. Previous studies have identified that a spatially correlated error generally exists in regional networks, caused by the reference frame, satellites orbits, ocean tide correction models, and other unknown errors. This is usually referred to as the common mode error (CME) [6–11]. The presence of this error affects the accuracy of the station coordinates and velocity solutions, and conceals many weak and transient signals in a coordinate time series (e.g., deep magma motion, fault motion).

**Citation:** Ma, X.; Liu, B.; Dai, W.; Kuang, C.; Xing, X. Potential Contributors to Common Mode Error in Array GPS Displacement Fields in Taiwan Island. *Remote Sens.* **2021**, *13*, 4221. https://doi.org/10.3390/ rs13214221

Academic Editors: Gino Dardanelli and Shuanggen Jin

Received: 19 September 2021 Accepted: 19 October 2021 Published: 21 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Spatial filtering techniques can be used to mitigate the CME. There are presently two types of spatial filtering techniques, stacking [6,12] and empirical orthogonal function (EOF) decomposition [7], which have been widely used in the field of geodesy to detect transient signals. Wdowinski et al. [6] initially introduced stacking to detect co-seismic displacements in a residual position time series. However, Nikolaidis et al. [12] illustrated that stacking works satisfactorily when the CME in a regional network is essentially spatially homogeneous. Larger areas and weaker spatial correlations generally result in differences between the spatial response of the different stations to the CME. Previous studies have illustrated that the CME spatial correlation decreases with increasing distance and reaches 95% within 1000 km. For networks that exceed 2000 km, the CMEs exhibit significant differences in various regions and the spatial correlation is essentially nonexistent beyond 6000 km [13,14]. To address this problem, Tian et al. [15] proposed a correlation weighting stacking scheme. However, this algorithm is still affected by local-scale or site-specific changes, and further research is necessary to obtain the appropriate weights of all the sites at each epoch.

In comparison, Dong et al. [7] and Serpelloni et al. [16] verified that the EOF can decompose residual coordinate time series into different modes to extract the CME, which provides a stricter mathematical framework than stacking. Principal component analysis (PCA) does not require CME spatial uniformity, but adopts a uniform temporal function that affects stations across a regional network. Dong et al. [7] reviewed the spatial consistency of stations using a visual appraisal of the decomposed principal components (PCs) to extract the CME from the Southern California Integrated GPS Network (SCIGN). PCA and its extensions have been subsequently used to eliminate the CME in local area networks [8,17–24].

Previous studies focused on the effective removal of CME from a GNSS coordinate time series to improve the accuracy of coordinate solutions from regional network measurements [25], which is conducive to the study of some geophysical applications using GNSS, such as the assessment of glacial isostatic adjustment models [26] and crustal movement [27,28]. Other studies have addressed the potential contributors of CME, particularly the potential geophysical signals [18,21,29–38]. This is particularly relevant for crustal deformation caused by the redistribution of surface mass, such as strongly seasonal water movements and atmospheric pressure changes. Kumar et al. [39] extracted the CME of Taiwan using the EOF, and quantitatively analyzed the CME in the U direction and atmospheric mass loading (ATML) of Taiwan, as calculated by the International Mass Loading Service (IMLS). They found a significantly high correlation and degree of concordance between the CME and ATML residuals for the vertical component. A further regression analysis revealed that the highest 90% of the non-seasonal ATML displacements in Taiwan are present in the CME variations.

However, the extraction of the CME using PCA still experiences some limitations for geophysical interpretation. PCA is based on the second-order statistics of variance and covariance, and therefore does not make full use of the higher-order CME statistics. PCA can work effectively in cases where there is only a single source in the CME. However, when the CME contains multiple competing sources, a single principal component (PC) may contain a mixture of contributions from different sources, which can easily generate artificial features. In contrast, Liu et al. [18,19,26], Ming et al. [38], and Yan et al. [29] applied an independent component analysis (ICA) to geodetic data sets and successfully separated seasonal and long-term trend signals.

Previous studies indicated the effectiveness of ICA in studying the potential contributors of CME in a regional GNSS coordinate time series. Thus, following the work of Kumar et al. [39], we apply a combination of PCA and ICA to the vertical coordinate time series of 44 stations from Taiwan island in China, and extract the CMEs of this network. The associated effects on the stations' coordinate residuals and observed values are analyzed and compared. We discuss the relationship between the CME sources and environmental loadings, and compare GPS CMEs with vertical surface displacement predictions from the

IMLS. We determine the cause of the Taiwan GNSS network CME, except for the atmosphere pressure loading (ATML), and separate the signals related to land water storage mass loading (LWS). The results further validate the accuracy of ICA applications for CME signal separation and the geophysical interpretation of regional networks. Changes in the optimal noise model of the GNSS stations before and after the ICA filtering are analyzed, as well as the differences between the time series parameters of the GNSS stations.

#### **2. Methods and Data Processing**

#### *2.1. GNSS Data Processing*

All of the stations in Taiwan's continuous GPS (CGPS) array are equipped with dualfrequency geodetic GNSS receivers (Trimble 4000 SST Geodetic II P and 4000 SSE/SSIGeodetic Surveyor) [40]. A station was usually occupied for more than two sessions, and each session included 6–14 h of GPS observations by tracking all available satellites that had risen higher than a 15◦ elevation angle. The sampling interval for the data logging is 15 s. The collected data are downloaded from the internal RAM of the receivers to a PC hard disk or floppy disk. The raw data of each station are then transferred to the Receiver INdependent EXchange (RINEX) format using a transfer program for post-processing. The GNSS data in RINEX format is converted to the Bernese format, compiled and processed, and solved with respect to the Penghu Structural Stabilization Station on the west coast [40–42]. Figure 1 displays the location distribution map of CGPS sites in Taiwan Province. In this study, the GPS measurements of 44 sites in Taiwan Province (red triangle in Figure 1) from 2006 to 2016 can be obtained from the GPS Lab web application (http://gps.earth.sinica.edu.tw, accessed on 2 August 2021) of Institute of Earth Sciences, Academia Sinica, Taiwan Province, which requires a simple registration process to access data. You can download threecomponent continuous GPS time series data for all sites through the "download time series" button.

In this study, we first selected data from the 283 available sites in the Taiwan CGPS array that have been observed for more than 10 years, and 109 sites that have been observed for less than 10 years. We then calculated the data length ratio of each site (*Ratioi* = *Li*/3652, *L* is the data length of the *i th* site, *<sup>i</sup>* = 1, 2, ··· , 392) based on the vertical coordinate time series data of each site from January 1, 2006 to December 31, 2015 (ten-year span, totaling 3652 days). We selected 188 sites with *Ratio* ≥ 90%. We adopted the approach of visual appraisal to screen and select 44 high-quality time series sites according to the data consistency and quality. Figure 2 shows the data span of the 44 sites used for the CME estimation, with nearly all of the sites having a data cycle longer than 9.5 years.

In the time series spatio-temporal analysis, GNSS coordinate time series require a series of pre-processing steps. Coordinate series processing mainly includes outlier detection and elimination, offsets detection, missing data interpolation, etc. Kurtosis is the non-Gaussian property of random variables that can be measured by the fourth-order moment of the variables, defined as:

$$kurt(x) = E\left\{\mathbf{x}^4\right\} - 3\left(E\left\{\mathbf{x}^2\right\}\right)^2\tag{1}$$

where *x* is an independent random variable and *E* is a mathematical expectation operator. In the ICA process, it is assumed that *x* is normalized, and *E* ! *y*2 " = 1. The kurtosis is simplified to *E* ! *y*4 " <sup>−</sup> 3. For Gaussian variables, the kurtosis *kurt*(*y*) <sup>=</sup> *<sup>E</sup>* ! *y*4 " <sup>−</sup> <sup>3</sup> <sup>=</sup> 3*E* ! *y*2 " − 3 = 0, which implies that the kurtosis of the Gaussian variable is 0, the super-Gaussian variable is positive, and the sub-Gaussian variable is negative. A stronger non-Gaussian expression of the variable is associated with a greater absolute value of kurtosis. It should be noted that the kurtosis calculation is highly sensitive to outliers [43–46]. Outliers in a GNSS time series must therefore be identified in advance. The triple standard difference method is used for outlier detection and elimination. An offset in a time series is most often caused by antenna replacements and co-seismic displacement deformation, which affects the analysis of periodic changes of the GNSS coordinate series and the estimation of some geophysical signals. Here, we use the most convenient manual visual inspection

method to eliminate the offsets [47]. The missing data in the time series are interpolated after eliminating the outliers and fitting the offsets. The initial values of the missing data are interpolated using a PCA iteration approach [26]. An initial interpolation result is obtained using an inverse distance weighting interpolation. PCA is then applied and the first several PCs for which the cumulative percentage exceeds 70% are used to reconstruct the series. The PCA reconstruction is repeated until the maximum absolute value of the difference between two iteration is less than 0.01 mm.

**Figure 1.** Taiwan continuous global positioning system (CGPS) network. The black, green, blue, and purple circles represent CGPS sites operated by the Central Weather Bureau (CWB), the Institute of Earth Sciences (IES), the Central Geological Survey (CGS) and the Water Resources Agency (WRA) and other research units. Among them, the sites represented by red triangles are used in this ICA/PCA study.


**Figure 2.** Data availability of the 44 GPS sites. Blanks represent the missing data.

In the data pre-processing of the 44 GNSS stations, we used Equation (2) to fit the time series of each station, and eliminated the linear term and offset term, to smoothen the observation, and input the missing data. The time series after unified treatment ranges from approximately −20 to 20 mm. We selected five stations with a notable linear term, offset term and missing data in the coordinate time series (CHKU, JUNA, S101, SFON and SHAN), as indicated on the left in Figure 3, and the pre-processed data are presented on the right. It can be seen that the linear term, offset term, and missing data are all corrected.

**Figure 3.** Original GPS time series (**left**) and preprocessing sequence of the CHKU, JUNA, S101, SFON, and SHAN stations (**right**).

The coordinate time series for each site can be modeled as follows to obtain the effects of seasonal variations and colored noise [12]:

$$y(t\_i) = a + bt\_i + c\sin(2\pi t\_i) + d\cos(2\pi t\_i) + e\sin(4\pi t\_i) + f\cos(4\pi t\_i) + \sum\_{j=1}^{n\_\mathcal{g}} g\_i H(t\_i - T\_{\bar{g}j}) + \varepsilon\_i \tag{2}$$

where *ti* for *i* = 1, 2, ··· , *N* is the daily solution epoch in units of years, *a* and *b* are the site position and linear rate, respectively, *c* and *d* are the annual term coefficients, *e* and *f* are the semi-annual coefficients, *H ti* <sup>−</sup> *Tgi* is a Heaviside step function (when *ti* ≥ *Tgi*, its value is 1; otherwise it is 0), *gi* is the change of rate from *ti* = *Tgi*, and *ε<sup>i</sup>* is the model residual and also represents noise. It is usually assumed that *ε<sup>i</sup>* is white noise; however, the analysis of a large amount of data shows that the noise in the coordinate time series obtained from GNSS continuous observation data contain both white noise and colored noise [5,48–54]. The velocity error in a GPS coordinate time series may be underestimated by factors of 5–11 if a pure white noise model is assumed [5]. If the influence of colored noise is not considered, a biased or inaccurate geophysical interpretation may be obtained.

#### *2.2. Principal Component Analysis*

PCA is an orthogonal decomposition method that combines a group of related data and decomposes the group into a group of linear uncorrelated orthogonal eigenvectors and corresponding PCs. Each of the vectors reflects some of the information from the original matrix to different degrees. The PCs are arranged in descending order of eigenvalues, and the first few PCs can generally capture most of the information from the original matrix. Therefore, the main functions of the PCA is to manage the statistical method of dimensionality reduction in mathematics, explore a small number of PCs to represent the signals of the original data, and reflect the characteristic modes contained within the original matrix as much as possible [7,20,55].

The daily station residual time series with a span of *m* days in a continuous GNSS network is constructed by *n* stations, *X ti*, *xj* (*i* = 1, 2, ··· , *m*; *j* = 1, 2, ··· , *n*). Firstly, we preprocess each data column in matrix *X* with trends and outliers, and obtain its covariance matrix *B* using the following equation:

$$b\_{ij} = \frac{1}{m-1} \sum\_{k=1}^{m} \varkappa\_{k,i} \varkappa\_{k,j} \tag{3}$$

where *xk*,*i*, *xk*,*<sup>j</sup>* respectively represent the time series corresponding to the *i th* and *j th* stations at epoch *k*. The symmetric matrix *B* can be decomposed as:

$$B = V\Lambda V^{\mathsf{T}} \tag{4}$$

where eigenvector matrix *<sup>V</sup>*<sup>T</sup> is an (*<sup>n</sup>* × *<sup>n</sup>*) matrix and <sup>Λ</sup> is a non-zero diagonal matrix of order *k*. For a complete rank covariance matrix *B*, *k* = *n*. The matrix *X* can be expressed as:

$$X = VP \tag{5}$$

where *<sup>P</sup>* is a (*<sup>m</sup>* × *<sup>n</sup>*) matrix. The *<sup>k</sup>*-th row vector in *<sup>P</sup>* is the *<sup>k</sup>th* PC for matrix *<sup>X</sup>* and represents a temporal signature of the mode. The *kth* column vector in *V* is the *kth* spatial response corresponding to the mapped PC, which represents the spatial distribution of the mode. The above decomposition method is known as an EOF analysis or a PCA. The matrix *P* can be denoted as:

$$P = V^{\mathrm{T}} X \tag{6}$$

We arrange the eigenvalues *λi*(*i* = 1, 2, ··· , *n*) in descending order, and the contribution of each PC to the original data set can be represented by the feature cumulative contribution rate *mk*(*k* = 1, 2, ··· , *n*):

$$m\_k = \frac{\sum\_{i=1}^k \lambda\_i}{\sum\_{i=1}^n \lambda\_i} \tag{7}$$

The first few PCs represent those that contribute the most to the variance of a participating time series, and are usually related to the common source of temporal function. Higher-order PCs usually relate to the influence of local or individual stations. The CME calculated by PCA can then be obtained according to:

$$\varepsilon^{\rm PCA} = \sum\_{k=1}^{q} v\_k p\_k (1 \le q < n) \tag{8}$$

where *q* is the number of PCs that define the CME.

#### *2.3. Independent Component Analysis*

ICA is a data-driven approach based on blind source separation (BSS). BSS only assumes that the signal sources are independent non-Gaussian signals, and the relationship between the information of the signal sources and the linear transformation is unknown. ICA is a digital signal processing algorithm developed to solve this problem. Compared to the model-driven approaches, ICA works without prior information about the underlying sources, which allows it to effectively detect some signals that cannot be obtained when using other methods [43–46].

It is assumed that signal sources are composed of several statistically independent signals that overlap in both temporal and spatial domains, ICA synchronously observes the overlapping signals using multiple channels and decomposes the observed signals into several ICs, after the unmixing, as a set of source signal estimates (Figure 4). The channels exert no influence on the signal, and the number of observed channels is the same as the number of signals. The standard noiseless ICA mathematical model is:

$$X(t) = AS(t) \tag{9}$$

where the random vectors *<sup>X</sup>*(*t*) <sup>=</sup> [*X*1(*t*), ··· , *XM*(*t*)]<sup>T</sup> represent the observed signals and random vectors *Si*(*t*), *i* = 1, 2, ··· , *N*; *M* ≥ *N* represent the source signals, and *A* is a mixing matrix. The ICA decomposition process can be regarded as the inverse transformation of Equation (9):

$$Y(t) = S(\hat{t}) = BX(t) \tag{10}$$

where *B* is the unmixing matrix. It is assumed that the row vectors in *S* are statistically independent of each other, their joint probability density function (pdf) is the product of their marginal probability density function, which means that, the joint entropy of each component is the sum of the entropy of each component, such that:

$$p(S) = \prod\_{i=1}^{n} p(S\_i) \tag{11}$$

where *p* is the pdf of *S*, where both *A* and the source signal *S* are unknown. However, as long as the output components separated by the unmixing matrix B are statistically independent of each other, this approach is equivalent to separating the source signals.

$$\begin{array}{ccccc} \mathbf{X}(t) - \mathbf{A}\mathbf{S}(t) & & & \mathbf{X} = \mathbf{A} \times \mathbf{S} \\ \mathbf{Y}(t) = \mathbf{B}\mathbf{X}(t) & & & \begin{matrix} \mathbf{Y} \\ \end{matrix} = \begin{matrix} \mathbf{A} \end{matrix} \times \mathbf{S} \\ \begin{matrix} \mathbf{Y} \\ \end{matrix} = \begin{matrix} \mathbf{B} \end{matrix} \times \mathbf{X} \\ \begin{matrix} \mathbf{Y} \\ \end{matrix} = \begin{matrix} \mathbf{B} \end{matrix} \times \mathbf{X} \\ \mathbf{X} = \mathbf{A} \end{matrix} \end{array}$$

**Figure 4.** Schematic illustration of independent component analysis.

The basic hypothesis of ICA is that the source signals are statistically independent, allowing for no more than one source with a Gaussian distribution while describing the remaining sources with non-Gaussian distributions. Most of the deformed signals in GNSS coordinate sequences are also non-Gaussian signals, and the non-Gaussian signals are quantized by negentropy. The ICA attempts to create various linear transformations on the observed signals, and the maximum negentropy of the transformed quantity is likely to indicate the source signals [43,46]. Various ICA algorithms have been derived based on this fundamental idea. In this study, we use the FastICA algorithm based on negentropy to estimate the original signals [56,57]. The negentropy is defined as:

$$f(X) = \int\_{-\infty}^{+\infty} p\_G(X) \ln p\_G(X) dX + \int\_{-\infty}^{+\infty} p(X) \ln p(X) dX \tag{12}$$

where *p*(*X*) is the pdf of *X*, and *pG*(*X*) is the pdf of the Gaussian distribution with the same covariance matrix as *p*(*X*). The greater negentropy values are associated with stronger non-Gaussian signals [57–59].

A detailed description of the FastICA algorithm procedure is as follows:


We can similarly reconstruct the original data set of a GNSS time series of *n* stations across a period of *m* days and calculate the CME using the source signals that have been separated by ICA. We can obtain ICs by *Y*(*t*) = *BX*(*t*), assuming that *A* = *B*<sup>−</sup>1, and *X* can be obtained by:

$$X(t) = B^{-1}Y(t) = AY(t) \tag{13}$$

The CME based on ICA can thus be given as:

$$\mathfrak{e}^{\text{ACA}} = \sum\_{k \in R}^{R} a\_k y\_k \tag{14}$$

where *R* is a group of temporal components in spatio-temporal filtering.

#### **3. Results**

#### *3.1. Common Mode Error Extraction Using PCA/ICA*

To intuitively compare the amplitude of each component, the corresponding spatial response is usually divided by the maximum absolute value and scaled to a variation interval of −100% ∼ 100%, where the scaled amount is multiplied by the corresponding temporal components. An upward movement is a positive response and a downward is a negative response. The results of the first five PCs are shown on the left side of Figure 5, and their corresponding spatial responses are provided in the first row of Figure 6.

**Figure 5.** Five components extracted from Taiwan PCA (**left**) and ICA (**right**).

We use the standard CME definition provided by Dong et al. [7], in which most stations (50%) for a certain PC mode exhibit a clearly normalized spatial response (> 25%) and the eigenvalue of the mode exceeds 1% of the collective eigenvalues; it can therefore be considered as a common mode. The spatial response corresponding to the first PC in the first line of Figure 6 displays notably good regional consistency, and its average normalized amplitude (absolute value) is 77.3%. However, the 2nd–5th. PCs do not satisfy this condition. In terms of eigenvalues, Table 1 provides the percentage of the first 10 PC eigenvalues for the total eigenvalues and the cumulative contribution rate (eigenvalues). The first PC accounts for 24.7%, the second PC accounts for 6.5%, and the cumulative contribution rate of the eigenvalues tends to stabilize after the fifth order, which indicates that the first PC is able to represent the most common-mode components. In summary, we define the displacement component caused by the first PC as the CME of the entire GPS network.

**Figure 6.** Spatial response of five PCs (first row) and ICs (second row) at each site in Taiwan. The red arrow represents the positive spatial response, and the blue arrow represents the negative spatial response.


**Table 1.** Individual contribution rate and cumulative contribution rate of the top ten PCs.

We use ICA to process the five components from PCA. The principle is realized using ICA's fast fixed-point algorithm [57,60]. The first step of the algorithm is to centralize and whiten the observed values before processing. The mixing matrix is orthogonal so as to reduce the number of free parameters. The mixed signals are then linearly transformed to express an unrelated variable with a variance that is equal to 1 (whitening or spheroidizing). This step comprises the pretreatment process of ICA, which is realized by PCA. This is conducted because the number of PCs reserved for the ICA analysis should be lower than the data dimension. The process of trial and error is therefore used to select the appropriate quantity [61]. In this study, the effect of retaining the first five PCs becomes apparent.

The ranking of the ICs obtained by ICA does not correspond to a decrease in the variance. We therefore calculate the average ratio of the GPS\_IC displacement in each GPS site calculated by the ICA filtering to the observed time series, in accordance with the procedure described by Liu et al. [19] and Milliner et al. [62]. Smaller ratios are associated with contributions that are more significant, and the ICs are reordered in ascending order. The corresponding spatial response of each IC is normalized according to the absolute

value of its maximum value so that the comparison is more accurate. The IC time series and the corresponding spatial response value, after the IC, is multiplied by the normalization factor, shown on the right in Figure 5 and in the second row of Figure 6. The first two ICs are evidently the largest of the original data according to the previous definition of the CME, and their spatial response values display the same sign. We can therefore locate the displacement caused by these two components as the CME of the entire GPS network in ICA filtering.

After calculating the CME, we subtract its influence on each station from the original time series. The filtering effect can be described by the reduction rate of the root mean square (RMS) of the residual time series (Figure 7). The results in Table 2 reveal that after PCA and ICA filtering, the average RMS values are reduced by approximately 14% and 17%, respectively. These two methods can therefore be used to effectively improve the signal-to-noise ratio of the residual time series.

**Figure 7.** RMS reduction percentage of residual time series after ICA filtering.



Figure 7 shows the decline of the RMS value of the residual time series at each station after ICA filtering. Figure 7 also presents notable differences in the ICA filtering effect between the east and west. The RMS reduction rate in the west is relatively large, whereas that in the east is small, which is partially due to topography. There are many mountains on the eastern side of Taiwan and the stations are sparsely distributed, whereas the western side is relatively flat and the stations are relatively dense. Another explanation is the orogenic processes, as there is topographic uplift in the eastern region [63] and sinking in the southwest region due to groundwater extraction [64]. This phenomenon can also be observed in the second-order PCs obtained via PCA filtering, and the station response indicates notable local variation characteristics, as shown in PC2 in the first line of Figure 6.

#### *3.2. Noise Analysis of GNSS Coordinate Time Series*

In the noise analysis of the Taiwan GNSS time series, for the original time series of each station, we fit the offsets in advance, interpolate the original data, and then analyze the power spectrum of the noise series. The results are presented in Figure 8. Due to the large volume of data, this study uses the GS39 station as an example to elaborate on an explanation. It can be seen that the power spectrum at low and high frequencies of the noise series presents different peaks and demonstrates a clear linear trend term. The spectral energy in the low-frequency band is higher than that in high-frequency band. The spectral exponent of this component is −0.67, which indicates that the noise of this station is influenced by colored noise. Through piecewise fitting, the slope of this component is approximately approached, yielding −1 at a low frequency and 0 at a higher frequency. The low-frequency result possesses the characteristics of flicker noise, and the high-frequency result possesses the characteristics of white noise. It can thus be estimated that the noise type of this point is a combination type of white noise and flicker noise.

**Figure 8.** U component power spectrum of the GS39 site noise sequence. The black line represents the power spectrum, and the red line represents the slope.

#### **4. Discussion**

#### *4.1. Potential Geophysical Interpretation of the CME*

We note that although the RMS reduction results are similar, the PCA and ICA filtering processes exhibit different spatial and temporal patterns. Previous studies have also revealed that PCA-derived CMEs are not completely unrelated to local effects or random noise. The decomposition of similar contribution components in an actual network residual time series is difficult to achieve using PCA, and we cannot identify potential geophysical mechanisms or study the subtle signals in GPS observations [29]. In contrast, there is a high correlation between the common mode components extracted by ICA and the simulated surface mass loading deformation. An important criterion to better understand and describe the differences between the ICA and PCA results is the extent to which the temporal and spatial patterns of annual crustal deformation, caused by mass loadings of various geophysical sources, are restored by using these two methods.

We submitted requests for real-time solutions for the ATML, and LWS 3-D displacement files from 44 stations in Taiwan between 1 January 2006 and 31 December 2015 using the online solution function of the IMLS (http://massloading.net/, accessed on 2 August 2021) [65–68]. For ATML and LWS, we used the GEOS-FPIT model developed by the Global Modeling and Assimilation Office at NASA Goddard Space Flight Center, which considers a time range from 1 January 2000 to present, that is updated several times per day. The model resolution is 0.50◦ × 0.625◦ × 72 layers × 3 h. We obtained the displacement time series in the center-of-figure frame.

We averaged all of the loading series displacements into the daily results, and we calculated the average GPS\_IC1 and GPS\_IC2,comparing it with the average ATML and LWS. As indicated in Figure 9, GPS\_IC is the outer product of the IC and its corresponding spatial response. We note that ATML and LWS are consistent in terms of their amplitude and phase with GPS\_IC1 and GPS\_IC2, respectively. The correlation coefficients of these two temporal patterns are 0.58 and 0.4, respectively. We therefore suggest that atmospheric and hydrological mass loading are the main components of the CME in the Taiwan GPS network, which can be reflected by the IC extracted by ICA spatio-temporal filtering. Although IC1 and IC2 can be interpreted by ATML and LWS, other temporal components are statistically independent, and their unretrieved information still requires further study.

**Figure 9.** Comparison of the GPS\_ICs (calculated with IC1 and IC2) and mass loading deformation signals (ATML and LWS). The upper figure shows GPS\_IC1 and the ATML, and the lower figure shows GPS\_IC2 and the LWS.

#### *4.2. Effect of Removing the CME*

We then use the Hector software [69] to analyze the 44 noise sequences with only the offsets corrected using the white-plus-power law (W + PL) model. We do not fix the spectral index, but treat it as an unknown factor, and solve it alongside the linear, annual, and semi-annual periodic terms [14]. On this basis, we subtract the CME after ICA filtering and analyze the obtained noise sequence, as shown in Table 3. We also analyze the power spectrum of the original noise sequence, PCA filtered noise sequence, and ICA filtered noise sequence of station GS39, as shown in Figure 10. The ICA filtering approach is observed is more useful than the PCA in terms of the amplitude of the annual term.

Table 3 lists the estimated parameters and noise terms of the 44 stations in Taiwan. We can see that the PL amplitude is reduced by an average of approximately 27.8% after ICA filtering, which indicates that the filtering has a certain suppression effect on noise. However, there is not an obvious change in the amplitude of the white noise (W), which indicates that it originated locally. The amplitude of annual and semi-annual periodic terms decreased by approximately 60% and 18%, respectively. Before and after the filtering, the average spectral index shifted from −0.77 to −0.92, which is consistent with the result of Rao et al. [70], indicating that the noise model of the Taiwan GPS data is closer to the white-plus-flicker noise (W+FN) model with a spectral index of −1.

**Table 3.** Estimated parameters using Hector software before and after spatiotemporal filtering (ICA).


**Figure 10.** Comparison of the power spectrum of the GS39 site noise sequence and noise sequence filtered by PCA and ICA.

#### **5. Conclusions**

We use PCA and ICA approaches to analyze the vertical coordinate time series of 44 sites in Taiwan island in China from 2006 to 2016. The results are summarized as follows.


**Author Contributions:** All authors contributed significantly to the manuscript. Conceptualization, X.M., B.L. and W.D.; methodology, X.M. and B.L.; validation, X.M. and B.L.; formal analysis, X.M.; investigation, X.M.; resources, B.L. and W.D.; data curation, X.X.; writing—original draft preparation, X.M.; writing—review and editing, B.L., W.D., C.K. and X.X.; supervision, W.D., C.K. and X.X.; funding acquisition, B.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is funded by the National Natural Science Foundation of China (Grant No. 41904003, 41774040, 41674011), the Natural Science Foundation of Hunan Province, China (Grant No. 2020JJ5571) and Open Fund of Engineering Laboratory of Spatial Information Technology of Highway Geological Disaster Early Warning in Hunan Province (Changsha University of Science & Technology) (Grant No. KFJ190602).

**Data Availability Statement:** The GPS data used in this work can be obtained from the GPS Lab web application (http://gps.earth.sinica.edu.tw, accessed on 2 August 2021) of the Institute of Earth Sciences, Academia Sinica, Taiwan. This application requires a simple registration process to access the data. You can download the three-component continuous GPS time series data of all the sites through the "download time series" button. The pre-calculated ATML and LWS displacement data sets were obtained from the International Mass Loading Service Web application (http://massloading. net/, accessed on 2 August 2021) Explore individual loadings, used in the "GEOS-FPIT" model, the "Center of figure" frame, and the time range is in "Between 2006.01.01\_00:00" and "2016.01.01\_00:00". Site file uses our solved XYZ coordinate file. The displacement is provided every 3 h and is stored in a total of 120 .eph files at all 44 sites. The file name contains the corresponding timestamp. We extract the U-component displacement from each file to form the original time series data set for our spatio-temporal analysis. The time series of the GPS stations used in this paper are uploaded to Figshare, and readers can download the data from: https://figshare.com/articles/dataset/data\_for\_ EES\_rar/13490442, accessed on 2 August 2021.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

