Introduction

Precise point positioning (PPP) has emerged more than two decades ago as an alternative method to differential carrier-phase positioning (Malys and Jensen 1990; Zumberge et al. 1997). While the differential real-time kinematic (RTK) approach relies on the elimination of measurement errors or nuisance parameters through differential positioning with respect to a real or virtual reference station, the PPP technique relies on accurate models and correction information to achieve a similar positioning accuracy (Héroux et al. 2004; Kouba and Héroux 2001).

Aiming at decimeter to centimeter level accuracy, both offline and real-time applications of PPP build on the use of the most accurate orbit and clock products. With limited exceptions (Gunning et al. 2019; Hadas et al. 2019), broadcast ephemerides have not been considered as a suitable option for PPP in view of their limited accuracy. However, two new GNSS constellations, Galileo and BeiDou-3, have advanced and are now close to or already in operational status. They provide broadcast ephemerides with significantly smaller orbit and clock errors than the legacy systems. Galileo, in particular, has a global average root mean square (RMS) signal-in-space range error (SISRE) of only 20 cm, which is approximately a factor of three smaller than for GPS (Montenbruck et al. 2018). The main driver behind the low SISRE are frequent uploads of updated broadcast ephemerides to the Galileo satellites, together with the use of very precise passive hydrogen maser clocks through the entire constellation. Modernized and more stable clocks have also been deployed in the GPS constellation with the Block IIF satellites. Similar and even better clocks are also used in the latest generation GPS III satellites, which promise a further decrease in the SISRE over the next years.

The low SISRE of Galileo makes this constellation already today the most promising candidate for PPP with broadcast ephemerides in applications aiming at accuracy from sub-meter to a few decimeters. The application of broadcast ephemerides in a PPP model appears of particular interest for real-time positioning, as it eliminates the dependency on a PPP correction data stream.

To minimize the impact of broadcast ephemeris errors in PPP, we study two different strategies. Following an analysis of SISREs of GPS and Galileo to characterize the current quality of the respective broadcast ephemerides, the relevant algorithms are presented. Subsequently, practical tests with permanent reference stations are presented to assess the achievable accuracy of PPP with broadcast ephemerides. Furthermore, these are complemented by a boat test to demonstrate the practical use of the method on a moving vehicle under realistic conditions.

Assessment of GPS and Galileo broadcast ephemeris errors

The achievable PPP accuracy with broadcast ephemeris (BCE) products mainly depends on the magnitude of orbit and clock errors in those products. For a comprehensive assessment of these errors, the SISRE is typically used as a metric (Montenbruck et al. 2018). A recent analysis yields SISRE RMS values of approximately 21 cm for Galileo and 35 cm for GPS Block IIF satellites (Wu et al. 2020). For a proper understanding of the PPP tests with permanent reference stations described in the next section, a dedicated SISRE analysis was performed for the entire month of December 2019, considering all healthy satellites in the GPS and Galileo constellation. GPS and Galileo broadcast ephemeris products from the International GNSS Service (IGS; Johnston et al. 2017) have been compared with precise products from the Center for Orbit Determination in Europe (CODE; Prange et al. 2020a, b). The results are depicted in Fig. 1, which shows the cumulative distribution of SISRE values for both constellations.

Fig. 1
figure 1

Cumulative distribution of SISRE values for GPS and Galileo for December 2019

Full size image

The plot confirms the finding that the Galileo constellation exhibits significantly smaller broadcast orbit and clock offset errors compared to GPS. The RMS SISRE for Galileo reaches a value of only 12 cm during this time period, while GPS has an RMS SISRE of 50 cm. The 95th-percentile values are roughly twice as high and amount to 24 cm for Galileo and 97 cm for GPS. The significantly lower errors for Galileo can be explained by the use of highly stable passive hydrogen maser clocks onboard most satellites, limiting the clock prediction error. Second, Galileo offers a much higher upload rate of the broadcast navigation data compared to GPS, which reduces orbit and clock extrapolation errors. It can be concluded from this analysis that the errors which need to be compensated by the additional SISRE parameter in the state vector have different magnitude for both constellations and also have different temporal behavior. Therefore, different initial standard deviations, as well as process noise settings, must be used.

