Precision Joint RF Measurement of Inter-Satellite Range and Time Difference and Scalable Clock Synchronization for Multi-Microsatellite Formations

The rapid development of multi-satellite formations requires inter-satellite radio frequency (RF) measurement to be both precise and scalable. The navigation estimation of multi-satellite formations using a unified time reference demands the simultaneous RF measurement of the inter-satellite range and time difference. However, high-precision inter-satellite RF ranging and time difference measurements are investigated separately in existing studies. Different from the conventional two-way ranging (TWR) method, which is limited by its reliance on a high-performance atomic clock and navigation ephemeris, asymmetric double-sided two-way ranging (ADS-TWR)-based inter-satellite measurement schemes can eliminate such reliance while ensuring measurement precision and scalability. However, ADS-TWR was originally proposed for ranging-only applications. In this study, by fully exploiting the time-division non-coherent measurement characteristic of ADS-TWR, a joint RF measurement method is proposed to obtain the inter-satellite range and time difference simultaneously. Moreover, a multi-satellite clock synchronization scheme is proposed based on the joint measurement method. The experimental results show that when inter-satellite ranges are hundreds of kilometers, the joint measurement system has a centimeter-level accuracy for ranging and a hundred-picosecond-level accuracy for time difference measurement, and the maximum clock synchronization error was only about 1 ns.


Introduction
As the importance of microsatellite formations has increased [1,2], high-precision intersatellite ranging and clock synchronization, two inseparable key supporting technologies for the relative navigation of microsatellite formations, have received increasing attention. Inter-satellite clock synchronization can be established based on the inter-satellite time difference measurement, thereby providing a unified time reference for inter-satellite range measurement-based navigation solutions.
The inter-satellite range and time difference can be measured using the navigation satellite common view method or the inter-satellite RF autonomous measurement method [3,4]. The inter-satellite RF autonomous measurement method can effectively overcome the problems of the navigation satellite common view method, which has a high measurement accuracy only for low earth orbits and thus a limited scope of application. In addition, the inter-satellite RF autonomous measurement method does not rely on any external system and has a high degree of autonomy. Hence, it has been applied widely. Typical with these studies. First, the model was proposed for application scenarios with very low dynamics and had only a sub-meter measurement accuracy. Second, these studies focused on the localization estimation algorithm and did not consider the sources of measurement errors and provided no means to implement the joint measurement. Obviously, there is still a long way to go for the realization and application of space-borne high-precision joint measurement.
The other problem of the existing RF measurement research resides in the lack of a precision RF measurement method for miniaturized space-borne navigation systems. The traditional time division multiple access (TDMA)-based TWR method can effectively overcome the disadvantage of frequency division multiple access (FDMA) and code division multiple access (CDMA) based methods in terms of poor scalability and has been successfully applied to the inter-satellite links of the global positioning system (GPS) and Beidou navigation constellations [14]. However, the high measurement accuracy achieved with TWR in these navigation constellations is heavily reliant on high-performance atomic clocks and the assistance of navigation ephemeris, which are not available on microsatellite platforms due to limited resources.
In a previous work [15], authors proposed an ADS-TWR-based multi-satellite intersatellite ranging scheme, and high-ranging accuracy could be achieved with only a common miniaturized frequency source and without any external assistance. However, ADS-TWR was originally proposed for ranging-only applications, so a separate time difference measurement equipment is still needed to constitute a space-borne navigation system [16,17]. To this end, by fully exploiting the non-coherent measurement characteristic of ADS-TWR that potentially supports joint measurement, this study proposes a multi-satellite joint inter-satellite range and time difference measurement method.
It is noteworthy that inter-satellite time difference cannot be obtained directly as an inter-satellite range. Therefore, a critical problem of establishing a higher precision time difference reference needs to be solved to validate the time difference measurement performance. To address this issue, a CDMA-based high-precision time difference reference design approach is proposed to evaluate the time difference measurement performance of the joint measurement method. This reference is implemented on the same hardware platform of the joint measurement system and thus has the advantages of simplicity and miniaturization.
Regarding the clock synchronization method, the network time protocol (NTP) structure is generally used to achieve time synchronization for terrestrial wireless sensor networks. It works following three steps. First, part of the sensors in the network synchronizes 34 with the root node which possesses the reference time. Then, these sensors become new reference nodes. This process is passed layer by layer until all nodes in the network complete time synchronization [18]. However, this scheme has a problem in that errors will accumulate layer by layer and the synchronization accuracy depends heavily on the time accuracy of the root node. The consensus network [19] collects the clock information of all nodes, calculates the virtual clock of the current formation, and then each node is synchronized to the virtual clock. The consensus network scheme effectively solves the problems of NTP. However, the synchronization moment of all nodes must be strictly consistent, otherwise, the process will fail. From this point of view, this scheme is not suitable for TDMA. This is because, in a TDMA system, satellites can only communicate with other satellites in their own time slots for clock synchronization and cannot synchronize all the time like FDMA and CDMA. Meanwhile, in a consensus network, the clock failure of any node in the network will cause all nodes to lose synchronization with the virtual clock.
Based on the work of the joint measurement method, this paper further proposed an inter-satellite clock synchronization scheme for multi-satellite formations. This scheme can effectively alleviate the influence of frequency deviation among different satellites in a formation. It has the advantages that the synchronization error is not accumulated, its synchronization accuracy is not limited to a single node, and the clock synchronization is minimally affected by the failure of one satellite's clock.
The joint measurement and clock synchronization methods proposed in this paper enable large-scale formation microsatellites to achieve precision inter-satellite range and time difference measurement and clock synchronization with only one set of equipment.

