Identiﬁcation of BDS Satellite Clock Periodic Signals Based on Lomb-Scargle Power Spectrum and Continuous Wavelet Transform

: Onboard satellite clocks are the basis of Global Navigation Satellite Systems (GNSS) operation, and their revolution periods are at the level of 2 per day (about 12 h) in the case of the Medium Earth Orbit (MEO) satellites. In this work, the authors analysed the entire BeiDou Navigation Satellite System (BDS) space segment (BDS-2 and BDS-3) in terms of the occurrence of periodic, repetitive signals in the clock products, and checked if they coincide with the orbital periods or their multiples. The Lomb-Scargle (L-S) power spectrum was used as a tool to determine the periods present in the BDS clock products, allowing for analyses based on incomplete input data; in this case, the incomplete data were the phase data with jumps and outliers removed. In addition, continuous wavelet transform (CWT) was used to produce a time − frequency representation showing the more complex behaviour of the satellite clock products. As shown in the case of geostationary and geosynchronous inclined orbit satellites, the main period was 23.935 h, while for the Medium Earth Orbit it was 12.887 h, with the BDS satellite orbital period being 12 h 53 m (12.883 h). Some effects connected with reference clock swapping are also visible in the power spectrum. The conducted analyses showed that the BDS-2 satellite clocks have much higher noise than the BDS-3 satellite clocks, meaning that the number of designated periods is greater, but their reliability is signiﬁcantly lower. BDS-3 satellites have only been in operation for a very short time, thus this is the ﬁrst analysis to include this type of data. Moreover, such a wide and complex analysis has not been carried out to date.

Since October 2020, WUM has been providing final clock correction data for C38-C46 BDS-3 satellites, and thus precise clock correction for all active BeiDou satellites has become available from IGS MGEX data centres ( Figure 1). Our analysis covers the last 88 days of 2020 (5 October 2020-31 December 2020) and was performed for BDS-2 and BDS-3 satellite clocks. We focused on periodic effects connected with the orbital period of different types of satellites (MEO, GEO, and IGSO). The resulting analysis presents the periods in the range of 0.8-36 h. The lower limit was set to 0.8 h in order to take into account the higher harmonics of the periodic changes of the clock correction. A summary of the analysed satellite, its orbit, and clock type information are provided in Appendix A.

Methods
The processing algorithm is shown in the diagram in Figure 2. The subject of analysis was raw daily clock correction data (phase data) obtained from WUM with 30 s sampling intervals and was decimated in the pre-processing phase to a 300 s sampling interval. As an example of data analysis, the processing stages for satellite no. C32 are presented in Figure 3. The analysis was performed on differentiated data (frequency data).

