Singular spectrum analysis for modeling seasonal signals from GPS time series

Journal of Geodynamics 72 (2013) 25–35

Contents lists available at ScienceDirect

Journal of Geodynamics journal homepage: http://www.elsevier.com/locate/jog

Singular spectrum analysis for modeling seasonal signals from GPS time series Q. Chen a,∗ , T. van Dam b , N. Sneeuw a , X. Collilieux c , M. Weigelt b , P. Rebischung c a b c

Institute of Geodesy, University of Stuttgart, Geschwister-Scholl-Str. 24D, 70174, Stuttgart, Germany Faculté des Science, de la Technologie et de la Communication, University of Luxembourg, 6 rue Richard Coudenhove-Kalergi L-1359, Luxembourg IGN/LAREG, Université Paris Diderot, GRGS, 35 rue Hélène Brion, 75013 Paris, France

a r t i c l e

i n f o

Article history: Received 23 December 2012 Received in revised form 21 May 2013 Accepted 24 May 2013 Available online 29 June 2013 Keywords: GPS time series Modulated seasonal signals Singular spectrum analysis Least-squares fitting Kalman filtering

a b s t r a c t Seasonal signals in GPS time series are of great importance for understanding the evolution of regional mass fluctuations, i.e., ice, hydrology, and ocean mass. Conventionally these signals (quasi-annual and semi-annual signals) are modeled by least-squares fitting harmonic terms with a constant amplitude and phase. In reality, however, such seasonal signals are modulated, i.e., they will have a time-variable amplitude and phase. Recently, Davis et al. (2012) proposed a Kalman filter based approach to capture the stochastic seasonal behavior of geodetic time series. Singular Spectrum Analysis (SSA) is a non-parametric method, which uses time domain data to extract information from short and noisy time series without a priori knowledge of the dynamics affecting the time series. A prominent benefit is that trends obtained in this way are not necessarily linear. Further, true oscillations can be amplitude and phase modulated. In this work, we will assess the value of SSA for extracting time-variable seasonal signals from GPS time series. We compare our SSA-based results to those obtained using (1) least-squares analysis and (2) Kalman filtering. Our results demonstrate that SSA is a viable and complementary tool for extracting modulated oscillations from GPS time series. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction In the last three decades, GPS has overwhelmingly demonstrated its value in monitoring the deformation and/or displacement of the Earth’s surface. Innumerable studies have demonstrated the use of GPS technology for observing tectonic motion, surface mass induced displacements, post-glacial rebound, and many other geodynamically driven surface displacements. Within any one GPS time series, annual, semi-annual, decadal, long-term trends and system noise may all be superimposed. Differentiating these different components is vital for understanding the underlying geophysical or environmental processes. Much focus has been given to understanding seasonal signals that are associated with many sources. For instance, Dong et al. (2002) outlined a number of potential contributors to seasonal variations of GPS site positions. Apart from systematic errors like local multi-path, satellite orbital models or troposphere mis-modeling (Steigenberger et al., 2009), seasonal signals are known to result

∗ Corresponding author. Tel.: +49 711 68584086. E-mail addresses: [email protected] (Q. Chen), [email protected] (T. van Dam), [email protected] (N. Sneeuw), [email protected] (X. Collilieux), [email protected] (M. Weigelt), [email protected] (P. Rebischung). 0264-3707/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jog.2013.05.005

from surface-mass loading (e.g., Blewitt et al., 2001; van Dam et al., 2001, 2007; Blewitt and Lavallée, 2002; Dong et al., 2002; Davis et al., 2004). Periodic site motions from the redistribution of environmental surface masses, such as continental water variations, atmospheric pressure changes, tidal and non-tidal fluctuations in the ocean, redistribution of ice and snow, have been routinely represented by sums of sinusoids with annual and semi-annual cycles (e.g., van Dam et al., 2001). To retrieve these harmonics together with trend, a linear model is resolved by least-squares fitting prior to or simultaneous with some noise assumption (e.g., Zhang et al., 1997; Mao et al., 1999; Langbein, 2004; Williams et al., 2004). However, in this case, only constant amplitudes and phases can be obtained. True seasonal variations are not constant from year to year. Several studies have suggested determining the time-varying periodic signals by relying, for instance, on non-parametric annual signal (Freymueller, 2009; Tesmer et al., 2009), on Kalman filter based techniques (e.g., Davis et al., 2012), on piecewise continuous linear polynomials (e.g., Davis et al., 2006), or on a flexible semi-parametric model by Bennett (2008). In the geodesy, data-driven methods, such as Empirical Orthogonal Function (EOF)/Principal Component Analysis (PCA) and Independent Component Analysis (ICA), are used to smooth or model geodetic time series (e.g., Dong et al., 2006; Rangelova et al., 2007, 2010, 2012; Forootan and Kusche, 2012). Among these contributions, Dong et al. (2006) applied PCA to remove

26

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

the common mode errors from regional spatial GPS time series. Rangelova et al. (2007, 2010, 2012) used PCA and MSSA (multichannel singular spectrum analysis) respectively to investigate the main time-variable mass variations from GRACE data. In this paper, we apply singular spectrum analysis (SSA), a special case of MSSA, to investigate the possibility of extracting modulated seasonal signals from GPS time series. This paper is organized as follows: in Section 2 we introduce the detailed methodology of SSA and shortly outline the methods of least-squares fitting and Kalman filtering. This is followed in Section 3 by a short description of the GPS data we used in our research. Section 4 is the simulation part, which demonstrates the abilities of SSA in separating time-variable seasonal signals from simulated GPS time series. In this part, we also investigate how to determine the lag-window size in the process of implementing SSA. To further demonstrate its abilities, SSA derived seasonal signals are compared to the signals obtained from the conventional method (least-squares fitting) and Kalman filtering in Section 5. Finally, we discuss the advantages and disadvantages of the SSA in Section 6. 2. Methodology 2.1. Singular spectrum analysis SSA and PCA are both based on the Karhunen-Loève theory of random fields and of stationary random processes (Ghil et al., 2002). Both approaches rely on identifying the main patterns in the data in order of decreasing order of variance. The technique is implemented by diagonalizing the corresponding empirical covariance matrix. SSA is designed for univariate time series and can be extended to spatial dimension. In that case the result is Multichannel SSA (MSSA). PCA is actually a special case of the MSSA when no time lags are introduced. SSA therefore is a data-driven technique that uses time domain data to extract information from short and noisy time series without prior knowledge of the dynamics affecting the time series. An important feature of SSA is that trends obtained in this way are not necessarily linear. More importantly oscillations can be modulated in bother amplitude and phase. A number of authors provide thorough descriptions of this technique (see e.g., Broomhead and King, 1986; Vautard and Ghil, 1989; Vautard et al., 1992; Plaut and Vautard, 1994; Allen and Robertson, 1996; Ghil et al., 2002). In the following discussion, we briefly summarize their contributions and outline the SSA algorithm. The starting point of the SSA algorithm is to embed a time series (xt , 1 ≤ t ≤ N, assuming with mean subtracted) into a trajectory matrix D by sliding a M-point window. The dimension of the trajectory matrix is N × M, where N = N − M + 1. D has the form of

