Time–frequency analysis of the Galileo satellite clocks: looking for the J2 relativistic effect and other periodic variations

Abstract

When observed from the ground, the frequency of the atomic clocks flying on the satellites of a Global Navigation Satellite System is referred to as apparent frequency, because it is observed through the on-board signal generation chain, the propagation path, the relativistic effects, the measurement system, and the clock estimation algorithm. As a consequence, the apparent clock frequency is affected by periodic variations of different origins such as, for example, the periodic component of the J2 relativistic effect, due to the oblateness of the earth, and the clock estimation errors induced by the orbital estimation errors. We present a detailed characterization of the periodic variations affecting the apparent frequency of the Galileo clocks, obtained by applying time–frequency analysis and other signal processing techniques on space clock data provided by the European Space Agency. In particular, we analyze one year of data from three Galileo Passive Hydrogen Masers, flying on two different orbital planes. Time–frequency analysis reveals how the spectral components of the apparent frequency change with time. For example, it confirms that the amplitude of the periodic signal due to the orbital estimation errors depends on the angle between the sun and the orbital plane. Moreover, it allows to find a more precise estimate of the amplitude of the J2 effect, in agreement with the prediction of the general theory of relativity, and it shows that such amplitude suddenly decreases when the corresponding relativistic correction is applied to the data, thus validating the analytical formula used for the correction.

Introduction

Space qualified atomic clocks are at the core of any Global Navigation Satellite System (GNSS), flying on board of each satellite of the system. When received and observed on ground, the frequency of such space clocks is referred to as apparent frequency, since it is the clock frequency observed through the signal generation chain of the satellite and the measurement system, and it is affected by environmental and relativistic effects.

The apparent frequency is affected by periodic variations of different origins, among which the periodic component of relativistic effects such as time dilation and gravitational redshift, whose amplitude is proportional to the eccentricity of the satellite orbit. A review of GNSS relativistic effects can be found in Ashby (2003) and in Nelson (2011). The periodic relativistic effects must be corrected at the user level, when the receiver software computes the navigation solution. However, the implemented relativistic corrections usually neglect the contribution due to the oblateness of the earth, which is hereafter referred to as the J2 relativistic effect. This tiny effect corresponds to the second order multipole expansion of the gravitational potential of the earth, where ${J}_{2}$ is the coefficient quantifying the oblateness of the earth, and it also includes a periodic component, with period equal to one half of the orbital period of the satellite. A review of the J2 relativistic effect can be found in Kouba (2004).

The J2 periodic effect is expected to be clearly visible on the apparent frequency of space clocks such as the Passive Hydrogen Masers (PHMs) flying on board the Galileo satellites, thanks to their impressive stability, as already shown in Formichella et al. (2019) and in Kouba (2019). For this reason, Steigenberger et al. (2015) correct the periodic component of the J2 effect before estimating the stability of the Galileo clocks. Moreover, the J2 periodic effect can be also detected on the Rubidium Atomic Frequency Standards (RAFSs) flying on board the GPS satellites, as shown in Formichella (2018) and in Senior et al. (2008).

As discussed in Senior et al. (2008) for GPS and in Waller et al. (2010) for Galileo, orbital estimation errors are another source of periodic frequency variations. Indeed, orbital and clock estimates are products of the same algorithm, so that an error in estimating the satellite position is translated into an error on the estimated time offset of its on-board clock. Periodic orbital estimation errors are due, for instance, to the implementation of non-optimal Solar Radiation Pressure (SRP) models, and their period is expected to correspond to the orbital period of the satellite. Such periodic orbital estimation errors are then translated into periodic variations of the estimated time offset, and hence into periodic variations of the apparent frequency. While the first harmonic of this periodic signal has a period corresponding to the orbital one, its second harmonic has the same period of the J2 signal, and this makes more difficult to obtain an unbiased estimate of the J2 signal amplitude. However, the magnitude of the orbital estimation errors, and hence the amplitude of the corresponding periodic frequency variations, changes in time as a function of the angle between the sun and the orbital plane of the satellite. Therefore, in order to lower the bias on the estimate of the J2 signal amplitude, one can find the period of the year during which the orbital estimation errors are smaller and then perform the J2 amplitude estimation over such period.

Another source of periodic variations is the sunlight thermal effect. Indeed, the frequency of any atomic clock is sensitive to temperature variations, as those experienced by the satellite depending on its orientation with respect to the sun. Such orientation varies periodically along the orbit, therefore the induced frequency variations are expected to have a period equal to the orbital one. The sunlight thermal effect is usually mitigated by the action of an on-board temperature control system, but we cannot exclude the presence of some residual periodic effect on the apparent clock frequency, also due to on-board phase delay variations caused by temperature variations, as conjectured in Waller et al. (2010).

