GPS + Galileo + BDS-3 medium to long-range single-baseline RTK: an alternative for network-based RTK?

Abstract

Thanks to the development of the real-time kinematic (RTK) algorithm and the emerging Global Navigation Satellite System (GNSS), especially for Galileo and BeiDou-3, reliable positioning accuracy for medium and long-baseline RTK became possible globally. Moreover, with the development of the GNSS receiver hardware, baseline length limitations due to radio-based communications are removed thanks to internet-based communication. In this work, single-baseline RTK, incorporated partial ambiguity resolution with troposphere and ionosphere weighting, using GPS (G), Galileo (E), BeiDou-3 (C3) and multi-GNSS (GE and GEC3), is conducted with real GNSS data of EUREF Permanent GNSS network under three different cutoff angles (10°, 20°, and 30°) for six different lengths of baselines (~50, ~150, ~250, ~350, ~450, and ~550 km). The results show that the multi-GNSS RTK solution significantly contributed to the positioning accuracy and convergence time of the single-system RTK solutions. Based on the results, non-available epoch-wise solutions for the high-degree cutoff angles are more obvious for the single-system RTK, whereas multi-GNSS solutions provide 100% solutions for each cutoff angle and baseline. The results indicate that instantaneous and a few epochs single-epoch ambiguity resolution is feasible for 50, 150, 250 and 350 km baseline lengths for multi-GNSS RTK. Based on the positioning results, horizontal–vertical positioning improvements of multi-GNSS RTK (GEC3) compared with the single-system GPS RTK are found as 50%–37%, 40%–35%, 55%–47%, 53%–54%, 57%–49% and 57%–49% for 50, 150, 250, 350, 450 and 550 km, respectively, under a 10° cutoff angle. For 20° and 30° cutoff angles, the accuracy improvements are much higher. The convergence time improvements (n/e/u) of multi-GNSS RTK (GEC3) compared with the single-system GPS RTK are found as 86/92/75%, 77/67/72%, 75/77/83%, 53/56/52%, 69/49/62%, and 52/45/39% for 50, 150, 250, 350, 450 and 550 km, respectively, under a 10° cutoff angle.

1. Introduction

Real-time kinematic (RTK) positioning is one of the most frequently used positioning techniques for achieving high-accuracy positioning from mm to cm level. With the development of the new Global Navigation Satellite System (GNSS) receiver hardware, data transmission from base to rover does not rely solely on the radio-based communications that significantly restrict the baseline length between rover and base sites. Thanks to internet communications, baseline restrictions and other disadvantages of radio-based communications were eliminated in RTK (Jones and Kelly, 2007). However, differential ionospheric and tropospheric delays that are not present in the short baseline are major limiting factors of the reliable ambiguity resolution (AR) and positioning accuracy for medium and long-baseline RTK (Li et al., 2010; Takasu and Yasuda, 2010).

To overcome this limitation, relative ionospheric and tropospheric biases must be considered in the functional and stochastic models of RTK. Since ionospheric effects are frequency dependent, the first-order ionospheric errors can be eliminated by using a dual-frequency ionospheric-free (IF) linear combination as a conventional method (Chen et al., 2022). However, fixing the integer ambiguities using the IF combination is difficult, and the noise amplification is higher than in the uncombined model. As a result, slow convergence time and loss of positioning accuracy are more evident in the IF combination for medium and long baselines (Takasu and Yasuda, 2010). The other approach is to estimate the relative ionospheric delay parameters for each double-differenced (DD) satellite pair using uncombined phase and code observations (Zhang et al., 2020b). This approach can be supported by wide-lane (WL) and narrow-lane (NL) ambiguity estimation to improve the AR reliability (Li et al., 2014).

Over the last several years, GNSS has been extensively expanded with the completion of Galileo and BDS-3 full operational capability (FOC). The multi-GNSS RTK solution became more reliable than a single-system RTK in terms of ambiguity resolution, convergence time and positioning accuracy, especially for medium and long baselines (Zhang et al., 2020b; Mirmohammadian et al., 2022; Xia et al., 2022). Some studies paid attention to the performance of multi-GNSS RTK. Zhao et al. (2014) investigated GPS + BDS-2 RTK single-frequency for short baselines. They found that time to first fix (TTFF) and ambiguity fixing rate is improved compared with GPS-only RTK. Odolinski et al. (2015) conducted GPS only and GPS + BDS-2 RTK for short and long baselines. They found that adding BDS-2 to GPS improves the TTFF, positioning accuracy and ambiguity dilution of precision (ADOP), especially for the high degree of cutoff angles. Zhang et al. (2020b) examined the contribution of QZSS to GPS-only and GPS + Galileo RTK for long baselines in Japan. They found that adding QZSS to GPS-only RTK and GPS + Galileo RTK improves positioning accuracy, ambiguity fixing rate and convergence time for long baselines in Japan. Shu et al. (2018) conducted GPS-only and GPS + BDS-2 RTK for ~60 km baselines using the troposphere constraint method. They found that adding BDS-2 to GPS-only RTK improves positioning accuracy and ambiguity fixing rate. Odijk et al. (2014) evaluated the Galileo contribution to GPS-only RTK for 250 km baseline using simulated GNSS data. According to the results, TTFF obtained from GPS + Galileo RTK using partial ambiguity resolution requires a few minutes compared to 15–25 min in GPS-only RTK using dual-frequency data. Liu et al. (2022) analysed the performance of single-frequency and dual-frequency RTK using GPS-only, BDS-2-only, BDS-3-only and BDS-2 + BDS-3 data for short baselines in China. In addition to the legacy B1I/B3I signals, new B1C/B2a signals were also used for BDS-3. The results show that BDS-3 RTK using B1C/B2a dual-frequency signals has the best ambiguity fixing rate and positioning accuracy among the other combinations. For the single-frequency RTK, the results indicate that BDS-2 + BDS-3 RTK produced the best results in terms of the ambiguity fixing rate and positioning accuracy. Wang et al. (2022) investigated BDS-2 + BDS-3 and GPS dual-frequency RTK performance for five different short to medium-length baselines (from 6 to 128 km). The results showed that BDS-2 + BDS-3 and GPS RTK performances are comparable in short baselines, but BDS-2 + BDS-3 RTK has slightly better performance than GPS in medium baselines. The results also indicate that single-frequency BDS-3 RTK using B1C, B1I and B3I produced similar accuracy, while B2a single-frequency BDS-3 RTK produced poorer accuracy.

