Differential Inter-System Biases Estimation and Initial Assessment of Instantaneous Tightly Combined RTK with BDS-3, GPS, and Galileo

Abstract

In this contribution, we assess, for the first time, the tightly combined real-time kinematic (RTK) with GPS, Galileo, and BDS-3 operational satellites using observations from their overlapping L1-E1-B1C/L5-E5a-B2a frequencies. First, the characteristics of B1C/B2a signals from BDS-3 operational satellites is evaluated compared to GPS/Galileo L1-E1/L5-E5a signals in terms of observed carrier-to-noise density ratio, pseudorange multipath and noise, as well as double-differenced carrier phase and code residuals using data collected with scientific geodetic iGMAS and commercial M300Pro receivers. It’s demonstrated that the observational quality of B1C/B2a signals from BDS-3 operational satellites is comparable to that of GPS/Galileo L1-E1/L5-E5a signals. Then, we investigate the size and stability of phase and code differential inter-system bias (ISB) between BDS-3/GPS/Galileo B1C-L1-E1/B2a-L5-E5a signals using short baseline data collected with both identical and different receiver types. It is verified that the BDS-3/GPS/Galileo ISBs are indeed close to zero when identical type of receivers are used at both ends of a baseline. Moreover, they are generally present and stable in the time domain for baselines with different receiver types, which can be easily calibrated and corrected in advance. Finally, we present initial assessment of single-epoch tightly combined BDS-3/GPS/Galileo RTK with single-frequency and dual-frequency observations using a formal and empirical analysis, consisting of ambiguity dilution of precision (ADOP), ratio values, the empirical ambiguity resolution success rate, and the positioning accuracy. Experimental results demonstrate that the tightly combined model can deliver much lower ADOP and higher ratio values with respect to the classical loosely combined model whether for GPS/BDS-3 or GPS/Galileo/BDS-3 solutions. The positioning accuracy and the empirical ambiguity resolution success rate are remarkably improved as well, which could reach up to approximately 10%∼60% under poor observational conditions.

1. Introduction

China's BeiDou Navigation Satellite System (BDS) has evolved from the first phase of a demonstration system (BDS-1) to the second phase of a regional navigation satellite system (BDS-2) and is currently in construction of the third phase of a global navigation system (BDS-3) [1]. Before formal deployment of BDS-3 constellation, a demonstration system with five BDS-3 experimental satellites, including two Inclined Geosynchronous Orbit (IGSO) satellites (i.e., I1-S and I2-S) and three Medium Earth Orbit (MEO) satellites (i.e., M1-S, M2-S, and M3-S), have been established from March 2015 to February 2016 to carry out in-orbit test and verification for the new payloads, new designed signals, and new techniques [1,2,3]. The first two operational satellites for BDS-3 were successfully launched on 5 November 2017, which marked the official start of the establishment of BDS-3 [2]. By the end of 2018, China has successfully launched a total of 19 BDS-3 operational satellites, including 18 MEO satellites in operation and one Geostationary Earth Orbit (GEO) satellite under in-orbit test, to complete a primary system for providing global services [4]. Currently (early May 2019), 20 operational BDS-3 satellites have been launched (as listed in

Table 1. The status of launched BDS-3 operational satellites (early May 2019) [5].

SVN	RRN	Common Name	COSPAR ID	NORAD ID	Launch Time	Orbit
C201	C19 1	BDS-3 MEO-1	2017-069A	43001	5.11.2017	MEO
C202	C20	BDS-3 MEO-2	2017-069B	43002	5.11.2017	MEO
C203	C27	BDS-3 MEO-7	2018-003A	43107	11.01.2018	MEO
C204	C28 2	BDS-3 MEO-8	2018-003B	43108	11.01.2018	MEO
C205	C22	BDS-3 MEO-4	2018-018A	43207	12.02.2018	MEO
C206	C21	BDS-3 MEO-3	2018-018B	43208	12.02.2018	MEO
C207	C29	BDS-3 MEO-9	2018-029A	43245	29.03.2018	MEO
C208	C30	BDS-3 MEO-10	2018-029B	43246	29.03,2018	MEO
C209	C23	BDS-3 MEO-5	2018-062A	43581	29.07.2018	MEO
C210	C24	BDS-3 MEO-6	2018-062B	43582	29.07.2018	MEO
C211	C26	BDS-3 MEO-11	2018-067A	43602	24.08.2018	MEO
C212	C25	BDS-3 MEO-12	2018-067B	43603	24.08.2018	MEO
C213	C32	BDS-3 MEO-13	2018-072A	43622	19.09.2018	MEO
C214	C33	BDS-3 MEO-14	2018-072B	43623	19.09.2018	MEO
C215	C35	BDS-3 MEO-16	2018-078A	43647	15.10.2018	MEO
C216	C34	BDS-3 MEO-15	2018-078B	43648	15.10.2018	MEO
C217	C59	BDS-3 GEO-1	2018-085A	43683	1.11.2018	GEO
C218	C36	BDS-3 MEO-17	2018-093A	43706	18.11.2018	MEO
C219	C37	BDS-3 MEO-18	2018-093B	43707	18.11.2018	MEO
C220	C38	BDS-3 IGSO-1	2019-023A	44204	20.04.2019	IGSO

^1,2 PRN C19 used until 12.06.2018 and since 15.11.2018. PRN C28 used until 11.06.2018 and since 20.12.2018. As no official space vehicle numbers (SVNs) are available, SVNs are assigned based on the satellite generation (C2?? with monotonically increasing number for BDS-3). The SVN used here follows NORAD ID. In case of dual launches, the smaller NORAD ID is assigned the smaller SVN and vice versa. Information on BDS-3 is preliminary and might be subject to change.

) [5]. BDS-3 is expected to complete its full constellation with 30 satellites, including 3 GEO, 3 IGSO, and 24 MEO, to provide positioning and navigation services for global users by 2020 [1,4].