We analyze Galileo space clock data provided by the European Space Agency (ESA) with two main objectives: first, to characterize the most relevant periodic variations affecting the apparent clock frequency, and to clearly show them by exploiting the effective visual rendering of time–frequency analysis; second, to estimate the amplitude of the J2 periodic signal for a comparison with its theoretical value. Apart from being a further test of general relativity, which, however, has been proven with much higher accuracy by Delva et al. (2018) and Herrmann et al. (2018) with another experiment using Galileo satellites, this second objective represents a specific validation of the J2 correction that could be applied at the user level, to improve the timekeeping and positioning performance of Galileo and the other GNSSs.

After a review of the J2 relativistic effect, we present the data analysis and the experimental results. A final discussion of the results and of the further developments of our work is given in the conclusions.

The J2 relativistic effect

After a short introduction on relativistic effects in GNSS, we review the theory of the J2 relativistic effect and we discuss and fix some ambiguities which can be found in the literature. Then, we compute the expected amplitude of the periodic component of the J2 effect, in the case of Galileo satellites on a nominal orbit, and we give some insights on what happens for satellites with a higher-eccentricity orbit.

According to the general theory of relativity, in the weak-field approximation, the ratio between the infinitesimal intervals of proper time measured by a ground and a satellite clock is

$$\frac{{\text{d}}{\tau }_{\mathrm{G}}}{{\text{d}}{\tau }_{\mathrm{S}}}=\frac{1}{{c}^{2}}\left(1+{\phi }_{\mathrm{G}}-{\phi }_{\mathrm{S}}-\frac{{v}_{\mathrm{G}}^{2}}{2}+\frac{{v}_{\mathrm{S}}^{2}}{2}\right)$$

(1)

where the subscripts $G$, $S$ stand for ground and satellite, respectively, $c$ is the speed of light, $\phi$ is the gravitational potential, and $v$ is the speed measured in the Earth Centered Inertial (ECI) frame. The gravitational potential at a distance $r$ from the center of the earth can be written as

$$\phi \left(\overrightarrow{r}\right)=-\frac{\mathrm{GM}}{r}+R\left(\overrightarrow{r}\right)$$

(2)

where $G$ is the universal gravitational constant, $M$ is the mass of the earth, and $R$ is a perturbing potential accounting for the fact that the earth is not a perfect sphere, and including the gravitational effect of other celestial bodies as the sun and the moon. According to Kouba (2004), the main contribution to the perturbing potential comes from the oblateness of the earth and hence, hereafter, we neglect the other contributions. With such an approximation and by integrating (1), following the derivation of Kouba (2004), we can obtain the time delay between a ground and a satellite clock, $\Delta {t}_{\mathrm{GS}}={\tau }_{G}-{\tau }_{S}$, accumulated during a time interval $\Delta {\tau }_{G}$ measured by the ground clock,

$$\Delta {t}_{\mathrm{GS}}=\left(\frac{3\mathrm{GM}}{2{\mathrm{a}}_{0}{\mathrm{c}}^{2}}-{L}_{\mathrm{G}}\right)\Delta {\tau }_{\mathrm{G}}+\frac{2\overrightarrow{\mathrm{r}}\bullet \overrightarrow{\mathrm{v}}}{{c}^{2}}+\Delta {t}_{{\mathrm{J}}_{2}}$$

(3)

where ${a}_{0}$ is the mean semi-major axis of the satellite orbit, ${L}_{\mathrm{G}}=$ 6.969290134 × 10^–10 by definition, assuming that the ground clock is at rest on the geoid and according to IERS (2010), $\overrightarrow{\mathrm{r}}$ and $\overrightarrow{\mathrm{v}}$ are the position and the velocity of the satellite, and $\Delta {t}_{{\mathrm{J}}_{2}}$ is the contribution of the J2 relativistic effect, given by

$$\Delta {t}_{{\mathrm{J}}_{2}}=\frac{{a}_{\mathrm{E}}^{2}{J}_{2}}{2{\mathrm{a}}_{0}^{2}{\mathrm{c}}^{2}}\left[\frac{7\mathrm{GM}}{{a}_{0}}\left(1-\frac{3}{2}{\mathrm{sin}}^{2}i\right)\Delta {\tau }_{\mathrm{G}}+3\sqrt{{\mathrm{GMa}}_{0}}{\mathrm{sin}}^{2}i\mathrm{sin}\left(2u\right)\right]$$

(4)

where ${a}_{\mathrm{E}}$ is the equatorial radius of the earth, $i$ is the inclination of the satellite orbit and $u$ is the angular position of the satellite along the orbit. The coefficient ${\mathrm{J}}_{2}$ quantifies the oblateness of the earth and, according to IERS (2010), its value is 1.0826359 × 10^–3 with uncertainty 1 × 10^–10. The first terms of (3) and (4), both proportional to $\Delta {\tau }_{\mathrm{G}}$, can be interpreted as the result of a constant frequency offset between the two clocks. The second term of (3) is instead a periodic component, with period equal to the orbital one and amplitude proportional to the orbital eccentricity. Finally, the second term of (4) is the periodic component of the J2 relativistic effect, with period equal to one half of the orbital one, and amplitude depending on the oblateness of the earth as well as on the inclination of the satellite orbit.