⎛

x1

x2

x3

·

···

xM

⎜x x x ··· · ⎜ 2 3 4 · ⎜ ⎜ x3 x4 x5 · ··· · ⎜ ⎜ ⎜ · x5 x6 x7 · · · · ⎜ D=⎜ .. .. .. .. .. ⎜ .. ⎜ . . . . . . ⎜ ⎜ · · · · · xN−1 ⎜ ⎜ ⎝ xN  · · · xN−1 xN

1 T D D N

⎛

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

c1

·

c2

···

cM−1

···

·

···

·

···

·

.. .

.. .

·

c1

c1

c0

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(3)

where entries cj , 0 ≤ j ≤ M − 1, are the covariance of x at lag j. Its unbiased estimates are: 1  xi xi+j , N−j N−j

cj =

0≤j ≤M−1

(4)

i=1

The similarities and differences between the BK and VG algorithm for computing the covariance matrix have been discussed by Allen and Smith (1996) and Ghil and Taricco (1997). The VG approach has the advantage of improved noise reduction when applied to short time series(Ghil et al., 2002). In this paper, we apply the VG approach to compute the covariance matrix C. In the second step, we apply eigenvalue decomposition to C in order to obtain the eigenvalues, k , and eigenvectors (also called EOFs), Ek , of this matrix. These are then sorted in descending order of k , where index k = 1, 2, . . ., M. The kth principal component (PC) is aki =

M 

xi+j Ejk ,

0≤i≤N−M

(5)

j=1

Finally, we can reconstruct each component of the original time series identified by the SSA using the kth reconstructed component (RC) series as given by Vautard et al. (1992)

xik = (1)

c0

⎜ c ⎜ 1 c0 c1 · ⎜ ⎜ c2 c1 c0 · ⎜ ⎜ ⎜ · c2 c1 c0 ⎜ CVG = ⎜ .. .. .. ⎜ .. ⎜ . . . . ⎜ ⎜ · · · · ⎜ ⎜ ⎝ cM−1 · · ·

⎞

⎧ i ⎪ 1 k k ⎪ ⎪ ai−j Ej ⎪ ⎪ i ⎪ ⎪ j=1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ M ⎨ 1

ak E k

1≤i≤M−1

M ≤i≤N−M+1

i−j j M ⎪ ⎪ j=1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ M ⎪  ⎪ 1 ⎪ ⎪ aki−j Ejk N − M + 2 ≤ i ≤ N. ⎪ ⎩ N−i+1

(6)

j=i−N+M

The covariance matrix is then formed by CBK =

This method of constructing the covariance matrix was originally proposed by Broomhead and King (1986), known for BK algorithm. We subsequently lable the covariance matrix as CBK . An alternative approach for computing the covariance matrix, C, is the VG algorithm (proposed by Vautard and Ghil (1989)). This algorithm is based on the lagged-covariance matrix of the process xt . With a maximum lag (or window size), M, the matrix CVG has a Toeplitz structure, i.e., constant diagonals corresponding to equal lags:

(2)

SSA typically decomposes a time series into RCs that are nearly periodic with periods less than M, whereas one or two RCs contain variations in the time series with periods greater than M. According to Plaut and Vautard (1994), harmonic oscillations can be identified in terms of the three fundamental properties of SSA: (1) two consecutive eigenvalues are nearly equal; (2) the two corresponding time

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

sequences described by EOFs are nearly periodic, with the same period and in quadrature; (3) the associated PCs are in quadrature. Sometimes large data gaps exist in a given GPS time series. For these series, one solution is to fill the gaps using interpolation techniques and apply SSA afterwards. An alternative approach was suggested by Schoellhamer (2001). He used a modified SSA algorithm to obtain spectral estimates from records with a large fraction of missing data. Another parameter f(0 ≤ f ≤ 1) is introduced that represents a specified fraction of allowable missing data points within a given window size M. When we encounter GPS time series with data gaps, we follow Schoellhamer’ s approach. We generally assign different f to different GPS time series, depending on how many gaps we have in a particular GPS time series. Big f value applies to the GPS time series with few data gaps and vice versa. For instance, Schoellhamer (2001) utilized f = 0.5 in his case [personal communication] and we normally apply f = 0.8 in our GPS time series.

27

term. When resolving a linear trend together with annual (f1 = 1cpy) and semi-annual (f2 = 2cpy) signals, the design matrix A becomes ai = [ 1 ti

sin(2ti )

cos(2ti )

sin(4ti )

cos(4ti ) ] .

(8)

The unknown vector x is x = [ y0

v a1 b1 a2 b2 ] .

(9)

The unknowns are solved for by minimizing the sum of the squares of the differences between the observations and the predictions for each epoch. For simplicity, the noise characteristics of the GPS time series are not accounted for and formal errors are used. The main advantage of the least squares is its ease of implementation and the straightforward interpretation of the estimates of the linear trend and seasonal amplitudes. Nevertheless, only constant amplitudes and/or phases are obtained. One disadvantage of the technique is that long-periodic variations can be mistakenly identified as a linear trend (Rangelova et al., 2012). In Section 5.2, we compare and contrast the SSA with least-squares technique.

