Spectral analysis of multi-year GNSS code multipath time-series

In the presented study multi-year time series of changes in the L1 pseudo-range multipath are analysed. Data from 8 stations of the EUREF Permanent Network (EPN) were used in the study. Periodic components present in the signal and their stability over time were analysed. Also, the type of background noise was determined, based on the spectral index. In some cases, the presence of weak components with a 1/2 and 1/3 of the Chandler period has also been found. Time-frequency analysis shows that periodic signals are not stationary in most of the examined cases, and particular signal components occur only temporarily. The analysed signals were characterised by pink noise in the lower frequency range and by white noise for higher frequencies, which is also characteristic for time series of coordinates obtained from GNSS measurements.


Introduction
GNSS multipath is studied mainly as a source of errors affecting code and phase observations in positioning. The impact of multipath on GNSS data is widely analysed in both time and frequency domain [1][2] [3]. The magnitude of the pseudo-range measurement error caused by multipath can be determined based on linear combinations of code and phase observations [4]. The impact of multipath on phase observations is estimated based on changes in signal power reaching the GNSS receiver (signal to noise ratio) [5][6] [7].
The analysis of the multipath phenomenon applies, among others in GNSS reflectometry, where it is often used to study the change in the environment of the GNSS antenna. Periodic oscillations in the multipath are the basis for determining the change in the position of the surface reflecting the satellite signal about the phase centre of the antenna. Multipath oscillations analyses allow assessing changes in snow cover thickness [8], ocean surface [9] [10], and soil moisture [11].
At EPN/IGS reference stations, indices of the multipath are routinely monitored [12] [13], with basic statistical parameters characterising the variability of the multipath for daily observations on individual carrier frequencies being mainly calculated. For this purpose, software such as G-Nut/Anubis [14] or teqc [15] is used.
Apart from short-term changes, the multipath is also characterised by seasonal variability, resulting from annual changes in vegetation, soil moisture, and snow or ice layer thickness. Therefore, when analysing the long time series of multipath, one should expect components with a period close to one year. Since the value of the error caused by multipath is also affected by the position of the antenna relative to reflecting obstacles, even small changes in the position of the antenna can also cause changes in indices describing the multipath.
In this work, multi-year time-series of indices characterising GNSS multipath for 8 European permanent stations belonging to the EPN/IGS network ( Fig. 1) were analysed in terms of the occurrence of periodic components. The analysed time-series of L1 code multipath were obtained using freely available teqc software.

Methodology
The multipath in code pseudoranges at L1 frequency is described by the value of the m p1 index, calculated as a linear combination of phase and code observations for a single epoch in accordance to the equation [15]: fected by the position of the antenna relative to reflecting obstacles, even small changes in ition of the antenna can also cause changes in indices describing the multipath. In this work, multi-year time-series of indices characterising GNSS multipath for 8 an permanent stations belonging to the EPN/IGS network (Fig. 1) were analysed in terms occurrence of periodic components. The analysed time-series of L1 code multipath were d using freely available teqc software.

dology
The multipath in code pseudoranges at L1 frequency is described by the value of the mp1 calculated as a linear combination of phase and code observations for a single epoch in ance to the equation [15]: ultiparty error, pseudorange, -L1 phase, -L2 phase, 1/L2) 2 .
The multi-year time series were created on the base of the root mean square of mp1 index aily observations. The calculated value of mp1 depends on many factors. In addition to the enon of interference of signals arriving in different ways, the mp1 is affected by changes by the replacement of the antenna, receiver or firmware upgrade. A change in the uration of the measuring equipment can cause a jump in the value of mp1 as well as a (1) where: m p1 -multiparty error, p1 -L1 pseudorange, also affected by the position of the antenna relative to reflecting obstacles, even small changes in the position of the antenna can also cause changes in indices describing the multipath.
In this work, multi-year time-series of indices characterising GNSS multipath for 8 European permanent stations belonging to the EPN/IGS network ( Fig. 1) were analysed in terms of the occurrence of periodic components. The analysed time-series of L1 code multipath were obtained using freely available teqc software.

Methodology
The multipath in code pseudoranges at L1 frequency is described by the value of the mp1 index, calculated as a linear combination of phase and code observations for a single epoch in accordance to the equation [15]: where: mp1 -multiparty error, p1 -L1 pseudorange, The multi-year time series were created on the base of the root mean square of mp1 index from daily observations. The calculated value of mp1 depends on many factors. In addition to the phenomenon of interference of signals arriving in different ways, the mp1 is affected by changes caused by the replacement of the antenna, receiver or firmware upgrade. A change in the configuration of the measuring equipment can cause a jump in the value of mp1 as well as a -L1 phase, also affected by the position of the antenna relative to reflecting obstacles, even small changes in the position of the antenna can also cause changes in indices describing the multipath.
In this work, multi-year time-series of indices characterising GNSS multipath for 8 European permanent stations belonging to the EPN/IGS network ( Fig. 1) were analysed in terms of the occurrence of periodic components. The analysed time-series of L1 code multipath were obtained using freely available teqc software.