Thanks to the completion of BDS-3 and Galileo FOC, GPS + Galileo + BDS-3 (G + E + C3) multi-GNSS RTK became more feasible due to the strength of the model and redundancy, especially for medium and long baselines, and it can be considered as an alternative technique to network-based RTK (Zhang et al., 2020a). Since GPS, Galileo and BDS-3 employ the code division multiple access (CDMA) technique, DD AR can be performed without the necessity of identical receivers between base and rover and the interfrequency biases (IFB) calibration file, as in the case of GLONASS RTK (Zaminpardaz et al., 2021). No literature has thoroughly addressed the evaluation of G + E + C3 RTK performance for medium and long baselines. Moreover, fully exploiting the performance of G + E + C3 RTK is critical for the feasibility analysis of an alternative option on behalf of network-based RTK applications. Therefore, this study assesses the performance of G + E + C3 RTK on medium and long baselines in terms of ambiguity fixing rate and positioning accuracy to fully exploit the potential of G + E + C3 RTK. The remainder of this paper is presented as follows. Section 2 gives the functional and stochastic models of G + E + C3 RTK . Section 3 presents the experimental setup. Section 4 gives the results and analyses, and the study is concluded in Section 5.

2. Multi-GNSS RTK functional and stochastic model

A medium to long baseline DD functional model for pseudo-range and carrier-phase can be expressed as (Shu et al., 2018; Liu et al., 2022):

(1)\begin{gather}\nabla \Delta \emptyset _{rb,f}^{ij} = \nabla \Delta \rho _{rb,f}^{ij} + {\lambda _f}\nabla \Delta N_{rb,f}^{ij} - {\gamma ^2}\nabla \Delta I_{rb,f}^{ij} + \nabla \Delta T_{rb}^{ij} + \nabla \Delta \varepsilon _{rb,f}^{ij}\end{gather}

(2)\begin{gather}\nabla \Delta P_{rb,f}^{ij} = \nabla \Delta \rho _{rb,f}^{ij} + {\gamma ^2}\nabla \Delta I_{rb,f}^{ij} + \nabla \Delta T_{rb}^{ij} + \nabla \Delta \varepsilon _{rb,f}^{ij}\end{gather}

where $\nabla \Delta$ is the double-differenced operator; $r$ and b represent the rover and base station; i and j are the reference satellite and the other satellite; $\emptyset$ and P are the carrier phase and pseudo-range observations; $\rho$ is the geometric range between the phase centre of satellite and receiver; f, ${\lambda _f}$ and N are the frequency channel, wavelength and the integer ambiguity on the corresponding frequency channel; ${\gamma ^2}$ is the frequency-depended ionospheric coefficient $({f_1^2/f_f^2} )$; $I$and T are the slant ionospheric and tropospheric delay; and $\varepsilon$ is the DD residuals, including observation noise, orbital error and multipath. The zenith DD hydrostatic delay can be modelled using empirical global tropospheric models. In this study, meteorological parameters were derived from the GPT2_5w model (Liu et al., 2017), and the zenith hydrostatic delay is estimated from the Saastamoinen model (Saastamoinen, 1973). The Vienna mapping function 1 (Boehm et al., 2006) was used for the conversion of slant tropospheric delay to zenith tropospheric delay. The zenith DD wet delay can be estimated as an unknown parameter instead of independent line-of-sight delays and it can be expressed using the mapping function ($MF$) as:

(3)\begin{equation}\nabla \Delta ZWD = (MFw_r^i(\varepsilon )ZW{D_r} - MFw_b^i(\varepsilon )ZW{D_b}) - (MFw_r^j(\varepsilon )ZW{D_r} - MFw_b^j(\varepsilon )ZW{D_b})\end{equation}

where $\varepsilon$ represents the elevation angle of the satellite. Generally, double-differenced tropospheric and ionospheric delays can be ignored for short baseline RTK (<10 km), and this model is the so-called ionosphere- and troposphere-fixed model. For medium and long baselines, double-differenced tropospheric and ionospheric delays cannot be cancelled, and they should be estimated as unknown parameters, along with the position of the rover and DD ambiguities. When the ionospheric slant delay is incorporated with the appropriate constraint, it is referred to as the ionosphere-weighting model (Odolinski et al., 2015). Tropospheric delay is mainly related to the length of baseline and height differences between the rover and base stations, while the ionospheric delay is additionally related to the latitude of the base and rover stations along with the baseline length. In this study, the constraints and random walk noises of zenith troposphere and ionosphere parameters are taken as (Zhang et al., 2020b):