In order to improve the interoperability with other satellite navigation systems, in addition to the backward BDS-2 compatible B1I (1561.098 MHz) and B3I (1268.520 MHz) signals, the BDS-3 satellites (both experimental and operational) are also capable of broadcasting new B1C (1575.42 MHz), B2a (1176.45 MHz), and B2b (1207.14 MHz) open service navigation signals (as listed in

Table 2. Available open service signals transmitted by BDS-3 operational satellites [1,4,6].

Band	Frequency (MHZ)	Signal Component	Modulation	Code Rate (Mcps)
B1	1575.42	B1CD	1 BOC (1,1)	1.023
		B1CP	2 QMBOC (6,1,4/33)	1.023
	1561.098	B1I	3 BPSK (2)	2.046
B2	1176.45	B2aD	4 QPSK (10)	10.23
		B2aP		10.23
	1207.14	B2bI	QPSK (10)	10.23
		B2bQ		10.23
B3	1268.52	B3I	BPSK (10)	10.23

^1,2,3,4 Abbreviations: BOC, binary offset carrier; QMBOC, quadrature multiplexed binary offset carrier; BPSK, binary phased shift keying; QPSK, quadrature phased shift keying. Subscript D represents data component and P represents pilot component.

) that overlap with the GPS L1/L5 and Galileo E1/E5a/E5b signals [1,4,6]. During the past two years, the characteristics of these new signals from BDS-3 experimental satellites has been intensively investigated in various aspects, including signal carrier-to-noise density ratio (C/N₀), pseudorange multipath effects, triple-frequency carrier phase ionosphere-free and geometry-free combination, together with precision of code and carrier phase observations, etc. [7,8,9,10,11]. Their results demonstrated that the observational quality of the new-generation B1C/B2a/B2b signals from BDS-3 experimental satellites is comparable to their GPS/Galileo counterparts. Moreover, the performance of precise orbit determination, satellite clocks, precise point positioning, and relative positioning were also evaluated using observations from BDS-3 experimental satellites or together with BDS-2 satellites [2,7,8,12,13,14,15]. Recently, several researchers have further presented their initial evaluation results for the new navigation signals transmitted by eight BDS-3 operational satellites, and similar results are also achieved as the experimental satellites [2,3,16].

The availability of overlapping frequencies between BDS-3 and GPS/Galileo brings opportunity for the tightly combined processing or inter-system double-differencing of observations from GPS/Galileo, and BDS-3 operational satellites. Different from the classical loosely combined double-differencing of observations from each system independently (i.e., choosing a single-reference satellite for each system separately), the tightly combined observational model utilizes a single-reference satellite for observations from all the involved Global Navigation Satellite Systems (GNSSs). As a result, additional double-differenced observations can be created, and the model strength is then improved [17,18]. However, when inter-system observations are created, additional between-receivers differential inter-system biases (ISBs), which originate from difference in the receiver hardware delays that affecting signals from different GNSS systems, are also introduced and should be carefully considered [17]. The phase and code differential ISBs are frequency dependent, which means that they may present with different values for different frequencies [17,18,19]. Previous studies have also revealed that the phase and code ISBs between overlapping frequencies that are shared by GPS, Galileo, QZSS, and IRNSS (i.e., L1-E1 and L5-E5a) are receiver-type dependent. They are close to zero when identical types of receivers are utilized at both ends of a baseline, but generally present with non-zero values for baselines with different receiver types. Fortunately, since the ISBs are stable in the time domain, they can be easily calibrated and corrected in advance [17,18,19,20]. If these ISBs are a-priori calibrated and introduced as corrections, the performance of precise relative positioning can be remarkably improved with respect to the loosely combined model in terms of ambiguity resolution and positioning accuracy [21,22,23,24,25].

Currently, most studies only focus on tightly combined relative positioning using observations from the overlapping frequencies between GPS, Galileo, QZSS, and IRNSS (i.e., L1-E1 and L5-E5a). Regarding BDS, existing contributions either focus on the overlapping frequency of BDS-2 B2I and Galileo E5b signals [19,24] or the non-overlapping frequencies between BDS-2 and GPS (i.e., B1I-L1 and B2I-L2) [26,27,28,29,30]. None have focused on the model and performance evaluation of tightly combined precise relative positioning of BDS-3 (especially operational satellites) with other GNSS systems such as GPS, Galileo, etc. As an aspect of evaluating the interoperability of BDS-3 operational satellites and their new navigation signals with other systems, for the first time this contribution will investigate the feasibility of conducting tightly combined real-time kinematic (RTK) positioning between BDS-3 B1C/B2a signals and GPS/Galileo L1-E1/L5-E5a signals, and how well it performs. The results could provide us with new insights into the interoperability of China’s BDS-3 with modernized GPS and Galileo.

In this contribution, we will provide initial assessment of tightly combined RTK with GPS, Galileo, and BDS-3 operational satellites using observations from their overlapping frequencies of B1C-L1-E1/B2a-L5-E5a. We first evaluate the characteristics of B1C/B2a signals from BDS-3 operational satellites in comparison to GPS/Galileo L1-E1/L5-E5a signals in terms of observed carrier-to-noise density ratio, pseudorange multipath and noise, together with double-differenced carrier phase and code residuals using data collected with scientific geodetic iGMAS and commercial M300Pro receivers. Then, we investigate the size and stability of phase and code differential ISBs between BDS-3/GPS/Galileo B1C-L1-E1/B2a-L5-E5a signals using short baseline data collected with both identical and different receiver types. Finally, the performance of single-epoch tightly combined BDS-3/GPS/Galileo RTK with single-frequency (B1C-L1-E1) and dual-frequency (B1C-L1-E1/B2a-L5-E5a) observations is evaluated by a formal as well as an empirical analysis, consisting of ambiguity dilution of precision (ADOP), ratio values, empirical ambiguity resolution success rate and the positioning accuracy. The emphasis is on instantaneous positioning since this is obviously the most challenging for integer ambiguity resolution. This approach has the advantage that the results will be insensitive to the occurrences of troublesome carrier phase cycle slips.

