Extracting Seasonal Signals in GNSS Coordinate Time Series via Weighted Nuclear Norm Minimization

: Global Navigation Satellite System (GNSS) coordinate time series contains obvious seasonal signals, which mainly manifest as a superposition of annual and semi-annual oscillations. Accurate extraction of seasonal signals is of great importance for understanding various geophysical phenomena. In this paper, a Weighted Nuclear Norm Minimization (WNNM) is proposed to extract the seasonal signals from the GNSS coordinate time series. WNNM assigns different weights to different singular values that enable us to estimate an approximate low rank matrix from its noisy matrix. To address this issue, the low rank characteristics of the Hankel matrix induced by GNSS coordinate time series was investigated ﬁrst, and then the WNNM is applied to extract the seasonal signals in the GNSS coordinate time series. Meanwhile, the residuals have been analyzed, obtaining the estimation of the uncertainty of velocity. To demonstrate the effectiveness of the proposed algorithm, a number of tests have been carried out on both simulated and real GNSS dataset. Experimental results indicate that the proposed scheme offers preferable performances compared with many state-of-the-art methods.


Introduction
Global Navigation Satellite System (GNSS) coordinate time series often appear in many geophysical and geodetic applications. In addition to long-term trends and noise, GNSS coordinate time series usually shows significant seasonal variations that mainly manifest as annual and semi-annual oscillations [1][2][3][4]. A common research task in GNSS time series post-processing is the extraction of the seasonal signal from a GNSS observable time series. Accurate extraction is helpful not only for improving the reliability of GNSS time series, but is also vital for the studies of various geophysical phenomena, such as the seasonal variations of crustal movement and the velocity of plate movement [5][6][7][8][9][10].
To retrieve the seasonal signals of GNSS coordinate time series, several methods have been investigated by numerous scholars during the past several decades [3,4,[11][12][13][14][15][16][17][18]. A detailed review can be found in [19]. The common method for devising a linear model is least-squares (LS) fitting [14], however, LS lacks random changes in the estimation of signals, and thus only seasonal signals with constant amplitude can be obtained [3], which is not consistent with the true seasonal variations. It has been suggested in several other studies to extract the periodic signals by data-driven methods, such as wavelet decomposition (WD) [13] and singular spectral analysis (SSA) [11,16,20]. Recently, Klos et al. compare and contrast these algorithms under different noise levels [17]. The results show that promising results have been achieved with these methods, however, most of these methods have several limitations. For instance, for WD, the wavelet basis function and the number of wavelet decomposition layers must be specified in advance, which will greatly affect the effect of wavelet denoising [13]. Although SSA can be used to effectively extract seasonal signals from univariate coordinate time series in the time domain without any prior information, it does not perform well when an inappropriate lag window is selected [11].
Weighted Nuclear Norm Minimization (WNNM) is a representative approach for solving low rank matrix approximation (LRMA) problems, which has been widely used in many fields, such as image denoising [21,22], image completion and repainting [23], seismic data denoising [24] and MIMO channel estimation [25]. WNNM outline several commonly used low-rank matrix methods, such as NNM, which can actually be viewed as a special counterpart of WNNM in which the weight vector of WNNM is set the same. Recently, Zha et al. compared and analyzed WNNM and NNM from the perspective of group sparse representation (GSR), mathematically proving that WNNM is more flexible than NNM for many practical problems [26]. In summary, rational weight will enhance the WNNM model for achieving better estimation. As is well known, seasonal signals in GNSS time series have been represented by sums of sinusoids with annual and semi-annual cycles [3]. More recently, Cai et al. proved that the Hankel matrix induced by the harmonic model has low-rank characteristics [27]. To the best of our knowledge, however, there has not been any work in which WNNM has been applied to GNSS time series. In this context, to differentiate these seasonal components in GNSS coordinate time series, the WNNM is applied in this study to investigate the possibility of extracting seasonal signals from GNSS coordinate time series. To circumvent this problem, the low-rank characteristics of the Hankel matrix induced by GNSS coordinate time series was investigated. Furthermore, the WNNM method was utilized to extract signals from the GNSS coordinate time series under different noise models and noise levels, and the results are compared with several state-of-the-art methods. Finally, the noise (i.e., the so-called residuals) was analyzed, and the corresponding parameters were calculated.
The rest of this paper is organized as follows. In Section 2, the GNSS time series model used in this study is briefly introduced, and the detailed WNNM methodology analyzed. The efficiency of WNNM in extracting the seasonal signals from simulated GNSS time series is tested in Section 3. Section 4 applies the WNNM to real GNSS data and compares the proposed algorithm with state-of-the-art algorithms. Section 5 concludes the paper.

Model
In general, the GNSS coordinate time series of a single component of a single station are routinely modeled as follows [2,4,11,[28][29][30]: where t(t = 0, 1, · · · , n − 1) represents the epoch, y(t), x(t), ε(t) represent the observation, signal and noise at the t epoch, written as y t , x t , ε t for simplicity, respectively. α, β are the initial displacement and velocity, respectively. a j (t), b j (t)(j = 1, 2) represent the instantaneous amplitudes, which are assumed to consist of a mean value (that is, a j , b j , j = 1, 2) plus a random component. f 1 , f 2 corresponding to the annual and semi-annual frequencies, respectively. The optimal stochastic model of the observation noise for most stations has been considered as a combination of white noise (WN) and flicker noise (FN) [1,4,[28][29][30][31][32].
Cai et al. proved that the Hankel matrix induced by the harmonic model has low-rank characteristics [27]. Inspired by this, the Hankel matrix induced by GNSS coordinate time series having a low-rank property was further investigated in this paper. A brief conclusion is given here via Theorem 1. Before stating the theorem, a linear operator H is presented to map the vector z ∈ n to the Hankel matrix as follows [27]: where Hz ∈ n 1 ×n 2 , the index of the vector and matrix is from 0 and [·] ij represents the (i, j)-th entry of matrix. In order to explore the low-rank characteristics of the Hankel matrix induced by the GNSS time series, let e t be the new noise term, including the observation noise ε t and the random variations of the seasonal signal. In this case, x t is still used to represent the signals, including the trend term and seasonal signals with amplitude a j , b j . Thus, the Hankel matrices induced by y t , x t and e t are easy to obtain via Equation (2), denoted as Y, X,E, then Equation (1) can be expressed as where , and n 1 = n − n 2 + 1 < n 2 . For the Hankel matrix X above, it is proved that it is of low rank using the following theorem. Theorem 1. Suppose GNSS coordinate time series is a regular sample in the time domain, i.e., t ∈ Z + , and the seasonal signal in x t is a mixture of r sinusoids with amplitude a j ,b j , the corresponding Hankel matrix is X ∈ n 1 ×n 2 , then the following inequality rank(X) ≤ 2r + 2, holds, where rank(X) denotes the rank of X, and r is the order of seasonal signals.

Proof. See Appendix A for details.
Theorem 1 indicates that the Hankel matrix X is low rank with a certain condition. Thanks to this useful conclusion, we can study the problem of seasonal signal extraction in GNSS coordinate time series from a new perspective. In this paper, a WNNM based method, which exploiting the low rank characteristic of the Hankel matrix, is proposed for extracting the seasonal signal from the GNSS coordinate time series.

Weighted Nuclear Norm Minimization
The low-rank characteristic of the Hankel matrix induced by the GNSS time series is proved in the preceding section. In this subsection, our intention is to use the WNNM to investigate the GNSS coordinate time series to obtain the seasonal signals therein, leading to the following WNNM-based model: where, X w, * = ∑ i=1 |w i σ i (X)| represents weighted nuclear norm, σ i (X) denotes the i-th singular value of X, w i ≥ 0 is a weight assigned to σ i (X), and w = [w 1 , · · · , w n ] T is the weight vector. Although the weight vector benefits the corresponding WNNM-based model for achieving better estimation. The difficulty of solving a WNNM model, however, is that it is non-convex in general cases, and thus the traditional approaches used for NNM problems are no longer applicable. To address this issue, a weighted nuclear norm proximal (WNNP) operator strategy has been investigated in detail in [23]. If Y ∈ n 1 ×n 2 , n 1 ≥ n 2 , let Y = UΣV T be the singular value decomposition (SVD) of , following Theorem 1 in [23], the global optimum solution of Equation (4) can be solved as follows: where and d 1 , d 2 ..., d n is the solution of the following convex optimization problem: More specifically, Equation (6) is a quadratic programming problem with linear constraints, which imply that Equation (6) can be solved by off-the-shelf convex optimization solvers. As introduced by Corollary 1 in [23], the global solution of Equation (6) iŝ provided the weights are sorted in a non-descending order. Following Equations (5)-(7), it can be seen that Equation (4) can be equivalently converted to a quadratic programming problem with a linear constraint that is a convex optimization problem, and thus the solution of Equation (4) can be obtained using off-the-shelf convex optimization tool [33]. Obviously, the weight vector, as the key parameter of WNNM-based models, plays a crucial role in the success of WNNM. Following the suggested setting of [23], w i = C √ n σ i (X)+ε is set herein, where C is a compromising constant that is set by experience. n is the size of the matrix, σ i (X) represent the i-th singular value of X and ε is a small positive constant to avoid dividing by zero.
Based on the relationship between the singular value of Y and the constant C, the singular value of estimated matrix X * can be obtained. If ε is small enough to make ε < min{ √ C, C σ 1 (Y) } hold (ε is set to 10 −6 in this paper), by using the weight vector w i = C √ n σ i (X)+ε , the closed-form solution of Equation (4) is where and The procedure of WNNM is summarized in Algorithm 1.

Results of Simulation
To validate the effectiveness of the proposed algorithm, it is tested on the simulated GNSS coordinate time series here. The WNNM-based algorithm is compared with several state-of-the-art methods, including moving ordinary least squares (MOLS) [17], WD [13] and SSA [11]. Following the previous studies, it can be found that these methods have achieved outstanding performance. Here, the similar strategies used in these studies are adopted. For MOLS, the time series is divided into small segments, each of which is 2-years long and overlap each other for 1-year. For WD, the seventh and eighth Meyer's wavelets are also used, as in [17], to model the seasonal signals. For SSA, Klos et al. analyzed the relationship between the lag window and the noise level, found that longer windows (3-years) produce better performance in the presence of higher noise levels, smaller windows (2-years) behave for low noise levels [17]. In this study, this method was also used to select the appropriate lag window.
To make the simulated time series more reasonable, GNSS time series is generated via Equation (1), where the seasonal signals are amplitude modulated. In all of the experiments, the parameters are set to α = 10 mm, β = 5 mm/yr. The mean amplitude of the annual and semi-annual signals is set as 3 and 2, respectively. The variations of the stochastic part of a j (t) and b j (t) are set to have standard deviations of 1.0 and 0.5, respectively. To quantitatively evaluate the performance, the reliability of WNNM is evaluated in terms of trend uncertainty, noise spectral index, noise amplitude and misfit. The trend uncertainty refers to the variance of the velocity estimation, which can be calculated by the general formulas (29)-(31) in [34]. k is the spectral index, representing the slope under the logarithm of the power spectral density (PSD) of the noise. Misfit represents the standard deviation between the signals estimated with a certain method and the signal that was simulated. To facilitate comparison of the performance of competing methods, first the effectiveness of the proposed WNNM-based method under the "FN" was verified. However, the optimal noise model of GNSS time series is well described as a combination of WN and FN, and thus the proposed WNNM-based method has also been verified under "WN+FN". To distinguish the different noises, σ w and σ f are used to represent the amplitude of WN and FN, respectively. In the following section, two cases will be discussed. The first case is that the noise model is pure FN, and in the other case a WN is added back into the synthetic time series. In each case, two different signal models, with or without a long-term trend, are considered separately.

Case 1: Pure FN
In this subsection, the performance of the proposed algorithm on the simulated time series, which was corrupted by the pure FN, is investigated. Following the method presented in the preceding section, the WNNM-based and competing methods were implemented, and the results are shown in Figure 1, which depicts the signals extracted by different methods for low noise levels, that is, The original GNSS time series is also plotted. It can be seen that the WNNM-based method, as well as the other competing methods, can achieve good performance in extracting the signals from the original GNSS coordinate time series, either in the case of without trend or with a linear trend. To compare the performance of the WNNM-based method with the other methods from the point of view of frequency, the residuals were also analyzed using maximum likelihood estimation (MLE) [7,28,29,35,36]. The PSD of residuals obtained by applying various methods are shown in Figure 2, and the PSD of original GNSS time series and synthetic noise are also plotted for comparison. As shown in Figure 2, the WNNM-based approach, as well as the other methods, can extract the signals from GNSS time series. The lower power-spectrum curve of WD may be due to the fact that the method absorbs some power, which is also shown in several previous studies.
To validate the results from a quantitative point of view, the corresponding parameters of the residuals are calculated and the results are summarized in Table 1, where the last row, "Actual", shows the actual spectral index, noise amplitude and trend uncertainty of the simulated noise used in this subsection. As can be seen in Table 1, for a lower noise level, the misfit of the WNNM-based method is significantly lower than that of other methods, especially when the long-term trend is removed. Furthermore, the spectral indices of other methods are generally closer to 0, which may be due to the absorption of some colored noise in the process of extracting seasonal signals. This is also verified in Figure 2. As can be seen from Figure 2, the residual of WNNM is closest to the synthetic noise, while the residuals of several other methods are lower than that of synthetic noise, especially at the annual and semi-annual epoch.  Although the WNNM-based method achieves a better estimation compared to the other methods under low noise level (i.e., σ f = 1 mm/yr − k 4 ), further work must be done, as the noise in a real GNSS series is usually not very small. To evaluate the performance of the WNNM-based method under different noise levels, its performance when a higher noise level is employed was further verified (i.e., σ f = 10 mm/yr − k 4 ) [17]. As with the higher noise levels, the signals extracted using various methods and the corresponding parameters of the residuals are given in Figure 3 and Table 2, respectively. The corresponding PSD were presented in Figure 4. It can be observed that the PSD of residuals show an obvious downward trend at the annual and semi-annual epoch. Among them, WNNM is the most gentle, followed by MOLS and SSA, which indicate that WNNM can achieve a better performance in extracting signals.
From the experimental results, it can be seen that, with increasing noise level, the performance of all methods used in this paper decreased. More specifically, when the long-term trend is removed, the performance of the WNNM-based approach is slightly better compared to that of MOLS, WD and SSA. Meanwhile, WNNM can achieve a good performance as the other competing methods, even when a linear trend is added back into the GNSS time series. In general, under pure FN model, the WNNM method achieves excellent estimation at different noise levels. However, it is well known that noise types in the GNSS coordinate time series are very complex, and it is not sufficient to consider the pure flicker noise. Based on this fact, the performance of the proposed WNNM-based method under "WN+FN" is further verified in the next section.

Case 2: WN+FN
As described above, for most of the stations the optimal noise model of the GNSS time series is well described as a combination of WN and FN [31]. To demonstrate the availability of WNNM in the application of GNSS time series, the performance of WNNM under the combination of WN and FN is discussed in this subsection. More specifically, Gaussian white noise with zero mean and variance 1 is added to these GNSS time series.
As shown in Figure 5 and Table 3, the WNNM-based approach can also achieve better performance under mixed noise. Compared with the other methods, the accuracy of WNNM-based method still has some advantages. Further more, for lower noise level, the spectral indices of other methods are generally closer to 0, which is consistent with the previous section. With the increasing of FN, the WNNM-based method also has a good ability to extract seasonal signals from the GNSS time series. The results are shown in Figure 6 and Table 4. Since the PSD of residuals is similar to those shown in Figures 2 and 4, we omit it here for brevity.
Experimental results indicate that the seasonal signals extracted by the WNNM method are more similar to the synthetic seasonal signals with smaller misfit under different noise models. Compared with the other methods, the power spectrum curve of residuals obtained by using WNNM method are also more similar to the power spectrum curve of simulated noise. In summary, it can be concluded that the WNNM-based method proposed in this paper can achieve good performance under different noise levels, whether the noise model is pure FN or a mixture of WN and FN. Compared with the case in which a linear trend is added, the WNNM-based method can achieve better estimation by removing the trend firstly.

Application to Real Data
In this section, the WNNM is applied to analyze the real GNSS coordinate time series. Nine GNSS stations in China are selected to assess the WNNM-based method as well as the competing methods. The locations and relevant information about each station are shown in Figure 7 and Table 5, respectively. Most of the information in Table 5, such as epochs of offset, are taken from information provided by Scrips Orbit and Permanent Array Center (SOPAC (http://sopac-ftp.ucsd.edu/pub/ timeseries/)). The last column in Table 5 shows the optimal noise model (ONM) of each site in the Up component, which follows from [37] (please refer to [37] for details).
The "Detrend" GNSS coordinate time series used in this study are available from SOPAC. The longest one is SHAO station (6049 d) and the shortest is LHAS station (2790 d). There are two main reasons for choosing these data. First, SOPAC uses "st_ f ilter" software to combine the solutions obtained by the Jet Propulsion Laboratory (JPL) (using GIPSY software) and by the SOPAC (using GAMIT/GLOBK software) analysis center, to eliminate errors due to software and processing strategies and finally obtain a unified, self-consistent solution (a detailed combination strategy is presented in Section 2.2 [19]). Second, the data do not contain long-term trends, and the ONM of most of these sites is "WN+FN" [37]. Although the ONM of some stations (such as URUM, GUAO and WUHN) are not "WN+FN", but they were still retained to verify the robustness of our method under different noise models. Although the Detrend time series was adopted, before the experiments were started, a standard pre-processing procedure, such as missing data interpolation [38] and offsets extraction [39], were also performed and the gross error, if existing, were removed using the 3σ criterion [29,40]. Owing to the fact that this is not the focus of this study, the data pre-processing strategies will not be discussed in detail, interested readers are referred to the relevant literature.  As can be seen from Table 5, most of these sites have approximately 30% missing data, and to minimize the impact of missing data, the analysis was limited to 2000 days for each site. Figure 8 shows the time series before and after pre-processing with BJFS and TNML sites as examples. It can be seen that before the pre-treatment there are missing and outliers in the coordinate time series, and after pre-processing these data are well suppressed. Following the methods introduced in the preceding section, the pre-processed data from the nine sites were analyzed. Here, the graphical results of three directions (North, East and Up) of BJFS station are shown. The time series extracted by applying MOLS, WD, SSA and WNNM are shown in Figure 9a-c, respectively, and the observed GNSS time series is also illustrated for comparison. It can be seen from Figure 9a-c that the WNNM method can accurately extract seasonal signals in all directions of this site. As shown in Figure 9b, there is still a slight linear trend in the East component, for which the LS should be first used to remove this trend. The results of removing the linear trend term are shown in Figure 9d. Obviously, the time-varying amplitude periodic signals of BJFS station obtained by the four methods are close to each other.
Results of trend uncertainty, spectral index and noise amplitude estimated from the GNSS time series of BJFS are provided in Table 6. Different from the simulation, the true signal in the real data is unknown, it is impossible to calculate the corresponding misfit as in the simulation experiment. Here, the validity is verified by calculating the correlation coefficients between the signals extracted by various methods. Results of BJFS station are shown in Table 7. As can be seen from Table 7, the highest correlation coefficient with WNNM was MOLS, followed by SSA. The correlation coefficients between the signals obtained by each method are relatively large. For the three components, the correlation coefficients are about 0.99. The results indicate that the four methods mentioned in this paper are all effective in estimating the actual periodic signals with time-varying amplitude.    Moreover, processing schemes of the other sites are the same as that of BJFS, which will be omitted due to space limitations. Similarly, the residuals obtained by various methods in all stations were also analyzed. The results are presented in Table 8. According to the results, one can see that the noise amplitude and noise spectral indices of different sites are different, and even different directions of the same site have certain differences, which have been demonstrated and verified by previous studies. The corresponding spectral indices for the North and East components were estimated to vary between −1.0636 and −1.6248, −1.0859 and −1.5159, respectively. For the vertical components, the spectral indices of the power-law process of the sites URUM, WUHN and GUAO were estimated to range from −1.3459 to −1.6447, with a median of −1.4872. The spectral indices of the remaining sites were estimated to be a mean of −1.1022, with the minimal and maximal being, respectively, −1.  Table 9. As shown in Tables 8 and 9, for most stations, results estimated by applying WNNM are not significantly different from the results obtained by the other methods. Correlation coefficients of the seasonal signals for the Up component of the other eight stations are also presented in Table 10. As shown in Table 10, the mean correlation coefficients between the signals extracted by WNNM and those by MOLS, WD and SSA were 0.9223,0.9641 and 0.9711, respectively. In addition, it is not hard to find that the correlation coefficients between MOLS and the other methods are small at the URUM station. This is primarily due to the large changes during 2018, which leads to the poor estimation of MOLS. Over all, the correlation between WNNM and the other methods was largely that verified the effectiveness of WNNM method in extracting the seasonal signals. The results of horizontal components are similar to these. We omit it here.

Discussion and Conclusions
GNSS coordinate time series contains obvious seasonal variations. In recent years, the previously published papers show that the harmonic model with time-varying amplitude can better model the seasonal oscillations. As a result of the time-varying amplitude, the traditional least square method is no longer tried. The data-driven methods, such as WD and SSA, are another main branch along this line of research. These methods can extract the seasonal signal accurately without knowing the prior information of the signal, but they also have some shortcomings. The WNNM method utilizes the low rank characteristics of Hankel matrix formed by GNSS time series. Different from the SSA that simply drops the smaller singular value, which will lose some useful information, WNNM shrinks the singular values according to their amplitude, which makes it more flexible in dealing with practical problems.
In this paper, seasonal signals were extracted from the GNSS coordinate time series by using the WNNM method. This is based on the fact that the Hankel matrix constructed by GNSS time series is low rank. To work around this issue, the Hankle matrix induced by the GNSS coordinate time series is investigated, the results demonstrate that it has low-rank characteristic under certain conditions. Then, a WNNM-based method was therefore proposed to extract the signals from the GNSS observations. The capability of WNNM approach to extract seasonal signals from GNSS time series was explored by extensive experiments that have been carried out on simulation time series. The results indicate that the seasonal signals extracted by WNNM are closer to the simulated signals, and smaller misfits were obtained under different noise models and noise levels. Whereas, seasonal signals extracted by WD method may be contaminated by colored noise more or less. In the application of real GNSS coordinate time series, the estimates of trend uncertainty, spectral indices and noise amplitude after applying WNNM are not significantly different from the results obtained by the other methods. This implies that the WNNM as well as other methods, can accurately extract seasonal signals in GNSS observations. As shown in Tables 7 and 10, it can be found that the highest correlation with WNNM is SSA, with a mean of 0.9734 for the Up component. For most of the nine sites, the correlation between WD and SSA is also larger, which may be due to the fact that they absorb some colored noise when extracting seasonal signals. The specific reasons are for further study. To conclude, the WNNM-based method is a reliable method for extracting the seasonal signals in GNSS coordinate time series. The WNNM method is not only beneficial to better study the noise in the geodetic time series, but also provides strong support for further improving the estimation accuracy of site velocities and their uncertainty estimates.
Sanchez of University of Extremadura, Spain, for providing the SSA code, which were used for comparison. We appreciate the SOPAC (http://sopac-ftp.ucsd.edu/pub/timeseries/) for providing the GNSS time series.

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

Appendix A. Proof of Theorem 1
Proof. The signals in GNSS coordinate time series is: x(t) = α + βt + r ∑ j=1 a j cos(2π f j t) + b j sin(2π f j t)), where α, β are the initial displacement and velocity, respectively. r is the order of seasonal signals, a j ,b j are the mean value of amplitudes. The seasonal signals can be expressed as a complex form by using Euler's formula, that is: d j e i2πλ j t , (t = 0, 1, · · · , n − 1), where, q = 2r, A j e iφ j , j = 1, 2, · · · , r, 1 2 A j−r e −iφ j−r , j = r + 1, r + 2, · · · , q, λ j = f j , j = 1, 2, · · · , r, − f j−r , j = r + 1, r + 2, · · · , q. Denote x(t) = x t (t = 0, 1, · · · , n − 1), the corresponding Hankel matrix X can be constructed by using Equation (2) where n 1 = n − n 2 + 1. The following matrices can be obtained by splitting the matrix X: where T, S are defined as trend matrix, seasonal matrix, respectively. For matrix T, owing to the GNSS coordinates time series are sampled uniformly in days, that is, the value of t is a positive integer, so the T can be expressed as follows: For matrix S, define y k = e i2πλ k for k = 1, ..., q, Cai et al. proved that S is low rank by using the vandermonde decomposition of Hankel matrix [27], that is S = E L DE T R , where,  If all x t is distinct, and q ≤ min(n 1 , n 2 ), the E L ,E R are both full rank matrixs, rank(S) = q, when all d j = 0. In summary, rank(X) = rank(T + S) ≤ rank(T) + rank(S) = q + 2 min{n 1 , n 2 } hold, thus, the hankel matrix induced by GNSS coordinate time series is low rank.