An Improved Multi ‐ Satellite Method for Evaluating Real ‐ Time BDS Satellite Clock Offset Products

: Two methods are widely used for evaluating the precision of satellite clock products, namely the single ‐ satellite method (SSM) and the multi ‐ satellite method (MSM). In the satellite clock product evaluation, an important issue is how to eliminate the timescale difference. The SSM selects a reference satellite to eliminate the timescale difference by between ‐ satellite differencing, but its evaluation results are susceptible to the gross errors in the referenced satellite clock offsets. In the MSM, the timescale difference is first estimated and then removed. Unlike the GPS, the BeiDou Navigation Satellite System (BDS) consists of three types of satellites, namely geosynchronous earth orbit (GEO), inclined geosynchronous orbit (IGSO), and medium earth orbit (MEO) satellites. The three types of satellites have uneven orbital accuracy. In the generation of satellite clock products, the orbital errors are partly assimilated into the clock offsets. If neglecting the orbital accuracy difference of the three types of BeiDou satellites, the MSM will obtain biased estimates of the timescale difference and finally affect the clock product evaluation. In this study, an improved multi ‐ satellite method (IMSM) is proposed for evaluating the real ‐ time BDS clock products by removing the assimilated orbital errors of the three types of BDS satellites when estimating the timescale difference. Three real ‐ time BDS clock products disseminated by three different International GNSS Service (IGS) analysis centers, namely CLK16, CLK20, and CLK93, over a period of two months are used to validate this method. The results indicate that the assimilated orbital errors have a significant impact on the estimation of the timescale difference. Subsequently, the IMSM is compared with the SSM in which the referenced satellite is rigorously chosen, and their RMS difference is only 0.08 ns, which suggests that the evaluation results obtained by the IMSM are accurate. Compared with the traditional MSM, the IMSM improves the RMS by 0.16, 0.11, and 0.07 ns for CLK16, CLK20, and CLK93, respectively. Finally, three real ‐ time BDS clock products are evaluated using the proposed method, and results reveal a significant precision difference among them. We the effect of the assimilated orbital errors on the precision evaluation of clock products. The effect directly caused by the assimilated orbital errors is significantly different for the GEO, IGSO, and MEO satellites with an averaged RMS of 2.265, 0.407, and 0.139 m, respectively. Besides, the assimilated orbital errors can cause some discontinuities in the clock products as the orbit updates. The effect of the discontinuities is comparable for the three types of satellites, with a mean value of 0.18 ns. The evaluation results using the IMSM method are compared with those of the existing methods, i.e., the multi ‐ satellite method and the single ‐ satellite method. For the latter, the referenced satellite is carefully selected to ensure the correctness of the results obtained by the single ‐ satellite method. The results show that the RMS differences between the IMSM method and the multi ‐ satellite method exceed 0.13 ns, indicating an obvious effect of the assimilated orbital errors on the multi ‐ satellite method. In contrast, the results derived from the single ‐ satellite method and the IMSM method are more consistent, with an RMS difference of smaller than 0.08 ns for CLK16, CLK20, and CLK93. Compared with the traditional MSM, the IMSM improves the RMS by 0.16, 0.11, and 0.07 ns for CLK16, CLK20, and CLK93, respectively. These results suggest that the evaluation results obtained by the IMSM are accurate and confirm the correctness of the IMSM method. Moreover, the IMSM method can overcome the problem of the single ‐ satellite method that is susceptible to the gross errors in the referenced satellite clock offsets. BDS and CLK93 nanoseconds.