System Model
The researchers [15,20] described a TDMA and ADS-TWR-based method for distributed multi-satellite measurement. Because the ADS-TWR method was originally proposed for ranging applications, they only analyzed and simulated the inter-satellite ranging performance. In this study, the feasibility of the method for inter-satellite time difference measurement will be analyzed. On this basis, a model for high-precision joint measurement will be established.
An ADS-TWR measurement involves three signals between any two satellites in a formation, as shown in Figure 1. Through pseudo-noise (PN) code measurement, satellite B can obtain the time messages of T A (t 1 ), T B (t 2 ), T A (t 5 ), and T B (t 6 ). Similarly, satellite A can obtain T B (t 3 ) and T A (t 4 ). T tof1 , T tof2 , and T tof3 are the times of flight. become new reference nodes. This process is passed layer by layer until all nodes in the network complete time synchronization [18]. However, this scheme has a problem in that errors will accumulate layer by layer and the synchronization accuracy depends heavily on the time accuracy of the root node. The consensus network [19] collects the clock information of all nodes, calculates the virtual clock of the current formation, and then each node is synchronized to the virtual clock. The consensus network scheme effectively solves the problems of NTP. However, the synchronization moment of all nodes must be strictly consistent, otherwise, the process will fail. From this point of view, this scheme is not suitable for TDMA. This is because, in a TDMA system, satellites can only communicate with other satellites in their own time slots for clock synchronization and cannot synchronize all the time like FDMA and CDMA. Meanwhile, in a consensus network, the clock failure of any node in the network will cause all nodes to lose synchronization with the virtual clock. Based on the work of the joint measurement method, this paper further proposed an inter-satellite clock synchronization scheme for multi-satellite formations. This scheme can effectively alleviate the influence of frequency deviation among different satellites in a formation. It has the advantages that the synchronization error is not accumulated, its synchronization accuracy is not limited to a single node, and the clock synchronization is minimally affected by the failure of one satellite's clock.
The joint measurement and clock synchronization methods proposed in this paper enable large-scale formation microsatellites to achieve precision inter-satellite range and time difference measurement and clock synchronization with only one set of equipment.

System Model
The researchers [15,20] described a TDMA and ADS-TWR-based method for distributed multi-satellite measurement. Because the ADS-TWR method was originally proposed for ranging applications, they only analyzed and simulated the inter-satellite ranging performance. In this study, the feasibility of the method for inter-satellite time difference measurement will be analyzed. On this basis, a model for high-precision joint measurement will be established.
An ADS-TWR measurement involves three signals between any two satellites in a formation, as shown in Figure 1. Through pseudo-noise (PN) code measurement, satellite B can obtain the time messages of T A (t 1 ), T B (t 2 ), T A (t 5 ), and T B (t 6 ). Similarly, satellite A can obtain T B (t 3 ) and T A (t 4 ). T tof1 , T tof2 , and T tof3 are the times of flight.  According to the above time messages, satellite B can construct the time intervals between the signal transmission and reception, namely T rdA , T rdB , T reA , and T reB as follows: The subscript rd represents the round-trip time and re represents the reply time.
Based on the ADS-TWR time measurement, the distance, R, between the two nodes can be obtained as follows [20]: where c is the light speed. Based on the range measurement and the times at the two nodes, the time difference, ∆T, between the two nodes at time t 3 can be obtained as follows: Equations (2) and (3) constitute the basic model for joint measurement. The high-precision joint measurement of inter-satellite range and time difference in multi-satellite formations requires in-depth analysis and modeling of the sources of measurement errors, and particular attention needs to be paid to the difference between the time difference measurement errors and ranging errors. The major sources of measurement errors of the multi-satellite joint measurement method proposed in this study include the frequency deviation of the frequency source and satellite motion-induced dynamics. In addition, because the measurement system uses PN code measurement signals, PN code phase jitter noise, receiver thermal noise, and dynamic stress error are inevitable, which are hereafter referred to collectively as phase tracking noise [11,21]. Moreover, the ionosphere delay and hardware delay need to be corrected. The influence of each error source on ranging and the corresponding compensation measures can be found in [20]. For the time difference measurement, the influence of each error source and the compensation measures are similar to ranging.
Since the first and third signals in the ADS-TWR method are symmetric with respect to the second signal, the range between the two satellites at time t 3 is used as the reference range R AB , with the corresponding time of flight designated as T tofAB . Considering the above sources of errors, with the ionosphere delay and hardware delay corrected (the detailed correction is shown in Appendix A), the total ADS-TWR ranging error and time difference measurement error, E R and E ∆T , can be numerically modeled as follows: where K A and K B represent the ratios of the onboard clock frequency to the nominal frequency for satellites A and B, respectively, n R is the phase tracking noise error of ranging, and n ∆T is the sum of the phase tracking noise errors at times T A (t 4 ) and T B (t 3 ).
The measurement parameters are listed in Table 1. These parameters are used in all experiments unless otherwise noted. Although longer coherent integration time can improve the carrier-to-noise ratio of the received signal, it will also increase the frequency error. Furthermore, the coherent integration time must be less than the data bit period to avoid the influence of code polarity change. Thus, a 1/4 data bit period is selected as the coherent integration time, which is 50 µs.