2.2. Choice of lag-window size M 2.4. Kalman filtering When implementing the SSA, a key problem is choosing the lag-window size, M, that defines the spectral resolution of the algorithm. The spectral resolution can be thought of as the ability to distinguish between two spectral peaks. According to Vautard et al. (1992), M should neither be too large nor too small. If M is too small, the coarse resolution may cause several neighboring peaks in the spectrum of signal to appear as one. On the other hand, large M values will split the peak into several components with neighboring frequencies. In addition, large M will weaken the ability of the algorithm to isolate peaks at frequencies lower than the resolution, 1/M, or equivalently periods larger than M. Vautard et al. (1992) therefore concluded that SSA is able to resolve oscillations with periods between M/5 and M. Moreover, Ghil et al. (2002) also proposed that when choosing the window-lag time one should balance two considerations: the quantity of information extracted versus the degree of statistical confidence in that information. The former consideration requires that the window lag should be as wide as possible, i.e. a large M. While the latter factor, requires as many repetitions of the features of interest as possible, i.e., as large a ratio N /M as possible. In the geodetic field, Rangelova et al. (2012) applied a three-year window (156 weeks) to retrieve long-term changes together with seasonal variations from weekly GRACE solutions. For the case of weekly GPS time series, annual and semi-annual signals are currently of most interests for us. As a result, a window size of two or three years should be appropriate. Section 4.3 confirms this choice. Recently, a mathematic way of window length selection which is based on statistical test (test the convergence of the autocovariance function) is proposed by Khan and Poskitt (2012). As our main interets focus on extracting annual and semi-annual signals, we therefore only follow the empirical rules rather than try the newly proposed mathematic method to determine the window length. 2.3. Least-squares fitting Least-squares (LS) fitting is a conventional method that is often applied to weekly GPS time series to simultaneously fit linear and periodic terms as well as periodic terms, e.g. y(ti ) = y0 + vt i +

q 

ak sin(2fk ti ) + bk cos(2fk ti ) + (ti ),

(7)

k=1

where ti is the epoch of observation i in decimal year, y0 is a constant initial offset, v is a constant velocity, ak and bk are the coefficients of periodic terms, fk are the periodic frequencies and (t) is the noise

Davis et al. (2012) developed an expression for the time-varying seasonal signal in GPS time series as the sum of a deterministic component with an annual period and a zero-mean stochastic component. In Davis et al. (2012), the authors implemented a model (see Eq. (10)) that included stochastic annual and semiannual terms using linear Kalman filtering. x(t) = x0 + v(t)(t − t0 ) + a1 (t) cos(2f0 (t − t0 )) + b1 (t) sin(2f0 (t − t0 )) + a2 (t) cos(4f0 (t − t0 )) + b2 (t) sin(4f0 (t − t0 )),

(10)

where f0 = 1cpy and t0 is a reference epoch, v(t) is time-variable velocity term and a1 (t), a2 (t), b1 (t) and b2 (t) are instantaneous time-variable amplitudes. We have implemented the Davis et al. (2012)’s Kalman filtering scheme to compare with SSA (see Section 5.2). We used the same state process (random constant for x0 and random walk for stochastic terms) and the same input variance (variance rate values of 1 mm2 yr−3 for the rate term and 0.5 mm2 yr−1 for sinusoidal amplitudes). Nevertheless, those stochastic parameters can be also modeled as other stochastic process, such as First Order Gaussian Markov (FOGM) process (e.g., Reilinger et al., 2006; Ji and Herring, 2011). 3. Data We used GPS combined weekly station position time series computed by the International GNSS Service (IGS) (Dow et al., 2009) from 1998.0 to 2011.3. These data include homogeneously reprocessed coordinates from the first IGS reprocessing campaign up to GPS week 1459 (December 2007). However, the reprocessed individual solutions each IGS Analysis Center have been recombined using the combination strategy of the new IGS combination center (Rebischung et al., 2012). For the period considered here, the origin and scale information from these GPS solutions were preserved and not projected into the International Terrestrial Reference Frame (ITRF) realizations as compared with previous combinations (Ferland and Piraszewski, 2009). We thus followed the approach of (Collilieux et al., 2012) to derive station position time series consistent with ITRF orientation but deprived of GPS apparent geocenter motion. Their approach was specifically developed in order to preserve crustal deformation information in the position time series, especially at the seasonal timescale. We first estimated reference station positions and

28

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

Amplitude/mm

15

(a)

(b)

Simulated data Deterministic part Linear trend Annual part Semi−annual part

Simulated annual signal constant part

2

deviation part

1

mm

20

0 −1

10

−2 Simulated semi−annual signal 5

constant part

2

deviation part

mm

1 0

0 −1

−5

−2 2000

2002

2004 Time/yr

2006

2000

2008

2002

2004

2006

2008

Time/yr

Fig. 1. (a) Simulated 641 weeks GPS height-coordinate time series consisting of (1) a linear trend, (2) annual and semi-annual cycles with constant amplitudes, (3) annual and semi-annual cycles with time-variable amplitudes, and (4) noise that includes white noise, flicker noise and random walk noise; (b) detailed description of simulated annual and semi-annual parts.

velocities as well as position offsets, in case of position discontinuities, for every station. We then estimated translation and rotation parameter time series between input weekly network station positions and these reference coordinates using the IGS08 core subset of stations (Rebischung et al., 2012). These transformation parameters were then applied to the original position time series in order to remove GPS apparent geocenter motion and possible orientation mis-alignment. The output positions were then corrected for the long-term positions and velocities to extract the residual non-linear coordinate variations.

where ti is the epoch of observation i in decimal year, y0 is a constant initial offset, v is a constant velocity and (ti ) represents the noise. Here we follow the suggestions of previous studies (e.g., Zhang et al., 1997; Mao et al., 1999) and adopt a noise model consisting of white noise, flicker noise, and random walk noise. However, we assign different variances to those three noise components because white noise and flicker noise are prevalent in GPS time series and random walk plays a secondary role. s(ti ) is our model for time-variable station motion. In contrast to Bennett (2008), a time-variable semiannual cycle is also added. Therefore s(ti ) is in the form of: s(ti )

4. Simulation

= a sin(2ti ) + b cos(2ti ) + c(ti ) sin(2ti ) + c(ti ) cos(2ti )

(12)

+d sin(4ti ) + e cos(4ti ) + f (ti ) sin(4ti ) + f (ti ) cos(4ti ),

As described before, SSA has been used widely and successfully in the field of climate studies. To demonstrate the application of SSA in the field of GPS, we simulate 641 epochs of weekly GPS time series, where the true signal is amplitude modulated. The intention of the simulation is to investigate the performance of SSA under the control of signal and noise composition. In reality, most GPS stations show characteristics with trend and periodic signals plus noise model. To make the simulation simple and reasonable, we use a kinematic model for the time-evolution of a GPS continuous coordinate time series y(ti ) given by Bennett (2008) y(ti ) = y0 + vt i + s(ti ) + (ti ),

(11)

where a, b, d and e are constant time-averaged wave amplitudes and c(ti ), f(ti ) are amplitude variations having the form of Eq. (13).

⎧ 0.3 sin(ti ) ⎪ ⎨ c(ti ) = e i

The reader should note that this option for the deviation function might not be more reasonable than any other potential options. The simulated data is shown in Fig. 1(a) and the detailed simulated annual and semi-annual parts are depicted in Fig. 1(b).

0.3

30%

(b)

0.25

