Investigation of Tightly Combined Single-Frequency and Single-Epoch Precise Positioning Using Multi-GNSS Data

Abstract

The loose combination (LC) and the tight combination (TC) are two different models in the combined processing of four global navigation satellite systems (GNSSs). The former is easy to implement but may be unusable with few satellites, while the latter should cope with the inter-system bias (ISB) and is applicable for few tracked satellites. Furthermore, in both models, the inter-frequency bias (IFB) in the GLObal NAvigation Satellite System (GLONASS) system should also be removed. In this study, we aimed to investigate the performance difference of ambiguity resolution and position estimation between these two models simultaneously using the single-frequency data of all four systems (GPS + GLONASS + Galileo + BeiDou Navigation Satellite System (BDS)) in three different environments, i.e., in an open area, with surrounding high buildings, and under a block of high buildings. For this purpose, we first provide the definition of ISB and IFB from the perspective of the hardware delays, and then propose practical algorithms to estimate the IFB rate and ISB. Thereafter, a comprehensive performance comparison was made between the TC and LC models. Experiments were conducted to simulate the above three observation environments: the typical situation and situations suffering from signal obstruction with high elevation angles and limited azimuths, respectively. The results show that in a typical situation, the TC and LC models achieve a similar performance. However, when the satellite signals are severely obstructed and few satellites are tracked, the float solution and ambiguity fixing rates in the LC model are dramatically decreased, while in the TC model, there are only minor declines and the difference in the ambiguity fixing rates can be as large as 30%. The correctly fixed ambiguity rates in the TC model also had an improvement of around 10%. Once the ambiguity was fixed, both models achieved a similar positioning accuracy.

1. Introduction

The Global Navigation Satellite System (GNSS) has entered into a new era in recent decades with four different systems operating simultaneously, i.e., American Global Positioning System (GPS), Russian GLObal NAvigation Satellite System (GLONASS), European Galileo, and Chinese BeiDou Navigation Satellite System (BDS), forming a joint multi-GNSS system. The increasing number of visible satellites have been confirmed to improve the satellite availability, convergence time, and reliability [1,2], and will bring great benefits to urban navigation services. Since satellite signals may be obstructed by the surrounding huge buildings or a block of high buildings, the multi-GNSS system can provide the potential ability to resist the risk of satellite deficiency and retain a satisfactory positioning performance.

How to integrate the observation data of different systems is a crucial issue toward realizing multi-GNSS precise positioning. In the case of popular relative positioning, such as real-time kinematic (RTK), there are two different models to realize the integration [3]. The first is to individually choose pivot satellites for different systems and undertake differencing between satellites within each system; this method is called the loose combination (LC). The other method is to choose only one pivot satellite for all systems and undertake differencing across different systems; this method is called the tight combination (TC) [4].

The data processing for hybrid systems began from the integration of GPS and GLONASS precise positioning. Due to the frequency division multiple access (FDMA) technique used in the current GLONASS system, the frequencies in the same frequency band for different satellites are different [5]. Dodson et al. [6], Wang et al. [7], Dai et al. [8], and Al-Shaery et al. [9] discussed the model of this hybrid system and the ambiguity resolution (AR) issues, and they focused on the techniques used to eliminate the relative receiver clock errors from GLONASS data. Considering that different hardware biases in GLONASS receiving channels cannot be removed by between-satellite differencing, the inter-frequency bias (IFB) was introduced and analyzed [10,11], and lots of attention has been paid to the processing of IFB. Considering the long-term stability, the IFB can be treated as a constant for a certain pair of receivers, and there are two ways to manage it. One way is to pre-estimate and calibrate it in the positioning process [12,13,14], and the other is to directly estimate it as an unknown parameter in the processing procedure [15,16,17]. All these studies for the GPS + GLONASS system used the LC model, and similar LC models can also be found in other hybrid systems, such as GPS + Galileo [3,18] and GPS + BDS [19,20].

For the TC model, the difference in hardware delays in different systems leads to the inter-system bias (ISB) [21]. However, at early stages, the ISB was not considered and this model was only applied in hybrid systems with overlapping frequencies, such as GPS + Galileo L1/E1 + L5/E5a [4,22]. Odijk and Teunissen [23] and Paziewski et al. [24] identified that for heterogeneous receivers, the ISB cannot be ignored even if overlapping frequencies were used. The subsequent studies put emphasis on the ISB estimation in hybrid systems without overlapping frequencies [25], and the RTK application of TC models can be found in GPS + BDS [26] and other hybrid systems using the code division multiple access (CDMA) technique [27,28]. The only contribution for the TC model between CDMA and FDMA systems was made by Gao et al. [29], where only GPS + GLONASS was used. Following their work, in this paper we will further make an investigation of tightly combined RTK positioning using the single-frequency signals of GPS L1 (1575.42 MHz), GLONASS L1 (1602+k·9/16 MHz) (k is the channel number), BDS B1I (1561.098 MHz), and Galileo E1 (1575.42 MHz). To achieve this goal, we will first propose algorithms to estimate the IFB in the GLONASS system and ISB between GPS and the other three systems.

Furthermore, the performance of tight combined single-frequency RTK positioning using multi-GNSS data should also be comprehensively investigated. Paziewski and Wielgosz [22] discussed the performance of a tightly combined GPS + Galileo under clear-sky satellite visibility (15° elevation mask) and obstructed satellite visibility (30° elevation mask). They also compared the performance using a TC model and an LC model in a hybrid GPS + Galileo system [30]. The results indicated that the TC model shows significant improvement in the obstructed environment compared with single systems, and both TC and LC models gave comparable results regarding the ambiguity resolution and coordinate domains. Combining these two experiments, in this paper, we will investigate the performance of the TC model by comparing it with that of LC model under three different environments. In this way, the advantage of the TC model using multi-GNSS data becomes clearer.

For these two reasons, this paper has launched an investigation of tightly combined single-frequency and single-epoch precise positioning using all four GNSSs. We first derive the double differenced mathematical model of tightly combined single-frequency multi-GNSS observations and achieve the definition of IFB and ISB from the perspective of hardware biases. In order to eliminate the IFB and ISB, we propose an easy-to-implement strategy to manage each of them. After the removal of hardware delays, we give the results of simulating three different measurement scenarios in the following experiments, and discuss the corresponding AR and positioning performance of tightly combined single-frequency single epoch multi-GNSS data and its improvements compared with the loose combined data. In the last section, we make some concluding remarks for users of multi-GNSS data processing.