(4)\begin{gather}{P_{\textrm{trop}({k - 1} )}} = \log (1 + D \,{\cdot}\, 5e - 4) \,{\cdot}\, 0 \,{\cdot}\, 05 + H \,{\cdot}\, 5e - 5\end{gather}

(5)\begin{gather}{P_{\textrm{iono}({k - 1} )}} = \left( {D \,{\cdot}\, 5e - 6 \,{\cdot}\, \exp \left( {\frac{{90 - \varphi }}{{50}} - 1} \right)} \right){/}\textrm{sin}(\varepsilon )\end{gather}

(6)\begin{gather}{Q_{\textrm{trop}}} = (\textrm{log}(1 + D \,{\cdot}\, e - 5) \,{\cdot}\, 0 \,{\cdot}\, 02 + H \,{\cdot}\, e - 5) / \sqrt {3600 \,{\cdot}\, \Delta t} \end{gather}

(7)\begin{gather}{Q_{\textrm{iono}}} = \left( {D \,{\cdot}\, 5e - 6 \,{\cdot}\, \exp \left( {\frac{{90 - \varphi }}{{50}} - 1} \right)} \right){/}\textrm{sin}(\varepsilon ) \,{\cdot}\, \sqrt {3600 \,{\cdot}\, \Delta t} \end{gather}

where ${P_{\textrm{trop}({k - 1} )}}$, ${Q_{\textrm{trop}}}$, ${P_{\textrm{iono}({k - 1} )}}$ and ${Q_{\textrm{iono}}}$ are the constraints and random walk noises of troposphere and ionosphere; $(k - 1)$ indicates the previous corresponding epoch; $D$ and H are the inter-station distance and height differences; $\varphi$ is the average latitude of rover and base stations; and $\varepsilon$ is the satellite elevation. As an example, for a baseline of 350 km with an interstation height difference of 308 m (BSCN-EGLT), the constraint and random walk noises of troposphere/ionosphere are set as 0 ⋅ 27 m/0 ⋅ 033 m/$\sqrt h$ and 1 ⋅ 9 m/1 ⋅ 9 m/$\sqrt h$ for 20° cutoff angle. This zenith ionosphere maps to 5 ⋅ 5 m and 5 ⋅ 5 m/$\sqrt h$ of slant ionosphere using the single-layer model (Schaer, 1999).

Design matrix ($H$) and vector of observed minus computed ($b$) based on linearly independent DD observations can be expressed as:

(8)\begin{gather}H({m \times n} )= \left[ {\begin{array}{*{20}{l}} {\dfrac{{\partial ({\nabla \Delta \emptyset_{rb,f}^{ij}} )}}{{\partial ({{X_r},{Y_r},{Z_r}} )}}}& {\dfrac{{\partial ({\nabla \Delta \emptyset_{rb,f}^{ij}} )}}{{\partial ({{\lambda_f}\nabla \Delta N_{rb,f}^{ij}} )}}}& {\dfrac{{\partial ({\nabla \Delta \emptyset_{rb,f}^{ij}} )}}{{\partial ({\nabla \Delta ZWD} )}}}& {\dfrac{{\partial ({\nabla \Delta \emptyset_{rb,f}^{ij}} )}}{{\partial ({{\gamma^2}\nabla \Delta I_{rb,f}^{ij}} )}}}\\ {\dfrac{{\partial ({\nabla \Delta P_{rb,f}^{ij}} )}}{{\partial ({{X_r},{Y_r},{Z_r}} )}}}& {\dfrac{{\partial ({\nabla \Delta P_{rb,f}^{ij}} )}}{{\partial ({{\lambda_f}\nabla \Delta N_{rb,f}^{ij}} )}}}& {\dfrac{{\partial ({\nabla \Delta P_{rb,f}^{ij}} )}}{{\partial ({\nabla \Delta ZWD} )}}}& {\dfrac{{\partial ({\nabla \Delta P_{rb,f}^{ij}} )}}{{\partial ({{\gamma^2}\nabla \Delta I_{rb,f}^{ij}} )}}} \end{array}} \right]\end{gather}

(9)\begin{gather}b({m \times 1} )= \begin{array}{*{20}{c}} {\nabla \Delta \emptyset _{rb,f}^{ij} - \nabla \Delta \rho _{rb,f}^{ij}}\\ {\nabla \Delta P_{rb,f}^{ij} - \nabla \Delta \rho _{rb,f}^{ij}} \end{array}\end{gather}

where ${X_r},{Y_r},{Z_r}$ are the rover's 3D coordinates and $\partial$ represents the partial derivate with respect to the unknown parameters. Assuming the base and rover observe four satellites for the corresponding epoch, the number of linearly independent DD observations equals 12 after choosing the reference satellite for the dual frequency phase and code observations. The epoch-wise reference (pivot) satellite was chosen for the satellite with the maximum elevation angle. The minimum number of common satellites should be five (equals 16 linearly independent DD observations for dual frequency) to estimate the receivers’ 3D position (xyz) along with the eight DD ambiguity (for dual frequency), one zenith DD wet tropospheric delay, and four slant ionosphere parameters for medium and long single-baseline RTK. As a result, the dimensions of H and b equal $({16 \times 16} )$ and$({16 \times 1} )$, respectively, for the five common satellites between the rover and base stations. However, for the epochs that include four common satellites, baseline processing can still be conducted without explicitly estimating the ionospheric and/or tropospheric parameters using the a priori information of the ionosphere and troposphere parameters (i.e. the initial value and its covariance). In this study, the minimum number of common satellites for baseline processing was set to four. A prior information constraint is used without epoch-wise estimation for the zenith DD wet tropospheric delay when the minimum number of common satellites is four. Variance-covariance ($P$) and process noise ($Q$) matrixes of the estimated parameters can be expressed as:

