Extracting Common Mode Errors of Regional GNSS Position Time Series in the Presence of Missing Data by Variational Bayesian Principal Component Analysis.

Removal of the common mode error (CME) is very important for the investigation of global navigation satellite systems’ (GNSS) error and the estimation of an accurate GNSS velocity field for geodynamic applications. The commonly used spatiotemporal filtering methods normally process the evenly spaced time series without missing data. In this article, we present the variational Bayesian principal component analysis (VBPCA) to estimate and extract CME from the incomplete GNSS position time series. The VBPCA method can naturally handle missing data in the Bayesian framework and utilizes the variational expectation-maximization iterative algorithm to search each principal subspace. Moreover, it could automatically select the optimal number of principal components for data reconstruction and avoid the overfitting problem. To evaluate the performance of the VBPCA algorithm for extracting CME, 44 continuous GNSS stations located in Southern California were selected. Compared to previous approaches, VBPCA could achieve better performance with lower CME relative errors when more missing data exists. Since the first principal component (PC) extracted by VBPCA is remarkably larger than the other components, and its corresponding spatial response presents nearly uniform distribution, we only use the first PC and its eigenvector to reconstruct the CME for each station. After filtering out CME, the interstation correlation coefficients are significantly reduced from 0.43, 0.46, and 0.38 to 0.11, 0.10, and 0.08, for the north, east, and up (NEU) components, respectively. The root mean square (RMS) values of the residual time series and the colored noise amplitudes for the NEU components are also greatly suppressed, with average reductions of 27.11%, 28.15%, and 23.28% for the former, and 49.90%, 54.56%, and 49.75% for the latter. Moreover, the velocity estimates are more reliable and precise after removing CME, with average uncertainty reductions of 51.95%, 57.31%, and 49.92% for the NEU components, respectively. All these results indicate that the VBPCA method is an alternative and efficient way to extract CME from regional GNSS position time series in the presence of missing data. Further work is still required to consider the effect of formal errors on the CME extraction during the VBPCA implementation.