The remainder of this contribution is organized as follows. Section 2 overviews the GNSS data collection. Section 3 evaluates the characteristics of BDS-3 B1C/B2a signals in comparison with GPS L1/L5 and Galileo E1/E5a signals. In Section 4, we first briefly introduce the observation model; and then investigate the size and stability of phase and code differential ISBs between BDS-3/GPS/Galileo B1C-L1-E1/B2a-L5-E5a signals. After that, we present performance assessment of single-epoch tightly combined BDS-3/GPS/Galileo RTK with single-frequency and dual-frequency observations. Finally, some conclusions are given in Section 5.

2. Data Collection

Compared with our previous research [7], the scientific geodetic iGMAS receivers (manufactured by one university in China) [7], commercial ComNav M300Pro (manufactured by ComNav Technology Ltd. in Shanghai, China) receivers, together with two Trimble Zephyr Geodetic 2 antennas are used in the experiments here. All the selected GNSS receivers can simultaneously track GPS L1/L2/L5, Galileo E1/E5a/E5b, BDS-2 B1I/B2I/B3I, and BDS-3 B1I/B3I/B1C/B2a/B2b signals, except that B2b cannot be tracked by M300Pro. The tracking modes of GPS/Galileo/BDS-3 L1-E1-B1C/L5-E5a-B2a signals, which are our major concerns in the following analysis, are listed in

Table 3. The observables including tracking mode for the iGMAS and M300Pro receiver.

Receiver	GPS		Galileo		BDS-3
iGMAS	L1 C/A	L5I	E1B	E5aI	B1C Data	B2a Data
M300Pro	L1 C/A	L5I	E1 (B + C)	E5a (I + Q)	B1CC 1	B2aA 2

^1,2 They do not conform to RINEX 3.04 and currently their tracking modes are not visible to the users.

.

The observations are also collected with 10 s sampling interval on the roof of the School of Geodesy and Geomatics (SGG) building at Wuhan University but in June 2018, which is different from our previous experiments in [7]. Figure 1 shows the observational environment and the GNSS antennas used in the experiments, while Figure 2 illustrates the sky plot and availability of BDS-3, Galileo, and GPS satellites with 10° elevation cut-off angle from 3–4 June 2018. 12 GPS Block IIF, 16 Galileo, and 8 BDS-3 operational satellites are tracked during the experiments.

3. Quality Assessment of the B1C/B2a Signals from BDS-3 Operational Satellites

Based on the raw data collected from 3–4 June 2018 with iGMAS and M300Pro receivers, this section evaluates the quality of the BDS-3 B1C/B2a signals by comparing their signal carrier-to-noise density ratio (C/N₀), pseudorange multipath and noise, as well as the double-differenced carrier phase and code residuals to those of the GPS/Galileo L1-E1/L5-E5a signals. The GPS satellites here only restrict to the Block IIF satellites that can simultaneously transmit L1/L5 signals, while the Galileo signals come from both In-Orbit-Validation (IOV) and Full-Operation-Capability (FOC) satellites, which is the same as our previous research [7]. However, since the signals transmitted by the Galileo IOV satellites show a lower level (approximately 5 dB-HZ) of C/N₀ in comparison with their FOC counterparts [31], and the C/N₀ of first two FOC satellites (E14 and E18) that were injected into incorrect orbital planes could be higher than regular Galileo FOC satellites [32], only regular FOC satellites are involved in the computation of average C/N₀ for Galileo.

3.1. Carrier-to-Noise Density Ratio (C/N₀)

The average C/N₀ values against elevation is used here to investigate the characteristics of the BDS-3 B1C/B2a signals. For definition of the C/N₀ and the procedures of how to compute the average C/N₀ values against elevation for each available signal, please refer to our previous work for detail [7]. The C/N₀ values against elevation for BDS-3, GPS, and Galileo with the iGMAS and M300Pro receivers are plotted in Figure 3. In general, the C/N₀ values increase with the increasing of satellite elevation angles for both receivers, from about 30 dB-Hz at low elevation angles to approximately 50 dB-Hz close to the zenith. Moreover, for the iGMAS receiver, it is obvious that the B1C signal exhibits approximately 2–8 dB-Hz lower C/N₀ values than other signals over all the elevation angles. This is reasonable considering that the data component is tracked for B1C signal [7]. With respect to the M300Pro receiver, however, insignificant differences are observed for the obtained C/N₀ values between B1C and other signals, which may due to the reason that the pilot component is tracked for B1C signal with M300Pro.

3.2. Pseudorange Multipath and Noise

The same as reference [7], we adopt the multipath (MP) combination to determine the pseudorange multipath and noise. For detailed definition of the MP combination, please refer to reference [7,33]. Figure 4 illustrates the MP combination against elevation for BDS-3 C20 (B1C/B2a/B1I/B3I signals) and BDS-2 C11 (B1I/B2I/B3I signals) with the iGMAS (top panels) and M300Pro (bottom panels) receiver. In line with our previous results for BDS-3 experimental satellites [7], the MP combination against elevation series from both iGMAS and M300Pro receivers show that the satellite-induced elevation-dependent code biases, which have been previously identified to exist in BDS-2 code observations [34], seem to be no more obviously present for BDS-3 operational satellites as well. Similar results are also observed in [3].

