*Article* **Measuring Height Difference Using Two-Way Satellite Time and Frequency Transfer**

**Peng Cheng 1, Wenbin Shen 1,2,\*, Xiao Sun 1, Chenghui Cai 1, Kuangchao Wu <sup>1</sup> and Ziyu Shen <sup>3</sup>**


**Abstract:** According to general relativity theory (GRT), the clock at a position with lower geopotential ticks slower than an identical one at a position with higher geopotential. Here, we provide a geopotential comparison using a non-transportable hydrogen clock and a transportable hydrogen clock for altitude transmission based on the two-way satellite time and frequency transfer (TWSTFT) technique. First, we set one hydrogen clock on the fifth floor and another hydrogen clock on the ground floor, with their height difference of 22.8 m measured by tape, and compared the time difference between these two clocks by TWSTFT for 13 days. Then, we set both clocks on the ground floor and compared the time difference between the two clocks for seven days for zero-baseline calibration (synchronization). Based on the measured time difference between the two clocks at different floors, we obtained the height difference 28.0 ± 5.4 m, which coincides well with the tapemeasured result. This experiment provides a method of height propagation using precise clocks based on the TWSTFT technique.

**Keywords:** TWSTFT; satellite; geopotential; altitude transmission

## **1. Introduction**

The gravity potential (geopotential) plays a significant role in geodesy. It is essential for defining the geoid and measuring orthometric height. The conventional method of determining the geopotential is combining leveling and gravimetry, but there are shortcomings: with the increase in measurement length, the error accumulates and becomes larger and larger, and it is impossible or difficult to transfer the orthometric height between two points separated by oceans.

To overcome the shortcomings existing in the conventional method, time–frequency comparison methods based on general relativity theory (GRT) [1] were proposed in recent decades [2–7]. The basic idea is that by comparing the time elapsed or frequency shift between two remote clocks, the geopotential difference between the two sites where the clocks are located could be determined.

In the clock-transportation method, the most critical conditions are the clock's precision and time transfer. Precise clocks generally include microwave-atomic clocks (MACs) and optical-atomic clocks (OACs). MACs play an important role in time service research. However, with the development of a high-precision clock, its precision could not match that of the OAC. The concept of OAC was first proposed by Nobel laureate Dehmelt (1973) in the late 1970s. He used the energy level transition of a single ion to realize an ultra-highprecision optical clock. In recent years, the precision of OAC has reached a total systematic uncertainty of 9.4 × <sup>10</sup>−<sup>19</sup> and frequency stability of 1.2 × <sup>10</sup>−<sup>15</sup> / <sup>√</sup>*<sup>τ</sup>* [8]. However, OACs

**Citation:** Cheng, P.; Shen, W.; Sun, X.; Cai, C.; Wu, K.; Shen, Z. Measuring Height Difference Using Two-Way Satellite Time and Frequency Transfer. *Remote Sens.* **2022**, *14*, 451. https://doi.org/10.3390/rs14030451

Academic Editor: Xiaogong Hu

Received: 15 December 2021 Accepted: 16 January 2022 Published: 18 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

are often bulky and can only work in laboratory environments, which significantly limits their application scope and makes it difficult to conduct a clock-transportation experiment. It has always been the wish of scientists to realize a transportable, reliable, and quasicontinuous high-precision OAC, but it is also challenging. To broaden the application scope of OAC, many research groups in the worldwide have been devoted to the development of transportable optical-atomic clocks (TOCs). In 2014, a group reported a TOC based on laser-cooled strontium atoms trapped in an optical lattice, and this TOC fits within a volume of <2 m3, and its relative uncertainty is 7 × <sup>10</sup>−<sup>15</sup> [9]. Three years later, a group from Physikalisch-Technische Bundesanstalt (PTB) reported a TOC with 87Sr, with its characterization against an OAC resulting in a systematic uncertainty of 7.4 × <sup>10</sup>−<sup>17</sup> [10]. With the development of TOC, many scientists started to conduct clock-transportation experiments. Grotti et al. (2018) reported the first field measurement campaign with a 87Sr TOC with an uncertainty of 1.8 × <sup>10</sup>−<sup>16</sup> and a 171Yb OAC with an uncertainty of 1.6 × <sup>10</sup><sup>−</sup>16. They used these clocks and fiber link to determine the geopotential difference between the middle of a mountain and a location 90 km away with a height difference of 1000 m. Their experimental result of potential difference of 10,034(174) m2s−<sup>2</sup> agrees well with value of 10,032.1(16) m2s−<sup>2</sup> determined independently by the conventional geodetic approach [3].