Introduction
In the last three decades, global navigation satellite system (GNSS) technology has provided abundant, high-accuracy position information for the Earth, which allows researchers to investigate many types of geophysical phenomena, such as geocenter motion [1,2], crustal deformation [3,4], seismic monitoring [5,6], and glacial isostatic adjustment [7,8]. However, the GNSS observations contain not only site-specific temporal noise, which is well quantified and described as a combination of white noise (WN) and power-law noises (PLN) using the maximum likelihood estimation (MLE) method [9], but also spatial correlation error known as the common mode error (CME) [10]. Interference between temporal noise and spatial correlation error would bias the station velocity estimation, increase velocity uncertainty, and even lead to incorrect geophysical interpretations; thus, it is highly recommended to mitigate the interference with some effective filtering techniques to enhance the GNSS signal-to-noise ratio.
For regional GNSS networks, strong CME, coming from the reference frame error, mis-modeling of satellite orbits and clocks, large-scale environmental effects, etc., presents in the position time series [11]. This makes it difficult to discern weak and transient tectonic signals, such as slow slip events and aseismic episodic tremor, from GNSS data. Until now, many methods have been proposed to filter out the CME, and Table 1 shows the overview of each mathematical method. The stacking approach was originally introduced to detect the seismic displacements [12]. After that, a weighted stacking method that takes individual position error into account was proposed [13]. This kind of stacking approach is based on the assumption of an equivalent CME for all stations; thus, it is not suitable for larger scale networks; e.g., continental networks, wherein CME exhibits inhomogeneous behavior. To overcome this restriction, Márquez-Azúa et al. took the length of the coordinate time series and site distance as weights to estimate the daily CME, and then removed it from the original GNSS time series to study the large-scale deformation of Mexico [14]. Tian et al. used the distances between neighboring sites and interstation correlations as a weighting scheme to extract the spatially correlated transients [15,16]. Nevertheless, the aforementioned methods are unable to identify stations with strong local effects that could affect the CME detection, and the suitable weights for all stations at each epoch still need to be investigated.
Another way to remove the CME is to transform the solutions into a regional reference frame [17,18]. In this way, the stations located at the edge of the region may exhibit distorted signals or some remnant CME [16]. To suppress the CME, more reliable, statistical signal decomposition techniques, such as principal component analysis (PCA) [19][20][21], multi-channel singular spectrum analysis (MSSA) [22,23], and independent component analysis (ICA) [24,25], have been introduced into GNSS time series analyses. Based on an assumed inhomogeneous distribution of CME and a more rigorous mathematical framework, these methods, in particular the PCA, have been widely used to eliminate CME from regional network. Nevertheless, the traditional PCA method normally requires a complete set of observations, indicating that interpolation should be implemented to recover the missing data beforehand, which may cause deviations in the extracted CME. Because of the limitation, Shen et al. [26] successfully developed an improved principal component analysis (IPCA) to extract CME from the discontinuous time series of 27 permanent GNSS stations in China. Then Li et al. [27,28] extended this approach by considering the formal errors as weights to extract CME. However, the IPCA approach is incapable of handling the condition that any two time series of a regional network have no, or only a few common epochs. Very recently, the probabilistic principal component analysis (PPCA) method was first introduced by Gruszczynski et al. [29] to filter out CME from time series with few common observational epochs, yet the PPCA method is sometimes sensitive to the initial parameters and prone to face the over-fitting problem due to the usage of maximum likelihood criterion to update model parameters [30]. To compensate for the above shortcomings, an alternative method called variational Bayesian principal component analysis (VBPCA) came into being [31]. The VBPCA method, as an extension of PPCA, imposed a priori distributions on the model parameters and estimated the hyperparameters by maximizing the evidence of observed signals instead of likelihood in the maximum-likelihood estimation paradigm. This property makes it more resistant against the overfitting problem in comparison with PPCA method. Moreover, the VBPCA method can recover the missing values in GNSS time series more conveniently, since it offers automatic dimensionality selection during the Bayesian inference procedure [30], while for other methods, i.e., PCA, IPCA, and PPCA, the appropriate number of principal components used to reconstruct the missing signal needs to be systematically quantified. In this paper, we utilize the VBPCA method to extract the CME from a regional GNSS position time series with missing values, and compare its performance with the other two approaches, PCA and IPCA. The PCA method is the most common algorithm with which to filter CME. The SOPAC uses this method to remove the CME in the PBO GNSS network. It can only handle a continuous residual coordinate matrix. Thus, it needs to fill the missing epoch using the spatially averaged values derived from other stations on these days. Then, traditional PCA analysis should be performed iteratively until convergence to refill in these missing points. In contrast, the IPCA algorithm can directly deal with an incomplete positional time series based on minimizing the weighted quadratic norm of principle component (PC) unknowns of an epoch. More details about PCA and IPCA algorithms can be found in Dong et al. [19] and Shen et al. [26], respectively. The formal errors of daily GNSS time series are not taken into account in this research, and its effect on CME extraction should be further explored. The remainder of this paper is organized as follows. Section 2 presents the VBPCA algorithm. Section 3 describes the data processing and demonstrates the effectiveness of this method by numerical synthetic and real time series analysis. Variations of the noise characteristics and the linear velocity uncertainty for the regional GNSS stations before and after removing CME are evaluated in Section 4; conclusions are summed up in Section 5.

Methodology
Since the GNSS measurement is often incomplete or missing in practice, the traditional approaches usually handle missing data by explicitly recovering them. The short-term missing data can be filled well by interpolation, yet long-term missing data are difficult to recover [26]. Moreover, most data interpolation methods focus only on single point time series, without regarding to the spatial correlation among points in a region [32]. This may result in incorrect extraction of CME. Here, we put forward the VBPCA algorithm to characterize and mitigate the impact of CME. The VBPCA method, including regularization, consideration of the noise term, and introduction of the a prior distribution over the model parameters, performs well in terms of finding the intrinsic dimension and the optimal number of clusters for the data with missing values [33].