Error Caused by Frequency Deviation
The dominant component of the error caused by frequency deviation in Equation (5) can be expressed as: The simulation result of the time difference measurement is shown in Figure 2a, where the (K A +K B )/2 are 0.2 ppm, 0.05 ppm, and 0.02 ppm. For comparison, Figure 2b depicts the result of the range measurement, which is normalized to time by being divided by the speed of light.  Considering that the probability density functions of (K A +K B )/2 − 1 and (K A − K B )/2 are exactly the same, the influence of the frequency deviation on ranging and the time difference measurement is the same. For a typical inter-satellite distance of 200 km (the inter-satellite distance of the GRAIL formation is 175-225 km, and that of GRACE-FO is about 220 km) and a frequency accuracy of 0.01 ppm (which is available for a state-of-theart miniaturized frequency source), the time difference measurement error caused by the frequency deviation can be controlled under 10 ps.

Motion Error
Due to the fast movement of LEO satellites, T tof1 , T tof2 , and T tof3 are not equal, which leads to the satellite motion error. The time difference measurement error caused by satellite motion in Equation (5) can be expressed as: Considering that the probability density functions of (K A +K B )/2 − 1 and (K A − K B )/2 are exactly the same, the influence of the frequency deviation on ranging and the time difference measurement is the same. For a typical inter-satellite distance of 200 km (the inter-satellite distance of the GRAIL formation is 175-225 km, and that of GRACE-FO is about 220 km) and a frequency accuracy of 0.01 ppm (which is available for a state-ofthe-art miniaturized frequency source), the time difference measurement error caused by the frequency deviation can be controlled under 10 ps.

Motion Error
Due to the fast movement of LEO satellites, T tof1 , T tof2 , and T tof3 are not equal, which leads to the satellite motion error. The time difference measurement error caused by satellite motion in Equation (5) can be expressed as: Correspondingly, the range measurement error caused by satellite motion in Equation (4) can be expressed as: In a similar way to the analysis in [20], the following derivation will show how the satellite motion error arises when ADS-TWR is used for the time difference measurement. The relationship among the inter-satellite range, the times of flight, and the component of the satellite's absolute speed in the direction of the line connecting satellite A to B can be deduced as follows: where r(t 12 ) represents the geometric distance between the position of satellite A at moment T A (t 1 ) and the position of satellite B at moment T B (t 2 ), r(t 1 ) represents the distance between satellite A and B at moment T A (t 1 ), and v B12 represents the average absolute velocity in the baseline direction between satellite A and B from time T A (t 1 ) to T B (t 2 ), as shown in Figure 1. Other variables are named in the same way. The relationship among r(t 1 ), r(t 3 ), and r(t 5 ) is: where v AB represents the instantaneous relative velocity between the two satellites. The relative motion between the satellites is visualized in Figure 3. This diagram uses r(t − τ, t) to represent r(t 12 ), r(t 34 ), and r(t 56 ) in Equation (9), and uses r(t) to stand for r(t 1 ), r(t 3 ), and r(t 5 ). Substituting Equations (9) and (10) into Equations (7) and (8), we obtain M ∆T and M AT as follows: where K is the number of time slots in one measurement period, N AB is the number of time slots in the interval between time slot A and time slot B, T s is the duration of a time slot, and v AB13 and v AB35 denote the average value of v AB in the time intervals (t 1 , t 3 ) and (t 3 , t 5 ), respectively (the detailed derivation is shown in Appendix B).
where v AB represents the instantaneous relative velocity between the two satellites. The relative motion between the satellites is visualized in Figure 3. This diagram uses r(t-τ,t) to represent r(t 12 ), r(t 34 ), and r(t 56 ) in Equation (9), and uses r(t) to stand for r(t 1 ), r(t 3 ), and r(t 5 ). Substituting Equations (9) and (10) into Equations (7) and (8), we obtain M ∆T and M AT as follows: where K is the number of time slots in one measurement period, N AB is the number of time slots in the interval between time slot A and time slot B, T s is the duration of a time slot, and v AB13 and v AB35 denote the average value of v AB in the time intervals (t 1 , t 3 ) and (t 3 , t 5 ), respectively (the detailed derivation is shown in Appendix B). In fact, the direction of v A34 is opposite to that of v B12 and v B56 , and for a short time v AB35 and v AB13 can be considered to be equal. Therefore, the combinations of v A34 +v B12 , v A34 +v B56 and v AB35 − v AB13 enable ADS-TWR to reduce the motion error to a certain extent.
Equations (11) and (12) show that, compared with the ranging error, the time difference measurement error caused by satellite motion adds one error item, which is related to the transmission times and the satellites' absolute speed.
The following takes a practical multi-satellite formation as an example to illustrate the influence of satellite motion. A four-satellite circular formation is simulated (S0 is a virtual node) as shown in Figure 4. The orbital elements are displayed in Table 2. ors 2023, 23, x FOR PEER REVIEW In fact, the direction of v A34 is opposite to that of v B12 and v B56 , and v AB35 and v AB13 can be considered to be equal. Therefore, the combination v A34 + v B56 and v AB35 − v AB13 enable ADS-TWR to reduce the motion error tent.
Equations (11) and (12) show that, compared with the ranging error, ence measurement error caused by satellite motion adds one error item, w to the transmission times and the satellites' absolute speed.
The following takes a practical multi-satellite formation as an exam the influence of satellite motion. A four-satellite circular formation is sim virtual node) as shown in Figure 4. The orbital elements are displayed in T    a, e, i, Ω, w, and M 0 represent the semi-major axis, eccentricity, inclination, longitude of ascending node, argument of periapsis, and mean anomaly, respectively.
The simulation of the joint measurement error caused by satellite motion is shown in Figure 5. For the convenience of comparison, the range measurement results are normalized to time. As shown in this figure, the motion error in the time difference measurement is about 100 ps larger than that in the range measurement. However, the time difference error in the joint measurement remains at a hundred-picosecond level, which can meet the requirements of general microsatellite formations.

