BDS Satellite Clock Prediction Considering Periodic Variations

Abstract

The periodic noise exists in BeiDou navigation satellite system (BDS) clock offsets. As a commonly used satellite clock prediction model, the spectral analysis model (SAM) typically detects and identifies the periodic terms by the Fast Fourier transform (FFT) according to long-term clock offset series. The FFT makes an aggregate assessment in frequency domain but cannot characterize the periodic noise in a time domain. Due to space environment changes, temperature variations, and various disturbances, the periodic noise is time-varying, and the spectral peaks vary over time, which will affect the prediction accuracy of the SAM. In this paper, we investigate the periodic noise and its variations present in BDS clock offsets, and improve the clock prediction model by considering the periodic variations. The periodic noise and its variations over time are analyzed and quantified by short time Fourier transform (STFT). The results show that both the amplitude and frequency of the main periodic term in BDS clock offsets vary with time. To minimize the impact of periodic variations on clock prediction, a time frequency analysis model (TFAM) based on STFT is constructed, in which the periodic term can be quantified and compensated accurately. The experiment results show that both the fitting and prediction accuracy of TFAM are better than SAM. Compared with SAM, the average improvement of the prediction accuracy using TFAM of the 6 h, 12 h, 18 h and 24 h is in the range of 6.4% to 10% for the GNSS Research Center of Wuhan University (WHU) clock offsets, and 11.1% to 14.4% for the Geo Forschungs Zentrum (GFZ) clock offsets. For the satellites C06, C14, and C32 with marked periodic variations, the prediction accuracy is improved by 26.7%, 16.2%, and 16.3% for WHU clock offsets, and 29.8%, 16.0%, 21.0%, and 9.0% of C06, C14, C28, and C32 for GFZ clock offsets.

1. Introduction

The real-time satellite orbit and clock products are indispensable for real-time precise point positioning (RT-PPP) [1]. The ultra-rapid orbit products with sufficient accuracy from MGEX Analysis Centers (ACs) can be used as real-time orbit correction directly [2]. However, the clock offsets are closely related to the satellite atomic clocks in space which can be influenced by many factors, such as external environmental effects that are not corrected by conventional models [3]. It is difficult to grasp the behavior of the satellite clocks accurately. Therefore, analyzing the characteristics of the clock offsets and improving the accuracy of real-time clock offsets is of great importance to RT-PPP [4,5].

For BDS satellites, the domestic atomic clocks are equipped to provide the on-board time reference. BDS-2 satellites are equipped with the rubidium (Rb) atomic clocks, while BeiDou-3 satellites are equipped with the new generation Rb atomic clocks and high-quality hydrogen (H) atomic clocks [6,7]. The BDS real-time clock products can be obtained through the ultra-rapid clock products or the real-time clock corrections. The ultra-rapid clock products are comprised of observation and prediction parts, and the prediction part is available for real-time applications. Limited by the accuracy of clock prediction model, the accuracy of the predicted clock offsets is lower than that of real-time clock corrections calculated based on observation data from global or regional network stations [8]. The real-time clock corrections enable the RT-PPP via internet. Due to a network outage, a disruption in data streaming is inevitable [9], which will interrupt the availability and stability of RT-PPP [10]. In this case, the predicted clock offsets are usually used to maintain the continuity of real-time applications [11]. Therefore, an accurate clock prediction model is of great importance. Currently, some available models, such as a quadratic polynomial model (QPM), spectral analysis model (SAM), grey mode (GM), and autoregressive integrated moving average model (ARIMA) [12,13,14,15], are given for clock prediction. Among them, the QPM is the most basic clock prediction model, modeled with the parameters of the clock offset, clock velocity, and clock drift, and intuitively reflects the physical characteristics of satellite clocks [16,17]. However, the QPM ignores the periodic noise, which presents in clock offsets [18,19]. Some studies attribute the periodic noise to the factors such as orbit determination errors, temperature variations, and perturbation errors [20,21]. The SAM, composed of quadratic polynomial and periodic terms, improves the prediction accuracy of the QPM by compensating the periodic noise, and is widely used to obtain ultra-rapid clock products [22,23,24]. To further improve the prediction accuracy of SAM, different numbers of periodic term and different length of fitting period used in SAM are evaluated. The evaluated results show that the optimal number and length vary with the clock type and its in-orbit performance [25]. The remaining nonlinear system errors in SAM residuals are modeled and compensated by the Back Propagation (BP) neural network, and the accuracy of 24 h predicted clock offsets is around 4 ns [26]. For the SAM, the compensation effect of periodic noise directly affects the fitting and prediction accuracy. Therefore, the accurate periodic terms are the prerequisite for the SAM. The FFT is usually used to extract periodic terms from the clock offsets. For good frequency resolution, long-term clock offset series are applied for FFT, such as 60-day, 100-day, and even one year [25,26,27]. However, due to space environment changes, temperature variations, and various disturbances, the periodic noise is time-varying [28]. The time-varying performance of GPS satellites based on the clock products of the International GNSS Service (IGS) is characterized, and the results show that the amplitudes of periodic terms vary with the eclipse seasons, and there is a relationship between the periodic noise and orbit dynamics [29]. The Galileo satellite clock offsets have a similar performance. A detailed time-varying analysis of the first two significant periodic terms is carried out for Galileo satellites based on the European Space Agency (ESA) clock products, and confirms that amplitude variations of the first periodic term affected by eclipse seasons, and the amplitude sinusoidal variations of second periodic term related to the J2 relativistic effect [30]. Studies show that the BDS periodic terms have some distinction between the GPS and Galileo, but similar variations of the periodic variations have been found in the BDS satellite clock offsets [31]. Considering that the FFT cannot characterize the periodic variations of the clock offsets in time domain [32], the periodic terms obtained by FFT may have certain deviations. In this case, the accuracy of the SAM will be affected, and it is necessary to quantify the periodic terms more accurately.

