Preliminary experimental results of determining the geopotential difference between two synchronized portable hydrogen clocks at different locations

As early as in 1980, it was proposed that the CVSTT technique could be used for comparing clock offsets between two remote clocks (Allan and Weiss,, 1980; Allan et al.,, 1985). This technique was adopted by the Bureau International des Poids et Mesures (BIPM) as one of the main methods for transferring international atomic time (TAI) signals (Allan and Thomas,, 1994; Imae et al.,, 2004). Using this method, the uncertainty of comparing remote time may achieve several nanoseconds (Lewandowski and Thomas,, 1991; Lewandowski et al.,, 1999; Ray and Senior,, 2003; Sun et al.,, 2010; Yang et al.,, 2012; Rose et al.,, 2016). The main advantages of CVSTT technique lie in that the satellite clock errors are cancelled, and various other errors, especially ionospheric and tropospheric errors, are largely reduced due to the simultaneous two-way observations (Allan and Weiss,, 1980; Yang et al.,, 2012; Defraigne and Baire,, 2011).

In this study, we use the CVSTT technique for comparing the clock offsets. The reasons are stated as follows. Firstly, the CVSTT technique is a stable and accurate method for comparing clock offsets between two remote clocks, and this technique has been adopted by the Bureau International des Poids et Mesures (BIPM) as one of the main methods for transferring international atomic time (TAI) signals. Therefore, the results of CVSTT technique are reliable and stable. Secondly, using CVSTT technique is relative cheap and convenient in practice. Though the experimental results using optical atomic clocks via fiber links are quite accurate and stable indeed, the fiber-link is limited to a relative short distance and fixed stations with fiber links, and for a long distance even intercontinental comparison, the fiber-link is extremely challengeable and costly. In addition, the Two-way satellite time and frequency transfer (TWSTFT) technique is more accurate and stable than CWSTFT technique, but it is very expensive in economy, which could be applied in the future when the condition is available. Lastly, though the carrier phase measurements in GNSS are more precise than the corresponding precise code measurements in time transfer, the absolute clock offsets cannot be determined from carrier phase measurements due to the ambiguities issues (Defraigne and Bruyninx,, 2007). Therefore, the code measurements are necessary for determining the absolute clock offsets between remote time and frequency standards for a long time. Hence, at present experiment condition, the CVSTT technique and precise code measurements are adopted here.

To better understand the text, here we briefly introduce the principle of the CVSTT technique. Suppose there are two ground stations ($A$, $B$) and a common-view satellite $S$. The coordinates of the two ground stations are denoted as $(x_{i},y_{i},z_{i})$ ($i=A,B$), and the satellite coordinates are denoted as ($x_{{}_{S}},y_{{}_{S}},z_{{}_{S}}$). At an appointed time, the $i$-th staion’s clock and satellite’s clock emits 1 purse per second (1 PPS) signals, simultaneously. The 1 PPS signal of the $i$-th station (denoted as 1PPS-local signal) arrives to the GNSS time transfer receiver (GNSS-TTR) via cables, and opens the door of the time interval counter (TIC) for starting timing. At the same time, the 1PPS signal of the satellite $S$ (denoted as 1PPS-S signal) is received by the GNSS antenna of the ground station. After signals demodulation the 1PPS-S signal arrives to the GNSS-TTR via cables. And the 1PPS-S signal is used for closing the TIC’s door for stopping timing. The $i$-th staion’s clock also generates the 10 MHz sinusoidal signals with extremely high accuracy. Therefore, for the $i$-th station and at a certain timepoint, by counting the number of the sinusoidal waves between the starting timepoint (the time of 1PPS-local signals arrive at TIC) and the stopping timepoint (the time of 1PPS-S signals arrive at TIC), the time interval between the two 1PPS signals can be determined as follows (Allan and Weiss,, 1980; Lee et al.,, 1990; Lewandowski and Thomas,, 1991):

