Multi-GNSS PPP-RTK: From Large- to Small-Scale Networks

Precise point positioning (PPP) and its integer ambiguity resolution-enabled variant, PPP-RTK (real-time kinematic), can benefit enormously from the integration of multiple global navigation satellite systems (GNSS). In such a multi-GNSS landscape, the positioning convergence time is expected to be reduced considerably as compared to the one obtained by a single-GNSS setup. It is therefore the goal of the present contribution to provide numerical insights into the role taken by the multi-GNSS integration in delivering fast and high-precision positioning solutions (sub-decimeter and centimeter levels) using PPP-RTK. To that end, we employ the Curtin PPP-RTK platform and process data-sets of GPS, BeiDou Navigation Satellite System (BDS) and Galileo in stand-alone and combined forms. The data-sets are collected by various receiver types, ranging from high-end multi-frequency geodetic receivers to low-cost single-frequency mass-market receivers. The corresponding stations form a large-scale (Australia-wide) network as well as a small-scale network with inter-station distances less than 30 km. In case of the Australia-wide GPS-only ambiguity-float setup, 90% of the horizontal positioning errors (kinematic mode) are shown to become less than five centimeters after 103 min. The stated required time is reduced to 66 min for the corresponding GPS + BDS + Galieo setup. The time is further reduced to 15 min by applying single-receiver ambiguity resolution. The outcomes are supported by the positioning results of the small-scale network.


Introduction
Integer ambiguity resolution-enabled precise point positioning, PPP-RTK, enables single-receiver users to recover the integerness of their carrier-phase ambiguities, thereby reducing the positioning convergence time as compared to that of PPP [1,2]. Apart from the satellite clocks, the PPP-RTK network platform also provides users with the satellite phase biases and (optionally) the ionospheric corrections for fast integer ambiguity resolution [3,4]. Compared to the standard PPP technique, which normally requires long convergence time from tens of minutes to hours to reach cm-level accuracy [5][6][7], the resolved integer ambiguities in PPP-RTK lead to shorter convergence time in user positioning results. Using 30-s single-frequency GPS data, the ambiguities can be resolved within several minutes with the rover positioning precision reduced to mm-and cm-level in the horizontal and vertical directions [3].

GNSS-Derived Estimable Parameters Output by Curtin's PPP-RTK Platform
In this section, the GNSS measurement and dynamic models on which Curtin's PPP-RTK platform is based are briefly reviewed. We commence with the main GNSS observations underlying the stated models, i.e., the undifferenced (UD) carrier-phase and pseudo-range (code) data collected on multiple frequencies.

Undifferenced Phase and Code Observation Equations
Let the observed-minus-computed carrier-phase and code observations of receiver r (r = 1, · · · , n), tracked by satellite s (s = 1, · · · , m), be denoted by ∆φ s r,j and ∆p s r,j , respectively. The subscript j (j = 1, · · · , f ) indicates the frequency on which the observations are collected. The corresponding observation equations read [19] ∆φ s r,j = ∆ρ s r + g s r τ r − µ j ι s r + dt r − dt s + δ r,j − δ s ,j + λ j a 1s r,j , ∆p s r,j = ∆ρ s r + g s r τ r + µ j ι s r + dt r − dt s + d r,j − d s ,j , where ∆ρ s r = c sT r ∆x r denotes the 'increment' of the geometric range between receiver r and satellite s, containing the receiver's position increment ∆x r . Thus, the 3 × 1 vector c s r represents the receiver-to-satellite unit direction vector. While the hydrostatic components of the zenith tropospheric delays (ZTDs) are a priori corrected using the Saastamoinen model, the wet components, denoted by τ r , are treated as unknowns and linked to the observations through the mapping functions g s r . The first-order slant ionospheric delays, experienced on the first frequency j = 1, are denoted by ι s r . Thus, the corresponding coefficient is given as µ j = (λ j /λ 1 ) 2 , with λ j being the carrier-phase wavelength on frequency j. The receiver and satellite clock parameters are denoted as dt r and dt s , respectively. The frequency-dependent receiver and satellite code biases are denoted by d r,j and d s ,j , respectively. Likewise, δ r,j and δ s ,j , respectively, denote the receiver and satellite phase biases. The integer-valued ambiguities are denoted as a s r,j . All quantities are expressed in units of range, except the ambiguities a s r,j , which are given in cycles. The precise satellite orbits are assumed to be included in ∆φ s r,j and ∆p s r,j .