Equation (4) is well known although some ambiguities can be found in the literature, which we now discuss. First of all, the result reported in Kouba (2004) has a sign error originating in  (28), regarding the constant frequency offset term. This is also confirmed by a note at page 155 of IERS (2010). Once the sign is fixed, the final result reported in  (31) of Kouba (2004) is fully in agreement with (4), apart from an overall minus sign due to the fact that we compute ${\tau }_{\mathrm{G}}-{\tau }_{\mathrm{S}}$ instead of ${\tau }_{\mathrm{S}}-{\tau }_{\mathrm{G}}$. The same sign error seems to be present in  (6) of Kouba (2019), which indeed refers to Kouba (2004). Finally, there is a coefficient error in  (25) of Nelson (2011), again regarding the constant frequency offset term, which should be multiplied by a factor 7 to give the correct result compatible with (4).

The first two terms of (3) are the main part of the total relativistic effect and must be corrected, in order to avoid positioning errors up to several kilometers, whereas the J2 relativistic effect is so small that it is usually neglected. In the case of GPS, according to the Interface Specification document IS-GPS-200K (2019), the constant frequency offset, first term of (3), is compensated at the satellite level by applying an hardware offset to the output of the on-board clock, whereas the periodic component, second term of (3), is corrected at the user level by the receiver software. In the case of Galileo, according to the Interface Control Document OS-SIS-ICD (2016) and to Mudrak et al. (2015), the periodic component is corrected by the receiver software as in GPS, whereas it is not specified if the constant frequency offset is corrected at the satellite-hardware level or not. However, even if no hardware correction is applied to a specific satellite clock, the constant frequency offset is automatically included in the clock correction parameters broadcasted to the users within the navigation message, so that it may be corrected at the user-receiver level. Finally, both GPS and Galileo do not apply any special correction for the J2 relativistic effect. This is actually not a problem for the constant frequency offset component of the J2 effect, as such an offset is automatically included in the clock correction parameters, whereas the periodic component shows up in the apparent clock frequency and affects the position estimation, unless compensated by the user in post-processing by using the second term of (4) as a time correction.

We now focus on the J2 periodic component only, $\Delta {t}_{{\mathrm{J}}_{2}}^{\mathrm{per}}$. From (4) we have

$$\Delta {t}_{{\mathrm{J}}_{2}}^{\mathrm{per}}=A\mathrm{sin}\left(2u\right)$$

(5)

which is a sinusoidal signal with amplitude

$$A=\frac{3\sqrt{\mathrm{GM}} {J}_{2}{a}_{\mathrm{E}}^{2}{\mathrm{sin}}^{2}i}{2{\mathrm{a}}_{0}^{3/2}{c}^{2}}$$

(6)

affecting the apparent time offset of the satellite clock measured with respect to a stable ground clock. The fractional frequency offset of a clock, $y\left(t\right)$, is the time derivative of its time offset. Therefore, the J2 signal affecting the apparent fractional frequency offset of the satellite clock can be obtained by differentiating (5), namely,

$${y}_{{\mathrm{J}}_{2}}^{\mathrm{per}}=B\mathrm{cos}\left(2u\right)$$

(7)

which is a sinusoid with amplitude

$$B=2A\frac{{\text{d}}u}{{\text{d}}t}$$

(8)

To the first-order approximation (circular orbit), the time-dependent angular position of the satellite along the orbit can be written as

$$u\left(t\right)=\frac{2\pi t}{T}$$

(9)

where the orbital period $T$, to the same order of approximation, is given by

$$T=\frac{2\pi {\mathrm{a}}_{0}^{3/2}}{\sqrt{\mathrm{GM}}}$$

(10)

leading to

$$B=\frac{3\mathrm{GM}{J}_{2}{a}_{\mathrm{E}}^{2}{\mathrm{sin}}^{2}i}{{\mathrm{a}}_{0}^{3}{c}^{2}}$$

(11)

which is the theoretical value of the amplitude of the J2 signal affecting the apparent satellite clock frequency.

Putting numbers in (6) and (11) gives the expected numerical value of the J2 signal amplitude. The values of the GNSS-independent numerical constants are reported in Table 1, taken from Mohr et al. (2016) for the fundamental physical constants and from The Astronomical Almanac 2018 (2017) for the earth parameters. The Galileo system orbital parameters are reported in Table 2, taken from Kouba (2019), where the column value gives the nominal parameters for the Galileo constellation. However, once sent into orbit, each satellite has its own parameters with a small deviation from the nominal ones, which can be treated as an uncertainty associated to the nominal parameters. Such an uncertainty is inferred from the data reported in Kouba (2019), where the mean semi-major axis and the inclination are estimated for each Galileo satellite by using MGEX products for November 23, 2018 and February 15, 2019. For each of these two epochs, we compute the average and the standard deviation of the differences between the nominal parameters and the measured once, including all the Galileo satellites except GSAT0201 and GSAT0202, which are on a higher-eccentricity orbit. Then, we give a conservative estimate of the uncertainty associated to the nominal parameters by directly summing the absolute value of the average and the standard deviation computed at the previous step. We do this for the two epochs, thus obtaining two conservative estimates of the uncertainty for each of the parameters and, finally, we select the largest one to be used in the following calculations. Such uncertainties are those reported in Table 2.