In addition to the SISRE values, the magnitude of orbit and clock offset discontinuities are relevant. The SISRE modeling in the Kalman filter assumes only small temporal variations within limits defined by the process noise. This assumption is, of course, violated when a handover of two consecutive batches of broadcast ephemeris sets happens. A change in the issue-of-data counter indicates such handovers for ephemeris (IODE) or clock (IODC) in case of GPS or in the IODnav counter in case of Galileo. The modeled range from satellite- to receiver-antenna may be affected by a discontinuity in either the projected orbit onto the user line-of-sight (LOS) vector or a discontinuity in the clock offset, or both. Depending on the magnitude of the discontinuity, it can either go unnoticed or lead to a rejection of measurements, re-initializing the SISRE state after data quality control.

Orbit and clock offset discontinuities on handovers of consecutive broadcast ephemeris records have been analyzed for the same time period for all satellites of GPS and Galileo constellations. The discontinuities have been determined from the difference between the current and the previous broadcast ephemeris data evaluated at the epoch at which a new broadcast ephemeris set has become available. The clock discontinuities are simply the difference between the clock offsets of the current and the previous broadcast ephemeris set. The orbit discontinuities are computed as the orbit difference between the current and the previous set mapped onto the LOS vector of a ground-based user at the worst user location (WUL) (Li et al. 2011). This user location corresponds to the position on the earth, within the visibility area of a satellite, at which the discontinuity has the largest impact on the range between satellite and receiver. It thus yields a conservative assessment of the impact of orbit discontinuities.

The corresponding results are depicted in Figs. 2 and 3, which show the frequency of occurrence of discontinuities of a certain size with a bin size of 5 cm. Figure 2 shows the statistics of orbit discontinuities. It can be observed that about 75% of all discontinuities for Galileo are less than 5.0 cm, whereas GPS is more often affected by larger values. The 95th-percentile is 9.9 cm for Galileo and 48.7 cm for GPS. Similar to the overall SISRE, this difference can be explained by the shorter update interval of the orbit information onboard the Galileo satellites.

Fig. 2
figure 2

Frequency of occurrence of 3D orbit discontinuities at the WUL for GPS and Galileo for December 2019

Full size image
Fig. 3
figure 3

Frequency of occurrence of clock offset discontinuities for GPS and Galileo for December 2019

Full size image

The statistics of the clock discontinuities are depicted in Fig. 3. Like the orbit, most of the clock discontinuities for Galileo are smaller than 5.0 cm. Almost 85% of all discontinuities fall into this bin. A higher percentage of larger discontinuities is present for GPS. The 95th-percentile is 10.2 cm for Galileo and 26.7 cm for GPS. It should be noted that orbit and clock discontinuities exceeding 0.5 m are summarized in the right-most bar. Only very few of these large discontinuities exist, which can reach magnitudes of a few meters.

All static and dynamic tests presented in our research are performed with a data interval of 30 s. At this time scale, the time variations in the clock errors are bigger than those of the orbit, dominating the rate of change. The consequence is that the time variation of the SISRE not only differs between constellations but also between different satellites. Those with more stable clocks, such as the GPS IIF and GPS III rubidium clocks, and the Galileo hydrogen masers, are affected by a smaller variation of the clock error and need smaller process noise. Satellites with less stable clocks, such as the IIR rubidium and IIF cesium atomic frequency standards, require higher process noise to account for these variations. For simplicity, no satellite-specific settings are applied in our tests. However, different process noise is assigned for the GPS and Galileo constellation to account for the fact that the GPS system has a larger percentage of satellites with higher clock noise, whereas Galileo uses predominantly highly stable passive hydrogen masers.