Dynamic Models: Temporal Behaviour of the Parameters
The UD formulation (1) represents a rank-deficient system of observation equations, that is, not all of the parameters involved in (1) are estimable, only combinations of them. A number of parameters, equal to the rank-deficiency, must therefore be chosen as the system's S-basis so as to form a set of minimum constraints that recovers the system of equations to one of full rank [17,18]. From an algorithmic point of view, this means that one has to remove those columns of the model's design matrix corresponding to the S-basis parameters [16]. The number of the parameters required to be chosen as S-basis is prone to be affected by any assumption placed on the underlying dynamic model as well as on the measurement model. For instance, the receivers may collect data on only one frequency or on multiple frequencies. The temporal behavior of the parameters, involved in (1), may be modelled or some (all) of the parameters may be assumed to be unlinked in time. Each of the stated scenarios changes the number of S-basis parameters. Clearly, by removing the corresponding columns of the model's design matrix, one would not estimate the remaining parameters in an unbiased manner. In other words, the interpretation of the remaining parameters does change. For instance, consider the case where all the parameters, except the clocks (and ionospheric delays), are assumed linked in time. The interpretations of the remaining parameters, referred to as the estimable parameters, are given in Table 1. In the table, the estimable parameters are distinguished from their original counterparts by the.-symbol. As shown, the estimable slant ionospheric delays at epoch i, i.e.,ι s r (i), do not represent their original versions, namely, the unbiased slant ionospheric delays ι s r (i). Instead, they represent biased slant ionospheric delays that are lumped with the (scaled) differential code biases (DCBs) d r,GF (1) and d s ,GF (1) at the first epoch i = 1. In this contribution, we present results under the following two different dynamic models: Table 1. Estimable GNSS (Global Navigation Satellite System) network-derived parameters formed by an S-basis [16]. The additional parameters ∆x 1 (1) and τ 1 (1) (within {.}) are taken as S-basis for the small-to-regional scale networks, i.e., when c s r ≈ c s 1 and g s r ≈ g s 1 (r = 2, . . . , n). The argument (i) stands for the epoch index. • Time-unlinked clocks: While both the receiver and satellite clocks are treated unlinked in time, the temporal behaviour of the remaining parameters (except the ambiguities) are modeled as a 'random-walk' process. The ambiguities are assumed constant in time, unless slips occur (cf. Section 3.1).

•
Time-linked satellite clocks: While the receiver clocks are assumed unlinked in time, the temporal behaviour of the satellite clocks is modeled by a constant-velocity dynamic model (i.e., the 2-state clock-plus-drift model) [15] (cf. Section 3.3).

Measurement Models
Next to the dynamic models that concern the temporal behaviour of the parameters, any assumption underlying the measurement model can also affect the rank-deficiencies involved in Equation (1). For instance, if one would parametrize the slant ionospheric delays ι s r into a fewer number of unknown parameters using, e.g., the single-layer models [20], the interpretation of the estimable parameters in Table 1 changes. The following measurement models are considered: -From small-to large-scale networks. Consider the case where the inter-station distances between the network receivers are so short such that all receivers view satellite s from almost the same direction angle. Thus, c s r ≈ c s 1 and g s r ≈ g s 1 (r = 2, . . . , n). As a consequence, the GNSS data cannot distinguish between the tropospheric delays g s 1 τ 1 (or the geometric delays ∆ρ s 1 ) and the delays caused by the satellite clocks dt s in Equation (1). To tackle such near rank-deficiency between the ZTDs τ r (position increments ∆x r ) and the satellite clocks dt s , one has to therefore choose the reference ZTD τ 1 (1) (reference position ∆x 1 (1)) as an additional S-basis parameter [16]. As shown in Table 1, the estimable ZTDsτ r do not represent the unbiased ZTDs τ r , but ZTDs relative to the reference ZTD τ 1 (1) at the first epoch i = 1. Similarly, the estimable satellite clocks dt s are further biased by the terms g s 1 τ 1 (1) and c sT 1 ∆x 1 (1). In case of large-scale networks, however, the stated rank-deficiency is absent. Numerical results concerning both the large-scale and small-scale networks are presented in Sections 3 and 4, respectively.
-Fully-and partially-known receivers' positions. In the system of observation Equation (1), the positions of the network receivers are assumed unknown. Note, however, that the positions of some (or all) of the network receivers can be a priori known. For instance, in case of continuously operating reference station (CORS) networks, the positions of all the network receivers are assumed known. Next to the CORS network setup, we also consider the case where the positions of some of the receivers are fully unknown. This is particularly useful when the network is extended by newly-established stations whose positions are required to be determined.
-Ionosphere-float and -weighted models. Due to the spatial correlation of the ionosphere, the slant ionospheric delays experienced by nearby receivers from a given satellite are almost identical. One can therefore make use of this property and augment the observation Equation (1) with the following extra observation Equation [21]: in which dι s r denote the ionospheric pseudo-observations taking zero samples values. The smaller the distance between receiver r and r = 1, the better the approximation 0 ≈ (ι s r − ι s 1 ) becomes. Therefore, one can assign weights to the pseudo-observations dι s r such that the weights increase when the inter-station distances decrease. Examples of such weights are presented in, e.g., (Appendix A, [15]). When the GNSS observation Equation (1) is augmented with Equation (2), we refer to the model as the ionosphere-weighted model. In the absence of the extra observation Equation (2), the underlying model is referred to as the ionosphere-float model. Both of these models are considered in the following.