Phase Tracking Error
The phase tracking error includes the receiver thermal noise, dynamic stre transmitter code-phase jitter noise. For general satellite formations, such as cartwh circular formations, the relative dynamics among the satellites in the formation are and the dynamic stress error can be neglected. Therefore, we focus mainly on r thermal noise and transmitter code-phase jitter noise. We used a delay-locked loop for code tracking, and the standard deviation σ DLL of the receiver thermal noise er be expressed as [20]: is rier-noise ratio (CNR, Carrier-to-Noise Ratio, typically in units of dBHz); and D early-to-late correlation spacing (in chips).

Phase Tracking Error
The phase tracking error includes the receiver thermal noise, dynamic stress, and transmitter code-phase jitter noise. For general satellite formations, such as cartwheel and circular formations, the relative dynamics among the satellites in the formation are small, and the dynamic stress error can be neglected. Therefore, we focus mainly on receiver thermal noise and transmitter code-phase jitter noise. We used a delay-locked loop (DLL) for code tracking, and the standard deviation σ DLL of the receiver thermal noise error can be expressed as [20]: where B fe is the double-side front-end bandwidth; T c is the chip period; B L is the single-side equivalent loop bandwidth; T coh is the coherent integration time; C/N 0 is the carrier-noise ratio (CNR, Carrier-to-Noise Ratio, typically in units of dBHz); and D is the early-to-late correlation spacing (in chips). The transmitter code-phase jitter noise error σ θ can be considered to be linear in short time slots as: where f pn is the PN code rate; T u is time slot; and σ Allan is the short-term Allan deviation (ADEV) of the satellite frequency source, see [11] for a detailed analysis of the equation. The receiver thermal noise at times t 1 , t 3 , and t 5 are denoted as σ DLL1 , σ DLL3 , and σ DLL5 , respectively. The transmitter code-phase jitter noise affects T tof2 with the effect denoted as σ θtof2 . n ∆T is the sum of the phase tracking noise errors at times T A (t 4 ) and T B (t 3 ), and it can be expressed as: The expressions of n R are available in Appendix C. The simulation result of the phase tracking error is shown in Figure 6, with the trend of the curve similar to that of the range measurement in the above reference. Under high CNR (greater than 70 dBHz), the tracking error is less than 200 ps.

Summary of Error Modelling
The overall time difference error is shown in Figure 7. In summary, the joint m urement system has a centimeter-level accuracy for ranging and a hundred-picoseco level accuracy for the time difference measurement.

Summary of Error Modelling
The overall time difference error is shown in Figure 7. In summary, the joint measurement system has a centimeter-level accuracy for ranging and a hundred-picosecond-level accuracy for the time difference measurement.

Summary of Error Modelling
The overall time difference error is shown in Figure 7. In summary, the joint measurement system has a centimeter-level accuracy for ranging and a hundred-picosecondlevel accuracy for the time difference measurement.

Design of Time Difference Reference
This section describes the design of the time difference reference. For a two-node formation, the data from the transmitters of the two nodes modulate different PN code sequences. Each receiver of the two nodes has two independent receiving channels which receive the TDMA signal and CDMA signal, respectively, as shown in Figure 8. The CDMA channel keeps tracking the received signal and calculates the time difference between the two nodes using the time difference measurement method proposed in [11].

Design of Time Difference Reference
This section describes the design of the time difference reference. For a two-node formation, the data from the transmitters of the two nodes modulate different PN code sequences. Each receiver of the two nodes has two independent receiving channels which receive the TDMA signal and CDMA signal, respectively, as shown in Figure 8. The CDMA channel keeps tracking the received signal and calculates the time difference between the two nodes using the time difference measurement method proposed in [11]. The TDMA channel works only in the time slot of the node and measures the range and time difference using ADS-TWR.  Different from TDMA, the CDMA-based time difference measurement is temporally continuous. Thus, the frequency deviation caused error and motion caused error are small and the total measurement error of time difference measurement using CDMA is much less than that using TDMA. In a ground static environment, the motion error is zero for both CDMA and TDMA; however, the phase tracking noise is related to the loop bandwidth of the delay-locked loop (DLL) [20], and the total measurement error of CDMA can be controlled much less than that of TDMA by decreasing the loop bandwidth of CDMA. Therefore, the reference designed as such can be used for the experimental verification of the time difference measurement performance of the proposed joint measurement method. Because the CDMA system is realized on the same hardware as the TDMA joint measurement system, the overall system had the prominent advantages of simplicity and miniaturization.