In this paper, we investigate the periodic variations of BDS satellites clock offsets based on the clock products from WHU and GFZ. To obtain more periodic variation details, the STFT is used, and the results are compared with FFT. The time-varying amplitude and frequency of the periodic terms in BDS satellite clock offsets are detected through STFT. Awareness of periodic variations can be used to improve the clock prediction model, we proposed a time frequency analysis model (TFAM), and the fitting performance as well as prediction performance are evaluated.

The following section is arranged as follows. First, the periodic variations of the BDS satellite clock offsets are characterized and analyzed, and the TFAM considering periodic variations is proposed in Section 2. Then, the fitting and prediction performance of the proposed TFAM are evaluated and compared with that of the SAM in Section 3. The discussions are given in Section 4. Finally, the conclusions are summarized in Section 5.

2. Materials and Methods

In this section, the periodic variations of BDS satellite clock offsets are characterized and analyzed, and the TFAM considering periodic variations is proposed. In order to obtain the clean clock offsets, a preprocessing step is implemented. This paper introduces a double median absolute deviation (MAD) detection method, which can distinguish the gross errors and phase jumps effectively. After preprocessing, the amplitude spectrums based on different FFT lengths are compared, and the limitation of FFT in analyzing the time-varying periodic noise is discussed. To fully analyze the periodic variations, the QPM fitting residuals in the time domain, the spectrum analysis results in frequency domain, and the time-frequency analysis results in time and frequency domain of three constellation satellites are given. By time-frequency analysis of STFT, the periodic variations of BDS satellite clock offsets are characterized, and the relationships between the periodic variations and sun elevation angle above the orbit plane (β angle) are investigated. After that, the frequency variations of the main periodic term are analyzed in detail. Finally, the clock prediction model is modified by taking into account the periodic variations of clock offsets, and the TFAM is proposed.

2.1. Preprocessing of Clock Offsets

The preprocessing of raw clock offsets is of great importance for periodic variation analysis and clock prediction. The gross errors and phase jumps are the main anomalies of the clock offsets, which cannot objectively reflect the characteristics of the satellite clocks and degrade the performance of clock prediction [33]. Such clock anomalies should be processed before periodic variation analysis and clock offset modeling. With high calculation efficiency and anti-error performance, the MAD method is usually used to detect gross errors [34]. When the clock frequency series satisfies Equation (1), the related clock offsets are considered abnormal:

$$\left|y_i\right|>k$$

(1)

where $y_i$ represents the clock frequency series, $i$ is the epoch, and $k$ is the threshold value, which can be calculated by $k=m+n\cdot MAD$. $m$ denotes the median of the clock frequency series, $MAD=Median\left\{\left|y_i-m\right|/0.6745\right\}$, and the factor $n$ can be set to 3–5.

The traditional MAD method is a preprocessing method corresponding to clock frequency series. Both gross errors and phase jumps of clock offsets in time domain will bring outliers in clock frequency series, which means the outliers detected by traditional MAD are caused by gross errors or phase jumps. If all the outliers are processed without distinguishing them, some effective clock offsets may be destroyed simultaneously. Therefore, it is necessary to identify the outliers caused by gross errors or phase jumps. In this paper, we use double MAD detection to distinguish the outliers. The detailed flowchart of the double MAD detection is shown in Figure 1.

As shown in Figure 1, the double MAD detection performs two MAD detection. In the first MAD detection, the raw clock offsets are converted to clock frequency series, and the threshold value $k$ is calculated. The outliers of clock frequency series can be easily detected by MAD detection. Then, store the clock offset outliers identified by the first MAD detection temporarily and eliminate the frequency outliers. After that, the clock offset series without gross errors are recovered by integral algorithm. In the second MAD detection, the new clock frequency series are recalculated, while the threshold $k$ still uses the value obtained by the first MAD detection. In this case, the outliers of clock frequency series absolutely correspond to the phase jumps, and the clock offsets, which are accidentally removed during the first MAD detection, can be recovered. After double MAD detection, the outliers caused by gross errors and phase jumps are marked respectively.

It should be noted that the gross errors and phase jumps identified by double MAD detection should be processed reasonably. In this paper, the gross errors are set to zero in the BDS satellite clock periodic variation analysis, while removed directly during modeling. The constant offset brought by phase jump is corrected to guarantee sufficient clock offset series for modeling. We use the mean of $f_2$ without outliers to substitute the outliers of the clock frequency series. After that, the clock offsets before and after the phase jump are aligned.

2.2. Analysis of Periodic Variations in BDS Satellite Clock Offsets