25% Variance(in percent)

Normalized Eigenvalues

(a)

0.2 0.15 0.1 0.05 0 0

(13)

⎪ ⎩ f (t ) = e0.3 sin(ti ) .

20% 15% 10% 5%

20

40 60 80 Eigenvalue indices

100

0%

1

2 4 Frequency/cpy

10

15 20

30

Fig. 2. (a) Eigenvalues from eigenvalue decomposition of the covariance matrix from simulated data; (b) eigenvalues versus the dominant frequency associated with their corresponding EOFs.

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

0.2

(a)

EOF1

EOF2

0.1

20

0

0

−0.1 −0.2 0.2

(c)

PC1

PC2

(d)

PC3

PC4

−20

(b)

EOF3

20

EOF4

0.1

10

0

0

−0.1

−10

−0.2 0

29

−20

20

40

60 Lag/weeks

80

100

120

2000

2002

2004 Time/yr

2006

2008

Fig. 3. Left: the first four empirical orthogonal functions (EOFs) from simulated data are grouped into two pairs, (1, 2) and (3, 4), representing annual and semi-annual periods respectively; right: the first four principal components (PCs) from the simulated data are grouped into two pairs as well.

In the following section, two cases will be discussed. The first case is that the simulated data is composed of periodic signals, i.e., we have omitted the linear trend (vti ). In the other case, we investigate the ability of SSA in the presence of a strong linear trend.

Annual signal 5

Original

Reconstructed

0 4.1. Case I: simulated data without a linear trend

−5 Semi−annual signal 4

Original

Reconstructed

2 0 −2 −4 2000

2002

2004

2006 Time/yr

2008

2010

Fig. 4. Reconstructed annual signal and semi-annual signal from the first four EOFs and PCs compared to the original annual and semi-annual signal.

6

True signal

Least−squares

SSA

4 Amplitude/mm

In this section, we analyze the simulated data without a linear trend. Following the methodology presented in Section 2.1, the SSA decomposes the simulated data into 105 temporal modes using a two-year lag-window, M = 105 weeks (the choice of lag-window M will be discussed in Section 4.3). This result is shown in Fig. 2(a). The eigenvalues are normalized in such a way that they represent fractions of the total data variance. Fig. 2(a) clearly shows the first two eigen-pairs that explain 53.95% and 31.21% of the data variance respectively and demonstrating the first property of SSA (two consecutive eigenvalues are nearly equal). Fig. 3(a–d) depicts the first four EOFs and PCs, which are in quadrature and perfectly fulfill the other two properties (corresponding EOFs and PCs are periodic and in quadrature). From the spectral point of view (Vautard et al., 1992; Allen and Smith, 1996), an alternative way to display the information is to plot the eigenvalues against the dominant frequency associated with their corresponding EOFs. This result is shown in Fig. 2(b). Fig. 2(b) depicts the prominent annual and semiannual cycles. According to the criteria of Plaut and Vautard (1994), the first four components are reconstructed and grouped into two components, annual signal and semi-annual signal, shown in Fig. 4. The reconstructed annual signal and semi-annual signal fit the original simulated time-variable annual cycles and semi-annual cycles quite well. It demonstrates that SSA has the ability to capture the time-variable component of periodic signals. A combined signal (annual plus semi-annual) is shown in Fig. 5. For comparison, we also plot a periodic signal extracted by least-squares. It is evident that least-squares fails to capture the time-variable peaks. RMS (Root-Mean-Square) of the difference between the true seasonal signal and the SSA-derived signal, LS-derived signal are calculated, which are 0.21 mm and 0.32 mm respectively. Both methods can achieve quite small RMS. However, SSA performs better than LS. This case study indicates that SSA is able to distinguish the amplitude variations that really exist in the simulated data and extract more seasonal signals.

2

0

−2

−4 2000

2002

2004

2006 Time/yr

2008

2010

Fig. 5. Seasonal signals from SSA and least squares are compared with the true simulated seasonal signals. It clearly shows least-squares fitting fails to capture the time variable part but SSA does.

30

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

EOF1

0.12

PC1

200

0.1

100

0.08 0.2

EOF2

Aliased!

0

0 50

PC2

0

−0.2 0.2

−50 50

EOF3

0

PC3

0

−0.2 0.2

EOF4

Aliased!

0

−50 50

PC4

0

−0.2 0.2

−50 50

EOF5

0

PC5

0

−0.2 0.2

EOF6

Aliased!

0

−50 50

PC6

0

−0.2 0

50 Lag/weeks

100

−50

2000 2002 2004 2006 2008 Time/yr

Fig. 6. When a linear trend exists in the time series, it will alias into EOF2, EOF4 and EOF6.

4.2. Case II: data with a linear trend

4.3. Simulation on determining lag-window M Recall from Section 2.2 that the lag-window size is important in the process of SSA. In this part, we will show a simple demonstration on determing the lag-window M. We use the simulation data from Case I again here. Lag-window sizes of 1-year (M = 53 weeks), 2-year (M = 105 weeks), 3-year (M = 157 weeks), 4-year (M = 209 weeks), 5-year (M = 261 weeks) are selected. Fig. 7 shows the eigenvalues (represented in variance) as a function of the dominant periodic signal in the EOFs and in terms of different lag-window sizes. According to Fig. 7, annual and semi-annual cycles are evidently obtained with all five different lag-window sizes. However,

25%

Variance(in percent)

In this section, we add the linear trend back into the data (see Fig. 1). Analogous to Case I, a window size of 105 weeks is applied here. Fig. 6 gives the first six EOFs and PCs, which represent longterm variation (EOF1&PC1), annual cycle (EOF2, EOF3&PC2, PC3), and semi-annual cycle (EOF4, EOF5, EOF6&PC4, PC5, PC6) respectively. Among these, spectral mixing can be observed in EOF2, EOF4 and EOF6 (they do not purely vary around zero). This aliasing is probably due to the presence of the strong linear trend and colored noise (Allen and Robertson, 1996; Rangelova et al., 2010, 2012). In this situation, the trend adversely affects the periodic modes. To handle this problem, a simple solution is to remove the linear trend using a least-squares estimation. Rangelova et al. (2010, 2012) suggested that the removal of the GIA signal largely reduces the variance of the long-term modes leaving the annual cycle as the dominant signal when analyzing the GRACE-derived water mass variability with MSSA (Multi-Channel SSA). In the case of continuous GPS time series, strong trends also often dominate the motion of GPS stations in most areas. Before applying the SSA technique, we advise detrending the data so that better results can be obtained. When the linear trend is removed, Case II is equivalent to Case I.