Introduction
On April 1, 2013, the International GNSS Service (IGS) Real-Time Pilot Project was officially launched to provide real-time service (RTS). The RTS products consist of GNSS satellite orbit and clock corrections to the broadcast ephemeris. Such corrections are formatted according to the Radio Technical Commission for Maritime Services (RTCM) state-space representation (SSR) standard and disseminated using the networked transport of RTCM via internet protocol (NTRIP). Based on the RTS products, the precise point positioning (PPP) technology can be processed in real-time mode [1][2][3], which has been applied to various real-time fields such as vehicle navigation [4], aerial triangulation [5], time transfer [6], and architecture deflection detection [7]. On December 27, 2018, the BeiDou Navigation Satellite System (BDS) started to provide global services. Currently, a few IGS analysis centers (ACs) and institutes provide free-accessed real-time precise BDS orbit and clock products, including Wuhan University (WHU), Centre National D 'etudes Spatiales (CNES), Deutsches Zentrum für Luft-und Raumfahrt (DLR), and Deutsches GeoForschungsZentrum (GFZ). Besides, some commercial PPP services, such as Trimble RTX, Fugro Starfix, and UniStrong Atlas/ChinaCM, also provide BDS real-time precise products. The availability of real-time BDS products enables users to perform real-time BDS PPP. For RTS providers, assessing the precision of real-time products is a routine job. For the users, the real-time BDS PPP performance relies on the quality of the precise products. The product precision is also a key indicator of satellite health status and an essential component of service integrity information. Therefore, it is vital to evaluate real-time BDS clock products.
In the precise clock estimation, the satellite and receiver clock offsets are linearly dependent. To make the clock offsets estimable, a referenced clock offset or the so-called timescale is usually introduced, which means that the estimated clock offsets are relative values [8]. Therefore, precision is a key quality indicator for the clock products. In the precision evaluation of clock products, it is necessary to eliminate the timescale difference, which is a time-varying system bias with respect to the referenced clock product [8,9]. Conventionally, there are two commonly used methods to evaluate the precision of GPS clock products, i.e., the single-satellite method (SSM) and multi-satellite method (MSM). The SSM selects a reference satellite to form satellite-differenced (SD) and productdifferenced (PD) clock offsets in which the timescale difference can be eliminated by the differencing operation [10,11]. However, in this method, the selection of the referenced satellite is vital since the gross errors in the clock offsets of the referenced satellite can be propagated into the SD and PD clock offsets. Unlike the SSM, the MSM estimates the timescale difference using the PD clock offsets of all satellites, and subsequently removes the timescale difference at every epoch [8,12]. Therefore, the evaluation results obtained by the MSM can be less affected by the clock gross errors included in the selected reference satellite [8].
Although many efforts have been made to assess the precision of IGS legacy or multi-GNSS experiment (MGEX) products [13][14][15][16][17][18][19], including real-time BDS clock products [11,12,20,21], few studies paid attention to the impact of orbital errors on the evaluation results. In the precise clock estimation, the orbital errors are partially assimilated in the clock offsets [22][23][24]. Due to the smaller precision difference between different GPS satellites, the effect of assimilated orbital errors on the precision evaluation of clock products is usually negligible. Dissimilar to the GPS, the BDS constellation consists of three types of orbits, namely geosynchronous earth orbit (GEO), inclined geosynchronous orbit (IGSO), and medium earth orbit (MEO). As the accuracy of GEO orbit products is much lower than that of IGSO and MEO satellites [25], the assimilated orbital errors for the GEO satellites are also much larger, which will bias the timescale difference estimates in the MSM since it is a common parameter for all three types of BDS satellites. In the process of the real-time clock offset estimation, the satellite coordinates are usually fixed using the predicted ultra-rapid orbit products [26][27][28]. When the ultra-rapid orbit updates, some discontinuities in the generated real-time orbit and clock products will occur, which means the clock discontinuity is essentially caused by the assimilated orbital errors. Similar to the assimilated orbital errors, the caused clock discontinuities also affect the estimation of the timescale difference. Based on the above analysis, an improved multisatellite method (IMSM) is proposed to evaluate the real-time BDS clock products by removing the assimilated orbital errors of the three types of BDS satellites when estimating the timescale difference. The idea of removing the assimilated orbital errors for the three types of BDS satellites in the IMSM is analogical to assign different weights for different types of BDS satellites in the estimation of the timescale difference.
The paper is outlined as follows. First, an extended model of real-time BDS clock products is depicted, and then the improved multi-satellite method is presented. Next, three kinds of available real-time precise BDS clock products, i.e., WHU CLK16, DLR CLK20, and CNES CLK93, are used to validate the IMSM. Afterward, three BDS clock products are evaluated using the IMSM. Finally, some conclusions are summarized.

Method
In the precision evaluation of real-time clock products, a more precise clock product is usually selected as a reference. Then, the timescale difference needs to be eliminated. In this section, an improved multi-satellite method is presented for evaluating the real-time BDS clock products by removing the assimilated orbital errors in the estimation of the timescale difference.

