Signal Extraction from GNSS Position Time Series Using Weighted Wavelet Analysis

The daily position time series derived by Global Navigation Satellite System (GNSS) contain nonlinear signals which are suitably extracted by using wavelet analysis. Considering formal errors are also provided in daily GNSS solutions, a weighted wavelet analysis is proposed in this contribution where the weight factors are constructed via the formal errors. The proposed approach is applied to process the position time series of 27 permanent stations from the Crustal Movement Observation Network of China (CMONOC), compared to traditional wavelet analysis. The results show that the proposed approach can extract more exact signals than traditional wavelet analysis, with the average error reductions are 13.24%, 13.53% and 9.35% in north, east and up coordinate components, respectively. The results from 500 simulations indicate that the signals extracted by proposed approach are closer to true signals than the traditional wavelet analysis.


Introduction
During the past three decades, the global or regional GNSS (Global Navigation Satellite System) station networks have been constructed for various applications [1][2][3], and nearly thirty years' of position time series of a larger number of permanent stations have been derived using GNSS tracking data. Thanks to these position time-series data, various geophysical phenomena, such as plate tectonics [4], crustal deformation [5], post-glacial rebound [6] and vertical land motion at tide gauges [7], have been investigated more clearly. Due to time variable characteristics of external environments and geophysical effects, GNSS position time series show significant seasonal changes, especially for the vertical component [8]. The main causes of seasonal oscillations are attributed to gravitational excitation, soil moisture and environmental loading variations [9][10][11]. If seasonal signals are not well characterized or there is non-random noise in the position time series, the velocity of a permanent station and its uncertainty will be misestimated [12,13]. In addition, in order to analyze the noise components using existing software CATS [14], Hector [15] and est_noise [16], all signals must be removed from the GNSS position time series beforehand. Therefore, there has been an increasing interest in extracting signals from GNSS position time series, especially seasonal signals. Since a seasonal signal is irregular and its amplitude varies over time [17,18], it cannot be characterized by using a combination of annual and semi-annual harmonic functions with constant amplitudes that are solved by least squares estimation (LSE) [19][20][21]. Moreover, besides the annual and semi-annual signals, there still exist some other periodic signals with periods of 350 days and their fractions [22][23][24]. Consequently, the nonlinear variation of a GNSS position time series cannot be reliably described via a harmonic model.
Several studies have been carried out to reliably describe and deal with time-varying seasonal signals, such as singular spectrum analysis (SSA) [25], Kalman filter (KF) [26,27] and adaptive Wiener

Binary Wavelet Analysis
Assuming that ϕ(t) is a mother wavelet function, then a family of wavelet functions can be generated by translating and scaling the ϕ(t) as [52]  where a and b denote the scaling factor and translation factor, respectively. If we perform binary discretization on scaling and translation factor, namely, taking a = 2 j , b = 2 j k, where j, k are integers, then Equation (1) is rewritten as [53] ϕ j,k (t) = 2 − j/2 ϕ 2 − j t − k The Equation (2) is a family of binary wavelet functions. For a discrete position time series of base station x = x 0 x 1 · · · x N−1 T , its j-th WT can be calculated by the following equation [54] w( j, k) = where S i represents i-th sampling interval and w( j, k) denotes k-th value of j-th wavelet coefficients. The matrix form of Equation (3) is where, w j represents the vector of j-th wavelet coefficient with dimension n j = N/2 j+1 , and the capital form W j represents the n j × N wavelet transform matrix of the j-th layer, which can be expressed as Stacking the wavelet coefficient of each layer and adding the scale coefficient to last layer v J−1 , the wavelet transform of x can be expressed as where, J denotes the number of layers to be decomposed, V J−1 is a unit row vector satisfying V J−1 ⊥W j , which can be uniquely determined according to the orthogonality condition. The wavelet transform matrix W is a standard orthogonal matrix satisfying W T W = WW T = I N , where I N denotes N-dimensional identity matrix. As the MAR process shown in Figure 1 [41], the original time series can be reconstructed by WT matrix and wavelet coefficient as follows where d j denotes the detail component of j-th layer and a J−1 denotes the appropriate component of time series. From the frequency domain perspective, the time series with frequency of f can be decomposed into several components with different frequencies. For daily GNSS position time series, its sampling frequency is one day. Considering that the annual and semi-annual signals are the two main periodic oscillations, with the periods of about 182 ∈ (2 7 , 2 8 ) days and 365 ∈ (2 8 , 2 9 ) days, respectively, the number of layers to be decomposed J is chosen as 8 and the reconstructed components of seventh and eighth level can be considered to be the annual and semi-annual signal, which is also confirmed by Klos et al. [18].