PPP algorithm description

The functional models of pseudorange and carrier-phase observations in PPP rely on using precise correction data or models to eliminate unknown terms in the measurement equations. As it is common practice, dual-frequency observation combinations are used to remove the ionospheric delay up to first order. We start the derivation of the observation equation with these standard models for code and phase measurements (Kouba et al. 2017)

$$p =\rho +\xi + c \left(\mathrm{ d}{t}_{r}-\mathrm{d}{t}^{s} \right)+\left(T+\mathrm{d}T\right)+e,$$
(1)

and

$$\varphi =\rho + \zeta + c \left(\mathrm{ d}{t}_{r}-\mathrm{d}{t}^{s} \right)+\left(T+\mathrm{d}T\right)+\lambda \left(A+\omega \right)+\epsilon ,$$
(2)

where \(p\) and \(\varphi\) are the ionosphere-free combination of pseudorange and carrier-phase measurements, \(\rho\) is the geometrical range between the satellite’s and the receiver’s antenna reference points, \(\xi\) and \(\zeta\) are the corrections for code and phase center offsets for transmitting and receiving antennas, \(c\) is the speed of light, \(\mathrm{d}{t}_{r}\) and \(\mathrm{d}{t}^{s}\) are the receiver and satellite clock offsets, \(T\) is the modeled tropospheric delay, \(\mathrm{d}T\) is an additional, estimated tropospheric delay correction, \(\lambda\) is the wavelength of the ionosphere-free combination, \(A\) is the float-valued ionosphere-free combination of carrier phase ambiguities, \(\omega\) is the carrier-phase wind-up, and where \(e\) and \(\epsilon\) are the combined noise and multipath errors for pseudorange and carrier-phase. The user position coordinates \({x}_{r}\), \({y}_{r}\) and \({z}_{r}\) are included in the geometric range

$$\rho =\sqrt{({x}^{s}-{x}_{r}{)}^{2}+({y}^{s}{-y}_{r}{)}^{2}+({z}^{s}-{z}_{r}{)}^{2}},$$
(3)

where \({x}^{\mathrm{s}}\), \({y}^{\mathrm{s}}\) and \({z}^{\mathrm{s}}\) are the coordinates of the satellite.

The observations are processed in a Kalman-filter, which estimates the position, receiver clock offset, and ambiguities \({A}_{0} ... {A}_{n}\) for \(n\) tracked satellites as part of the filter’s state vector

$$\underline{x} =\left[ {x}_{r} \ {y}_{r} \ {z}_{r} \ \mathrm{ d}{t}_{r} \ \mathrm{ d}T \ {A}_{0} ... {A}_{n}\right].$$
(4)

When computing dual-constellation solutions, a receiver clock offset for each constellation is estimated. The functional model for pseudorange and carrier-phase observations in (1) and (2) and the state vector in (4) represent a typical formulation for PPP with float ambiguities and ionosphere-free dual-frequency observations. It would normally be used with precise orbit and clock products, and the ambiguities are treated as constant parameters. The typical PPP methodology is here referred to as PRE. The same formulation is also used in a second method, where precise orbit and clock products are replaced with broadcast ephemerides. This method does not include any strategy for SISRE compensation. It is here referred to as BCE.

For simplicity, no motion model is used, and all states are predicted as constants in the time update. For state vector components treated as random walk parameters, white process noise with variance \(q={\sigma }_{P}^{2} \cdot \Delta t/{\tau }_{P}\) at a sampling interval \(\Delta t\) is used in most cases. For the clock offset, process noise with a very large variance \(q={\sigma }_{\mathrm{P}}^{2}\) is applied irrespective of the filter step size. In this way, clock offsets are essentially estimated as free parameters at each epoch. The filter settings for the initial standard deviation \({\sigma }_{0}\), process noise standard deviation \({\sigma }_{\mathrm{P}}\) and time constant \({\tau }_{P}\) are shown in Table 1.