(10)\begin{gather}P({n \times n} )= \begin{array}{*{20}{c}} {\sigma _x^2}& .& .& .& .& .& .& .\\ .& {\sigma _y^2}& .& .& .& .& .& .\\ .& .& {\sigma _z^2}& .& .& .& .& .\\ .& .& .& {{P_{\textrm{iono}({k - 1} )}}}& .& .& .& .\\ .& .& .& .& {{P_{\textrm{trop}({k - 1} )}}}& .& .& .\\ .& .& .& .& .& {\sigma ({{\lambda_f}\nabla \Delta N_{rb,f}^{i1}} )}& .& .\\ .& .& .& .& .& .& \ddots & .\\ .& .& .& .& .& .& .& {\sigma ({{\lambda_f}\nabla \Delta N_{rb,f}^{ik}} )} \end{array}\end{gather}

(11)\begin{gather}Q({n \times n} )= \begin{array}{*{20}{c}} {{Q_x}}& .& .& .& .& .& .& .\\ .& {{Q_y}}& .& .& .& .& .& .\\ .& .& {{Q_z}}& .& .& .& .& .\\ .& .& .& {{Q_{iono}}}& .& .& .& .\\ .& .& .& .& {{Q_{trop}}}& .& .& .\\ .& .& .& .& .& 0& .& .\\ .& .& .& .& .& .& \ddots & .\\ .& .& .& .& .& .& .& 0 \end{array}\end{gather}

where ${Q_x},$ ${Q_y}$, and ${Q_z}$ are the a priori variances of the receiver's 3D position, and they were set to 3 m. ${P_{\textrm{iono}({k - 1} )}}$ and ${P_{\textrm{trop}({k - 1} )}}$ are the constraints of the ionosphere and troposphere which are computed from Equations (4) and (5), ${Q_{\textrm{iono}}}$ and ${Q_{\textrm{trop}}}$ are the process noises of ionosphere and troposphere which are computed from Equations (6) and (7), $\sigma ({{\lambda_f}\nabla \Delta N_{rb,f}^{i1}} )$ to $\sigma ({{\lambda_f}\nabla \Delta N_{rb,f}^{ik}} )$ are the a priori variance of DD ambiguities using i reference satellite and other satellites (from $1\textrm{st}$ to $k\textrm{th}$ satellites) and they were computed as the difference between pseudorange and phase observations on the corresponding signals. The dots in the matrixes represent 0. The initial position of the receiver was taken from the RINEX header. If the computed position and the initial position differ 100 m from each other, single-point positioning was conducted as an initial position for the corresponding epoch. The dimension of P and Q can be expressed as:

(12)\begin{equation}n = 5 + [\kern-0.15em[{\nabla \Delta \emptyset_{rb,f}^{ij}} ]\kern-0.15em]\end{equation}

where $[\kern-0.15em[{\nabla \Delta \emptyset_{rb,f}^{ij}} ]\kern-0.15em]$ denotes the number of linearly independent DD ambiguities. For four common satellites between the base and rover, n equals 8.

An elevation dependent stochastic model for variance–covariance ($R$) matrix of the phase and code measurements can be expressed as (Zhang et al., 2020b):

(13)\begin{gather}{R_k}({m \times m} )= \begin{array}{*{20}{c}} {{\sigma _{\emptyset _{rb,f}^{ij}}}}& .\\ .& {{\sigma _{P_{rb,f}^{ij}}}} \end{array}\end{gather}

(14)\begin{gather}{\sigma _\emptyset } = 0 \,{\cdot}\, 003\ m \,{\cdot}\, \left( {0 \,{\cdot}\, 5 + \frac{{0 \,{\cdot}\, 5}}{{\sin (\varepsilon )}}} \right)\end{gather}

(15)\begin{gather}{\sigma _P} = 0 \,{\cdot}\, 3\ m \,{\cdot}\, \left( {0 \,{\cdot}\, \le 5 + \frac{{0 \,{\cdot}\, 5}}{{\sin (\varepsilon )}}} \right)\end{gather}

State transition matrix ($A$) between previous $(k - 1)$ and current ($k$) epochs can be described as:

(16)\begin{equation}A({n \times n} )= \begin{array}{*{20}{c}} 1& .& .& .& .& .& .& .\\ .& 1& .& .& .& .& .& .\\ .& .& 1& .& .& .& .& .\\ .& .& .& 1& .& .& .& .\\ .& .& .& .& 1& .& .& .\\ .& .& .& .& .& 1& .& .\\ .& .& .& .& .& .& \ddots & .\\ .& .& .& .& .& .& .& 1 \end{array}\end{equation}

In this study, the receiver dynamic model was not considered. Epoch-wise conventional Kalman filtering (Kim, 2011) was applied to estimate the unknown parameters as:

(17)\begin{gather}\hat{x}_k^ - = A \,{\cdot}\, {\hat{x}_{k - 1}}\end{gather}

(18)\begin{gather}P_k^ - = A \,{\cdot}\, {P_{k - 1}} \,{\cdot}\, {A^T} + Q\end{gather}