Table 1 Values of fundamental physical constants and earth parameters required to compute the J2 signal amplitude, along with their 1-sigma uncertainties

Table 2 Nominal values of Galileo system orbital parameters required to compute the J2 signal amplitude, along with uncertainties inferred from data reported in Kouba (2019)

In order to associate an uncertainty to the expected amplitude of the J2 signal, we apply the law of uncertainty propagation to (6) and (11), as recommended by the guide to the expression of uncertainty in measurement, GUM (2008), and by assuming uncorrelated uncertainties for the quantities listed in Tables 1 and 2. The results are reported in Table 3. Note that, despite the conservative uncertainty associated to the orbital parameters, the final relative uncertainty affecting the amplitude of the J2 signal is as low as 2%, so that we can consider the values reported in Table 3 to hold for the whole Galileo constellation, except GSAT0201 and GSAT0202, which we discuss at the end of this section. We also stress the fact that the uncertainties reported in Table 3 are completely dominated by the uncertainty on the inclination angle, accounting for approximately the 99.999% of the total uncertainty on the theoretical value of the J2 amplitude.

Table 3 Amplitude of the J2 periodic signal, along with 1-sigma uncertainties, valid for all the Galileo satellites except GSAT0201 and GSAT0202

According to the values reported in Table 3, the J2 signal affecting the Galileo apparent clocks gives a time error up to about 60 ps. Remembering that, as a rule of thumb, the position error induced by a time error ${\varepsilon }_{t}$ is at the level of $c\bullet {\varepsilon }_{\mathrm{t}}$, we find that the position error induced by the J2 signal is up to about 2 cm, which could be relevant for GNSS high-accuracy applications. Therefore, it is worth studying the J2 periodic component and validating its theoretical formula, which could then be used to correct the J2 periodic effect at the user level, improving the navigation solution.

For a Galileo satellite with nominal parameters, the orbital period $T$ is of about 14 h, so that the period of the J2 signal, ${T}_{{\mathrm{J}}_{2}}$, is of about 7 h. According to the specifications of the Galileo clocks reported in Waller et al. (2010), the frequency stability of a PHM at an averaging time equal to the period of the J2 signal, ${\sigma }_{y}\left({T}_{{\mathrm{J}}_{2}}\right)$, is about 1 × 10^–14. Compared to the value of $B$, this number tells us that the J2 periodic signal should be visible on the apparent frequency of the Galileo PHM clocks, as we are going to demonstrate in the following section.

Among all the in-orbit Galileo satellites, GSAT0201 and GSAT0202 have been deployed on a high-eccentricity non-nominal orbit due to a problem upon launch, as described in Steigenberger and Montenbruck (2017). In principle, Eqs. (6) and (11) cannot be applied to these satellites, as they derive from (4), which has been computed under a set of assumptions. Specifically, following the derivation of Kouba (2004), the eccentricity of the satellite orbit is neglected consistently with the declared order of approximation. This can be safely done for Galileo satellites on a nominal orbit, with a typical eccentricity $e\hspace{0.17em}<\hspace{0.17em}$0.0005, but not for GSAT0201 and GSAT0202, whose orbits have an eccentricity of about 0.166, as reported in Kouba (2019). However, in the same paper, Kouba compares the results obtained from the analytical formula (6) and from the numerical integration of (1), computed at the specific epoch of November 23, 2018, founding a good agreement between the two results despite the large eccentricity. Nonetheless, in the absence of an analytical formula for the satellites with high eccentricity, we cannot say if such an agreement is always verified, especially considering that the argument of perigee of the orbits of GSAT0201 and GSAT0202 changes in time at a rate of about 0.034 degrees per day, as reported by the European GNSS Agency (GSA) at https://www.gsc-europa.eu/system-service-status/orbital-and-technical-parameters. Therefore, we think it is worth deriving an analytical formula for satellites with high eccentricity too. In case of need, such an analytical formula would be much easier than a numerical integration to be implemented as a relativistic correction at the user-receiver level.

Data analysis and experimental results

In this section, we analyze Galileo space clock data, provided by ESA in form of phase offsets between the satellite clocks and a ground reference clock. In particular, we analyze one year of data from three PHMs acting as master clock on three Galileo satellites, flying on two different orbital planes. The chosen year is 2017 and the chosen satellites are listed in Table 4, along with relevant information retrieved from the GSA web page. Note that, hereafter, we will refer to a satellite and/or its on-board clock by its Pseudo-Random Noise (PRN) number, that is, by the satellite ID as listed in Table 4.

Table 4 List of analyzed Galileo satellites and on-board clocks, along with relevant information

