Sub-meter-level navigation with an enhanced multi-GNSS single-point positioning algorithm using iGMAS ultra-rapid products

Abstract

Ultra-rapid products have the advantage of being used in real-time positioning with no external connections. In this study, these products provided by the international GNSS Monitoring and Assessment System (iGMAS) for four global constellations (GPS, GLONASS, Galileo and BDS-3) were assessed in terms of service rate and accuracy in navigation. In this regard, a MATLAB-based in-house code solving the problem was developed for all possible combinations of the constellations. To explore the effectiveness of the iGMAS products, the same dataset has been also processed using GFZ rapid products. The results demonstrate that the GPS and Galileo solutions were substantially comparable to the rapid products concerning service rate and accuracy, but that the GLONASS and BDS-3 iGMAS products require some enhancements. In addition, the Galileo solution produced remarkably good results both individually and in combination. The GPS/GLONASS/Galileo/BDS-3 SPP solution generated a mean root mean square (RMS) error of 0 ⋅ 54 m horizontally and 0 ⋅ 89 m vertically. Thus, GPS-only, GLONASS-only, Galileo-only and BDS-3-only solutions were improved by 42%, 79%, 28% and 74% in 3D mean RMS error with the quad system solutions, respectively.

1. Introduction

In the wake of GPS, operational usage of Global Navigation Satellite Systems (GNSS) technology, together with recently developed projects, such as GLONASS, Galileo and the BeiDou Satellite System (BDS), has become more powerful and reliable in satellite-based positioning (Li et al., 2015). Especially after 2010, modernisations and innovations in these systems rendered each of them operative in worldwide use. In the renovation stage, the revisions made in the Russian satellites system GLONASS, nominal 24 Medium Earth Orbit (MEO) operational satellites have been in service since 2012 (Revnivykh et al., 2017). Today, the constellation consists of 18 GLONASS-M satellites, two GLONASS-M+ satellites and one GLONASS-K satellite (Bury et al., 2022). The Chinese BDS was initially realised as a demonstration system BDS-1, followed by the regional BDS-2, which specifically serves the Asia-Pacific region, was completed at the end of 2012. The third-generation navigation system BDS-3, which aims to serve the whole world, was announced on 27 December 2018. Currently, this system includes 24 MEO, three Geostationary Earth Orbit (GEO), and three Inclined Geosynchronous Orbit (IGSO) satellites (Yang et al., 2019; CSNO-TARC, 2022). The Galileo constellation project was started by the European Union towards the end of 2005. The success of the only-Galileo positioning solution was tested by launching two pairs of In-Orbit Validation (IOV) satellites in 2011 and 2012. In 2014, the constellation was developed with Full Operational Capability (FOC) satellites. As of February 2022, Galileo contains usable 22 MEO satellites (EUSPA, 2022). GPS, GLONASS, Galileo and BDS-3 are the basic components of Positioning Navigation and Timing (PNT) services capable to provide information to users independently.

Each constellation has a comparable structure, mathematical modelling and observation type. Therefore, an integrated positioning has become possible, which is named multi-GNSS, together with the new generation receivers (Montenbruck et al., 2014). Obtaining precise orbit and clock products of these systems has become a requirement for multi-GNSS precise point positioning (PPP) applications. International GNSS Service (IGS) began providing these products for GPS satellites in 1994 (Beutler et al., 1999). The orbits of all operational GLONASS satellites were included in the IGS precise products in 2004 (Beutler et al., 2009). With the advent of Galileo, BDS and the other regional navigation systems within the framework of GNSS, the IGS Multi-GNSS Experiment (MGEX) project was put on the agenda with the modern receivers deployed in 2013. Thus, different analysis centres (ACs) have started to calculate precise products utilising MGEX network data for GPS, GLONASS, Galileo and BDS (Montenbruck et al., 2017). According to the publishing timeliness, IGS precise products are classified as final, rapid and ultra-rapid, with the final product having a delay of 12–18 days, the rapid product having a delay of 17 h and the ultra-rapid product (observed half) having a delay of 3–9 h, although another type of ultra-rapid product (predicted half) can be acquired in real-time (Dow et al., 2009; Shi et al., 2015).

Thanks to the MGEX products, especially multi-GNSS, PPP has become applicable. Therefore, many PPP studies in the literature started using GPS-GLONASS and continued with the integration of the Galileo and BDS systems (Cai and Gao, 2013; Jiao et al., 2019; Kiliszek and Kroszczyński, 2020). Parallel to the post-processing applications, IGS Real-Time Service (RTS) was officially initiated in 2013 which provides orbit and clock correction for the broadcast ephemeris to the real-time PPP users (Hadaś and Bosy, 2015; Wang et al., 2018; Ge et al., 2021). A strong external connection is required to acquire RTS orbit and clock corrections, which are broadcasted via an internet stream. During real-time positioning, an unexpected communication failure or network connection latency will cause a gap of several seconds to hours in receiving orbital and clock data. Hence, an alternate approach to the RTS might be considered by utilising ultra-rapid products to ensure positioning uninterruptedly (El-Mowafy et al., 2017). The ultra-rapid products are used in a variety of studies, including near-time troposphere estimation (Satirapod et al., 2011; Hadaś et al., 2013), precise timing (Cerretto et al., 2012), clock estimation (Cao et al., 2022), examination of dynamic displacements (Yigit et al., 2020), earthquake analysis (Shu et al., 2020), time transfer (Chen et al., 2022), real-time PPP and SPP applications (Elsobeiey and Al-Harbi, 2016; Bahadur and Nohutcu, 2021; Bahadur, 2022; Geng et al., 2022; Jiang et al., 2022; Jiao and Song, 2022; Ogutcu and Farhan, 2022).