where $T_{i}$ is time interval measurement of the two 1 PPS signals at station $i$; $t_{i}$ and $t_{s}$ are the emission timepoints of 1PPS signals of the atomic clocks from $i$-th station and satellite $S$, respectively; $\delta\rho_{is}$ is the error source (e.g. ionospheric delay, tropospheric delay, and Sagnac effect); $t_{r_{i}}$ is the receiver delay which is usually calibrated by the manufacturing company using simulated signals (Lewandowski and Thomas,, 1991); and $\rho_{is}$ is the geometric distance between $i$-th station and satellite $S$, expressed as:

By applying equation (7), the clock offsets between the $i$-th station and satellite S can be expressed as follows:

Finally, the clock offsets between two ground stations A and B can be determined by taking the common-view satellite’s clock as common reference:

We note that the cable delays from station’s atomic clock to TTR are measured beforehand with the value of 50.2 ns, as well as the cable delays from GNSS antenna to TTR with the value of 204.5 ns. And the cable delays are input into the TTR as a parameter. In fact, the cable delays are not significant in our study due to the fact that we compare the time elapse between two remote clocks.

Previous studies (Shen and Ding,, 2014) demonstrates that the ensemble empirical mode decomposition (EEMD) (Wu and Huang,, 2009) is an effective technique for isolating target signals from environmental noises. The EEMD technique is developed from the empirical mode decomposition (EMD) (Huang et al.,, 1998) (https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1998.0193). The signals series are decomposed into a series of Intrinsic Mode Functions (IMFs) and a residual trend $r$ after EEMD decomposition. These IMFs series that are sifted stage by stage reflect local characteristics of the signals, while the residual trend $r$ series reflects slow change of the signals (Zhu et al.,, 2013). The specific steps of EEMD are stated as follows (Wu and Huang,, 2009) (https://www.worldscientific.com/doi/abs/10.1142/S1793536909000047):

1) The overall signal series $X(t)$ is determined after adding Gaussian white noise $\omega(t)$ into original signal series $x(t)$:

2) The EMD method is applied to implement decomposition for $X(t)$ subsequently, these IMFs are determined as:

where $c_{j}$ denotes IMF_j and $r_{n}$ denotes the residual trend.

3) Adding different Gaussian white noise $\omega_{i}(t)$ into $x(t)$ and repeating step 1) and step 2), a series of $c_{ij}$ and $r_{in}$ are determined:

4) According to zero mean principle of Gaussian white noise, the interference cuased by Gaussian white noise can be cancelled after averaging.Therefore, the IMF series $c_{n}(t)$ can be determined as:

Finally, the original signal series $x(t)$ is decomposed into a series of IMFs series $c_{n}(t)$ and a residual trend series $r(t)$ by EEMD technique:

After EEMD decomposition of the signal series $x(t)$, the index of the orthogonality (IO) is calculated to check the completeness of the decomposition, which is defined as follows (Huang et al.,, 1998):

where $X^{2}(t)$ is defined as:

If the decomposition is completely orthogonal, the cross terms in the right-hand side of expression (29) should be zero, namely IO$=0$; and for the worst case IO$=1$.

Here, using a simulation experiments we explain the advantages of the EEMD technique. Suppose we have a synthetic series $S_{base}(t)$ that consists of 3 periodic signals series,$S_{i}(t)=A_{i}\cdot\sin(2\pi f_{i}t)\cdot\exp(-10^{-7}t)$(units:ns), where $A_{1}=10$, $A_{2}=A_{3}=5$, $f_{1}=2/86400$ Hz, $f_{2}=1/86400$ Hz, $f_{3}=0.5/86400$ Hz, respectively, and a linear signal series, $S_{4}(t)=3\cdot 10^{-4}t-10^{-7}$, with a data length of 10 days and sampling interval 960 seconds. The $S_{base}(t)$ can be expressed as:

The results of the constructed series are shown in Fig. 4. We added noise signals $N(t)$ into $S_{base}(t)$ with 2$\%$ (Case 1), 5$\%$ (Case 2) and 10$\%$ (Case 3) of the standard deviation (STD) of $S_{base}(t)$, and then construct a signal series $S(t)$ which contains noise. The $S(t)$ can be expressed as:

where the $N(t)$ is added noise, which consists of five types of noises, expressed as (Li and Sun,, 2011; Zhai et al.,, 2012; Ashby,, 2015):