M=53 M=105 M=157 M=209 M=261

20%

15%

10%

5%

0%

1

2 Cycles/yr

5

10

15 20 2530

Fig. 7. Eigenspectrum for five different lag-window sizes.

more variance is reconstructed with smaller lag-windows, whereas larger lag-windows split part of main variance into other artificial long-term variations. Table 1 presents a comparison of the variances from the true signal and the SSA extracted signal. In the simulated data, the annual and semi-annual signal make up 53.42% Table 1 Variance comporison of five different lag-window sizes. True variance

Annual Semi-annual

53.42% 30.24%

M (weeks) 53

105

157

209

261

56.01% 32.57%

53.95% 31.21%

52.13% 30.09%

50.36% 29.09%

48.72% 28.16%

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

0.2 (a) 0 −0.2

EOF 1

EOF 2

0.2 (b) 0 −0.2

EOF 3

EOF 4

0.2 (c) 0 −0.2

EOF 5

50 (e) 0 −50 50 (f)

PC 1

PC 2

PC 3

PC 4

PC 5

PC 6

PC 7

PC 8

0 −50

EOF 6

50 (g) 0 −50

0.2 (d)

EOF 7

EOF 8

50 (h) 0

0 −0.2 0

31

20

40

60

80

100

120

140

160

−50 1998

1999

2000

2001 2002 Time/yr

Lag/weeks

2003

2004

2005

Fig. 8. Left: first eight EOFs of the height-coordinate time series from the WIS1 station; These EOF’s represent an annual cycle, a long-term variation, a quarterly cycle and a semi-annual cycle; right: first eight corresponding PCs of the height coordinate of WIS1 station.

and 30.24% of the total variance respectively. With the 1-year lagwindow (53 weeks), the annual and semi-annual signal make up 56.01% and 32.57% of the true variance. The 1-year lag then overestimates the annual and semi-annual signals in the time series. With the lag-window size of 2-years (105 weeks) and 3-years (157 weeks), the variance from SSA reconstruction varies slightly around the true variance. Large under-estimation of the variance happens with larger lag-window sizes, e.g. 4 years and 5 years. Additionally, unrealistic long-term changes can be observed with the lag-windows set larger than 2 years. To a certain extent, this analysis confirms the findings of Vautard et al. (1992) that large M values will split the peak into several components with neighboring frequencies and weaken the capability of SSA to isolate peaks at frequencies lower than the resolution 1/M. Based on this simulation, we determine that 2-year or 3-year lag-windows should be appropriate for extracting the annual and semi-annual periods from weekly GPS time series. However, if long-term variations are of interest, we recommend that larger lag-window sizes to be used (Rangelova et al., 2012). In the following real data analysis part, a 2-year or 3-year lag-window will be applied.

Table 2 Variance contribution of each signal in the height-coordinate time series of WIS1 station as a percent of the total variance.

Annual Long-term Semi-annual Total variance

Modes

Variance in % of the total variance

1, 2 3, 4 7, 8

17.49, 17.30 4.50, 3.91 1.89, 1.87 46.96

which explain 46.96% of the total variance. This is close to half of the energy within the whole signal. Fig. 9 presents the EOFs spectrum against eigenvalues (represented in variance). It is obvious that the the annual cycle stands out in the whole spectrum. This is probably due to seasonal water-storage effects around Great Lakes. Besides, it is interesting to find that semiannual and quarterly periods, 2.09 cycles/yr and 4.13 cycles/yr respectively, probably correspond to the second (2.08 cycles/yr)and fourth draconitic year (4.16 cycles/yr) (Collilieux et al., 2007; Ray et al., 2008). Fig. 10 displays the original time series and the reconstructed series (including a long-term variation, annual and semi-annual

5. Real data analysis

20%

10%

4.16 cycles line

15%

2.08 cycles line

1.04 cycles line

To assess the performance of SSA in the case of real data, we select the weekly GPS height-coordinate time series from WIS1, Wisconsin Point (USA). WIS1 station is located in the Great Lakes ◦ ◦ region at a latitude of 46 . 7 N and 92 . 0 W. The whole time series spans almost 10 years, from 1998.021 to 2007.814. Within the time series, an annual cycle dominates the signal. The original time series is presented in Fig. 10. Prior to building the covariance matrix C in Eq. 3, the GPS time series is centered and normalized by its standard deviation. Based on the simulation in Section 4.3, a lag-window of 157 weeks is chosen. Following the steps in Section 2.1, Fig. 8 presents the decomposed EOFs and PCs. Periodic pairs and long-term variation pairs are clearly observed. Fig. 8(a) and (e) describe the annual behavior that does not show much variations with time. Fig. 8(b) and (f) show the long-term variations. Fig. 8(c) and (g) correspond to quarterly signal pairs. Semi-annual pairs are displayed in Fig. 8(d) and (h). In terms of the three criteria of Plaut and Vautard (1994), those pairs are grouped and reconstructed (the figure is not shown here). Table 2 shows the reconstructed variance of the main signals,

Variance(in percent)

5.1. Extracting seasonal signals

5%

0%

1

2 4 Cycles/yr

10

20 30

Fig. 9. Variance contribution versus the dominant frequency associated with corresponding EOFs. The dominant frequencies in the semi-annual and quarterly cycle are close to the second and fourth multiplies of the frequency of the draconitic year.

32

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

15 Original time series

Long−term+Annual+Semi−annual

10

5

Table 3 Correlation coefficients of original time series, least-squares filtered seasonal signal, SSA filtered seasonal signal and Kalman filtering filtered seasonal signal for the height-coordinate of WIS1.

Data LS SSA KF

mm

0

Data

LS

SSA

KF

1 0.57 0.63 0.63

1 0.94 0.96

1 0.98

1

−5

−10

−15

−20 1998

2000

2002

2004

2006

2008

Time/yr Fig. 10. Original height-coordinate time series of WIS1 and the SSA extracted seasonal signals.

cycles). The filtered signal fits the original time series very well, capturing all the peaks of the data. 5.2. Comparison with least-squares fitting and Kalman fitering In the previous sections, we have demonstrated that SSA can be used to extract the time-variable periods from GPS time series. To assess the performances further, in this section we compare SSA with conventional least-squares fitting and Kalman filtering by Davis et al. (2012) to extract a trend plus periodic oscillations from GPS time series. We discuss two examples and only seasonal signals are considered. For the first example, periodic signals dominate the time series. In the other example, the periodic signals are small and non-linear signals are significant. The height-coordinate time series from WIS1, with a strong annual signal, is used again. A lag-window of 157 weeks is chosen for the SSA analysis. Settings for LS fitting and Kalman filtering are described in the methodology section (see Section 2.3 and 2.4). Fig. 11(a) shows the annual and semi-annual cycles derived using LS fitting, SSA and Kalman filtering respectively, together with the original time series. In the LS fitting series, the seasonal signal peaks during the beginning of September and reaches a minimum during the beginning of April every year. The seasonal signals retrieved from SSA and Kalman filtering, however, vary with time. We do not discuss further whether those peak shifts are geophysically meaningful or not because this requires further analysis of geophysical cause of the motion of the station. We only focus on the methodology and demonstrate that SSA has the ability to capture those shifted peaks. In addition, we observe that the time-vaying periodic signals obtained from SSA and Kalman filtering fit each other quite well (Fig. 11(a)). The RMS of difference between SSA and Kalman filtering derived signals is 0.64 mm which is smaller than that between SSA and LS (1.1 mm), and that between Kalman filtering and LS (0.86 mm). We do not observe much difference in terms of RMS because this GPS time series is dominating by seasonal signals. We calculate the correlation between the original time series and the (1) least-squares fit, (2) SSA and (3) Kalman filter derived series. Table 3 displays the result. With respect to the original time series, a high correlation is observed between the SSA and Kalman filtering series. The LS series has the lowest correlation. We can conclude that SSA and Kalman filtering perform similarly and that both are superior to LS fitting.