Methodology
The multipath in code pseudoranges at L1 frequency is described by the value of the mp1 index, calculated as a linear combination of phase and code observations for a single epoch in accordance to the equation [15]: where: mp1 -multiparty error, p1 -L1 pseudorange, The multi-year time series were created on the base of the root mean square of mp1 index from daily observations. The calculated value of mp1 depends on many factors. In addition to the phenomenon of interference of signals arriving in different ways, the mp1 is affected by changes caused by the replacement of the antenna, receiver or firmware upgrade. A change in the configuration of the measuring equipment can cause a jump in the value of m as well as a -L2 phase, also affected by the position of the antenna relative to reflecting obstacles, even small changes in the position of the antenna can also cause changes in indices describing the multipath.
In this work, multi-year time-series of indices characterising GNSS multipath for 8 European permanent stations belonging to the EPN/IGS network ( Fig. 1) were analysed in terms of the occurrence of periodic components. The analysed time-series of L1 code multipath were obtained using freely available teqc software.

Methodology
The multipath in code pseudoranges at L1 frequency is described by the value of the mp1 index, calculated as a linear combination of phase and code observations for a single epoch in accordance to the equation [15]: where: mp1 -multiparty error, p1 -L1 pseudorange, The multi-year time series were created on the base of the root mean square of mp1 index from daily observations. The calculated value of mp1 depends on many factors. In addition to the phenomenon of interference of signals arriving in different ways, the mp1 is affected by changes caused by the replacement of the antenna, receiver or firmware upgrade. A change in the The multi-year time series were created on the base of the root mean square of m p1 index from daily observations. The calculated value of m p1 depends on many factors. In addition to the phenomenon of interference of signals arriving in different ways, the m p1 is affected by changes caused by the replacement of the antenna, receiver or firmware upgrade. A change in the configuration of the measuring equipment can cause a jump in the value of m p1 as well as a change in the noise level or amplitude of the examined signal. There may also be data gaps in the time series due to the temporary inactivity of the reference station [12] [13].
For research purposes, 8 permanent stations were selected from Europe ( Fig. 1) with the longest continuous observation possible and the smallest number of changes in equipment in the analysed period. Data in RINEX format comes from the EUREF Permanent Network data centres, and daily RMS m p1 values have been calculated using teqc software.
The main goal of the research was to verify the occurrence of the characteristic periodic components in the m p1 times series. Before the spectral analysis, data preprocessing had to be carried out. For this purpose, the jumps in value, the linear trend and outliers (> 5 ise level or amplitude of the examined signal. There may also be data gaps in the the temporary inactivity of the reference station [12][13]. poses, 8 permanent stations were selected from Europe ( Fig. 1) with the longest rvation possible and the smallest number of changes in equipment in the Data in RINEX format comes from the EUREF Permanent Network data RMS mp1 values have been calculated using teqc software. goal of the research was to verify the occurrence of the characteristic periodic he mp1 times series. Before the spectral analysis, data preprocessing had to be this purpose, the jumps in value, the linear trend and outliers (> 5) were e obtained mp1 time series, which caused additional data gaps in the time series. uity in the observational data, the method of Lomb-Scargle periodogram was d is effective for spectrum analysis in the case of non-evenly data sampling [16] ce of predetermined periodic components in the signal has not been assumed in on the spectrum obtained, the presence of components with the tropical year ays) and the Chandler wobble period (433 days) and their three subsequent verified. S observations, the time series of coordinates are most often analysed [18].
wer-law P∝ f k , where P is power spectral density, f is the frequency and k it can be shown that in this type of observation, background noise is a ink noise and white noise [19] [20]. ork, for the background noise determination, the spectral index k was estimated a straight line to the power spectral density expressed on a logarithmic scale. of stationarity in the frequency domain was performed with wavelet transform, tool for the characterisation of frequency variations in the signal [21] that ation on both the time and frequency of the signal. Individual harmonics in the eries were analysed using the continuous wavelet time-frequency spectrum with t wavelet. ) were removed from the obtained m p1 time series, which caused additional data gaps in the time series. Due to discontinuity in the observational data, the method of Lomb-Scargle periodogram was used. This method is effective for spectrum analysis in the case of non-evenly data sampling [16] [17]. The existence of predetermined periodic components in the signal has not been assumed in the paper. Based on the spectrum obtained, the presence of components with the tropical year period (356.25 days) and the Chandler wobble period (433 days) and their three subsequent harmonics were verified.
For GNSS observations, the time series of coordinates are most often analysed [18]. Based on the power-law P change in the noise level or amplitude of the examined signal. There may also be data gaps in the time series due to the temporary inactivity of the reference station [12] [13]. For research purposes, 8 permanent stations were selected from Europe ( Fig. 1) with the longest continuous observation possible and the smallest number of changes in equipment in the analysed period. Data in RINEX format comes from the EUREF Permanent Network data centres, and daily RMS mp1 values have been calculated using teqc software.
The main goal of the research was to verify the occurrence of the characteristic periodic components in the mp1 times series. Before the spectral analysis, data preprocessing had to be carried out. For this purpose, the jumps in value, the linear trend and outliers (> 5) were removed from the obtained mp1 time series, which caused additional data gaps in the time series. Due to discontinuity in the observational data, the method of Lomb-Scargle periodogram was used. This method is effective for spectrum analysis in the case of non-evenly data sampling [16] [17]. The existence of predetermined periodic components in the signal has not been assumed in the paper. Based on the spectrum obtained, the presence of components with the tropical year period (356.25 days) and the Chandler wobble period (433 days) and their three subsequent harmonics were verified.
For GNSS observations, the time series of coordinates are most often analysed [18]. Based on the power-law P∝ f k , where P is power spectral density, f is the frequency and k spectral index, it can be shown that in this type of observation, background noise is a combination of pink noise and white noise [19] [20].
In this work, for the background noise determination, the spectral index k was estimated based on fitting a straight line to the power spectral density expressed on a logarithmic scale. The verification of stationarity in the frequency domain was performed with wavelet transform, a widely used tool for the characterisation of frequency variations in the signal [21] that preserves information on both the time and frequency of the signal. Individual harmonics in the examined time series were analysed using the continuous wavelet time-frequency spectrum with a complex Morlet wavelet. f k , where P is power spectral density, f is the frequency and k spectral index, it can be shown that in this type of observation, background noise is a combination of pink noise and white noise [19] [20].
In this work, for the background noise determination, the spectral index k was estimated based on fitting a straight line to the power spectral density expressed on a logarithmic scale. The verification of stationarity in the frequency domain was performed with wavelet transform, a widely used tool for the characterisation of frequency variations in the signal [21] that preserves information on both the time and frequency of the signal. Individual harmonics in the examined time series were analysed using the continuous wavelet time-frequency spectrum with a complex Morlet wavelet.