Table 1 Kalman-filter noise settings for state vector elements
Full size table

As a first approach for SISRE compensation, named BCE1, we use the same model with broadcast orbits and clocks. In contrast to the standard PPP modeling, process noise is applied in our case to the ambiguities to account for variations of the SISRE over time. Any unmodeled biases, like the SISRE, would be partly compensated by the float ambiguities in the carrier phase model but essentially be neglected in the pseudorange model. This approach has earlier been suggested in Montenbruck and Ramos-Bosch (2008) for real-time orbit determination of earth orbiting satellites and found to be effective for the compensation of broadcast ephemeris errors in that application.

The second approach, named BCE2, is to include the projected orbit and clock errors as additional parameter \(s\) into the pseudorange and carrier-phase observation equations as suggested by Gunning et al. (2019)

$$p =\rho +\xi + c \left(\mathrm{ d}{t}_{r}-\mathrm{d}{t}^{s} \right)+\left(T+\mathrm{d}T\right)+s+e,$$
(5)

and

$$\varphi =\rho + \zeta + c \left(\mathrm{ d}{t}_{r}-\mathrm{d}{t}^{s} \right)+\left(T+\mathrm{d}T\right)+s + \lambda \left(A+\omega \right)+\epsilon .$$
(6)

The following parameters of (5) and (6) are now part of the filter’s state vector

$$\underline{x} =\left[ {x}_{r} \ {y}_{r} \ {z}_{r} \ {\text{d}}{t}_{r} \ {\text{d}}T \ {s}_{0} ... {s}_{n} \ {A}_{0} ... {A}_{n}\right]$$
(7)

where \({s}_{0} ... {s}_{n}\) are the projected SISRE values for \(n\) tracked satellites. In this approach, SISREs are explicitly estimated using a separate parameter for each tracked satellite and process noise is applied to account for their temporal variations. The initial standard deviation for the SISRE parameters in Table 1 corresponds to the RMS value from the broadcast ephemerides assessment of Fig. 1.

The process noise values for methods BCE1 and BCE2 were defined from the results of a sensitivity analysis. The basic assumption of such analysis is that there exists a value of the process noise that yields the best accuracy. This minimum can, therefore, be selected as the best process noise. The same test case was thus repeated using different values for the process noise standard deviation, over a range from 0 to 15 mm. The process noise was applied at a 30 s measurement update interval. Observations from a pool of 11 stations, chosen from the IGS network and depicted in Fig. 4, were used to compute 24 h single-constellation solutions over a month (December 2019). The solution was computed in static conditions for both methods BCE1 and BCE2. The three-dimensional (3D) RMS position errors were averaged over all tests to obtain a single accuracy metric for each process noise value. The sensitivity analysis was performed separately for BCE1 and BCE2 and for both constellations to obtain a process noise value for each case. The results are depicted in Fig. 5.

Fig. 4
figure 4

Map with the locations of the IGS stations used for the IGS monitoring stations tests

Full size image
Fig. 5
figure 5

Sensitivity of 3D RMS position error on process noise at 30 s update steps for GPS-only and Galileo-only static PPP solutions

Full size image

For Galileo, the process noise values for which the minima are reached correspond to 1 mm (BCE1) and 3 mm (BCE2), while for GPS, the minima are found at 3 mm (BCE1) and 10 mm (BCE2). This is in line with the fact that Galileo is characterized by a smaller overall SISRE than GPS. The results indicate that the BCE2 method is not very sensitive to process noise for both constellations. However, the pronounced dips around the minima indicate high sensitivity to the process noise for BCE1. The figure also denotes how the absence of process noise causes an important deterioration of the position accuracy.

