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, Reference Kouba and Héroux2001). 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, Reference Kouba and Héroux2001; Zhou et al., Reference Zhou, Dong, Li, Li and Schuh2019). 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., Reference Lou, Zheng, Gu, Wang, Guo and Feng2016; Liu et al., Reference Liu, Jiang, Chen, Zhao, Huo, Huang and Chen2019). 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, Reference Niell1996; Leandro et al., Reference Leandro, Langley and Santos2008). 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, Reference Koch1999)

(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., Reference Pan, Cai, Santerre and Zhang2017)

(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, Reference Wanninger2012).

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, Reference Cai and Gao2013). 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, Reference Zhang and Pan2022). 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., Reference Angrisano, Gaglione, Crocetto and Vultaggio2020). 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, Reference Rousseeuw and Leroy1987; Hekimoglu et al., Reference Hekimoglu, Erdogan, Soycan and Durdag2014)

(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, Reference Pearson2001; Chiang et al., Reference Chiang, Pell and Seasholtz2003). 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., Reference Yang, Song and Xu2002; Guo and Zhang, Reference Guo and Zhang2014)

(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.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. (Reference Zhou, Cai, Chen, Jiao, He and Yang2022) and (Zhou et al. (Reference Zhou, Cai, Chen, Jiao, He and Yang2022) 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