Multi-Satellite Clock Synchronization
Considering the advantages and disadvantages of the NTP and the consensus clock Different from TDMA, the CDMA-based time difference measurement is temporally continuous. Thus, the frequency deviation caused error and motion caused error are small and the total measurement error of time difference measurement using CDMA is much less than that using TDMA. In a ground static environment, the motion error is zero for both CDMA and TDMA; however, the phase tracking noise is related to the loop bandwidth of the delay-locked loop (DLL) [20], and the total measurement error of CDMA can be controlled much less than that of TDMA by decreasing the loop bandwidth of CDMA. Therefore, the reference designed as such can be used for the experimental verification of the time difference measurement performance of the proposed joint measurement method. Because the CDMA system is realized on the same hardware as the TDMA joint measurement system, the overall system had the prominent advantages of simplicity and miniaturization.

Multi-Satellite Clock Synchronization
Considering the advantages and disadvantages of the NTP and the consensus clock synchronization schemes, as well as the distributed multi-satellite measurement scheme in [10], we propose an innovative clock synchronization scheme which is shown in Figure 9. The synchronization process can be described as follows: each satellite only transmits a signal in its own time slot and receives other satellites' signals in the other time slots. Take three satellite (numbered S1, S2, S3) formations as an example. When S1 is in the transmission state, the remaining two satellites in the formation measure the time difference between themselves and S1 based on the joint measurement method, thus completing clock synchronization with S1. When the time slot of S1 ends, satellite S2 switches into the transmission state, so S2 becomes the new root node, and the remaining two satellites synchronize with S2. In the end, clock synchronization is established and maintained among all satellites in the formation.
be controlled much less than that of TDMA by decreasing the loop bandwidth of CDMA. Therefore, the reference designed as such can be used for the experimental verification of the time difference measurement performance of the proposed joint measurement method. Because the CDMA system is realized on the same hardware as the TDMA joint measurement system, the overall system had the prominent advantages of simplicity and miniaturization.

Multi-Satellite Clock Synchronization
Considering the advantages and disadvantages of the NTP and the consensus clock synchronization schemes, as well as the distributed multi-satellite measurement scheme in [10], we propose an innovative clock synchronization scheme which is shown in Figure  9. The synchronization process can be described as follows: each satellite only transmits a signal in its own time slot and receives other satellites' signals in the other time slots. Take three satellite (numbered S1, S2, S3) formations as an example. When S1 is in the transmission state, the remaining two satellites in the formation measure the time difference between themselves and S1 based on the joint measurement method, thus completing clock synchronization with S1. When the time slot of S1 ends, satellite S2 switches into the transmission state, so S2 becomes the new root node, and the remaining two satellites synchronize with S2. In the end, clock synchronization is established and maintained among all satellites in the formation. S1 S3 S2 S1 is in transmission state S2 is in transmission state S3 is in transmission state S1 S3 S2 S1 S3 S2 Figure 9. Clock synchronization scheme.
Similar to the consensus clock synchronization scheme, the new clock synchronization is divided into two steps as shown in Figure 10. First, each satellite measures the time difference between itself and the satellite in the transmission state and corrects the clock phase by modifying the NCO (numerically controlled oscillator) phase register. Second, the frequency difference between the two satellites is obtained by calculating the difference between the current and last time difference measurements and dividing it by the measurement cycle time, which can be compensated by modifying the frequency control word of the NCO. Similar to the consensus clock synchronization scheme, the new clock synchronization is divided into two steps as shown in Figure 10. First, each satellite measures the time difference between itself and the satellite in the transmission state and corrects the clock phase by modifying the NCO (numerically controlled oscillator) phase register. Second, the frequency difference between the two satellites is obtained by calculating the difference between the current and last time difference measurements and dividing it by the measurement cycle time, which can be compensated by modifying the frequency control word of the NCO.

Clock Synchronization Performance Analysis
The phase correction error depends on the time difference measurement error. Since a hundred-picosecond-level accuracy of time difference measurement can be achieved, the same level of accuracy can also be expected for phase correction.
We start with the case that the clock frequency is not corrected. The maximum synchronization error can be expressed as:

Clock Synchronization Performance Analysis
The phase correction error depends on the time difference measurement error. Since a hundred-picosecond-level accuracy of time difference measurement can be achieved, the same level of accuracy can also be expected for phase correction.
We start with the case that the clock frequency is not corrected. The maximum synchronization error can be expressed as: where T update is the synchronization cycle equal to 5 s and the frequency accuracy is 0.01 ppm. According to Equation (16), when satellite A is in the transmission state, the maximum synchronization error T syn_error of satellite B will reach hundreds of nanoseconds. Thus, frequency correction is necessary in order to achieve high-precision inter-satellite clock synchronization. Frequency correction is performed by calculating the ratio of the frequencies of the two satellites and feeding the result to the frequency control word of the NCO. The modified parameters are calculated as in Equation (17).
where ∆ f (i) is the frequency difference at time i, ∆T(i) is the time difference measurement at time i, f local is the local clock frequency of the satellite in the receiving state, and ∆ f NCO is the NCO frequency control word error. Combining Equations (16) and (17), we can obtain the expression of the synchronization error after frequency correction as shown in Equation (18).
where f 0 is the nominal frequency equal to 40 MHz. The frequency correction error calculated in Equation (18) is affected by three factors. The first one is the frequency source stability, and its effect is only tens of picoseconds when the frequency stability is 0.01 ppb, according to Equation (16). The second one is the time difference measurement error of the joint measurement. When the parameters are used as in Section 2, the difference of the time difference measurement error at time i and i − 1 is about 317 ps, which will bring a frequency correction error of 2.5 mHz according to Equation (17). Substituting ∆ f (i) in Equation (18) with 2.5 mHz will result in a synchronization error of 317 ps as well, since f local and f 0 are almost the same. The third one is the NCO frequency word error. The worst error of the frequency control word in this paper is 9.313 MHz. By substituting it into Equation (18), we find that the synchronization error incurred by the NCO frequency word error is 1.1641 ns. In summary, the worst synchronization error can reach about 2 ns after frequency correction.

Simulation of Multi-Satellite Clock Synchronization Scheme
The simulation parameters are displayed in Table 3, where K i represents the frequency accuracy of satellite i; eNCO i represents the NCO frequency word error of satellite i; and eTD ij represents the standard deviation of the time difference measurement error between satellite i and j. Take the clock synchronization between satellites 1 and 2 as an example, K 1 and K 2 in Table 3 represent K A and K B in Equation (18), respectively. Additionally, (eNCO 1 − eNCO 2 ) is equivalent to ∆ f NCO in Equation (17), and eTD 12 corresponds to the standard deviation of the errors of ∆T(i) and ∆T(i − 1). The simulation results of the performance of the clock synchronization scheme are shown in Figure 11. During the first time slot (0~5 s), clock synchronization has not been established, so the clock deviation is more obvious. The initial clock time of satellite 1, satellite 2, and satellite 3 are 5.12 + 2 × 10 −7 s, 5.12 s, and 5.12 − 2 × 10 −7 s, respectively.

Experimental Platform
An inter-satellite RF measurement transceiver was designed and implemented, based on which the measurement system was constructed. The S-band RF measurement transceiver consisted of a digital signal processing module and two RF front-ends for reception and transmission, respectively. Since the system is based on the TDMA mechanism, the receiver and transmitter operate on the same carrier frequency and can be freely switched on and off. Figure 12 shows a prototype of the RF measurement transceiver.  Figure 11a shows the clock deviation between Satellites 1 and 3 (noted as CD 13 ), Satellites 2 and 3 (noted as CD 23 ), and Satellites 1 and 2 (noted as CD 12 ) without phase correction and frequency correction for each satellite. Figure 11b shows the clock deviation with phase correction only. It can be seen that there is still a hundred-ns level gap among the satellite clocks due to the absence of frequency error correction. Figure 11c shows the clock deviation with both phase correction and frequency correction, which remains at only an ns level, demonstrating the best synchronization performance.

Experimental Platform
An inter-satellite RF measurement transceiver was designed and implemented, based on which the measurement system was constructed. The S-band RF measurement transceiver consisted of a digital signal processing module and two RF front-ends for reception and transmission, respectively. Since the system is based on the TDMA mechanism, the receiver and transmitter operate on the same carrier frequency and can be freely switched on and off. Figure 12 shows a prototype of the RF measurement transceiver. The measurement processing module of the transceiver was developed based on an FPGA (field programmable gate array). The structure of the processing module is displayed in Figure 13 and its main specifications are listed in Table 4. The local PN code count and phase can be extracted from the local PN code generator in the transmitter to form the signal reception time. The TDMA receiving channel obtains the received PN code chip count and phase through the tracking loop to form the signal transmission time and then obtains the range and time difference measurements by using the ADS-TWR method. The measurement processing module of the transceiver was developed based on an FPGA (field programmable gate array). The structure of the processing module is displayed in Figure 13 and its main specifications are listed in Table 4. The local PN code count and phase can be extracted from the local PN code generator in the transmitter to form the signal reception time. The TDMA receiving channel obtains the received PN code chip count and phase through the tracking loop to form the signal transmission time and then obtains the range and time difference measurements by using the ADS-TWR method. Sensors 2023, 23, x FOR PEER REVIEW 16 of 26 Figure 13. Joint measurement module.

Two-Node Experiment
The test platform consists of two RF transceivers and other auxiliary devices including miniaturized oven-controlled crystal oscillators (OCXOs), RF cables, attenuators, an RS232 cable, and a computer. As shown in Figure 14, the attenuator was connected to both transceivers via RF cables to adjust the signal strength and protect the transceivers. The results of the joint measurement and clock synchronization were sent to the computer via the RS232. The transmission speed in the RF cable needs to be calibrated before performance evaluation, and the calibration result was 20,660,000 m/s. At the same time, the hardware delay could be calculated as 3.294 μs.