The estimated receiver position is corrected for by the displacement due to solid Earth tides and pole tides following the IERS 2010 conventions (Petit and Luzum 2010). Furthermore, an eccentricity correction accounting for the offset between antenna reference point and marker position is applied when necessary. In case, the clock reference signals of the broadcast or precise clock product differ from that of the pseudorange observations processed in the PPP, and the corresponding differential code bias (DCB) corrections are applied. Precise clock products from the International GNSS Service (IGS; Johnston et al. 2017) typically refer to the L1 and P(Y)-code for GPS and the E1 and E5a signals for Galileo. In the case of broadcast ephemerides, the clock offsets refer to the same signals for GPS but can refer to either E1 and E5a for FNAV messages or E1 and E5b for INAV. The Galileo FNAV messages have been used exclusively for all analyses in this paper. The satellite clock offsets are furthermore corrected for periodic effects of special and general relativity, and the modeled range is corrected for the Shapiro effect (Ashby 2003).

The satellite position is corrected for the rotation of the earth during signal flight time, also known as the Sagnac effect (Enge and Mira 2006). The algorithm uses the global mapping function (GMF, Boehm et al. 2006) with atmospheric parameters based on the global pressure and temperature (GPT) model (Boehm et al. 2007). An additional tropospheric correction is estimated based on the non-hydrostatic mapping function of the GMF model. Antenna phase-center offsets and phase variations are modeled using the igs14.atx offsets and patterns (Rebischung and Schmid 2016). In accordance with the current MGEX practice, GPS L1/L2 calibrations are used in place of missing Galileo E1/E5a. In the case of PPP with precise products, the corrections are applied for transmitter and receiver antennas. For broadcast ephemerides, in contrast, no satellite antenna offset needs to be applied, since broadcast orbits are already referred to the satellites’ antenna reference point rather than the center of mass. Finally, the carrier-phase wind-up effect is modeled using satellite-type dependent attitude models for the different generations of satellites. GPS Block IIA/IIR satellites were modeled according to Kouba (2009), while Block IIF was modeled according to Dilssner (2010). Galileo was instead modeled following GSC (2019, Sect. 3.1.1.) for IOV satellites and GSC (2019, Sect. 3.1.2) for FOV satellites.

The ionospheric-free combinations of pseudorange and carrier phase observations are processed in the measurement update. An elevation-dependent weighting function is used for both observation types. An elevation cut-off angle of 10° is used. The assumed measurement standard deviations are summarized in Table 2. The individual values have been chosen based on the average measurement residuals over the set of stations used in the test case.

Table 2 Overview of observation standard deviations for Kalman filter update
Full size table

For the test performed with IGS monitoring stations, daily observation data at 30 s sampling and corresponding broadcast ephemerides received by the permanent reference station data were obtained from the IGS in the receiver independent exchange format (RINEX). For the PRE test, CODE Final MGEX precise ephemeris products (Prange et al. 2020a, b) with 5 min orbit step size and 30 s clock sampling were used. The respective clock solutions are referenced to semi-codeless L1/L2 P-code signal tracking, identified by RINEX observation codes C1W and C2W, for GPS and E1/E5a pilot tracking, identified as C1C and C1W, for Galileo. For DCBs, the daily products provided by the Chinese Academy of Sciences (CAS) were used (Wang et al. 2016). For all other tests, LNAV broadcast ephemerides of GPS and FNAV of Galileo were used. The IGS weekly global solutions for the station positions were used as a reference against which the 24 h solutions were compared.

PPP results

The models described in this study were characterized by a series of tests. The algorithm was applied to a series of different cases in terms of location, time, and conditions. Data from a pool of IGS stations over multiple days were used for the tests. These analyses are described in the following sub-section. A kinematic boat test was performed with a receiver on a motorboat on Lake Ammer, in southern Bavaria. The data were processed a-posteriori, and the results of this test are described in the second sub-section.

Tests with IGS monitoring stations