Results
The studied RMS m p1 time series covered a period from 13.2 to 20.6 years. At the preprocessing stage, some outliers were removed from the time series obtained directly as a result of the calculations (eq. 1). Depending on the permanent station, between 3% and 26% of the original data was deleted. Based on spectral analysis, the occurrence of annual components was determined at the significance level of 95%. In the case of the BACA, BOLG, DELF and WSRT stations, in addition to the annual component and its higher harmonics, weak signals of 1/2 or 1/3 of the Chandler period are visible (Fig. 3). These weak signals of small amplitude, which appeared in the m p1 spectrum require further research. The summary of data regarding the time series for each permanent station tested are presented in Table 1. All the analysed stations show power-law behaviour in the entire spectrum. It was estimated a corner frequency equals 6 cycles per year (cpy) between the low-frequency part where pink noise appears and the higher frequency with dominating uncorrelated white noise. Spectral indices k1 and k2 were calculated, by dividing the frequency range into two parts, relative to the corner frequency: k1 -below 6 cpy and k2 -above 6 cpy (middle column in Figures 2a and 2b). In most cases, the background noise in the low-frequency part of the spectrum is close to pink noise, and as the frequency increases, the noise character becomes white noise (-1.5 < k1 < -0.5 < k2). For the AQUI and JOEN stations, pink noise is dominant in the entire frequency range. Table 1 summarises the values of spectral indices k1 and k2 for the permanent stations analysed.
Time-frequency analysis clearly shows the non-stationarity of long-term multipath changes. In most of the analysed cases, components with an annual period and its higher harmonics appear only in a part of the examined time interval (Fig. 2a and 2b -

Conclusions
The time series of the m p1 parameter describing the signal's multipath in the observed L1 code was analysed. For the analysis 8 GNSS permanent stations from Europe were chosen. The frequency spectrum obtained with the Lomb-Scargle periodogram indicates the presence of 1 cpy frequency components and some higher harmonics. In the signal spectrum from several reference stations, the second and third harmonic of the Chandler frequency are also visible. The annual component has been expected because it reflects seasonal changes in the vicinity of the GNSS antenna. The occurrence of the Chandler wobble period in the m p1 time series requires detailed research. Background noise in the obtained spectrum can be described as pink noise for lower frequencies and white noise for frequencies above 6 cpy. A similar distribution of noise in the spectrum is also characteristic for GNSS coordinate time series. Time-frequency analysis allows confirming non-stationarity and high variability of signal parameters over time. It is caused mainly by the equipment replacement in the reference stations and modification of the operating parameters of the receivers. Orange (dash-dot) vertical lines correspond to the yearly period and its first three harmonics. Green (dash) vertical lines correspond to the Chandler wobble period and its first three harmonics