In the second example we use the east-coordinate GPS time series from MIL1. MIL1 is located on the west coast of Michigan ◦ ◦ Lake with coordinates 43 . 0 N and 87 . 9 W. The time series is 9 years long. A large west-ward motion can be observed from 1998 to 2002. This trend is probably caused by a long-term decrease in the water level of the Great Lakes. The station receiver was changed during this period in June 1999 but no obvious offset occurs at that time. In contrast to the first example, a window size of 105 weeks is chosen. Fig. 11(b) presents the comparison between the filtered periodic signals derived from LS fitting, SSA and Kalman filtering. It should be pointed out that SSA is able to resolve non-linear longterm variations buried in the GPS time series (see curve in magenta in Fig. 11(b)). Nevertheless, as our main interests focus on annual and semi-annual periods, we do not include the long-term variation in this comparison. Because a non-linear long-term variation dominates the time series, LS fitting only extracts relatively weak periodic information from the data. However, SSA has the ability to retrieve the nonlinear information and considerable annual and semi-annual cycles are identified in the presence of other strong signals. Kalman filtering can retrieve the time-varying annual and semi-annual cycles as well. As in the previous example, the Kalman signal fits the SSA very well, although some disagreements appear at the beginning of time series that might be due to a problem with the convergence of the Kalman filter. The RMS of difference between the seasonal signals received by SSA and Kalman filtering is 0.30 mm which is only half of that between SSA and LS (0.62 mm) and that between Kalman filtering and LS (0.58 mm). Significant differences are received in terms of the RMS of difference when the GPS time series is not dominating by periodic signals. Moreover, it is interesting to note that SSA and Kalman signals peak after the LS signals in the first half of the time series (up until 2002). After 2002 and until the end of the time-series, the SSA and Kalman signals peak earlier than the LS signal. This effect might be caused by local environment effects, e.g., hydrology. Power spectral densities (PSD) (calculated using the LombScargle algorithm, see Scargle (1982), Press et al. (2001)) of the original time series and the residuals of the various techniques are shown in Fig. 12(a) for WIS1 and Fig. 12(b) for MIL1. In the case of WIS1, the power of the annual cycles are reduced with respect to the original power by all three approaches. However, as the semiannual period deviates from exactly 2 cycles per year, LS fitting fails to capture the shift and power remains in the LS residuals around the semi-annual cycles. On the contrary, SSA and Kalman filtering both succeed in fitting the semi-annual peaks and are able to considerably reduce the power around the semi-annual period. In the case of MIL1, LS fitting does not reduce the power around annual and semi-annual cycles because they both deviate from their exact frequencies. However, SSA and Kalman filtering can both approximate the original signal around the frequencies of interest and can reduce the power at those frequencies accordingly. We repeat the correlation analysis as was used in the first example. Table 4 presents the result. The correlation coefficient between SSA filtered seasonal signal and the Kalman filtering filtered signal is 0.92. Both SSA and Kalman filtering have a very low correlation

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

20

5

(a)

Original time series LS: Annual+Semi−annual KF: Annual+Semi−annual SSA: Annual+Semi−annual

15 10

33

(b)

Original time series LS: Annual+Semi−annual KF: Annual+Semi−annual SSA: Annual+Semi−annual SSA: Long−term variation

4 3 2 1

mm

mm

5 0

0 −5

−1

−10

−2

−15

−3

−20 1998

2000

2002

2004

2006

−4

2008

2000

Time/yr

2002 Time/yr

2004

2006

Fig. 11. Seasonal signals estimated by Least-squares fitting, SSA and Kalman filtering for the height-coordinate of WIS1 (a) and the east-coordinate of MIL1 (b).

(a)

Original time series Residuals after LS

2

(b)

Residuals after KF Residuals after SSA

10

1

1

10 Power

Power

Residuals after KF Residuals after SSA

10

10

0

10

−1

0

10

−1

10

10

−2

−2

10

10

−3

10

Original time series Residuals after LS

2

−3

0.1

1 2 Frequency/cpy

4

10

20

10

0.1

1 2 Frequency/cpy

4

10

20

Fig. 12. Power spectral density of residuals obtained by Least-squares fitting, SSA and Kalman filtering for the height-coordinate of WIS1 (a) and for the east-coordinate of MIL1 (b). Residuals here refer to original time series less the seasonal signals derived from the three techniques.

with LS fitting demonstrating that LS fitting fails to retrieve the annual and semi-annual signals in the presence of other strong signals. The SSA filtered series correlates to the original time series slightly better than the Kalman filtered series. This shows that SSA performs slightly better than Kalman filtering in the presence of long-period signals. The above two comparisons between LS fitting, SSA and Kalman filtering shows SSA and Kalman filtering outperform LS fitting in both of the cases analyzed. SSA and Kalman filtering are superior to least-squares in their ability to capture signals with modulated amplitudes and phases. LS fitting can not fit the GPS time series very well when non-linear signals dominate the data. SSA Table 4 Correlation coefficients of the original time series, least-squares filtered seasonal signal, SSA filtered seasonal signal and Kalman filtering seasonal signal for the eastcoordinate of MIL1.

Data LS SSA KF

Data

LS

SSA

KF

1 0.26 0.58 0.55

1 0.57 0.56