Model of Real-time BDS Clock Product
The generation methods of the real-time clock offset product can be classified into two categories, i.e., the un-differenced methods [29] and mixed-differenced methods [30]. For both types of methods, a timescale must be introduced to separate the satellite and receiver clock offsets due to their linear dependency. The introduced timescale may differ among different ACs due to their different data processing schemes, which will cause a timescale difference between different clock products. Besides, an initial clock is introduced to convert the epoch-differenced clock offsets into the un-differenced ones in the mixed-differenced method. Overall, the traditional model of the clock product can be described as [8,31] , where superscript and subscript denote the satellite and AC, respectively, is the clock product in seconds, is the timescale in seconds, is the initial clock bias in seconds, and is the phase estimation clock correction in seconds, which can be expressed as , (2) where is the true value of the clock correction, and is the noise. The timescale is a timevarying offset for all satellites, and thus it can be absorbed by the receiver clock offset. The initial clock bias is a constant offset on a continuous arc, which only affects the PPP convergence time [31]. The final PPP accuracy is dependent on the clock correction [8]. In summary, the PPP performance is only affected by and . The accuracy of BDS GEO orbit products is over five times worse than that of IGSO and MEO orbit products [21,25]. In the real-time clock estimation, the ultra-rapid orbit products are usually used to fix the satellite coordinates, and the orbital errors can be partially assimilated by the clock offsets [22][23][24]. Besides, the discontinuities of the orbital errors will also cause clock discontinuities due to the high correlation between the radial orbit and clock offset [20,21]. In other words, the clock discontinuity can be regarded as another effect caused by the assimilated orbital errors. Thus, the derived clock products for different BDS satellite types contain different magnitudes of orbital errors. The time-varying part of the assimilated orbital errors can be absorbed into the estimated timescale difference, resulting in deviated evaluation results for all satellites. Therefore, the traditional model described by Equation (1) is not suitable for BDS clock products. An extended model of real-time BDS clock products can be expressed as , where is the effect caused by the assimilated orbital error in the unit of second. In the precision evaluation, a post-processed clock product is usually employed as the referenced product. Due to its higher precision, the initial clock bias, assimilated orbital error, and noise in the referenced clock product can be neglected. Thus, the reference product provided at the analysis center can be expressed as .
By subtracting Equation (4) from Equation (3), the derived PD clock offset can be formulated as follows: , where • is the between-product differencing operator for the products provided by the analysis center and , is the PD clock offset, and is the timescale difference.

Effect of the Assimilated Orbital Errors
As stated in the previous section, the effect of the assimilated orbital errors includes not only the effect directly caused by the orbital errors but also the influence caused by the clock discontinuity. Therefore, the effect of the assimilated orbital errors in Equation (3) can be further expressed as where is the effect directly caused by the orbital errors in the unit of second, and , is the effect of the clock discontinuity at epoch in the unit of second.
Due to the correlation of the orbit and clock offset, about 96% of radial orbital errors can be absorbed by the clock offset in the estimation of satellite clock offsets [23]. The orbital error along the line-of-sight direction can be computed and corrected from the observations [23,24]. A practical and straightforward formula can be expressed as where is the orbital error on the line-of-sight direction in meters, , , is the unit lineof-sight vector from receiver to satellite , and Δ , , is the PD coordinate vector for satellite in the Earth-Centered Earth-Fixed system in meters. Since post-processed orbit products can provide satellite coordinates with negligible errors, real-time orbit products can be compared with the post-processed orbit products to obtain the PD coordinate vector. Thus, Equation (7) can be used to estimate the radial orbital errors. As the satellite clock offsets are estimated in a GNSS reference network, the orbital error can be approximately described as where is the elevation angle-dependent weight, is the number of stations in the reference network, and is the speed of light in vacuum in meters per second.
The orbit update can easily result in clock discontinuity, which is another effect caused by the assimilated orbital errors, and also should be considered in the precision evaluation. After correcting the orbital errors and timescale difference , the corrected product-differenced and epochdifferenced (PDED) clock offset can be used to detect and correct the clock discontinuities in an iterative manner. The initial value of the timescale difference can be calculated based on the PD clock offsets of MEO satellites using the MSM method. The criteria for identifying the discontinuities are provided as follows: where ∇ • , is the between-epoch differencing operator for epoch and 1, ∇Δ , , and ∇Δ , , are the corrected PDED clock offsets, and ∇ and ∇ are the mean and standard deviation (STD) of the corrected PDED clock offsets, respectively. If the conditions in the inequations (9) are satisfied, a discontinuity is considered to occur at epoch as Once the discontinuity is determined, it can be removed from the clock offset series since epoch .
After correcting the effect of the assimilated orbital errors, the corrected PD clock offset Δ can be expressed as Equation (11) following Equation (5).
where is the noise of the corrected PD clock offset Δ , which includes the noise of , i.e., , and the noise introduced by . In the proposed method, the corrected PD clock offset series obtained from Equation (11) are used to estimate the timescale difference. In contrast, if the effect of the assimilated orbital errors is not considered, as done in the multi-satellite method, the timescale difference is estimated as [8] , 1 , where , is the timescale difference estimate obtained by the multi-satellite method, and is the number of satellites. After placing Equation (5) into Equation (12) under the zero-mean assumption ∑ 0, we can obtain where ∑ is the noise, and ∑ is the bias caused by the assimilated orbital errors. Unfortunately, the bias is not a constant, meaning that it will affect the calculation of the STD for the clock offsets. Therefore, the effect of the assimilated orbital errors must be carefully considered when assessing the precision of BDS real-time clock products.