(19)\begin{gather}{K_k} = P_k^ -\,{\cdot }\,{H^T} \,{\cdot}\, {({H \,{\cdot}\, P_k^ -\,{\cdot }\,{H^T} + R} )^{ - 1}}\end{gather}

(20)\begin{gather}{\hat{x}_k} = \hat{x}_k^ -{+} {K_k} \,{\cdot}\, ({{b_k} - H \,{\cdot}\, \hat{x}_k^ - } )\end{gather}

(21)\begin{gather}{P_k} = P_k^ - - {K_k} \,{\cdot}\, H \,{\cdot}\, P_k^ - \end{gather}

where the superscript $^{\prime} - ^{\prime}$ means predicted value, ${K_k}$ is the Kalman gain, ${b_k}$ is the vector of observed minus computed and ${\hat{x}_k}$ is the estimated parameters using the float solution. After the float solution is obtained, float DD ambiguities are fixed to the integers using the LAMBDA (Least-squares AMBiguity Decorrelation Adjustment) method (Teunissen, 1995; Li et al., 2018) as:

(22)\begin{equation}\check{a} = \textrm{min}{(a - \hat{a})^T} \,{\cdot}\, Q_{\hat{a}}^{ - 1} \,{\cdot}\, ({a - \hat{a}} )\end{equation}

where a is the vector of unknown integer DD ambiguities, $\hat{a}$ is the vector integer least squares estimates of the ambiguities and ${Q_{\hat{a}}}$ is the variance–covariance matrix of $\hat{a}$, and $\textrm{min}$ means the integer least square (ILS) estimates solution of the minimization problem (Chang et al., 2005). After the ambiguity validation, using something such as a fixed failure ratio test (Hou et al., 2016), the float parameters estimates and the corresponding variance–covariance matrix can be updated to ambiguity constrained solutions by solving the following equations (Odolinski et al., 2014; Li et al., 2018):

(23)\begin{gather}\check{b} = \hat{b} - {Q_{\hat{b}\hat{a}}} \,{\cdot}\, Q_{\hat{a}\hat{a}}^{ - 1} \,{\cdot}\, \left( {\hat{a} - \breve{a} } \right)\end{gather}

(24)\begin{gather}Q_{\check{b}\check{b}} = {Q_{\hat{b}\hat{a}}} - ({Q_{\hat{b}\hat{a}}} \,{\cdot}\, Q_{\hat{a}\hat{a}}^{ - 1} \,{\cdot}\, {Q_{a\hat{b}}})\end{gather}

where $\check{b}$ is the fixed solution of the estimated parameters and ${Q_{\check{b} \check{b}}}$ is the variance–covariance matrix of $\check{b}$. The vector of unknown parameters from ambiguity constrained solutions based on the uncombined DD observations can be expressed as:

(25)\begin{equation}\check{b} = [{{X_r},{Y_r},{Z_r},{\gamma^2}\nabla \Delta I_{rb,f}^{ij}, \nabla \Delta ZWD} ]\end{equation}

When the integer ambiguities of all satellites cannot be fixed using the LAMBDA method, partial ambiguity resolution (PAR) can be considered (Li et al., 2018). In this study, AR for BDS-3 geostationary orbit (GEO) satellites was not conducted due to the significantly poor orbit accuracy compared to GPS and Galileo.

In this study, DD wide-lane AR with the geometry-based model (Gao et al., 2015) was used to assist the medium- to long-range RTK. In this approach, the baseline parameters and DD wide-lane ambiguities are estimated together.

3. Experimental setup

RTK processes were conducted for 10 days (DOY: 060-069) in 2022 for approximately 50, 150, 250, 350, 450 and 550 km baseline lengths. Three different baselines were chosen for each length of baselines. A total of 18 baselines were formed by the real GNSS data from the European Reference Frame (EUREF) permanent GNSS network. The length and elevation differences of all the baselines are listed in Table 1. The average results obtained from three baselines for each baseline length were given (46-57-59 km for 50 km, 149-150-153 for 150 km, and so on). The distributions of the GNSS stations are given in Figure 1.

Table 1. Interstation distances and heights for the baselines

Figure 1. Distribution of the EUREF GNSS stations used in this study

In the selection of GNSS stations, it was confirmed that all chosen GNSS stations can track all operational GPS, Galileo and BDS-3 satellites. This confirmation is particularly important for BDS-3, because a few types of GNSS receivers can track all available BDS-3 satellites (Cao et al., 2021; Steigenberger et al., 2022). Since data loss in the RINEX files can cause the convergence state reinitialization, data availability in terms of epoch and corresponding phase/code frequencies was checked using in-house software. The RINEX data below 100% epoch and 98% corresponding phase/code frequency availability, respectively, were removed from the processing to avoid the impact of spurious reconvergence due to the data loss.

Since the IGS SINEX file does not contain each chosen station reference coordinates, precise point positioning with ambiguity resolution (PPP-AR) using the undifferenced uncombined observations was employed to obtain the reference coordinates of the stations. CSRS-PPP online positioning service was used to conduct PPP-AR (Banville et al., 2021). GPS + GLONASS combination (with integer ambiguity and real-value ambiguity for GPS and GLONASS) using static adjustment was employed for each station's daily data for 10 days. The average of 10 days’ results was taken as a reference coordinate for each station. The PPP-AR fix rate was higher than 98% for each station and each day, which indicates the reliability of the solutions.