The analysis is based on three different techniques, namely frequency analysis, time–frequency analysis and data fitting, which are applied to properly preprocessed data. In particular, we use time–frequency analysis to characterize the time variations of the observed spectral components, whereas data fitting is used to estimate the amplitude of the J2 effect for a comparison with its theoretical value.

A brief introduction to time–frequency analysis

Time–frequency analysis is a multidisciplinary field for the representation of signals whose frequency content changes with time, see Cohen (1995) for a thorough introduction to the subject, Amin et al. (2017) for applications to GNSSs, and Galleani (2018) for applications to atomic timing. Contrary to frequency analysis, where the frequency content of a signal is obtained through its Fourier transform, in time–frequency analysis there are infinite ways to obtain the time–frequency spectrum of a signal, referred to as time–frequency distribution or time–frequency representation.

An effective and widely used time–frequency distribution is the spectrogram, defined as

$${P}_{y}(t,f)={\left|\underset{-\infty }{\overset{+\infty }{\int }}h\left({t}^{^{\prime}}-t\right)y\left({t}^{^{\prime}}\right){e}^{-i2\pi f{t}^{^{\prime}}}{\text{d}}{t}^{^{\prime}}\right|}^{2}$$

(12)

where $y(t)$ is the analyzed signal, such as the fractional frequency offset of a clock, and $h(t)$ is the analysis window, typically a Hamming, Hann, or rectangular window. The spectrogram is the squared magnitude of the short-time Fourier transform, and it therefore describes how the energy of the time-varying frequencies of the signal changes with time.

To better illustrate the basic ideas and the effectiveness of time–frequency analysis, we consider the test signal $y(t)$ shown in Fig. 1, where only a generic noisy behavior can be seen. Conversely, the spectrogram in Fig. 2 clearly reveals that the test signal has a rich nonstationary structure, made by four components. The component 1 is a sinusoid at the normalized frequency $f\hspace{0.17em}=\hspace{0.17em}$0.05, with constant amplitude and lasting for the entire duration of the signal. Note that its period is 1/0.05 $=\hspace{0.17em}$20 samples. The component 2 is a short-duration sinusoid with a Gaussian amplitude. The component 3, namely, the noisy pattern between $f\hspace{0.17em}=\hspace{0.17em}$0.1 and $f\hspace{0.17em}=\hspace{0.17em}$0.3, is a bandpass noise, obtained by bandpass filtering a white Gaussian noise. The component 4 is a linear chirp, namely, a sinusoid whose frequency changes linearly with time. The amplitude of this chirp varies with time, with a sinusoidal behavior. Since this chirp is corrupted by the bandpass noise, the amplitude of its time–frequency representation fluctuates with respect to time.

Fig. 1

Test signal

Fig. 2

Spectrogram of the test signal shown in Fig. 1

Frequency and time–frequency analysis of the Galileo satellite clocks

We now investigate the frequency and time–frequency structure of the 2017 apparent frequency of the Galileo satellite clocks listed in Table 4.

We start by preprocessing the phase offset data provided by ESA with a five-step procedure. First, we select a stable ground reference clock, possibly an Active Hydrogen Maser (AHM), trying to minimize the changes of reference during the analysis period. Second, to obtain the frequency offset data that we will analyze in the frequency and time–frequency domains after the preprocessing operations described in the subsequent steps, we follow the standard definition (see, for example, Kartaschoff (1978)) and differentiate the phase offset data with a sampling time ${T}_{\mathrm{s}}=$ 300 s. Note that the differentiation process translates the random walk noise affecting the phase data into a white frequency noise, and the observed day-boundary phase jumps, resulting from the data processing applied by the analysis centers generating the clock products, into frequency outliers. Third, we filter out these and the other possible frequency outliers according to the approach described in Sesia et al. (2016). Fourth, we estimate the time-varying mean $\widehat{\mu }(t)$ of the frequency offset $y(t)$ by sliding a window made by N_W = 169 samples. By multiplying by the sampling time ${T}_{\mathrm{s}}$, we find that this window length corresponds to 50,700 s, that is, to 14 h and 5 min, a good approximation of the 14 h, 4 min and 45 s orbital period declared by Subirana et al. (2013) for the Galileo satellites on a nominal orbit. Therefore, this window length guarantees that the periodic signals due to the orbital estimation errors and to the J2 relativistic effect are filtered out in $\widehat{\mu }(t)$. Then, we detrend the frequency offset $y(t)$ by subtracting $\widehat{\mu }(t)$. Note that the detrended frequency still contains the periodic signals of interest, because they are not present in $\widehat{\mu }(t)$, whereas the linear frequency drift typical of the PHM and any other possible long-term non-stationarity of the clock frequency are removed from $y\left(t\right)$. Finally, in the last step, we replace with zero values the data blocks with two or more missing values, and we interpolate the individual missing data.

