**1. Introduction**

The inter-satellite link (ISL) is not only an important feature and technological innovation of BDS, but also an important means to optimize the distribution of regional station, improve the accuracy of spatial signals and increase the frequency of message injection [1–3]. BDS-3 has provided the global services officially on 31 July 2020, and all of the satellites are equipped with ISL [4]. Based on the Time Division Multiple Access (TDMA) technique and the Ka-band technology of BDS ISL strategy, any visible satellites can achieve real-time inter-satellite measurement and communication by the multipoint-to-multipoint time-division two-way links [5,6].

In contrast with GPS and GALILEO systems that update ephemeris and clock offset parameters mainly through the global uniform distribution of stations [7,8], BDS adopts the combined time synchronization strategy which determines satellite clock offset parameters by the satellite-ground link (SGL) for visible satellite and inter-satellite link (ISL) for invisible satellite [9,10]. This strategy was first published in 2018, and usually called as the "one-hop" reduction method [11]. According to this method, when the satellite is visible by the Master Control Center (MCC), the clock offset is

**Citation:** Wang, D.; Guo, R.; Liu, L.; Yuan, H.; Li, X.; Pan, J.; Tang, C. A Method of Whole-Network Adjustment for Clock Offset Based on Satellite-Ground and Inter-Satellite Link Observations. *Remote Sens.* **2022**, *14*, 5073. https://doi.org/ 10.3390/rs14205073

Academic Editor: Yunbin Yuan

Received: 11 August 2022 Accepted: 6 October 2022 Published: 11 October 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/).

measured through the SGL L-band satellite-ground two-way time synchronization; when the satellite is invisible by MCC, the clock offset is calculated by one node satellite (which is visible to both MCC and the invisible node satellite) of the SGL L-band twoway time synchronization and these two satellites of the ISL Ka-band two-way time synchronization, thus clock offsets of invisible satellite can be calculated through the "one-hop" reduction mode by one node visible satellite [11,12]. Compared with only SGL, the tracking coverage is extended by 40% by using ISL [11], the prediction error can be reduced from 3 ns to 1 ns [13], and the orbit result is also greatly improved [14–16].

Some achievements have been published to estimate the affections of the systematic error. Yan [17] and Wang [18] estimate the satellite antenna phase center offsets (PCOs), phase variations (PVs) and the differential code bias (DCB) of BeiDou-3 satellites, and has the results that PCOs and PVs is different with satellites, and DCB value is related to the receiver type. Moreover, Liu [19] proposes a network adjustment model by the correction of inter-satellite clock offsets to detect and analyze the closed residuals, and has the results that the random noise of the inter-satellite clock corrections is reduced by 30% to 50%. These achievements have evaluated the performance of current system and given the application conclusions that the signal in space accuracy has reached up to 0.5 m [2]. However, systematic errors of the measured SGL clock offsets and the calculated ISL clock offsets using different node satellite exist, which affect the accuracy of the time synchronization [12,20].

This paper further studies the whole-network adjustment method by multi-source clock offset of the satellite-ground and inter-satellite observation. Taking SGL and ISL data of all orbiting satellites as a whole, this new method is used to realize the satelliteground and inter-satellite time synchronization in the entire constellation. The propose is to optimize the time synchronization strategy of BDS, minimize the impact of different devices, and improve the accuracy of BDS broadcast clock offset.

### **2. Materials and Methods**

#### *2.1. Time Synchronization Principle*

#### 2.1.1. L-Band Satellite-Ground Two-Way Time Synchronization

BDS achieves satellite-ground time synchronization using L-band two-way time comparison. The satellite and ground station, respectively, generate and broadcast pseudocode ranging signals based on the control of local clocks. Specifically, the ground station obtains the downlink pseudorange at the time corresponding to the ground 1 pulse per second (pps), and the satellite obtains the uplink pseudorange at the time corresponding to the satellite 1 pps. Meanwhile, the satellite sends the observed value of uplink pseudorange to the ground station, and then the ground station calculates the difference between the locally measured downlink pseudorange and the received uplink pseudorange to obtain the clock offset of satellite relative to the ground station. It lays a foundation for completing the satellite-ground two-way time comparison. The principle of L-band satellite-ground two-way time comparison is shown in Figure 1.