where $N_{W-PM}(t)$ is the white noise phase modulation (W-PM), $N_{F-PM}(t)$ is the flicker noise phase modulation (F-PM), $N_{W-FM}(t)$ is the white noise frequency modulation (W-FM), $N_{F-FM}(t)$ is the flicker noise frequency modulation (F-FM), $N_{RW-FM}(t)$ is the random walk noise frequency modulation (RW-FM).

The constructed signal series $S(t)$ with different noise magnitudes are shown in Fig. 5. For our present purpose, the linear signal series $S_{4}(t)$ is the target signal which need to be extracted from the synthetic signals series $S(t)$.

The EEMD technique is applied to identify the periodic signals $S_{1}(t)$, $S_{2}(t)$ and $S_{3}(t)$ from $S(t)$. After EEMD decomposition, we use the index IO to check completeness of the decomposition. In our simulation experiments, the values of IO are 0.0015 (Case 1), 0.0017 (Case 2) and 0.0017 (Case 3), respectively, suggesting that signal series $S(t)$ are effectively decomposed in three different cases. The decomposed IMFs are shown in Figs. 6-8, and we found that the periodic signals $S_{1}(t)$, $S_{2}(t)$ and $S_{3}(t)$ can be clearly identified, respectively. The red dotted curves in Figs. 6-8 denote the original signals $S_{i}(t)$ ($i$ = 1, 2, 3) for comparison purpose.

Further, the Hilbert transform (HT) is executed to display the variety of instantaneous frequencies for each IMF component. The corresponding Hilbert spectra and marginal spectra for IMFs (expect for r) in three cases are shown in Figs. 9-11, where each of the subfigures (a)s shows the frequency variations of each IMF, and each of the subfigures (b)s shows measures of the total amplitude (or energy) contribution from each frequency value, representing the cumulated amplitude over the entire data span in a probabilistic sense (Huang et al.,, 1998).

From Figs. 9-11, we see that the mode-mixing problem exists in EEMD decomposition, and it becomes more obvious as noise increases from 2$\%$ to 10$\%$. However, the three periodic signals $S_{1}(t)$, $S_{2}(t)$, and $S_{3}(t)$ can be identified clearly. The set frequencies of $S_{1}(t)$, $S_{2}(t)$, and $S_{3}(t)$ are 2 circle per day (cpd), 1 cpd and 0.5 cpd, respectively, denoted as $f_{1-set}=2$ cpd, $f_{2-set}=2$ cpd, and $f_{3-set}=0.5$ cpd. After EEMD decomposition, the marginal spectra show that the corresponding values are $f_{1-case1}=2.148$ cpd, $f_{2-case1}=1.049$ cpd, $f_{3-case1}=0.599$ cpd (Case 1); $f_{1-case2}=2.248$ cpd, $f_{2-case2}=1.099$ cpd, $f_{3-case2}=0.599$ cpd (Case 2); $f_{1-case3}=2.347$ cpd, $f_{2-case3}=1.049$ cpd, $f_{3-case3}=0.599$ cpd (Case 3), respectively. And the detected signals corresponding to the original set signals $S_{1}(t)$, $S_{2}(t)$ and $S_{3}(t)$ are shown by the peaks denoted in green circles. In addition, there appear not-real signals with frequencies around 0.2 cpd after EEMD decomposition, which are denoted as $F(t)$ and shown by the peaks denoted in red rectangle. The $F(t)$ might be meaningless in physical explaining, and it will be an interference if we focus on periodic signals. In this study, however, the target is to detect and identify the linear signal $S_{4}(t)$, which means that all periodic signals are useless and should be removed. After all these periodic signals are removed, we reconstruct a new signal series by summing the residual IMFs and the residual trend r. The reconstructed signal series is denoted as $S^{{}^{\prime}}(t)$.

We take the signal series $S_{4}(t)$ as real signal series, and try to recovery it by two different methods. The first method is to perform a least squares linear fitting on the signal series $S(t)$ directly; and the second method is to perform EEMD decomposition on $S(t)$ and then reconstruct the new signal series $S^{{}^{\prime}}(t)$ by removing periodic series IMFs, and the least squares linear fitting is performed on $S^{{}^{\prime}}(t)$. The results are shown in Fig. 12.

