Broadcast Ephemeris with Centimetric Accuracy: Test Results for GPS, Galileo, Beidou and Glonass

: Here we test the capability of the Broadcast Ephemeris Message, in both its GPS-like (Keplerian ellipse with secular and periodic perturbations) and Glonass-like (numerical integration of a 9D state vector) formats, to reproduce a corresponding precise ephemeris. We start from a daily Rinex 3.04 navigation ﬁle for multiple GNSS and the corresponding SP3 precise orbits computed by CNES (Centre National d’Etudes Spatiales) for GPS, Glonass, Galileo and CODE (Center for Orbit Determination in Europe) for Beidou, and compute broadcast ECEF coordinates and clocks. The pre-ﬁt discrepancies are converted by least squares to corrections to the broadcast ephemeris parameters in two-hour consecutive arcs (for GPS, Galileo and Beidou) and to a set of seven Helmert parameters for the entire day, to align in origin, orientation and scale to the common GNSS IGS14 Reference Frame. The test cases suggest that the Broadcast Ephemeris Message, complemented with Reference Frame information, can reproduce the precise ephemeris and clocks with centimetric accuracy for intervals at least equal to the respective validity times, typically 2 h. The broadcast ephemeris of Glonass consists of three initial positions and velocities at epoch, three constant Lunisolar accelerations for the satellite position, and of three polynomial coefﬁcients for the satellite clock. The 9D vector of state is numerically integrated to generate position and velocity data within the validity time (0.5 h) of the message. To test the capability of this model to reproduce the corresponding values of a precise ephemeris, the 9D vector of state and clock polynomials are adjusted until the rms (root mean squared spread) of the post-ﬁt residuals relative to a precise orbit (CNES’s in our case) is minimum. We show in one example (one satellite for one day) that the Glonass type of message can reproduce a precise ephemeris and clock with a rms spread of 0.025 m over one-hour arcs. Volume computations on one month of data with all available satellites conﬁrm the test results. For GPS, Glonass, Galileo and Beidou, the best ﬁtting clock values predicted by our second order polynomials, based on a 15 min sampling, are shown to ﬁt the corresponding high rate clocks (30 s sampling) of MGEX with zero bias and a rms spread of 0.062 ns (GPS G01), 0.023 ns (Galileo E01), 0.43 ns (Glonass R01), 0.086 ns (Beidou C07) and 0.086 ns (Beidou C12). Modiﬁcations to the GPS-like message structure and Glonass algorithm are proposed to increase the validity time by including the effect of the 3rd zonal harmonic of the Earth’s gravity ﬁeld. The potential of the RTCM messages for broadcasting the improved navigation message is reviewed.


Introduction
Is the limited accuracy of the broadcast ephemeris caused by inherent limitations of the model to compute the spacecraft coordinates and clock offset, or by the model coefficients in the broadcast message being insufficiently accurate? In this paper we provide arguments based first on detailed test cases and then on volume computations that the broadcast ephemeris can in fact provide spacecraft coordinates and clock offsets as accurate as the corresponding precise ephemeris provided that the coefficients of the broadcast ephemeris are appropriately tuned. addresses the improvement of the Glonass navigation message using SP3 orbits and clocks. Then, specific issues are addressed: comparison between high rate clock products of IGS/MGEX and the polynomial model (Section 4), and possible improvements in the message content for increased accuracy (Section 5). In Section 6, we apply the approach described above to volume computations involving one month of data and all the satellites of the four constellations. We show that the seven Helmert parameters averaged over the 30 days and all the satellites of each constellations average to zero, within one standard deviation, with the only exceptions of the rotations Rx and Rz of Galileo slightly above the one sigma threshold. This suggests that, on average, the GNSS-specific reference frame is aligned to the ITRF2014, but that for centimetric accuracy the seven Helmert parameters appropriate to the day and satellite should be used. Then we compare spectra of the pre-fit and post-fit residuals of the Broadcast-SP3 position differences averaged over one month and all satellites of a given constellation. We show that the former has significant power (1 m typically) at low frequencies (<1 cycle/6 h), whereas the latter has a spectrum very close to flat. This supports the idea that the proposed approach has effectively accounted for most of the signal, at the centimetric level. We remark in the Conclusions that if the broadcast ephemeris, when appropriately tuned, is capable of centimetric accuracy, then there are some important consequences. One is that the broadcast format, which is optimized in bandwidth requirements, could be used in real time for precision applications, taking advantage of the existing RTCM messages. Another inference is that ephemeris and clocks in the broadcast format could be used in place of the same data in the SP3 format, at least for MEO satellites.