Weighted Wavelet Analysis by Considering Formal Errors
For the traditional WA approach, the Equations (1) ~ (7) are generally employed to analyze the time series without considering the prior formal errors. To introduce the weight factors according to the formal errors, the position time series is first decomposed into several components with different frequency by WA approach, and the signal and noise are separated by using correlation coefficient method [55], by which the layer with local minimum correlation coefficient is the boundary layer between signal and noise. The correlation coefficient of j-th layer is computed as follows 1 , 0 x d (8) where, Considering that the temporal correlation is not provided in GNSS time series and the precision of observations is only related to noise, not the signal [50], the weight of time series is adjusted by multiplying the noise with a weight factor i p as follows ( 0,1 , 1) where i x denotes the adjusted time series. According to the law of error propagation, the formal error of i x can be derived by where i  is the formal error of the original time series at i-th epoch. According to the definition of weight, we have

Weighted Wavelet Analysis by Considering Formal Errors
For the traditional WA approach, the Equations (1)~(7) are generally employed to analyze the time series without considering the prior formal errors. To introduce the weight factors according to the formal errors, the position time series is first decomposed into several components with different frequency by WA approach, and the signal and noise are separated by using correlation coefficient method [55], by which the layer with local minimum correlation coefficient is the boundary layer between signal and noise. The correlation coefficient of j-th layer is computed as follows where, x and d j denote the vectors of original time series and its j-th layer elements, x and d j are their correspondent average values, d i,j is the i-th element of d j and i represents the epoch index of the time series.
Consequently, both the signal and noise are reconstructed according to the boundary layer from the decomposed coefficients. At i-th epoch of the time series, the signal s i and noise e i are represented as Considering that the temporal correlation is not provided in GNSS time series and the precision of observations is only related to noise, not the signal [50], the weight of time series is adjusted by multiplying the noise with a weight factor p i as follows where x i denotes the adjusted time series. According to the law of error propagation, the formal error of x i can be derived by where σ i is the formal error of the original time series at i-th epoch. According to the definition of weight, we have where m is a constant denoting the formal error of unit weight. It is obvious from Equations (11) and (12) that each component of the adjusted time series has the same formal error m. Therefore, processing the adjusted time series with WA is equivalent to processing the original time series by introducing a Remote Sens. 2020, 12, 992 5 of 16 weight factor, for this reason the proposed approach is called weight wavelet analysis (WWA). In order to keep the total energy of the adjusted time series unchanged [50,51], the sum of the weights for all epochs should be equal to the number of epochs, i.e., From Equations (12) and (13), the weight factor p i and constant m can be derived as follows Once the weight factor is determined by Equation (15), the adjusted time series can be derived by Equation (10), where initial signal and noise are derived from the original time series by WA approach. Then the adjusted time series is iteratively processed by using WA approach to extract signal and noise until the difference ∆s i of the signals between two adjacent iterations satisfies the following threshold [50] |∆s where ε is a small value used to terminate the iterative procedure, which is suggested to be chosen as 0.001 [50]. The algorithm of WWA is summarized in Figure 2.