Variational Bayesian PCA
In the context of VBPCA, the GNSS observed datum X = [x 1 , x 2 , · · · , x n ] ∈ R d×n is the stacking of all unfiltered residual time series x j that can be formulated as: where d indicates the number of epochs, n is the number of GNSS stations, W ∈ R d×k denotes the loading matrix, p j ∈ R k×1 are the latent variables or principal components, µ ∈ R d×1 is the bias vector, and e i ∈ R d×1 represents the noise term.
In addition to the assumption that the latent variables p j and noise e j are normally Gaussian distributed, the VBPCA also imposes a Gaussian a priori distribution on the parameters of µ and W to overcome the overfitting problem: where ε x is the noise variance, ε µ and ε W,j represent the a priori variance for µ and W :j separately, W :j is the j−th column of the loading matrix W, and I denotes the Identity Matrices. In the VBPCA framework, ε= {ε x , ε µ , ε W,j is regarded as the hyperparameter set, and θ= W, p, µ is treated as a hidden variable. Their values can be calculated by the variational expectation-maximum algorithm. Since the true posterior p(θ x, ε) is analytically intractable, it is usually approximated by using a simple probabilistic density function q(θ). In this way, the calculation of the posterior p(θ x, ε) is modified to update the approximation q(θ), so as to minimize the cost function in the E-step: Then the hyperparameters ε can be optimized with the obtained q(θ) to maximize the likelihood p(x ε) in the M-step.
For the sake of calculation convenience, q(θ) is used in the following form: where µ i is the i th element of the bias vector µ; w i is the column vector corresponding to the i th row of W; p j are the latent variables of the jth observations; and µ i and µ i are the posterior mean and variance of µ. w i and Σ w i are the posterior mean and variance of w i , and p j and Σ p j are the posterior mean and variance of p j , respectively. Thus, the cost function (Equation (3)) can be minimized by update alternatively one factor of q(θ i ), while other factors are held invariant. Finally, the model parameters and hyperparameters can be carried out by calculating the Equations (5)-(7) iteratively until convergence. The main steps and the corresponding flowchart of this algorithm ( Figure 1) are as follows: Step1: Update of the latent variables: Step2: Estimation of the bias vector and loading matrix: where x i,j represents the element i, j of matrix X, O represents the set of indices i, j for which x i,j is observed, O i is the set of indices j for which x i,j is observed, and |O i | denotes the number of elements in O i .
Step3: Calculation of variance parameters: where w ik stands for the kth element on the diagonal of Σ w i .

GNSS Data Processing
We selected 44 continuous GNSS stations located in Southern California ( Figure 2) to evaluate the VBPCA algorithm, and then compare its performance with the other two methods, including PCA and IPCA. The latest Scripps Orbit and Permanent Array Center (SOPAC) provides position time  The VBPCA can automatically select the optimal number of principal components for data reconstruct based on the hyperparameter ε W,i , which tends to be zero if the evidence of the corresponding principal component for reliable data modeling is weak [30]. This characteristic is called the automatic relevance determination (ARD) and makes it easy to recover the missing values in the GNSS time series. Furthermore, the unreliable reconstruction of missing data can also be detected since the VBPCA method provides uncertainty information for unknown quantities.

GNSS Data Processing
We selected 44 continuous GNSS stations located in Southern California ( Figure 2) to evaluate the VBPCA algorithm, and then compare its performance with the other two methods, including PCA and IPCA. The latest Scripps Orbit and Permanent Array Center (SOPAC) provides position time series spanning 14-years from 2005 to 2018 for the selected GNSS stations; they are used here during the experiment, so as to avoid the effects of mismodeling seasonal signals and noises on identifying CME [34]. Any position value with large formal error (>10, >10, >20 mm for the north, east, and up directions, respectively) is considered an outlier and first discarded from the raw time series. Then the linear trend, annual and semiannual terms, and the offsets are separately modeled and subtracted to obtain the residual position time series [35]. After that, the post-seismic displacements due to strong earthquakes are corrected, and the residuals that exceed three interquartile range (IQR) thresholds are also removed. Finally, we obtain the "clean" and unfiltered residual time series with an average missing rate of 5.25% for the experiment.

Comparison of CME Relative Errors from Different Methods
We first filled the data gaps of the 44 GNSS residual time series to obtain fully complete coordinate matrix X using the regularized expectation-maximization (RegEM) algorithm [36].
Then we randomly picked up 5-30% observations with 5% increments and removed them from the time series to simulate the incomplete data. In this way, we can evaluate the capability of the VBPCA algorithm to extract CMEs with different amounts of missing data. Note that the gaps in the real GNSS observations are usually lost over a successive time span rather than a completely random absence [32]; thus, in the synthetic experiment, some of the deleted elements were distributed continuously, while the others were lost completely at random. Finally, we separately employed the PCA, IPCA, and VBPCA methods to extract the CMEs from our simulated complete and incomplete data. Since the CMEs for a regional GNSS network should be constant over a time span no matter whether the data are missing or not, we treat the CMEs from the simulated complete time series as

Comparison of CME Relative Errors from Different Methods
We first filled the data gaps of the 44 GNSS residual time series to obtain fully complete coordinate matrix X using the regularized expectation-maximization (RegEM) algorithm [36]. Then we randomly picked up 5-30% observations with 5% increments and removed them from the time series to simulate the incomplete data. In this way, we can evaluate the capability of the VBPCA algorithm to extract CMEs with different amounts of missing data. Note that the gaps in the real GNSS observations are usually lost over a successive time span rather than a completely random absence [32]; thus, in the Sensors 2020, 20, 2298 7 of 19 synthetic experiment, some of the deleted elements were distributed continuously, while the others were lost completely at random. Finally, we separately employed the PCA, IPCA, and VBPCA methods to extract the CMEs from our simulated complete and incomplete data. Since the CMEs for a regional GNSS network should be constant over a time span no matter whether the data are missing or not, we treat the CMEs from the simulated complete time series as the "true" reference for comparison purpose, and then calculate the relative errors of CME as [29]: where E which equals 100 here, and indicates the repeated number of every experiment; CME 0 and CME k are the common mode errors before and after introducing gaps, respectively. Figure 3 presents the average relative errors of CME calculated by PCA, IPCA, and VBPCA methods with different percentages of deleted data for the north component. We can see that the relative errors of CME computed by all the three methods become larger with increasing percentages of missing data. Within all experiments, the performance of VBPCA in estimating the CME is better than that of PCA and IPCA. When there is less missing data (i.e., 5%), the relative errors of CME calculated by the three methods are not very significant with the maximum distinction of less than 1%. However, with increasing missing data, the three results become much different. In particular, when 30% of the observations are removed from the time series, the relative errors of CME using VBPCA reach the maximum of only 10.9%, indicating improvements of 2.7% and 5.1% compared with IPCA (13.6%) and PCA (16.0%), respectively. Similar results were obtained for the east and up components. Therefore, we conclude that the VBPCA method is an alternative, effective way to obtain the CME from incomplete GNSS position time series, especially for larger number of stations with a bigger missing data rate.
Sensors 2020, 20, x FOR PEER REVIEW 8 of 22 and PCA (16.0%), respectively. Similar results were obtained for the east and up components. Therefore, we conclude that the VBPCA method is an alternative, effective way to obtain the CME from incomplete GNSS position time series, especially for larger number of stations with a bigger missing data rate.

The Performance of Missing Value Estimation
An advanced feature of VBPCA is that it can automatically select the right number of principal components for missing data reconstruction. This property makes it convenient for GNSS missing value imputation. In order to evaluate the performance of missing value estimation, we selected twelve stations with missing data of less than 2% for artificial data deleting, and used the normalized root mean squared error (NRMSE) to assess the estimation accuracy. The introduced missing data are ranged from 5% to 30% of the total observation data with 5% increments, and the NRMSE is defined as:

The Performance of Missing Value Estimation
An advanced feature of VBPCA is that it can automatically select the right number of principal components for missing data reconstruction. This property makes it convenient for GNSS missing value imputation. In order to evaluate the performance of missing value estimation, we selected twelve stations with missing data of less than 2% for artificial data deleting, and used the normalized root mean squared error (NRMSE) to assess the estimation accuracy. The introduced missing data are ranged from 5% to 30% of the total observation data with 5% increments, and the NRMSE is defined as: where E is the same as Equation (5), N is the number of deleted data, x est k represents the estimation of the artificial deleted values, x obs indicates the real GNSS observed values corresponding to x est k , and var(·) defines the variance. The NRMSE value will come to 0 when the estimation is accurate, and 1 when the estimation is poor or the noise involved is too large. Figure 4 shows the NRMSE results of the three algorithms for the north component with various percentages of missing data. With rising missing data, the NRMSE values also present an increasing order for all the three methods, among which the VBPCA algorithm provides the most accurate results, followed by the IPCA method, while PCA exhibits the poorest. However, the difference between VBPCA and IPCA is not significant. Similar results are obtained for the east and up components; hence, we do not show the graphs here due to the limited space. Therefore, we conclude that the VBPCA method is more suitable to fill the missing values for GNSS time series, especially in the case of large amount of missing data.   Tables 2-4 show the a priori variance W, j ε corresponding to the first five principal components of the north, east, and up (NEU) components obtained by the VBPCA algorithm, respectively. We notice that the W,1 ε values, which separately interpret 54.71%, 60.84%, and 54.84% of the total variance for the NEU components, are all significantly larger than the others. W,2 ε contributes 10.00%, 6.43%, and 7.74% of the total variance, while the upper W, j ε accounts for no higher than 5.01% for each direction. We then scale the PCs/eigenvectors by the normalization factor similar to [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution.   Tables 2-4 show the a priori variance ε W, j corresponding to the first five principal components of the north, east, and up (NEU) components obtained by the VBPCA algorithm, respectively. We notice that the ε W,1 values, which separately interpret 54.71%, 60.84%, and 54.84% of the total variance for the NEU components, are all significantly larger than the others. ε W,2 contributes 10.00%, 6.43%, and 7.74% of the total variance, while the upper ε W, j accounts for no higher than 5.01% for each direction. We then scale the PCs/eigenvectors by the normalization factor similar to [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution. The mean responses are 89.32%, 83.34%, and 82.56% for the north, east, and up components, respectively, while the minimum values are 70.92% at station WWMT, 62.68% at station AZU1, and 71.32% at station GNPS. Quite different from the 1st PC, the 2nd and 3rd PCs show both positive and negative spatial responses in all directions, and the spatial patterns do not exhibit strong coherence over the network. In addition, some stations, e.g., MSOB in the north, WRHS in the east, and HOLP upward, show significantly larger responses compared with the other stations. We think that this may be due to local effects, unmodeled signals or errors, etc., which are uncommon for the entire regional stations. Therefore, in the following analysis, we will only use the 1st PC and its corresponding eigenvectors to recover the CME, and then subtract the CME from our regional GNSS residual time series to eliminate the spatial correlation of inter-stations.  [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution. The mean responses are 89.32%, 83.34%, and 82.56% for the north, east, and up components, respectively, while the minimum values are 70.92% at station WWMT, 62.68% at station AZU1, and 71.32% at station GNPS. Quite different from the 1st PC, the 2nd and 3rd PCs show both positive and negative spatial responses in all directions, and the spatial patterns do not exhibit strong coherence over the network. In addition, some stations, e.g., MSOB in the north, WRHS in the east, and HOLP upward, show significantly larger responses compared with the other stations. We think that this may be due to local effects, unmodeled signals or errors, etc., which are uncommon for the entire regional stations. Therefore, in the following analysis, we will only use the 1st PC and its corresponding eigenvectors to recover the CME, and then subtract the CME from our regional GNSS residual time series to eliminate the spatial correlation of inter-stations. 5.01% for each direction. We then scale the PCs/eigenvectors by the normalization factor similar to [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution. The mean responses are 89.32%, 83.34%, and 82.56% for the north, east, and up components, respectively, while the minimum values are 70.92% at station WWMT, 62.68% at station AZU1, and 71.32% at station GNPS. Quite different from the 1st PC, the 2nd and 3rd PCs show both positive and negative spatial responses in all directions, and the spatial patterns do not exhibit strong coherence over the network. In addition, some stations, e.g., MSOB in the north, WRHS in the east, and HOLP upward, show significantly larger responses compared with the other stations. We think that this may be due to local effects, unmodeled signals or errors, etc., which are uncommon for the entire regional stations. Therefore, in the following analysis, we will only use the 1st PC and its corresponding eigenvectors to recover the CME, and then subtract the CME from our regional GNSS residual time series to eliminate the spatial correlation of inter-stations. 5.01% for each direction. We then scale the PCs/eigenvectors by the normalization factor similar to [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution. The mean responses are 89.32%, 83.34%, and 82.56% for the north, east, and up components, respectively, while the minimum values are 70.92% at station WWMT, 62.68% at station AZU1, and 71.32% at station GNPS. Quite different from the 1st PC, the 2nd and 3rd PCs show both positive and negative spatial responses in all directions, and the spatial patterns do not exhibit strong coherence over the network. In addition, some stations, e.g., MSOB in the north, WRHS in the east, and HOLP upward, show significantly larger responses compared with the other stations. We think that this may be due to local effects, unmodeled signals or errors, etc., which are uncommon for the entire regional stations. Therefore, in the following analysis, we will only use the 1st PC and its corresponding eigenvectors to recover the CME, and then subtract the CME from our regional GNSS residual time series to eliminate the spatial correlation of inter-stations. 5.01% for each direction. We then scale the PCs/eigenvectors by the normalization factor similar to [19], and show the first three scaled PCs and its normalized eigenvectors of the VBPCA solution for the NEU components in Figure 5. It can be clearly seen that the normalized amplitudes of the first eigenvectors all have the same sign and show a positive spatial response for all the three components with uniform distribution. The mean responses are 89.32%, 83.34%, and 82.56% for the north, east, and up components, respectively, while the minimum values are 70.92% at station WWMT, 62.68% at station AZU1, and 71.32% at station GNPS. Quite different from the 1st PC, the 2nd and 3rd PCs show both positive and negative spatial responses in all directions, and the spatial patterns do not exhibit strong coherence over the network. In addition, some stations, e.g., MSOB in the north, WRHS in the east, and HOLP upward, show significantly larger responses compared with the other stations. We think that this may be due to local effects, unmodeled signals or errors, etc., which are uncommon for the entire regional stations. Therefore, in the following analysis, we will only use the 1st PC and its corresponding eigenvectors to recover the CME, and then subtract the CME from our regional GNSS residual time series to eliminate the spatial correlation of inter-stations.

Interstation Correlation Analysis
The distance correlation coefficient is applied in this research to measure the dependency of two GNSS stations before and after reducing CME from the residual time series, due to the reason that this method can detect both linear and nonlinear correlation between two variables, whereas the classic Pearson's correlation can only measure the linear relationship [37]. The values of distance correlate to the range between 0 and 1 (0 implies independence, 1 represents similarity). Figure 6 illustrates the interstation correlations among the 44 stations and their corresponding histograms before and after removing CME from the NEU residual time series. We notice that the 44 GNSS stations show a significant correlation before filtering, with average values of 0.43, 0.46, and 0.38 for the north, east, and up components, respectively, among which the maximum correlation reaches up to 0.71 between station SBCC and SCIA (about 120 km apart) for the east component. These results confirm that there are strong spatial correlations in the regional GNSS network, which need to be eliminated before scientific research and application. After removing CME with VBPCA, the average correlation coefficients are remarkably reduced to 0.11, 0.10, and 0.08, representing average interstation correlation reductions of 74.42%, 78.26%, and 78.95% for the NEU components, respectively. Hence, we conclude that the VBPCA filtering is an effective technique to reconstruct the CME for regional GNSS network.

Interstation Correlation Analysis
The distance correlation coefficient is applied in this research to measure the dependency of two GNSS stations before and after reducing CME from the residual time series, due to the reason that this method can detect both linear and nonlinear correlation between two variables, whereas the classic Pearson's correlation can only measure the linear relationship [37]. The values of distance correlate to the range between 0 and 1 (0 implies independence, 1 represents similarity). Figure 6 illustrates the interstation correlations among the 44 stations and their corresponding histograms before and after removing CME from the NEU residual time series. We notice that the 44 GNSS stations show a significant correlation before filtering, with average values of 0.43, 0.46, and 0.38 for the north, east, and up components, respectively, among which the maximum correlation reaches up to 0.71 between station SBCC and SCIA (about 120 km apart) for the east component. These results confirm that there are strong spatial correlations in the regional GNSS network, which need to be eliminated before scientific research and application. After removing CME with VBPCA, the average correlation coefficients are remarkably reduced to 0.11, 0.10, and 0.08, representing average interstation correlation reductions of 74.42%, 78.26%, and 78.95% for the NEU components, respectively. Hence, we conclude that the VBPCA filtering is an effective technique to reconstruct the CME for regional GNSS network. Sensors 2020, 20, x FOR PEER REVIEW 13 of 22 Figure 6. Interstation correlation coefficients (left panels) and the histograms (right panels) before and after filtering CME from residual time series for the north (top panels), east (middle panels), and up (bottom panels) components. The unfiltered and filtered interstation correlation coefficients are separately shown at the upper triangle and lower triangular of the left panels, while their corresponding histograms are shown as red and blue colors in the right panels, respectively. Table 5 lists the root mean square (RMS) values of the 44 GNSS residual time series before and after filtering CME with VBPCA for the NEU components, respectively. We can see that the RMS reduction ranges between 7.63% and 44.00% for all stations, with an average value of 27.11%, 28.15%, and 23.28% for the North, East, and Up components, respectively. In particular, station SBCC, PPBF and SCIA exhibits the maximum RMS reduction of 44.00% (1.15 vs. 0.64), 42.86% (1.20 vs. 0.69), and 33.35% (3.92 vs. 2.61) for the north, east, and up component.

Time Series Analysis
To further analyze the periodicity of the residual time series, we take station SCIA as an example and calculate its power spectra using the Lomb-Scargle periodogram before and after removing CME for the NEU components. The CME and its power spectra are plotted in the left panels of Figure 7, while the unfiltered and filtered time series, together with their corresponding power spectra are shown in the right panels. It can be clearly seen that the 1.04 cycle per year (cpy), which is called the draconitic harmonic, shows high power with significant peaks for all the three components in the CME power spectra, and the noise characteristics of the CME are proximate to pure flicker noise. The spectrum of the unfiltered time series shows a combination of colored noise at low frequencies and white noise at high frequencies for all three directions, while there is a clear decrease in the power spectra at low frequencies compared with that of the high frequencies, especially for several peaks of the draconitic harmonic frequencies after eliminating CME with VBPCA. This indicates that the CME of the selected region contains part of flicker noise and draconitic signals, which can be significantly reduced with VBPCA filtering. Table 5. RMS values of the 44 GNSS residual time series before and after removing CME (unit: mm).

No.
Station Figure 6. Interstation correlation coefficients (left panels) and the histograms (right panels) before and after filtering CME from residual time series for the north (top panels), east (middle panels), and up (bottom panels) components. The unfiltered and filtered interstation correlation coefficients are separately shown at the upper triangle and lower triangular of the left panels, while their corresponding histograms are shown as red and blue colors in the right panels, respectively. Table 5 lists the root mean square (RMS) values of the 44 GNSS residual time series before and after filtering CME with VBPCA for the NEU components, respectively. We can see that the RMS reduction ranges between 7.63% and 44.00% for all stations, with an average value of 27.11%, 28.15%, and 23.28% for the North, East, and Up components, respectively. In particular, station SBCC, PPBF and SCIA exhibits the maximum RMS reduction of 44.00% (1.15 vs. 0.64), 42.86% (1.20 vs. 0.69), and 33.35% (3.92 vs. 2.61) for the north, east, and up component.

Time Series Analysis
To further analyze the periodicity of the residual time series, we take station SCIA as an example and calculate its power spectra using the Lomb-Scargle periodogram before and after removing CME for the NEU components. The CME and its power spectra are plotted in the left panels of Figure 7, while the unfiltered and filtered time series, together with their corresponding power spectra are shown in the right panels. It can be clearly seen that the 1.04 cycle per year (cpy), which is called the draconitic harmonic, shows high power with significant peaks for all the three components in the CME power spectra, and the noise characteristics of the CME are proximate to pure flicker noise. The spectrum of the unfiltered time series shows a combination of colored noise at low frequencies and white noise at high frequencies for all three directions, while there is a clear decrease in the power spectra at low frequencies compared with that of the high frequencies, especially for several peaks of the draconitic harmonic frequencies after eliminating CME with VBPCA. This indicates that the CME of the selected region contains part of flicker noise and draconitic signals, which can be significantly reduced with VBPCA filtering.

Effect of CME on Noise Amplitude and Velocity Estimation
To quantitatively analyze the effects of CME filtering on the final velocity field estimation, we estimate the velocity together with the PLN+WN noise model using Hector software [38] for the total 44 stations. Figure 8 illustrates the variation (unfiltered minus filtered) of the amplitudes of PLN and WN, and the velocity, after removing CME by VBPCA algorithm for the NEU components. Figure 9 presents the corresponding velocity uncertainty before and after CME filtering. We observe that the PLN amplitudes exhibit a significant decrease after filtering out CME, with average values of 49.90%, 54.56%, and 49.75% in the NEU directions, respectively. However, the WN amplitudes only vary with average decrease ratios of 1.8% and 5.82% for the north and east components, but an increase ratio of Figure 7. Power spectra of CME, unfiltered and filtered residual time series of station SCIA for the north (top panels), east (middle panels), and up (bottom panels) components. Left panels: CME and its power spectra (blue). Right panels: unfiltered (gray) and filtered time series (red), together with their corresponding power spectra. The vertical black-solid lines indicate harmonics of 1.0 cpy, and brown-dash lines indicate harmonics of 1.04 cpy.

Effect of CME on Noise Amplitude and Velocity Estimation
To quantitatively analyze the effects of CME filtering on the final velocity field estimation, we estimate the velocity together with the PLN+WN noise model using Hector software [38] for the total 44 stations. Figure 8 illustrates the variation (unfiltered minus filtered) of the amplitudes of PLN and WN, and the velocity, after removing CME by VBPCA algorithm for the NEU components. Figure 9 presents the corresponding velocity uncertainty before and after CME filtering. We observe that the PLN amplitudes exhibit a significant decrease after filtering out CME, with average values of 49.90%, 54.56%, and 49.75% in the NEU directions, respectively. However, the WN amplitudes only vary with average decrease ratios of 1.8% and 5.82% for the north and east components, but an increase ratio of 27.65% for the up component. This is due to the fact that flicker noise is dominated in the CME and can be remarkably suppressed using VBPCA filtering, while white noise is likely related to site-specific errors that cannot be removed by CME filtering.
can be remarkably suppressed using VBPCA filtering, while white noise is likely related to sitespecific errors that cannot be removed by CME filtering.
The velocity differences vary between −0.27 mm/year and 0.16 mm/year for the horizontal components, and −0.36-0.21 mm/year for the up component (right panels of Figure 8), indicating that CME has a certain influence on the velocity estimation of the GNSS stations. Compared with velocity difference, the velocity uncertainties decline obviously with averages of 51.95%, 57.31%, and 49.92% after filtering out CME for the NEU components, among which the maximum improvement ratio equals 79.07% at station SRS1 for the east component, while the minimum improvement is 12.02% at station MSOB for the up component. Therefore, we conclude that CME filtering would be essential to improve the precision of velocity estimates. Figure 8. The variation of amplitudes of PLN (left panels) and WN (middle panels) before and after CME filtering with VBPCA, together with the velocity differences (right panels). From top to bottom panels represent the North, East and Up components.

Conclusions
In this paper, we introduce an alternative VBPCA approach to extract CME from regional GNSS position time series with missing values. Without establishing the covariance matrix and computing eigenvalue decomposition, the VBPCA method considers PCA from a probabilistic point of view, and utilizes the variational expectation-maximization iterative algorithm to search principal subspace. Moreover, The VBPCA can automatically select the optimal number of principal components for data reconstruction based on the hyperparameter W,i ε , which makes it convenient CME has a certain influence on the velocity estimation of the GNSS stations. Compared with velocity difference, the velocity uncertainties decline obviously with averages of 51.95%, 57.31%, and 49.92% after filtering out CME for the NEU components, among which the maximum improvement ratio equals 79.07% at station SRS1 for the east component, while the minimum improvement is 12.02% at station MSOB for the up component. Therefore, we conclude that CME filtering would be essential to improve the precision of velocity estimates.

Conclusions
In this paper, we introduce an alternative VBPCA approach to extract CME from regional GNSS position time series with missing values. Without establishing the covariance matrix and computing eigenvalue decomposition, the VBPCA method considers PCA from a probabilistic point of view, and utilizes the variational expectation-maximization iterative algorithm to search principal subspace. Moreover, The VBPCA can automatically select the optimal number of principal components for data reconstruction based on the hyperparameter ε W,i , which makes it convenient to recover the missing values in the GNSS time series.
The daily position time series of 44 continuous GNSS stations located in Southern California from SOPAC are selected to simulate and investigate the performance of filtering out CME with VBPCA. Our results show that the VBPCA algorithm achieves lower relative errors of CME compared with the other two widely used approaches; namely, PCA and IPCA. In particular, it shows maximum improvements of 2.7% and 5.1% compared to IPCA and PCA in the case of 30% missing data. With respect to the performance of missing value estimation, the VBPCA algorithm always exhibits the lowest NRMSE; thus, it is more suitable for filling the gaps of the GNSS time series, especially for a large amount of missing data.
We then apply the VBPCA algorithm to the 44 incomplete residual time series, and successfully recover the CME using the first PC and its eigenvector for the selected regional GNSS network. After removing CME, both the interstation correlation coefficients and the RMS values of the residual time series show remarkable decreases. The draconitic harmonics, especially for the top eight frequencies, also exhibit significant reductions. Therefore, we conclude that the CME indeed exists in the regional GNSS network and the VBPCA method performs quite well to reconstruct the CME for the selected region.
Meanwhile, the velocities of unfiltered and filtered position time series, together with their corresponding velocity uncertainties are simultaneously estimated with the PLN+WN noise model. We observe that after CME filtering, the PLN amplitudes are significantly reduced, with average values of 49.90%, 54.56%, and 49.75% in the NEU directions, respectively. However, the WN amplitudes exhibit only slight decreases in their horizontal components, but increases in the up direction. The velocity differences can reach 0.36 mm/year, indicating that CME has a certain influence on the velocity estimation of the GNSS stations. On the other hand, the velocity uncertainty represents an average reduction of 51.95%, 57.31%, and 49.92% for the three components; thus, we conclude that the CME filtering with VBPCA could not only attenuate the colored noise, but also improve the accuracy of the estimated station velocity. Further work is still required to take the formal errors into account during the VBPCA implementation, which may also have some effect on the CME extraction.
Finally, with the fast expansion of GNSS stations and the accumulated longer-time, or higher-rate GNSS observations, we are facing a challenge to find an optimal and objective approach to extract CME for regional GNSS network. The deep leaning method that is regarded as a learning technique on artificial neural networks can provide fast and effective solutions, especially in the analysis of big data. Considering the advantages of deep learning, its application in extracting the CME from GNSS position time series will be further explored.
Author Contributions: W.J. and W.L. developed the idea that led to this paper. Z.L. and W.L. researched the theory of VBPCA algorithm, and wrote the original draft. H.C. and Q.C. performed the GNSS data processing and analyzed the results. J.W. collected the GNSS data and visualized the results. G.Z. and W.L. wrote the source code of different methods to compare their extracted CME. In addition, W.L. and Z.L. further revised the manuscript according two reviewer's comments. All authors have read and agreed to the published version of the manuscript.