Two-Node Experiment
The test platform consists of two RF transceivers and other auxiliary devices including miniaturized oven-controlled crystal oscillators (OCXOs), RF cables, attenuators, an RS232 cable, and a computer. As shown in Figure 14, the attenuator was connected to both transceivers via RF cables to adjust the signal strength and protect the transceivers. The results of the joint measurement and clock synchronization were sent to the computer via the RS232. The transmission speed in the RF cable needs to be calibrated before performance evaluation, and the calibration result was 20,660,000 m/s. At the same time, the hardware delay could be calculated as 3.294 µs.

Joint Measurement Experimental Verification
After compensating for the hardware delays, we obtained the overall measurement errors compared with the theoretical values as shown in Figure 15. It can be seen that the experimental accuracy is well consistent with the theoretical analysis. At a high carrier-tonoise ratio of 70 dBHz, an overall ranging error of 4.52 cm and a time difference measurement error of 233 ps can be achieved. It is noteworthy that, as revealed by the preceding analysis, the error of time difference measurement should be equal to that of ranging in a ground static environment. However, the former seems larger than the latter, as shown in Figure 15. For example, the time difference measurement error and normalized ranging error are 233 ps and 151 ps, respectively, at a carrier-to-noise ratio of 70 dBHz. This is because the measurement signals were transmitted in RF cables rather than in a vacuum. They are actually equivalent when the calculation is performed with the transmission speed in RF cables.

Joint Measurement Experimental Verification
After compensating for the hardware delays, we obtained the overall measurement errors compared with the theoretical values as shown in Figure 15. It can be seen that the experimental accuracy is well consistent with the theoretical analysis. At a high carrier-to-noise ratio of 70 dBHz, an overall ranging error of 4.52 cm and a time difference measurement error of 233 ps can be achieved.

Joint Measurement Experimental Verification
After compensating for the hardware delays, we obtained the overall measurement errors compared with the theoretical values as shown in Figure 15. It can be seen that the experimental accuracy is well consistent with the theoretical analysis. At a high carrier-tonoise ratio of 70 dBHz, an overall ranging error of 4.52 cm and a time difference measurement error of 233 ps can be achieved. It is noteworthy that, as revealed by the preceding analysis, the error of time difference measurement should be equal to that of ranging in a ground static environment. However, the former seems larger than the latter, as shown in Figure 15. For example, the time difference measurement error and normalized ranging error are 233 ps and 151 ps, respectively, at a carrier-to-noise ratio of 70 dBHz. This is because the measurement signals were transmitted in RF cables rather than in a vacuum. They are actually equivalent when the calculation is performed with the transmission speed in RF cables.  It is noteworthy that, as revealed by the preceding analysis, the error of time difference measurement should be equal to that of ranging in a ground static environment. However, the former seems larger than the latter, as shown in Figure 15. For example, the time difference measurement error and normalized ranging error are 233 ps and 151 ps, respectively, at a carrier-to-noise ratio of 70 dBHz. This is because the measurement signals were transmitted in RF cables rather than in a vacuum. They are actually equivalent when the calculation is performed with the transmission speed in RF cables. Figure 16a compares the clock synchronization error with and without frequency correction. Apparently, the synchronization error between two phase corrections without frequency correction can reach hundreds of nanoseconds, which is much larger than that with frequency correction. Figure 16b is a magnified display of the clock synchronization results with frequency correction. It can be seen that the maximum error between two phase corrections is only about 1 ns, which is in accordance with the theoretical analysis.  Figure 16a compares the clock synchronization error with and without frequency correction. Apparently, the synchronization error between two phase corrections without frequency correction can reach hundreds of nanoseconds, which is much larger than that with frequency correction. Figure 16b is a magnified display of the clock synchronization results with frequency correction. It can be seen that the maximum error between two phase corrections is only about 1 ns, which is in accordance with the theoretical analysis.

Multi-Node Experiment
In [10] It was demonstrated that the measurement accuracy hardly deteriorated as the number of nodes increased, even when the number of nodes reached tens. A threenode experimental platform was established, as shown in Figure 17, to verify the performance and scalability of the proposed joint RF measurement method and clock synchronization scheme for multi-microsatellite formations.
The experimental platform was constructed as follows: the three nodes were connected with fixed attenuators and splitters. Node A and node B were connected with 3m-long RF cables, node A and node C were connected with 1-m-long RF cables, and node B and node C were connected with 5-m-long RF cables.

Multi-Node Experiment
In [10] It was demonstrated that the measurement accuracy hardly deteriorated as the number of nodes increased, even when the number of nodes reached tens. A three-node experimental platform was established, as shown in Figure 17, to verify the performance and scalability of the proposed joint RF measurement method and clock synchronization scheme for multi-microsatellite formations.  Figure 16a compares the clock synchronization error with and without frequen correction. Apparently, the synchronization error between two phase corrections witho frequency correction can reach hundreds of nanoseconds, which is much larger than th with frequency correction. Figure 16b is a magnified display of the clock synchronizati results with frequency correction. It can be seen that the maximum error between tw phase corrections is only about 1 ns, which is in accordance with the theoretical analys