Methods
The processing algorithm is shown in the diagram in Figure 2. The subject of analysis was raw daily clock correction data (phase data) obtained from WUM with 30 s sampling intervals and was decimated in the pre-processing phase to a 300 s sampling interval. As an example of data analysis, the processing stages for satellite no. C32 are presented in Figure 3. The analysis was performed on differentiated data (frequency data). Clock metrology uses a form of phase or fractional frequency measurements relative to a reference clock. The phase reading δt is instantaneous (timestamp), while the fractional frequency reading is the difference between the two phases' values divided by their interval, according to [44]:  Clock metrology uses a form of phase or fractional frequency measurements relative to a reference clock. The phase reading δt is instantaneous (timestamp), while the fractional frequency reading is the difference between the two phases' values divided by their interval, according to [44]: The differentiation of the periodic signal changes its amplitude a and shifts phase ϕ by π/2, but has no impact on signal component frequency f. The amplitudes are scaled by 2πf, according to the following equations: (3) The first step of our analysis was to differentiate clock correction data δt (Figure 3a) to produce a time series of clock correction rate dδt/dt (Figure 3b). The outliers present in the time series were removed two-fold. First, values of dδt/dt larger than 10 −10 s/s (arbitrarily chosen) were treated as outliers and were removed from the data set ( Figure 3b). The threshold value was chosen based on the analysis of the whole data set and is adopted as 10 −10 s/s for proper detection of the jumps in the signal (phase data). The jumps that existed in the data series were mainly caused by the day-boundary effect [18]. The maximum number of removed data points at this stage was 0.8%. Such a small amount of removed data has no impact on further data processing. After removal of the primary outliers, jumps in clock correction rate data remained identifiable for several satellites, as is visible in Figure 3C. These discontinuities were removed ( Figure 3d) from the data. Next, the third order trend (polynomial) and secondary outlier detection and removal were performed. At this stage, outliers were detected with the aid of the MAD (median absolute deviation) method commonly used in clock stability analysis [44] (Figure 3e). The MAD coefficient is calculated according to: In Equation (4), x i is the frequency data and med(x) is a median of x, and then outliers in the time series are identified if: In our calculations, the typical value for constant k = 5 was adopted [44], and after this stage, the maximum removed data point was 2.2%. Lastly, the Lomb-Scargle [45,46] power spectrum was calculated ( Figure 3f).
Herein, we used the Lomb-Scargle method to obtain an approximation of the power spectrum of the analysed signal when gaps in the data occurred due to the removal of outliers or missing observations from the data source. In the resulting L-S spectrum, only period values for which the probability of false identification lower than P (spectrum peak | noise) = 0.00005 were selected as significant. A relatively low value of P was chosen in order to filter out spurious periods present in the analysed clock data. In Figure 3f, meaningful periods are those for which the peaks are above the orange line, i.e., 6.4 h and four peaks indicating periods over 30 h. The relevant periods for each analysed BeiDou satellite are presented as a blue mark in Appendix B.
In our calculations, the typical value for constant k = 5 was adopted [44], and after this stage, the maximum removed data point was 2.2%. Lastly, the Lomb-Scargle [45,46] power spectrum was calculated (Figure 3f).  As a last step of our analysis the time−frequency spectrum has been calculated using Continuous Wavelet Transform. A CWT is used for time series stationarity analysis. A comprehensive description of the application of the wavelet transform can be found in [47].
The Morlet wavelet was adopted as the mother wavelet. This kind of wavelet is commonly used in the time−frequency analysis and provides results with equal variance in time and frequency domains [48]. CWT allows for the examined signal to be presented in the form of a scalogram-a time−frequency spectrum, which permit for the identification of both frequencies present in the signal and their location in time. CWT is defined as a convolution of the analysed signal s(t) with scaled and shifted wavelet function ψ(t) [49]: where, a represents the scale b denotes the translation ψ represents the mother wavelet, and * denotes the complex conjugate The normalization factor 1/ √ a is used to preserve the direct comparability of the transform for all of the scales. The Morlet wavelet used here to calculate the time−frequency spectrum is expressed as follows: where f denotes the center frequency and f b represents the bandwidth. The continuous wavelet transform method does not optimize the time and frequency resolution simultaneously. In the CWT analysis, a trade-off must be found between these quantities [41]. In the present study, the center frequency was taken as f = 1.5 Hz and the bandwidth parameter was f b = 2.0 for all pf the clock data time series analyzed. The parameter values were selected empirically to ensure a satisfactory balance between the time and frequency resolution for the analyzed data.

Results and Discussion
In the analysed data (clock correction rates) after removal of the outliers and trend, numerous discontinuities are clearly visible (e.g., Figure 3e). These jumps are less visible for the older BDS-2 than for the BDS-3 satellites because of the relatively higher noise level in the clock correction time series. As the offsets are the same for all satellites (see Appendix C), their origin is also common and must be connected with some maintenance works or with a processing strategy adopted for clock correction determination. Moreover, as is indicated in the clock data files obtained, during the analysed time span, the reference clock for the WUM clock products was switched 44 times ( Figure 4). The clock switch occurs at the end/beginning of the day and eight different clocks-hydrogen masers-worked as a reference clock. The switching of the reference clock leads to jumps in analysed clock data, visible for all BDS-3 satellites (e.g., Figure 3e).  For verification, if the presence of this effect had an impact on the resulting power spectrum, we performed an analysis of the simulated data. Simulated data were Gaussian noise of amplitude 2 × 10 −13 s (similar to the BDS-3 satellite clock rate signal) with jumps of 0.2 ns added in the same epochs as the jumps occurring in the real data ( Figure 5). These For verification, if the presence of this effect had an impact on the resulting power spectrum, we performed an analysis of the simulated data. Simulated data were Gaussian noise of amplitude 2 × 10 −13 s (similar to the BDS-3 satellite clock rate signal) with jumps of 0.2 ns added in the same epochs as the jumps occurring in the real data ( Figure 5). These discontinuities are visible in the BDS-3 satellite clock data only. Simulated data processed the same way as real clock data (described in Section 2.2) allows us to confirm that the signals of periods of 17 h and 30-40 h are caused by frequent reference clock swapping. An L-S periodogram of the simulated data is presented in Figure 5. It shows a similar pattern for higher periods as periodograms for all of the analysed satellite clocks (Appendix C). The same effect is visible on the time−frequency spectrograms (see Appendix D). For verification, if the presence of this effect had an impact on the resulting power spectrum, we performed an analysis of the simulated data. Simulated data were Gaussian noise of amplitude 2 × 10 −13 s (similar to the BDS-3 satellite clock rate signal) with jumps of 0.2 ns added in the same epochs as the jumps occurring in the real data ( Figure 5). These discontinuities are visible in the BDS-3 satellite clock data only. Simulated data processed the same way as real clock data (described in Section 2.2) allows us to confirm that the signals of periods of 17 h and 30-40 h are caused by frequent reference clock swapping. An L-S periodogram of the simulated data is presented in Figure 5. It shows a similar pattern for higher periods as periodograms for all of the analysed satellite clocks (Appendix C). The same effect is visible on the time−frequency spectrograms (see Appendix D). An example spectrogram for satellite C32 is shown in Figure 6. Visible is a component with a period of about 6.4 h, which appears from about 7200 h and continues until the end of the period analyzed. An example spectrogram for satellite C32 is shown in Figure 6. Visible is a component with a period of about 6.4 h, which appears from about 7200 h and continues until the end of the period analyzed.    Figures 3e and 5). Spectrograms for all of the analyzed satellites are shown in Appendix D.
The identified periods in the clock correction rate depend on the satellite orbit type. For GEO and IGSO satellites, the main period of 23.935 h corresponding to the sidereal day is present. Moreover, higher harmonics are visible-up to the ninth in the case of the C01 satellite-which shows that the periodic changes are not sinusoidal. The orbital period for the BeiDou MEO satellites is about 12.887 h. All BDS-3 MEO satellites have dominant second harmonics for this period (Appendices B-D).
The IGSO BDS-3 satellites (C38, C39, and C40) equipped with H-maser clock exhibit many low power spectrum peaks around 24 h and 12 h. In the case of these satellites, orbital periods close to the sidereal day and its harmonics are of course natural. In the BDS-3 segment, only these three satellites are IGSO and they all have an H-maser clock. The presence of signals with periods between 12 and 24 h, which are not visible in the other BDS-3 satellites data, may depend on these factors. Several BDS-2 satellites' clock corrections oscillate with very low periods, e.g., C10, 1.02 h; C11, 1.21 h; and C14, around 1.8 h, not connected with their orbital periods. These effects require further investigation.