Precision Evaluation of BDS Real-time Clock Products
Before evaluating the clock products, the timescale difference needs to be estimated and removed. In the IMSM, the assimilated orbital errors are corrected from the clock offsets when estimating the timescale difference. Since the assimilated orbital errors for the BDS GEO satellites are larger than the IGSO and MEO satellites, the corrections for the GEO satellites are also larger. This strategy is similar to assigning different weights for different types of BDS satellites in the estimation of the timescale difference.
To estimate the timescale difference, a robust least-squares method is employed. The error equation is expressed as where satellite 1,2, … , and epoch =1,2,…,n, Δ , Δ , , … , Δ , is the -by-1 vector of the PD clock offsets after correcting the effect of assimilated orbital errors, is the residual vector, is the -by-identity matrix, is the -by-matrix whose elements at the column are one while the others are zero, , … , is the -by-1 vector of the estimated initial clock bias for the satellites, Δ , , … , , is the -by-1 vector of the estimated timescale difference for the epochs, and is the ( * )-by-( * ) weight matrix, which can be determined according to a robust weighting function [32] with an assumption of no correlation between satellites or epochs. Since Equation (14) is under-determined, a zero-mean constraint is added to . Once the timescale difference , is determined by Equation (14), the estimated timescale difference can be removed from the PD clock offsets in Equation (5) as follows.
where Δ , is the PD clock offset after subtracting the timescale difference, and it is also called the double-differenced (DD) clock offset for simplicity. The items on the right side of Equation (15) are the errors to affect PPP performance [8,24,31]. Therefore, the precision of the clock products can be evaluated using the STD of the DD clock offset as where and are STD and mean values of the DD clock offset series, respectively, and is the element number in the clock offset series. For instance, for a 24-h dataset at a sampling interval of 30 seconds, the value of is 2880. The entire process for implementing the IMSM is depicted in Figure  1.