Similarly, the Root Mean Square (RMS) values of the MP combinations are also calculated with elevation angle interval of 5°, so as to further evaluate the quality of code observations among the various signals, and the results are shown in Figure 5. Consistent with our previous results [7], we observe that the RMS of MP combinations are clearly dependent on the satellite elevation angles for all signals. Using the iGMAS receiver, the RMS values for L1/E1/B1C are approximately 0.3 m close to the zenith and increase to approximately 1.0 m at low elevation angles. Nevertheless, the values are about 0.1 m close to the zenith and increase to approximately 0.5 m at low elevation angles for L5/E5a/B2a. Therefore, we conclude that the RMS values for L1/E1/B1C are obviously larger than L5/E5a/B2a. Similar phenomenon is observed for M300Pro receiver, but with larger RMS values. Note that the RMS values for L1 is nearly identical with L5/E5a/B2a and much smaller than B1C/E1, which is quite different from using the iGMAS receiver. We speculate that this may due to the specific hardware design of the M300Pro.

3.3. Zero-Baseline Double-Differenced Carrier Phase and Code Residuals

It has already been proved that the zero-baseline double-differenced residuals could be applied to evaluate the noise characteristics of the GNSS observations [7,35]. In this experiment, we also form a zero-baseline with two iGMAS receivers and calculate its double-differenced carrier phase and code residuals. The corresponding results for representative satellite pairs of C19-C22, G08-G01, and E09-E24 are shown in Figure 6. For the definition of zero-baseline double-differenced carrier phase and code observations together with their reduction, please refer to reference [7] for details. We observe that overall, the noise of phase and code observations from the B1C/B2a signals of BDS-3 operational satellites is comparable to GPS L1/L5 and Galileo E1/E5a signals. Most double-differenced carrier phase residuals for all available signals vary within a range of −0.005 to 0.005 cycles (−1 to 1 mm), while −0.1 to 0.1 m for the code residuals. Similar results are also obtained for M300Pro receivers, which are not shown here due to the limited space.

4. Differential ISB Estimation and Performance Evaluation of Tightly Combined RTK with BDS-3, GPS, and Galileo

In this section, we present observation models, differential ISB estimation results, and performance evaluation of tightly combined RTK with GPS, Galileo, and BDS-3 operational satellites. All GPS (including block IIF, IIA, IIR, and IIR-M) and Galileo satellites available are involved in the data analysis. However, Galileo E14 and E18 are not used in tightly combined RTK since they are not tracked by the M300Pro receiver.

4.1. Observation Models for Combined BDS-3, GPS, and Galileo RTK

In this section, we will briefly introduce the observation models corresponding to loosely and tightly combined RTK with BDS-3 B1C/B2a, GPS L1/L5, and Galileo E1/E5a observations. Since our research is restricted to zero and short baselines, the atmospheric delays are ignored.

Assuming receivers b and r simultaneously track satellites s_G = 1_G, 2_G...m_G and s_* = 1_*, 2_*...m_* on overlapping frequency $i=1,2...f$, where $s_G$ and $s_{*}$ are the number of satellites from GPS and constellation * (The symbol * is B for BDS-3 and E for Galileo) respectively, the undifferenced phase and code observation equation for receiver b and GPS satellite $s_G$ on frequency i is given by [17,18]:

$$\begin{array}{l}{\varphi }_{b,i}^{s_G}={\rho }_b^{s_G}+d{t}_b-d{t}^{s_G}+{\lambda }_i\left({\delta }_{b,i}^G+{\phi }_{b,i}-{\delta }_i^{s_G}-{\phi }_i^{s_G}+N_{b,i}^{s_G}\right)+{\epsilon }_{b,i}^{s_G}\\ P_{b,i}^{s_G}={\rho }_b^{s_G}+d{t}_b-d{t}^{s_G}+d_{b,i}^G-d_i^{s_G}+e_{b,i}^{s_G}\end{array}$$

(1)

where ϕ and P are the phase and code measurements in meters, respectively; ρ is the geometric range from satellite to receiver; $d{t}_b$ is the receiver clock bias; $d{t}^{s_G}$ is the satellite clock bias; ${\lambda }_i$ is the wavelength; ${\delta }_{b,i}^G$ and ${\delta }_i^{s_G}$ are the phase hardware delay in the receiver and satellite, respectively; ${\phi }_{b,i}$ and ${\phi }_i^{s_G}$ are the initial phase bias in the receiver and satellite, respectively; $d_{b_i}^G$ and $d_i^{s_G}$ are the code hardware delays in the receiver and satellite, respectively; N is the ambiguity; $\epsilon$ and $e$ are all other unmodelled errors including the measurement noise. Note that the subscript G means that ${\delta }_{b,i}^G$ and $d_{b,i}^G$ are related to the GPS system rather than a certain GPS satellite.

Similarly, the undifferenced phase and code observation equation for receiver b and satellite $s_{*}$ on frequency i is given by:

$$\begin{array}{l}{\varphi }_{b,i}^{s_{*}}={\rho }_b^{s_{*}}+d{t}_b-d{t}^{s_{*}}+{\lambda }_i\left({\delta }_{b,i}^{*}+{\phi }_{b,i}-{\delta }_i^{s_{*}}-{\phi }_i^{s_{*}}+N_{b,i}^{s_{*}}\right)+d_{TO}+{\epsilon }_{b,i}^{s_{*}}\\ P_{b,i}^{s_{*}}={\rho }_b^{s_{*}}+d{t}_b-d{t}^{s_{*}}+d_{b,i}^{*}-d_i^{s_{*}}+d_{TO}+e_{b,i}^{s_{*}}\end{array}$$