Mathematical Model and Results for GPS, Galileo, BeiDou
To compute the ECEF coordinates of an SV as a function of time, the GPS navigation message adopts an algorithm based on a Keplerian ellipse with secular and periodic perturbations on selected orbital parameters. The validity time of the model is nominally two hours. Outside the validity time interval, the accuracy of the coordinates tends to degrade progressively. The orbital arc of nominal duration two hours has an orientation defined by three angles: Inclination I 0 , Longitude of Ascending Node Ω 0 and Longitude of Perigee ω. The size of the ellipse is given by the semimajor axis a and the eccentricity e. The angular position of the SV at a Reference Time Toe (Time of Ephemeris) is M 0 and is counted from the perigee (for a pictorial view and definitions see https://www.gsc-europa. eu/system-service-status/orbital-and-technical-parameters, accessed on 10 October 2021).
These six orbital parameters, which would be all constant if the total force were that of a point mass, are allowed to vary with time to account for deviations from the point mass force model. Inclination and Ascending Node are allowed to increase with a rate dI 0 /dt and dΩ 0 /dt constant during the validity interval; the orbital angular velocity n, which by Kepler's third law is inversely proportional to the semimajor axis to the 3/2 power, is perturbed by a constant term ∆n. The radial, along track and cross track position of the satellite are assigned periodic perturbations with a period equal to one-half the orbital period. These perturbations are described by sinusoids with separate amplitudes for the sine and cosine components. There is a total of 15 parameters which refer to the Time of Ephemeris Toe. Three additional clock parameters (offset a0, drift a1 and drift rate a2) refer to the Time of Clock Toc. Hence, every two-hour arc is described by 18 parameters and two reference epochs, Toe and Toc, which often coincide.
The Reference Precise Ephemeris are computed with rigorous constraints on the origin, orientation and scale of the Terrestrial Reference Frame (TRF), which in our case is ITRF2014 or, specifically, IGS14. In the GPS-like algorithm, only the shape and size of the orbital arc are defined. It follows that to generate spatial coordinates as close as possible to those in the SP3 file, it is appropriate to solve for a set of Helmert parameters. This can be done once per day, i.e., about two complete revolutions, because their global character requires sampling over the entire orbital arc. The mathematical model used for the GNSSs with a GPS-like message to express precise ephemeris and clocks in a broadcast form is described by Equation (1): [∆(XYZT)] 2 = f (Tx, Ty, Tz, Rx, Ry, Rz, Sc; a0, a1, a2, Crs, ∆n, M 0 , Cuc, e, Cus, sqrt(a), Cic, Ω 0 , Cis, I 0 , Crc, ω, dΩ dt , dI dt ) = min (1) We seek to minimize the sum, over the 96 epochs (384 equations) in one daily SP3 file, of the squared differences ∆(XYZT) of XYZT(SP3) and XYZT(BRDM), that is, the difference between the SP3 and broadcast ECEF coordinates and clock correction of a given Space Vehicle. This difference is parametrized by one set of seven Helmert parameters (three translations T, three rotations R and one scale Sc) valid for the entire day and 12 sets each of 18 parameters, one for each two-hour arc. The function f can be linearized by series expansion in a neighborhood of the broadcast values of the parameters (the Helmert parameters have zero a priori value). A partial derivative matrix H is computed numerically by approximating the partial derivatives with finite central differences. The best fitting corrections to the a priori values are computed by solving the linear system of 223 normal equations. The matrix H T H ( T stands for transposed), however, is poorly conditioned, as it may be expected. The rotational Helmert parameters tend to be highly correlated with the angles of orientation of the orbit. A well-conditioned matrix is obtained by increasing the weight of the matrix elements corresponding the parameters in red in Equation (1), which are therefore fixed to the broadcast values. In practice we solve for seven global Helmert parameters (i.e., one set for one day) and 12 sets each of seven parameters (the mean anomaly M 0 at the reference epoch Toe and three pairs of amplitudes of cosine and sine components of periodic perturbations along track (Cuc, Cus), radial (Crc, Crs) and across track (Cic, Cis) ( Figure 1). This is equivalent to constrain the Helmert rotations and solve for the angles defining the spatial orientation of the orbit.
The mathematical model used for the GNSSs with a GPS-like message to express precise ephemeris and clocks in a broadcast form is described by Equation (1) (1) We seek to minimize the sum, over the 96 epochs (384 equations) in one daily SP3 file, of the squared differences ∆(XYZT) of XYZT(SP3) and XYZT(BRDM), that is, the difference between the SP3 and broadcast ECEF coordinates and clock correction of a given Space Vehicle. This difference is parametrized by one set of seven Helmert parameters (three translations T, three rotations R and one scale Sc) valid for the entire day and 12 sets each of 18 parameters, one for each two-hour arc. The function f can be linearized by series expansion in a neighborhood of the broadcast values of the parameters (the Helmert parameters have zero a priori value). A partial derivative matrix H is computed numerically by approximating the partial derivatives with finite central differences. The best fitting corrections to the a priori values are computed by solving the linear system of 223 normal equations. The matrix H T H ( T stands for transposed), however, is poorly conditioned, as it may be expected. The rotational Helmert parameters tend to be highly correlated with the angles of orientation of the orbit. A well-conditioned matrix is obtained by increasing the weight of the matrix elements corresponding the parameters in red in Equation (1), which are therefore fixed to the broadcast values. In practice we solve for seven global Helmert parameters (i.e., one set for one day) and 12 sets each of seven parameters (the mean anomaly at the reference epoch Toe and three pairs of amplitudes of cosine and sine components of periodic perturbations along track ( , ), radial ( , ) and across track ( , ) (Figure 1). This is equivalent to constrain the Helmert rotations and solve for the angles defining the spatial orientation of the orbit. . schematic structure of the partial derivative matrix for the GPS-like model. The first vertical block represents the partials of XYZT pre-fit residuals relative to the Helmert parameters, appropriately scaled. The block diagonal part contains 12 blocks each of 4 rows of 18 partial derivatives, using as a priori the broadcast values of the parameters which are indexed with Toe for the coordinates and Toc for the clock. Schematic structure of the partial derivative matrix for the GPS-like model. The first vertical block represents the partials of XYZT pre-fit residuals relative to the Helmert parameters, appropriately scaled. The block diagonal part contains 12 blocks each of 4 rows of 18 partial derivatives, using as a priori the broadcast values of the parameters which are indexed with Toe for the coordinates and Toc for the clock.
In the following, we use this model for data referring to 2 January 2020, as an arbitrary choice. We consider G01 for GPS, E01 for Galileo, C07 for Beidou IGSO and C12 for Beidou MEO. The reference SP3 file is available at the IGS/MGEX servers. We have arbitrarily chosen the CNES SP3 file for GPS and Galileo (https://cddis.nasa.gov/archive/gnss/ products/mgex/2086/GRG0MGXFIN_20200020000_01D_15M_ORB.SP3.gz, accessed on 10 October 2021), and the CODE SP3 file for Beidou (https://cddis.nasa.gov/archive/gnss/ products/mgex/2086/COD0MGXFIN_20200020000_01D_05M_ORB.SP3.gz, accessed on 19 September 2021). The Mixed broadcast ephemeris file was downloaded from the BKG server (https://igs.bkg.bund.de/root_ftp/MGEX/BRDC_v3/2020/002/brdm0020.20p.Z, accessed on 10 October 2021). As mentioned in the Introduction, our results will be valid for the specific satellites and epochs only. To extend the validity of the conclusions drawn in this paper beyond the specific satellites and epoch requires volume calculations. In this paper we concentrate on the mathematical approach as a prerequisite for further massive processing. Epoch-dependent departures from the results described here could for example be caused by a different angle between the Sun's geocentric vector and the satellite orbital plane (beta angle), which controls the force caused by the solar radiation pressure. Figure 2 shows the pre-fit and post-fit residuals for G01 (Block IIF). The pre-fit residuals are computed using the Broadcast Ephemeris blocks of appropriate Toe and Toc, and the post-fit residuals are computed with the tuned parameters computed with Equation (1). In Figure 2, the pre-fit differences are represented with colored dots, one color for each twohour ephemeris block, with the exception of the first block where the broadcast ephemeris with Toe = 2 h has been used for an interval of 4 h, from 00:00 to 04:00, for test purposes. The vertical lines have a two-hour spacing. The diamonds represent the post-fit residuals, after the adjustment of the seven Helmert parameters and of the 11 sets of seven arc parameters. In Section 5, we will discuss in more detail the oscillatory pattern of the post-fit residuals, in connection to an improvement of the model to include the contribution of the third zonal harmonic of the gravity field. Likewise, it will be pointed out that the higher noise in the first 4 h of the post-fit residuals of the clock corrections suggests that the second order clock polynomial is unable to track the high frequency noise of the satellite clock, if the fit interval is as long as 4 h, as it may be expected. This will be further discussed in Section 4 in relation to the comparison with high rate clocks.

Results for GPS
The overall improvement in doing the fit described by Equation (1) is summarized as follows: the mean and rms (root mean square spread) of the pre-fit residuals is −0.013 ± 0.917 m, and 0.000 ± 0.051 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The adjusted clock polynomials fit the SP3 clock data with 0.007 m rms, whereas for the spatial coordinates the corresponding figure is a factor of seven larger. Figure 3 describes the epoch-wise corrections to the 11 sets of seven arc parameters, and the corrections to two of the three estimated clock parameters (a2 is not shown but computed). The correction to the mean anomaly M 0 means that at the epoch Toe the satellite is moved back or forth along the orbit, relative to the orbital phase predicted by the broadcast ephemeris. Table 1 finally provides the seven Helmert parameters for this daily set, and their estimated formal error 1 sigma.  The overall improvement in doing the fit described by Equation (1) is summarized as follows: the mean and rms (root mean square spread) of the pre-fit residuals is −0.013 ± 0.917 m, and 0.000 ± 0.051 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The adjusted clock polynomials fit the SP3 clock data with 0.007 m rms, whereas for the spatial coordinates the corresponding figure is a factor of seven larger.

Results for Galileo
For Galileo we have chosen the E01, an FOC (Full Operational Capability) SV, as an example. The repetition rate of the broadcast ephemeris is for Galileo usually 10 min, but not always at a regular rate. The rate of update of the message is considerably higher than GPS, so that a tailored ephemeris/clock is available at each epoch of observation. We seek broadcast model coefficients which make this high repletion rate unnecessary both for position and clocks. We expect that the validity time of the Galileo navigation message should be in fact comparable with that of GPS, for example. We have chosen those ephemeris blocks with a Toe/Toc closest to the even hours of day 2 January 2020. Figure 4 shows the difference between the raw residuals (original broadcast coordinates and clock minus corresponding SP3 values) described by colored dots (one color for each Toe); the post-fit residuals are shown by diamonds. The clock correction is expressed in meters. By comparison with GPS, it can be noted that the original ephemeris has, relative to the reference precise orbit, a pattern suggesting periodic resets, to keep the broadcast position and time as close as possible to the precise value, in the neighborhood of the respective Toe. In this way the effects of the third zonal harmonic (J3) are probably less visible than with a longer sampling time, as discussed for G01 in Figure 2. Moving away from this reference epoch, the departure from the reference values increases, so that the ephemeris needs to be replaced by a more recent block near the edge of the validity period. The broadcast ephemeris with best-fitting parameters (post-fit residuals, y axis to the right) tracks instead the precise ephemeris very closely and with no divergence at the edges of the validity interval.

Results for Galileo
For Galileo we have chosen the E01, an FOC (Full Operational Capability) SV, as an example. The repetition rate of the broadcast ephemeris is for Galileo usually 10 min, but not always at a regular rate. The rate of update of the message is considerably higher than GPS, so that a tailored ephemeris/clock is available at each epoch of observation. We seek broadcast model coefficients which make this high repletion rate unnecessary both for position and clocks. We expect that the validity time of the Galileo navigation message should be in fact comparable with that of GPS, for example. We have chosen those ephemeris blocks with a Toe/Toc closest to the even hours of day 2 January 2020. Figure 4 shows the difference between the raw residuals (original broadcast coordinates and clock minus corresponding SP3 values) described by colored dots (one color for each Toe); the post-fit residuals are shown by diamonds. The clock correction is expressed in meters. By comparison with GPS, it can be noted that the original ephemeris has, relative to the reference precise orbit, a pattern suggesting periodic resets, to keep the broadcast position and time as close as possible to the precise value, in the neighborhood of the respective Toe. In this way the effects of the third zonal harmonic (J 3 ) are probably less visible than with a longer sampling time, as discussed for G01 in Figure 2. Moving away from this reference epoch, the departure from the reference values increases, so that the ephemeris needs to be replaced by a more recent block near the edge of the validity period. The broadcast ephemeris with best-fitting parameters (post-fit residuals, y axis to the right) tracks instead the precise ephemeris very closely and with no divergence at the edges of the validity interval. The overall improvement for Galileo's E01 in doing the fit described by Equation (1) is summarized as follows: the mean and rms of the pre-fit residuals is −1.224 ± 1.992 m, and 0.000 ± 0.014 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The adjusted clock polynomials fit the SP3 clock data with 0.002 m rms, whereas for the spatial coordinates the corresponding figure is a factor of seven larger. Figure 5 shows the corrections to two clock parameters (a), to the Mean anomaly at epoch Toe and cosine and sine amplitude of the corrections to second harmonic periodic perturbations along track (b). In the bottom left and bottom right plots, we have the cosine and sine amplitudes of the corrections for the radial (c) and across track (d) components, respectively. Table 2 gives the Helmert parameters of the broadcast reference frame relative to the IGS14 frame of the precise ephemeris. Here it is noteworthy that the broadcast frame of E01 has a negligibly small misalignment to IGS14. To account for the CoM-APC offset the offsets in https://www.gsc-europa.eu/support-to-developers/galileo-satellitemetadata (accessed on 13 October 2021) have been used. This calibration implies that the scale factor must be negligibly small, as it results from Table 2. Table 2. Helmert parameters relating the origin, orientation and scale of the Broadcast reference frame of E01 relative to the IGS14 frame of the SP3 precise ephemeris, 2 January 2020. The overall improvement for Galileo's E01 in doing the fit described by Equation (1) is summarized as follows: the mean and rms of the pre-fit residuals is −1.224 ± 1.992 m, and 0.000 ± 0.014 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The adjusted clock polynomials fit the SP3 clock data with 0.002 m rms, whereas for the spatial coordinates the corresponding figure is a factor of seven larger. Figure 5 shows the corrections to two clock parameters (a), to the Mean anomaly at epoch Toe and cosine and sine amplitude of the corrections to second harmonic periodic perturbations along track (b). In the bottom left and bottom right plots, we have the cosine and sine amplitudes of the corrections for the radial (c) and across track (d) components, respectively. Table 2 gives the Helmert parameters of the broadcast reference frame relative to the IGS14 frame of the precise ephemeris. Here it is noteworthy that the broadcast frame of E01 has a negligibly small misalignment to IGS14. To account for the CoM-APC offset the offsets in https://www.gsc-europa.eu/support-to-developers/galileo-satellite-metadata (accessed on 10 October 2021) have been used. This calibration implies that the scale factor must be negligibly small, as it results from Table 2. Table 2. Helmert parameters relating the origin, orientation and scale of the Broadcast reference frame of E01 relative to the IGS14 frame of the SP3 precise ephemeris, 2 January 2020.

Results for Beidou C07
We address here C07, an example of Beidou GNSS, that is an SV in an IGSO orbit: this is a geosynchronous orbit with an inclination comparable to that of GPS or Galileo. The orbit pattern in ECEF coordinates is such that the X and Y coordinates sample a limited range of values. Consequently, one expects a rank deficiency in the normal equations and that the translational Helmert parameters are poorly constrained by the lack of geometry. We have therefore constrained the translational Helmert parameters to zero and solved for Helmert rotations and scale (Table 3), and the seven arc parameters and three clock parameters at intervals of 2 h. For Beidou, the precise SP3 ephemeris computed by CODE were used, being those of CNES not yet available. This also provides a way of using data from more than one Center.
One more aspect to consider is that this SV is modeled with the GPS type of ephemeris which is optimized on an orbital radius of ca. 26,000 km, whereas the geosynchronous orbit has a radius of roughly 42,000 km. It follows that the perturbations due to the Earth's gravity field will be attenuated, and those due to the Sun and the Moon will tend to increase.

Results for Beidou C07
We address here C07, an example of Beidou GNSS, that is an SV in an IGSO orbit: this is a geosynchronous orbit with an inclination comparable to that of GPS or Galileo.  Figures 6 and 7 summarize our results for C07. The broadcast ephemeris shows a discrepancy relative to the precise orbit of the order of the meter, and there is a constant bias in the clock corrections ( Figure 6). After the fit we observe that the post-fit residuals for both coordinates and clock agree very closely with the corresponding precise ephemeris, suggesting that the correction to the 12 sets of seven parameters each (mean anomaly at epoch Toe; amplitudes and phases of the second zonal components) have successfully accounted for the modeling imperfections of the broadcast message.  is summarized as follows: the mean and rms of the pre-fit residuals is −6.753 ± 13.310 m, and 0.000 ± 0.029 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The large bias in the pre-fit residuals is primarily due to the clock offset of the broadcast message relative to the SP3 values, equivalent to nearly 30 m or 10 −7 s (Figure 7a). The rms of the post-fit residuals is dominated by the spatial component: the clock polynomials fit with a typical rms of 0.002 m, about a factor of eight smaller than the rms of the spatial components. The overall improvement for Beidou's C07 in doing the fit described by Equation (1) is summarized as follows: the mean and rms of the pre-fit residuals is −6.753 ± 13.310 m, and 0.000 ± 0.029 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The large bias in the pre-fit residuals is primarily due to the clock offset of the broadcast message relative to the SP3 values, equivalent to nearly 30 m or 10 −7 s (Figure 7a). The rms of the post-fit residuals is dominated by the spatial component: the clock polynomials fit with a typical rms of 0.002 m, about a factor of eight smaller than the rms of the spatial components.

Results for Beidou C12
The SV C12 belongs to the part of the Beidou GNSS which is in Medium Earth Orbit (MEO), like GPS. Therefore, it is meaningful to use Equation (1) and model Helmert transformations on a full day arc, and arc parameters (mean anomaly, cosine and sine amplitudes of the periodic perturbations along track, cross track and radial) on two-hour arcs. Figures 8 and 9 summarize our results for C12. The broadcast ephemeris shows, as for C07, a discrepancy relative to the precise orbit of the order of the meter, and there is a constant bias in the clock corrections ( Figure 8). After the fit we observe that the post-fit residuals for both coordinates and clock agree very closely with the corresponding precise ephemeris, suggesting that the correction to the 12 sets of seven parameters each (mean anomaly at epoch Toe; amplitudes and phases of the second zonal components) plus the seven Helmert parameters have successfully accounted for the modeling imperfections of the broadcast message.
The overall improvement for Beidou's C12 in doing the fit described by Equation (1) is summarized as follows: the mean and rms of the pre-fit residuals is −5.263 ± 9.351 m, and 0.000 ± 0.046 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The large bias in the pre-fit residuals is, similarly to C07, primarily due to the clock offset of the broadcast message relative to the SP3 values.

Results for Beidou C12
The SV C12 belongs to the part of the Beidou GNSS which is in Medium Earth Orbit (MEO), like GPS. Therefore, it is meaningful to use Equation (1) and model Helmert transformations on a full day arc, and arc parameters (mean anomaly, cosine and sine amplitudes of the periodic perturbations along track, cross track and radial) on two-hour arcs. Figures 8 and 9 summarize our results for C12. The broadcast ephemeris shows, as for C07, a discrepancy relative to the precise orbit of the order of the meter, and there is a constant bias in the clock corrections ( Figure 8). After the fit we observe that the post-fit residuals for both coordinates and clock agree very closely with the corresponding precise ephemeris, suggesting that the correction to the 12 sets of seven parameters each (mean anomaly at epoch Toe; amplitudes and phases of the second zonal components) plus the seven Helmert parameters have successfully accounted for the modeling imperfections of the broadcast message.   The overall improvement for Beidou's C12 in doing the fit described by Equation (1) is summarized as follows: the mean and rms of the pre-fit residuals is −5.263 ± 9.351 m, and 0.000 ± 0.046 m for the post-fit residuals. This statistic includes the spatial and temporal coordinates. The large bias in the pre-fit residuals is, similarly to C07, primarily due to the clock offset of the broadcast message relative to the SP3 values. Table 4 summarizes the seven Helmert parameters estimated for C12 on a daily basis. There is a clear shift of the origin in the z direction of ca. 0.7 m, and a smaller one in the x direction. Again, this can be interpreted as a good indication that our optimization process does require an adjustment of the reference frame parameters, if centimetric rms of the post-fit residuals relative to a precise ephemeris is the goal.

Mathematical Model and Results for Glonass
In the Glonass navigation message, the clock corrections of the satellite time to the Glonass time is, as is for GPS, expressed in terms of a second order polynomial. The instantaneous ECEF position and velocity of a Glonass satellite is instead computed by numerically integrating the equations of motion. These are conveniently formulated in terms of six ordinary differential equations of first order. The vector of state of the satellite (po-  Table 4 summarizes the seven Helmert parameters estimated for C12 on a daily basis. There is a clear shift of the origin in the z direction of ca. 0.7 m, and a smaller one in the x direction. Again, this can be interpreted as a good indication that our optimization process does require an adjustment of the reference frame parameters, if centimetric rms of the post-fit residuals relative to a precise ephemeris is the goal. Table 4. Helmert parameters relating the origin, orientation and scale of the Broadcast reference frame of C12 relative to the IGS14 frame of the SP3 precise ephemeris for SV C12, 2 January 2020. The scale factor accounts for the Center of Mass-Antenna Phase Center correction in the radial direction.

Mathematical Model and Results for Glonass
In the Glonass navigation message, the clock corrections of the satellite time to the Glonass time is, as is for GPS, expressed in terms of a second order polynomial. The instantaneous ECEF position and velocity of a Glonass satellite is instead computed by numerically integrating the equations of motion. These are conveniently formulated in terms of six ordinary differential equations of first order. The vector of state of the satellite (position and velocity) is therefore broadcast at a reference epoch Toe. The numerical integrator, normally a 4th order Runge-Kutta, maps this vector of state from time Toe to any other epoch within the validity interval of the message. The force field consists of the gradient of the Earth's gravitational potential truncated to the second zonal harmonic J 2 plus the centrifugal and Coriolis terms, since the equations of motion are integrated in a rotating frame. Glonass broadcasts three additional terms, the Lunisolar accelerations, which are constant accelerations during the validity time of the broadcast message, normally 30 min.
Tuning the broadcast parameters on a precise orbit requires therefore the adjustment of the clock and vector of state in arcs of at least 30 min. In the sample broadcast ephemeris data set analyzed here (https://igs.bkg.bund.de/root_ftp/MGEX/BRDC_v3/2020/002/ brdm0020.20p.Z, accessed on 10 October 2021), the clock drift and clock drift rate (a1, a2) were zero, and the lunisolar accelerations were allowed to change in fixed increments. It is also important to note that the Glonass message constrains the terrestrial reference frame by providing ECEF initial coordinates and velocities. Therefore, for Glonass it is not necessary to estimate Reference Frame parameters. For the CoM to APC correction the z-bias of 2.450 m appropriate for Glonass M was used (https://files.igs.org/pub/station/ general/pcv_archive/igs14_2056.atx, accessed on 10 October 2021).
To tailor the broadcast parameters on a precise ephemeris and clock, we formulate the minimum variance algorithm as follows: We minimize the sum of the 96 × 4 = 384 square discrepancies ∆(XYZT) between the SP3 precise values and those computed with the broadcast message, using as adjustable parameters the three clock terms a0, a1 and a2, and a 9D vector of state containing position and velocities at epoch Toe, and three constant accelerations. If the 24-h interval is broken into 24 consecutive arcs each of one hour, then we have to adjust 288 parameters. These become 144 if we test arcs of a 2-h duration. The Time of clock in the clock polynomial and the Time of ephemeris are taken coincident.
Contrary to the GPS-like approach, where the arc terms were coupled by the global terms, i.e., the Helmert parameters (Figure 1), the partial derivative matrix set up to linearize Equation (2) is strictly block diagonal, with no correlation between arc parameters of different arcs ( Figure 10).
Tuning the broadcast parameters on a precise orbit requires therefore the a of the clock and vector of state in arcs of at least 30 min. In the sample broadcast data set analyzed (https://igs.bkg.bund.de/root_ftp/MGEX/BRDC_v3/2020/002/brdm0020.20p.Z on 13 October 2021), the clock drift and clock drift rate ( 1, 2) were zero, and t lar accelerations were allowed to change in fixed increments. It is also importa that the Glonass message constrains the terrestrial reference frame by provid initial coordinates and velocities. Therefore, for Glonass it is not necessary t Reference Frame parameters. For the CoM to APC correction the z-bias of 2.450 priate for Glonass M was used (https://files.igs.org/pub/station/gener chive/igs14_2056.atx, accessed on 13 October 2021).
To tailor the broadcast parameters on a precise ephemeris and clock, we the minimum variance algorithm as follows: Contrary to the GPS-like approach, where the arc terms were coupled by terms, i.e., the Helmert parameters (Figure 1), the partial derivative matrix set arize Equation (2) is strictly block diagonal, with no correlation between arc p of different arcs ( Figure 10). We consider on 2 January 2020, satellite R01 and precise ephemeris and clock computed by CNES (https://cddis.nasa.gov/archive/gnss/products/mgex/2086/GRG0MGXFIN_20 200020000_01D_15M_ORB.SP3.gz, accessed on 10 October 2021). The results of the parameter optimization for this test case are shown in Figures 11 and 12. Figure 11 suggests that also in the example of Glonass it is possible to fine tune the broadcast model parameters so that the positions and clock corrections computed with the adjusted parameters very closely fit the reference CNES final orbit of R01. The corrections to the broadcast parameters are shown graphically in Figure 12. Besides the rms of 0.024 m as an indication of the final accuracy, it should be mentioned as an additional benefit that the interval of validity of the corrected broadcast ephemeris has been doubled to 1 h. Tests with 2-h arcs resulted in rms spreads considerably larger. We propose 1 h as a reasonable compromise between accuracy and refresh rate.

Polynomial Clock Model vs. IGS/MGEX High Rate Clocks
The clock polynomials we have computed for GPS, Galileo, Beidou and Glonass fit with centimetric rms the clock values tabulated at 15 min intervals in the SP3 files. For most applications, the polynomial clock corrections are used to interpolate the values at rates of the order of 1 Hertz or higher. In [2], it is pointed out that noise in rubidium, cesium and hydrogen maser clocks is characterized by a random walk phase modulation. For lag times of up to 10 s, the Allan variance is between 10 −11 and 10 −12 for the Cesium or Rubidium clocks. For Galileo's Passive Hydrogen Maser, the frequency stability is somewhat smaller than 10 −12 [18,19,24,26]. This means that two consecutive, non-overlapping time segments of 100 s nominal length will have a 1 sigma difference of 1 ns to 0.1 ns for a two-sample Allan stability of 10 −11 and 10 −12 respectively. It follows that our clock polynomials need to be tested against high rate clock estimates. The IGS/MGEX makes available such files with sampling rate of 30 s. In the rest of this section, we will therefore compute differences in the clock corrections, between our polynomial values based upon 15 min sampling and high rate IGS clocks based on 30 s sampling. Figure 13 gives a comparative example of IGS/MGEX high rate clocks (30 s sampling) and the predictions of our polynomials based on best fit to 15 min data (http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/2020/COM20864.CLK.Z for Beidou and https://cddis.nasa.gov/archive/gnss/products/mgex/2086/GRG0MGXFIN_20200020000_01D_30S_CLK.CLK.gz for GPS, Glonass and Galileo, accessed on 13 October 2021). Figure 13 indicates that the rms discrepancy is smaller than 0.1 ns for GPS, Galileo and Beidou, while for Glonass we have a factor 10 worse rms. Several departures from random noise are well visible, suggesting that our clock polynomials smooth the high frequency part of the noise. The differences are on average of up to few equivalent cm or less. This implies that an overall figure of merit for our improved broadcast model of a few cm is the accuracy level we can expect, even at high sampling rates. These rms estimates are within the consistencies among the various MGEX solutions for clocks and orbits, as discussed in [16].

Polynomial Clock Model vs. IGS/MGEX High Rate Clocks
The clock polynomials we have computed for GPS, Galileo, Beidou and Glonass fit with centimetric rms the clock values tabulated at 15 min intervals in the SP3 files. For most applications, the polynomial clock corrections are used to interpolate the values at rates of the order of 1 Hertz or higher. In [2], it is pointed out that noise in rubidium, cesium and hydrogen maser clocks is characterized by a random walk phase modulation. For lag times of up to 10 s, the Allan variance is between 10 −11 and 10 −12 for the Cesium or Rubidium clocks. For Galileo's Passive Hydrogen Maser, the frequency stability is somewhat smaller than 10 −12 [18,19,[24][25][26]. This means that two consecutive, non-overlapping time segments of 100 s nominal length will have a 1 sigma difference of 1 ns to 0.1 ns for a two-sample Allan stability of 10 −11 and 10 −12 respectively. It follows that our clock polynomials need to be tested against high rate clock estimates. The IGS/MGEX makes available such files with sampling rate of 30 s. In the rest of this section, we will therefore compute differences in the clock corrections, between our polynomial values based upon 15 min sampling and high rate IGS clocks based on 30 s sampling. Figure 13 gives a comparative example of IGS/MGEX high rate clocks (30 s sampling) and the predictions of our polynomials based on best fit to 15 min data (http://ftp.aiub.unibe. ch/CODE_MGEX/CODE/2020/COM20864.CLK.Z for Beidou and https://cddis.nasa.gov/ archive/gnss/products/mgex/2086/GRG0MGXFIN_20200020000_01D_30S_CLK.CLK.gz for GPS, Glonass and Galileo, accessed on 10 October 2021). Figure 13 indicates that the rms discrepancy is smaller than 0.1 ns for GPS, Galileo and Beidou, while for Glonass we have a factor 10 worse rms. Several departures from random noise are well visible, suggesting that our clock polynomials smooth the high frequency part of the noise. The differences are on average of up to few equivalent cm or less. This implies that an overall figure of merit for our improved broadcast model of a few cm is the accuracy level we can expect, even at high sampling rates. These rms estimates are within the consistencies among the various MGEX solutions for clocks and orbits, as discussed in [16].

Improving the Broadcast Model for GPS-like and Glonass-like Messages
We have seen indications that both the GPS-like and Glonass-like navigation messages have a potential for high accuracy, in the sense that the model coefficients can be tuned to a precise ephemeris with centimetric rms spread relative to a precise SP3 ephemeris. At these levels of accuracy, it is meaningful to ask if the model implemented in both types of Broadcast Ephemeris is adequate to such accuracies.
For GPS-like messages, we have secular perturbations in inclination and node, and periodic perturbations in mean anomaly, in the radial direction and inclination. The period is one half the Keplerian period, roughly 6 h, implying that the periodic perturbations caused by the part of the Earth's gravity field are modeled. For Glonass-like messages, the Glonass ICD (Interface Control Document) prescribes a force field which includes, for the Earth's gravity, the monopole and quadrupole ( ) components.

Improving the Broadcast Model for GPS-like and Glonass-like Messages
We have seen indications that both the GPS-like and Glonass-like navigation messages have a potential for high accuracy, in the sense that the model coefficients can be tuned to a precise ephemeris with centimetric rms spread relative to a precise SP3 ephemeris. At these levels of accuracy, it is meaningful to ask if the model implemented in both types of Broadcast Ephemeris is adequate to such accuracies.
For GPS-like messages, we have secular perturbations in inclination and node, and periodic perturbations in mean anomaly, in the radial direction and inclination. The period is one half the Keplerian period, roughly 6 h, implying that the periodic perturbations caused by the J 2 part of the Earth's gravity field are modeled.
For Glonass-like messages, the Glonass ICD (Interface Control Document) prescribes a force field which includes, for the Earth's gravity, the monopole and quadrupole (J 2 ) components.
To investigate the effect of higher order terms of the gravitational potential on the orbital elements, we refer to [27]. The differential equations relating a disturbing potential U to the temporal changes of the orbital elements ε k (k = 1:6) have the general form: where L k (a, e, I) is a linear differential operator dependent on the semimajor axis a eccentricity e and inclination I. If we assume that interaction of perturbations can be ignored (this is appropriate for all the harmonics except J 2 ), then we can write: where δε k lm is the perturbation of the k-th orbital element due to the harmonic C lm and ε k 0 is the unperturbed orbital element. The general form of the perturbed orbital elements is .
The coefficients A lmpq scale as (a e /a) l , e q and sin p (I), cos p (I), and abs(q) < 5 typically. θ are respectively the mean anomaly, the rate of the perigee, the mean motion, the rate of the ascending node and the Earth rotation rate. Secular terms arise when m = 0 and (l − 2p) = q = 0 that is for even zonal harmonics. Because the rate of perigee and node are of the order of 10 −3 n and the period is n ∼ = 2 . θ, the period P of a perturbation is primarily determined by: 2π It is observed that the amplitude of the harmonics C lm decrease with the degree l according to the Kaula's rule: Equation (9) implies that, except for resonant terms, we have to consider small degree terms. For example, considering the J 3 zonal (l = 3, m = 0), for p = q = 0 we have perturbations with period close to 4 h. The J 3 term, accounts for the pear-shaped figure of equilibrium of the Earth, that is, the lack of symmetry between north and south hemispheres. The conventional values from the JGM-3 gravity model [28] are: We can estimate a maximum value of the perturbation caused by J 3 for an orbit of roughly 26,000 km by scaling the maximum periodic perturbation due to J 2 , which along track has a maximum amplitude of ca. 200 m: The period of the associated perturbation is 1/3 of the orbital frequency, which is close to 4 h. We suggest that both GPS-like and Glonass-like messages should account for such term. The changes in the respective messages which could accommodate the improved gravity model are described by the following equations for a GPS-like ephemeris and a Glonass-like model. For GPS, denoting by uk the true anomaly from the node, we need to introduce six additional coefficients c3 uc , c3 us , c3 ic , c3 is , c3 rc , c3 rs which define the amplitude of the cosine and sine component of the third harmonic periodic perturbation in the along track, cross track and radial directions, respectively. This notation is similar to that adopted (c uc , . . . , c is ) to represent the amplitudes of the perturbations caused by the second harmonic. This approach was first proposed by Zhou et al. [29]. The total periodic perturbations due to J 2 and J 3 are then given by: duk = c uc * cos(2uk) + c us * sin(2uk) + c3 uc * cos(3uk) + c3 us * sin(3uk) (12) drk = c rc * cos(2uk) + c rs * sin(2uk) + c3 rc * cos(3uk) + c3 rs * sin(3uk) For Glonass, we simply need to modify the force function used in the numerical integration of the equations of motion in a co-rotating frame. The additional terms are marked in red.
For the GPS-like approach, the amplitude of the 3rd harmonic perturbations need to be estimated by including them in the set of solve for parameters (Equation (1)) and the relevant partials in the partial derivative matrix. These six additional parameters would justify an increase in length of the validity interval to 4 h or perhaps more, so that in a daily fit the total number of parameters to be estimated with 4-h arcs is 7 (Helmert) + 6 (arcs) * [3 (clocks) + 6 (orbit)] = 103 parameters, smaller than the 127 parameters of Section 2. Appropriate words would have to be defined in the Navigation message to allocate space to the six additional terms. For Glonass, there would be no change in the vector of state, nor in the navigation message. Only the function which describes the force needs to be modified.
We have seen in Figure 2 that the fit interval of the first arc of G01 was extended to 4 h to test the possibility of extending to 4 h the nominal validity time of 2 h. Figure 14 shows the detail of the post-fit residuals in this 4-h arc. The blue dots represent the best fitting signal we computed by modifying the GPS navigation message according to Equations (12)- (14). The similarity is evident, but for this short arc it is difficult to draw a sufficiently motivated conclusion.

Results for Volume Calculations
In the previous sections, we presented results based on limited data, in order to privilege the detail of the analysis more than the volume calculations on a larger number of days and satellites. In this section, we apply the same approach systematically to a larger body of data. We concentrate on the full month of January 2020.
In the Appendix, Tables A1-A4 contain the Helmert parameters and their standard deviations for GPS, Glonass, Galileo and Beidou, respectively, averaged over one month. The Helmert parameters were estimated by modeling the differences between the ECEF coordinates computed with the broadcast ephemeris and those available in the SP3 file of CNES, on a satellite-by-satellite basis. For BeiDou, the CODE SP3 ephemeris were used.

Results for Volume Calculations
In the previous sections, we presented results based on limited data, in order to privilege the detail of the analysis more than the volume calculations on a larger number of days and satellites. In this section, we apply the same approach systematically to a larger body of data. We concentrate on the full month of January 2020.
In the Appendix A, Tables A1-A4 contain the Helmert parameters and their standard deviations for GPS, Glonass, Galileo and Beidou, respectively, averaged over one month. The Helmert parameters were estimated by modeling the differences between the ECEF coordinates computed with the broadcast ephemeris and those available in the SP3 file of CNES, on a satellite-by-satellite basis. For BeiDou, the CODE SP3 ephemeris were used. Then the 31 daily values were averaged. For GPS, BeiDou and Glonass, the mean of each of the seven parameters is zero within one standard deviation, implying that on average the broadcast and SP3 reference frames are aligned. For Galileo, the only parameters which are nonzero within one standard deviation are Rx and Rz. In general, the broadcast and SP3 frames are aligned. However, individual values (e.g., Rz of G01 or G11) can be large enough to generate significant departures of the broadcast predictions relative to the reference SP3 positions. Consequently, it is advisable, for maximum accuracy, to make available the Helmert parameters in conjunction with the corrections to the broadcast ephemeris, for each satellite and day.
We now turn to evaluate the effectiveness of our proposed approach to account for unmodeled reference frame adjustments and orbit effects in the broadcast ephemeris on all the satellites and days of January 2020. To this purpose, we have the average spectra of the pre-fit and post-fit residuals for each constellation. Figures 15-18 show the results for GPS, Glonass, Galileo and Beidou (IGSO and MEO), respectively. The spectra are computed in the Radial, Tangential and Across track triad. The upper plot shows the spectrum of the pre-fit residuals, i.e., the differences between broadcast and SP3 coordinates prior to the adjustment. The lower plot shows the spectra of the post-fit residuals, i.e., after the fit. The plots clearly show that the pre-fit spectra are dominated by low frequency signals (periods larger than 6 h, or frequency less than 0.05 mHz) with amplitudes of the order of 1 m. The lower plots show that the parameter adjustment proposed in this paper effectively removes the low frequency signals, so that the spectra of the differences between the positions computed with the corrected broadcast model and the SP3 positions are very nearly flat. The scales of the y axis are different. In particular for Galileo, it is worth noting that the proposed algorithm removes from the original broadcast ephemeris most of the unmodeled perturbations, whereas for Glonass we observe a lower efficiency. 1 m. The lower plots show that the parameter adjustment proposed in this paper effectively removes the low frequency signals, so that the spectra of the differences between the positions computed with the corrected broadcast model and the SP3 positions are very nearly flat. The scales of the y axis are different. In particular for Galileo, it is worth noting that the proposed algorithm removes from the original broadcast ephemeris most of the unmodeled perturbations, whereas for Glonass we observe a lower efficiency.  tively removes the low frequency signals, so that the spectra of the differences between the positions computed with the corrected broadcast model and the SP3 positions are very nearly flat. The scales of the y axis are different. In particular for Galileo, it is worth noting that the proposed algorithm removes from the original broadcast ephemeris most of the unmodeled perturbations, whereas for Glonass we observe a lower efficiency.

Discussion
The sample of data we have tested suggests that the GPS-like (Keplerian orbit plus perturbations) navigation message can be tuned to a precise ephemeris and clock with an rms spread ranging from 1.4 cm (Galileo) to 5.2 cm (GPS), provided that the corrections to the broadcast parameters are complemented by daily adjustments in origin and orientation of the reference frame parameters. These adjustments are negligibly small for Galileo but can be of several tens of cm for GPS and Beidou. For an IGSO orbit like Beidou C07, the translation parameters are undefined due to lack of geometry. It may be expected that similar results can be obtained for Japan's QZSS which are also in a IGSO orbit, and for India's NAVIC/IRNSS.
The validity time of the ephemeris is two hours. For GPS, we tested a fit interval of 4 h. In such a case, the post-fit residuals indicate the presence of an oscillation of an approximately 4-h period, which is the signal expected to be caused primarily by the third zonal (J 3 ) of the Earth's gravity field. Including this perturbation in the navigation message would help in increasing the validity time and have a more random spread of the post-fit residuals.
The inference is that the format of the broadcast ephemeris adopted by GPS, Galileo, Beidou, and possibly NAVIC/IRNSS and QZSS, can represent the SV position to very high accuracy, comparable with that of the precise ephemeris. Comparison of power spectra of the Broadcast-SP3 residuals before and after the adjustment computed for all the satellites and one month support this hypothesis. Satellites in the GEO orbit have not been tested because a reference precise ephemeris is at this time unavailable.
For the clock corrections, a similar reasoning applies. The three-parameter model can likewise be tuned on precise clock corrections at 15 min sampling. It is to be noted, however, that flicker noise at sampling rates of up to 100 s can increase considerably the Allan variance for Cesium or Rubidium clocks. For the Galileo's Passive Hydrogen Maser this is also true but not so marked as for Cesium or Rubidium clocks. Therefore, the effective validity of the polynomial clock model at sampling rates of the order of 1 Hz needs to be tested with IGS's high rate clocks.
One of the most important challenges in satellite navigation precise positioning is to be able to broadcast corrections to the navigation message so that the user has positions and clocks of the used GNSS's of sufficient accuracy to apply for example a PPP (Precision Point Positioning) algorithm.
Currently, the International GNSS Service (IGS) agency and various analysis centers (ACs) provide users with ultra-rapid precise satellite orbit and clock products for GPS and Glonass with 6-h updates and a 3-h latency (igs.org/acc/gps-only/#ultra-rapid, accessed on 10 October 2021). Corresponding products for Galileo and Beidou are discussed by [30]. The accuracy of the predicted orbits is not good enough for high-precision PPP [25]. A widely used approach to generate high accuracy orbits and clocks usable in Real Time PPP is the SSR (State Space Representation). In the SSR concept (https://www.igscb.org/wpcontent/uploads/2020/10/igs_ssr_-v1_00_20201005.pdf, accessed on 10 October 2021), the broadcast ephemeris data are complemented with a set of corrections along track, across track and radial, which are described by a piecewise linear function of time. Likewise, the clock corrections to the broadcast values are described by a piecewise quadratic function of time. Tests done on several Real Time and Orbit and Clock products indicate an update rate of typically 5 s, with a comparable latency. Specific RTCM messages 1060 and 1066 have been defined for GPS and Glonass SSR, respectively.
We suggest that if sufficiently accurate orbit and clock predictions are available, the navigation message could be broadcast in a corrected form as a streamed RTCM (Radio Technical Commission for Maritime Services, www.rtcm.org, accessed on 10 October 2021) message as described in this paper, with additional reference frame information for GPSlike messages. Given the validity time of two hours, perhaps extendable to four hours with a refined gravity model, for GPS-like messages, and one hour for Glonass, the refresh rate could be considerably longer than SSR. If the SSR corrections are based on the improved navigation message, there would be a twofold advantage: (a) the SSR corrections would address the finest details of the orbit and clocks, for improved overall accuracy; (b) the user would have a redundancy in case of unavailable SSR data, as he would still be able to compute its position with a high accuracy ephemeris in a broadcast rather than SP3 format.
Applicable RTCM messages could be 1019, 1020, 1042, 1045/6 for GPS, Glonass, Beidou and Galileo I/NAV and F/NAV, respectively (the two Galileo navigation messages refer to the two different ionospheric free linear combinations and corresponding clocks), and 1021 for the Helmert parameters [31].
In conclusion, we suggest that a corrected broadcasted message as it has been presented here, rather than the message broadcast by the various GNSSs, could be the basis for an even more accurate set of SSR corrections for both position and clocks.

Conclusions
The limited accuracy of the broadcast data has been discussed in detail in several publications. Possible actions towards an improvement of the broadcast message are the focus of this paper. We have investigated the potential of the GPS and Glonass navigation messages to represent satellite position and clock with an accuracy comparable with that of the precise ephemeris delivered by the IGS within the MGEX project. The reference precise products (ephemeris and clocks) were generated by CNES for GPS, Glonass and Galileo, and by CODE for Beidou (MEO and IGSO). We have examined one satellite per constellation in one specific day, and obtained indication that the broadcast message can reproduce the precise ephemeris with a centimetric rms, provided that: (a) for GPS-like messages (GPS, Galileo and Beidou), arcs of two hours are used for the fit, complemented by one set of Helmert parameters for the day; and (b) for Glonass-like messages arcs of one hour are used for the fit. For Glonass the Helmert parameters are unnecessary because the Glonass message is based directly on Cartesian ECEF coordinates and velocities. Volume calculations on all satellites and for one month indicate that the proposed approach is applicable in general.
Corrections to the clock parameters are also computed, and similar accuracy has been demonstrated. High frequency clock files from IGS have also been tested against the estimated clock polynomials. The results suggest that the rms spread is limited to fractions of nanoseconds for GPS, Galileo and Beidou, and a factor of 10 higher for Glonass. Systematic drifts in the clock differences are clearly visible, implying that our clock polynomial is likely to smooth the high frequency noise of the on board clocks, particularly for the Rubidium and Cesium clocks. Order of magnitude estimates of the perturbations caused by the higher order terms of the gravity field suggest the opportunity to modify the message to include the effect of the J 3 component of the Earth's gravity field.
Implications of our results for real time applications very much depend on the availability of predicted precise orbits and clocks which can be represented in a broadcast form. In such case, appropriate RTCM messages are available for the examined constellations, both for the broadcast message and the complementary Helmert parameters.