For each baseline, GPS-only (G), Galileo-only (E), BDS-3-only (C3), GPS + Galileo (GE) and GPS + Galileo + BDS-3 (GEC3) combinations were used under 10°, 20° and 30° cutoff angles. Since most of the GLONASS satellites are broadcasting signals in frequency division multiple access (FDMA) structure, AR for GLONASS satellites is challenging. Interfrequency bias (IFB) parameter can be added for each frequency in the Kalman filter, but this significantly reduces the redundancy. If the same receiver is used for the base and rover, estimation of the IFB parameter is not necessary (Brack et al., 2021). But this condition significantly restricts the flexibility of the users. Therefore, GLONASS satellites were not used in this study. RTK processing strategy was summarized in Table 2.

Table 2. RTK processing parameters

4. Results

The following results give redundancy, AR fixing rate, positioning accuracy and convergence time analysis.

4.1 Redundancy analysis

All epoch-wise solutions from G, E, C3, GE and GEC3 RTK processes were investigated in terms of DD satellite number and non-available epoch-wise solutions under 10°, 20° and 30° cutoff angles. For each baseline length, Figure 2 shows the mean DD satellite number of each epoch within 24|h for 10 days.

Figure 2. Mean DD satellite number for each baseline length and cutoff angle

For solving the unknown parameters from the uncombined DD model, the minimum number of DD satellite numbers between the base and rover should be equal to three for the corresponding epoch. As seen in Figure 2, the redundancy significantly decreases under high-degree cutoff angles. Particularly for the single-system RTK, restricted satellite visibility can cause non-available solutions for the corresponding epoch. Figure 2 also shows that Galileo satellite visibility is low compared with GPS and Galileo for most baseline lengths, especially for the 10° and 20° cutoff angles. All non-available epoch-wise solutions were computed with respect to each baseline and cutoff angle for G, E, C3, GE and GEC3 RTK processes. The results show that non-available epoch-wise solutions did not exist for multi-GNSS RTK modes (GE and GEC3) for each cutoff angle. Multi-GNSS RTK modes provided 100% solutions for each cutoff angle thanks to the increased redundancy versus the single-system RTK modes. Non-available epoch-wise solutions were generally accumulated under 20° and 30° cutoff angles for single-system RTK solutions. Table 3 shows the ratio of the non-available epoch-wise solutions for single-system RTK solutions for the 20° and 30° cutoff angles.

Table 3. The ratio of the non-available epoch-wise solutions (unit: %)

The results of non-available epoch-wise solutions show that the satellite visibility of Galileo is significantly low compared with GPS and BDS-3 for the 20° and 30° cutoff angles. Two representative stations (HERS and SOD3) were taken to investigate the satellite visibility of Galileo and BDS-3 under a 30° cutoff angle for DOY 060 (Figure 3). TRIMBLE GNSS mission planning software was used to produce the satellite visibility time series.

Figure 3. Mean DD satellite numbers

As seen in Figure 3, Galileo satellite visibility is significantly low compared with BDS-3 for HERS and SOD3 stations under a 30° cutoff angle. When common satellite number between the base and rover drops to three for the corresponding epoch, no RTK solution can be conducted. When the satellite visibility time series is investigated, it can be seen that the Galileo time series contain more sessions with fewer than four satellites compared with BDS-3 for the representative HERS and SOD3 stations. The other stations have similar satellite visibility patterns. Therefore, the non-available epoch-wise solutions of Galileo-only RTK are much higher compared with BDS-3-only RTK under the high-degree cutoff angles.

4.2 Fixing rate analysis

The fixing rate of full AR and partial AR and also the rate of the float solutions are given in Tables 4–6 for each baseline length, RTK mode and cutoff angle. Figure 4 also shows the all fixing rates (full + partial) for each RTK mode.

Table 4. Ambiguity fixing rate for 10° cutoff angle (unit: %, full/partial/float)

Table 5. Ambiguity fixing rate for 20° cutoff angle (unit: %, full/partial/float)

Table 6. Ambiguity fixing rate for 30° cutoff angle (unit: %, full/partial/float)

Figure 4. Total fixing rates for 10° cutoff angle (full + partial)

The results of the fixing rate show that the fixing rate of GPS-only RTK is best among the Galileo-only and BDS-3-only RTK for each baseline and cutoff angle. This can be attributed to the better GPS satellite visibility compared to Galileo and BDS-3 (Figure 2). As expected when the cutoff angle increases, the fixed rates decrease for the single-system RTK solutions due to poor satellite visibility. The results also indicate that partial AR plays an important role, especially for long baselines, high-degree cutoff angles and multi-GNSS RTK. In multi-GNSS RTK, the ambiguity dimensions increase significantly, resulting in a more difficult estimation of all ambiguities compared to the single-system RTK; therefore, the fixing rate of AR is reduced especially for long baselines (Wang, 2012). Instead of forcing to solve all ambiguities simultaneously, resolving a subset of ambiguities can be more effective, especially for long baselines and multi-GNSS RTK. For example, 2% full set and 74% partial ambiguities were resolved for 550 km using GEC3 RTK under a 20° cutoff angle. Ambiguity validation can be performed using the model-driven approach or data-driven approach. The success rate and fixed failure rate ratio test (FFRT) ambiguity validations should be implemented for model-driven and data-driven approaches, respectively (Hou and Verhagen, 2014). AR success rate only effective when the partial AR mode contains a model-driven approach because ambiguity subsets are determined based on the satellite geometry instead of measurements in the model-driven approach (Hou et al., 2022). In this study, the data-driven approach (ambiguity subsets are determined based on the measurements) was used with the FFRT while taking the threshold as 2.0. Therefore, AR success was not considered and set to 0 (Table 2).