1 0.92

1

produces the similar results as Kalman filtering in the both examples. As for which method is better, further comparisons and analysis are required. 6. Discussion and conclusion Investigating time variable seasonal signals associated with time varying amplitudes and phases for GPS time series is challenging. Some promising research has already been carried out, e.g., Davis et al. (2012). The main objective of our study is to try to address the problem in an alternative way and to extract the modulated periodic cycles from the original GPS time series. From the simulation and real data analysis, we demonstrate that singular spectrum analysis (SSA), as a data-driven method, has the ability to extract amplitude and phase varying periodic variabilities from GPS coordinate time series. Bennett (2008) and Davis et al. (2012) proposed model based approaches to investigate the time-varying signals in GPS time series. These model based methods are implemented under certain assumptions. For example, Bennett’ s semi-parametric assumed time-variable component changes smoothly with time and approximately assumed the phase  is constant to reduce the computation

34

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35

burden. In Davis’s approach, a random walk process was assumed as the state model in the course of Kalman filtering. In this case, the noise process (input variance) should be assumed in advance or should be estimated (Collilieux et al., 2007), which significantly increases the total computational cost. SSA, as a data-driven method, can avoid the complication of assuming the noise processes that affect the time series. In SSA, the only variable parameter is the lag-window size M, which depends on the length of time series and the periodic cycles of interest. In the simulation, we demonstrated that a two-year or three-year lag-window should be appropriate for extracting time-varying annual and semi-annual cycles. If other signal having different frequencies are of interest, a different lag-window size needs to be assigned. However, there are some drawbacks to SSA. Like other datadriven methods, such as EOF/PCA, no error propagation is carried out during the analysis. This means that we cannot obtain a quantitative evaluation of the analysis. A significance test is required after the analysis, especially in the presence of colored noise. Allen and Smith (1996) proposed to adopt Monte Carlo SSA to distinguish a real geophysical signal from an artificial signal driven by colored noise. In addition, to a certain extent, SSA is not a globally efficient analysis tool with respect to those model based methods because analysis needs to be performed station by station and component by component. One potential solution to this drawback is to apply multi-channel singular spectrum analysis (MSSA) to time series in a regional or global GPS network. In conclusion, we have introduced an alternative approach to extract time-variable seasonal signals from GPS time series. SSA is a viable and complementary tool for exploring modulated oscillations from GPS time series. Acknowledgements We would like to thank two anonymous reviewers for helpful comments and suggestions on the manuscript. Q. Chen is supported by China Scholarship Council (CSC) for his PhD study. References Allen, M., Smith, L., 1996. Monte Carlo SSA: detecting irregular oscillations in the presence of colored noise. J. Climate 9, 3373–3404, doi: 10.1175/15200442(1996)009<3373:MCSDIO>2.0.CO;2. Allen, M.R., Robertson, A.W., 1996. Distinguishing modulated oscillations from coloured noise in multivariate datasets. Climate Dyn. 12, 775–784, doi: 10.1007/s003820050142. Bennett, R.A., 2008. Instantaneous deformation from continuous GPS: contributions from quasi-periodic loads. Geophys. J. Int. 174, 1052–1064, doi: 10.1111/j.1365246X.2008.03846.x. Blewitt, G., Lavallée, D., 2002. Effect of annual signals on geodetic velocity. J. Geophys. Res. 107, 2145, doi: 10.1029/2001JB000570. Blewitt, G., Lavallée, D., Clarke, P., Nurutdinov, K., 2001. A new global mode of earth deformation: seasonal cycle detected. Science 294, 2342–2345, doi: 10.1126/science.1065328. Broomhead, D., King, G.P., 1986. Extracting qualitative dynamics from experimental data. Phys. D: Nonlinear Phenom. 20, 217–236, doi: 10.1016/01672789(86)90031-X. Collilieux, X., Altamimi, Z., Coulot, D., Ray, J., Sillard, P., 2007. Comparison of very long baseline interferometry, GPS, and satellite laser ranging height residuals from ITRF2005 using spectral and correlation methods. J. Geophys. Res. 112, B12403, doi: 10.1029/2007JB004933. Collilieux, X., van Dam, T., Ray, J., Coulot, D., Métivier, L., Altamimi, Z., 2012. Strategies to mitigate aliasing of loading signals while estimating GPS frame parameters. J. Geodesy, 1–14, doi: 10.1007/s00190-011-0487-6. van Dam, T., Wahr, J., Lavallée, D., 2007. A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe. J. Geophys. Res. 112, B03404, doi: 10.1029/2006JB004335. van Dam, T., Wahr, J., Milly, P.C.D., Shmakin, A.B., Blewitt, G., Lavallée, D., Larson, K.M., 2001. Crustal displacements due to continental water loading. Geophys. Res. Lett. 28, 651–654, doi: 10.1029/2000GL012120. Davis, J.L., Elósegui, P., Mitrovica, J.X., Tamisiea, M.E., 2004. Climate-driven deformation of the solid earth from GRACE and GPS. Geophys. Res. Lett 31, L24605, doi: 10.1029/2004GL021435.