Only two different ground references have been used for the whole analysis period, both of them being Galileo Experimental Sensor Stations (GESSs) connected to high-precision atomic clocks. In particular, for the first six months of the year 2017, the reference is GIEN, that is, the GESS hosted at the Italian National Institute of Metrological Research (INRiM), which is connected to the local realization of the Universal Time Coordinated (UTC), namely UTC(IT), a time scale generated by an AHM steered toward UTC. The data provided by ESA are corrected for the first order relativistic effect for the whole analysis period, but not for the J2 relativistic effect until the month of November 2017.

Figure 3 shows the pre-processed frequency offset for all of the three considered satellites. As it can be observed from the small data gaps, the data coverage is very good, being approximately 94% for E24 and 98% for E30 and E09.

Fig. 3

Preprocessed fractional frequency offset for 2017 of the PHMs on board the E24 (top), E30 (middle), and E09 (bottom) Galileo satellites

Subsequently, we perform a frequency analysis of the clock data. Figure 4 shows the energy spectrum ${S}_{\mathrm{y}}(f)$ of the frequency offset for the three satellites, obtained as the squared magnitude of the Fourier transform of $y\left(t\right)$. All of the three spectra clearly highlight the presence of two peaks at low frequencies, corresponding to periods of approximately 14 and 7 h, that is, to the first and second harmonics of a signal with period equal to the orbital one. Such periodic signal is expected to be due to the orbital estimation errors and, possibly, to residual sunlight thermal effects. However, the peak corresponding to the second harmonic is also due to the sinusoidal component of the J2 relativistic effect. Only for E24 the peak corresponding to the first harmonic is smaller than that of the second harmonic. In the middle frequencies, and particularly in the bandwidth 0.2 ≤ f ≤ 0.3, for all of the three satellites we also notice several peaks whose amplitudes seem to follow a noisy pattern.

Fig. 4

Frequency spectrum for the E24 (top), E30 (middle), and E09 (bottom) clock data shown in Fig. 3

We finally perform a time–frequency analysis of the PHM clock data for the three satellites. Figures 5, 6, 7 show the spectrogram of the data reported in Fig. 3, computed with a window length of 1 month. The first and second harmonics are clearly visible in the low frequencies. We immediately notice that the amplitude of the harmonics varies with time, as well as that this variation is not identical for each of the three clocks. In the first six months of the year and in the middle frequencies, we also notice several components with time-varying amplitude. These components correspond to the peaks with random amplitude located in the middle frequencies of the energy spectrum shown in Fig. 4. In the second half of the year, more precisely in the month of October and approximately for $f=0.42$, we see a component similar to a short-duration sinusoid with an amplitude modulation. We believe that both the noise-like components in the middle frequencies and this short-duration sinusoid are due to the ground reference used to measure the clock data, which actually changes on July 1, exactly when the noise-like components in the middle frequencies disappear.

Fig. 5

Spectrogram of the E24 clock data shown in Fig. 3, for the year 2017

Fig. 6

Spectrogram of the E30 clock data shown in Fig. 3, for the year 2017

Fig. 7

Spectrogram of the E09 clock data shown in Fig. 3, for the year 2017

To better investigate the behavior of the first and second harmonics, as well as their relationship with the satellite orbit, in Figs. 8, 9, 10 we show the spectrogram of the 2017 clock data for the three satellites in the frequency interval $0\le f\le 0.017$. The variation with time of the amplitude of the two harmonics is now clearly visible for each of them. The blue curve in the left side panel is the $\beta$-angle of the satellite as a function of time, that is, the angle between the sun and the orbital plane of the satellite. When $\beta$ is in the ± 12.4° range, whose bounds are indicated by the red dashed curves, the satellite is in an eclipse season, during which the incident sunlight is eclipsed by the earth for some time each day. Eclipse seasons are indicated by the solid red curves. We see that during the eclipse seasons the amplitude of the first harmonic reaches its maximum values. In the case of E30, a local maximum is also reached before the second eclipse season. Conversely, the amplitude of the first harmonic reaches its minimum values when the absolute value of $\beta$ is at its maximum. This is in agreement with what has already been reported by Steigenberger et al. (2015) and Steigenberger and Montenbruck (2017). Such behavior is likely due to the use of non-optimal SRP models, to modeling deficiencies linked to the limited knowledge of the attitude of the satellite during eclipses, and to possible thermal issues when the satellite passes from the sunlight to the shadow of the earth or the other way around.

Fig. 8

Spectrogram of the E24 clock data shown in Fig. 3, zoom on the low-frequency region of Fig. 5

Fig. 9

Spectrogram of the E30 clock data shown in Fig. 3, zoom on the low-frequency region of Fig. 6

Fig. 10

Spectrogram of the E09 clock data shown in Fig. 3, zoom on the low-frequency region of Fig. 7