The satellite-ground two-way time comparison method is used to obtain the observed satellite-ground clock offset, which can be expressed by a second-order polynomial:

$$
\Delta T\_{\text{\tiny\text\text\text\text\text\textheadtext\textheadtext\text\text\text\text\text\textheadtext\text\text\text\text\text\text\text\text\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbullet}{2c}\!(\rho\_{\text{\tiny\text\texttext\texthead}}-\rho\_{\text{d}})-\frac{1}{2c}(R\_{\text{\tiny\text\textheadtext\textheadtext\text\text\textheadtext\text\text\texthead\text\text\text\text\text\texthead\text\text\text\text\text\text\text\text\text\textbig\textbig\textbig\textbig\textbig\textbig\textbig\textbullet\}{}^{\text\text\text\text\texthead\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\textqu\}}{} \tag{1}
$$

where Δ*Ts*, Δ*Tg*, respectively, denote the satellite clock offset and the ground station clock offset; *ρu*, *ρ<sup>d</sup>* represent the uplink and downlink pseudorange, respectively; *Ru*, *Rd* indicate the theoretical values of uplink and downlink geometric distances, respectively; *τrel* represents the relativistic effect delay in the signal propagation path; *ε* indicates the synthetic measurement noise of the satellite-ground two-way link; and *c* represents the velocity of light.

The satellite-ground two-way time comparison method is used to obtain the observed value of satellite-ground clock offset. The satellite-ground clock offset of satellite *i* can be expressed by a second-order polynomial [21]:

$$
\Delta T\_i(t) = a\_0^i + a\_1^i(t - t\_0) + a\_2^i(t - t\_0)^2 \tag{2}
$$

where Δ*Ti*(*t*) represents the discrete point of satellite-ground clock offset of satellite *i*; *t* denotes the observation time of satellite-ground clock offset; *t*<sup>0</sup> refers to the reference time of clock offset parameter; *a<sup>i</sup>* <sup>0</sup>, *<sup>a</sup><sup>i</sup>* <sup>1</sup>, *<sup>a</sup><sup>i</sup>* <sup>2</sup> represent the clock offset, clock rate, and clock drift, respectively.

#### 2.1.2. Ka-Band Inter-Satellite Two-Way Time Synchronization

It is similar to the satellite-ground two-way time comparison method. If satellite *i* and satellite *j* send time signals to each other at the clock *Ti* and *Tj*, the signal sent by satellite *i* is received by satellite *j* after the delay from satellite *i* to satellite *j*. On this basis, the pseudorange from satellite *i* to satellite *j* can be measured. Similarly, the signal sent by satellite *j* is received by satellite *i* after the delay from satellite *j* to satellite *i*, so that the pseudorange from satellite *j* to satellite *i* can be measured. Next, satellite *i* and satellite *j* exchange their respective observation results, and finally calculate the relative clock offset between them. The detailed principle of Ka-band inter-satellite two-way time comparison is shown in Figure 2.

**Figure 2.** Schematic diagram of inter-satellite time synchronization.

Similar with the satellite-ground two-way time synchronization, the inter-satellite time synchronization of navigation satellites is based on pseudorange measurement of Ka-band two-way time synchronization, and the inter-satellite clock offset can be obtained as:

$$
\Delta T\_{\dot{j}} - \Delta T\_{\dot{i}} = \frac{1}{2c} (\rho\_{\dot{i}\dot{j}} - \rho\_{\dot{j}\dot{i}}) - \frac{1}{2c} (R\_{\dot{i}\dot{j}} - R\_{\dot{j}\dot{i}}) + \delta \tag{3}
$$

where Δ*Ti*, Δ*Tj* represent the clock offsets of satellite *i* and satellite *j*, respectively; *ρij* refers to the pseudorange of satellite *j* receiving from satellite *i*; *ρji* refers to the pseudorange of satellite *i* receiving from satellite *j*; *Rij*, *Rji* respectively, represent the theoretical values of two-way geometric distances; *δ* indicates the synthetic measurement noise of the intersatellite two-way link; and *c* represents the velocity of light.

The inter-satellite two-way time comparison method is used to obtain the observed inter-satellite clock offset, which can be expressed by a second-order polynomial [21]:

$$
\Delta T\_{ij}(t) + v\_{ij} = a\_0^j + a\_1^j(t - t\_0) + a\_2^i(t - t\_0)^2 - \\\lambda \tag{4}
$$

$$
\begin{pmatrix}
a\_0^i + a\_1^i(t - t\_0) + a\_2^i(t - t\_0)^2 \\
\end{pmatrix} \tag{5}
$$

where Δ*Tij*(*t*) represents the discrete point of inter-satellite clock offset of satellite *j* relative to satellite *i*; *t* denotes the observation time of inter-satellite clock offset; *t*<sup>0</sup> refers to the reference time of clock offset parameter; *vij* indicates the Ka-band inter-satellite link delay calibrated based on the L-band satellite-ground link; *a<sup>i</sup>* <sup>0</sup>, *<sup>a</sup><sup>i</sup>* <sup>1</sup>, *<sup>a</sup><sup>i</sup>* <sup>2</sup> represent the clock offset, clock rate and clock drift, respectively.

#### *2.2. Whole-Network Adjustment Method*

The Whole-network adjustment method uses multi-source clock offset data to obtain the clock offset parameters of satellite-ground and inter-satellite of all satellites. Assuming there is *m* satellites and *n* inter-satellite links, the *m* satellite-ground clock offsets and *n* inter-satellite clock offsets can be performed by whole-network adjustment method, and the equation can be abbreviated as:

$$L = \mathbb{C}\mathfrak{x} - \mathfrak{v} \tag{5}$$

where *L* = Δ*Ti*(*t*) Δ*Tij*(*t*) (*m*+*n*)∗1 represents the discrete point of *m* satellite-ground and *n*

inter-satellite clock offset; *C* = *Ci Cij* (*m*+*n*)∗3*m* refers to the coefficient matrix of observation

$$\text{time}; \mathbf{x} = \begin{pmatrix} \mathbf{x}^1 \\ \vdots \\ \mathbf{x}^i \\ \vdots \\ \vdots \\ \mathbf{x}^m \end{pmatrix}\_{3m+1} \\ \text{in indicates the column vector of clock offset parameters;} \\ \mathbf{x}^i = \begin{pmatrix} a\_0^i \\ a\_1^i \\ \vdots \\ a\_2^i \end{pmatrix} \\ \tag{2}$$

denotes the clock offset parameter of each satellite; *v* = 0 *vij* (*m*+*n*)∗1 refers to the link delay error [22]. The coefficient matrix is expressed as follows:

$$\mathbf{C}\_{i} = \begin{bmatrix} 1 & t - t\_{0} & (t - t\_{0})^{2} & 0 & 0 & 0 & \cdots & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & t - t\_{0} & (t - t\_{0})^{2} & \cdots & 0 & 0 & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 1 & t - t\_{0} & (t - t\_{0})^{2} \end{bmatrix}\_{m \times 3m}$$

$$\mathbf{C}\_{ij} = \begin{bmatrix} 0 & \cdots & -1 & -(t - t\_{0}) & -(t - t\_{0})^{2} & 0 & \cdots & 1 & (t - t\_{0}) & (t - t\_{0})^{2} & 0 & \cdots\\ 0 & \cdots & -1 & -(t - t\_{0}) & -(t - t\_{0})^{2} & 0 & \cdots & 1 & (t - t\_{0}) & (t - t\_{0})^{2} & 0 & \cdots\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & \cdots & -1 & -(t - t\_{0}) & -(t - t\_{0})^{2} & 0 & \cdots & 1 & (t - t\_{0}) & (t - t\_{0})^{2} & 0 & \cdots \end{bmatrix}\_{n \times 3m}$$

The principle of least squares is used to obtain the normal equation corresponding to the above equation:

$$\mathbf{C}^{\mathsf{T}}P(L+\boldsymbol{\sigma}) = \mathbf{C}^{\mathsf{T}}P\mathbf{C}\mathbf{x} \tag{6}$$

where *P* refers to the weight matrix of the observed column vector, and it can be assigned as the identity matrix when all satellites and atomic clocks are stable according to the classical least square method [23]. The normal equation is solved to obtain the clock offset parameter *x* of satellite-ground and inter-satellite multi-source data of all satellites.

#### *2.3. Accuracy Evaluation Method*

The accuracy evaluation of satellite clock offset is effective to verify the reliability of satellite clock offset data processing results. The optimal clock offset accuracy evaluation uses the theoretical real clock offset to evaluate the calculated clock offset accuracy. However, it is very difficult to obtain the theoretically real clock offset. The clock offset accuracy can be evaluated through different methods such as the inner coincidence accuracy comparison method and the outer coincidence accuracy comparison method. As for the internal coincidence accuracy comparison method, the accuracy is judged based on internal errors that generally refer to the fitting residual or prediction error of clock offset. According to the external coincidence accuracy comparison method, the accuracy is judged based on external information or other link information, which generally refers to the closure error of clock offset, mainly including the inter-satellite station closure error and three-satellite closure error.

#### 2.3.1. Fitting Residual

The discrete points of clock offset are fitted to obtain the fitted clock offset, and the interval between discrete points is 1 s. The differences between the fitted clock offset and the actual clock offset is calculated to obtain the fitting residual of clock offset. The root mean square (RMS) is usually used to represent the fitting residual of clock offset, which can be obtained as follows:

$$
\sigma = \sqrt{\frac{\sum\_{i=1}^{n} \left(\Delta T\_i^{\mathbb{C}} - \Delta T\_i^{\mathbb{O}}\right)^2}{n-1}} \tag{7}
$$

where *n* represents the number of observation data, Δ*T<sup>C</sup> <sup>i</sup>* represents the calculated value of the *i*th fitted clock offset, and Δ*T<sup>O</sup> <sup>i</sup>* denotes the observed value of the *i*th actual clock offset. Generally speaking, the smaller RMS of the fitting residual indicates the better stability of clock offset, and the higher accuracy of time synchronization. it is noted that the actual clock offset is the post precision satellite clock offset [24].

#### 2.3.2. Prediction Error

The time scale is set to improve the accuracy of satellite broadcast clock offset parameters, so as to ultimately serve users well. The clock offset parameter is calculated by the whole-network adjustment method and used for clock offset prediction. The difference

between prediction result and actual clock offset is calculated to obtain the prediction error of clock offset.

$$\begin{aligned} \Delta T\_{\text{error}}(t) &= \Delta T\_{\text{pred}}(t) - \Delta T\_{\text{actual}}(t) \\ &= a\_0 + a\_1(t - t\_0) + a\_2(t - t\_0)^2 - \Delta T\_{\text{actual}}(t) \end{aligned} \tag{8}$$

where Δ*Tactual*(*t*) refers to the actual clock offset, Δ*Tpred*(*t*) denotes the predicted clock offset, and *a*0, *a*1, *a*<sup>2</sup> indicate the time difference, frequency difference and frequency drift, respectively, which are the clock offset parameters calculated through the whole-network adjustment.

The satellite clock offset prediction is divided into short-term prediction and long-term prediction. Specifically, first-order polynomial fitting is performed by 2-h clock offset data to predict 1-h result, and statistics the RMSs of 1-h prediction error. In contrast, secondorder polynomial fitting is performed by 48-h clock offset data to predict 24-h result, and statistics the RMSs of 24-h prediction error. This paper adopts the short-term prediction method for evaluation.

#### 2.3.3. Closure Error

In order to further verify the effectiveness of the whole-network adjustment method, time synchronization performance is tested according to the closure errors of clock offset. The commonly used closure errors of clock offset include the inter-satellite station closure error and three-satellite closure error. The inter-satellite station closure error is shown in Figure 3.

**Figure 3.** Schematic diagram of inter-satellite station closure error.

Supposing the satellite *i* is the node satellite, the clock offset of satellite *j* can be observed indirectly. If there exist the two-way clock offset Δ*Ti* between satellite *i* and ground station, and Δ*Tij* between satellite *i* and satellite *j* in the same time period, the clock offset of satellite *j* can be obtained as:

$$
\Delta T\_j' = \Delta T\_{i\bar{j}} + \Delta T\_{\bar{i}} \tag{9}
$$

The difference between Δ*T <sup>j</sup>* and satellite-ground two-way time synchronization observed value Δ*Tj* of satellite *j* is calculated to obtain the residual of the indirect clock offset. Since it is difficult to find the observed values of Δ*Ti*, Δ*Tij* and Δ*Tj* in the same time period, one day of data are acquired. Δ*Ti* and Δ*Tij* are fitted to obtain the indirect fitting parameter of Δ*T <sup>j</sup>* and then calculate the fitting residual of Δ*Tj* within the direct observation period.

### **3. Results**

#### *3.1. Data Situation*

In order to obtain comprehensive and sufficient calculation results, this paper calculates the whole-network adjustment results of multi-source data according to the L-band satellite-ground two-way time synchronization and Ka-band inter-satellite two-way time synchronization of BDS on 1 March 2022. In addition, this paper makes a comparative analysis of "one-hop" reduction results for the same period. It is worth noting that GEO satellites mainly deliver SBAS, RDSS and other services, so they are without involvement in inter-satellite calculations. Therefore, this paper selects the data of 24 MEO satellites and 3 IGSO satellites for calculation. The data of all the 27 satellites are used for the intersatellite calculation. The satellite-ground clock offset curve is shown in Figure 4, and ISL establishment during this period is shown in Figure 5.

**Figure 4.** Curve of SGL clock offset.

**Figure 5.** ISL between C38 and other satellites.

As shown in Figure 4, the satellite-ground clock offset curve changes steadily, and we can see from Figure 5 that satellite C38 and other 26 satellites all have established links. Based on such data, we verify the validity of the whole-network adjustment method.

#### *3.2. Fitting Residual*

The fitting residuals are calculated based on the clock offset time series. Taking satellites C35, C36 and C37 as example, the fitting residuals using "one-hop" reduction method and the whole-network adjustment method are shown in Figure 6.

**Figure 6.** Clock residuals.

As shown in Figure 6, the fitting residuals of "one-hop" reduction method are fluctuating between ±3.5 ns, and the RMSs of fitting residuals of satellites C35, C36 and C37 are 0.25 ns, 0.37 ns, and 0.38 ns, respectively. The fitting residuals of the whole-network adjustment method are fluctuating between ±1 ns, and the RMSs of fitting residuals of satellites C35, C36 and C37 are 0.10 ns, 0.14 ns, and 0.13 ns, respectively.

The RMSs of fitting residuals of all satellites are further calculated, which reflect the internal coincidence accuracy of fitting results. The RMS statistical results of clock offset fitting residuals of 27 satellites are detailed in Table 1.


**Table 1.** RMS statistical results of fitting residuals (unit: ns).


**Table 1.** *Cont.*

As shown in Table 1, the RMSs of clock offset fitting residuals of all satellites subject to "one-hop" reduction method between 0.19 ns and 0.42 ns, with a mean of 0.32 ns. The RMSs of clock offset fitting residuals of all satellites subject to whole-network adjustment range between 0.10 ns and 0.19 ns, with a mean of 0.13 ns. Compared with the "one-hop" reduction, whole-network adjustment improves the clock offset fitting residual accuracy by nearly 45.06%.

It is noted that the improvement rate differs widely, the reason is that the value of all satellites subject to "one-hop" reduction method differs widely, and the value subject to the whole-network adjustment method is relatively small and close. This result further proves that the whole network adjustment method eliminates the equipment errors of different satellites.

#### *3.3. Prediction Error*

In order to further analyze the effectiveness of whole-network adjustment, the 2-h clock offset data are used for slip prediction in the next one hour. The slip prediction result is compared with the actual observation value to obtain the prediction error. In order to better analyze the effectiveness of different methods, we also use only the satellite-ground data for comparison. Taking satellites C28, C29 and C30 as example, the clock offset prediction errors of only satellite-ground method, "one-hop" reduction method and the whole-network adjustment method are shown in Figure 7.

**Figure 7.** Prediction error.

As shown in Figure 7, the prediction error of only satellite-ground is fluctuating between ±6 ns, and the RMSs of clock offset prediction errors of satellites C28, C29 and C30 are 0.54 ns, 0.68 ns and 1.12 ns, respectively. The clock offset prediction error of "one-hop" reduction method are fluctuating between ±4 ns, and the RMSs of clock offset prediction errors of satellites C28, C29 and C30 are 0.46 ns, 0.50 ns, and 0.91 ns, respectively. The clock offset prediction error of the whole-network adjustment method are fluctuating between ±1 ns, and the RMSs of clock offset prediction errors of satellites C28, C29 and C30 are 0.15 ns, 0.26 ns, 0.23 ns, respectively. The RMSs of clock offset prediction errors of all satellites are further calculated and the statistical results are shown in Table 2.


**Table 2.** RMS statistical results of prediction error (unit: ns).

As shown in Table 2, the RMSs of predicted only satellite-ground method range between 0.33 ns and 1.66 ns, with a mean of 0.72 ns. The RMSs of predicted clock offset of "one-hop" reduction method range between 0.32 ns and 0.98 ns, with a mean of 0.54 ns. The RMSs of predicted clock offset of the whole-network adjustment method range between 0.15 ns and 0.43 ns, with a mean of 0.25 ns. Compared with the only satellite-ground method and "one-hop" reduction method, the prediction accuracy of the whole-network adjustment method improves about 62.13% and 52.15%, respectively.

#### *3.4. Closure Error*

Analyzing the data on 1 March 2022, it can be seen that there are 226 inter-satellite station closure errors and 1710 three-satellite closure errors concerning 27 satellites. Some inter-satellite station closure errors and three-satellite closure errors are detailed in Table 3, and statistical results of all closure errors are itemized in Table 4.


**Table 3.** RMS statistical results of some inter-satellite station and three-satellite closure errors (unit: ns).

**Table 4.** RMS statistical results of closure errors (unit: ns).


It can be seen from Tables 3 and 4, among 226 inter-satellite station closure errors and 1710 three-satellite closure errors, the closure error of inter-satellite station and threesatellite of "one-hop" reduction method are 0.69 ns and 0.23 ns, respectively. In addition, the closure error of the whole-network adjustment method is 1.34 × <sup>10</sup>−<sup>10</sup> and 5.54 × <sup>10</sup>−11, respectively, which both approximate 0. The results show that the closure error of inter-satellite station is subject to SGL reduction, while the closure error of threesatellite is subject to ISL measurement only. Therefore, three-satellite closure error is superior to the inter-satellite closure error, indirectly indicating that the SGL measurement error is a limiting factor that improves in the system time synchronization accuracy.

#### **4. Conclusions**

Based on the L-band satellite-ground clock offset data and the Ka-band inter-satellite clock offset data, this paper provides a method of whole-network adjustment for clock offset. In addition, the fitting residual, prediction error and closure error are used to evaluate and verify the effectiveness of this new method. This new method makes the RMS of the clock offset fitting residual reach up to 0.13 ns, which is lower than that of "one-hop" reduction by about 45.06%. As well, the RMS of clock offset prediction error hits about 0.25 ns, which is lower than those of only satellite-ground and "one-hop" reduction by about 62.13% and 52.15%, respectively. In addition, the closure error of inter-satellite station and three-satellite are approximate 0, which is much lower than 0.69 ns and 0.23 ns for "one-hop" reduction. This paper enriches the satellite navigation time synchronization system, and can be applied to clock offset measurement, systematic difference calculation and so forth. Furthermore, it provides a theoretical basis for effectively resolving systematic deviations in satellite-ground links and inter-satellite links. Most importantly, it will ultimately contribute to an improvement in the overall performance of BDS in navigation, positioning and timing.

**Author Contributions:** Conceptualization, D.W. and J.P.; methodology, L.L.; software, R.G. and X.L.; validation, D.W. and R.G.; formal analysis, R.G. and H.Y.; resources, D.W. and C.T.; data curation, D.W. and R.G.; writing—original draft preparation, R.G. and D.W.; writing—review and editing, D.W.; funding acquisition, R.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 41874043.

**Data Availability Statement:** The data sources are supported by Beijing Satellite Navigation Center.

**Acknowledgments:** Thanks for the data sources support of Beijing Satellite Navigation Center, and thanks for the fund support of the National Natural Science Foundation of China (No. 41874043).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