The first part of the study focuses on a series of tests that were carried out using stations and data from the IGS network, given their high availability in terms of time and geographic locations. The set of 11 stations chosen for the analyses is depicted in Fig. 4. All tests are based on a month of observations and ephemeris data covering December 2019. The same data have been used to compute 24-h solutions with the four different methods listed in Table 1. All tests were performed in both static and kinematic mode, even though the reference station antenna positions were static.

In order to characterize the capabilities of Galileo and GPS, single-constellation solutions were computed for each system. Along with those, dual-constellation GPS + Galileo solutions were also estimated. This gives the possibility to study the application of our approach to multi-GNSS positioning, along with the effects that a larger amount of tracked satellites can have on the accuracy. This is particularly relevant for Galileo, for which, in certain situations, only a reduced number of healthy satellites could be tracked simultaneously in the test period. When computing the statistics of each solution, the first 120 epochs, i.e., one hour, were removed to exclude the convergence phase of the Kalman filter. For each set of tests with similar conditions, the RMS position errors were averaged to obtain a single value that characterizes the test. Solutions that showed noticeable deviations from the main observed distribution were removed from the statistics. Tables 3 and 4 list the results of the four different methods for static and kinematic conditions, respectively. Along with the three-dimensional RMS, the horizontal two-dimensional (RM) and the vertical RMS are shown, giving the possibility to assess the two separately. The 3D RMS position error statistics of the stations are depicted in Fig. 6 for the static case and Fig. 7 for the kinematic case with a box-and-whiskers diagram. The results of the test with precise products have been omitted in the figures because their small size made a visual comparison difficult.

Table 3 Positioning performance for IGS monitoring stations in static mode
Full size table
Table 4 Positioning performance for IGS monitoring stations in kinematic mode
Full size table
Fig. 6
figure 6

Accuracy (in terms of 3D RMS position error) for the different tests in static mode, divided by station. The three horizontal lines of the colored rectangles represent the first, second and third quartile of the distribution. The vertical lines extend to the maximum and minimum value. Outliers are not plotted

Full size image
Fig. 7
figure 7

Time evolution of the position error projected on the three components East (dE), North (dN) and Up (dU) for the four processing variants on December 3, 2019 for station CHOF. The Galileo-only solution in static conditions is shown

Full size image

The tests with precise products, yielding a horizontal accuracy in the order of 1 cm for static and 3 cm for kinematic conditions, demonstrate the validity of the employed PPP models and algorithm. The values are very similar for all three system cases, with only subtle differences that are mostly due to the number of satellites available for each case.

When it comes to the three solutions with broadcast ephemerides, the results indicate that GPS has the worst performance in all cases. This confirms the expectations, given the larger SISRE for GPS. The Galileo-only and the dual-constellation solutions, on the other hand, show better accuracy and similar behavior between them. These facts already indicate that the dual-constellation case seems to be able to bring together the robustness of GPS in terms of the number of satellites and the smaller SISRE which characterizes the modern Galileo.

The two methods applying a SISRE compensation technique, BCE1 and BCE2, yield similar results. In both single-constellation solutions, the 3D RMS values are within 5% of each other and within 15% in the dual-constellation case. In every case, the improvement compared to BCE is substantial, from a minimum of 35% for Galileo up to 50% for GPS. The achievable 3D accuracy is in the order of 10 cm for Galileo and dual-constellation and 25 cm for GPS. The horizontal 2D accuracy is as small as 7 cm for Galileo and 18 cm for GPS. As a way of example, the 24-h coordinates time series obtained with the four different methods in static conditions are depicted together in Fig. 8 for a selected station. It is possible to observe how the SISREs induce a deviation from the reference position in BCE, and how this deviation is mitigated in BCE1 and BCE2. The effect is particularly noticeable in the East and Up component.

Fig. 8
figure 8

Accuracy (in terms of 3D RMS position error) for the different tests in kinematic mode, divided by station. The three horizontal lines of the colored rectangles represent the first, second and third quartile of the distribution. The vertical lines extend to the maximum and minimum value. Outliers are not plotted

Full size image