Nowadays, ultra-rapid products are estimated by both IGS and International GNSS Monitoring and Assessment System (iGMAS) independently. IGS uses eight ACs: Centre for Orbit Determination in Europe (CODE), Center National d'Etudes Spatiales and Collecte Localisation Satellites (CNES/ CLS), Wuhan University (WHU), Deutsches GeoForschungs Zentrum (GFZ), Japan Aerospace Exploration Agency (JAXA), Technische Universität München (TUM), Shanghai Astronomical Observatory (SHAO), and Information and Analysis Centre (IAC) (Shen et al., 2021). In addition to the IGS products, the Chinese project called iGMAS provides products obtained from combinations of its own ACs. Currently, 12 iGMAS ACs process the multi-GNSS data collected from about 30 iGMAS tracking stations and generate precise GNSS products. Each ACs has the same mission and operates independently (iGMAS, 2022; Zhou et al., 2022). When comparing IGS and iGMAS, the main difference is that they are based on two different networks. In addition, iGMAS offers users an important opportunity to provide BDS ultra-rapid products with the contribution of 12 ACs, whereas the IGS offers this only by WHU and CNES. Therefore, the iGMAS ultra-rapid products availability for BDS looks stronger than for IGS. iGMAS also provides ultra-rapid clock products for four GNSS systems. Consequently, iGMAS ultra-rapid products make positioning possible individually or in different combinations of GNSS systems in real-time.

The current status of iGMAS products, particularly in the scope of BDS-3, directed us to use them in navigation problems. Hence, the single-point positioning (SPP) technique was applied to investigate the contribution of the real-time products to the navigation solution. For the processing of data, a comprehensive MATLAB code was developed in this study. The dataset used was taken from 16 IGS MGEX stations. A detailed analysis was performed in terms of positioning accuracy with single, dual, triple and quad systems. The service rates provided by iGMAS products for GPS, GLONASS, Galileo and BDS-3 were investigated. Furthermore, for evaluating the quality of iGMAS ultra-rapid products, GFZ rapid products were also employed, and their performance was compared with those obtained from the iGMAS products.

This study is divided into five parts. In Section 2, the multi-GNSS SPP functional and stochastic model as well as the outlier detection technique are explained in detail. In Section 3, the introduction of the dataset and the evaluation strategy of the developed code is given. In Section 4, an extensive accuracy analysis has been performed for all combinations with iGMAS products. Finally, Section 5 summarises the conclusions obtained from the study.

2. Methodology

2.1 Multi-GNSS SPP model

The multi-GNSS SPP technique provides users with time and position information based on code observations from multiple satellite systems. The equation of code ($P_{i,r}^{s,j}$) measurements on the $i\textrm{th}$ ($i = \,1,\,2$) frequency can be written as:

(1)\begin{equation}P_{i,r}^{s,j} = \rho _r^{s,j} + cdt_r^s - cd{T^{s,j}} + T{r^{s,j}} + {I_i}^{s,j} + b_{i,r}^s - b_i^{s,j} + {\varepsilon _{P_{i,r}^{s,j}}}\end{equation}

where superscript s represents the GNSS constellations ($G$: GPS, $R$: GLONASS, $E$: Galileo, $C$: BDS-3) and $j$ denotes the evaluated satellite of each constellation; subscript r and variable i refer to GNSS receiver and frequency of GNSS signal, respectively; $\rho _r^{s,j}$ is the geometric range in meters between the satellite and the receiver; c is the speed of light in a vacuum in meters per second; $cd{T^{s,j}}$ and $cdt_r^s$ are satellite and receiver clock offsets in seconds, respectively; $T{r^{s,j}}$ is the tropospheric delay in meters; ${I_i}^{s,j}$ is the first-order ionospheric delay depending on $i\textrm{th}$ frequency in meters; $b_{i,r}^s$ and $b_i^{s,j}$ are the receiver and satellite hardware code delay in code observations at $i\textrm{th}$ frequency in meters, respectively; ${\varepsilon _{P_{i,r}^{s,j}}}$ consists of the measurement noise, multipath, orbital error, etc., in meters. It should be outlined that some systematic errors, such as satellite phase offset, Earth rotation effects, relativistic effects and solid tide effect, were removed by using common models (Kouba and Héroux, 2001). The first-order ionospheric delay was eliminated by using ionosphere-free ($IF$) combination of the dual-frequency data.

Precise products are provided by many ACs to remove satellite orbit and clock offset from the equation. The satellite clock offset products, which have a high correlation with the satellite hardware code delays, are generated by ionosphere-free code observations. Thereby, these estimated satellite clock offset products include the $IF$ combination of satellite hardware code delays. Furthermore, the receiver hardware code bias which is strongly correlated with the receiver clock offset is absorbed by the receiver clock offset parameter in the estimation (Kouba and Héroux, 2001; Zhou et al., 2019). Thus, these two parameters were completely disregarded in the model during the evaluation process. This can be expressed with the following equations:

(2)\begin{gather}\tilde{P}_{IF,r}^{s,j} = P_{IF,r}^{s,j} + cd{\tilde{T}^{s,j}}\end{gather}

(3)\begin{gather}cd{\tilde{T}^{s,j}} = cd{T^{s,j}} + b_{IF}^{s,j}\end{gather}

(4)\begin{gather}cd\tilde{t}_r^s = cdt_r^s + b_{IF,r}^s\end{gather}