PPP-RTK Results: A Large-Scale Network
This section presents PPP-RTK network and user results based on processing of multi-GNSS data from Australia-wide network and user stations tracking at least GPS, BDS, and Galileo satellites. As depicted in Figure 1, the network is formed by 22 multi-GNSS continuous operating reference stations consisting of various receiver types including Trimble NetR9 (Sunnyvale, CA, USA), Septentrio PolRx5 (Leuven, Belgium), Septentrio PolaRx4TR, Septentrio PolaRx4TR, Septentrio PolaRxS, and Leica GR30 (St. Gallen, Switzerland). Results discussed in the following sections are based on data collected from these network stations and 50 user stations across Australia for ten days across three months (every seventh day) starting from day of year 340 of year 2016 with 30 s sampling interval. Station 1 (UCLA, Eucla in Western Australia) is arbitrarily selected as the reference station for PPP-RTK network processing.

Network Results
Multi-GNSS network data was processed with Curtin's PPP-RTK Network platform. Satellite positions were determined using precise IGS (International GNSS Service) Multi-GNSS Experiment (IGS-MGEX) orbits [22][23][24] while station coordinates were extracted from Geoscience Australia's Solution Independent Exchange format (SINEX) files. The estimable satellite clocks of multi-GNSS observables were aligned to respective reference observables using IGS-MGEX satellite DCB files. Stochastic, dynamic, and ambiguity resolution parameter settings for network processing are summarized in Table 2. Figure 2 depicts satellite visibility and Local Overall Model (LOM) test statistics [25]. The top panel indicates the time period(s) for which any given satellite of GPS, BDS, and Galileo is in common view of at least two stations. The middle panel shows the number of satellites visible above 10 degrees for individual systems as well as combined systems. As of the time of data collection, the number of satellites of combined (GPS, BDS, and Galileo) systems varies from 23 to 32. The bottom panel shows the time-series of LOM test statistics, along with their critical values. These LOM test values are obtained in the Detection step, which is the first step of the Main Detection-Identification-Adaptation (DIA) procedure. The purpose of the DIA procedure, implemented in Curtin PPP-RTK platform, is to identify mis-modeled effects (e.g., phase slips and code outliers) and adapt the model accordingly. After identifying mis-modeled effects and successful model's adaptation, all LOM values are found to be less than their critical values during the day, indicating the validity of the network model used. Figure 3 shows the time series of the satellite clock estimates (δt s ), together with their formal standard deviations for an arbitrary satellite per system. It usually takes about 2 h to have satellite clock estimates with formal standard deviations lower than 2 cm. It can be seen that, for each satellite shown, its clock estimates behave almost linearly as a function of time, which is explored for satellite clock modeling in Section 3.3. Figures 4 and 5

User Results
With availability of network-derived satellite clock and satellite phase bias products from network processing above, the convergence behavior of user position is analyzed next using multi-GNSS data from 50 CORS receivers evenly distributed across Australia. The associated parameter settings for user processing were identical to that of network processing above, except that user position was assumed to be unknown and unlinked in time. To get enough samples for generating reliable convergence curves, user processing was repeated every hour and across ten days resulting in 3800 data windows. The earliest user processing window each day started an hour after the initialization of respective network processing allowing network derived products to converge. The user processing was repeated for four different GNSS combinations: GPS-only, GPS + BDS, GPS + Galileo, GPS + BDS + Galileo.