(2)

where $d_{TO}$ is the time offset between GPS and constellation *.

4.1.1. Intra-System Double-Differenced Observations with Loosely Combined Model (LCM)

The LCM assumes one satellite as the reference satellite for each constellation separately. All systematic errors, including the satellite-specific errors, the receiver-specific errors, and the time offset, are absent when double-differenced observations are created. The double-differenced observation equation corresponding to LCM wit $1_G$ as reference satellite for GPS and $1_{*}$ as reference satellite for constellation * is given by:

$$\begin{array}{l}{\varphi }_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+{\lambda }_i{N}_{br,i}^{1_G{s}_G}+{\epsilon }_{br,i}^{1_G{s}_G}\\ {\varphi }_{br,i}^{1_{*}{s}_{*}}={\rho }_{br}^{1_{*}{s}_{*}}+{\lambda }_i{N}_{br,i}^{1_{*}{s}_{*}}+{\epsilon }_{br,i}^{1_{*}{s}_{*}}\\ P_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+e_{br,i}^{1_G{s}_G}\\ P_{br,i}^{1_{*}{s}_{*}}={\rho }_{br}^{1_{*}{s}_{*}}+e_{br,i}^{1_{*}{s}_{*}}\end{array}$$

(3)

where $s_G\ge 2_G$ and $s_{*}\ge 2_{*}$.

4.1.2. Differential ISBs Estimation with Tightly Combined Model (TCM)

The TCM assumes a single-reference satellite for all the tracked constellations. When inter-system double-differenced observations are created between different constellations, additional phase and code differential ISBs are present and should be carefully considered. The double-differenced observation equation corresponding to TCM with ISBs estimation is given by [17,18,36]:

$$\begin{array}{l}{\varphi }_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+{\lambda }_i{N}_{br,i}^{1_G{s}_G}+{\epsilon }_{br,i}^{1_G{s}_G}\\ {\varphi }_{br,i}^{1_G{s}_{*}}={\rho }_{br}^{1_G{s}_{*}}+{\lambda }_i(N_{br,i}^{1_{*}{s}_{*}}+{\overline{\delta }}_{br,i}^{G*})+{\epsilon }_{br,i}^{1_G{s}_{*}}\\ P_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+e_{br,i}^{1_G{s}_G}\\ P_{br,i}^{1_G{s}_{*}}={\rho }_{br}^{1_G{s}_{*}}+d_{br,i}^{G*}+e_{br,i}^{1_G{s}_{*}}\end{array}$$

(4)

where $s_G\ge 2_G$ and $s_{*}\ge 1_{*}$, ${\overline{\delta }}_{br,i}^{G*}={\delta }_{br,i}^{G*}+N_{br,i}^{1_G{1}_{*}}$ is the lumped phase ISB parameter, which is composed of the phase ISB parameter and double-differenced ambiguity between reference satellites of GPS and constellation *, and $d_{br,i}^{G*}$ is the code ISB parameter. According to Equation (4), the phase and code ISBs can be estimated separately for each overlapping frequency and each system pair for single-epoch or multi-epoch solutions.

4.1.3. Differential ISBs Correction with TCM

Since the differential ISBs between overlapping frequencies are generally very stable in the time domain, they can be calibrated in advance and introduced as corrections to strengthen the adjustment model. The double-differenced observation equation corresponding to the correction of ISB parameters is then given by [17,18,36]:

$$\begin{array}{l}{\varphi }_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+{\lambda }_i{N}_{br,i}^{1_G{s}_G}+{\epsilon }_{br,i}^{1_G{s}_G}\\ {\widehat{\varphi }}_{br,i}^{1_G{s}_{*}}={\varphi }_{br,i}^{1_G{s}_{*}}-{\lambda }_i{\widehat{\delta }}_{br,i}^{G*}={\rho }_{br}^{1_G{s}_{*}}+{\lambda }_i{N}_{br,i}^{1_G{s}_{*}}+{\epsilon }_{br,i}^{1_G{s}_{*}}\\ P_{br,i}^{1_G{s}_G}={\rho }_{br}^{1_G{s}_G}+e_{br,i}^{1_G{s}_G}\\ {\widehat{P}}_{br,i}^{1_G{s}_{*}}=P_{br,i}^{1_G{s}_{*}}-{\widehat{d}}_{br,i}^{G*}={\rho }_{br}^{1_G{s}_{*}}+e_{br,i}^{1_G{s}_{*}}\end{array}$$

(5)

where $s_G\ge 2_G$ and $s_{*}\ge 1_{*}$, ${\widehat{\varphi }}_{br,i}^{1_G{s}_{*}}$ and ${\widehat{P}}_{br,i}^{1_G{s}_{*}}$ are ISBs corrected double-differenced phase and code observations, respectively. Compare to the loosely combined model in Equation (3), the redundancy is increased by 1 per overlapping frequency per system pair. Consequently, improved performance of ambiguity resolution and position estimation is expected with this tightly combined model.

Similar observation equations can also be obtained if a Galileo or BDS-3 satellite is chosen as reference satellite, which are not presented here.

4.1.4. Stochastic Model

During data processing procedure of ISB estimation and RTK positioning with Equations (3)–(5), the observational model is resolved with least squares adjustment and the exponential elevation weighting function is adopted as stochastic model [37,38]:

$$\sigma \left(\theta \right)={\sigma }_0\left(1+10\exp \left(\frac{-\theta }{10}\right)\right)$$

(6)