Here we use the standard deviation (STD) of the difference series between the real signal series $S_{4}(t)$ and that determined by Method 1 or 2 to evaluate the reliability of the two methods. Using Method 1, namely the direct least squared linear fitting of the series $S(t)$, the STDs of the results are 1.31 ns (Case 1), 1.63 ns (Case 2), and 2.16 ns (Case 3), respectively. Using Method 2, namely the least squared linear fitting of the reconstructed series $S^{{}^{\prime}}(t)$, which is obtained after removing the periodic series via EEMD technique, the STDs of the results are 0.82 ns (Case 1), 1.16 ns (Case 2), and 1.87 ns (Case 3), respectively. The comparative results clearly suggest that the EEMD technique is effective for extracting the linear signals of interest by a priori removing the contaminated periodic signals from the original observations.

In our experiments, we used two hydrogen time and frequency standards, with one is a fixed reference clock $C_{A}$ (iMaser3000) and the other is a portable clock $C_{B}$ (BM2101-02). The experiments are conducted at Beijing 203 Institute Laboratory (BIL) and Luojiashan Time-Frequency Station (LTS), and the locations of the two stations are shown in Fig. 13 and Table 2. And we used the Modified Allan deviation (MDEV)(Allan et al.,, 1991) to evaluate the frequency stability of the two hydrogen atomic clocks, which is estimated from a set of frequency measurements for time averaging. The MDEV is expressed as follows (Allan et al.,, 1991; Bregni and Tavella,, 1997; Lesage and Ayi,, 2007):

where $N$ denotes the total number of the time series, $\tau_{0}$ denotes the basic time interval, and $\tau=n\tau_{0}$.

The frequency stabilities of the two clocks on the different average time intervals are given in Table 3. The MDEV with shorter interval represents random noise of the clock, while with longer interval represents secular frequency drift (Kopeikin et al.,, 2016).

The experiments are divided into two Periods (shown in Fig. 14). In Period 1 that spans from January 09, 2018 to January 29, 2018, the zero-baseline measurement is implenented at BIL with clocks $C_{A}$ and $C_{B}$ are placed in the same shielding room with the same geopotentials. According to the International GNSS Tracking Schedule (Allan and Thomas,, 1994; Defraigne and Petit,, 2015), the clock comparisons between $C_{A}$ and $C_{B}$, $\Delta C_{AB}$, are conducted via CVSTT technique, and a series of clock offsets $\Delta C_{AB}(t)$ can be obtained, where $\Delta C_{AB}=C_{B}-C_{A}$, and $\Delta C_{AB}(t)=C_{B}(t)-C_{A}(t)$.

When zero-baseline measurement is finished, the portable clock $C_{B}$ is transported to LTS from BIL while keeps clock $C_{A}$ unmoved at the BIL, and during transportation the portable clock $C_{B}$ continues operating. After installation and adjustment of $C_{B}$ at the LTS, the clock comparisons between $C_{A}$ and $C_{B}$ are conducted in geopotential difference measurement, similarly. And this period is denoted as Period 2, which spans from Febrary 01, 2018 to April 07, 2018. It should be noted that there is no difference except for the location diversity of clock $C_{B}$, when comparing the setup difference of Period1 and Period 2. In addition, there are two days of missing data (from January 30, 2018 to January 31, 2018) due to the fact that during transportation (around 8 h) there are no observations and the observations from the day after the re-observation are not stable and consequently not used.

Due to the fact that temperature is a major factor which may disturb the performance of atomic clocks as well as the cable delays (Weiss et al.,, 1999; Lisowiec and Nafalski,, 2004),we controll the environment with a relatively constant temperature in whole experiments. The temperature is held nearly constant ( $24\pm 0.5^{\circ}C$) for the fixed clock $C_{A}$ in Period 1 and Period 2 at BIL. For the portable clock $C_{B}$ , it is in the same laboratory with clock $C_{A}$ during Period 1, and the temperature is locked at $24\pm 1.5^{\circ}C$ in Period 2 at LTS.