Impact of Latency and Satellite Clock Modelling
With the dynamic satellite clock model incorporated in the network Kalman filter (Table 6), the satellite clock rate dṫ s is able to be estimated at each epoch [15]. In case of latency of the network products, denoted by ∆i, the estimable satellite phase biases and satellite clocks can be time-predicted with dť s (i + ∆i) = dt s (i) + (∆i) dˆṫ s (i),  Table 2.

Dynamic Model
Satellite clocks q c = 0.001 m/ √ s In this study, the satellite clocks are modeled with white frequency noise (WFN) with the process noise set to be 0.001 m/ √ s [15]. With 1 Hz GPS dual-frequency data of the Australia wide network processed from 13:00:00 to 14:59:59 in GPS Time (GPST) on 12 October 2017, the network products were estimated without latency and predicted with latencies of 10 and 15 s. The starting time of the user processing varies from 13:40:00 to 14:00:00 with a time shift of 1 min and continues for 1 h for each station. Using 51 user stations, more than 1000 data samples with 1 h coordinate time series were generated and used for computing the convergence curves for different latencies in ambiguity-fixed and -float cases. Figure 9 shows the 75% convergence curves of the horizontal radial components and the absolute up components of the user coordinates deviated from the ground truth. The network products were estimated without latency (see the green lines in Figure 9) and predicted with latencies of 10 and 15 s (see the blue and magenta lines in Figure 9). The red lines mark the y-values of 5 cm and 1 dm. We see that the convergence time to 1 dm increases with the increasing latencies. For ambiguity-float case, more time is required to reduce the coordinate increments to 1 dm compared to the ambiguity-fixed case. With the ambiguities fixed and the network products estimated without latency, it takes around 10 min to let the horizontal and up coordinate increments reduce to below 1 dm in 75% of the cases (see also the green lines in Figure 9a,b). As the latency of the network products increases, more time is needed. With a latency of 15 s (see also the magenta lines in Figure 9a,b), it takes around 15 min to reduce the coordinate increments to below 1 dm. In ambiguity-float case, more time is required for the coordinate convergence. As shown by the green lines in Figure 9c

PPP-RTK Results: A Small-Scale Network
This section presents user positioning performance results of a single user (UWA0) in a small-scale network (30 km) formed by three stations (shown in green in Figure 10). All three network stations are equipped with Trimble NetR9 receivers, while a user station is equipped with a Septentrio PolRx5 receiver (Leuven, Belgium). Results discussed in this section are based on data collected from these four stations for ten days across three months (every seventh day) starting from day of year 190 of year 2017 with 30 s sampling interval. CUT0, which is Curtin University's IGS reference station, is arbitrarily selected as the reference station for PPP-RTK network processing. For this small scale network, satellite positions were determined using broadcast navigation messages. Ground truths of station coordinates were computed using Geoscience Australia's AUSPOS service. The estimable satellite clocks of multi-GNSS observables were aligned to respective reference observables using IGS-MGEX satellite DCB files. For this small network, differential ionosphere delays were weighted as informed in Table 7 to strengthen the underlying model, while other parameter settings are as in Table 2. In addition to processing PPP-RTK user mode, user position was also determined using an in-the-loop mode, whereby the user position is assumed to be unknown and unlinked in time, is estimated together with other network parameters. The associated parameter settings for PPP-RTK user processing were identical to that of network processing above, except that user position was assumed to be unknown and unlinked in time. To get enough samples for generating reliable convergence curves, PPP-RTK user processing and in-the-loop user positioning were repeated every 40 min and across ten days resulting in 340 data windows. Unlike processing in previous sections, network and user processing started at the same epoch for a fair comparison between PPP-RTK user and in-the-loop user positioning performance. Table 7. Parameter settings for small network processing; Other parameter settings are as in Table 2.