where $\sigma \left(\theta \right)$ is the standard deviation (STD) of the undifferenced observation; θ is the satellite elevation in degrees; σ₀ is the a-priori zenith-referenced code and phase STD, which are set as 3 mm and 0.3 m for phase and code respectively in this contribution.

4.2. Differential ISB estimation between BDS-3/GPS/Galileo B1C-L1-E1 and B2a-L5-E5a signals

In this section, we estimate phase and code differential ISBs between BDS-3 B1C/B2a and GPS/Galileo L1-E1/L5-E5a signals using short baseline (separated by approximately 3.85 m) data collected from 3–4 June 2018 at SGG of Wuhan University with the iGMAS and M300Pro receivers. According to Equation (4), the phase and code ISBs are estimated epoch-by-epoch with known receiver positions under 10° elevation cut-off angle.

Figure 7 shows estimated phase and code differential ISBs between identical receiver types (two iGMAS receivers) for GPS/BDS-3 L1-B1C/L5-B2a and Galileo/BDS-3 E1-B1C/E5a-B2a signals, respectively. As expected, the ISB series are indeed stable in the time domain and the average ISB values are close to zero. For GPS/BDS-3 L1-B1C/L5-B2a signals, the phase and code ISBs are −0.004/−0.001 cycles and −0.25/0.02 m, respectively, while for Galileo/BDS-3 E1-B1C and E5a-B2a signals, the values are 0/0 cycles and −0.25/0.03 m, respectively. Similar results are also observed for ISBs between two M300Pro receivers, although not shown here. The STDs of the mean phase ISB obtained in the single-epoch solutions varied from ∼0.015 to 0.018 cycle, which is at the same level of phase noise. With respect to the code ISBs, the STDs for L1-B1C and E1-B1C are present with, to some extent, relatively large values (0.83/0.9 m), which is much larger than L5-B2a and E5a-B2a (0.27/0.24 m). This is reasonable considering that the L1-E1-B1C code observations are susceptible to be affected by larger multipath effect and noise as shown in Section 3. Such results allow for neglecting these ISBs in the tightly combined RTK with BDS-3/GPS/Galileo observations.

Figure 8 shows the estimated phase and code ISBs with different receiver types (iGMAS and M300Pro receivers) for L1-B1C/L5-B2a and E1-B1C/E5a-B2a signals, respectively. We observe that the ISB series are indeed present with non-zero values but stable in the time domain. Compared with the above identical receiver case, the STDs for the phase and code ISBs are slightly larger than that in Figure 7. For GPS/BDS-3 L1-B1C/L5-B2a signals, the average phase and code ISBs are −0.272/0.500 cycles and 14.51/−90.28 m, respectively, while for Galileo/BDS-3 E1-B1C/E5a-B2a signals, the values are −0.252/0 cycles and −0.17/−89.65 m, respectively. Such results allow for calibrating and correcting these ISBs in advance.

4.3. Initial Assessment of Instantaneous Tightly Combined RTK with GPS, Galileo, and BDS-3 Operational Satellites

In this section, we present initial assessment of the tightly combined BDS-3/GPS/Galileo B1C-L1-E1/B2a-L5-E5a RTK using short baseline data in terms of the ADOP, ratio values, the empirical success rate, and positioning accuracy. The performance of the TCM with ISBs calibration and correction (cf., Equation (5)) will be evaluated in comparison with the LCM solution (cf., Equation (3)).

4.3.1. Evaluation Statistics

ADOP is an easy-to-compute scalar diagnostic to measure the intrinsic precision of the ambiguities, which is defined as [39]:

$$ADOP={\sqrt{\left|Q_{\widehat{a}}{}_{\widehat{a}}\right|}}^{\frac{1}{n}}\left(\mathrm{cy}cle\right)$$

(7)

where n the dimension of the ambiguity vector, $\left|Q_{\widehat{a}}{}_{\widehat{a}}\right|$ the determinant of the variance-covariance (vc)-matrix of float ambiguities. A lower ADOP corresponds to a higher precision of the ambiguities, and an ADOP smaller than 0.12 cycles corresponds to an ambiguity success rate larger than 99.9% [40].

Besides ADOP, the empirical success rate can also be used to verify the performance of ambiguity resolution, which is defined as the percentage of accepted and correctly-fixed epochs in total number of epochs [41]. The ambiguities are confirmed to be correct by comparing them with the true ambiguities, which are estimated based on all the available GPS/Galileo/BDS observations over the observation period and resolved using LAMBDA method [42]. The ratio, defined as the quadratic norm of a suboptimal solution divided by that of the optimal solution, is utilized to validate the correctness of fixed ambiguity solution. If the ratio value is larger than a given threshold, e.g., 2.0 or 3.0, then the fixed ambiguity solution is accepted. The Fixed Failure Rate Ratio Test with user defined failure rate (P_f) of 1% and 0.1% are used as well [43,44].

4.3.2. Test Results

Another short baseline experiment with different receiver types (iGMAS and M300Pro) was conducted on 15 June 2018 (from 07:00 to 21:00 GPST) to investigate the impact of ISB correction on ambiguity resolution and position estimation. The baseline is also separated by approximately 3.85 m as previous experiment that conducted from 3–4 June 2018 for ISB estimation in Section 4.2, but the involved receivers have been rebooted between data collection of the two experiments. The observations were collected with 10 s sampling interval and 10° elevation cut-off angle. During data processing, the position estimation and integer ambiguity estimation is conducted epoch-by-epoch. The phase and code ISBs are corrected in advance with their estimated values in Section 4.2.