Based on the observations via CVSTT technique, the clock offsets series $\Delta C_{AB}(t)$ in each Period are determined. Here, the clock offsets series $\Delta C_{AB}(t)$ of the observations for all common-view satellites during Period 1 and Period 2 are denoted by the black dotted curve (simply denoted as raw data) in Fig. 15. And take the satellite elevation angles into consideration, the initial CVSTT observations (initial observations) are obtained, which are denoted by the yellow curve in Fig. 15. After that, we performe data preprocessing on the initial observations by removing gross error, correcting clock jumps and inserting missing data, and the preprocessed CVSTT data (preprocessed data) are finally obtained which denoted by the blue curve in Fig. 15. After data preprocessing, the clock offsets series $\Delta C_{AB}(t)$ become more stable. In the following data processing, all operations are based on the preprocessed data sets.

At the same time, it should be noted that, however, the observations from MJD 58137 to MJD 57141 are not used due to a systematic malfunction, and the observations from MJD 58142 to MJD 57147 are also not used due to random movement of the GNSS antenna (the two clocks are not transported, but one antenna is replaced). Therefore, the valid observations in Period 1 spans from MJD 58127 to MJD 57136 (seen in Fig. 15(a)). In Period 2, since the stability of clock $C_{B}$ could be significantly influenced by various environment factors during its transportation to LTS from BIL, the initial observations at LTS contains large fluctuations and noises (from MJD 58150 to 58161, seen in Fig. 15(b) and Fig. 20(b)), which will contaminate final determination of the geopotentials. Hence, as for preprocessed data in Period 2, we adopt two different preprocessed data sets by whether contains the observations of the first 12 days, and the details are shown later.

Then, the EEMD technique is applied to the preprocessed data sets for removing the uninteresting periodic components and extracting the geopotential-related signals. The preprocessed data sets in Period 1 and Period 2 are decomposed into a series of IMFs (with frequencies from higher to lower) and a long trend component $r$ after EEMD decomposition, as shown in Figs. 16 and 17. Generally, the components of the EEMD convey a physical meaning due to the fact that the characteristic scales are physical (Huang et al.,, 1998). However, the first several high frequency components might be fictitious due to the fact that the sampling interval (960s) is too large to capture the high-frequency variations. As a result, the data are jagged at these frequencies (e.g. $c_{1}$ in Fig. 16; $c_{1}$, $c_{2}$, $c_{3}$ in Fig. 17).

After using EEMD technique for preprocessed data sets in Period 1 and Period 2, we examine completeness and orthogonality of the EEMD decomposition. The IO values corresponding to Period 1 and Period 2 are determined with values of $-0.0085$ and $-0.0094$, respectively, demonstrating that the signals are effectively decomposed. It should be noted that the IO value during Period 2 is worse than that in Period 1. This might be due to the environmental difference of $C_{B}$ and $C_{A}$ during Period 2, or due to the difference of the signal propagation paths during Period 2.

The corresponding Hilbert spectra in Period 1 and Period 2 are given in Fig. 18. In order to measure the energy contributions from each frequency value, the marginal spectra of all IMFs are given in Fig. 19, which represent the cumulated amplitude over the entire data span in a probabilistic sense. From Figs. 18 and 19, the periodic signals with frequencies of around 0.2 cpd, 0.5 cpd, 1 cpd and 3.3 cpd of zero-baseline measurement are detected; and the periodic signals with frequencies of around 0.05 cpd, 0.36 cpd, 1 cpd and 2 cpd of geopotential difference measurement are detected. Finally, these periodic signals included in the original CVSTT observation series are removed, and we reconstructed the residual clock offsets series $\Delta C_{AB-re}(t)$ by summing the residual components. The reconstructed clock offsets series $\Delta C_{AB-re}(t)$ are regarded as the geopotential-related signals, based on which we can determine the geopotential difference between $C_{B}$ (LTS) and $C_{A}$ (BIL). The reconstructed clock offsets series $\Delta C_{AB-re}(t)$ are denoted as EEMD results, as described in section 4.2.