2. Mathematical Model

The original single-frequency pseudorange and carrier phase observation equations for a certain satellite system M, which can be rewritten as G for GPS, R for GLONASS, C for BDS, and E for Galileo, can be expressed as [31]:

$$\{\begin{array}{c}{P}_r^{s,M}={\rho }_r^{s,M}+c(d{t}_r-d{t}^{s,M})+T_r^{s,M}+I_r^{s,M}+B_r^{s,M}-B^{s,M}+{\epsilon }_P,\\ {\varphi }_r^{s.M}={\rho }_r^{s,M}+c(d{t}_r-d{t}^{s,M})+T_r^{s,M}-I_r^{s,M}-{\lambda }^{s,M}\left(N_r^{s,M}-b_r^{s,M}+b^{s,M}\right)+{\epsilon }_{\varphi },\end{array}$$

(1)

where P and $\varphi$ are the original code and phase observables in meters. The superscript s and the subscript r are the satellite and the receiver, and $\rho$ is the geometric distance between them. The symbols c, $d{t}_r$, and $d{t}^{s,M}$ are the speed of light in a vacuum, the receiver clock bias, and the satellite clock bias, respectively. T and I are the tropospheric and ionospheric delays, respectively. $B_r^{s,M}$ and $B^{s,M}$ are the receiver and satellite hardware code delays in meters, respectively, and $b_r^{s,M}$ and $b^{s,M}$ are the combination of the hardware phase delays and the initial phase offsets for the receiver and the satellite in cycles, respectively. $\lambda$ is the signal wavelength, and N is the unknown integer ambiguity. ${\epsilon }_P$ and ${\epsilon }_{\varphi }$. are the residuals for the code and the phase, respectively.

After differencing between the two receivers r and b, we can obtain the single-differenced (SD) pseudorange and carrier phase observation equations as follows:

$$\{\begin{array}{c}\Delta P_{rb}^{s,M}=\Delta {\rho }_{rb}^{s,M}+c·\Delta d{t}_{rb}+\Delta T_{rb}^{s,M}+\Delta I_{rb}^{s,M}+\Delta B_{rb}^{s,M}+\Delta {\epsilon }_P,\\ \Delta {\varphi }_{rb}^{s.M}=\Delta {\rho }_{rb}^{s,M}+c·\Delta d{t}_{rb}+\Delta T_{rb}^{s,M}-\Delta I_{rb}^{s,M}-{\lambda }^{s,M}\left(\Delta N_{rb}^{s,M}-\Delta b_{rb}^{s,M}\right)+\Delta {\epsilon }_{\varphi },\end{array}$$

(2)

where the symbol $\Delta (·)$ is the single-differencing operator between two receivers. The satellite-dependent terms, i.e., the satellite clock bias, the satellite code, and phase hardware delays, and the satellite initial phase offsets, are all removed. The atmospheric delay and the receiver-dependent terms are reduced. For short baselines, the atmospheric delays can be ignored due to their strong spatial correlation and Equation (2) can be simplified to:

$$\{\begin{array}{c}\Delta P_{rb}^{s,M}=\Delta {\rho }_{rb}^{s,M}+c·\Delta d{t}_{rb}+\Delta B_{rb}^{s,M}+\Delta {\epsilon }_P,\\ \Delta {\varphi }_{rb}^{s.M}=\Delta {\rho }_{rb}^{s,M}+c·\Delta d{t}_{rb}-{\lambda }^{s,M}\left(\Delta N_{rb}^{s,M}-\Delta b_{rb}^{s,M}\right)+\Delta {\epsilon }_{\varphi }.\end{array}$$

(3)

For medium or long baselines, the residuals of the ionospheric and tropospheric delay after double-differencing cannot be neglected and should be corrected in advance using external information, such as regional correction products, or estimated together with the baseline and ambiguity using the ionosphere and troposphere weighted model [30].

To manage the relative receiver clock bias, a further differencing between two satellites is made. Here, both LC and TC models can be used to eliminate the relative receiver clock bias. Before this operation, for a joint multi-GNSS system, the time systems and coordinate systems should be aligned. Usually, the GPS Time system and World Geodetic System 1984 are adopted as the reference, and the biases of the time system and the coordinate system from other systems to the GPS should be corrected.

In the loose combined model, the receiver code and phase hardware delays and the receiver initial phase offsets are also canceled, except those for the GLONASS system because of its IFB. As a result, the loose combined double-differenced (DD) pseudorange and carrier phase observation equations for satellites p and q in a single GPS/BDS/Galileo system can be expressed as:

$$\{\begin{array}{c}\Delta \nabla P_{rb}^{pq, M}=\Delta \nabla {\rho }_{rb}^{pq,M}+\Delta \nabla {\epsilon }_P,\\ \Delta \nabla {\varphi }_{rb}^{pq.M}=\Delta \nabla {\rho }_{rb}^{pq,M}-{\lambda }^M·\Delta \nabla N_{rb}^{pq,M}+\Delta \nabla {\epsilon }_{\varphi },\end{array}$$

(4)

where p is the reference satellite. For different observation types, the IFB in the GLONASS system can be divided into the inter-frequency code bias (IFCB) and inter-frequency phase bias (IFPB), and the corresponding observation equations can be expressed as:

$$\{\begin{array}{c}\Delta \nabla P_{rb}^{pq, R}=\Delta \nabla {\rho }_{rb}^{pq, R}+IFC{B}_{rb}^{pq}+\Delta \nabla {\epsilon }_P,\\ \Delta \nabla {\varphi }_{rb}^{pq.R}=\Delta \nabla {\rho }_{rb}^{pq, R}-\left({\lambda }^{q,R}·\Delta N_{rb}^{q,R}-{\lambda }^{p,R}·\Delta N_{rb}^{p,R}\right)+IFP{B}_{rb}^{pq}+\Delta \nabla {\epsilon }_{\varphi },\end{array}$$

(5)

where the symbol $\Delta \nabla (·)$ is the double-differencing operator, and:

$$\{\begin{array}{c}IFC{B}_{rb}^{pq}=\Delta B_{rb}^{q,R}-\Delta B_{rb}^{p,R},\\ IFP{B}_{rb}^{pq}={\lambda }^{q,R}·\Delta b_{rb}^{q,R}-{\lambda }^{p,R}·\Delta b_{rb}^{p,R}.\end{array}$$

(6)

In the tightly combined model, when we undertake differencing across different systems, neither the code and phase hardware delays nor the receiver initial phase offsets can be eliminated. For the code observation equation, the inter-system code bias (ISCB) arises; for the carrier phase observation equation, the differenced receiver initial phase offsets are absorbed by the inter-system phase bias (ISPB). On the other hand, if the GLONASS system is involved, the IFCB and IFPB must come out. Taking GPS and GLONASS systems as example, when a GPS satellite is chosen as the reference satellite, from Equation (3), the tightly combined observation equations can be derived to produce:

$$\{\begin{array}{c}\Delta \nabla P_{rb}^{GR}=\Delta \nabla {\rho }_{rb}^{GR}+ISC{B}_{rb}^{GR}+IFC{B}_{rb}^R+\Delta \nabla {\epsilon }_P,\\ \Delta \nabla {\varphi }_{rb}^{GR}=\Delta \nabla {\rho }_{rb}^{GR}-\left({\lambda }^R·\Delta N_{rb}^R-{\lambda }^G·\Delta N_{rb}^G\right)+ISP{B}_{rb}^{GR}+IFP{B}_{rb}^R+\Delta \nabla {\epsilon }_{\varphi },\end{array}$$

(7)

where

$$\{\begin{array}{ll}ISC{B}_{rb}^{GR}& +IFC{B}_{rb}^R=\Delta B_{rb}^R-\Delta B_{rb}^G=(\Delta B_{rb}^{R0}-\Delta B_{rb}^G)+(\Delta B_{rb}^R-\Delta B_{rb}^{R0})\\ & ISP{B}_{rb}^{GR}+IFP{B}_{rb}^R={\lambda }^R·\Delta b_{rb}^R-{\lambda }^G·\Delta b_{rb}^G\\ & =({\lambda }^{R0}·\Delta b_{rb}^{R0}-{\lambda }^G·\Delta b_{rb}^G)+({\lambda }^R·\Delta b_{rb}^R-{\lambda }^{R0}·\Delta b_{rb}^{R0}).\end{array}$$

(8)

Herein, $R0$ refers to a virtual signal of the constant central frequency in the GLONASS system with a channel number of zero, i.e., 1602 MHz for L1 and 1246 MHz for L2. From Equation (8), we can obtain the definition expressions for ISB and IFB as follows:

$$\{\begin{array}{c}ISC{B}_{rb}^{GR}=\Delta B_{rb}^{R0}-\Delta B_{rb}^G,\\ ISP{B}_{rb}^{GR}={\lambda }^{R0}·\Delta b_{rb}^{R0}-{\lambda }^G·\Delta b_{rb}^G,\\ IFC{B}_{rb}^R=\Delta B_{rb}^R-\Delta B_{rb}^{R0},\\ IFP{B}_{rb}^R={\lambda }^R·\Delta b_{rb}^R-{\lambda }^{R0}·\Delta b_{rb}^{R0}.\end{array}$$

(9)

This indicates that ISCB and ISPB describe the hardware delay difference between the GLONASS system with the central frequency and GPS system, and IFCB and IFPB describe the hardware delay difference of the GLONASS system between the actual frequency and the central frequency.

For a combined system without GLONASS, the terms IFCB and IFPB can be left out and only ISCB and ISPB need to be considered. As a result, in the tight combination of four systems, we highly suggest the users select a non-GLONASS satellite as the pivot satellite for easy processing.

Derived from the standard deviation of the original pseudorange and carrier phase, and taking the correlation into consideration, we can determine the corresponding stochastic model. For the loose combination, the weight matrix of each system can be determined separately in a traditional way [7], and the weight matrix of the combined multi-GNSS system can be expressed as $P=\mathrm{diag}\left(n_G{P}_G,n_R{P}_R,n_C{P}_C,n_E{P}_E\right)$, where $\left(n_G,n_R,n_C,n_E\right)$ is the ratio of the weights for different systems; in this paper, the four systems are equally weighted, i.e., $n_G=n_R=n_C=n_E=1$.

For the tight combination, given the multi-GNSS original code variance matrix D_P, carrier phase variance matrix D_φ, and the DD transition matrix H, the final stochastic model can be expressed as follows:

$$\mathit{P}={\mathit{D}}_{\Delta \nabla }^{-1}={\left[\begin{array}{cc}\mathit{H}{\mathit{D}}_P{\mathit{H}}^T& 0\\ 0& \mathit{H}{\mathit{D}}_{\phi }{\mathit{H}}^T\end{array}\right]}^{-1},$$

(10)

where

$${\mathit{D}}_P=\mathrm{diag}\left({\mathit{D}}_P^G,{\mathit{D}}_P^C,{\mathit{D}}_P^R,{\mathit{D}}_P^E\right), {\mathit{D}}_{\phi }=\mathrm{diag}\left({\mathit{D}}_{\phi }^G,{\mathit{D}}_{\phi }^C,{\mathit{D}}_{\phi }^R,{\mathit{D}}_{\phi }^E\right),$$

(11)

and

$$\mathit{H}=\left[\begin{array}{ccccccc}1& -1& -1& 1& \cdots & 0& 0\\ 1& -1& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& -1& 0& 0& \cdots & -1& 1\end{array}\right].$$

(12)

Herein, it is assumed that the first two columns denote the reference satellite for the base station and the rover station. In this way, the definition of ISB and IFB is functionally expressed, and the weight of the multi-constellation relative positioning is determined.

3. IFB and ISB Processing

The two kinds of relative hardware delays, i.e., IFB and ISB, are mathematically defined in Equation (9). Considering the temporal stability of IFB [11] and ISB [32], they can be estimated in the preprocessing step before the determination of the ambiguity and the position. Taking six pairs of available multi-GNSS data as an example, we propose a IFPB calibration and the subsequent ISB estimation method. Before the estimation of IFB and ISB using the above algorithms, the satellites whose elevation angles were lower than 10° were excluded, and a further outlier detection step was also applied.

The used data were all collected on a whole day in 2019 with 30-s intervals. They could be divided into two categories: the first three baselines were provided by the GNSS Research Centre of Curtin University in Australia and were all ultra-short baselines, while the other three baselines were collected from Hong Kong SatRef and were all short baselines. Herein, the ultra-short and short baselines are defined when their lengths are shorter than 100 m and 15 km, respectively. The detailed data configuration is listed in

Table 1. Configuration of used static baselines.

Tag	Baseline Name	Receiver Type		Date	Length (m)
		Base Station	Rover Station
A	CUT2-CUTB	TRIMBLE NETR9	TRIMBLE NETR9	1 January	4.27
B	CUT3-CUCC	JAVAD TRE_G3T	JAVAD TRE_G3T	1 March	7.99
C	CUT3-CUTA	JAVAD TRE_G3T	TRIMBLE NETR9	1 March	8.42
D	HKWS-HKSS	LEICA GR50	LEICA GR50	1 January	6812.47
E	HKQT-HKLM	TRIMBLE NETR9	TRIMBLE NETR9	1 January	12,486.21
F	HKQT-HKSC	TRIMBLE NETR9	LEICA GR50	1 January	8183.52

, where the three ultra-short baselines are marked as A, B, and C, and the three short baselines are marked as D, E, and F.

From Table 1, we can find that baselines A, B, D, and E used the homogeneous receivers on both the base and rover stations, while baselines C and F used the heterogeneous receivers. The three ultra-short baselines were all shorter than 10 m, and the lengths of the three short baselines varied from 5 to 15 km.

3.1. IFB Calibration

According to Equations (5) and (7), the pseudorange and carrier phase observation equations are biased by IFCB and IFPB, respectively, and these biases must be removed before the subsequent AR.

For the IFCB, Shi et al. [33] analyzed their characteristics and employed both the linear and quadratic frequency functions for different receivers. Zhou et al. [34] used data from 132 International GNSS Service stations, estimated the IFCB for each satellite with a datum constraint, and compared its performance with linear or quadratic polynomial functions. To simplify the model, in this paper only the linear function was adopted, i.e., $\Delta B_{rb}^s=k^s·c\Delta H_{rb}$, where $c\Delta H_{rb}$ is the IFCB rate in meters. Thus, the IFCB in Equation (6) can be re-written as:

$$IFC{B}_{rb}^{pq}=k^q·c\Delta H_{rb}-k^p·c\Delta H_{rb}=k^{pq}·c\Delta H_{rb}.$$

(13)

In this model, the IFCB rate can be easily determined using the DD pseudorange observation expressions in Equation (5). In the following procedure, the code observables will be significantly down-weighted, and this will hardly cause noticeable bad effects on the solution [16].

For the IFPB, some researchers used the calibration strategy before the positioning process [12,13,14], and others adopted the estimation strategy in the processing procedure [15,16,17]. Since the IFPB has a great impact on the AR, it should be precisely estimated. Herein, we will propose an easy-to-implement method to quickly determine the accurate IFPB value.

As derived from the SD carrier phase observation equation in meters in Equation (3), we can obtain the DD carrier phase observation equation in cycles for the GLONASS as follows:

$$\Delta \nabla {\phi }_{rb}^{pq}=\frac{\Delta \nabla {\rho }_{rb}^q}{{\lambda }^q}-\frac{\Delta \nabla {\rho }_{rb}^p}{{\lambda }^p}-\Delta \nabla N_{rb}^{pq}+f^{pq}·\Delta d{t}_{rb}+\Delta \nabla b_{rb}^{pq}+\Delta \nabla {\epsilon }_{\phi },$$

(14)

where $f^{pq}$ is the frequency difference between two satellites. Furthermore, the frequency can be expressed as a linear function of $k^s$, i.e.:

$$f^s=f^0+a{k}^s,$$

(15)

where $f^0$ is the central frequency and a is the constant coefficient. For the L1 signal, they are 1602 MHz and 9/16 MHz, respectively. $f^{pq}$ is then transformed to $a{k}^{pq}$. Since we know that $k^s\in \left[-7, 6\right]$, if the precision of the estimated $\Delta d{t}_{rb}$ can reach 10 ns, the largest bias caused by $\Delta d{t}_{rb}$ will be smaller than 0.08 cycles and the ambiguity rounding can only be affected by the IFPB.

Furthermore, the linear relationship between the IFPB and $k^s$ is given as [11]:

$${\lambda }^s·\Delta b_{rb}^s=k^s·c\Delta h_{rb},$$

(16)

where $c\Delta h_{rb}$ is the IFPB rate in meters. Combining Equations (15) and (16), $\Delta \nabla b_{rb}^{pq}$ can be derived:

$$\Delta \nabla b_{rb}^{pq}=k^q{f}^q\Delta h_{rb}-k^p{f}^p\Delta h_{rb}=f^0{k}^{pq}\Delta h_{rb}+a({(k^q)}^2-{(k^p)}^2)\Delta h_{rb}.$$

(17)

Benefitting from the previous research by Wanninger [11], an empirical range of $c\Delta h_{rb}$ that is less than 6 cm can be assumed as the tolerance, i.e., $\Delta h_{rb}$< 0.2 ns. For the single-frequency case, when the condition $k^{pq}=\pm 1$ is satisfied, the bias $\Delta \nabla b_{rb}^{pq}$ is smaller than 0.322 cycles, and the bias caused by $\Delta d{t}_{rb}$ is further reduced to no more than 0.006 cycles. As a result, for baselines with known coordinates for the base station and rover station, the total bias in Equation (14) will be no more than 0.33 cycles. If we temporarily ignore the effect of the hardware delay, the DD ambiguity can be estimated by directly rounding:

$$\Delta \nabla N_{rb}^{pq}=round\left(\frac{\Delta \nabla {\rho }_{rb}^q}{{\lambda }^q}-\frac{\Delta \nabla {\rho }_{rb}^p}{{\lambda }^p}-\Delta \nabla {\phi }_{rb}^{pq}+k^{pq}a\Delta d{t}_{rb}\right),$$

(18)

where round(*) means the rounding off operation. In this way, the DD integer ambiguity is determined, and the next step is to calculate the IFPB rate.

It should be mentioned that the efficiency of this algorithm is affected by two factors: the precision of $\Delta d{t}_{rb}$ and the limit of the IFPB rate. If the above limit of 6 cm for IFPB rate is used, the precision of $\Delta d{t}_{rb}$ can be as tolerant as 200 ns. On the other hand, considering there is no guarantee that the above critical value of the IFPB rate is the absolute limit, the IFPB rate that can be supported by our algorithm should be no larger than 9 cm, where the precision of $\Delta d{t}_{rb}$ should not exceed 10 ns. Otherwise, the DD integer ambiguity cannot be fixed to its right value and the algorithm is inapplicable.

According to Equation (16) and the definition of IFPB in Equation (6), the IFPB is expressed in a similar model as the IFCB in Equation (13), i.e., $IFP{B}_{rb}^{pq}=k^{pq}c\Delta h_{rb}$. From the DD carrier phase observation equation in Equation (5), we can obtain the IFPB rate as follows:

$$c\Delta h_{rb}=\left(\Delta \nabla {\varphi }_{rb}^{pq}-\Delta \nabla {\rho }_{rb}^{pq}+{\lambda }^q\Delta \nabla N_{rb}^{pq}+\left({\lambda }^q-{\lambda }^p\right)\Delta N_{rb}^p\right)/k^{pq},$$

(19)

where the float DD ambiguity is converted to the combination of an integer DD ambiguity $\Delta \nabla N_{rb}^{pq}$ and the SD ambiguity for the reference satellite $\Delta N_{rb}^p$. The initial value of $\Delta {\tilde{N}}_{rb}^p$ can be easily calculated via Equation (3) when the hardware delays are ignored. The bias in the calculated $\Delta {\tilde{N}}_{rb}^p$ will be dramatically reduced due to the mm-level’s coefficient $\left({\lambda }^q-{\lambda }^p\right)$. If no cycle slip occurs, $\Delta N_{rb}^p$ remains the same and the IFPB rate can be determined.

Since there may be only a very few DD observables meeting $k^{pq}=\pm 1$ at a single epoch, the obtained IFPB rates may be very noisy. Therefore, a sufficient duration is required to smooth the noisy IFPB rates and achieve a relatively accurate initial value. After the initial IFPB rate is calculated, the bias $\Delta \nabla b_{rb}^{pq}$ in Equation (17) can thereby be calculated, and all DD ambiguity can then be re-estimated using the corresponding DD carrier phase observables:

$$\Delta \nabla N_{rb}^{pq}=round\left(\frac{\Delta \nabla {\rho }_{rb}^q}{{\lambda }^q}-\frac{\Delta \nabla {\rho }_{rb}^p}{{\lambda }^p}-\Delta \nabla {\phi }_{rb}^{pq}+k^{pq}a\Delta d{t}_{rb}+\Delta \nabla b_{rb}^{pq}\right).$$

(20)

Similarly, from Equation (19) we can get the IFPB rate again at each epoch using all DD observables, and further obtain a much more accurate value through smoothing. In other word, the IFPB rate is iterated from an initial value to a much more accurate value. The whole procedure is listed in Figure 1.

Before applying the above procedure, some remarks should be made. The first is the inputs, where only the pair of observation data files with known coordinates can be used in our algorithm. The second is the condition difference to calculate the initial and the final IFPB rate. When we use DD observables to determine the DD ambiguity the first time, the condition $k^{pq}=\pm 1$ must be satisfied. However, when we determine the DD ambiguity again after we obtained the initial IFPB rate, all DD observables were utilized to get more calculated IFPB rates and thus a more precise result.

Taking the data in Figure 2 as examples and using the above-proposed method, we obtained the epoch-by-epoch IFCB rates and IFPB rates, which are plotted in Figure 2 The final smoothed IFCB rates, IFPB rates, and corresponding standard deviations (STDs) are listed in Figure 2.

From Figure 2, we can see that the IFB rates were very stable during one day’s observation session for all six baselines; thus, we can calibrate this bias using a pre-estimated value. Meanwhile, the IFPB rates for all baselines with homogeneous receivers were almost normally distributed and the statistical results were close to zero, while the IFB rates for baselines with heterogeneous receivers were obviously biased.

Table 2. Estimated IFCB rates, IFPB rates, and the corresponding standard deviations (STDs) for different baselines.

Baseline	A	B	C	D	E	F
IFCB rate (m)	−0.0048	−0.0544	−0.3336	0.0148	−0.0297	0.1576
STD (m)	0.3318	0.2796	0.3202	0.1367	0.1930	0.1764
IFPB rate (m)	−0.0001	−0.0001	−0.0040	0.0000	−0.0002	0.0312
STD (m)	0.0018	0.0005	0.0006	0.0017	0.0030	0.0029

shows that the IFPB rate for baseline F reached up to 0.03 m. If $k^{pq}$ is larger than 3, the IFPB for this baseline will exceed half cycle, and severely obstruct the subsequent AR. Furthermore, the STD values reflected that by adopting the proposed algorithm, we could obtain a 0.3-m-level IFCB rate and a 0.003-m-level IFPB rate.

After the IFCB and IFPB rates were determined, the IFCB and IFPB can be calculated and expressed as:

$$\{\begin{array}{c}IFC{B}_{rb}^R=\Delta B_{rb}^R-\Delta B_{rb}^{R0}=k^R·c\Delta H_{rb},\\ IFP{B}_{rb}^R={\lambda }^R·\Delta b_{rb}^R-{\lambda }^{R0}·\Delta b_{rb}^{R0}=k^R·c\Delta h_{rb}.\end{array}$$

(21)

The calculated IFBs were then used for calibration in the tightly combined model.

3.2. ISB Estimation

Taking a tight combination between the GPS and GLONASS as an example, after the IFB is calibrated, the tightly combined observation equations can be re-written as:

$$\{\begin{array}{c}\Delta \nabla P_{rb}^{GR}=\Delta \nabla {\rho }_{rb}^{GR}+ISC{B}_{rb}^{GR}+\Delta \nabla {\epsilon }_P,\\ \Delta \nabla {\varphi }_{rb}^{GR}=\Delta \nabla {\rho }_{rb}^{GR}-\left({\lambda }^R·\Delta N_{rb}^R-{\lambda }^G·\Delta N_{rb}^G\right)+ISP{B}_{rb}^{GR}+\Delta \nabla {\epsilon }_{\varphi }.\end{array}$$

(22)

From Equation (9), and considering the consistency of the hardware delays for all satellite pairs, we can regard both the ISCB and ISPB as short-term constants, which means their values remain unchanged for at least several hours of continuous operation [23,32]. The ISCB can be easily estimated using a least-squares adjustment or Kalman filtering. As for the ISPB, from Equation (22), we find the DD ambiguity has lost its integer nature. Therefore, a transformation similar to Equation (19) is employed to recover the integer nature, and the tightly combined DD carrier phase observation equation is re-expressed as follows:

$$\Delta \nabla {\varphi }_{rb}^{GR}=\Delta \nabla {\rho }_{rb}^{GR}-{\lambda }^R\Delta \nabla N_{rb}^{GR}-\left({\lambda }^R-{\lambda }^G\right)\Delta N_{rb}^G+ISP{B}_{rb}^{GR}+\Delta \nabla {\epsilon }_{\varphi }.$$

(23)

Similarly, the float DD ambiguity is separated into the terms of $\Delta \nabla N_{rb}^{GR}$ and $\Delta N_{rb}^G$, and the initial value of $\Delta {\tilde{N}}_{rb}^G$ can be calculated. If no cycle slip occurs, the reduced bias in the calculated $\Delta {\tilde{N}}_{rb}^G$ will be absorbed into the ISPB.

Since the ISPB and the DD ambiguity terms are also very difficult to separate, they will be regarded as one parameter and estimated together with the coordinates, or directly calculated using the following formula:

$$\Delta N_{rb}^{GR}-ISP{B}_{rb}^{GR}/{\lambda }^R=\left(\Delta \nabla {\rho }_{rb}^{GR}-\Delta \nabla {\varphi }_{rb}^{GR}-\left({\lambda }^R-{\lambda }^G\right)\Delta {\tilde{N}}_{rb}^G\right)/{\lambda }^R.$$

(24)

Herein, the geometry can be obtained when we estimate the ISCB. After rounding off, the integer part of the real ISPB is absorbed into the integer DD ambiguity, and the residual is changed to a constant value by multiplying the coefficient ${\lambda }^R/{\lambda }^{R0}$ and then taken as the ISPB in cycles. Furthermore, if the determined ISPB is around ±0.5 cycles, one cycle is added to the ISPB around −0.5 cycles to make all results consistent.

For a tight combination between GPS and BDS, there may be a 0.5-cycles inter-satellite type bias in BDS between the geostationary orbit (GEO) satellites and the non-GEO satellites [35], and this bias is corrected in advance. For a tight combination between GPS and Galileo, the overlapping frequency is used and the SD ambiguity term is eliminated, thus the model becomes simple.

The determined ISPBs for each pair of satellites are averaged together to one ISPB value at one epoch. For a short-term observation period, the ISPB can be further smoothed due to its stability to achieve its final value.

Using data in Table 1 and the above method, we obtained the epoch-by-epoch ISCB and ISPB, which are plotted in Figure 3. The final smoothed ISCBs, ISPBs, and corresponding STD values are listed in

Table 3. Estimated ISCB, ISPB, and the corresponding STDs between GLONASS/BDS/Galileo and GPS. The symbols “G-R”, “G-C”, and “G-E” in the table mean operation between GPS and GLONASS, between GPS and BDS, and between GPS and Galileo, respectively.

Baseline	ISCB (m)			ISPB (cycle)
	G-R	G-C	G-E	G-R	G-C	G-E
A	−0.1624	−1.2582	−0.0323	−0.0515	0.4706	−0.0013
STD-A	0.5022	0.3825	0.3503	0.0378	0.0201	0.0087
B	0.7834	0.2206	−0.0404	−0.1605	−0.0095	−0.0012
STD-B	0.6134	0.3476	0.5469	0.0233	0.0132	0.0143
C	0.6460	3.4453	0.4662	0.3141	0.5612	−0.0007
STD-C	0.7422	0.3678	0.4579	0.0414	0.0158	0.0143
D	0.2339	0.0017	−0.0181	−0.1089	0.4182	−0.0037
STD-D	0.3655	0.1964	0.2480	0.0266	0.0169	0.0216
E	−1.1997	0.4515	0.3265	−0.0354	0.4248	−0.0163
STD-E	0.5001	0.1702	0.1736	0.0365	0.0395	0.0347
F	−0.9774	−1.5025	−4.1108	0.1078	0.1831	0.5008
STD-F	0.4391	0.2043	0.2138	0.0483	0.0331	0.0369

.

From Figure 3 and Table 3, we can see that the ISB had a diurnal stability and could be corrected in advance. In addition, the ISCB and ISPB differ from one another, whether homogeneous or heterogeneous receivers are used. The ISCB could be as large as 4 m and the ISPB could reach up to 0.5 cycles. Such large biases will no doubt hamper the subsequent AR and must be removed. In the meantime, the given STD values show that the precision of the estimated ISCB was almost at the same level as the code measurement noise, and the ISCB between GPS and GLONASS was slightly larger than other inter-system pairs. This was reasonable due to the large code noise of GLONASS. As for ISPB, for ultra-short baselines, the precision could be up to 0.015 cycles, and a general precision of 0.03 cycles was achieved.

What is interesting is that for most baselines, the ISPB values between GPS and Galileo were close to zero because of their identical frequency, with the only exception being for baseline F. This indicates that the ISPB may exist even if overlapping frequencies are used in tightly combined systems.

4. Numerical Analysis and Discussion

After the removal of the IFB and ISB, the DD ambiguity can be resolved in both loosely and tightly combined models. In this part, we will conduct some experiments and set three different scenarios to investigate and compare the performance of TC with that of LC using multi-GNSS single-frequency data.

Three different scenarios were set according to the satellite elevation mask and azimuth range to simulate different observation environments, which are shown in

Table 4. Three scenarios with different elevation and azimuth masks.

Scenario	Case 1	Case 2	Case 3
Elevation mask (°)	10	40	10
Azimuth range (°)	0–360	0–360	180–360

. Case 1 was the typical situation in an open area, case 2 simulated the ambient occlusion with surrounding high buildings, and case 3 simulated the unilateral shade under a block of high buildings.

In the following discussion, we first give the epoch-by-epoch results of baseline A in case 3 as an example. A final comparison of the float solution and fixed solution between TC and LC in all three scenarios will be made subsequently.

Figure 4 gives the usable number of satellites of each system and the position dilution of precision (PDOP) for all four systems. From the figure, we can see that because of signal obstructions, for most of the time, fewer than five satellites can be used for each system, which means efficient precise positioning is difficult to realize if only a single system is used. However, when all satellites are combined, this problem will be solved. This can be inferred from the corresponding PDOP, whose average value was only 2.3.

4.1. Float Solution

In the data processing, we assumed all code precision was 0.3 m, except for GLONASS with 0.45 m, and all carrier phase precision was 0.003 m. Setting the coordinate and ambiguity as unknown parameters and ignoring the integer nature of the ambiguity, the single-epoch float solution was first obtained using a standard least-squares adjustment ahead of the ambiguity fix. The horizontal and vertical positioning biases for the LC model and TC model, along with the corresponding 95% confidence intervals and standard deviations, are shown in Figure 5. The results indicate that the float solution for TC model in both the horizontal and vertical directions had a moderate improvement compared with those for the LC model.

The position root-mean-square errors (RMSEs) in the TC model and LC model using all six baselines under three scenarios and the corresponding improvements were calculated and are listed in

Table 5. Comparison of RMSEs in the float solution between LC and TC under different scenarios and the corresponding improvements.

Tag	Case 1			Case 2			Case 3
	LC (m)	TC (m)	Pi (%)	LC (m)	TC (m)	Pi (%)	LC (m)	TC (m)	Pi (%)
A	0.898	0.826	8.02	2.487	1.874	24.65	2.317	1.983	14.42
B	1.049	1.018	2.96	2.906	2.533	12.84	2.872	2.323	19.12
C	0.960	0.895	6.77	5.275	2.171	58.84	2.797	2.268	18.91
D	0.444	0.429	3.38	1.260	1.110	11.90	1.125	1.044	7.200
E	0.501	0.474	5.39	1.228	0.975	20.60	1.295	1.095	15.44
F	0.587	0.556	5.28	1.619	1.278	21.06	2.080	1.709	17.84
Mean	0.740	0.700	5.43	2.463	1.657	32.72	2.081	1.737	16.53

. The results show that in typical situations (case 1), there were few improvements in the positioning errors in the TC model compared with those in the LC model, but in other two situations (cases 2 and 3) with obstructed satellites, obvious decreases in positioning errors were found. The improving percentage is represented as P_i, and the average values for all six baselines exceeded 30% and 15% for cases 2 and 3, respectively.

The results indicate that the performance of the TC model was similar to that of the LC model in open areas, but was superior to the LC model when there were signal obstructions around. In these complex environments, the TC model was a better choice. In the meantime, we also found that cases 2 and 3 had different degrees of influence on the performance of the TC and LC models, where the TC model could reduce the negative impacts more in situations with a high elevation mask.

4.2. Fixed Solution

For single-epoch precise positioning, the unknown integer ambiguity must be determined to acquire the centimeter level’s fixed solution. After the float ambiguity and its variance were obtained in the float solution, the least-squares ambiguity decorrelation adjustment method [36] was employed to search for the integer ambiguity. During the ambiguity validation procedure, both the data-driven R-ratio test and the model-driven success rate [37] were considered, whose empirical thresholds were set as 3.0 and 99%, respectively, in this paper. When both these two values reached or exceeded their thresholds, the ambiguity was regarded as being fixed. Meanwhile, we realized that the ambiguity may be incorrectly fixed. Among the epochs with fixed ambiguities, only when the positioning biases in all three directions were smaller than 0.1 m could the ambiguity be regarded as being correctly fixed.

Based on the above consideration, we defined the ambiguity fixing rate $P_{fix}$ and ambiguity correctly-fixed rate $P_c$ as:

$$\begin{array}{cc}{P}_{fix}=\frac{N_{fix}}{N_{All}}\times 100\%,& P_c=\frac{N_c}{N_{fix}}\times 100\%\end{array},$$

(25)

where $N_{All}$ is the total epoch number, $N_{fix}$ is the epoch number of the fixed ambiguity, and $N_c$ is the epoch number of the correctly fixed ambiguity.

The test results with the reciprocal values of the R-ratio are plotted in Figure 6. It is clear that more ratio values exceeded the threshold in the TC model than those in the LC model. In other words, in the TC model, more fixed epochs were achieved.

Combining the ratio-test and the AR success rate, we could calculate $P_{fix}$ for both the LC and TC models under different situations; the results for all six baselines are listed in

Table 6. Comparison of the ambiguity fixing rates between LC and TC under different scenarios and the corresponding differences.

Tag	Case 1			Case 2			Case 3
	LC (%)	TC (%)	${\Delta }_{\mathit{f}\mathit{i}\mathit{x}}(\%)$	LC (%)	TC (%)	${\Delta }_{\mathit{f}\mathit{i}\mathit{x}}(\%)$	LC (%)	TC (%)	${\Delta }_{\mathit{f}\mathit{i}\mathit{x}}(\%)$
A	100.00	100.00	0.00	71.08	99.55	28.47	41.91	74.65	32.74
B	99.86	99.86	0.00	48.96	78.92	29.97	21.91	58.65	36.74
C	99.97	99.97	0.00	61.49	83.40	21.91	23.23	69.17	45.94
D	99.58	99.58	0.00	66.35	99.27	32.92	54.34	95.73	41.39
E	74.58	72.26	−2.33	67.85	95.73	27.88	35.66	68.40	32.74
F	89.51	86.04	−3.47	65.28	97.85	32.57	38.47	74.20	35.73
Mean	93.92	92.95	−0.97	63.50	92.45	28.95	35.92	73.47	37.55

. The difference in $P_{fix}$ between TC and LC is denoted as ${\Delta }_{fix}$, i.e., ${\Delta }_{fix}=P_{fix}\left(TC\right)-P_{fix}\left(LC\right)$. It shows that $P_{fix}$ in the LC and TC models were nearly the same in case 1, which means they had a similar performance in normal situations. However, in cases 2 and 3, the TC model outperformed the LC model, where the difference of $P_{fix}$ could be more than 30%. It is clear that the TC model increased the performance of ambiguity fixing.

Furthermore, by judging whether the ambiguity was correctly fixed, we could calculate $P_c$ for both LC and TC models under different situations, whose results for all six baselines are listed in

Table 7. Comparison of correctly fixed ambiguity rates between LC and TC under different scenarios and the corresponding differences.

Tag	Case 1			Case 2			Case 3
	LC (%)	TC (%)	Δc (%)	LC (%)	TC (%)	Δc (%)	LC (%)	TC (%)	Δc (%)
A	96.70	96.15	−0.56	82.85	93.41	10.55	73.74	85.44	11.71
B	97.64	98.12	0.49	82.27	86.05	3.78	68.78	79.04	10.26
C	97.78	96.56	−1.22	82.50	90.17	7.68	61.14	79.12	17.98
D	93.93	94.11	0.17	91.99	95.49	3.49	79.81	90.46	10.65
E	78.31	80.54	2.23	84.24	90.71	6.48	61.73	75.63	13.90
F	72.19	68.72	−3.46	79.84	84.32	4.47	42.33	54.28	11.95
Mean	89.42	89.03	−0.39	83.95	90.03	6.08	64.59	77.33	12.74

. The difference in $P_c$ between TC and LC is written as Δ_c. From the table, we can see that just like the results of $P_{fix}$, in case 1, we obtained similar correctly fixed ambiguity rates for TC and LC. In cases 2 and 3, the correctly fixed ambiguity rates in the TC model showed a moderate improvement. Considering the remarkable improvement of the ambiguity fixing rate, many more ambiguities were reliably fixed to their correct values. The TC model demonstrated significant advantages in AR in special situations with signal occlusion.

It should be pointed out that in cases 1, 2, and 3 of baseline F, there were six, five, and six epochs when the LC model was unavailable because there was only one usable satellite in the GLONASS and Galileo systems. However, this issue was resolved in the TC model. The only exception was for three epochs in case 3 when there was only one satellite in each system, thus the formed three-DD observation equations could not be used for positioning.

The fixed ambiguity could be taken as a known parameter and the fixed solution could thereby be calculated. Excluding the results with an incorrectly fixed ambiguity, we display the positioning biases in the horizontal and vertical directions for baseline A in case 3 with both LC and TC models in Figure 7. The 95% confidence intervals and standard deviations are also plotted in the figures.

We could hardly find any significant differences between the TC and LC models. Both achieved nearly 1-cm-level positioning results. The results for other baselines and other situations were almost the same and are not given in this paper. Therefore, from the perspective of positioning results, the TC and LC models showed the same level of performance.

4.3. Discussion

The above experimental results indicate an obvious advantage of the TC model in situations with obstructed satellites. Although the TC and LC models gave comparable results in open sky, this was usually not the case for urban residents. The pedestrian on the streets must face a situation under high buildings, where severe signal occlusion inevitably happens and hampers efficient navigation. From the functional models of TC and LC, if the ISBs are calibrated ahead of the positioning, the TC model can provide one more observation equation for each pair of systems, and the degree of freedom becomes large [30]. When the four systems are used together, the redundancy increases to three. In situations with limited tracked satellites, the additional equations can make a positive impact on the AR and positioning performance. That is why better results were achieved in TC model, as shown in Table 5 and Table 6.

Meanwhile, the functional model of TC infers that many more advantages can be acquired if more systems are involved compared with LC. If the regional satellite systems, such as Japanese QZSS and Indian IRNSS, are not taken into account, the integration of GPS, GLONASS, Galileo, and BDS maximizes the performance of a tightly combined solution. Although lots of contributions have been made using the hybrid systems with CDMA signals, it is a pity that in their work, the GLONASS signals are abandoned, maybe due to its FDMA signals and the corresponding IFB. This paper provides a further investigation of precise positioning with all four global navigation satellite systems by considering both IFB and ISB.

Furthermore, as more and more satellites are available, we do not need to resolve all tracked satellites since fixing the full set of ambiguities of multiple constellations may possibly decrease the AR success rate and computational efficiency, and the partial ambiguity resolution (PAR) technique, which resolves a subset of the ambiguities, was suggested to achieve a high success rate [38]. After the selection of an ambiguity set based on a certain strategy, the number of remaining satellites may be not that sufficient. At this moment, the TC model can be applied to ensure the position estimation. Therefore, the combination of the PAR technique and the tightly combined multi-GNSS system will bring about a high AR success rate and ambiguity fixing rate, along with time-saving data processing, regardless of whether the observation environment is in clear-sky or suffers from signal occlusion. This is expected to become a more efficient work pattern in the coming multi-GNSS era.

5. Conclusions

In this paper, we provided a rigorous functional model for tightly combined single-frequency multi-GNSS data. The virtual signal of the constant central frequency in the GLONASS system with a channel number of zero was introduced, and the IFPB and ISPB in GLONASS involving a tightly combined system were re-defined from the hardware delays.

Based on this functional model, we proposed a practical algorithm to estimate the IFB rate and the ISB, and their characteristics were analyzed using six pairs of available single-frequency multi-GNSS data of ultra-short and short baselines. Numerical results confirmed that the IFPB rates for baselines with homogeneous receivers were negligible, while those with heterogeneous receivers should be calibrated. Results also indicated that there may be considerable ISPB even if overlapping frequencies are used in tightly combined systems.

After the removal of the IFB and ISB, a comprehensive performance comparison was made from the float solution to the fixed solution between the TC and LC models. Experiments were conducted to simulate three different observation environments: one is the typical situation and the other two suffer from signal obstruction with high elevation angles or limited azimuths, respectively. From the results, we found that in typical situations, the TC model and LC model had a similar performance in the float solution and AR. However, when the satellite signals were severely obstructed and few satellites could be tracked in a single system, the float solution and ambiguity fixing rates in the LC model were dramatically decreased, while those in the TC model had relatively smaller declines, where the difference in ambiguity fixing rates could be as large as 30%. The correctly fixed ambiguity rates in the TC model also had a moderate improvement of around 10%. Once the ambiguity was fixed, both models achieved a similar positioning accuracy. Summarily, the TC model significantly outperformed the LC model in terms of AR in cases of signal obstruction, though in other cases, they had a similar performance.