A method comparable to BCE was used in Hadas et al. (2019), obtaining an accuracy that compared to our results is up to 40% better for Galileo and up to 50% for GPS. The reason for this difference is attributed to the different filter settings and mostly to the smaller cut-off angle, which is 5° in the mentioned work and 10° in ours. However, the same results show an accuracy similar to that found with BCE1 and BCE2, suggesting that SISRE compensation techniques can counteract the adverse effects of a higher cut-off angle.

The results yield a similar picture in kinematic conditions, where all three constellation cases show similar behavior with all methods. Both BCE1 and BCE2 bring a marked improvement with respect to BCE and are characterized by similar accuracies. For the Galileo-only and the dual-constellation cases, the proposed strategies reach a 3D RMS position error within 30 cm and a horizontal one below 20 cm. GPS shows values that are roughly twice as big. In kinematic conditions, compensating the SISRE with either model improves the solution with respect to the BCE model by 50% for Galileo and dual-constellation and by 65% for GPS. Like it is the case for the static tests, the best accuracies show one order of magnitude difference compared to the standard PPP approach. For kinematic conditions as well, the 24-h coordinates time series of the different solutions are depicted in Fig. 9 for a selected station. In this case, as well, it can be observed how the amplitude of the deviation from the reference position is much bigger for BCE and for the other two approaches with broadcast ephemerides. In particular, the last few hours of the solution depict a sudden deterioration of the BCE solution, which cannot be observed for BCE1 and BCE2.

Fig. 9
figure 9

Time evolution of the position error projected on the three components East (dE), North (dN) and Up (dU) for the four processing variants on December 3, 2019 for station CHOF. The Galileo-only solution in kinematic conditions is shown

Full size image

In kinematic conditions, the difference between the results found with BCE and those from Hadas et al. (2019) increases. In particular, the BCE GPS-only solution of Table 4 if worse by a factor of four. The difference is attributed to the different filter settings, mainly to the cut-off angle. As in the static case, the improvement brought by the compensation techniques in BCE1 and BCE2 bring the values close together. In the Galileo-only solution of Table 4, the 2D RMS position error is actually almost 40% smaller than Hadas, while for GPS the same value is roughly 20% bigger.

Test with kinematic boat measurements

One of the reasons why PPP is such an attractive approach is the ability to perform absolute positioning with an accuracy that is virtually independent of the location. Since the proposed strategies aim at understanding the capabilities of PPP when precise ephemerides are not available, the studying of a real-life kinematic scenario is of particular interest to us. On September 11, 2019, a motorboat was driven on Lake Ammer, in southern Bavaria, to record GNSS observations on a moving vehicle for a period of 1 h and 20 min. Dual-frequency measurements of GPS (L1 + L2 P(Y)) and Galileo (E1/E5a) were collected with an AsteRx SB receiver connected to a Trimble Zephyr 3 Geodetic antenna. The setup is depicted in Fig. 10. In addition to the observations, the receiver recorded an RTK solution, which was later used as a reference for assessing the accuracy of the different approaches. The RTK reference antenna, located at the DLR site in Oberpfaffenhofen, was at a distance of approximately 15 km.

Fig. 10
figure 10

Picture of the Trimble Zephyr 3 Geodetic antenna used for recording the observation during the test at the Lake Ammer (southern Bavaria)

Full size image

The methods used for the boat tests are the same four listed in Table 1 and used in the previous analysis. Similar to the tests with IGS monitoring stations, one dual-constellation (Galileo + GPS) solution and one single-constellation solution for each system were computed. The tests were performed in kinematic mode. In these boat tests, the filter was updated at 1 s steps. Given the short period of the test of less than 1.5 h, this was necessary to increase the number of epochs processed. The first 5 min of data were removed from the statistics.