Conclusions
GNSS satellite clocks are a key factor in GNSS time transfer and precise positioning. Periodic effects present in satellite clock data may transfer to the GNSS products, giving false signals in areas such as geodynamical analysis and time transfer. This is of great importance, especially in Precise Point Positioning (PPP) applications with zero-difference data processing. In this paper, the authors prove the existence of identifiable effects connected with satellite orbital periods and common for the same satellite type (MEO, GEO, and IGSO). Other effects, present in all of the analysed signals, and giving a similar spectral structure, may be connected with the processing strategy applied for clock correction determination or maintenance works in the time scale-reference clock switch. Although the jumps caused by switching the reference clock are relatively small (ca. 2 ns), it is worth noting that such effects are not clearly visible and are not documented in final products in any way. Moreover, the usage of the L-S method allows for a determination of the periods not connected with satellites orbital movement nor with the time scale maintenance in several BDS-2 satellite clocks. This kind of phenomenon needs further, extended investigation. The time−frequency analysis shows that the satellite clock behaviour may vary in time, and determined periods have variable amplitudes over the analysed time span. The analysed data may also include periodic components related to the sampling interval, but this phenomenon was not analysed in this study. The presented analysis was carried out only for the BeiDou system. In future research, it is necessary to study this type of phenomenon also for other GNSS systems, possibly with the use of special software packages developed for analysis unequally spaced data (e.g., [50]) to verify potential limitations of the methodology used.  Appendix B Figure A1. BeiDou BDS-2 (shaded grey) and BDS-3 satellites clock correction periodicity. The green lines correspond to BeiDou MEO satellites orbital period and its 2nd, 3rd and 4th harmonics. The orange lines correspond to GEO and IGSO satellites orbital periods and their harmonics up to 9. Spurious periods of 17 h and 30-40 h visible in the clock correction are caused by frequent jumps in the analysed signal (same as for simulated SIM data). The letters on the right side of diagram provide information on the particular satellite orbit and its clock: G-geostationary orbit, I-inclined geosynchronous orbit, M-Medium Earth orbit, R-rubidium clock, H-hydrogen maser. Figure A1. BeiDou BDS-2 (shaded grey) and BDS-3 satellites clock correction periodicity. The green lines correspond to BeiDou MEO satellites orbital period and its 2nd, 3rd and 4th harmonics. The orange lines correspond to GEO and IGSO satellites orbital periods and their harmonics up to 9. Spurious periods of 17 h and 30-40 h visible in the clock correction are caused by frequent jumps in the analysed signal (same as for simulated SIM data). The letters on the right side of diagram