Figure 9 shows the number of satellites tracked above 10° and 20° elevation cut-off angles together with their corresponding Position Dilution of Precision (PDOP) series. We observe that the tracked GPS Block IIF and BDS-3 operational satellites above 10° and 20° elevation cut-off angles are no less than one over the time span, and the available Galileo satellites are two to seven. The average tracked GPS, Galileo, and BDS-3 operational satellites during the day are 13.1 and 11.6, while the corresponding average PDOPs are 1.60 and 2.05, respectively.

Figure 10 and Figure 11 show the epoch-by-epoch ADOP of LCM and TCM with single-frequency (L1-E1-B1C) and dual-frequency (L1-E1-B1C/L5-E5a-B2a) observations under 10° elevation cut-off angle for GPS/BDS-3 and GPS/Galileo/BDS-3 solutions, respectively. We find that the single-epoch ADOP time series of TCM is obviously smaller than LCM for both GPS/BDS-3 and GPS/Galileo/BDS-3 solutions. In the case of GPS/BDS-3 combination, the mean ADOP values decrease from 0.397, 0.260 to 0.284, 0.179 cycles for single-frequency and dual-frequency, respectively, while for GPS/Galileo/BDS-3 solutions, the values decrease from 0.221, 0.128 to 0.155, 0.093 cycles. The improvement is more significant with single-frequency observations due to its relatively weaker model strength.

Figure 12 shows the epoch-by-epoch ratio values for LCM and TCM with single-frequency and dual-frequency observations under 10° elevation cut-off angle for GPS/BDS-3 (top panels) and GPS/Galileo/BDS-3 (bottom panels), respectively. The corresponding statistics is listed in

Table 4. Statistics of ratio values for GPS/BDS-3 and GPS/Galileo/BDS-3 solutions with 10° elevation cut-off angle. Number of epochs 5041.

Freq.	Ratio	GPS/BDS-3 (Epoch)			GPS/Galileo/BDS-3 (Epoch)
		LCM	TCM	Imp. 3	LCM	TCM	Imp.
SF 1	≥2.0	1941	3084	58.9%	4004	4805	20.0%
	≥3.0	1087	2114	94.5%	3108	4523	45.5%
	≥5.0	450	1133	151.8%	1782	3616	102.9%
DF 2	≥2.0	3607	4387	21.6%	4922	4978	1.1%
	≥3.0	3005	4015	33.6%	4800	4888	1.8%
	≥5.0	2088	3337	59.8%	4334	4602	6.2%

^1,2 “SF” and “DF” represent single-frequency and dual-frequency, respectively; “Imp.” means increase rate of epochs for TCM with respect to LCM.

. We find that the achieved ratio values of TCM are much higher than LCM. For example, the epochs with ratio values larger than 3.0 are increased from 1087, 3005 to 2114, 4015 with GPS/BDS-3 single-frequency and dual-frequency observations, respectively, while from 3108, 4800 to 4523, 4888 with GPS/Galileo/BDS-3 observations. Such results imply a higher empirical ambiguity resolution success rate for TCM.

Table 5. Empirical ambiguity resolution success rates (percent) for GPS/BDS-3 and GPS/Galileo/BDS-3 solutions with 10° and 20° elevation cut-off angles.

Freq.	Elev.	AVE. SAT.	AR Validation	GPS/BDS-3			GPS/Galileo/BDS-3
				LCM	TCM	Imp. 2	LCM	TCM	Imp.
SF	10°	(9.3, 13.1)	Th 1 = 2.0	33.2	58.7	25.5	78.4	95.1	16.7
			Th = 3.0	19.7	40.8	21.1	60.6	89.4	28.8
			Pf = 1%	9.8	39.5	29.7	71.8	97.5	25.7
			Pf = 0.1%	0.1	14.5	14.4	43.5	93.9	50.4
	20°	(8.0, 11.6)	Th = 2.0	24.1	47.3	23.2	66.8	91.4	24.6
			Th = 3.0	14.7	34.5	19.8	49.8	84.5	34.7
			Pf = 1%	3.3	22.3	19.0	51.0	93.2	42.2
			Pf = 0.1%	0.0	5.8	5.8	23.6	84.5	60.9
DF	10°	(9.3, 13.1)	Th = 2.0	68.7	85.9	17.2	97.6	98.7	1.1
			Th = 3.0	58.2	79.2	21.0	95.1	96.9	1.8
			Pf = 1%	63.2	84.8	21.6	99.2	100.0	0.8
			Pf = 0.1%	53.6	79.0	25.4	97.2	100.0	2.8
	20°	(8.0, 11.6)	Th = 2.0	59.1	80.1	21.0	95.1	96.9	1.8
			Th = 3.0	48.8	73.4	24.6	92.8	96.4	3.6
			Pf = 1%	50.6	78.9	28.3	96.3	97.5	1.2
			Pf = 0.1%	35.5	70.1	34.6	93.2	97.5	4.3

^1,2 “Th” represents threshold of ratio for integer ambiguity validation; “Imp.” means empirical success rate improvement of TCM with respect to LCM.