In addition, Figs. 8–10 show that, for all of the three satellites, the amplitude of the second harmonic becomes very small in the last two months of the year. To better highlight this phenomenon, in Fig. 11 we show the spectrogram of E24 starting from October 15 only, computed with a window length of 3 days. The picture reveals that the second harmonic basically disappears from the date of November 1, indicated by the solid red line. To understand this behavior, we have to recall two facts: first, that what we generically indicated as “second harmonic” in Figs. 8, 9, 10, it is not just the second harmonic of the periodic signal due to the orbital estimation errors and the possible residual thermal effects, but it is rather the superposition of such second harmonic and of the sinusoidal signal due to the J2 relativistic effect; second, that starting from November 2017 the clock data provided by ESA are corrected for the J2 effect. Therefore, the reason for the observed behavior is that the J2 correction applied at the beginning of November actually removed the J2 signal from the clock data, whereas any residual signal visible at a period of 7 h is probably the “pure” second harmonic of the periodic signal at 14 h. Such an observation is a first validation of the J2 correction computed according to (5) and (6), which could also be applied at the user level. A second validation is given by the analysis results presented in the next section, based on data fitting for a direct estimation of the J2 signal amplitude.

Fig. 11

Spectrogram of the E24 clock data shown in Fig. 3 starting from October 15, 2017, zoom on the low-frequency region, highlighting the disappearance of the second harmonic

Data fitting and J2 signal amplitude estimation

We now estimate the amplitude of the J2 signal by fitting the available clock data with a model function and compare the obtained estimates with the theoretical value reported in Table 3.

We start by preprocessing the phase offset data provided by ESA with a similar procedure with respect to the frequency and time–frequency analysis. As in the previous case, we select a stable ground reference clock, we differentiate the phase offset data to obtain frequency offset data, and we filter out the frequency outliers. The next steps of the preprocessing are a little bit different with respect to the previous analysis. First of all, in this case we work with selected non-overlapping data batches with a typical duration of one month, discarding those batches showing severe apparent clock non-stationarities. Then, for each selected batch, we independently estimate and remove a linear frequency drift. Finally, only if it is strictly necessary because of the presence of stochastic long-term fluctuations, we estimate and remove the time-varying mean $\widehat{\mu }(t)$, computed as in the previous analysis, by using a sliding window of ${N}_{\mathrm{w}}=169$ samples or integer multiples of such value.

The model function is the sum of two independent sinusoids, namely,

$$F\left(t\right)={A}_{1}\mathrm{sin}\left(\frac{2\pi }{{T}_{1}}t+{\varphi }_{1}\right)+{A}_{2}\mathrm{sin}\left(\frac{2\pi }{{T}_{2}}t+{\varphi }_{2}\right)$$

(13)

where the two amplitudes ${A}_{\mathrm{1,2}}$, the two periods ${T}_{\mathrm{1,2}}$, and the two phases ${\varphi }_{\mathrm{1,2}}$ are all free parameters. The objective is to fit the two sinusoids corresponding to the two harmonics discussed in the previous section, with periods $T$ and ${T}_{{\mathrm{J}}_{2}}$, so that we can estimate their amplitudes. By fitting the frequency data in the selected batches, we obtain a set of independent estimates of ${T}_{\mathrm{1,2}}$ and ${A}_{\mathrm{1,2}}$, and we further discard those batches for which none of the estimated periods is compatible with ${T}_{{\mathrm{J}}_{2}}$. Furthermore, in case both the estimated periods are compatible with ${T}_{{\mathrm{J}}_{2}}$, we repeat the fitting with a simplified 1-sinusoid model function. Therefore, taking only the results from the remaining batches, we obtain a set of independent and biased estimates of the amplitude of the J2 signal, where the bias is due to the “pure” second harmonic of the periodic signal caused by the orbital estimation errors and the residual thermal effects. We assume the presence of such bias to be the main reason why, in Figs. 8, 9, 10, the so called “second harmonic” is not a signal with constant power, as it should have been if it was only due to the J2 periodic effect.

An example of a successfully fitted 1-month data batch is given in Fig. 12 for the satellite E30, where we only show the first 5 days of the batch, in order to make the periodic signals visible by eye.

Fig. 12

Preprocessed fractional frequency offset of E30 for the first five days of June 2017 (blue), along with the regression curve obtained by using the model function specified in (13) (red)

After analyzing the three satellites for the whole year 2017, we end up with a set of 29 biased estimates of the J2 signal amplitude $B$, specifically, 11 for E24, 8 for E30, and 10 for E09. Figure 13 shows the results for E24 and E30, which are plotted together as these two satellites are on the same orbital plane and, therefore, they share the same $\beta$-angle and eclipse seasons. Analogously, we show in Fig. 14 the results for E09, which is on a different orbital plane. In both the figures, we highlight the eclipse seasons as shaded regions, indicate with a vertical dashed line the epoch from which the J2 effect is corrected in the input data, and indicate with horizontal lines the theoretical value of $B$ with a 3-sigma confidence interval. The estimates are also plotted with 3-sigma uncertainty bars, as obtained from the fit. However, note that such bars represent the statistical uncertainty only, that is, they do not include the unknown bias and are therefore underestimated. Moreover, note that the lack of estimates for the months of November and December, for all of the three satellites, is a direct consequence of the J2 correction applied to the input data.