The clock-transportation experiments mentioned above used optical fiber to transfer frequency signals. Although optical fiber has very high accuracy, distance still limits its application. By comparison, though GNSS common-view technique and TWSTFT have lower accuracy than optical fiber frequency signal transfer [11–15], they can realize longdistance time–frequency signal transmission. In TWSTFT, one uses a geostationary satellite as 'bridge' to transmit time–frequency signals from one station to another one, and the time elapse recorded at one station is compared with that at another one. Because using TWSTFT only requires the stations in a place where the geostationary satellite can receive and transmit signals, there is hardly any limit on the positions of stations. Most errors in the signal-transmission process are offset because of the approximate symmetry of the signal transmission path of TWSTFT technology. Therefore, this symmetry causes the high precision of this technology [16]. In August 1962, the USNO (U.S. Naval Observatory) launched the first communication satellite, Telstar I, to transmit telephone and high-speed data signals. Then, the USNO and the NPL (National Physical Laboratory) collaborated in an experiment using this satellite to relate the precise clocks at the USNO and the RGO (Royal Greenwich Observatory). This is considered to be the first two-way satellite time transfer experiment, and the accuracy of the experiment was 20 μs [17]. In 1992, some satellite systems and modems adapted for TWSTFT were commercialized. About ten coordinated universal time (UTC) laboratories have equipped with the modems and other equipment for the clock comparisons with TWSTFT [18]. Later, many TWSTFT experiments were conducted, and they almost exchanged timing information via the communication satellite; paired ground stations transmit and receive pseudo-random noise (PN) coded signals in TWSTFT links. The TWSTFT technique became promising using geostationary satellites for high-accuracy time and frequency transfer [16,19–21]. The accuracy of TWSTFT has been further improved to around 0.2 ns, and its improvement has been seriously limited by the chip rate of the coded signal [22]. Therefore, further improvements should come from the use of carrier phase information, because the resolution of the carrier phase is 100 to 1000 times more accurate than that of the code [23]. In 2016, an experiment using carriedphase TWSTFT was performed between the two stations of NICT (National Institute of Information and Communications Technology) and KRISS (Korea Research Institute of Standards and Science) with Sr and Yb OAC, and the instability for a frequency transfer at the 10−<sup>16</sup> level after 12 h was achieved [24]. Riedel et al. (2020) conducted a 26-day comparison of five simultaneously operated OACs and six MACs located at SYRTE, NPL, INRIM, LNE, and PTB by using TWSTFT and GPSPPP. Considering the correlations and gaps of measurement data, they improved the statistical analysis procedure; combined overall uncertainties in the range of 1.8 × <sup>10</sup>−<sup>16</sup> to 3.5 × <sup>10</sup>−<sup>16</sup> for the OAC comparisons

were found [25]. To investigate the feasibility of transportable atomic clock comparison using TWSTFT, we conducted a MAC comparison experiment at the Beijing Institute of Radio Metrology and Measurement (BIRMM), Beijing [26,27].

In Section 2, we introduce the principle of measuring the height difference by the TW-STFT, and in Section 3, we discuss the error corrections of TWSTFT. Section 4 demonstrates our experiment and data processing. In Section 5, we provide the results. Conclusions and discussions are placed in Section 6.

#### **2. Methods**

#### *2.1. Height Measurement Based on Time Difference*

According to general relativity theory (GRT), a clock at a position with lower height (stronger geopotential) ticks slower than an identical at a position with higher height (weaker geopotential) [5,28]. Inversely, we can determine the geopotential difference between two points A and B by measuring the elapsed time difference between two clocks located at A and B, expressed as (accurate to 1/c2) [7,29–31]

$$\frac{\Delta t\_{AB}}{T} = \frac{t\_B - t\_A}{T} = -\frac{W\_B - W\_A}{c^2} = -\frac{\Delta W\_{AB}}{c^2} \tag{1}$$

