Spatiotemporal analysis of GPS time series in vertical direction using independent component analysis

GPS has been widely used in the field of geodesy and geodynamics thanks to its technology development and the improvement of positioning accuracy. A time series observed by GPS in vertical direction usually contains tectonic signals, non-tectonic signals, residual atmospheric delay, measurement noise, etc. Analyzing these information is the basis of crustal deformation research. Furthermore, analyzing the GPS time series and extracting the non-tectonic information are helpful to study the effect of various geophysical events. Principal component analysis (PCA) is an effective tool for spatiotemporal filtering and GPS time series analysis. But as it is unable to extract statistically independent components, PCA is unfavorable for achieving the implicit information in time series. Independent component analysis (ICA) is a statistical method of blind source separation (BSS) and can separate original signals from mixed observations. In this paper, ICA is used as a spatiotemporal filtering method to analyze the spatial and temporal features of vertical GPS coordinate time series in the UK and Sichuan-Yunnan region in China. Meanwhile, the contributions from atmospheric and soil moisture mass loading are evaluated. The analysis of the relevance between the independent components and mass loading with their spatial distribution shows that the signals extracted by ICA have a strong correlation with the non-tectonic deformation, indicating that ICA has a better performance in spatiotemporal analysis.


Introduction
With the improvement of positioning accuracy and the growing number of CORS, GPS has found an increasingly wide utilization in the field of geodesy and geodynamics, such as earthquake studies (Segall and Davis 1997), volcano deformation monitoring (Ueda et al. 2013), crustal movement research (Niu et al. 2005;Teferle et al. 2008), and so on. The GPS time series in the horizontal direction are featured with the linear velocity caused by plate movement, and in the vertical direction with more cyclical variation caused by geophysical phenomena such as non-tectonic deformation. Some researches show that the mass loading can explain most seasonal variations estimated in vertical GPS time series Wang et al. 2005;Zhang et al. 2002).
The analysis of GPS data in the regional GPS monitoring networks indicates that there is a spatial correlation among different stations that is usually caused by the socalled common-mode error (CME) (Wdowinski et al. 1997;Nikolaidis 2002;Dong et al. 2006). So, it is important to reduce the CME and improve the accuracy of the GPS coordinates using spatiotemporal filtering technique. Wdowinski et al. (1997) proposed using the stacking, a regional spatial filtering method, to calculate the CME in a GPS network and to improve the accuracy of the observation. Considering the non-conformance of CME in different stations, Nikolaidis (2002) added weight factor to the stations in the calculation. Although widely used, regional spatial filtering is performed in a somewhat empirical fashion, and can hardly reflect the spatial distribution of CME in the whole region. Dong et al. (2006) proposed a spatiotemporal filtering method, principal component analysis (PCA), to remove the CME. In the spatiotemporal filtering using PCA, CME can be got with varying spatial responses, providing a more solid numerical framework for analyzing the physical source in the time series. Yuan et al. (2008) found that the cycle items of CME got by PCA filtering in vertical direction decreased significantly after a load mass correction, suggesting that non-tectonic deformation is the main cause of CME in vertical direction. In addition, as an effective spatiotemporal analysis method, PCA has also been used in the field of earth geophysics for various purposes (Kawamura and Yamaoka 2006;Tiampo et al. 2004).
However, when extracting the geophysical events from the GPS coordinate time series, PCA may cause contamination between modes, i.e., a mode obtained by PCA can be contaminated by other modes (Kawamura and Yamaoka 2006). Independent component analysis (ICA) is a blind source separation method proposed in the 1990s, which can transform the observed mixed signals into a series of signals with mutually independent components in the statistical sense (Comon 1994). Huang et al. (2012) adopted ICA as a spatiotemporal filtering method to remove the CME from GPS monitoring networks, but the CME got by ICA had not been interpreted reasonably. Dai et al. (2014) used ICA as spatiotemporal modeling method in dam deformation analysis and extracted components with clear physical interpretation.
In this paper, the data of vertical time series of 80 GPS stations in the UK and 40 GPS stations in Sichuan-Yunnan region in China were used as ICA input signals. Some independent temporal components without contamination could then be extracted from these data. Then their spatial responses to each monitoring point can be obtained from the mixing matrix in the process of ICA. After that, we computed the main non-tectonic signals of each station using Quasi-Observation Combination Analysis (QOCA) software (available at http://qoca.jpl.nasa.gov/). By analyzing the independent components with spatial responses and the non-tectonic signals, we could precisely define the CME of GPS networks and make right interpretations to the common mode independent components (ICs) extracted by ICA.

Independent component analysis and spatiotemporal filtering
Basic model of ICA ICA is a useful method for blind source separation and is considered as an extension of PCA that can only make use of the second order statistical characteristics of the source (Comon 1994). Its fundamental principle is illustrated in Fig. 1. Assuming M observations X, X(t) = [X 1 (t), ⋯, X M (t)] T , from N independent components S i (t), i = 1, 2, ⋯, N, we have: Without any other priori information about matrix A or source signals, we can use ICA to obtain a separating matrix W to separate the original signals S(t) in Eq. 1, based on some optimization criteria. Generally, we can get W through two steps: (1) Whiten the observed signals X(t) by a whitening matrix B, let Z = BX and E(ZZ T ) = I (I is a unit matrix). (2) Calculate the transformation matrix by the specific independence optimize rule, let Y(t) = UZ, where Y(t) is the best approximation vector of S(t).

FastICA algorithms
ICA algorithms are usually based on the non-Gaussianity and independence of the source signals. By making full use of the high-order statistical characteristics of the source, reflecting the non-Gaussianity and independence of the source signals, ICA can effectively extract the independent components from the measured mixed signals. FastICA is a fast optimization iterative algorithm with good stability (Hyvärinen 1999;Hyvärinen and Oja 2000). In this paper, a FastICA algorithms based on the negentropy, which is a common quantitative measure of the non-Gaussianity of a random variable, is used in the independent component analysis. The stronger the non-Gaussianity of a random variable is, the greater the negentropy will be. The detailed steps are as follows: 1. Centralize and whiten the observed data. 2. Choose an initial weight vector of unit norm (random) w.
ð Þ k k : 5. Go back to step 3 if not converged. Usually, there is a spatial correlation in different stations in a GPS network, and the common deformation features may cover the internal deformation feature. So, it is necessary to remove the so-called CME using the spatiotemporal filtering. PCA can serve as an effective tool for spatiotemporal filtering and can improve the signal-to-noise ratio (SNR) of GPS time series in regional network analysis. But PCA also causes contamination between modes, and cannot interpret the CME reasonably (Kawamura and Yamaoka 2006;Dong et al. 2006). Therefore, ICA is proposed as a spatiotemporal filtering method in this paper.
Assuming the GPS displacement series are expressed as the product of a finite number of spatial and temporal components, the time series can be presented as follows: where n is the number of stations, a k represents the kth temporal components (ICs or PCs), and w jk is the spatial response of a k to x j . Effect of CME in original time sires is given as follows: where ε i is the CME in x j and R is a set of temporal components in spatiotemporal filtering. Another important question is how to identify the common mode to calculate the CME in the spatiotemporal filtering. In PCA filtering, the CME may include few top principal components (PCs), whose spatial responses are generally identical. Similarly, when using ICA, we also define the common mode through the spatial responses of the ICs.

Case study
Vertical GPS displacement series in two representative regions, the UK and Sichuan-Yunnan region in China, are applied in the spatiotemporal filtering using PCA and ICA. Two selected regions are located in different areas and have different terrain features. The UK is a coastal nation with high latitudes. Most areas in the UK have gentle terrain. Sichuan-Yunnan region, with a very complex terrain, is located in the southeast of Tibetan Plateau, and is a region prone to frequent strong earthquakes in China. In this section, spatiotemporal filtering using PCA and ICA were applied in these two regions to examine the filtering effect in different areas with different terrain features. We compared the filtering effect of PCA and ICA by analyzing the reductions of root mean square (RMS) error of time series in these two regions. Meanwhile, the physical source of CME in vertical time series was discussed by analyzing the temporal components with their spatial responses in filtering and the main non-tectonic signals with their spatial distribution.
Spatiotemporal filtering using PCA/ICA in the UK Vertical GPS displacement series provided by the services of the NERC British Isles continuous GNSS Facility (BIGF) between 2008 and 2012 were used in the spatiotemporal analysis. Daily double-differenced solutions were computed using an in-house modified version of Bernese GNSS software version 5.0. Before analyzing the time series using PCA/ICA, the GPS data time series were preprocessed by three steps. First, remove the jumps caused by a changing GPS antenna and linear term in the data series. Second, remove the stations with large gaps within the observation span 2008-2012. Finally, interpolate the data and make sure the integrity of the data.
The time series of the displacement data recorded at four representative GPS stations (ABEP, ASAP, EASN, and POOL) are shown in Fig. 2 (left), and the preprocessed ones are shown on the right. Clearly, the jumps, discontinuities and the linear velocity term are removed during the preprocessing. In this paper, time series data of 80 GPS stations were preprocessed before the spatiotemporal analysis.
The GPS vertical displacement time series were processed with PCA firstly. The top five PCs are shown in Fig. 3 on the left and their corresponding spatial responses to GPS stations are shown in the first row in Fig. 4. The PCA spatiotemporal filtering results show that the first principal component (PC1), contributing~61.6 % energy to the original data, usually represents the CME. Its corresponding spatial response values to each point have the same sign (positive or negative) and similar size. However, for a large network, the PCs with smoothed transition pattern of the spatial responses are also likely CME (Serpelloni et al. 2013), as the PC2 shows in Fig. 4. Therefore, we defined the displacements caused by the first two principal components as the CME of the GPS networks. Then, ICA was used to process the five principal components from PCA, with independent components showing in Fig. 3 on the right. Their spatial responses to each station are shown in the second row in Fig. 4. According to the FastICA algorithms, variances of all independent components are the same and equal to one. But their contributions to the original data can be reflected by their spatial response values. As Fig. 4 shows, the third and fifth independent component contribute most to the original data, and their spatial response values are of the same sign. Therefore, we defined the displacements caused by them as the CME of the GPS networks in ICA filtering. CME can be obtained using Eq. (3).
After defining the CME, we subtracted its effect to each station from the GPS vertical displacement time series. We used the root mean square (RMS) error to describe the SNR of the time series. The results in Table 1 show that the RMS value reduced by 40 and 37 % after PCA and ICA filtering, respectively. So, both these two methods can improve the SNR of residual error time series effectively.

Analysis of the result in the UK
A crustal deformation field observed by GPS usually comes with non-tectonic deformation , especially in the vertical direction. And those deformation signals should be discarded when we use the GPS observations to study tectonic deformation. Non-tectonic deformation mainly caused by the effect of tides (solid earth tides, pole tides, and ocean tides) and the earth's surface mass loading. In the data processing, tide loading (including solid tides, pole tides, and ocean tides) have been corrected. As a result, non-tectonic deformation is mainly contributed by the mass loading on the earth's surface, such as atmosphere, non-tidal ocean, snow, and soil moisture.
The deformation caused by atmospheric and soil moisture mass loading are considered as the major non-tectonic deformation signals in vertical direction. Usually, the Green's function approach (Farrell 1972) is adopted to calculate the deformation caused by the mass loading. A sub-module "mload" in the QOCA software was used in this paper to calculate the deformation caused by atmospheric and soil moisture mass loading. The input data include a station list (with the name and positions of the stations), a user-defined driver file, and the load data. We adopted the atmospheric pressure data and soil moisture data from the National Center for Environmental Prediction (NCEP). The space resolution of atmospheric pressure data is 2.5°× 2.5°and the time resolution is 6-h. The space resolution of soil moisture data is 1.875°× 1.875°and time resolution is 24-h. Since the GPS time series used in this study are daily solutions, we calculated the daily averaged surface pressure data to get the corresponding site displacements. We calculated them with QOCA software for all the 80 stations with the time resolution as 24 h. These 80 GPS stations have differences in mass loading due to their different locations, but they enjoy generally the same temporal features. The mean time series caused by atmospheric and soil moisture mass loading in the 80 stations are shown in Fig. 5.
To study the relationship between the CME and the non-tectonic deformations, we compared the ICs with the main non-tectonic deformation from spatial and temporal domain. The spatial distribution of atmospheric and soil moisture mass loading are shown in  We can see that the spatial distribution of atmospheric mass loading is quite similar to that of IC3. And the spatial distribution of soil moisture mass loading is very similar to that of IC5, which decreasing gradually from the southeast to northwest direction. The signs are opposite because of the uncertainty about the signs of ICs in the process of ICA. The CME in the UK GPS networks are mainly reflected by IC3 and IC5 according to the last section. Then, the mean CME caused by IC3 and IC5 were calculated through Eq. (3) and compared with the atmospheric and soil moisture mass loading. The results (Fig. 7) show that the CME caused by IC3 (CME3) and IC5 (CME5) and their spatial responses share similar feature and amplitude with the atmospheric mass loading and soil moisture mass loading, respectively. So, we suggest that the non-tectonic deformation, mainly including the atmospheric and soil moisture mass loading, is the main component for the CME in GPS networks, and it can be reflected by the ICs extracted by ICA spatiotemporal filtering.
Aside from IC3 and IC5, which are considered to reflect the CME, other temporal components, being statistically independent, may represent the implicit information which needs a further research. Fig. 6 Spatial distribution of the atmospheric mass loading (left) and soil moisture mass loading (right) in the UK Fig. 7 Comparision between CME (calculated with IC3 and IC5) and mass loading deformation signals (atmospheric and soil moisture) in the UK. The top one shows the CME3 and atmospheric mass loading, the bottom one shows the CME5 and soil moisture mass loading. Red lines represent the CME, black lines represent the mass loading, and blue lines show the residuals (minus 20 and 10 mm, respectively) between mass loading and CME Spatiotemporal filtering using PCA/ICA in Sichuan-Yunnan region In order to further verify the effectiveness of ICA spatiotemporal filtering, another experiment with GPS data of crustal deformation surveying network established in Sichuan-Yunnan region in China was carried out. The daily vertical time series of GNSS reference stations we used are generated at the Crustal Movement Observation Network of China (CMONOC), and processed with GAMIT/GLOBK software. The same method was used in the time series preprocess, and altogether 40 stations were chosen in this case study.
Then, PCA and ICA were used in the spatiotemporal analysis and the top five PCs and ICs that are shown in Figs. 8 and 9 show the spatial responses of PCs and ICs to each station, respectively. Considering its complex terrain, we added the height information to the figure, which may associate with the non-tectonic deformation.
In the PCA spatiotemporal filtering, the spatial response values of PC1 are with the same sign and similar size, which means that the displacement caused by the first principal component is the CME of the GPS networks. The spatial responses in the second row in Fig. 9 show that in the ICA spatiotemporal filtering, the second independent component and the third independent component contribute most to the original data, and almost all their spatial response values are of the same sign. Therefore, we defined the displacements caused by the second  and third independent components as the CME of the GPS networks in ICA filtering, and the CME can be obtained by Eq. (3).
In the PCA spatiotemporal filtering, the CME can be calculated with PC1, and it can be calculated with IC2 and IC3 in the ICA spatiotemporal filtering. The results are shown in Table 2. The RMS value reduced by 46 and 48 % after PCA and ICA filtering, respectively. Similar to the experiment in the UK, both the filtering methods can improve the SNR of residual error time series effectively.

Analysis of the result in Sichuan-Yunnan region
In the data processing in Sichuan-Yunnan Region, IERS2003 model are used in the correction for the displacements of solid earth tides and pole tides, while FES2004 model are used in the correction for the displacements from the ocean loading. So, non-tectonic deformation is also mainly contributed by the mass loading on the earth's surface. However, in China, the impact of non-tide ocean mass loading is negligibly small compared with atmospheric and soil moisture mass loading Jiang et al., 2013). Besides, the snow mass loading was also ignored as there is a little snowfall in Sichuan-Yunnan region. Thus, the effect of atmospheric and soil moisture mass loading are the major factors in the analysis of nontectonic deformation.
To investigate the physical interpretation of CME obtained by ICA spatiotemporal filtering, we calculated the deformation time series in the 40 stations caused by atmospheric and soil moisture mass loading. As the results (Fig. 10) show, these effects on the vertical components of station positions are significant.
The spatiotemporal analysis results in the UK show that the CME resulted from calculating ICs have strong correlation with the non-tectonic deformation. Similarly, we validated the conclusion of Sichuan-Yunnan region both in temporal and spatial domain. The spatial distribution of atmospheric and soil moisture mass loading in the 40 stations are shown in Fig. 11. The spatial distribution of atmospheric mass loading resembles the distribution of IC2 (expect for the northwest area with a very complex terrain), and that of soil moisture mass loading resembles the distribution of IC3. Then, we made a comparison between the mean CME resulted from calculating ICs and the two main deformation of mass loading in the 40 stations. As the results (in Fig. 12) show, the CME (CME2 and CME3) has a very good consistency with the soil moisture and atmospheric mass loading, both in amplitude and features. The spatiotemporal analysis in Sichuan-Yunnan region further verified the conclusion that ICA can extract the CME with clear physical interpretation.
It should be noted that the distribution of the mass loading does not have a very good consistency with the spatial responses of ICs in the northwest area. Dissimilar with the experiment in the UK, in which the altitude of the stations is below 500 m, the northwest zone of Sichuan-Yunnan region has a complex terrain with an altitude above 2000 m. And the altitude differences are significant in the northwest and the southeast area.
During the calculation of the non-tectonic deformation, we adopted the atmospheric pressure data with space resolution of 2.5°× 2.5°. The accuracy of atmospheric pressure data is not high enough to reflect the complexity

Conclusion
A spatiotemporal analysis method using ICA was proposed in this paper and used as a filtering method to analyze the GPS time series in vertical direction of the UK and Sichuan-Yunnan region. The results show that ICA is an effective spatiotemporal filtering method as PCA to remove to CME in GPS networks and improve the SNR of observation data.
Then, taking into account the ability of signal separation of ICA, further analysis was conducted to the CME extracted by ICA and PCA. Compared with PCA, the spatiotemporal filtering using ICA can reflect more details of the CME, and the ICs extracted by ICA have strong correlation with the deformation caused by mass loadings such as atmosphere and soil moisture. The common displacements resulted from calculating ICs have the same features and amplitude with the mass loading, though the spatial distribution did not correspond well enough, which may be caused by the inaccuracy of the global model in local areas with complex terrain. Further, the ICA spatiotemporal method may be able to provide a new way to get the mass loading more accurately in local complex area using statistical methods. Fig. 11 Spatial distribution of the atmospheric mass loading (left) and soil moisture mass loading (right) in Sichuan-Yunnan region Fig. 12 Comparison between CME (calculated with IC2 and IC3) and mass loading deformation signals (atmospheric and soil moisture) in the UK. The top one shows the CME2 and atmospheric mass loading, the bottom one shows the CME3 and soil moisture mass loading. Red lines represent the CME, black lines represent the mass loading, and blue lines show the residuals (minus 10 and 30 mm, respectively) between mass loading and CME