The atomic clocks are crucial payloads of BDS satellites and can be easily influenced by environment and temperature [16,17]. Many researchers have confirmed that the periodic noise is significant in BDS satellite clock offsets, which are probably caused by some factors such as the orbit determination errors, temperature variations, and perturbation errors [20,21]. The FFT is a traditional method used to identify and extract periodic terms from the clock offsets. To improve the frequency resolution of amplitude spectrum, the length of clock offset series used for FFT is usually long, such as 60-day, 100-day, and even one year [25,26,27]. There are no uniform criteria for selecting the length of clock offset series in published literature. The FFT based on different lengths of clock offset series may give different periodic terms, especially the main periodic term. Here, the amplitude spectrums of different lengths of clock offset series are given first.

The raw BDS satellite clock offset series used here are WHU precise clock products, from day-of-year (DOY) 001 to 270 of the year 2019 (https://cddis.nasa.gov/archive/gnss/products/mgex, accessed on 14 August 2021). Considering that the BDS precise clock products are obtained by Orbit Determination and Time Synchronization (ODTS), some orbit determination errors (particularly for the radial component) are absorbed by the clock offsets [31]. Therefore, the orbit determination strategies are essential for the analysis of precise clock offsets [35].

Table 1. Processing strategies of WHU BDS precise orbit determination [35,37].

Items	Models
Satellite attitude	BDS-2 (GEO: orbit-normal; IGSO/MEO: yaw-steering and orbit-normal); BDS-3: continuous yaw-steering
Satellite antenna PCO/PCV	PCO according to igs14.atx; ignoring PCV
Reference frame	IGb14
Geopotential	12 × 12 EGM2008 model [38]
N-body	Jet Propulsion Laboratory (JPL) DE405
ocean tides	FES2004 [39]
Solid earth tides	IERS conventions 2010 [40]
pole tides	IERS conventions 2010 [40]
Relativity	IERS conventions 2010 [40]
Solar radiation	The 5-parameter ECOM1 model [35,41]
Antenna thrust	BDS-2 (GEO: none; IGSO: 185 W; MEO: 130 W [42]); BDS-3 (MEO: 130 W)
Earth Rotation Parameters	Static Earth’s Albedo model [43]

gives crucial information of the precise orbit determination strategies with which the used WHU precise clock products are estimated. Solar radiation pressure (SRP) is the dominant non-gravitational perturbation for satellites [36], which can affect the accuracy of the satellite orbit directly. For WHU orbit strategies, the 5-parameter ECOM1 model with an empirical a priori model is used as BDS SRP model [35,37].

To analyze the influence of different FFT lengths on spectral analysis, Geosynchronous Orbit (GEO) satellite C01, Inclined Geosynchronous Satellite Orbit (IGSO) satellite C10, and Medium Earth Orbit (MEO) satellite C25 are chosen randomly as an example. First, daily clock offsets are fitted by QPM after preprocessing. Second, the different length residuals of the QPM aligned day by day are utilized for FFT analysis. The amplitude spectrums of the 90-day, 180-day, and 270-day are shown in Figure 2. The horizontal axis represents frequency with the unit cycles per day (cpd), and the vertical axis represents the amplitude with the unit ns.

From Figure 2, the periodic terms in clock offsets are detected at multiple frequencies. For C01 and C10, the 1 cpd periodic term approximately equal to their orbit period of 23.9 h. Apart from four significant periodic terms, the amplitude of other periodic terms decreases with n (n is an integer) increases. The periodic terms with low amplitude are ignored in modeling and prediction. For C25, there are only four significant periodic terms in the amplitude spectrums. Among them, the 1.86 cpd periodic term is approximately equal to MEO orbit period of 12.9 h. It can be found that the BDS satellite clock offset periodic terms closely related to the orbit period. Similar periodic terms have been found in other satellites from the same constellation.

As shown in Figure 2, it can be seen that the main periodic term obtained by FFT based on different lengths of clock offset series are distinct. For C01, the periodic term of 1 cpd is the most significant in the amplitude spectrum of 90-day, while 3 cpd in the amplitude spectrum of 180-day and 270-day. For C10, the main periodic term in 90-day amplitude spectrum is 2 cpd, while 1 cpd for the 180-day and 270-day amplitude spectrums. The spectral peak in the 90-day and 270-day amplitude spectrums for C25 is approximately 3.72 cpd, while 1.87 cpd in the 180-day amplitude spectrum. Distinct amplitude spectrum results of different FFT lengths are probably led by periodic variations. Performing FFT on clock offset series gives the global spectrum analysis result for the entire time series. However, it cannot give the variations of frequency components in the time domain. The satellite atomic clocks in space are affected by space environment changes, temperature variations, and various disturbances, which make the periodic noise time-varying. Therefore, the main periodic term obtained by FFT may be inaccurate. The STFT has a unique advantage in processing non-stationary signals. To perform STFT, long-term clock offset series are broken up into chunks, and to reduce artifacts at chunk boundary, the chunks overlap each other. Each chunk is Fourier transformed, and the magnitude for each point in time and frequency domain is recorded. Therefore, the STFT can effectively detect the variations of amplitude spectrum. In this section, the STFT is used to analyze the periodic variations of BDS satellite clock offsets.

Although the constellation deployment of BDS-3 has been completed, some satellites of BDS-3 are in orbit only for a short time, and there are no sufficient precise clock products. Therefore, the discussion in this paper is limited to 15 satellites of BDS-2 and 15 satellites of BDS-3. The BDS-3 satellites are manufactured by China Academy of Space Technology (CAST) and Shanghai Engineering Center for Microsatellites (SECM). The satellites from the two manufacturers have different body sizes and attitude modes, and employ two different satellite platforms. Distinct attitude control models employed by CAST satellites and the SECM satellites affect the accuracy of orbit determination and clock offsets [19]. Therefore, it is necessary to classify the satellite clock offsets that from the CAST or the SECM. The detailed information of the BDS satellites is listed in

Table 2. The types of BDS satellites (http://www.csno-tarc.cn/, accessed on 14 August 2021).

System	Satellite Type	Clock Type	Manuf	PRN
BDS-2	GEO	Rb	CAST	C01, C02, C03, C04, C05
	IGSO	Rb	CAST	C06, C07, C08, C09, C10, C13, C16
	MEO	Rb	CAST	C11, C12, C14
BDS-3	MEO	Rb	CAST	C19, C20, C21, C22, C23, C24, C32, C33
		H	SECM	C25, C26, C27, C28, C29, C30, C34

.

To fully analyze the periodic variations, we investigate the periodic noise of the BDS satellite clock offsets in frequency and the time-frequency domain. First, we collected the raw clock products (DOY 001 to 360, 2019), and calculated the fitting residuals of the QPM after preprocessing. Then, we perform a spectrum analysis based on FFT and a time-frequency analysis based on STFT for the residuals. Limited by paper length, one satellite from each constellation is randomly selected to show the results. Figure 3 gives the residuals, amplitude spectrum, and spectrogram of the GEO satellite C01, the IGSO satellite C07, and the MEO satellite C27.

The figure of residuals shows small gaps, meaning that the raw clock offsets of C01, C07 and C27 have good continuity. The spectrograms in Figure 3 give the information of the time domain and the frequency domain simultaneously. The horizontal axis represents the time domain, and the vertical axis represents the frequency domain. The third dimension, which is represented by the intensity of the color, indicates the Power Spectral Density (PSD) of the frequency at that time.

From the results of residuals in Figure 3, it can be seen that the amplitude of the residuals varies over time irregularly, and there is no intuitive information about periodic variations. The amplitude spectrums show the holistic periodic terms of 360-day clock offset series. The spectrograms give the variations of the periodic terms through color variations. Two types of variations can be found through the spectrograms. First, the amplitude of the periodic terms varies with time, and the variations are diverse for different periodic terms and different satellites. Second, the frequencies of periodic terms vary over time. Therefore, the main periodic term is time-varying, and the detected results by FFT and STFT are different. For C01, the highest spectral peak in frame ① is 3 cpd, while 2 cpd in frame ②. For C07, the main periodic term in frame ① and frame ② of the spectrogram is 1 cpd, and 2 cpd, respectively. For C27, there are several significant peaks around 1.86 cpd in the amplitude spectrum, and the amplitude of them has little difference. The picture-in-picture of spectrogram shows that these significant periodic terms around 1.86 cpd play the main role at different times. Such frequency fluctuations of the main periodic term are easily ignored by FFT in a short time.

Considering that the accuracy of orbits and clocks during eclipse seasons are much lower than that in non-eclipse seasons, we characterize and quantify the periodic variations in WHU precise clock products during the eclipse seasons. A correlation analysis between the PSD and β angle is performed. The PSD of four significant periodic terms from Figure 3 are extracted and analyzed. The results of BDS-2 GEO satellites C01 and C02, BDS-2 IGSO satellites C06 and C08, BDS-2 MEO satellites C12 and C14, BDS-3 MEO satellites C19 and C21 made by CAST and BDS-3 MEO satellites C25 and C27 made by SECM are shown in Figure 4. The vertical axis on the left represents the PSD of periodic terms, and the vertical axis on the right represents the β angle. When the absolute value of β angle is less than a critical value, the satellite enters the eclipse seasons [44]. For GEO and IGSO satellites, the value is 8.5 deg, while for MEO satellites the value is 14.5 deg [45,46]. The shaded regions in Figure 4 represent the eclipse seasons.

Figure 4 shows the variations of four significant periodic terms. The temporal variations of PSD can be found for each satellite in eclipse seasons and non-eclipse seasons. For GEO and IGSO satellites, the PSD of four periodic terms vary seasonally. When the eclipse seasons come, the values of PSD are gradually increasing. The 1.86 cpd periodic term in the PSD of MEO satellites varies significantly, and it tracks the eclipse seasons. Other periodic terms show limited correlation with eclipse seasons. The 1 cpd periodic term of GEO and IGSO satellites and 1.86 cpd periodic term of MEO satellites, which are closely related to their orbit period, are significantly affected by the eclipse seasons. During the eclipse seasons, the actual attitude of the satellite is inconsistent with the nominal attitude. The yaw attitude will deteriorate the accuracy of non-conservative perturbation model and the correction of phase winding and PCO, thus affecting the accuracy of the satellite orbits [44]. Apart from the yaw attitude, the thermal radiation which is complicated to be accurately modeled also affects the orbit accuracy in eclipse seasons [46]. Considering that up to 97% of the radial orbit errors can be absorbed by the clock offsets [47], the PSD variations of the periodic term related to the orbit period can be affected by low orbit accuracy caused by yaw bias, thermal radiation, and other factors. The PSD variations of other periodic terms may be affected by perturbation errors and environment variations. Compared with BDS-2 MEO satellites, the PSD variations of BDS-3 MEO satellites in eclipse seasons are more outstanding. For BDS-3 MEO satellites, the PSD variations of C19 and C21 are higher than those of C25 and C27. It may be affected by the distinct orbit accuracy of CAST satellites and SECM satellites.

From Figure 4, it can be concluded that the PSD of periodic terms vary with time, and the variations of each periodic term are inconsistent. The inconsistent variations lead to the main periodic term varying over time. For GEO and IGSO satellites, the main periodic term varies irregularly, while the variations of the main periodic term generally occur around the β angle peak value for MEO satellites. Considering that the International GNSS Monitoring and Assessment System (IGMAS) and the GFZ only add one periodic term to the SAM for the BDS ultra-rapid clock products [26], we focus on the variations of the main periodic term. The frequency variations of the main periodic term for all satellites considered are shown in Figure 5.

Figure 5 gives the frequency variations of the main periodic term. The horizontal axis and the vertical axis represent the time and the PRN, respectively. The third dimension indicates the frequency of the main periodic term. The left subfigure contains the frequency variations of the main periodic term for 30 satellites. As can be seen from the figure, the result of C11 is obviously anomalous from the others. The right subfigure gives the results of all satellites except C11. From the right subfigure, it can be seen that the main periodic terms of 29 satellites vary over time, and the frequency range is between 0 and 6 cpd. Therefore, the main periodic term obtained by FFT according to long-term clock offset series cannot be used as the main periodic term simply.

As different MGEX ACs have different orbit and clock solutions, apart from the WHU, we also analyze periodic variations based on BDS clock products from GFZ. Since GFZ start to release BDS-3 satellite precise clock products on the website of Crustal Dynamics Data Information System (CDDIS) from DOY 166 in 2020, the 360-day clock products (DOY 166, 2020 to 159, 2021) provided by WHU and GFZ are used to analyze the periodic variations. The frequency variations of the main periodic term for WHU and GFZ clock offsets are shown in Figure 6.

From Figure 6, it can be seen that the frequency variations exist both in WHU and GFZ clock products, and the inconsistent variations may be caused by the different orbit and clock solutions. In addition to C11, special frequencies are also detected in C14 clock offsets. In order to analyze the special results of C11 and C14, the 360-day (DOY 166, 2020 to 159, 2021) spectrograms are given in Figure 7.

As shown in Figure 7, C11 has an obvious periodic term around 20 cpd, and the frequency of this periodic term increases with time. For C14, special periodic terms with decreasing frequency can be found in the spectrogram. A special periodic term can be identified both from WHU and GFZ clock offsets. It may be caused by anomalies of the satellite atomic clock itself.

Based on the analysis of this section, it can be found that the main periodic term of BDS satellite clock offsets is time-varying, and the variations occur frequently in all satellites. Therefore, the frequencies of the periodic terms obtained by STFT are more accurate for compensating the periodic noise than FFT. These time-varying periodic terms should be considered to improve the clock prediction model, especially for satellites with marked periodic variations.

2.3. BDS Satellite Clock Prediction Model Considering Periodic Variations

As a widely used and representative satellite clock prediction model, the SAM is expressed as:

$$\begin{array}{ll}{x}_i=& a_0+a_1\left(t_i-t_0\right)+a_2{\left(t_i-t_0\right)}^2+{\displaystyle \sum _{r=1}^l{A}_{r}sin\left(2\pi f_r\left(t_i-t_0\right)\right)}\\ & + B_{r}sin\left(2\pi f_r\left(t_i-t_0\right)\right)+{\mathsf{\Delta }}_i \left(0\le i\le n\right), \overline{P}\end{array}$$

(2)

where $x_i$ is clock offset of epoch $t_i$, $t_0$ represents the reference epoch, and $a_0$, $a_1$ and $a_2$ are clock offset (phase), clock velocity (frequency), and clock drift (frequency drift), respectively. $l$ is the number of periodic terms, and $A_r$, $B_r$, $f_r$ denote the amplitude and frequency of the periodic terms, respectively. ${\mathsf{\Delta }}_i$ is the residual of the model at epoch $t_i$. $\overline{P} = diag\left(i(t_i-t_0)\right), i=\left(1,2,\dots ,n\right)$ is the weight matrix [26].

Generally, parameters $l$ and $f_r$ in Equation (2) are determined first by spectrum analysis, and then other parameters such as $a_0$, $a_1$, $a_2$, $A_r$ and $B_r$ are estimated using LS. From Section 2, there are multiple frequency components in BDS satellite clock offset periodic noise. Inappropriate $l$ may easily cause overfitting, which adversely affects the prediction accuracy [3]. The current research generally adds one or two periodic terms in clock prediction model to compensate for the periodic noise [26,27,28]. Considering the SAM with one periodic term is commonly used in BDS ultra-rapid clock products [26], we set $l$ to 1 in this paper. The accurate frequency parameter $f_r$ is critical for the fitting and prediction accuracy of the SAM. From the previous analysis, the amplitude spectrum calculated by FFT is a global spectrum analysis result and cannot give the periodic variations in time dimension. Awareness of time-varying periodic term can be used to improve the clock prediction model, in this paper, the STFT is utilized to detect and quantify the frequency parameter $f_r$ of the main periodic term. Based on constantly updated clock offsets, the time-varying amplitude spectrum can be obtained by STFT, and the time-varying parameter $f_r$ can be determined accurately. Since the periodic term coefficient is obtained by STFT, the improved clock prediction model is referred to as the time-frequency analysis model (TFAM).

3. Results

In this section, based on the clock products of the WHU and the GFZ, the performance of the TFAM is analyzed. First, the fitting performance of the TFAM is evaluated and compared with that of the SAM. Then, the prediction results of the TFAM and the SAM are discussed. Finally, the prediction results of the satellites with marked periodic variations are investigated to show the performance of TFAM fully.

3.1. Accuracy Analysis of Clock Offset Fitting

In order to evaluate the fitting accuracy of the improved model, 30-day (DOY 200 to 229, 2020) WHU and GFZ clock offset fitting are analyzed statistically. The root mean square (RMS) error is used as the evaluation criterion of the fitting accuracy. Figure 8 gives the RMS of the QPM, the SAM, and the TFAM fitting residuals after removing the anomalies, respectively.

Figure 8 shows the mean fitting RMS of the WHU and GFZ clock offsets, respectively. It can be seen that both the SAM and the TFAM have better fitting accuracy than the QPM. For some satellites, the TFAM performs better than the SAM, and, for others, the fitting accuracy is basically the same. Adding a time-varying main periodic term to the clock prediction model indeed improves the fitting accuracy, and for the satellites with marked periodic variations, the TFAM can further improve the fitting accuracy. The BDS-3 MEO satellites have better fitting accuracy than BDS-2 MEO satellites, which implies that the BDS-3 MEO satellite clocks have better performance in modeling. This is mainly because the BDS-3 satellites are equipped with higher quality on-board atomic clocks than BDS-2. MEO satellites present better fitting performance than the GEO and IGSO satellites, which may be caused by different orbit model performance of MEO satellites and other satellites.

In order to further evaluate the fitting accuracy, the improvement rate of the TFAM relative to the SAM is analyzed. The improvement is given in Figure 9.

As shown in Figure 9, the improvement of fitting accuracy is presented by considering the time-varying periodic term. The improvement is not only in WHU clock offsets, but also in GFZ clock offsets. Nearly half of the satellites have clear improvements. Affected by inconsistent periodic variations of satellites, the improvement is different for different satellites. The maximum improvement of a single satellite for WHU clock offsets is 36.4%, while it is 50.0% for GFZ clock offsets. For WHU clock offsets, the fitting accuracy improves 10.5% and 12.2% for GFZ clock offsets. It is clear that the fitting residuals are further mitigated by considering the time-varying periodic term in the clock prediction model.

3.2. BDS Satellite Clock Prediction Analysis

In order to evaluate the prediction performance of the TFAM, we give and compare the prediction results of the TFAM and the SAM. The fitting and prediction are processed on a daily basis. To evaluate the prediction accuracy with different prediction time, 6 h, 12 h, 18 h, and 24 h prediction results are given in this section. As the data length of precise clock products publicly released by GFZ is limited, 300-day (DOY 166, 2020 to 099, 2021) clock offsets are used as the FFT length in this section. As shown in Table 2, the BDS-3 satellites are manufactured by CAST and SECM, and they have been partitioned to analyze the prediction performance, respectively. RMS after quadratic difference processes is used as the evaluation standard of the prediction accuracy. Clock prediction results of 30-day (DOY 140 to 169, 2021) are analyzed statistically. The mean RMS of 30 days for 30 satellites are listed in

Table 3. The mean prediction RMS based on the SAM and TFAM (unit: ns).

Clock Products	System	Satellite Type	6 h			12 h			18 h			24 h
			SAM	TFAM	IMP	SAM	TFAM	IMP	SAM	TFAM	IMP	SAM	TFAM	IMP
WHU	BDS-2	GEO (CASC)	1.52	1.46	3.9%	2.52	2.52	0%	3.73	3.68	1.3%	5.16	5.09	1.4%
		IGSO (CASC)	1.75	1.50	14.3%	3.52	3.00	14.7%	5.31	4.63	12.8%	6.97	6.36	8.8%
		MEO (CASC)	0.91	0.80	12.1%	1.56	1.43	8.3%	2.32	2.16	6.9%	2.77	2.57	7.2%
	BDS-3	MEO (CASC)	0.48	0.43	10.4%	0.85	0.79	7.1%	1.28	1.19	7.0%	1.77	1.66	6.2%
		MEO (SECM)	0.35	0.31	11.4%	0.66	0.58	12.1%	1.04	0.91	12.5%	1.43	1.27	11.2%
		Total	1.00	0.90	10.0%	1.82	1.66	8.8%	2.74	2.51	8.4%	3.62	3.39	6.4%
GFZ	BDS-2	GEO (CASC)	1.83	1.81	1.1%	2.74	2.74	0%	3.76	3.75	0.2%	5.06	5.05	0.1%
		IGSO (CASC)	1.79	1.38	22.9%	3.76	2.81	25.3%	5.78	4.42	23.5%	7.64	6.15	19.5%
		MEO (CASC)	0.86	0.76	11.6%	1.50	1.34	10.7%	2.23	2.04	8.5%	2.75	2.54	7.6%
	BDS-3	MEO (CASC)	0.45	0.41	8.9%	0.87	0.76	12.6%	1.33	1.18	11.2%	1.83	1.67	8.7%
		MEO (SECM)	0.41	0.36	12.2%	0.80	0.66	17.5%	1.23	1.06	13.8%	1.70	1.49	12.3%
		Total	1.07	0. 94	12.1%	1.94	1.66	14.4%	2.87	2.49	13.2%	3.8	3.38	11.1%

.

Table 3 shows the mean prediction RMS of the WHU and GFZ clock offsets using the SAM and TFAM. It can be seen that the TFAM improves the mean prediction accuracy of the SAM, which indicates that the effect of the periodic noise can be further weakened by considering the time-varying periodic term. For different satellite types, IGSO satellites have maximum improvement, and the 24 h accuracy improvement of GFZ and WHU clock offsets achieve 19.5% and 8.8%, respectively. For the GEO satellites, there is little accuracy improvement both for WHU and GFZ clock offsets. Figure 10 gives the frequency variations of the main periodic term for GFZ and WHU clock offsets of 30 days (DOY 140 to 169, 2021). As shown in Figure 10, there is little periodic variations for GEO satellites. Therefore, little accuracy improvement is achieved for GEO satellites. The periodic variations of IGSO satellites are more frequent than MEO satellites. The improvement of MEO satellites is lower than that of IGSO, which is mainly because the MEO periodic variations generally only occur around the β angle peak value. Compared with WHU, the GFZ improvement is more obvious. It is caused by the different periodic variations between WHU and GFZ clock offsets. As shown in Figure 10, the variations of the main periodic term are mostly consistent for GFZ and WHU clock offsets. The prediction accuracy improvement of the 6 h, 12 h, 18 h, and 24 h are in the range of 6.4% to 10.0% for WHU clock offsets, while 11.1% to 14.4% for GFZ clock offsets.

In addition, the prediction accuracy of different satellites is discussed. Compared with the BDS-2 MEO satellites, the prediction accuracy of BDS-3 MEO satellites is better. Because the BDS-3 satellites are equipped with the new generation high-quality atomic clocks, they have better frequency stability than the atomic clocks on BDS-2. For different constellation satellites, the MEO satellites have the best prediction accuracy, followed by the GEO satellites and the IGSO satellites. Comparing the prediction accuracy of satellites from different manufacturers, it can be found that the prediction accuracy of SECM satellites is higher than CAST satellites. For different prediction lengths, the prediction accuracy is generally reduced with increasing of the prediction length. The 6 h prediction accuracy of the TFAM for all satellites is approximately 0.9 ns.

Considering that the time-varying periodic terms of different satellites occur at different times, the overall improvement of TFAM is limited. For the TFAM, it performs better for the satellites with marked periodic variations. In order to evaluate the performance of the improved model for the satellites with marked periodic variations, the BDS-2 IGSO satellite C06, the BDS-2 MEO satellite C14, the BDS-3 SECM satellite C28, and the BDS-3 CAST satellite C32 are chosen to show improved accuracy. Considering that there are little periodic variations of BDS-2 GEO satellites from DOY 140 to 169, in 2021, four satellites from IGSO and MEO constellations are chosen as an example. Figure 11 and Figure 12 give the amplitude spectrums and spectrograms of WHU and GFZ clock offsets, respectively.

It can be seen from the amplitude spectrums and spectrograms in Figure 11 and Figure 12, the main periodic terms detected by FFT and STFT are distinct for C06, C14, and C32. For C28, the main periodic term in the amplitude spectrum and spectrogram of Figure 11 is the same, while different in Figure 12, which may be caused by different clock offset continuity of WHU and GFZ. For C06, C14, and C32, the amplitude spectrums and spectrograms of WHU and GFZ clock offsets have slight discrepancy, while the variation trend is basically the same. The spectral peak at the amplitude spectrum of the C06 is 1 cpd, while 1 cpd during DOY 140 to 145 and 2 cpd during DOY 145 to 170 in the corresponding spectrograms. For the C14, the main periodic term at the amplitude spectrum is 1.86 cpd. A periodic variation can be found in the spectrogram of the C14 satellite. From DOY 140 to 145, the main periodic term is 5.94 cpd, while the main periodic term is changed to 3.72 cpd after DOY 145. The main periodic term at the amplitude spectrum and the spectrogram of C32 are 1.86 cpd and 3.72 cpd, respectively. Therefore, for C06, C14, and C32 in Figure 11 and all of the satellites in Figure 12, the main periodic term has marked variations, and the frequency parameter of the TFAM and the SAM is different. Such significant periodic variation could affect prediction accuracy of the clock prediction model. The 30-day prediction errors (24 h prediction) of C06, C14, C28, and C32 based on the SAM and TFAM are given in Figure 13 and Figure 14.

As shown in Figure 13 and Figure 14, the prediction errors of the TFAM and the SAM increase with prediction time, and the increase trend is similar for C06, C14, and C32 using WHU and GFZ clock offsets. For the satellites with marked periodic variations, the divergence trend of prediction errors is significantly improved using the TFAM than SAM, which indicates that the prediction errors accumulated with time can be improved by considering the time-varying periodic terms. Apparently, the TFAM can enhance the clock prediction accuracy. In order to intuitively demonstrate the prediction accuracy of the SAM and the TFAM, the mean value of the 30-day prediction errors is calculated, and the mean RMS values of the 30-day are remarked in Figure 15.

As shown in Figure 15, the prediction accuracy is significantly improved using the TFAM compared to the SAM. For WHU clock offsets, the RMS of the C06, C14, and C32 are improved by 26.7%, 16.2%, and 16.3%, respectively. The GFZ clock offset prediction accuracy improvements for C06, C14, C28, and C32 are 29.8%, 16.0%, 21.0%, and 9.0%. The TFAM performs better for the satellites with marked periodic variations.

4. Discussion

This research focused on the characteristics of the BDS satellite clock offset periodic variations, and proposed TFAM to analyze the impact of the periodic variations on the performance of clock prediction model. According to previous researchers, the amplitude of the periodic terms is time-varying. From the FFT results of different clock offset lengths, it can be seen that the FFT is not suitable to detect and quantify the periodic variations. Considering that the STFT has a unique advantage in processing non-stationary signals, it is used to analyze the periodic variations of BDS satellite clocks in this paper. Comparing the amplitude spectrums obtained by FFT and spectrograms obtained by STFT, it can be seen that both the amplitude and frequency of periodic terms vary with time, and the variations are diverse for different satellites. By analyzing the PSD variations of four significant periodic terms with respect to the β angle of BDS satellite, it is found that the amplitude variations of periodic terms due to the orbit determination errors vary significantly, and it tracks the eclipse season. For other periodic terms, the PSD variations may be caused by perturbation errors and environment variations. As the PSD variations of each periodic term are inconsistent, the frequency of main periodic term varies over time, which is found both in the WHU and GFZ clock offsets. Awareness of periodic variations may be used to improve the clock prediction model, and the TFAM considering the time-varying periodic term is established. Using STFT instead of FFT in obtaining the main periodic term, the periodic variations of the clock offsets are properly considered in the clock prediction model.

Compared with the fitting and prediction performance of the TFAM and SAM, it can be concluded that adding a time-varying periodic term to the TFAM indeed improves the fitting and prediction accuracy. Since the periodic variations of different satellites occur at different times, the overall improvement of TFAM is limited. For the satellites with marked periodic variations, the TFAM can significantly improve the fitting and prediction accuracy. Although there are some differences between the periodic variations of the WHU and GFZ clock offsets, the TFAM can improve the fitting and prediction accuracy for both WHU and GFZ clock offsets.

Due to the inconsistency performance of different satellites, the fitting and prediction results are discussed, and some conclusions have been drawn. The fitting and prediction accuracy of BDS-3 MEO satellites is better than that of BDS-2 MEO satellites. This may be because the BDS-3 satellites are equipped with high-quality atomic clocks, which have more stable frequency than atomic clocks on BDS-2 satellites. MEO satellites present better fitting and prediction performance than the GEO and IGSO satellites, and this is mainly because the orbit model performs better for MEO satellites than other satellites. Comparing the prediction accuracy of satellites from different manufacturers, the prediction accuracy of SECM satellites is higher than that of CAST satellites. It may be affected by the distinct orbit accuracy of CAST satellites and SECM satellites, which is caused by different attitude control models.

It is noted that the TFAM performs well for satellites with marked periodic variations. Except for the periodic noise, the characteristics of residual noise will be further studied in future work.

5. Conclusions

Due to space environment changes, temperature variations, and various disturbances, the periodic noise of the BDS satellite clocks is time-varying. In this paper, the periodic variations of BDS satellite clock offsets are analyzed by STFT. Then, the time-varying periodic terms are used to improve the performance of the clock prediction. Some conclusions can be summarized from this paper:

• Both the amplitude and frequency of periodic terms in clock offsets vary with time for all BDS satellites discussed in this paper.
• For satellites with marked periodic variations, the fitting accuracy of TFAM is significantly improved. Compared with SAM, the mean RMS of the TFAM fitting residuals increases by 10.5% for WHU clock offsets and 12.2% for GFZ clock offsets.
• TFAM can improve the average prediction accuracy of the SAM. Compared with SAM, the improvement rate of the prediction accuracy using TFAM of the 6 h, 12 h, 18 h and 24 h is in the range of 6.4% to 10% for WHU clock offsets, while 11.1% to 14.4% for GFZ clock offsets.
• Compared with SAM, for C06, C14, and C32 with marked periodic variations, the prediction accuracy using TFAM is improved by 26.7%, 16.2%, and 16.3% for WHU clock offsets, and 29.8%, 16.0%, 21.0%, and 9.0% of C06, C14, C28, and C32 for GFZ clock offsets.