Multi-Node Experiment
In [10] It was demonstrated that the measurement accuracy hardly deteriorated the number of nodes increased, even when the number of nodes reached tens. A thre node experimental platform was established, as shown in Figure 17, to verify the perf mance and scalability of the proposed joint RF measurement method and clock synchr nization scheme for multi-microsatellite formations.
The experimental platform was constructed as follows: the three nodes were co nected with fixed attenuators and splitters. Node A and node B were connected with m-long RF cables, node A and node C were connected with 1-m-long RF cables, and no B and node C were connected with 5-m-long RF cables. The experimental platform was constructed as follows: the three nodes were connected with fixed attenuators and splitters. Node A and node B were connected with 3-m-long RF cables, node A and node C were connected with 1-m-long RF cables, and node B and node C were connected with 5-m-long RF cables.

Joint Measurement Experimental Verification
The test results of the joint measurement experiment are shown in Figure 18. After the circuit delay was compensated, the standard deviations of the errors of the range and time difference measurements between node A and node B were 4.49 cm and 224.29 ps, respectively. The standard deviations of the errors of the range and time difference measurements between node A and node C were 4.64 cm and 217.03 ps, respectively. The standard deviations of the errors of the range and time difference measurements between node B and node C were 4.49 cm and 226.59 ps, respectively. The above results show that with the increase in the number of nodes, there is no significant degradation of the measurement performance, which accords with the theoretical analysis and validates the scalability of the proposed joint measurement method. The test results of the joint measurement experiment are shown in Figure 18. After the circuit delay was compensated, the standard deviations of the errors of the range and time difference measurements between node A and node B were 4.49 cm and 224.29 ps, respectively. The standard deviations of the errors of the range and time difference measurements between node A and node C were 4.64 cm and 217.03 ps, respectively. The standard deviations of the errors of the range and time difference measurements between node B and node C were 4.49 cm and 226.59 ps, respectively. The above results show that with the increase in the number of nodes, there is no significant degradation of the measurement performance, which accords with the theoretical analysis and validates the scalability of the proposed joint measurement method.
(a) (b) (c) Figure 18. Multi-node joint measurement experimental results: (a) measurements between node A and B; (b) measurements between node A and C; and (c) measurements between node B and C.

Clock Synchronization Experimental Verification
Since each satellite performs clock synchronization independently, the results of multi-node clock synchronization are similar to that of two-node clock synchronization. Figure 19 shows the clock synchronization errors when different satellite nodes acted as the master nodes, i.e., when different satellites were in the transmission state. The sampling period was 10 s. It can be seen that in the multi-node case, the system can still work well and maintain synchronization errors of less than 1 nanosecond. Figure 18. Multi-node joint measurement experimental results: (a) measurements between node A and B; (b) measurements between node A and C; and (c) measurements between node B and C.

Clock Synchronization Experimental Verification
Since each satellite performs clock synchronization independently, the results of multinode clock synchronization are similar to that of two-node clock synchronization. Figure 19 shows the clock synchronization errors when different satellite nodes acted as the master nodes, i.e., when different satellites were in the transmission state. The sampling period was 10 s. It can be seen that in the multi-node case, the system can still work well and maintain synchronization errors of less than 1 nanosecond. Sensors 2023, 23, x FOR PEER REVIEW 20 of 26 (a) (b) (c) Figure 19. Multi-node clock synchronization experimental results: (a) results of node B and C when node A is the master node; (b) results of node A and C when node B is the master node; and (c) results of node A and B when node C is the master node.

Conclusions
In this study, a method was proposed for the joint RF measurement of inter-satellite range and time difference for multi-satellite formations. Based on this method, a new clock synchronization scheme was proposed. The experimental results showed that the joint measurement system had a centimeter-level accuracy for ranging and a hundred-picosecond-level accuracy for time difference measurement, thus verifying the correctness of the theoretical measurement model. The experimental results also demonstrated that the maximum clock synchronization error was only about 1 ns. This research effectively solves the existing problem in the literature that a highly accurate inter-satellite range and time difference measurement cannot be obtained simultaneously using only one piece of equipment and provides a foundation for the miniaturized and high-precision RF navigation of multi-satellite formations.   . Multi-node clock synchronization experimental results: (a) results of node B and C when node A is the master node; (b) results of node A and C when node B is the master node; and (c) results of node A and B when node C is the master node.

Conclusions
In this study, a method was proposed for the joint RF measurement of inter-satellite range and time difference for multi-satellite formations. Based on this method, a new clock synchronization scheme was proposed. The experimental results showed that the joint measurement system had a centimeter-level accuracy for ranging and a hundredpicosecond-level accuracy for time difference measurement, thus verifying the correctness of the theoretical measurement model. The experimental results also demonstrated that the maximum clock synchronization error was only about 1 ns. This research effectively solves the existing problem in the literature that a highly accurate inter-satellite range and time difference measurement cannot be obtained simultaneously using only one piece of equipment and provides a foundation for the miniaturized and high-precision RF navigation of multi-satellite formations.