4.3 Analysis of positioning accuracy

For analysing the positioning accuracy, coordinates derived from the fixed solutions were considered. The number of the fixed solutions was not the same among the RTK modes due to the different numbers of fixed/float and non-available solutions. For the statistically robust comparison, the fixed solutions that share the same epochs within each day among the RTK modes were taken for computing the positioning accuracy. In this way, all fixed coordinates from each RTK mode were derived from the same conditions, such as satellite geometry, multipath and atmospheric effects. Figures 5–7 show the RMSEs of the north, east and up components for each RTK mode and cutoff angle. Since the number of the satellite was significantly low under the 30° cutoff angle, especially for Galileo and BDS-3, dm-level positioning error was generally obtained for a single-system RTK. Therefore, only GE and GEC3 combinations were considered for the 30° cutoff angle.

Figure 5. RMSEs for each baseline length under 10° cutoff angle

Figure 6. RMSEs for each baseline length under 20° cutoff angle

Figure 7. RMSEs for each baseline length under 30° cutoff angle

To increase the visibility of the results, the horizontal (n/e) and vertical axes of the plots were scaled differently. In terms of single-system RTK solution accuracy, the results show that the lowest positioning accuracy was obtained from Galileo-only RTK among the GPS and BeiDou-3 for each baseline length and cutoff angle. Poor satellite visibility (Figures 2 and 3) was mainly responsible for this phenomenon. For a 20° degree cutoff, GPS produced the best positioning accuracy compared with BDS-3 for each baseline length. For a 10° degree cutoff angle, the positioning accuracy of BDS-3 is comparable to GPS for 50 km, 150 km, 250 km and 350 km baseline lengths, but for 450 km and 550 km baseline lengths, the positioning accuracy of BDS-3 is much better than GPS. While satellite clock errors are eliminated from DD observations, satellite orbit errors cannot be fully eliminated from DD observations. The rule of thumb for the satellite error impact on the estimated receiver position in a single-difference observation as a baseline length function can be described as (Teunissen and Montenbruck, 2017):

(26)\begin{equation}{\Delta _{\textrm{receiver}}} \le \frac{L}{\rho } \,{\cdot}\, {\Delta _{\textrm{satellite}}} \end{equation}

where ${\Delta _{\textrm{receiver}}}$ is the satellite error impact on the estimated receiver position, L is the baseline length, $\rho$ is the distance between the satellite and receiver, and ${\Delta _{\textrm{satellite}}}$ is the satellite orbit error. As seen from Equation (26), the satellite error impact on the estimated receiver position depends on the baseline length. In according with the GNSS orbit signal-in-space range error (SISRE) (Liu et al., 2023), orbit SISRE of GPS, Galileo and BDS-3 on a baseline of 550 km project maximum 1 ⋅ 6 cm, 0 ⋅ 7 cm, and 0 ⋅ 3 cm error, respectively, on the estimated receiver position, using the nominal altitude of the GNSS satellites. The low positioning accuracy of GPS compared with BDS-3 for long baselines can be attributed to the minimum SISRE orbit error of BDS-3. It can be seen that the projected orbit SISRE errors on the estimated receiver position for 550 km baseline length are insignificant, considering the overall accuracy obtained from multi-GNSS RTK for 550 km. For single- and double-difference observations, the influence of the along-track and cross-track orbital errors is stronger than the radial error. For double-difference observation, the magnitude of these orbital errors varies, depending on the mutual position of the satellites.

When multi-GNSS RTK solutions are considered, the results show that the GEC3 combination has the minimum horizontal and vertical RMSE compared with the GE combination for each baseline and cutoff angle. RMSEs of horizontal–vertical components of the GEC3 RTK solutions under 10° degree cutoff angle were computed as 0 ⋅ 7–1 ⋅ 5, 0 ⋅ 9–2 ⋅ 4, 1–3 ⋅ 2, 1 ⋅ 4–3 ⋅ 5, 1 ⋅ 3–3 ⋅ 5 and 1 ⋅ 6–4 ⋅ 1 cm, for 50, 150, 250, 350, 450 and 550 km length of baselines, respectively. The results also indicate that as baseline length increases, the accuracy of the single-system and multi-GNSS RTK solutions degrades to a certain extent, but the accuracy degradation of multi-GNSS RTK solutions is generally less than the single-system RTK solutions under a 10° degree cutoff angle. The results of RMSEs obtained from the GEC3 RTK solutions under a 10° degree cutoff angle are comparable to the overall accuracy of network-based RTK solutions (Edwards et al., 2010; Cina et al., 2015; Laoniphon et al., 2021; Yeh et al., 2022).

4 ⋅ 4 Analysis of the convergence time

The convergence time is an important factor for medium- and long-baseline RTK applications. Despite the NRTK applications, instantaneous single epoch AR is generally more difficult for medium- and long-baseline RTK, and more than single epoch is generally required for ambiguity dilution of precision (ADOP) to convergence (Odolinski et al., 2015). The convergence time was estimated every 3 h within one day by re-initializing the Kalman filter. The convergence time was defined as when the positioning error is less than 10 cm at the current epoch and the following 20 epochs. The convergence time was computed separately for the north, east and up components. The convergence time was computed only for a 10° degree cutoff angle. Figure 8 shows the mean convergence time for each RTK mode and baseline length. Similarly to the positioning accuracy analysis, the horizontal (n/e) and vertical axes of the plots were scaled differently to increase visibility.