Ionosphere spatial model
Ionosphere weighted Applicable inter-station distance (l 0 ) 2 km Standard deviation of undifference ionosphere observables 0.01 m/l 0 Figures 11 and 12 depict user convergence curves for GPS-only and triple system (GPS + BDS + Galileo) processing comparing performance of both PPP-RTK user and in-the-loop user mode processing, while Tables 8 and 9 summarize corresponding convergence time for 1 dm and 5 cm (in brackets). Slightly better performance of in-the-loop user especially during the initialization is due to that in-the-loop user enjoys the use of full (variance-covariance) information, while a PPP-RTK user uses network products as deterministic corrections even at the initialization period. Hence, it is recommended to allow sometime (e.g., an hour) for network processing to converge before using its products. By using multi-GNSS, one can achieve almost instantaneous ambiguity resolved precise user position.

PPP-RTK Results: Low-Cost Single-Frequency Receivers
This section presents low-cost user positioning performance results of a single user (SPR1, shown in red in Figure 13) using single station (UWAU, shown in green in Figure 13). In this analysis, a real data campaign was set up with the reference station on top the Physics Building of University of Western Australia and user on top of Building 207 of Curtin University. Both stations are equipped with U-Blox M8 multi-GNSS receivers connected to U-Blox multi-GNSS patch antennas. Results discussed in the this section are based on multi-GNSS, single-frequency data collected from these two stations for three consecutive days (day of year 300, 301, and 302 of year 2017) with 30 s sampling interval. For this small scale network, satellite positions were determined using broadcast navigation messages. Ground truths of station coordinates were computed using batch processing of data form these stations and a nearby IGS reference station CUT0 (see Section 4). The estimable satellite clocks of multi-GNSS observables were aligned to respective reference observables using IGS-MGEX satellite DCB files. Except parameter settings listed in Table 10, the parameter settings are as in Table 2. Similar to Section 4, user position was determined using both PPP-RTK user and in-the-loop-user modes, where differential ionosphere delays were weighted as informed in Table 10 strengthening underlying model. Large code noise standard deviation reflects the quality of code observations from low cost receivers connected to patch antennas. To get enough samples for generating reliable convergence curves, PPP-RTK user processing and in-the-loop user positioning were repeated every 10 min and across three days resulting in 356 data windows. Similar to Section 4, network and user processing started at the same epoch for a fair comparison between PPP-RTK user and in-the-loop user positioning performance. Table 10. Parameter settings for low-cost network processing; other parameters are as in Table 2.

Ionosphere
Ionosphere spatial model Ionosphere weighted Applicable inter-station distance (l 0 ) 2 km Standard deviation of undifference ionosphere observables 0.1 m/l 0 Figure 14 depicts user convergence curves triple system (GPS+BDS+Galileo), single-frequency processing comparing performance of both the PPP-RTK user and in-the-loop user mode processing, while Table 11 summarizes corresponding convergence time for 1 dm and 5 cm (in brackets). Better performance of in-the-loop user especially during the initialization is due to that in-the-loop user enjoys the use of full (variance-covariance) information, while the PPP-RTK user uses network products as deterministic corrections even at the initialization period.

Conclusions
In this contribution, we provided numerical insights into the role taken by the multi-GNSS integration in delivering fast and high-precision positioning solutions using PPP-RTK. With the aid of Curtin PPP-RTK platform, data-sets of GPS, BDS and Galileo were processed in stand-alone and combined forms. Given the ground-truth coordinates of several single-receiver users, a large number of samples of the positioning errors (∼3800 samples) were collected so as to compute representative positioning convergence curves (cf. Figures 6-9). In case of the Australia-wide GPS-only ambiguity-float setup, 90% of the horizontal positioning errors (kinematic mode) were shown to become less than five centimeters after 103 min. The stated required time is reduced to 66 min for the corresponding GPS + BDS + Galieo setup. The time is further reduced to 15 min by applying single-receiver ambiguity resolution. We also showed the impact of the latency in sending the network-derived corrections on the user positioning performance (cf. Figure 9). With a latency of 15 s, it takes around 15 min to have multi-GNSS ambiguity-fixed positioning errors less than 1 dm. For the corresponding ambiguity-float case, more time is required. In that case, the required time increases to 45 min.
We also presented multi-GNSS PPP-RTK results obtained by single-frequency low-cost receivers for which a 'small-scale' network is considered. The PPP-RTK user results were compared with the so-called 'in-the-loop' user, that is, the user's data are simultaneously processed together with those of the network. While the PPP-RTK ambiguity-fixed positioning errors become less than 1 dm after 14 min, the corresponding in-the-loop counterparts only require 2 min to reach 1 dm.