Results
To verify the IMSM method, experiments were conducted based on real-time BDS clock products provided by the IGS RTS. Currently, these products are generated based on the RTS global network consisting of 256 stations (http://www.igs.org/network). After registering online, the users can access these real-time products from the RTS NTRIP casters. For the caster located in the GNSS Research Center of WHU, available BDS products include WHU CLK15/16, DLR CLK20/21, and CNES CLK92/93. Among them, the orbit products of CLK16, CLK20, and CLK93 refer to the satellite antenna phase center (APC) and the others refer to the center of mass (CoM), which is the only difference for these product couples. Table 1 presents the data details and processing strategies used in this study. The products CLK16, CLK20, and CLK93 were collected using the software BNC (BKG NTRIP caster) in real-time on the day of year (DOY) 213-243 and 255-284, 2019. The GFZ MGEX post-processed clock products provide clock offsets for BDS GEO, IGSO, and MEO satellites at a sample interval of 30 seconds. Therefore, the GFZ MGEX clock products were used as the reference products in the following experiments. The sample intervals of the real-time clock products are 5 seconds. In order to facilitate comparison and avoid interpolation errors, the real-time clock products were resampled to 30 seconds.

Effect Analysis of Assimilated Orbital Errors
This section analyzes the effect of the assimilated orbital errors on the precision evaluation of the clock products. The assimilated orbital errors can be computed according to Equations (6)  The computed assimilated orbital errors are presented in Figure 3 with a representative case of CLK16 on DOY 256, 2019. As shown in Figure 3(a), the effect directly caused by the assimilated orbital errors is significantly different for the GEO, IGSO, and MEO satellites with an averaged RMS of 2.265, 0.407, and 0.139 m, respectively. The assimilated orbital errors are significantly larger for the GEO satellites. This is reasonable because the GEO satellites have the poorest orbit accuracy due to their geostationary characteristics. In addition, the clock discontinuity is another effect caused by the assimilated orbital errors, and the detected discontinuities are shown in Figure 3(b). As seen in this subfigure, some discontinuities occur approximately every three hours when the orbits are updated. Figure 3(b) also shows that the discontinuities have comparable magnitudes for GEO, IGSO, and MEO satellites with a mean value of 0.18 ns.    To analyze the effect of the assimilated orbital errors on the estimation of the timescale difference, the timescale differences are estimated using the MSM and IMSM. Then, they are compared. Figure 4 presents the timescale differences estimated by both methods and their discrepancies for CLK16 on DOY 256, 2019. In Figure 4(a), the timescale difference series derived from both methods show similar trends. After examining the discrepancies of the two series, a varying trend can be found, and some discontinuities occur at about 6 and 12 h, as shown in the blue circles of Figure 4(b). The former reflects a systematic effect of the assimilated orbital errors on the timescale difference estimates, and the latter is caused by the discontinuities.  Furthermore, we analyzed the impact of the assimilated orbital errors on the precision evaluation of BDS clock products. Figure 5 presents the PD radial orbits and DD clock offsets derived from the MSM and the IMSM for a few typical satellites. For the MSM, the DD GEO clock offsets are highly correlated with the PD radial orbits, as shown in the left panel. However, the MEO and IGSO clock offsets do not exhibit a strong correlation with their orbit errors but show an opposite variation to the GEO clock offset errors. This is because part of the assimilated orbital errors in the GEO clock offsets has been propagated into the timescale difference parameter. Consequently, these errors are further propagated into the DD clock offsets of the IGSO and MEO satellites because their DD clock offsets are derived by subtracting the common timescale difference from the PD clock offsets. For the IMSM, the correlation of the clock offsets between the GEO and IGSO/MEO satellites is largely reduced. This is because the effect of the assimilated orbital errors, especially from the GEO satellites, is removed when estimating the timescale difference. As a result, the evaluation of the BDS IGSO and MEO clock offsets will not be subject to the substantial effect of the assimilated orbital errors from the GEO satellites.

Comparison of Evaluation Results
In order to validate the IMSM method, its evaluation results are compared with those derived from the SSM and MSM. Since the SSM eliminates the timescale difference by between-satellite differencing operation, the timescale difference does not need to be estimated. Nevertheless, this method depends on the data quality of the selected reference satellite clock. In addition, the real-time clock offsets are often influenced by gross errors and data stream interruptions, which make the SSM performance unstable and unreliable [8]. Therefore, the referenced satellite clock is carefully selected to ensure the correctness of the SSM evaluation results by abiding by the following two conditions. First, the data availability rate is higher than the threshold. Second, daily STD is within the threshold limit. The thresholds for the data availability rate and daily STD are empirically set to 85%-95% and 0.2-0.3 ns, respectively. It is noted that the MEO orbits and clock products generally have the highest accuracy among the three kinds of BDS satellites. Therefore, the MEO satellites are chosen as referenced satellites since they are better at satisfying the above conditions. After selecting the referenced satellite cautiously, the SSM is used to validate the IMSM method.
The DD clock offsets are computed using the SSM, MSM, and the IMSM. Figure 6 shows the DD clock offsets of a few typical satellites for CLK16 on DOY 256, 2019. In the implementation of the SSM, MEO C14 is selected as the referenced satellite following the above two principles. In the MSM, a few discontinuities are found in the blue circles, which agree with the discontinuities in the estimated timescale differences shown in Figure 4. By contrast, the DD clock offsets obtained by the IMSM do not show noticeable discontinuities since they have been removed when estimating the timescale difference. Moreover, the DD clock offsets derived from the IMSM are similar to those obtained from the SSM. In this figure, although the improvement for the DD clock offsets of the GEO satellites is not obvious due to a larger axis range, the IMSM-derived DD clock offset series are more consistent with those derived from the SSM, as seen in Table 2. The STDs of the DD clock offsets are computed for all satellites and listed in Table 2. It appears that the STDs derived from the IMSM show good agreement with those of the SSM. However, the STDs between the MSM and the IMSM have significant differences for all three types of satellites, indicating a significant effect of the assimilated orbital errors.  We further evaluate the IMSM by comparing the results with the SSM and MSM using the data from DOY 213 to 243 and from 255 to 284, 2019. In the SSM, the referenced satellites are selected IGSO MEO carefully following the aforementioned principles. The data are excluded for the days on which no referenced satellite satisfies the required principles. In order to have sufficient samples and obtain accurate evaluation results in the SSM simultaneously, the threshold for the data availability rate and STD limit is empirically set to 90% and 0.2 ns, respectively. Due to the relatively lower precision of CLK16 and CLK20, the STD threshold is set to 0.3 ns to ensure adequate samples. As a result, the number of days used for the comparison is 9, 11, and 31 for CLK16, CLK20, and CLK93, respectively. In the SSM, the selected reference clock offsets have an averaged STD of 0.18, 0.26, and 0.10 ns for CLK16, CLK20, and CLK93, respectively. By differencing the evaluation results of all satellites derived from the IMSM and an existing method, i.e., the SSM or MSM, Figure 7 shows the result differences of the IMSM with respect to the MSM and SSM as well as their RMS statistical values. The differences are displayed in a form of boxplots to comprehensively describe the error distribution of all satellites by minimum/maximum values, lower/upper quartiles, and median. As seen in Figure 7(a), the result differences between the IMSM and the MSM exceed 0.4 ns on some days for CLK16, CLK20, and CLK93, and their RMS differences are between 0.13 and 0.23 ns. In contrast, Figure 7(b) shows that the evaluation results derived from the SSM and the IMSM are more consistent, with RMS differences smaller than 0.08 ns for the three products. Compared with the traditional MSM, the IMSM improves the RMS by 0.16, 0.11, and 0.07 ns for CLK16, CLK20, and CLK93, respectively. These results indicate that the evaluation results obtained by the IMSM are accurate and that the IMSM is feasible to assess realtime BDS clock products. Moreover, the IMSM can avoid the selection of the reference satellites.    Furthermore, the result differences are classified according to their satellite types, and the averaged result differences are shown in Figure 8 for the GEO, IGSO, and MEO satellites. The right panel illustrates that the result differences between the SSM and the IMSM are small and within 0.06 ns for the three types of satellites. As the left panel indicates, the evaluation results of the IMSM tend to be larger and smaller for the GEO and MEO satellites, respectively, compared with the MSM. This can be attributed to the fact that the assimilated orbital errors of the GEO satellites have been removed before the timescale difference estimation in the IMSM. Besides, as seen in the left figure, the result differences exceed 0.16 ns for the GEO and MEO satellites. Since the nominal precision of the realtime precise BDS clock product is about 0.3 ns [20], such a difference is considered to be significant. The larger differences indicate that the MSM is subject to an effect from the assimilated orbital errors.

Evaluation of Three Real-time BDS Satellite Clock Products
In this section, the IMSM method is used to evaluate the real-time BDS clock products of WHU CLK16, DLR CLK20, and CNES CLK93. Figure 9 presents the daily averaged STDs of the GEO, IGSO, and MEO DD clock offsets for the two periods, namely DOY 213-243 and 255-284, 2019. Statistics of CLK16 are not calculated on some days due to its low data availability rate. In addition, the DD GEO clock offsets of CLK20 are not shown due to the data absence. As seen from these figures, the MEO and GEO clock products exhibit the best and worst precision, respectively. In terms of daily stability, the IGSO and MEO satellites also mostly outperform the GEO satellites. The precision stability can be affected by various factors such as data stream interruptions, quality control in real-time clock estimation, and an unhealthy state of satellites. Besides, the degraded precision for the CLK20 IGSO clock products in the second session could be caused by the data stream instability. Table 3 lists the precision of the three real-time BDS clock products for each satellite, which can be reflected by the averaged STDs of the DD clock offsets. It is seen that the precision difference among different products can reach 1.0 ns, while the precision difference for different satellite types is smaller than 1.0 ns.

Discussion
Based on the IMSM, we can objectively evaluate the precision of BDS real-time clock products. In the precision assessment, an important issue is how to eliminate the timescale difference. The existing single-satellite method removes the timescale difference by between-satellite differencing, while the existing multi-satellite method and the proposed IMSM method firstly estimate the timescale difference and then correct it. However, the single-satellite method needs to select a referenced satellite, which makes it inconvenient for automatic processing and may pose a risk that the evaluation results can be easily affected by the gross errors in the referenced clock offset. To overcome these problems, the multi-satellite method puts forward to use the multi-satellite reference. The main distinction between the proposed IMSM and the multi-satellite method is that when estimating the timescale difference, the effect of the assimilated orbital errors is corrected. The idea is similar to assigning different weights for three types of BDS satellites in the estimation of the timescale difference. As a result, the estimated timescale difference is more accurate.
The principles of the IMSM are mainly based on the fact that the orbital error can be absorbed by the clock offset to some extent [22][23][24], and the effect of the assimilated orbital errors cannot be ignored for BDS GEO and IGSO satellites, as shown in Figure 3. In this study, only BDS-2 real-time clock products are used to validate the proposed IMSM method. On 31 July 2020, the BDS-3 was announced formally commissioned. The BDS-3 constellation is also composed of GEO, IGSO, and MEO satellites. Therefore, the models and principles proposed in this study are also applicable to the BDS-3. As more ground stations begin to support BDS-3 signals, more BDS-3 real-time precise products from different ACs will become available and their quality will be further refined, which will facilitate the subsequent investigation and application of the IMSM method for evaluating BDS-3 real-time clock products.

Conclusions
An improved multi-satellite (IMSM) approach is proposed for evaluating real-time BDS clock products. In this approach, the effect of assimilated orbital errors is corrected when estimating the timescale difference. To validate this IMSM method, the real-time BDS clock products WHU CLK16, DLR CLK20, and CNES CLK93 were collected in real-time using the software BNC at two periods, i.e., DOY 213-243 and 255-284, 2019. We first analyzed the effect of the assimilated orbital errors on the precision evaluation of clock products. The effect directly caused by the assimilated orbital errors is significantly different for the GEO, IGSO, and MEO satellites with an averaged RMS of 2.265, 0.407, and 0.139 m, respectively. Besides, the assimilated orbital errors can cause some discontinuities in the clock products as the orbit updates. The effect of the discontinuities is comparable for the three types of satellites, with a mean value of 0.18 ns. The evaluation results using the IMSM method are compared with those of the existing methods, i.e., the multi-satellite method and the single-satellite method. For the latter, the referenced satellite is carefully selected to ensure the correctness of the results obtained by the single-satellite method. The results show that the RMS differences between the IMSM method and the multi-satellite method exceed 0.13 ns, indicating an obvious effect of the assimilated orbital errors on the multi-satellite method. In contrast, the results derived from the single-satellite method and the IMSM method are more consistent, with an RMS difference of smaller than 0.08 ns for CLK16, CLK20, and CLK93. Compared with the traditional MSM, the IMSM improves the RMS by 0.16, 0.11, and 0.07 ns for CLK16, CLK20, and CLK93, respectively. These results suggest that the evaluation results obtained by the IMSM are accurate and confirm the correctness of the IMSM method. Moreover, the IMSM method can overcome the problem of the single-satellite method that is susceptible to the gross errors in the referenced satellite clock offsets.
Based on the two-month dataset, the real-time BDS clock products CLK16, CLK20, and CLK93 are evaluated using the IMSM method. The results show that the MEO clock products exhibit the best precision among the three kinds of BDS satellites. By contrast, the GEO clock products have the poorest precision. The precision difference among different products can reach 1.0 ns, while the precision differs from different types of satellites with a few sub-nanoseconds.