where $\tilde{P}_{IF,r}^{s,j}$ is the pseudorange observation corrected by the satellite clock offset; and $cd{\tilde{T}^{s,j}}$ and $cd\tilde{t}_r^s$ are redefined satellite and receiver clock offset, respectively. When using multi-GNSS, each constellation has different hardware delays and time scales. For this reason, two techniques are available to evaluate receiver clock offset in multi-GNSS. The first one is to independently obtain the receiver clock offset for each constellation, which is preferred in this study. The second method is to estimate the receiver clock offset for a selected reference system and the inter-system bias (ISB) parameters between the reference system and the others (Lou et al., 2016; Liu et al., 2019). For tropospheric delay, the UNB3m model, developed by the University of New Brunswick, using the Niell mapping function and known as a hybrid model, was used in this study (Niell, 1996; Leandro et al., 2008). Using dual-frequency GNSS data, the $IF$ pseudo-range observations of the multi-constellation integration may be stated as follows:

(5)\begin{equation}\left\{\begin{array}{@{}l} \tilde{P}_{IF,r}^{G,j} = \rho_r^{G,j} + cd\tilde{t}_r^G + {\varepsilon_{P_{IF,r}^{G,j}}}\\ \tilde{P}_{IF,r}^{R,j} = \rho_r^{R,j} + cd\tilde{t}_r^R + {\varepsilon_{P_{IF,r}^{R,j}}}\\ \tilde{P}_{IF,r}^{E,j} = \rho_r^{E,j} + cd\tilde{t}_r^E + {\varepsilon_{P_{IF,r}^{E,j}}}\\ \tilde{P}_{IF,r}^{C,j} = \rho_r^{C,j} + cd\tilde{t}_r^C + {\varepsilon_{P_{IF,r}^{C,j}}} \end{array} \right.\end{equation}

The observation model [Equation (5)] should be linearised around the approximate receiver position. For the estimation of receiver coordinates and clock parameters, the common Gauss–Markov model together with the weighted least squares method is used. The equations are given by (Koch, 1999)

(6)\begin{gather} E(\boldsymbol{l}) = \boldsymbol{l} + \boldsymbol{v} = \boldsymbol{Ax};\;E(\boldsymbol{v}) = 0\;\textrm{and}\;D(\boldsymbol{l}) = \sigma _0^2{\boldsymbol{Q}_{\boldsymbol{l}}} = \sigma _0^2{\boldsymbol{P}^{-1}}\end{gather}

(7)\begin{gather}\boldsymbol{x} = {({\boldsymbol{A}^T}\boldsymbol{PA})^{ - 1}}{\boldsymbol{A}^T}\boldsymbol{Pl}\end{gather}

where $E(\boldsymbol{l})$ is the expected value of observations; $D(\boldsymbol{l})$ is the variance–covariance matrix of the observations; $\boldsymbol{v}$ is the posterior residuals of measurement vector; $\sigma _0^2$ is the a priori variance factor; $\boldsymbol{x}$ is the estimated unknown parameter vector; $\boldsymbol{A}$ is the coefficients matrix of the linearised equation; $\boldsymbol{l}$ is the measurement vector, which is the difference between the corrected pseudo-range measurement and the distance calculated using the satellite coordinates and the receiver approximate coordinates; $\boldsymbol{P}$ is the weight matrix of the observations; and ${\boldsymbol{Q}_{\boldsymbol{l}}}$ is the cofactor matrix of observations. The multi-GNSS SPP observation model can be explicitly expressed as follows:

(8)\begin{equation}\boldsymbol{v} = \boldsymbol{Ax} - \boldsymbol{l} = \left[ \begin{array}{@{}l@{}} {\boldsymbol{v}^G}\\ {\boldsymbol{v}^R}\\ {\boldsymbol{v}^E}\\ {\boldsymbol{v}^C} \end{array} \right] = \left[ {\begin{array}{@{}lllll@{}} {{\boldsymbol{A}^G}}& {{\boldsymbol{B}^G}}& 0& 0& 0\\ {{\boldsymbol{A}^R}}& 0& {{\boldsymbol{B}^R}}& 0& 0\\ {{\boldsymbol{A}^E}}& 0& 0& {{\boldsymbol{B}^E}}& 0\\ {{\boldsymbol{A}^C}}& 0& 0& 0& {{\boldsymbol{B}^C}} \end{array}} \right]\,\left[ \begin{array}{@{}l@{}} \,\,\boldsymbol{\varDelta r}\\ cd{t^G}\\ cd{t^R}\\ cd{t^E}\\ cd{t^C} \end{array} \right] - \left[ \begin{array}{@{}l@{}} {\boldsymbol{l}^G}\\ {\boldsymbol{l}^R}\\ {\boldsymbol{l}^E}\\ {\boldsymbol{l}^C} \end{array} \right]\end{equation}

where $\boldsymbol{\varDelta r} = {\left[ {\begin{array}{*{20}{l}} {\varDelta x}& {\varDelta y}& {\varDelta z} \end{array}} \right]^T}$ parameters vector for receiver coordinates; $\boldsymbol{B}$ is the clock offsets coefficients vector in that each element equals one; and ${\boldsymbol{A}^G},{\boldsymbol{A}^R},{\boldsymbol{A}^E},{\boldsymbol{A}^C},$ in the model are sub-coefficients matrices of coordinate unknowns for different GNSS. These sub-matrices can be expressed for GPS observations as follows:

(9)\begin{gather} {[}{\boldsymbol{v}^G}] = [\begin{array}{@{}ll@{}} {{\boldsymbol{A}^G}}& {{\boldsymbol{B}^G}} \end{array}]\,\,\left[\begin{array}{@{}l@{}} \,\,\boldsymbol{\varDelta r}\\ cd{t^G} \end{array} \right] - [{{\boldsymbol{l}^G}} ]\end{gather}

(10)\begin{gather} {\boldsymbol{A}^G} = \left[ \begin{array}{@{}l@{}} {\boldsymbol{a}_1}\\ {\boldsymbol{a}_{2\,}}\\ \,\, \vdots \\ {\boldsymbol{a}_{n\,}} \end{array} \right]\quad ;\quad {\boldsymbol{a}_i} = \left[ {\frac{{ - ({X^i} - {X_r})}}{{||{{\boldsymbol{r}^s} - {\boldsymbol{r}_r}} ||}}\quad \frac{{ - ({Y^i} - {Y_r})}}{{||{{\boldsymbol{r}^s} - {\boldsymbol{r}_r}} ||}}\quad \frac{{ - ({Z^i} - {Z_r})}}{{||{{\boldsymbol{r}^s} - {\boldsymbol{r}_r}} ||}}} \right]\end{gather}

where ${\boldsymbol{r}^s} = ({X^i},\,{Y^i},{Z^i})$is the ith satellite position vector; ${\boldsymbol{r}_r} = ({X_r},\,{Y_r},\,{Z_r})$is the initial position vector of the receiver; and n is the observed number of satellites.

In the stage of parameter estimation, the stochastic model should be formed in the most accurate way to get optimum results. This issue requires more attention, especially in multi-GNSS applications. The variance–covariance matrix should represent actual stochastic conditions of observations. In multi-GNSS applications, this situation can be outlined as the variance of the observations, depending on the GNSS system and the function of the satellite elevation angle, which can be written as (Pan et al., 2017)

(11)\begin{equation}{\sigma ^2} = \frac{{\sigma _0^2}}{{{{\sin }^2}(Ele)}}\end{equation}

where ${\sigma ^2}$ is the variance value of the code observation; $Ele$ is the satellite elevation angle; and $\sigma _0^2$ is the a priori variance of the related navigation system. Unlike GPS, Galileo and BDS, which use Code Division Multiple Access (CDMA) signals, GLONASS uses Frequency Division Multiple Access (FDMA) signals, except for the recently modernised ones. It should be noted that different hardware delays exist on GLONASS receiver channels as different frequencies are generated for each satellite. Therefore, the GLONASS code observations contain inter-frequency biases (Wanninger, 2012).

Since it is usually not desirable to determine these values for each satellite, this term is neglected in the stochastic model by assigning a higher variance value to the GLONASS code observations. As a result of the solution, the effect of this term is seen in the residuals of the GLONASS code observations (Cai and Gao, 2013). In addition, the BDS GEO satellites have lower performance than do the BDS MEO and BDS IGSO satellites. The process should be advanced by reducing the weight of these satellites in the stochastic model (Zhang and Pan, 2022). In the construction of the stochastic model, the chosen standard deviation (${\sigma _0}$) of observations for GPS, GLONASS, Galileo, BDS-3 (MEO/IGSO) and BDS-3 (GEO) were taken as 0 ⋅ 3, 0 ⋅ 6, 0 ⋅ 3, 0 ⋅ 6 and 1 ⋅ 2 m, respectively.

2.2 Outlier detection

The receiver coordinates and clock offset parameters were estimated epoch-by-epoch using the least squares method without any constraints between observation epochs. It is necessary to pay attention to anomalous measurements in order to obtain optimal and reliable results. These measurements (i.e. outliers) may be seen in the dataset due to both hardware and environmental or external reasons. Thus, outliers should be identified and their impact on the solution should be minimised. Otherwise, undesired observations can degrade the positioning performance (Angrisano et al., 2020). The median absolute deviation (MAD) estimator is adopted because it is highly resistant to the outliers, with a 50% breakdown point. The MAD estimator is formulated as (Rousseeuw and Leroy, 1987; Hekimoglu et al., 2014)

(12)\begin{gather}\textrm{MAD} = 1.4826 \times \textrm{median}\{{|{{l_1} - \textrm{med}} |,\,|{{l_2} - \textrm{med}} |,\, \cdots ,\,|{{l_n} - \textrm{med}} |} \}\end{gather}

(13)\begin{gather}\textrm{med} = \textrm{median}\{{{l_1},\,{l_2},\, \cdots ,{l_n}} \}\end{gather}

where ${l_i}$ is the value of the ith observation in the measurement vector ($\boldsymbol{l}$) and n is the length of the vector. The $k \times \textrm{MAD}$ threshold is compared with each $({|{l - \textrm{med}} |} )$value. If a value is greater than the threshold value, the observation is marked as an outlier and excluded from the evaluation. Herein, the k value was set to 3. $\textrm{MAD}$ represents the standard deviation for normally distributed data. With the $3\sigma$ edit rule, which is generally called the Hampel identifier, outliers can be detected (Pearson, 2001; Chiang et al., 2003). Note that the approach performs inaccurately in multi-GNSS processing since different systems have different receiver clock offsets. Therefore, this method can be used to detect gross errors in pre-analysing each system independently. As a result, it has functioned very efficiently in epochs where the degrees of freedom are low and outliers cannot be detected by common statistical methods.

Following the processing stage of gross errors detection, for SPP results, an appropriate weight model is a crucial requirement. A robust weighting method was implemented to mitigate the impact of outliers that cannot be detected with the MAD technique in this study. The method creates an equivalent weighting matrix according to certain conditions by using the information obtained from the least squares solution. The Institute of Geodesy and Geophysics (IGG)-III function was employed to calculate the matrix. The weighting of observations was adjusted with IGG-III as (Yang et al., 2002; Guo and Zhang, 2014)

(14)\begin{gather}{\gamma _i} = \left\{ {\begin{array}{@{}ll} 1& {|{{{\tilde{v}}_i}} |\le {k_0}}\\ {\dfrac{{{k_0}}}{{|{{{\tilde{v}}_i}} |}}{{\left( {\dfrac{{{k_1} - |{{{\tilde{v}}_i}} |}}{{{k_1} - {k_0}}}} \right)}^2}}& {{k_0} < |{{{\tilde{v}}_i}} |\le {k_1}}\\ 0& {|{{{\tilde{v}}_i}} |> {k_1}\,} \end{array}} \right.\end{gather}

(15)\begin{gather}{\bar{P}_i} = {P_i}\,{\gamma _{ii}}\quad ;\quad \overline {\boldsymbol{P} }\left[ { = \begin{array}{*{20}{c}} {\dfrac{{{\gamma_{11}}}}{{\sigma_1^2}}}& {}& {}& {}\\ {}& {\dfrac{{{\gamma_{22}}}}{{\sigma_2^2}}}& {}& {}\\ {}& {}& \ddots & {}\\ {}& {}& {}& {\dfrac{{{\gamma_{nn}}}}{{\sigma_{nn}^2}}} \end{array}} \right]\end{gather}

where $\overline {\boldsymbol{P} }$ is the equivalent weighting matrix; $\gamma$ is the inflation factor; subscript i represents each observation; ${k_0}$ and ${k_1}$ are two threshold values and usually set ${k_0}$ = 1 ⋅ 5–3, ${k_1}$ = 3–8; and ${\tilde{v}_i}$ is the standardised residual, expressed by

(16)\begin{equation}{\tilde{v}_i} = \frac{{{v_i}}}{{\sqrt {\tilde{\sigma }_0^2{Q_{{v_i}}}} }}\end{equation}

where $\tilde{\sigma }_0^2$ is the estimate of unit weighted variance and ${Q_{{v_i}}}$ is the cofactor matrix of the residual:

(17)\begin{gather}\tilde{\sigma }_0^2 = \frac{{{\boldsymbol{v}^T}\boldsymbol{Pv}}}{{n - u}}\end{gather}

(18)\begin{gather}{\boldsymbol{Q}_v} = {\boldsymbol{P}^{ - 1}} - \boldsymbol{A}{({{\boldsymbol{A}^T}\boldsymbol{PA}} )^{ - 1}}{\boldsymbol{A}^T}\end{gather}

where n is the number of observations and u is the number of unknowns in the equation. With the IGG-III function, the equivalent weighting matrix is constructed iteratively. The maximum standardised residual value of observations is taken into account. Equation (15) is used to reweight this observation. The processing steps are repeated with the new matrix. The procedure is repeated until there is no change in the standardised residual values obtained in the two iterations.

3. Experimental design

3.1 Data description

To investigate the SPP performances of all combinations of GPS, GLONASS, Galileo and BeiDou navigation systems (BDS)-3, 16 stations of the multi-GNSS experimental (MGEX) monitoring network were selected. The five-day observation files, between Day of Year (DOY) 137–141 (17–21 May 2021) with a sample interval of 30 s, were obtained from the Crustal Dynamics Data Information System (CDDIS) in Receiver-Independent Exchange (RINEX) format (available at: https://cddis.nasa.gov/archive/gnss/). The main purpose for the selection of these stations is to measure worldwide performance of GNSS systems. The geographic distribution of these MGEX stations is shown in Figure 1.

Figure 1. Geographic distribution of 16 selected IGS MGEX stations in this study

The main objective of this research is to assess the real-time positioning performance, using pseudo-range observations, with external orbit and clock products offered by iGMAS. iGMAS, which is similar to the International GNSS Service (IGS) in terms of mission, supplies these products free of charge, calculating them with its global monitoring station's data. Zhu et al. (2022) and (Zhou et al. (2022) can be examined for more information about iGMAS's tracking network, analysis centres, and mission. Both orbit and clock products are accessible with a sampling interval of 15 min. The ultra-rapid (so-called super-fast) orbit and clock offset files provided by iGMAS for six-hourly include GPS, GLONASS, Galileo and BDS (available at: http://www.igmas.org/Product/). Thus, these products enable near-real-time and real-time applications with quad constellation integration. During the processing stage, iGMAS products were sorted through BDS System Time (BDT), whereas the observations referred to the GPS Time (GPST). The time difference between BDT and GPST is 14 s (iGMAS, 2022). Resolving this inconsistency, requires getting them into alignment with the chosen reference time (GPST) system. For the improvement of SPP results, the differential code bias (DCB) values should be applied to the pseudo-range observations. The Bias Solution Independent Exchange Format (Bias-SINEX) file contains the bias values, which vary depending on the satellites and type of code measurements, provided by the Centre for Orbit Determination in Europe (CODE) on a monthly basis. In our study, C1C code observations at GPS L1 frequency at BRST, GANP, KRGG, BOAV, MAR7 and JFNG stations were turned into C1W observations using differential code biases (C1C–C1W). The current antenna file (igs14 2163.atx) was applied to convert the satellite coordinates to their phase centres. Furthermore, the precise coordinates of the stations provided by IGS weekly were considered as the ground truth. Accuracy performance analyses were carried out by comparing the ground truth coordinates with the epoch-by-epoch coordinates acquired from the SPP solution.

3.2 Processing strategies

Within the context of the study, MATLAB-based in-house code was developed to evaluate the SPP solution of single, dual, triple and quad combinations of GPS, GLONASS, Galileo and BDS-3 constellations. Figure 2 shows the flowchart of the code.

Figure 2. Flowchart of SPP processing code. The code includes three components: data preparation component, data processing component and output component (SPP: single-point positioning; MAD: mean absolute deviation; IGG-III: institute of geodesy and geophysics-III; PDOP: position dilution of precision)

The orbit and clock offsets of BDS satellites are not included in the ultra-rapid products provided by IGS ACs except for WHU and CNES. However, the iGMAS ultra-rapid products, which incorporate orbit and clock offsets of BDS-3 satellites, give the advantage of a solution for real-time applications with the integration of the quad constellation. In addition, the performance of BDS-3 with both single and other navigation system combinations can also be evaluated. Furthermore, the developed program is capable of performing SPP solutions with both the iGMAS and IGS products. The same observation data was analysed again using the GeoForschungs Zentrum (GFZ) rapid orbit and clock offsets in order to compare the accuracy performance of all combinations. Table 1 describes the methodologies used in these evaluation stages.

Table 1. Processing settings for developed SPP-based program

4. Results

4.1 Assessment of iGMAS ultra-rapid orbit and clock offset products

The performance of iGMAS ultra-rapid (predicted part) products was compared with the Centre for Orbit Determination in Europe (CODE) final precise products throughout the solution days. Figure 3 shows the orbit and clock offset accuracies of the available satellites. While the RMS values were considered for orbital accuracies, standard deviation (STD) values were calculated for clock offsets due to the different clock datums of products. The average 3D orbit accuracy of available satellites was estimated as 0 ⋅ 07, 0 ⋅ 15, 0 ⋅ 13 and 0 ⋅ 25 m for GPS, GLONASS, Galileo and BDS-3, respectively. BDS-3 is distinct from the others because of the low-orbit accuracy of IGSO satellites. The clock offset STDs for GPS, GLONASS, Galileo and BDS-3 were found as 2 ⋅ 03, 2 ⋅ 27, 1 ⋅ 92 and 1 ⋅ 85 ns, respectively. In addition, note that the situations with no orbit and clock offset information, no clock offset information while the orbit information is available, or vice versa, have been observed in iGMAS products for some satellites. Providing the products for all satellites in the future will increase the effectiveness of the iGMAS in positioning.

Figure 3. The average RMS values for the orbit and the average STD values for the clock offset according to the final products produced by IGS CODE for available satellites

4.2 Availability and PDOP analyses

The number of the available satellite and position dilution of precision (PDOP) values were obtained for all combinations of GPS, GLONASS, Galileo and BDS-3 in SPP solutions for 16 IGS stations with iGMAS ultra-rapid products. Positioning performance was negatively affected by the satellite availability of the iGMAS ultra-rapid products. To be more specific, the satellites associated with the observed session might not be present in the precise products. Therefore, the satellites not having clock offset and/or orbit information were ignored in the solution. In the developed code, the solution outcomes were examined according to various cases. In this context, the SPP solutions are possible for cases containing at least 4, 5, 6 and 7 satellites on epochs for single, dual, triple and quad systems, respectively. In addition, if the PDOP value was less than 20 in any epoch, the solution were accepted as a valid result.

Figure 4 depicts the maximum, mean and minimum number of satellites employed in the solution at the stations, as well as the PDOP values. The average number of satellites in GPS, GLONASS, Galileo and BDS-3 were 9–10, 5–6, 7–8 and 5–6, respectively, while the PDOPs of these systems were 1 ⋅ 81, 3 ⋅ 30, 2 ⋅ 30 and 3 ⋅ 11, respectively. The results showed that GLONASS and BDS-3 performed worse than other systems in terms of the minimum satellite and maximum PDOP criteria. It is an expected result that this issue affects the positioning performance of GLONASS and BDS-3, which will be explained in more detail in the next section. On the other hand, it was observed that Galileo was the most stable system after GPS in terms of PDOP values and the number of satellites.

Figure 4. Maximum, average and minimum number of satellites and PDOP values for all combinations

Within the dual constellations, on average, the PDOP values ranged between 1 ⋅ 29-and 1 ⋅ 91. Thus, the dual-system solutions provided significant advantages in positioning accuracy compared with single-system solutions. However, note that GLONASS- BDS-3 integration results are not very compatible with the other dual constellations. Redundancy in SPP solutions for a single epoch increased dramatically, with an average of 17–23 and 27–28 visible satellites in all triple- and quad-constellation integrations. The PDOP values were at an extremely ideal level in these combinations. Furthermore, the narrower range of maximum and minimum PDOP values indicates more reliable positioning. The positioning accuracy, reliability and redundancy provided by the quad constellation were especially remarkable.

4.3 Single-system SPP performance

The positioning accuracy of navigation systems was assessed for the solution of single-system cases (GPS-only, GLONASS-only, Galileo-only and BDS-3-only) using the developed SPP code. However, as indicated in the previous section, the number of visible satellites is a critical condition for all positioning techniques, not just SPP. When there are insufficient satellites, or even when there are enough satellites, the orbit–clock products of some of these satellites are unavailable, and position acquisition is impossible. Theoretically, in single-system situations, orbit and clock information of at least four satellites must be available so that a solution can be obtained. However, due to satellite health or other factors, some satellites’ orbit and clock products may not be available. This circumstance particularly occurs in ultra-rapid (in the predicted-part) products. For this reason, some epochs could not be solved and were marked as unresolved. In other words, the success of providing solutions throughout all epochs is referred to as the service rate for navigation systems, and it is defined by valid solutions. Figure 5 shows the average service rate of independent navigation systems in five-day solutions employing both iGMAS ultra-rapid and GFZ rapid products. The results reveal that the service rates of GPS and Galileo are very similar. The GPS service rate for both products is more than 98% for all stations. Ultra-rapid products in Galileo solutions have a service rate of more than 94% at all stations, except for the USUD station, which has an 88% service rate. GLONASS had the worst service rate among all navigation systems. While the average GLONASS service rate for the iGMAS products was 76%, the service rate for all stations increased with GFZ products and reached an average of 95%. Furthermore, the GLONASS service rate is slightly better in the regions close to the pole due to its high orbital inclination. Except for the BOAV station, the BDS-3 showed an average service rate performance of 87% with ultra-rapid products. Excluding the same station as the GFZ rapid products, the service rate has exceeded 98%. The service rate of this station has not changed and has remained at 62%. This situation is considered to be caused by geographical location. When the service rates for GLONASS and BDS-3 are assessed, it is clear that the quality of iGMAS ultra-rapid products for the systems should be upgraded.

Figure 5. Service rate achieved by iGMAS ultra-rapid and GFZ rapid products in the event of a GPS, GLONASS, Galileo and BDS-3 single-system for 16 MGEX stations

Focusing on the individual GNSS SPP (GPS: G, GLONASS: R, Galileo: E, BDS-3: C) performances, Figure 6 shows the distribution of positioning errors for the north, east and up components for 16 stations. The distributions, containing a total of 230,400 epoch (the five-day) SPP solutions, are expressed as probability percentages. In addition, the average error and root mean square (RMS) statics were obtained for each component (Figure 6). When the RMS values for the north, east and up components of the generated error distributions are analysed in detail, Galileo emerges as the single-system solution with the best positioning performance, with RMS values of 0 ⋅ 56, 0 ⋅ 53 and 1 ⋅ 23 m, respectively. All of Galileo's error components have lower RMS values than the others. Thus, it can be said that the quality of the Galileo signals used made a significant contribution to the positioning accuracy. This resulted in agreement with the findings reported in a study with multi-GNSS SPP (Zhang and Pan, 2022). Then, the GPS showed higher positioning accuracy than GLONASS and BDS-3, with the RMS statistic being 0 ⋅ 75, 0 ⋅ 62 and 1 ⋅ 54 m in the north, east, and up components, respectively. GLONASS had the worst single-system positioning performance with the RMS statistics of 1 ⋅ 79, 1 ⋅ 86,and 4 ⋅ 04 m in the three components, respectively. BDS-3 offered better positioning performance than GLONASS, with the RMS statistics for three components being 1 ⋅ 62, 1 ⋅ 49 and 3 ⋅ 37 m, respectively. GPS and Galileo outperformed the positioning accuracy of GLONASS and BDS-3. Although BDS-3 can compete with GPS and Galileo in terms of service rate and PDOP values with the iGMAS ultra-rapid products, it performed poorly for positioning accuracy. In the mean error, the horizontal components have a maximum error of 7 cm and the vertical component has a maximum error of 27 cm. These mean errors, especially in single solutions, provide evidence that systematic errors are well modelled. It should be also noted that in a single-system solution, low redundancy might result in a potential outlier going undetected, causing problems with the solutions. This occurred, particularly in some GLONASS and BDS-3 solutions. The $\textrm{MAD}$ method prevented this issue, whereas the IGG-III robust approach did not operate properly in epochs with low redundancy.

Figure 6. Error distributions of epoch-by-epoch SPP solutions in the north, east and up components with iGMAS ultra-rapid products for GPS-only, GLONASS-only, Galileo-only and BDS-3-only

4.4 Dual-System SPP performance

Figure 7 demonstrates the distributions of positioning errors, RMS statistics and mean statistics for all dual combinations of the navigation systems utilised in the study. First of all, dual-constellation combination solutions considerably improved the positioning accuracy as compared with single-GNSS systems. The best dual-combination performance was produced with the GPS/Galileo solution. RMS statistics of the north, east and up components were 0 ⋅ 46, 0 ⋅ 41 and 1 ⋅ 00 m, respectively. Furthermore, as a significant finding, the Galileo-based combinations outperformed the GPS-based combinations in terms of all metrics. Considering the GLONASS and BDS-3 combinations of GPS and Galileo, BDS-3 further improved positioning accuracy according to GLONASS. BDS-3 enhanced the 3D RMS errors of GPS-only, Galileo-only and GLONASS-only solutions by about 12%, 12% and 44%, respectively. Among the dual combinations, GLONASS-BDS-3 had the worst outcomes. In addition, combining BDS-3 and GLONASS did not yield the same quality results as using GPS and Galileo alone.

Figure 7. Error distributions of epoch-by-epoch SPP solutions in the north, east and up components with iGMAS ultra-rapid products for GPS/GLONASS, GPS/Galileo, GPS/BDS-3, GLONASS/Galileo, GLONASS/BDS-3 and Galileo/BDS-3

4.5 Triple-system and quad-system SPP performance

The error distributions, RMS and error values for triple and quad combinations are shown in Figure 8. The quad combination produced higher accurate results than all the combinations from the five-day epoch-by-epoch solution of the 16 stations dataset. The RMS values of 0 ⋅ 54 m horizontally and 0 ⋅ 89 m vertically were estimated for the quad constellation. The benefit of employing the quad constellation in real-time applications, with its RMS value of 1 ⋅ 02 m in 3D, should be highlighted. In addition, note that there is no significant difference between the quad constellation results and the GPS/GLONASS/Galileo and GPS/Galileo/BDS-3 results. Because they have almost nearly the same accuracy as the quad constellation, these triple combinations can be employed in instances where quad combinations are not attainable. However, the accuracy of the GLONASS/Galileo/BDS-3 was better than the GPS/GLONASS/BDS-3 results. When GLONASS and BDS-3 were used simultaneously in the triple combination, the accuracy was considerably superior to that of its dual combination.

Figure 8. Error distributions of epoch-by-epoch SPP solutions in the north, east and up components with iGMAS ultra-rapid products for GPS/GLONASS/Galileo, GPS/GLONASS/BDS-3/, GPS/Galileo/BDS-3, GLONASS/Galileo/BDS-3 and GPS/Galileo/GLONASS/BDS-3

4.6 iGMAS ultra-rapid products versus GFZ rapid products

SPP solutions were performed with the same dataset also utilising GFZ rapid orbit and clock products to assess the quality of the iGMAS ultra-rapid orbit and clock products. Table 2 details the five-day mean RMS statistics for all combinations of each station based on both GFZ and iGMAS products. The most important point to highlight is that while the GLONASS-only and BDS-3-only results in the two different solutions differed dramatically, the GPS-only and Galileo-only results did not differ significantly. Additionally, the performance of GLONASS/BDS-3 combined improved by 35% with GFZ products. However, when GPS/Galileo was taken into consideration, it showed by 9% improvement. It can be concluded that GLONASS and BDS-3 orbit and clock products produced by the iGMAS lagged behind GPS and Galileo in terms of availability and accuracy. The positioning performance of the systems, both individually and in combination with other systems, is predicted to improve if the orbit and clock products of BDS-3 and GLONASS are enhanced. In triple-system solutions, the accuracies obtained with GFZ and iGMAS products were consistent with each other. Nevertheless, the GRC combination including GLONASS and BDS-3 produced poorer results. Finally, the 3D accuracy obtained with real-time positioning in the quad-constellation integration was only 0 ⋅ 13 m less than the result obtained with the rapid products. Namely, precise navigation needs can be met using the combination of GREC with the ultra-rapid products provided for GPS, GLONASS, Galileo and BDS-3.

Table 2. Five-day average 3D RMS statistics based on iGMAS ultra-rapid and GFZ rapid products for all stations (unit: m)

5. Conclusion

With the advancement of GNSS technology, space-based positioning has become more common in real-time applications, particularly navigation, guiding and surveying. Currently, there are four independent global systems, namely GPS, GLONASS, Galileo and BDS. For precise positioning, many analysis centres (ACs) provide basic information about the satellites as ultra-rapid products to real-time users going through the IGS MGEX network. Parallel to IGS, iGMAS offers ultra-rapid products for the four global constellations using its own network computed by 12 ACs. This study aimed to assess the performance of iGMAS ultra-rapid products in navigations problems. To the fulfilment of the objective, SPP solutions were performed for all 15 combinations (single, dual, triple and quad) of the constellations using a five-day dataset of 16 MGEX stations. In this context, an in-house code was developed in MATLAB for the SPP solutions. The MAD approach was successful in removing gross errors in the SPP solutions. The approach also supported reweighting robust-technique IGG-III performance. Additionally, all datasets were also processed with GFZ rapid products to assess the availability of the iGMAS ultra-rapid products, service rate and positioning accuracy. In comparison to GPS and Galileo, the results indicated that GLONASS and BDS-3 had poorer service rates. This situation can be interpreted as being dependent on the absence of some satellites data in ultra-rapid products. When using the ultra-rapid products against the rapid products, on average, the service rate decreased from 95% to 76% for GLONASS and 97% to 86% for BDS-3. Moreover, there were no noticeable degradations in the service rates for GPS and Galileo when employing ultra-rapid products.

For accuracy validation of the iGMAS products, the RMS values were calculated using all epoch-wise solution results. The single-system solutions showed that Galileo produced the best results, with RMS values of 0 ⋅ 56, 0 ⋅ 53 and 1 ⋅ 23 m in the north, east and up components, respectively. The accuracy achieved with Galileo-only was even better than some dual combinations. This can be explained by Galileo's observations being less sensitive to the multipath effect, including less noise, and being less influenced by non-modellable errors. The worst solutions were generated using the GLONASS, with its RMS values of 1 ⋅ 79, 1 ⋅ 86 and 4 ⋅ 04 m. The dual constellation solutions demonstrated that combinations with Galileo produced better results than did the GR, GC and RC solutions. In particular, the RC solution differed from the other dual solutions negatively. In the triple constellation, the other results, except for the GRC, varied in the range 0 ⋅ 43–0 ⋅ 45 m, 0 ⋅ 38–0 ⋅ 44 m and 0 ⋅ 93–1 ⋅ 02 m in the north, east and up components, respectively. The quad-solution results had the lowest RMS values of 0 ⋅ 40, 0 ⋅ 37 and 0 ⋅ 89 m in the three components. The combination results with respect to their RMS values from the worst to the best can be listed as R, C, RC, G, GR, GC, GRC, E, RE, EC, REC, GE, GRE, GEC and GREC. Results produced in this study indicated that Galileo and its combinations exhibited remarkable performance.

Finally, the accuracy level obtained in triple and quad combination (approximately 0 ⋅ 65 m horizontal, 1 m vertical component) with the proposed algorithm without using any augmentation systems can meet the requirements of many applications, such as civil aviation, smart agriculture practices, ship navigation, pedestrian and vehicle tracking, autonomous systems such as unmanned aerial vehicles (UAV), and for some road and railway applications. Positioning accuracy and reliability can be increased by expanding iGMAS's network and enhancing the availability of ultra-rapid products. For future study, the developed SPP algorithm is planned to be tested in harsh environments, such as urban canyons and forest areas, with multi-GNSS.