where *t*<sup>A</sup> and *t*<sup>B</sup> denote the times at sites *A* and *B*, respectively, after a standard time period of T; *W*<sup>A</sup> and *W*<sup>B</sup> are the geopotential at sites *A* and *B*, respectively (note that we apply the definition of geopotential given by geodetic community); and *c* is the speed of light in the vacuum. From Equation (1), we can determine the geopotential difference Δ*WAB* = *WB* − *WA* based on Δ*tAB*/*T*.

As shown in Figure 1, given the height of point A and the geopotential difference Δ*WAB* between *A* and *B*, one can measure the orthometric height of point *B*, expressed as [32,33]

$$H\_B = H\_A - \frac{\Delta W\_{AB}}{\overline{\mathcal{g}}} = H\_A + \frac{\Delta t\_{AB}}{T} \cdot \frac{c^2}{\overline{\mathcal{g}}} \tag{2}$$

where *HA* and *HB* are the orthometric heights of point A and B, respectively, and *g* is the 'mean value' between *g*<sup>B</sup> at point B and *g*OB at point OB

**Figure 1.** Principle of determining the orthometric heights. WA and WB are the geopotentials at point A and B, OA and OB are the projection points on the geoid (bold dashed curve) corresponding to point A and B along the plumb lines (light-blue dashed curves), red dashed curve denotes equipotential surface passing through point A, W0 is the geopotential on the geoid.

#### *2.2. One Pulse Per Second (1 PPS) Signal*

The 1 PPS signal provides precise clock synchronization [34]. The original signal is the frequency signal (with typical frequency10 MHz) output by the atomic clock. As an analog signal, in distant transmission, the information carried by the original frequency signal is easily distorted by various interferences [35]. Therefore, it needs to be converted into a digital signal, 1 PPS, to complete the time-information transmission. The process of converting frequency signal (say 10 MHz) into 1 PPS signal is shown in Figure 2. Suppose the frequency of the reference signal is 10 MHz, 107 cycles of the reference signal are one second. The generated 1 PPS signal rises from a low electrical level at the beginning of one reference signal cycle, and keeps the high electrical level for a short time (generally the pulse width is 20 μs), then declines to a low electrical level and keeps the site until the end of the 107 cycles (calculated from the first rise). This is a complete cycle (1 s) of 1 PPS signal. When the 1 PPS signal is used for time signal comparison, the rising edge of the signal will be used as the point to trigger the timer. Usually, the electrical level of the trigger is a predetermined value between the lowest electrical level and the highest electrical level. Ideally, when the electrical level rises to this predetermined value, the switch is triggered immediately, and the timing starts. However, there is a trigger delay during the process of the trigger switch. In fact, the switch can only be triggered when the actual electrical level is slightly higher than the predetermined value, so the trigger time will be within the time corresponding to the red oblique line [36]. Hence, the rising time duration *δt* is very important for precise time synchronization. Usually, *δt* should be smaller than 10 ns, since if *δt* is too large, the rising edge's slope (the blue slope line at the bottom of Figure 2) will become too small, which will increase the uncontrollable duration (red line) and reduce the time synchronization precision.

**Figure 2.** Generation of the one pulse per second (1 PPS) signal.

### *2.3. Transmission of 1 PPS Signal*

Figure 3 shows the transmission process of 1 PPS signal in the TWSTFT. When the 1 PPS signal is generated by a clock, a part of the 1 PPS signal enters time interval counter (TIC) as trigger gate open pulse. Another part of the 1 PPS signal is transmitted to a satellite, through the modulation and emitter. The satellite uses a transparent transponder to transmit it to another ground station. After the signal from the satellite is received, the receiver and demodulation recover the 1 PPS signal from the received signal. The recovered 1 PPS signal enters the TIC as a pulse to trigger the gate close. The above processes are conducted between each of the two stations.

**Figure 3.** Principle of the two-way satellite time and frequency transfer (TWSTFT). There are two same TWSTFT observation systems located at two sites A and B. Every system include clock, emitter, receiver, etc. (see text) (modified after ITU-R 2015).

If the 1 PPS signal is directly transmitted to another station, it is difficult to compare the 1 PPS signals generated by two clocks at two sites synchronously. Hence, using the PN code to accompany the 1 PPS signal for time marking is necessary. The maximal length linear feedback shift register sequences (m-sequences) are the PN code usually used in TWSTFT. The m-sequence is the largest code that can be generated by a given shift register or a delay element of a given length [37]. The code has good autocorrelation characteristics and can be easily generated (Enge et al., 1987). At the same time, to reduce the impact of the bit error on the final time comparison, the frame synchronization bits (FSB) will be inserted at the rising edge of the 1 PPS signal. The bit error means that in signal transmission, decay causes the signal to be damaged, so the originally transmitted signal is '1', but the received signal is '0' (it should also be '1'), and the FSB is a specific set of bits at the beginning of each frame of signal [38]. The spreading spectrum consists of 1 PPS and PN, which is completed when the 1 PPS signal combines with the synchronized PN code by way of modulo-2 sum in the code generator (CG). This spread spectrum signal (SSS) is the baseband signal with wide bandwidth and decays rapidly during transmission. If it is directly used for long-distance baseband transmission, due to quick attenuation, the signal at the receiver will have a too-low signal-to-noise ratio (SNR) to be identified [39]. Therefore, the SSS needs to be modulated on a proper carrier. Generally, the carrier is a cosine wave, with its amplitude, frequency, and phase being known, and is used to transmit the information (signal) by changing the carrier's amplitude, or frequency, or phase, or the combinations of these entities.

The SSS is modulated on the carrier in the modulator with the binary phase-shift keying (BPSK) [40]. The BPSK changes the initial phase of the carrier to transmit the binary digital information of the SSS, while the carrier's amplitude and frequency remain unchanged (in the modulation process). In this modulation technology, the initial phase of the carrier has only two values of '0' and 'π', which correspond to the digits '1' and '0' of the 1 PPS signal, respectively. When transmitting signal '1', the modem generates a carrier with initial phase 0, and when transmitting signal '0', the modem generates a carrier with the initial phase π. The modulated signal is an intermediate frequency (IF) signal, which is named 1 PPSTX.

Referring to Figure 3, the frequency of the 1 PPSTX signal generated at station A is about 70 MHz, which is not suitable for long-distance signal transmission. Therefore, the frequency of the 1 PPSTX signal is amplified to 14 GHz through the multiplier in the up-converter (UC) and transmitted to the geostationary satellite through the transmitter.

After it receives the up-link signal, the satellite transparent transponder (STT) converts the signal's frequency to 12 GHz and transmits the signal (as a down-link signal) to station B. When the down-link signal arrives at station B, the signal's frequency is reduced to about 70 MHz through the multiplier in the down-converter (DC) (the function of which is similar to UC, but in the opposite way). At the same time, the received signal will also be digitized to facilitate subsequent processing. Usually, when this signal is converted to the IF range, a frequency shift will impact the signal acquisition. The reason is that the carrier used to demodulate the IF signal to the baseband signal needs to have the same frequency as the IF signal. When capturing the received signal, the synchronization information of the signal is searched, and the known PN code is used to synchronize the received signal and PN code. Then, the local carrier frequency and code phase are almost identical to the received signal. After the capture is implemented, the signal will be tracked in the delay-locked loop (DDL), which can fully synchronize the local carrier and the PN code with the received signal, thus removing the carrier and PN code from the received signal and recovering the 1 PPS signal. Thus far, we have completed the signal demodulating and despreading and obtained the transmitting 1 PPS signal generated by station A [41]. However, due to the existence of bit error, it is necessary to detect the synchronization code in the 1 PPS signal to ensure that the correct synchronization point of the 1 PPS signal is generated [38]. Since it has been synchronized by PN code, the recovered 1 PPS signal needs to be input to TIC to compare with the local 1 PPS signal. Each rising edge of the local 1 PPS signal will be used as the opening point of the timing gate, and each corresponding rising edge of the recovered 1 PPS signal will be used as the closing point of the timing gate. In the TIC, the 10 MHz signal output by the local clock is used as the time base reference, measuring the time length between the opening and the closing each time, and then the TIC outputs the time difference as the observed value. Therefore, at both clock sites, the time signals are transmitted nominally at the same instant. Each clock site receives the signal from the other clock site, and its arrival time is measured. After exchanging the measured data, the time–frequency difference between the two clocks is calculated.

#### *2.4. Time Difference Calculation in TWSTFT*

The TWSTFT is based on the exchange of timing signals through geostationary telecommunication satellites. Any one of the clocks at site A and B generates a pulse signal every second (1 PPS signal). If the two clocks run at the same rate, the time difference between the two pulses holds the same at every second. Therefore, measuring the time difference between the two pulses, the running rate difference between the two clocks can be calculated. The task of TWSTFT technology is how to accurately compare the two time pulses from two stations.

It is assumed that the clocks of the two stations A and B have been synchronized in advance. At the appointed time, the clock *CA* and *CB* at A and B simultaneously generate the pulse signals *PA* and *PB*. A part of *PA* enters TIC(A) as trigger gate open pulse. Another part of *PA* is transmitted to station B. After a short time, the pulse *PB* is coming and entering the TIC(A) as a pulse to trigger the gate close. The TIC(A) will calculate and record the time from gate opening to gate closing. The above process is repeated every second.

The TIC reading at station A is expressed as:

$$
\pi\_I^A = \pi^A - \pi^B + \pi\_t^B + \pi\_{pu}^B + \pi\_{su}^B + \pi\_s^B + \pi\_{pd}^A + \pi\_{sd}^A + \pi\_r^A \tag{3}
$$

and that at station B is expressed as:

$$
\tau\_I^B = \tau^B - \tau^A + \tau\_t^A + \tau\_{pu}^A + \tau\_{su}^A + \tau\_s^A + \tau\_{pd}^B + \tau\_{sd}^B + \tau\_r^B \tag{4}
$$

Combining Expressions (3) and (4), the time scale difference could be expressed as:

$$\begin{aligned} \tau^A - \tau^B &= \frac{1}{2} [ \left( \tau\_I^A - \tau\_I^B \right) + \left( \tau\_s^A - \tau\_s^B \right) - \left( \tau\_{sd}^A - \tau\_{su}^A \right) + \left( \tau\_{sd}^B - \tau\_{su}^B \right) + \left( \tau\_{pu}^A - \tau\_{pd}^A \right) \\ &- \left( \tau\_{pu}^B - \tau\_{pd}^B \right) + \left( \tau\_t^A - \tau\_r^A \right) - \left( \tau\_t^B - \tau\_r^B \right) \end{aligned} \tag{5}$$

where *τ<sup>k</sup>* denotes the local time-scale, and k means station k (k = A, B); *τ<sup>k</sup> <sup>I</sup>* is the TIC reading; *τk su* and *τ<sup>k</sup> sd* are the Sagnac effect delay in the up-link and down-link, respectively; *<sup>τ</sup><sup>k</sup> pu* and *τk pd* the signal path up-link and down-link delays, respectively; *<sup>τ</sup><sup>k</sup> <sup>s</sup>* is the satellite path delay through the transponder, *τ<sup>k</sup> <sup>t</sup>* the transmitter delay, and *τ<sup>k</sup> <sup>r</sup>* the receiver delay [42].

#### **3. Error Analysis and Corrections**

The error sources in TWSTFT observations mainly come from three aspects: (1) equipment delay errors; (2) signal-propagation-path-delay errors; and (3) Sagnac effect errors.

#### *3.1. Equipment Delay Error*

Equipment delay error mainly includes TIC measurement error, modem errors, satellite transparent transponder delay error, and the emitting and receiving delay error of Earth's station equipment. The TIC measurement error and modems error, which are caused by their own measurement accuracy errors, are about a dozens of picoseconds and 100 ps [42], respectively. The satellite transparent transponder delay error, mainly due to the signal from A to B and the signal from B to A through different forwarding channels of the transparent transponder, is hard to control, and is generally in the range of 100 ps [43]. The emitting and receiving equipment delay error of the ground station is around 0.2~0.5 ns, including cable delay, the transmission and receiving system error, and temperature variation.

It may not be accurate enough to measure the delay error of the signal after passing through each piece of equipment alone. Therefore, the zero-baseline measurement could be a better way to measure all equipment-delay errors. Zero-baseline measurement is when two atomic clocks are placed at the same place simultaneously, and then a TWSTFT experiment is conducted.

#### *3.2. Propagation Path Delay Error*

The signal's propagation path delay error is caused by the satellite's signal path delay errors of up-link and down-link, mainly including the delay propagation path geometry error and tropospheric delay error, ionospheric delay error, and delay errors of the different distances between two stations and the satellite.

The propagation path geometry delay is related to the coordinates of satellites and ground stations. For the signal's arriving and receiving time, if there is an error for ground stations or satellite position, the error will directly affect the time delay between the satellite and ground station. However, the clock difference calculation formula includes the difference between the up-link *τ<sup>k</sup> pu* and down-link path *τ<sup>k</sup> pd*; therefore, through the difference between two paths, one can eliminate the error with only a few picoseconds left. That means, at the same time, coordinate errors of satellite and ground stations do not affect the calculation of the clock difference.

Tropospheric delay is also called the tropospheric refraction error. The troposphere, through changing the propagation path of the signal, causes time delay. Many factors can affect tropospheric delay, such as ground climate, atmospheric pressure, temperature, and humidity TEC. Tropospheric delay can use the tropospheric delay model to correct; common models are the Hopfield model, Saastamoinen model, EGNOS model, etc. [44–46]. The differences between each model are mainly in the low-altitude angle; at the zenith direction, the difference is very small. Moreover, the up-link and down-link paths are symmetrical, hence, the tropospheric delay can be mostly cancelled. The remaining time delay is less than 10ps, and it can be neglected in the present study [42].

The ionospheric delay error is mainly caused by charged particles in the ionosphere. The charged particles can change the speed and path of the signal. The ionospheric delay depends on the total electron content (TEC) and the signal's frequency. Because the up-link and down-link signals at each station differ in carrier frequency, the following formula is used to correct the ionospheric delay error [47]:

$$\left. \tau\_{pu}^{k} \right|\_{\rm ion} - \left. \tau\_{pd}^{k} \right|\_{\rm ion} = \frac{40.28 \,\epsilon\_{ks}}{c} \left( \frac{1}{f\_{\rm U}^{2}} - \frac{1}{f\_{\rm D}^{2}} \right) \tag{6}$$

where *fu* and *fD* are the up-link and down-link frequencies of the signals, respectively; *eks* is the TEC along the signal-propagation path between ground station *k* and satellite S, and *c* is the speed of light in vacuum. Thus, we have the following equation to correct the ionospheric delay:

$$\frac{1}{2} \left[ \left. \pi\_{pu}^{A} \right|\_{ion} - \left. \pi\_{pd}^{A} \right|\_{ion} \right] - \frac{1}{2} \left[ \left. \pi\_{pu}^{B} \right|\_{ion} - \left. \pi\_{pd}^{B} \right|\_{ion} \right] = \frac{20.14 (e\_{As} - e\_{Bs})}{c} \left( \frac{1}{f\_{\mathcal{U}}^{2}} - \frac{1}{f\_{\mathcal{D}}^{2}} \right) \tag{7}$$

Obviously, after the two-way signal transfer, the total effect of ionospheric delay could be further reduced. However, to better eliminate the error, we need to know the Δe= *eAs* − *eBs*. Usually, the global ionospheric TEC map provided by the IGS (International GNSS Service) every two hours can be used to eliminate or reduce the delay error. Taking the typical value of TEC, it can be estimated that the ionospheric delay is about 0.1 ns.

There is an error caused by the distances between the two stations and the satellite. The distance between the two stations and satellite is different, which means the signal arrival time from A to S and B to S is different. There is the problem that, when the signal of one station arrives at the satellite, another signal is on the way. Because the geostationary satellite is not strictly relative static with the Earth, the satellite's location is different between the first-arriving and the late-arriving. It will cause the signal path from A to B and B to A to be asymmetrical, seriously influencing the advantages of symmetrical paths of TWSTFT in eliminating errors. Approximately, when the distance converted from the longitude difference between the two ground stations and the satellite longitude is 300 km, the maximum error is about 30 ps. To reduce signal arrival time difference, we can slightly adjust the transmission delay to make them arrive at the satellite simultaneously. If the time difference between the two signals arriving at the satellite is less than 5 ms, the effect in TWSTFT will be less than 1 ps [43].

#### *3.3. Sagnac Effect Error*

When we transmitted the signal from the station to the satellite, the signal's route is not time-varying [45]. Since the satellite and ground stations are moving, this motion causes the signal-propagation change; the corresponding influence is referred to as the Sagnac effect [48]. After the two-way link difference, the error caused by the Sagnac effect is about 100-200 ns and needs further correction. The Sagnac effect correction for a one-way route from the ground station to satellite is given in a model that provides sufficient accuracy by the following expression [42,48,49]:

$$\tau\_{sd}^k = \frac{\Omega}{c^2} \mathcal{R} \left\{ a \cos \left[ \tan^{-1} \left( \tan \left. \varphi^k - f \tan \left. \varphi^k \right) \right| \right] + H^k \cos \varphi^k \right\} \sin \left( \lambda^k - \lambda^k \right) \right. \tag{8}$$

where Ω is the Earth's rotation rate, R is the distance from the satellite to the geocenter, *a* is Earth's equatorial average radius, *f* is the flattening of the Earth ellipsoid, *H<sup>k</sup>* (k = A or B) is the height of the station above the ellipsoid, *λ<sup>s</sup>* is the longitude of the satellite, and *ϕ<sup>k</sup>* and *λ<sup>k</sup>* are the latitude and longitude of the station, respectively.

After the model correction, only about a 10–100 picosecond error will be left. The reason for this residual is that, for the ground observer, the position of a geostationary satellite is not completely fixed. There will be a slight periodic movement with a daily period around a central point. As shown in Figure 4, in terrestrial reference frame, the position of satellite relative to station changes over time. That means the *λ<sup>s</sup>* and R in Equation (8) will change over time. This leads to Sagnac effect error, which is no longer a constant but varies with the orbital period of the satellite in a Earth-fixed system, and its maximum peak to peak amplitude is hundreds of ps. In our study, the influence magnitude is 10–100 ps. At the current accuracy level, it has almost no effect on our experimental results, but it may be considered if higher accuracy is needed. Various error sources and their magnitudes of influence are listed in Table 1.

**Figure 4.** The influence of small periodic movement of satellite on TWSTFT. Two antennas *TA* and *TB* located, respectively, at station A and B, are connected with clocks. S is the geostationary satellite.


### **4. Experiments and Data Processing**

### *4.1. Experiments*

To verify the feasibility of this method, we need clocks with high stability to maintain a stable frequency standard and a reliable time-transfer technique to measure the time difference between two positions. Therefore, we used two hydrogen atomic clocks in the experiment; a non-transportable clock *CA* (H-MASER VCH-1003A) and a transportable clock *CB* (H-MASER BM2101-01); both of their relative nominal frequency stabilities are <sup>5</sup> × <sup>10</sup>−<sup>15</sup> in one day. The TWSTFT was used as the time-transfer technique, one of the most accurate time-transfer methods than other GNSS-related techniques [50]. Time transfer by satellite does not have higher stability than fiber link, but it has better performance for long-distance time transfer.

We conducted the experiment in a building at the BIRMM, Beijing. The experiment we conducted was a geopotential difference measurement from 15 December 2016 to 27 December 2016 for a total of 13 days. During this process, *CB* was moved to the fifth floor, while *CA* was placed on the ground floor (see Figure 5a). The height difference between the fifth floor and ground floor is 22.8 m. After the experiment was conducted, we carried out a zero-baseline measurement from 27 December 2016 to 3 January 2017 for a total of seven days to calculate the clock drift and equipment error for calibration. As Figure 5b shows, clocks *CA* and *CB* were both placed on the ground floor. The equipment was temperature-stabilized and controlled during the experiment. At every site (the fifth floor or ground floor), through connection cables, one hydrogen clock, integrated control cabinet (ICC), and a satellite antenna were connected (see Figure 5).

**Figure 5.** Time difference measurement. (**a**) Geopotential difference measurement. (**b**) Zero-baseline measurement. *TA* and *TB* are the antennas connected to clock *CA* and clock *CB*, respectively. ICC is integrated control cabinet, which includes modulation, upconvertor, downconvertor, demodulation, time interval counter, etc. S is the geostationary satellite used in the experiments.

#### *4.2. Data Processing*

In the experiments, the data we collected were the time difference between clock *CA* and *CB* corresponding to time. The sampling rate of data is 1 Hz. According to Equation (1), we just want to calculate <sup>Δ</sup>*tAB <sup>T</sup>* , which is the slope of data, so it is insignificant whether the two clocks are synchronized in the beginning. What we are concerned with is the running rates of the two clocks.

There are some defects because of an equipment sampling error in the row data (Figure 6a,b). There are accidental data jumps, outliers, and some missing data. To improve the quality of the data, the following procedures are adopted:


**Figure 6.** Comparison between two clocks. (**a**) Row data of two clocks on different floors. (**b**) Row data of two clocks on the same ground floor. (**c**) Residual data of two clocks on different floors after processing. (**d**) Residual data of two clocks on the same ground floor after processing.

After these procedures, we obtained the 'valid' data. Then, we corrected some errors in time transfer, as discussed in Section 3. In the experiment, we just want to calculate the time interval difference (Δ*tAB <sup>T</sup>* ) caused by geopotential difference. In other words, we do not need to consider the synchronization accuracy of time, and we only need to consider the change rate of the time difference between the two sites. Therefore, the errors we need to consider are the changes over time rather than the fixed values over time in equipment delay, ionospheric delay, and the Sagnac effect. Among equipment delays, the time-varying part is the frequency drift of the clock. The output frequency of the atomic clock has a linear drift with time, which can be found in the zero-baseline experiment (Figure 6c), and could be corrected by the Zero-baseline experiment. In the ionospheric delay correction model, the TEC value changes in real-time. However, because the TEC values on the signal path are all the same (two antennas in one building), the influences on the result can be ignored. The time delay caused by the Sagnac effect is a fixed value because both the variation of the Earth's rotation speed and satellite orbits is minimal, which can be considered as constants under the current experimental accuracy [48]. Hence, the time delay caused by the Sagnac effect will not influence the experimental result at the present accuracy requirement.

Then, we obtained the data with an apparent linear trend and used a linear function to fit the data. As shown in Figure 6c and d, the blue curves denote original observations, the purple curves denote the data after removing the diurnal terms from the original data noted as residual data, and orange lines are linear fittings of the purple curves. In Figure 6c, the data show relatively intense jumps in the geopotential measurement. The reason for the jumps is that the temperature-control equipment was not particularly stable, which influenced the accuracy of the observations. Nevertheless, on the whole, this is in good agreement with the linear fitting line. This trend mainly includes the frequency drift of the clock, equipment error, and the time interval difference caused by geopotential differences. Therefore, we only need to eliminate the frequency drift of the clock and equipment error by zero-baseline measurement to obtain the results. In Figure 6d, in the whole zero-baseline comparison observation period, the data are relatively smooth and have a strong trend, which is conducive to the calculation of results.

#### **5. Results**

After calculation, the slopes of the geopotential comparison experiment and zerobaseline experiment are *kgeo* = 2.11639 × <sup>10</sup>−<sup>15</sup> and *kzero* = −0.93617 × <sup>10</sup><sup>−</sup>15, respectively. The slope of the zero-baseline measurement is the constant system shift; therefore, subtracting it from that of the geopotential comparison experiment could determine the difference of the clock running rates between CA and CB. To calculate the difference of the clock running rates ( <sup>Δ</sup>*tAB <sup>T</sup>* ), we differenced the two slopes:

$$\frac{\Delta t\_{AB}}{T} = k\_{\text{gco}} - k\_{zaro} = 3.0526 \times 10^{-15} \tag{9}$$

where *kremote* and *kzero* denote the slopes of the geopotential comparison experiment and zero-baseline experiment, respectively.

Based on Formulas (1) and (2), the measured height difference is Δ*H* = <sup>Δ</sup>*tAB <sup>T</sup> <sup>c</sup>*<sup>2</sup> *<sup>g</sup>* = 28.00 m. The residual standard deviation of the regression equation is expressed as [52]:

$$S = \sqrt{\frac{\sum\_{i=1}^{n} (\mathcal{Y}\_i - \mathcal{Y}\_i)^2}{n-2}} \tag{10}$$

where *n* is the total number of the sampling interval of the whole time series, *yi* denotes the *i*th observation, and *y*ˆ*<sup>i</sup>* denotes the *i*th fitting value. The uncertainty of the slope is given by [49]

$$\mu = \frac{\mathbb{S}}{\sqrt{\sum\_{i=1}^{n} (\boldsymbol{x}\_i - \overline{\boldsymbol{x}})^2}} \tag{11}$$

where S is the residual standard deviation of the regression equation, *xi* is the time of the *i*th observation, and *x* is the mean value of *x*.

Based on Formulas (10) and (11), we obtained the uncertainties of the slopes of the geopotential comparison experiment and zero-baseline experiment, written as *ugeo* = 0.26 × <sup>10</sup>−<sup>15</sup> and *uzero* = 0.52 × <sup>10</sup>−15, respectively. On the basis of the propagation law of errors, the uncertainty of the difference of the clock running rates ( <sup>Δ</sup>*tAB <sup>T</sup>* ) is

$$
\mu\_D = \sqrt{\mu\_{remove}^2 + \mu\_{zero}^2} = 0.58 \times 10^{-15} \tag{12}
$$

Finally, substituting the uncertainty values into Equations (1) and (2), we obtained the accuracy of the measurements, 5.4 m.

The relevant results are shown in Table 2.

**Table 2.** Results of geopotential comparison and zero-baseline comparison.


#### **6. Conclusions**

The experimental results in period 1 (geopotential comparison measurement) provide an uncertainty of the slope 0.26 × <sup>10</sup>−15, which has better accuracy than period 2 (zerobaseline comparison measurement) 0.52 × <sup>10</sup><sup>−</sup>15. The reason for this result may be that the observation time of period 1 is longer than period 2; a longer observation time is helpful to weaken the influence of observation noise on experimental results.

Based on TWSTFT observations, we determined the height difference between A and B as 28.0 ± 5.4 m. The bias between the measured value and the corresponding true value of 22.8 m is 5.2 m. Our current experimental accuracy is 5.4 m, which is not very high. The main reason is due to the accuracy limit of the hydrogen atomic clocks used, the stability of which is 5 × <sup>10</sup>−<sup>15</sup> and matches the final measurement accuracy. The current result is only exploratory research on this method, which may not be directly applied, but this result proves the feasibility of the method. Although the results reached expectations, there are still some problems worth further exploration. We cannot explain why there is strong periodicity in the observed data, but it is conjectured that these are connected with the atomic clock's performances and the relative motion of the satellite. The temperature could significantly influence the running rate of the atomic clock [53,54]. Due to the limitation of the experiment conditions, the ambient temperature of the atomic clock is not completely constant, and there is a certain fluctuation, which may be the reason for the jump in the observed data. In addition, the change of satellite orbit may also have an impact on the results. As described in Section 3.3, the satellite's orbit with a daily period, which is not fixed and changes every day.

Here, we employed two hydrogen MACs and TWSTFT technology to measure the height difference. The results of the experiment indicate that the accuracy of TWSTFT used in this experiment is sufficient for future applications. Our experimental results are preliminary. Further studies and experiments are needed to advance this research.

**Author Contributions:** Initiation and experiment design, W.S.; Conducting experiments, X.S., C.C., K.W., Z.S., and W.S.; Methodology, P.C. and W.S.; Software, P.C.; Data processing, P.C.; Validation, P.C., X.S., and C.C.; Formal analysis, P.C., W.S., and K.W.; Investigation, P.C. and W.S.; Original draft preparation, P.C.; Writing—review and editing, P.C., W.S., X.S., C.C., K.W. and Z.S.; Visualization, P.C.; Supervision, W.S.; Project administration, W.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (Grant nos. 41721003, 42030105, 41631072, 41804012, 41874023, 41974034), Space Station Project (2020)228 and Natural Science Foundation of Hubei Province (grant No. 2019CFB611).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data obtained in the study are available from the corresponding author upon reasonable request.

**Acknowledgments:** Thanks are given to the Beijing Institute of Radio Metrology and Measurement (BIRMM) for providing the experiment platform, and especially to the Frequency Measurement Laboratory of BIRMM for providing satellite two-way time transfer system. We would like to express our sincere thanks to three anonymous reviewers for their valuable comments and suggestions, which greatly improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