Davis, J.L., Wernicke, B.P., Bisnath, S., Niemi, N.A., Elósegui, P., 2006. Subcontinentalscale crustal velocity changes along the Pacific-North America plate boundary. Nature 441, 1131–1134, doi: 10.1038/nature04781. Davis, J.L., Wernicke, B.P., Tamisiea, M.E., 2012. On seasonal signals in geodetic time series. J. Geophys. Res. 117, B01403, doi: 10.1029/2011JB008690. Dong, D., Fang, P., Bock, Y., Cheng, M.K., Miyazaki, S., 2002. Anatomy of apparent seasonal variations from GPS-derived site position time series. J. Geophys. Res. 107, 2075, doi: 10.1029/2001JB000573. Dong, D., Fang, P., Bock, Y., Webb, F., Prawirodirdjo, L., Kedar, S., Jamason, P., 2006. Spatiotemporal filtering using principal component analysis and KarhunenLoeve expansion approaches for regional GPS network analysis. J. Geophys. Res. 111, B03405, doi: 10.1029/2005JB003806. Dow, J., Neilan, R., Rizos, C., 2009. The international GNSS Service in a changing landscape of Global Navigation Satellite Systems. J. Geodesy 83, 191–198, doi: 10.1007/s00190-008-0300-3. Ferland, R., Piraszewski, M., 2009. The IGS-combined station coordinates, earth rotation parameters and apparent geocenter. J. Geodesy 83, 385–392, doi: 10.1007/s00190-008-0295-9. Forootan, E., Kusche, J., 2012. Separation of global time-variable gravity signals into maximally independent components. J. Geodesy 86, 477–497, doi: 10. 1007/s00190-011-0532-5. Freymueller, J., 2009. Seasonal position variations and regional reference frame realization. In: Drewes, H. (Ed.), Geodetic Reference Frames. Springer Berlin Heidelberg. volume 134 of International Association of Geodesy Symposia, pp. 191–196, doi: 10.1007/978-3-642-00860-3 30. Ghil, M., Allen, M.R., Dettinger, M.D., Ide, K., Kondrashov, D., Mann, M.E., Robertson, A.W., Saunders, A., Tian, Y., Varadi, F., Yiou, P., 2002. Advanced spectral methods for climatic time series. Rev. Geophys. 40, 1003, doi: 10.1029/ 2000RG000092. Ghil, M., Taricco, C., 1997. Advanced spectral analysis methods. In: Castagnoli, G., Provenzale, A. (Eds.), Past and Present Variability of Solar-Terrestrial System: Measurement, Data Analysis and Theoretical Models, Soc. Ital. di Fis. Bologna, Italy, pp. 137–159. Ji, K.H., Herring, T.A., 2011. Transient signal detection using GPS measurements: Transient inflation at Akutan volcano, Alaska, during early 2008. Geophys. Res. Lett. 38, L06307, doi: 10.1029/2011GL046904. Khan, M.A.R., Poskitt, D.S., 2012. Moment tests for window length selection in singular spectrum analysis of short- and long-memory processes. J. Time Ser. Anal., doi: 10.1111/j. 1467-9892.2012.00820.x. Langbein, J., 2004. Noise in two-color electronic distance meter measurements revisited. J. Geophys. Res. 109, B04406, doi: 10.1029/2003JB002819. Mao, A., Harrison, C.G.A., Dixon, T.H., 1999. Noise in GPS coordinate time series. J. Geophys. Res. 104, 2797–2816, doi: 10.1029/1998JB900033. Plaut, G., Vautard, R., 1994. Spells of low-frequency oscillations and weather regimes in the northern hemisphere. J. Atmos. Sci. 51, 210–236, doi: 10.1175/15200469(1994)051<0210:SOLFOA>2.0.CO;2. Press, W., Teukolsky, S., Vetterling, W., Flannery, B., 2001. Numerical Recipes in Fortran 77:the art of scientific computing (2nd edn), Vol. 1 of Fortran Numerical Recipes. Cambridge University Press, Cambridge. Rangelova, E., Sideris, M., Kim, J., 2012. On the capabilities of the multichannel singular spectrum method for extracting the main periodic and non-periodic variability from weekly GRACE data. J. Geodyn. 54, 64–78, doi: 10.1016/j.jog.2011.10.006. Rangelova, E., van der Wal, W., Braun, A., Sideris, M.G., Wu, P., 2007. Analysis of Gravity Recovery and Climate Experiment time-variable mass redistribution signals over north america by means of principal component analysis. J. Geophys. Res. 112, F03002, doi: 10.1029/2006JF000615. Rangelova, E., Wal, W., Sideris, M.G., Wu, P., 2010. Spatiotemporal analysis of the GRACE-derived mass variations in North America by means of multi-channel singular spectrum analysis: Gravity. Geoid Earth Observ. 135, 539–546, doi: 10.1007/978-3-642-10634-7 72. Ray, J., Altamimi, Z., Collilieux, X., van Dam, T., 2008. Anomalous harmonics in the spectra of GPS position estimates. GPS Solut. 12, 55–64, doi: 10. 1007/s10291007-0067-7. Rebischung, P., Garayt, B., Collilieux, X., Altamimi, Z., 2012. IGS Reference Frame Working Group Coordinator Report 2011. IGS Technical Reports 2011. Astronomical Institute of the University of Bern. Rebischung, P., Griffiths, J., Ray, J., Schmid, R., Collilieux, X., Garayt, B., 2012. IGS08: the IGS realization of ITRF2008. GPS Solut. 16, 483–494, doi: 10.1007/s10291011-0248-2. Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., Ozener, H., Kadirov, F., Guliev, I., Stepanyan, R., Nadariya, M., Hahubia, G., Mahmoud, S.M., Saker, K., ArRajehi, A., Paradissis, D., Al-Aydrus, A., Prilepin, M., Guseva, T., Evren, E., Dmitrotsa, A., Filikov, S.V., Gomez, F., Al-Ghazzi, R., Karam, G., 2006. GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. J. Geophys. Res. (Solid Earth) 111, 5411, doi: 10.1029/2005JB004051. provided by the SAO/NASA Astrophysics Data System. Scargle, J.D., 1982. Studies in astronomical time series analysis. ii – statistical aspects of spectral analysis of unevenly spaced data. Astrophys. J. 263, 835–853, doi: 10.1086/160554. Schoellhamer, D.H., 2001. Singular spectrum analysis for time series with missing data. Geophys. Res. Lett. 28, 3187–3190, doi: 10.1029/2000GL012698. Steigenberger, P., Boehm, J., Tesmer, V., 2009. Comparison of GMF/GPT with VMF1/ECMWF and implications for atmospheric loading. J. Geodesy 83, 943–951, doi: 10.1007/s00190-009-0311-8.

Q. Chen et al. / Journal of Geodynamics 72 (2013) 25–35 Tesmer, V., Steigenberger, P., Rothacher, M., Boehm, J., Meisel, B., 2009. Annual deformation signals from homogeneously reprocessed VLBI and GPS height time series. J. Geodesy 83, 973–988, doi: 10. 1007/s00190-009-0316-3. Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Phys. D: Nonlinear Phenom. 35, 395–424, doi: 10.1016/0167-2789(89)90077-8. Vautard, R., Yiou, P., Ghil, M., 1992. Singular-spectrum analysis: A toolkit for short, noisy chaotic signals. Phys. D: Nonlinear Phenom. 58, 95–126, doi: 10.1016/0167-2789(92)90103-T.

35

Williams, S.D.P., Bock, Y., Fang, P., Jamason, P., Nikolaidis, R.M., Prawirodirdjo, L., Miller, M., Johnson, D.J., 2004. Error analysis of continuous GPS position time series. J. Geophys. Res. 109, B03412, doi: 10.1029/2003JB002741. Zhang, J., Bock, Y., Johnson, H., Fang, P., Williams, S., Genrich, J., Wdowinski, S., Behr, J., 1997. Southern California Permanent GPS Geodetic Array: Error analysis of daily position estimates and site velocities. J. Geophys. Res. 102, 18035–18055, doi: 10.1029/97JB01380.