Real GNSS Position Time Series Analysis
In order to verify the performance of our new approach, the position time series of 27 permanent stations of CMONOC were analyzed with the comparison of traditional WA and their corresponding locations are presented in Figure 3. The dataset of position time series was downloaded from the GNSS data product service platform of China Earthquake Administration (http://www.cgps.ac.cn/). The CMONOC provided the data solved by Bernese and GAMIT, respectively. Since the time series

Real GNSS Position Time Series Analysis
In order to verify the performance of our new approach, the position time series of 27 permanent stations of CMONOC were analyzed with the comparison of traditional WA and their corresponding locations are presented in Figure 3. The dataset of position time series was downloaded from the GNSS data product service platform of China Earthquake Administration (http://www.cgps.ac.cn/). The CMONOC provided the data solved by Bernese and GAMIT, respectively. Since the time series derived by Bernese is too short, we adopted the GAMIT solution; for specific processing details refer to ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf. Outliers were eliminated by interquartile range (IQR) criterion [56] and observations with the formal errors larger than 10 mm and 20 mm respectively for horizontal and vertical components were discarded [50]. Offsets were detected and corrected by sigseg software [57] and manual visual inspection. The gaps in the data were supplemented with zero values.

Real GNSS Position Time Series Analysis
In order to verify the performance of our new approach, the position time series of 27 permanent stations of CMONOC were analyzed with the comparison of traditional WA and their corresponding locations are presented in Figure 3. The dataset of position time series was downloaded from the GNSS data product service platform of China Earthquake Administration (http://www.cgps.ac.cn/). The CMONOC provided the data solved by Bernese and GAMIT, respectively. Since the time series derived by Bernese is too short, we adopted the GAMIT solution; for specific processing details refer to ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf. Outliers were eliminated by interquartile range (IQR) criterion [56] and observations with the formal errors larger than 10 mm and 20 mm respectively for horizontal and vertical components were discarded [50]. Offsets were detected and corrected by sigseg software [57] and manual visual inspection. The gaps in the data were supplemented with zero values. Here we take BJFS station as an example to specifically show the processing details of WWA and WA approach. Figures 4 and 5 present the position time series and formal errors of BJFS station during the period from year 1999 to 2018, respectively, which clearly indicate that the precision varies with the epoch and vertical component is obviously worse than horizontal component. Considering that the Coiflet wavelet is near symmetric and has shown to be excellent for the sampling approximation of smooth functions [53,59], thus we adopted the coif-5 wavelet as the mother wavelet. The coif-5 wavelet was employed to perform 8-layer decomposition on the detrended position time series; the reconstructed component of each layer are shown in Figure 6. To distinguish the frequency (period) of reconstructed components, the Lomb-Scargle algorithm [60] was adopted; the corresponding results are presented in Figure 7. It is obvious we can conclude that WA can decompose the position time series into several components with different frequencies efficiently. The correlation coefficients between original time series and reconstructed components are presented Here we take BJFS station as an example to specifically show the processing details of WWA and WA approach. Figures 4 and 5 present the position time series and formal errors of BJFS station during the period from year 1999 to 2018, respectively, which clearly indicate that the precision varies with the epoch and vertical component is obviously worse than horizontal component. Considering that the Coiflet wavelet is near symmetric and has shown to be excellent for the sampling approximation of smooth functions [53,59], thus we adopted the coif-5 wavelet as the mother wavelet. The coif-5 wavelet was employed to perform 8-layer decomposition on the detrended position time series; the reconstructed component of each layer are shown in Figure 6. To distinguish the frequency (period) of reconstructed components, the Lomb-Scargle algorithm [60] was adopted; the corresponding results are presented in Figure 7. It is obvious we can conclude that WA can decompose the position time series into several components with different frequencies efficiently. The correlation coefficients between original time series and reconstructed components are presented in Table 1, from which we can find that the boundary layers between signal and noise are d 4 , d 5 and d 5 for north, east and up components, respectively. According to Equation (14), the formal errors of unit weight for north, east and up components were 1.01 mm, 1.05 mm and 3.03 mm; the corresponding weight factors are presented in Figure 8. Then we analyzed and extracted the signals from the position time series of BJFS station with WWA and WA approaches; the results are presented in Figure 9, where red and blue lines denote the results from WA and WWA and green line denotes the difference between them. From Figure 9, we can find that both two approaches can better describe the non-linear variation of the GNSS position time Remote Sens. 2020, 12, 992 7 of 16 series and the extracted signals have slight difference, with mean absolute differences are 0.09 mm, 0.10 mm and 0.23 mm for north, east and up components, respectively. After the signals are removed from time series, the average error of residual time series can be calculated with Equation (17) for WA and Equation (18) for WWA as followsσ whereê WA andê WWA are the residuals by WA and WWA, S denotes the epoch set of observations and M is the number of dataset. Since the total weights are same for WA and WWA, the average error can be used to characterize the quality of extracted signal [50]. For BJFS station, the average errors of north, east and up components are 0.74 mm, 0.88 mm and 3.04 mm for WWA and 0.97 mm, 1.11 mm and 3.64 mm for WA, respectively. Therefore, we can conclude that WWA can extract more signals than the WA approach.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 17 in Table 1, from which we can find that the boundary layers between signal and noise are 4 5 , d d and 5 d for north, east and up components, respectively. According to Equation (14), the formal errors of unit weight for north, east and up components were 1.01 mm, 1.05 mm and 3.03 mm; the corresponding weight factors are presented in Figure 8. Then we analyzed and extracted the signals from the position time series of BJFS station with WWA and WA approaches; the results are presented in Figure 9, where red and blue lines denote the results from WA and WWA and green line denotes the difference between them. From Figure 9, we can find that both two approaches can better describe the non-linear variation of the GNSS position time series and the extracted signals have slight difference, with mean absolute differences are 0.09 mm, 0.10 mm and 0.23 mm for north, east and up components, respectively. After the signals are removed from time series, the average error of residual time series can be calculated with Equation (17) for WA and Equation (18) for WWA as follows 2 ( ) (18) where ˆW A e and ˆW WA e are the residuals by WA and WWA, S denotes the epoch set of observations and M is the number of dataset. Since the total weights are same for WA and WWA, the average error can be used to characterize the quality of extracted signal [50]. For BJFS station, the average errors of north, east and up components are 0.74 mm, 0.88 mm and 3.04 mm for WWA and 0.97 mm, 1.11 mm and 3.64 mm for WA, respectively. Therefore, we can conclude that WWA can extract more signals than the WA approach.     To further testify the performance of WWA with the comparison of traditional WA, we analyzed all 27 permanent stations series of CMONOC. Figure 10 presents the formal errors of unit weight derived by Equation (14)     Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 17 0.8213, 0.9002 and 0.8322 for north, east and up components, respectively. This indicates that the more the weight factor varies, the higher the relative improvement will be.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 17 0.8213, 0.9002 and 0.8322 for north, east and up components, respectively. This indicates that the more the weight factor varies, the higher the relative improvement will be.   As can be seen from Figure 11, the average errors of WWA are smaller than WA for all stations, which indicate that WWA performs better than WA on signals extraction from GNSS position time series. The mean improvements of WWA relative to WA are 13.24%, 13.53% and 9.35% for north, east and up components, respectively. Similar to the formal error of unit weight and the standard deviation of weight factor, the improvements of three components are also strong correlated, with the correlation coefficient of 0.8730 between north and east, 0.8051 between north and up and 0.8845 between east and up. Moreover, the relative improvements for three components are consistent with those of the standard deviation of weight factors in Figure 10 and their correlation coefficients are 0.8213, 0.9002 and 0.8322 for north, east and up components, respectively. This indicates that the more the weight factor varies, the higher the relative improvement will be.

Synthetic Time Series Analysis
To further evaluate the performance of the WWA approach, the simulation experiments were carried out in this section, since true signals are known in simulation cases. The simulated strategy was the same as the literature [50]. Here, the extracted signals of BJFS station with WA are treated as true signal S 0 and the noise e 0 is generated based on the formal errors using the Matlab function mvnrnd. Therefore, the simulated time series X 0 generated by adding S 0 and e 0 as follows Considering that the true signal is known, the root mean square error (RMSE) and mean absolute error (MAE) of the extracted signal can be calculated with where, S 0 (t) andŜ(t) denote the true signal and the extracted signal at the t epoch, M is the length of dataset. The smaller the RMSE and MAE, the closer the extracted signal will be to the true signal. In order to obtain statistically reliable results, we repeat the simulation experiments 500 times. For each experiment, the signal was extracted by WA and WWA, and then the RMSE and MAE were computed based on Equations (21) and (22); the specific results are presented in Figures 12 and 13. Obviously, all RMSEs and MAEs of WWA were smaller than those of WA, which indicates that WWA performs better than WA in extracting signal from GNSS time series. The mean RMSE and MAE of signal extracted by WA and WWA are shown in Table 2, as well as the improvements of WWA relative to WA, which are calculated as follows From

Discussion and Conclusions
Various seasonal signal and noise are contained in GNSS position time series. An effective separation of signal and noise not only helps to obtain a more accurate velocity estimate of permanent station, but also lays the foundation for noise analysis. Initially, the seasonal variation of GNSS position time series was modelled by harmonic model [19][20][21], which assumed that the amplitude of signal was constant. However, it can be seen from Figure 9 that the amplitude of nonlinear seasonal signals varies with epoch to epoch, which is also demonstrated by others [17,18]. Moreover, it has been verified that the WA can decompose the position time series into different scales ( Figure 6) and the reconstructed component of each layer has different frequencies ( Figure 7). Obviously, the reconstructed component of seventh (d7) and eighth level (d8) represent semi-annual and annual signals, respectively, which consists with Klos et al. [18]. It seems that the annual signal is more correlated with the original time series (Table 1) than semi-annual signal, which indicates that the annual signal is more significant.
This finding confirms the assumption of Kaczmarek and Kontny [17], in which the wavelet power spectrum analysis has been used to measure the amplitudes of periodic signals.
Since the formal errors of position time series of GNSS permanent station are also provided in daily solution, WWA is proposed in this paper for extracting the time-varying seasonal signal from a GNSS position time series. The proposed approach adjusts the weight of position time series by constructing weight factors via formal errors. When formal errors are homogeneous, the WWA will be reduced to WA. The real position time series of 27 permanent stations from CMONOC over the period from 1999 to 2018 were analyzed using WWA and WA. Both methods can extract the non-linear seasonal variation from noisy position time series. After removing the signals from original time series, the average errors of residual time series by WWA were all smaller than those by WA for all stations and the mean improvements are 13.24%, 13.53% and 9.35% for north, east and up components, respectively. Furthermore, the improvement seemed to be correlated with the standard deviation of weight factors. In order to further illustrate the superiority of WWA relative to WA, the 500 simulations were carried out. The results show that the signal extracted by WWA was closer to true signal compared with WA. For the north component, the mean RMSE and MAE of signal extracted by WWA relative to WA reduced from 0.5294 mm to 0.4770 mm and 0.4184 mm to 0.3354 mm, with improvements of 9.90% and 19.84%, respectively. For the east component, the corresponding mean RMSE and MAE reduced from 0.3708 mm to 0.3601 mm and 0.2543 mm to 0.2476 mm, with improvements of 2.89% and 2.63%. For the up component, the corresponding mean RMSE and MAE reduced from 0.9990 mm to 0.9684 mm and 0.7180 mm to 0.6964 mm, with improvements of 3.06% and 3.01%. Both the results of real and synthetic time series demonstrated that WWA performs better than WA. Therefore, the formal errors should be considered when extracting the signals from GNSS position time series using wavelet analysis.
Author Contributions: Y.S. proposed the key idea, supervised the research and designed the experiments; K.J. carried out the data analysis and wrote the paper; Y.S. and F.W. revised the manuscript. All authors have read and agreed to the published version of the manuscript.