Fig. 13

Biased estimates of the amplitude of the J2 signal for E24 (blue dots) and E30 (green triangles), along with the theoretical value of the J2 signal amplitude (red line)

Fig. 14

Biased estimates of the amplitude of the J2 signal for E09 (blue dots), along with the theoretical value of the J2 signal amplitude (red line)

We now focus on the bias affecting the estimates of $B$. Disentangling the effects of the J2 signal and of the second harmonic due to orbital estimation errors and residual thermal effects, the latter being the bias, is very difficult, and hence it is difficult to estimate the bias itself. However, assuming that the regions in which the first harmonic is at its minimum are those in which the orbital estimation errors and the thermal effects are also at their minimum, we can infer that the estimates of $B$ obtained in the same regions are those with the smallest bias. For instance, by observing Fig. 8, we find that the first harmonic of E24 reaches its minimum in the months of May, June, September and October. Actually, if we look at Fig. 13, we find that the estimates corresponding to the months of June and October are among the closest to the theoretical value. Conversely, the estimate corresponding to the month of July is one of the most far from the theoretical value, consistently with the fact that the first harmonic has a local maximum in the month of July, during one of the eclipse seasons. The estimates obtained from E30 follow a similar behavior with respect to that of E24, as it can be expected since the two satellites are on the same orbit. In particular, from Fig. 9 we find that the first harmonic of E30 reaches its minimum in the months of March, April, September and October, and indeed, by looking at Fig. 13, we find that the estimates corresponding to the months of March and September are among the closest to the theoretical value. Finally, looking at Fig. 10, we find that the first harmonic of E09 reaches its minimum in the month of June, which is exactly the one giving the estimate of $B$ closest to the theoretical value, as it can be observed from Fig. 14.

As a consequence of such observations, we try to give a single, improved estimate of $B$ for each of the satellites, by computing the weighted average of the estimates, which we inferred to have the lowest bias. The final results are graphically illustrated in Fig. 15 and are all well compatible with the theoretical value, the most far point being less than 2 sigma apart. Note that the uncertainty of the estimates is plotted as 1-sigma error bars, and the shaded region around the theoretical value (middle horizontal line) also represents a 1-sigma confidence interval. Moreover, the relative uncertainty of the three estimates is also pretty low, being approximately the 3% for E24 and E30, and the 6% for E09. Such results represent our final validation of the J2 signal amplitude computed according to (6), and hence of the J2 relativistic correction computed according to (5) and (6).

Fig. 15

Final estimates of the J2 signal amplitude for the three satellites (blue dots), along with its theoretical value (red line)

Conclusions

We have analyzed one year of clock phase offset data from three PHMs flying on board the Galileo satellites in 2017. The power spectrum of the apparent clock frequency reveals the presence of a periodic signal with a period equal to the orbital one, that is, about 14 h, primarily due to orbital estimation errors and possible residual thermal effects, and of which the first two harmonics are clearly visible. The peak corresponding to the second harmonic is also due to the presence of another sinusoidal signal, that is, the periodic component of the J2 relativistic effect, which has a period equal to one half of the orbital one.

Then, we applied time–frequency analysis to characterize the time variations of the power spectrum, and we verified that the amplitude of the first harmonic is, in general, lower for higher values of the $\beta$-angle, whereas it is particularly high during the eclipse seasons, as expected. We also estimated the amplitude of the J2 signal, and we found it in good agreement with the theoretical value, as the best estimate is less than 1 sigma apart from it and with a relative uncertainty as low as 3%. This represents a validation of (5) and (6), which can be used as a further relativistic correction to be implemented at the user level, for improving GNSS positioning and timing solutions, especially when a centimeter-level position accuracy is required. Such correction has been actually applied to the data provided by ESA, for the last two months of the analyzed period, and indeed our analysis confirmed that the J2 signal had been successfully removed from the observed clock frequency during such months. This is particularly evident from the results of our time–frequency analysis.

Regarding the future developments of our work, we aim at extending the analysis to all the other Galileo satellites and to other time periods. Moreover, a separate study will be devoted to the Galileo satellites GSAT0201 and GSAT0202, flying on a higher-eccentricity orbit and hence requiring a dedicated assessment of the J2 relativistic effect.

Finally, we comment on the possible application of our analysis as an alternative method for estimating the ${J}_{2}$ coefficient. Indeed, (11) can be inverted to find ${J}_{2}$ as a function of $B$, the latter being experimentally estimated as described in the previous section. However, considering as an example the value of $B$ obtained from E24, which has the lower uncertainty, even if we neglect the uncertainties of all the other parameters, the resulting uncertainty associated to the estimate of ${J}_{2}$ is orders of magnitude worse than that obtained with classical geophysical methods. Therefore, such kind of application is currently not of any practical interest, as it is limited by the uncertainty associated to the experimental estimate of $B$.