lists the instantaneous empirical ambiguity resolution success rates for GPS/BDS-3 and GPS/Galileo/BDS-3 solutions with different integer ambiguity validation strategies under 10° and 20° elevation cut-off angles. It’s shown that the observed satellites under 10° elevation cut-off angle are 9.3 and 13.1 for GPS/BDS-3 and GPS/Galileo/BDS-3, respectively, which is, to some extent, not promising. Therefore, the empirical ambiguity resolution success rates are not very satisfactory (e.g., only 33.2% and 78.4% with single-frequency observations and ratio threshold of 2.0). If TCM is adopted, improved empirical ambiguity resolution success rates are achieved with respect to the LCM whether for GPS/BDS-3 or GPS/Galileo/BDS-3 solutions. For positioning based on GPS/BDS-3 single-frequency observations, the empirical ambiguity resolution success rates of TCM are increased by approximately 25.5%, 21.1%, 29.7%, 14.4% and 23.2%, 19.8%, 19.0%, 5.8% with respect to LCM under 10° and 20° elevation cut-off angles, respectively. The improvement become approximately 16.7%, 28.8%, 25.7%, 50.4% and 24.6%, 34.7%, 42.2%, 60.9% if GPS/Galileo/BDS-3 single-frequency observations are adopted. Similar results are also obtained with GPS/BDS-3 dual-frequency observations. Nevertheless, the improvement is marginal (less than 5.0%) for GPS/Galileo/BDS-3 solutions with dual-frequency observations. This is reasonable since under this case, the LCM can already deliver excellent (more than 92.8%) ambiguity resolution success rates due to the relatively good observational geometry and strong model strength.

Table 6. STDs of float/correctly-fixed single-epoch RTK errors in the North, East, and Up.

Freq.	Model	10°			20°
		E (cm)	N (cm)	U (cm)	E (cm)	N (cm)	U (cm)
SF	GB 1 LCM	62.5/0.3	86.3/0.3	157.0/0.9	63.6/0.3	89.7/0.3	212.3/1.0
	GB TCM	58.2/0.3	82.6/0.3	156.0/0.9	60.6/0.3	85.6/0.4	199.0/1.0
	GEB 2 LCM	53.6/0.2	70.3/0.3	134.8/0.8	55.5/0.2	71.4/0.3	160.6/0.8
	GEB TCM	49.5/0.2	67.3/0.3	130.7/0.9	52.0/0.2	68.6/0.3	154.0/0.9
DF	GB LCM	54.5/0.2	75.0/0.3	137.6/0.8	55.2/0.2	76.4/0.3	177.3/0.8
	GB TCM	48.3/0.3	65.9/0.3	130.9/1.0	48.9/0.3	66.7/0.3	158.9/1.0
	GEB LCM	39.5/0.3	48.9/0.3	102.9/0.9	40.7/0.3	49.5/0.3	118.0/0.9
	GEB TCM	35.6/0.2	45.7/0.3	93.4/0.8	36.3/0.2	46.3/0.3	105.8/0.8

^1,2 “GB” and “GEB” represent GPS/BDS-3 and GPS/Galileo/BDS-3, respectively.

lists the empirical STDs of float and correctly-fixed (integer ambiguity validation with ratio threshold of 3.0) single-epoch positioning errors in the North, East, and Up under 10° and 20° elevation cut-off angles. We can observe that comparable results are achieved for correctly-fixed solutions whereas obvious improvements are obtained for float solutions if TCM is adopted. For example, in case of GPS/BDS-3 dual-frequency observations with 10° elevation cut-off angle, the three-dimension STDs are reduced from (54.5, 75.0, 137.6) to (48.3, 65.9, 130.9) cm, while for GPS/Galileo/BDS-3 solutions, they reduced from (39.5, 48.9, 102.9) to (35.6, 45.7, 93.4) cm. Although higher accuracy of the float solutions does not imply higher final positioning accuracy, they have an undeniable impact on the performance of the ambiguity resolution [23].

Figure 13 shows the horizontal and vertical positioning errors of different solutions with only single-frequency observations under 10° elevation cut-off angle. Compared with LCM solutions, it is obvious that the adoption of TCM can provide fewer float and wrongly-fixed solutions while more correctly-fixed solutions. Such results undoubtedly confirm the benefits of tightly combined RTK with observations from GPS, Galileo and BDS-3 operational satellites.

5. Conclusions

In this contribution, for the first time we provide initial assessment of tightly combined RTK with GPS, Galileo, and BDS-3 operational satellites using raw B1C-L1-E1/B2a-L5-E5a data collected with scientific geodetic iGMAS and commercial M300Pro receivers. The carrier-to-noise density ratio, pseudorange multipath and noise, as well as double-differenced carrier phase and code residuals of B1C/B2a signals are first characterized in comparison with L1-E1/L5-E5a signals. Our experimental results indicate that the observational quality of B1C/B2a signals that transmitted by BDS-3 operational satellites is comparable to their corresponding GPS/Galileo L1-E1/L5-E5a signals. Then, we investigate the size and stability of phase and code differential ISBs between B1C/B2a and L1-E1/L5-E5a signals using short baseline data collected with both identical and different receiver types. It is verified that the BDS-3/GPS/Galileo differential ISBs are indeed close to zero with two receivers of the same type. In contrast, they are generally present but stable in the time domain if receivers of different types are involved, which allows for calibration and correction in advance. Finally, a formal as well as empirical analysis is given to evaluate performance of instantaneous tightly combined BDS-3/GPS or BDS-3/GPS/Galileo RTK with single-frequency (B1C-L1-E1) and dual-frequency (B1C-L1-E1/B2a-L5-E5a) observations. Compared with the loosely combined model, the tightly combined model can provide us with much lower ADOP and higher ratio values. Meanwhile, the empirical ambiguity resolution success rate and positioning accuracy are remarkably improved as well, especially under constrained observational conditions with limited observed satellites and only single-frequency observations. In such cases, the improvement of the empirical ambiguity resolution success rate can reach up to approximately 10%∼60%, which further confirm the feasibility and benefits of tightly combined real-time kinematic (RTK) positioning with GPS, Galileo and BDS-3 operational satellites. In-depth investigation on the phase and code differential ISBs between BDS-3 operational satellites and GPS/Galileo for more types of receivers are still required in the future. Moreover, the ISBs between BDS-3 operational satellites with respect to QZSS, IRNSS, and SBAS will also need to be estimated so as to investigate the interoperability of BDS-3 with these systems.