The results, divided in vertical, 2D horizontal and 3D RMS position errors, are listed in Table 5. Compared to the kinematic solutions of Table 4, the standard PPP solutions show a deterioration by a factor of five for all cases. A degradation is expected given the less rigorous conditions of the boat test compared to those of the monitoring stations, like strong multipath due to surrounding water.

Table 5 Accuracy results of the tests performed with the boat data
Full size table

Concerning the BCE case, the accuracy obtained in the boat tests is in line with the accuracies obtained with individual IGS monitoring stations. GPS is once more characterized by the worst 3D accuracy, in the order of 1 m in this case. Galileo yields the best accuracy of 0.5 m, while the large SISRE of GPS causes the dual-constellation solution to fall between the two single-constellation cases. The coordinates time series of the dual-constellation solutions are depicted in Fig. 11, where it can be observed how the East component causes an important deterioration of the horizontal 2D accuracy. The solution with precise orbit and clock products seems to be affected as well, unlike BCE1 and BCE2. It should be noted that the gaps in the solution reflect interruptions of the RTK reference position and not of any of the post-processed solutions.

Fig. 11
figure 11

Time evolution of the position error projected on the three components East (dE), North (dN) and Up (dU) for the four processing variants of the boat test. Dual constellation GPS + Galileo solution in kinematic conditions

Full size image

The two proposed approaches show similar accuracies in all three cases, with 3D RMS values within 10% of each other for Galileo and dual-constellation, and within 20% for GPS. Compared to BCE, the improvement brought by the SISRE compensation techniques is up to 30% for Galileo, and 60% for GPS and dual-constellation. Overall, BCE1 and BCE2 show a 2D horizontal accuracy in the order of 30 cm, 20 cm and 10 cm for GPS, Galileo, and dual-constellation, respectively.

Summary and Conclusions

Among the current GNSSs, the European Galileo system is characterized by the best broadcast ephemerides, with a typical SISRE of one to two decimeters. This accuracy suggests the possibility of performing precise point positioning (PPP) without the need for precise orbit and clock products or additional real-time corrections. Our study assesses two functional models as strategies for performing PPP with broadcast ephemerides in Kalman-filter-based algorithms. The first approach lumps unmodeled orbit and clock errors within the float ambiguities states, adding proper process noise to allow for the time variation of these errors. On the other hand, the second model introduces a dedicated SISRE parameter in the observation equations and filter states. Here, process noise is applied to the SISRE state instead of the float ambiguities.

Compensation of SISRE in the two models allows for a reduction in positioning errors by 40–70% when working with broadcast ephemerides compared to established PPP algorithms. In general, the improvement is most pronounced for constellations characterized by larger SISRE values and for positioning performed in kinematic conditions. Among the two methods for SISRE compensation considered here, the use of process noise in ambiguity states allows for a particularly simple implementation. On the other hand, the explicit incorporation of SISRE states in the estimation vector shows slightly better performance and robustness.

Overall, positioning errors at the few-decimeter level could be achieved in kinematic PPP solutions with broadcast ephemerides in our study. As expected from the very low SISRE values of Galileo broadcast ephemerides, the use of Galileo offers the best performance with horizontal errors of down to 0.2 m and 3D position errors of 0.3–0.4 m. Roughly two times larger errors were obtained in GPS-only processing of a globally distributed IGS monitoring station set. Dual-constellation solutions achieve similar accuracy as Galileo-only solutions but offer increased robustness due to the larger number of tracked satellites.

The tests demonstrate that applying broadcast ephemerides to a PPP model is a viable approach for positioning with dual-frequency code and phase observations. While not competitive with established PPP concepts based on precise ephemerides or real-time correction data, it can still offer a ten-times accuracy improvement over code-based single-point positioning (SPP). This makes PPP with broadcast ephemerides an interesting alternative for applications aiming at sub-meter accuracies such as personal navigation or traffic management. Similar to SPP, PPP with broadcast ephemerides can be performed in real-time with exclusive use of data transmitted by the GNSS satellites themselves and does not require access to external correction services.