Figure 8. Mean convergence time for each baseline length under 10° degree cutoff angle

When the single-system RTK solutions are considered, the convergence time results show that GPS has the fastest convergence time compared with Galileo and BDS-3. No distinct pattern was found for Galileo and BDS-3 in terms of the convergence time performance. Similar to the positioning accuracy results, the GEC3 RTK combination produced the best convergence time results among all combinations. The mean convergence time of the horizontal–vertical components of the GEC3 RTK solutions under a 10° degree cutoff angle was computed as 0 ⋅ 9–1 ⋅ 9, 3 ⋅ 4–3 ⋅ 4, 5 ⋅ 8–4, 9 ⋅ 7–10 ⋅ 3, 12 ⋅ 1–10 ⋅ 7 and 13–14 ⋅ 1 min for 50, 150, 250, 350, 450 and 550 km length of baselines, respectively.

Instantaneous AR was also investigated in this study. Instantaneous AR was defined as when the convergence time for each component and TTFF were zero. In this situation, AR is derived from one single epoch. Instantaneous AR is seen for GE and GEC3 RTK solutions for 50, 150, 250 and 350 km baseline lengths. The ratio of instantaneous AR for all data sets of the GE and GEC3 RTK solutions was computed for the corresponding baseline lengths and is given in Table 7.

Table 7. The ratio of instantaneous AR for 10° cutoff angle (unit: %)

As seen in Table 7, instantaneous AR is more likely to occur for short baselines. The ratio of the instantaneous AR from GEC3 RTK solutions is significantly higher than the GE RTK solutions. It proves that BDS-3 significantly contributed to GPS + Galileo RTK in terms of instantaneous AR.

5. Conclusions

Since BDS-3 has reached the FOC, RTK positioning with AR using GPS, Galileo and BDS-3 has become possible. In this study, single-system GPS (G), Galileo (E), BeiDou-3 (C3) and multi-GNSS (GE and GEC3) single-baseline RTK solutions under 10°, 20° and 30° cutoff angles were conducted using six different lengths of baselines (~50, ~150, ~250, ~350, ~450 and ~550 km) within the EUREF permanent GNSS network for 10 days of 24 h data in 2022. The troposphere and ionosphere weighted model with wide-lane partial AR geometry-based model was used for each RTK processing. The main conclusions can be summarized as follows:

1. (1) For the corresponding baselines, Galileo satellite visibility was significantly low compared with GPS and BDS-3, especially for 10° and 20° cutoff angles, and it caused the poor positioning performance of Galileo among the single-system RTK solutions. For a 10° cutoff angle, GPS and BDS-3 horizontal positioning accuracies are comparable up toa 350 km baseline length. For longer baselines, the accuracy of BDS-3 horizontal positioning is much better than GPS. For the vertical component, the accuracy of BDS-3 is much better than GPS for baselines longer than 150 km. This phenomenon can be explained by the lowest orbit SISRE error of BDS-3 among the other GNSS constellations. For a 20° cutoff angle, GPS horizontal and vertical positioning accuracy are much better than BDS-3 for each baseline length due to the higher satellite visibility of GPS (Figure 2).

2. (2) In terms of multi-GNSS RTK positioning accuracy (GE and GEC3), the GEC3 combination produced the best positioning accuracy for each baseline length and cutoff angle. Even under a 30° cutoff angle, the GEC3 RTK solutions produced cm-level positioning accuracy with 100% epoch availability solutions, thanks to the high redundancy.

3. (3) The convergence time analysis showed that instantaneous AR was feasible for multi-GNSS RTK solutions for baselines up to 350 km under a 10° degree cutoff angle. The results indicated that the ratio of instantaneous AR from the GEC3 RTK solutions is significantly higher than the GE RTK solutions for baselines up to 350 km. No instantaneous AR was observed for the baseline lengths longer than 350 km. Similarly to the positioning accuracy results, the GEC3 RTK combination has the fastest convergence time compared with the other RTK combinations. As expected, when the baseline length increases, the convergence time increases for each RTK mode. The minimum and maximum mean horizontal–vertical convergence time of the GEC3 RTK solutions were computed as 0 ⋅ 9–1 ⋅ 9 and 13–14 ⋅ 1 min for 50 and 550 km baseline lengths, respectively.

4. (4) Orbit SISRE errors of GNSS constellations proved that the influence of the broadcast orbit error can be safely ignored for Galileo and BDS-3 for baseline lengths up to 550 km. For GPS, a maximum 1 ⋅ 6 cm error on the estimated receiver position was computed for a 550 km baseline length. When considering the multi-GNSS RTK, orbit SISRE of GPS can be largely compensated by using Galileo and BDS-3 for baseline length within 550 km.

The results of this study showed that the GEC3 RTK solutions for medium and long baselines are promising, thanks to the completed FOC of Galileo and BDS-3. The accuracy of the converged GEC3 RTK solutions (a few cm in horizontal and cm-level in vertical) is comparable to NRTK solutions even for a 550 km baseline length. However, a significantly longer convergence time is required by single-baseline multi-GNSS RTK solutions for baselines longer than 250 km. This is the major disadvantage of the single-baseline multi-GNSS RTK for long baselines compared with the NRTK. But when considering the costs of the NRTK for large countries or large offshore areas, single-baseline multi-GNSS RTK can dramatically reduce the costs of NRTK infrastructure maintenance.