Kinematic Absolute Positioning with Quad-Constellation GNSS

Abstract

The absolute positioning technique is based on a point positioning mode with a single Global Navigation Satellite System (GNSS) receiver, which has been widely used in many fields such as vehicle navigation and kinematic surveying. For a long period, this positioning technique mainly relies on a single GPS system. With the revitalization of Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, it is now feasible to carry out the absolute positioning with quad-constellation of GPS, GLONASS, BDS, and Galileo. A combination of multi-constellation observations can offer improved reliability, availability, and accuracy for position solutions. In this chapter, combined GPS/GLONASS/BDS/Galileo point positioning models for both traditional single point positioning (SPP) and precise point positioning (PPP) are presented, including their functional and stochastic components. The traditional SPP technique has a positioning accuracy at a meter level, whereas the PPP technique can reach an accuracy of a centimeter level. However, the later relies on the availability of precise ephemeris and needs a long convergence time. Experiments were carried out to assess the kinematic positioning performance in the two different modes. The positioning results are compared among different constellation combinations to demonstrate the advantages of quad-constellation GNSS.

1. Introduction

Position services have become an inevitable demand for the human activities. Advanced technologies of the position services can significantly improve human’s manufacturing efficiency, life quality, and resource utilization. Along with the development of human society, there is an increasing need of kinematic position services, such as automatic drive, intelligent transportation, precision agriculture, and so on. The Global Navigation Satellite System (GNSS), which rose in the 1980s of the last century, is an optimal infrastructure to realize the outdoor kinematic position services. The GNSS-based absolute positioning technologies have many advantages, such as no restriction by the inter-station distance, low cost, and simple data processing. The kinematic positions can be derived globally in all weather and any time based on the GNSS absolute positioning technique with a single receiver.

According to different performance demands, two kinds of kinematic absolute positioning technologies can be employed, namely, single point positioning (SPP) and precise point positioning (PPP). The SPP technology can provide meter-level positioning, while the PPP technology has a positioning accuracy at a centimeter level. As satellite-based positioning technologies, the performances of the SPP and PPP are quite dependent on the observed satellites. For a long period, the kinematic absolute positioning technologies are mainly based on a single GPS system. With the recent revitalization of the Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, the quad-constellation integrated absolute positioning has become feasible. Multi-constellation combination is expected to improve the reliability, availability, and accuracy of the SPP and PPP solutions due to the increased measurement redundancy and enhanced satellite geometry, especially when they are performed in areas with GNSS signal blockages.

2. Global Navigation Satellite System

A rapid development has been undergone for the satellite-based global navigation systems in recent years. The GNSS family has expanded from a single GPS constellation to four constellations of Galileo, BDS, GLONASS, and GPS. An overview of the four GNSS systems is conducted in terms of their space segment status and navigation signals.

2.1 Space segment status

A nearly circular orbit with an altitude of about 20,200 km is employed for GPS satellites. The GPS satellites pass a same place twice a day. All GPS satellites are located on six equally spaced orbital planes surrounding the earth. For each plane, there are four slots occupied by baseline satellites. The ascending nodes of the orbital planes are equally spaced 60° apart, and they are inclined at 55°. From virtually any point on the earth, users can view at least four GPS satellites, attributing to the 24-slot arrangement. Currently, a 27-slot constellation with improved coverage in most parts of the world is effectively operated for the GPS after expanding the constellation. The United States is committed to maintaining the availability of at least 24 operational GPS satellites 95% of the time. For the past few years, a total of 31 operational GPS satellites have been flying so as to ensure this commitment. The GPS constellation is a mix of old and new satellites, including 1 Block IIA, 11 Block IIR, 7 Block IIR-M, and 12 Block IIF satellites, as of March 2019 [1].

A complete revitalization of GLONASS with a full constellation including 24 operational satellites arranged into 3 orbital planes has been achieved since 2012. A nearly circular orbit is operated for each GLONASS satellite at an altitude of about 19,100 km, and approximately 11 h 16 min is needed for the GLONASS satellites to complete the orbit. The ascending nodes of orbital planes are separated by 120°, and each orbital plane has an inclination angle of 64.8°. The satellites within the same orbital plane are equally spaced 45° apart, while the difference in argument of latitude for satellites in equivalent slots in two different orbital planes is 15°. A continuous global navigation and position service can be provided due to the reasonable spacing of GLONASS satellites. Currently, there are 26 GLONASS satellites in orbit, but only 23 of them are in full operation, including 1 GLONASS-K1 and 22 GLONASS-M satellites. The GLONASS-M satellite SVN 716 is spare, and the SVN 720 satellite of the same series is in maintenance. The GLONASS-K1 satellite SVN 701 K is in the phase of flight tests [2].

Following the deployment timeline of BDS, its implementation has been performed in three steps: BeiDou navigation demonstration system (BDS-1) by 2000, regional BDS (BDS-2) by 2012, and global BDS (BDS-3) by 2020. Although the BDS-3 is under construction based on the “three-step” strategy, it has been providing basic services of positioning, navigation, and timing (PNT) for global customers since 27 December 2018. The nominal space constellation of BDS-2 consists of five geostationary earth orbit (GEO) satellites, five inclined geosynchronous orbit (IGSO) satellites, and four medium earth orbit (MEO) satellites. The nominal space constellation of BDS-3 consists of 3 GEO satellites, 3 IGSO satellites, and 24 MEO satellites. The five BDS-2 GEO satellites are located at 58.75°E, 80°E, 110.5°E, 140°E, and 160°E respectively, while the three BDS-3 GEO satellites are placed at 80°E, 110.5°E, and 140°E above the earth’s equator, respectively. Actually, spare satellites may be deployed in orbit, according to actual situation. An altitude of 35,786 km with a period of revolution of 23 h 56 min is adopted for the operation of the GEO satellites. The GEO satellites exhibit a non-zero inclination of 0.7−1.7°, as they are actively controlled in longitudes rather than latitudes. The altitude and orbital period of the IGSO satellites are the same as those of the GEO satellites, and they have an inclination of 55°. Each MEO satellite operates in a nearly circular orbit at an orbit inclination of 55° and an altitude of 21,528 km. The MEO satellites are arranged into three orbital planes, and an angle of 120° is used for the spacing between ascending nodes of different orbital planes. For MEO satellites, the period of the revolution is 12 h 53 min. As of March 2019, there are 15 BDS-2 satellites (5 GEO/7 IGSO/3 MEO), 5 BDS-3 demonstration system (BDS-3S) satellites (2 IGSO/3 MEO), and 19 BDS-3 satellites (1 GEO/18 MEO) in orbit. All BDS-3S satellites and the BDS-3 GEO satellite C59 are in the flight test phase, while the other BDS satellites are fully operational [3].

When fully deployed, the Galileo constellation will contain 30 satellites in 3 orbital planes. There are one inactive spare satellite and nine equally spaced operational satellites in each plane. The ascending nodes of the three planes are equally separated by 120°, and all of them are inclined at an angle of 56°. With a period of about 14 h 7 min and a semimajor axis of 29,600 km, all Galileo satellites are in a nearly circular orbit. The current Galileo space segment is composed of 26 satellites of 2 different generations. Two respective dual launches of four in-orbit validation (IOV) satellites in 2011 and 2012 initiated the buildup of the operational Galileo constellation. A permanent failure of the E5 and E6 signal transmission in May 2014 happened for the IOV-4 satellite due to a sudden power loss. Since then, the IOV-4 satellite can only transmit E1 signal. On 22 August 2014, the first pair of Full Operational Capability (FOC) satellites was launched. However, there was an “orbital injection anomaly” for the two FOC satellites 1 day later, which results in an elliptical orbit with an inclination roughly 5° smaller than planned. Because of the lack of broadcast ephemerides and single-frequency transmission, the IOV-4 satellite currently has the status “not available.” As of June 2016, the other three IOV satellites are declared “available,” namely, providing broadcast ephemerides and healthy signals that can be used in real-time navigation. Although the two FOC satellites FOC-1 and FOC-2 are not listed in the constellation status and in the eccentric orbit, they are generally transmitting broadcast ephemerides and navigation signals. The other 20 FOC satellites were successively declared “available” over the past 4 years [4].

Figure 1 shows a 24-h ground track of quad-constellation GNSS satellites on 18 March 2019. The quad-constellation mixed precise satellite orbit file is used to derive the satellite coordinates for all GNSS satellites. Different satellites are identified by different colors. The coordinate transformation from reference frame “IGS08” (three-dimensional positions x, y, and z), namely, the realization of international terrestrial reference frame 2008, to geodetic coordinate system (geodetic latitude B, longitude L, and height H) is first performed. Then, the tracks of GNSS satellites are projected to the grounds in the map based on their geodetic coordinates of latitudes and longitudes. The GNSS satellites can provide better PNT services for the stations near their ground tracks. All MEO satellites offer complete global coverage, including Galileo satellites, BDS MEO satellites, GLONASS satellites, and GPS satellites. The ground tracks of the four types of MEO satellites are confined from 57.2°S to 57.2°N latitude, from 56.1°S to 56.1°N latitude, from 65.6°S to 65.6°N latitude, and from 56.8°S to 56.8°N latitude, respectively. Due to the higher inclination angle of GLONASS satellites, they have the most extensive latitude coverage. More satellites at high-latitude areas will be visible by this extensive latitude coverage. The ground tracks of two Galileo FOC satellites in eccentric orbit (FOC-1 and FOC-2), which are represented by the thick lines, show notable asymmetric shape. For the two Galileo satellites, the orbital inclination is 5° lower than nominal. As a result, the peak values of the latitude coverage are reduced. Nevertheless, the scientific, geodetic, and surveying applications can still use the FOC-1 and FOC-2 satellites. The ground tracks of the BDS IGSO satellites are restricted from about 76.2°E to 138.0°E longitude and 57.5°S to 57.5°N latitude. Two figure-of-eight loops can be used to describe the IGSO satellite tracks. The two loops have an average longitude difference of around 30°, so as to effectively cover the western and eastern parts of China as well as the adjoining regions. The intersection points of the two loops are at the longitude of approximately 118°E and 95°E, respectively. The availability of satellites with high elevation angles can be improved by the employment of the inclined geosynchronous orbit. Thus, for users in densely populated areas, the “urban canyon” problems can be alleviated. As shown in Figure 1, the BDS GEO satellites have a movement within a range of 1.8°S–1.8°N latitude, but they are fixed in longitude. The south-north movement can be attributed to the non-zero inclination. To ensure enough visible satellites for users in Asian-Pacific regions, five GEO satellites are distributed in the Indian and Pacific oceans over the equator as supplements for the IGSO satellites [5].

Figure 1.

Ground tracks of quad-constellation navigation systems on 18 March 2019.

2.2 Navigation signals

The earlier navigation satellites provide signals on two frequencies so that the users can form dual-frequency observation combination to remove first-order effects of ionospheric delay based on the dispersive nature of the ionosphere. The modernized GNSS satellites have the capability of transmitting multi-frequency signals. Table 1 details the navigation signals of the four GNSS constellations [6]. k denotes the frequency factor of GLONASS satellites. The GPS Block IIF satellites can provide signals on L1, L2, and L5 frequencies, while the other GPS satellites are still transmitting L1/L2 signals. The last seven satellites of GLONASS-M series and all GLONASS-K satellites operate with three frequency bands, namely G1, G2, and G3, while the other GLONASS satellites can only offer G1 and G2 signals. All Galileo satellites are able to broadcast E1, E5A, E5B, E5 (A + B), and E6 signals. The BDS-2 satellites are capable of providing B1, B2, and B3 triple-frequency signals. In addition to the B1 and B3 signals, the BDS-3 satellites can transmit four new navigation signals, namely, B1C, B2a, B2b, and B2a + B2b [6].

Frequency (MHz)	GPS	GLONASS	Galileo	BDS-2	BDS-3
1602 + k × 9/16		G1
1575.42	L1		E1		B1C
1561.098				B1	B1
1278.75			E6
1268.52				B3	B3
1246 + k × 7/16		G2
1227.60	L2
1207.14			E5B	B2	B2b
1202.025		G3
1191.795			E5(A + B)		B2a + B2b
1176.45	L5		E5A		B2a

Table 1.

Navigation signals of quad-constellations.

3. Kinematic single point positioning with quad-constellations

With the use of single-frequency code observations and broadcast ephemeris, the SPP technology can provide meter-level positioning accuracy. Many researchers have focused on error mitigations to improve the SPP performance. The emerging multi-GNSS integration opens new prospects. In this section, the quad-constellation integrated SPP (QISPP) model with GPS, GLONASS, BDS, and Galileo measurements is developed, and its performance in the kinematic mode is evaluated.

3.1 QISPP model

Alignment of the coordinate and time references of the four GNSS systems is a key issue for the QISPP. With respect to the coordinate references, the coordinate systems of Galileo, BDS, GLONASS, and GPS satellites adopt the broadcast orbits of GTRF, CGCS2000, PZ90.11, and WGS-84, respectively. Although different coordinate references are employed for the four GNSS systems, the differences among them are only at a level of several centimeters [7, 8]. In view that the code-based positioning solutions using broadcast ephemeris can only achieve an accuracy at a meter level, such a small difference is negligible. In other words, the four GNSS systems can directly use their satellite coordinates without coordinate transformations in the QISPP. On the other hand, it is not the case for the time scales employed by the four GNSS systems. The GPS Master Control Station establishes the GPS time, which refers to the US Naval Observatory (USNO) Coordinated Universal Time (UTC) with a small difference of <1 μs. Besides, the UTC (USNO) is periodically corrected with integer leap seconds, and thus GPS time differs from it [1]. An atomic time scale UTC (Soviet Union, SU), which is maintained by Russia with an integer difference of 3 h and a fractional difference of <1 ms, is adopted by the GLONASS system [2]. Therefore, in addition to a tiny fractional difference, the GPS time differs from the GLONASS time by leap seconds. The BDS time system (BDT) was synchronized with UTC within 100 ns at 00:00:00 on 1 January 2006, and there exists a constant offset of 14 s between the GPS time and BDT [3]. Apart from a difference of 10 nanoseconds, the Galileo System Time (GST) is nearly identical to the GPS time [4]. The differences among the time references of the four GNSS systems will significantly affect the positioning solutions. Thus, unlike the coordinate reference frames, the inconsistent time scales cannot be ignored and must be properly handled in the QISPP.

Although there is only a physical clock in the multi-GNSS receiver, receiver clock parameters with respect to their respective time scales have to be estimated for each satellite system, since different time scales are adopted by the four GNSS systems. Alternatively, instead of adding a receiver clock parameter, a system time difference parameter with respect to a reference time scale can also be introduced [9]. We can directly estimate the GPS receiver clock offset as an unknown parameter, and the receiver clock offsets of the other satellite systems are regarded as a sum of the system time difference parameter and GPS receiver clock, provided that the GPS time scale is chosen as the reference. Following Pan et al. [10], the QISPP observation model reads

$$P^g={\rho }^g+{\mathit{cdt}}^g-{\mathit{cdT}}^g+d_{\mathit{orb}}^g+d_{\mathit{\text{trop}}}^g+d_{\mathit{ion}}^g+{\epsilon }_P^g$$ E1

$$P^r={\rho }^r+{\mathit{cdt}}^g+{\mathit{cdt}}_{\mathit{sys}}^{r,g}-{\mathit{cdT}}^r+d_{\mathit{orb}}^r+d_{\mathit{\text{trop}}}^r+d_{\mathit{ion}}^r+{\epsilon }_P^r$$ E2

$$P^b={\rho }^b+{\mathit{cdt}}^g+{\mathit{cdt}}_{\mathit{sys}}^{b,g}-{\mathit{cdT}}^b+d_{\mathit{orb}}^b+d_{\mathit{\text{trop}}}^b+d_{\mathit{ion}}^b+{\epsilon }_P^b$$ E3

$$P^e={\rho }^e+{\mathit{cdt}}^g+{\mathit{cdt}}_{\mathit{sys}}^{e,g}-{\mathit{cdT}}^e+d_{\mathit{orb}}^e+d_{\mathit{\text{trop}}}^e+d_{\mathit{ion}}^e+{\epsilon }_P^e$$ E4

where the superscripts e, b, r, and g refer to Galileo, BDS, GLONASS, and GPS satellites, respectively. P is the measured pseudorange; ρ is the geometric range; c is the speed of light in vacuum; dt^g is the GPS receiver clock offset; ${\mathit{dt}}_{\mathit{sys}}^{r,g}$ , ${\mathit{dt}}_{\mathit{sys}}^{b,g}$ and ${\mathit{dt}}_{\mathit{sys}}^{e,g}$ are the GPS-GLONASS, GPS-BDS, and GPS-Galileo system time differences, respectively; dT is the satellite clock offset; d_orb is the satellite orbit error; d_trop is the tropospheric delay; d_ion is the ionospheric delay; and ɛ is the measurement noise including multipath.

The broadcast ephemeris data is used to compute the clock offset and satellite position, as given in Eqs. (1)–(4). The Saastamoinen model is used to correct the tropospheric delay errors [11]. For the GPS, GLONASS, and BDS systems, the Klobuchar model is used to correct the ionospheric delay errors [12]. The second version of the NeQuick model is employed to perform the ionospheric error correction for the Galileo observations [13]. Regarding the Galileo ionospheric error corrections, the NeQuick model is better suited than the Klobuchar model [14]. Therefore, in the QISPP model, the unknown parameters to be estimated include three system time differences, one GPS receiver clock offset, and three receiver coordinates. Due to the low cost of single-frequency receivers, most SPP users are using them. Hence, in this study, we only use the single-frequency code observations of each system on the L1/G1/B1/E1 frequencies. The broadcast clock offsets and satellite orbits are referred to the ionosphere-free code combination on two frequencies. Consequently, the hardware delay biases in a form of ionosphere-free combination are contained in the satellite clock offsets derived from the broadcast ephemeris. When the ionosphere-free combined code observables are used, the hardware delay biases can be canceled out for dual-frequency users. But for single-frequency users, the hardware delay biases must be corrected. Fortunately, the broadcast navigation messages on a satellite-by-satellite basis have provided the time group delays, which can be employed to carry out the hardware delay bias corrections in the single-frequency pseudorange-based positioning.

Following Hoque et al. [15], the Klobuchar model for ionospheric error corrections of single-frequency users can be described as

$$d_{\mathit{ion}}=\left\{\begin{array}{ll}F\times \left[5.0\times {10}^{-9}\times \mathit{AMP}\times \left(1-\frac{x^2}{2}+\frac{x^4}{24}\right)\right]& \left|x\right|<1.57\\ F\times \left(5.0\times {10}^{-9}\right)& \left|x\right|\ge 1.57\end{array}\right.$$ E5

$$x=\frac{2\pi \left(t-50400\right)}{\mathit{PER}}$$ E6

$$\mathit{AMP}=\left\{\begin{array}{ll}\sum _{n=0}^3{\alpha }_n{\varphi }_m^n& \mathit{AMP}\ge 0\\ 0& \mathit{AMP}<0\end{array}\right.$$ E7

$$\mathit{PER}=\left\{\begin{array}{ll}\sum _{n=0}^3{\beta }_n{\varphi }_m^n& \mathit{PER}\ge 72,000\\ 72,000& \mathit{PER}<72,000\end{array}\right.$$ E8

where α and β are coefficients included as part of satellite message, ϕ_m is the geomagnetic latitude, t is the local time, and E is the elevation angle. For more details, refer to Klobuchar [12].

The second version of the NeQuick model can be depicted by several main equations provided below:

$$N_{\mathrm{bot}}\left(h\right)=N_{\mathrm{E}}\left(h\right)+N_{\mathrm{F}1}\left(h\right)+N_{\mathrm{F}2}\left(h\right)$$ E9

$$N_{\mathrm{E}}\left(h\right)=\frac{4\mathit{NmE}}{{\left(1+\text{exp}\left(\frac{h-\mathit{hmE}}{\mathit{BE}}\xi \left(h\right)\right)\right)}^2}\times \text{exp}\left(\frac{h-\mathit{hmE}}{\mathit{BE}}\xi \left(h\right)\right)$$ E10

$$N_{\mathrm{F}1}\left(h\right)=\frac{4\mathit{NmF}1}{{\left(1+\text{exp}\left(\frac{h-\mathit{hmF}1}{B1}\xi \left(h\right)\right)\right)}^2}\times \text{exp}\left(\frac{h-\mathit{hmF}1}{B1}\xi \left(h\right)\right)$$ E11

$$N_{\mathrm{F}2}\left(h\right)=\frac{4\mathit{NmF}2}{{\left(1+\text{exp}\left(\frac{h-\mathit{hmF}2}{B2}\right)\right)}^2}\times \text{exp}\left(\frac{h-\mathit{hmF}2}{B2}\right)$$ E12

$$N\left(h\right)=\frac{4\mathit{NmF}2}{{\left(1+\text{exp}\left(z\right)\right)}^2}\text{exp}\left(z\right)$$ E13

For more details and specific meaning of the above parameters, refer to Nava et al. [13].

As to the stochastic model for the QISPP, we can use the following covariance matrix of observations:

$$\mathit{Cov}=\left[\begin{array}{llll}{Q}_g& 0& 0& 0\\ 0& Q_r& 0& 0\\ 0& 0& Q_b& 0\\ 0& 0& 0& Q_e\end{array}\right]$$ E14

$$Q=\left[\begin{array}{llll}{\sigma }_1^2& 0& \cdots & 0\\ 0& {\sigma }_2^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_n^2\end{array}\right]$$ E15

where the subscript n denotes the number of satellites for each satellite system and σ² denotes the variance of code observations, which can be written as

$${\sigma }^2={\sigma }_0^2/{\left(\text{sin}E\right)}^2$$ E16

where σ₀ denotes the standard deviation (STD) of code observations, which differs among different satellite systems, and E denotes the satellite elevation angle.

The code observation precision is set to 0.3 m for GPS satellites [16]. Following Cai et al. [17], an initial weight ratio of 1:1 is appropriate for BDS and GPS code observations. Thus, the code observation precision is also set to 0.3 m for BDS satellites. The code observation precision is set to 0.6 m for GLONASS satellites due to their twice lower code chipping rate than the GPS code observations. The code observations of Galileo satellites are down-weighted by a factor of four, considering that the broadcast ephemeris has relatively lower accuracies [8]. That is, the precision of the Galileo code observations is also set to 0.6 m.

3.2 Performance analysis of QISPP solutions

Four different constellation combinations are employed for the purpose of comparison, namely, GPS/GLONASS/BDS/Galileo, GPS/GLONASS/BDS, GPS/GLONASS, and GPS-only. In the data processing, the receiver coordinates as well as other unknown parameters are estimated epoch-by-epoch without imposing any constraints between the epochs in order to analyze the single-epoch SPP performance. For brevity, in the following figures and tables, “GLO” and “GAL” are used to represent GLONASS and Galileo systems, respectively.

Figure 2 shows the epoch-wise SPP positioning errors for the four different combination cases at station NNOR on 8 April 2015. In the east, north, and up directions, the variations of positioning errors are consistent for the four cases, but the series of position errors show less fluctuation for the triple- and quad-constellation cases. It is seen that the GPS-only SPP achieves larger positioning errors than the GPS/GLONASS case. The positioning errors of GPS/GLONASS/BDS SPP at almost all epochs are further reduced by combining with BDS. The further introduction of Galileo observations does not exhibit significant change, since the blue lines are almost completely covered by the orange ones.

Figure 2.

Kinematic positioning errors of four different combination cases for the SPP processing at NNOR.

Figure 3 presents the position dilution of precision (PDOP) and the number of visible satellites for the four cases. It is clear that the number of visible satellites is obviously increased by the multi-constellation combination and the PDOP value is simultaneously decreased. The quad-constellation integrated case increases the average number of visible satellites from 6.9 to 27.8 in comparison to the GPS-only SPP, and thus the average PDOP values are significantly decreased from 2.3 to 1.1. The increased number of available satellites and decreased PDOP values explain why the positioning accuracies can be improved in the multi-constellation integrated cases.

Figure 3.

Number of satellites and PDOP values for the SPP processing at NNOR.

The datasets from 47 International GNSS Service (IGS) stations are processed for further analysis, and the average root mean square (RMS) statistics of epoch-wise SPP positioning errors, number of satellites, and PDOP are obtained. The average satellite numbers of the above four processing cases are 8.6, 15.6, 21.1, and 23.4, respectively, while the corresponding average PDOPs are 2.0, 1.4, 1.2, and 1.1, respectively. The combination of GPS and GLONASS improves the positioning accuracy over the GPS-only case by 7, 5, and 5% from 0.81, 2.05, and 3.13 m to 0.75, 1.94, and 2.96 m in the east, north, and up directions, respectively. In the triple-constellation SPP, an accuracy improvement of 9, 6, and 7% over the GPS/GLONASS case to 0.68, 1.82, and 2.75 m in the three directions is achieved, respectively. After a further integration with Galileo, the positioning accuracy is only improved by 3, 2, and 2% to 0.66, 1.78, and 2.70 m in the three directions, respectively [10].

4. Kinematic precise point positioning with quad-constellations

The PPP technique adopts an absolute positioning approach to achieve centimeter-level positioning accuracy using code and carrier phase observations as well as precise satellite orbit and clock offset corrections [18]. Both industrial applications and scientific research widely use GPS-based PPP. But the position solutions of GPS-based PPP require a long time to converge. Due to high measurement redundancy, significantly reduced convergence time and improved positioning accuracy can be expected by using multi-constellation GNSS PPP. The positioning model and processing method of quad-constellation integrated PPP (QIPPP) with GPS, GLONASS, BDS, and Galileo measurements are developed, and then the improvement in precision and convergence time from quad-constellations is verified by comparing the solutions of different constellation combinations.

4.1 QIPPP model

In view that the code observation equation is detailed in Eqs. (1)–(4), only the carrier phase observations on ith (i = 1, 2) frequency are formulated, that is:

$${\mathrm{\Phi }}_i=\rho +\mathit{cdt}-\mathit{cdT}+d_{\mathit{orb}}+d_{\mathit{\text{trop}}}-d_{\mathit{ion}/L_i}+{\lambda }_i{N}_i+{\epsilon }_{{\mathrm{\Phi }}_i}$$ E17

where Φ is the measured carrier phase, d_ion/Li is the ionospheric delay on ith frequency, λ is the wavelength, and N is the phase ambiguity term grouped with hardware delays. It should be noted that the wavelength λ is different for different GLONASS satellites because GLONASS employs the frequency-division multiple access (FDMA) technique.

To remove the first-order effects of ionospheric delays, the PPP normally utilizes the ionosphere-free (IF) combined observables, that is:

$$P_{\mathit{\text{IF}}}=\left(f_1^2\cdot P_1-f_2^2\cdot P_2\right)/\left(f_1^2-f_2^2\right)$$ E18

$${\mathrm{\Phi }}_{\mathit{\text{IF}}}=\left(f_1^2\cdot {\mathrm{\Phi }}_1-f_2^2\cdot {\mathrm{\Phi }}_2\right)/\left(f_1^2-f_2^2\right)$$ E19

where Φ_IF and P_IF denote the IF combined carrier phase and code observables, respectively, and f denotes the carrier phase frequency, which differs among GLONASS satellites.

In order to investigate emerging new satellite systems such as Galileo and BDS, the Multi-GNSS Experiment (MGEX) has been established by the IGS [19]. The correction of satellite clock and orbit errors for QIPPP is conducted with the use of precise satellite clock and orbit products provided by the MGEX. The first-order effects of ionospheric delays are removed using the IF combined observables, as shown in Eqs. (18) and (19). The tropospheric delays can be divided into a wet part and a dry part [20]. The wet part is estimated from the measurements, while the Hopfield tropospheric model is employed to correct the dry part. The projection from slant delays to zenith delays adopts the Niell mapping functions [21]. As some literatures such as Kouba and Héroux [18] have well described other error mitigations, they are not provided here. In PPP, the code-specific hardware delays at the receiver and the receiver clock offsets are usually estimated as a lumped term, as they are linearly correlated with each other. For different navigation systems, both the frequency and signal structures differ. Consequently, within a multi-GNSS receiver, the receiver-dependent code hardware delays are different for the four navigation systems. For the purpose of solving this issue, we should design a receiver clock offset parameter for each satellite system. Alternatively, the differences between receiver clock estimates of different satellite systems can be compensated by introducing an inter-system bias (ISB). The QIPPP observation model can be written as follows, provided that the GPS system is chosen as the reference [22]:

$$P_{\mathit{\text{IF}}}^g={\rho }^g+\mathit{cdt}+M^g{d}_{\mathit{zwd}}+{\epsilon }_{P_{\mathit{\text{IF}}}}^g$$ E20

$${\mathrm{\Phi }}_{\mathit{\text{IF}}}^g={\rho }^g+\mathit{cdt}+M^g{d}_{\mathit{zwd}}+N_{\mathit{\text{IF}}}^g+{\epsilon }_{{\mathrm{\Phi }}_{\mathit{\text{IF}}}}^g$$ E21

$$P_{\mathit{\text{IF}}}^b={\rho }^b+\mathit{cdt}+{\mathrm{ISB}}_{b,g}+M^b{d}_{\mathit{zwd}}+{\epsilon }_{P_{\mathit{\text{IF}}}}^b$$ E22

$${\mathrm{\Phi }}_{\mathit{\text{IF}}}^b={\rho }^b+\mathit{cdt}+{\mathrm{ISB}}_{b,g}+M^b{d}_{\mathit{zwd}}+N_{\mathit{\text{IF}}}^b+{\epsilon }_{{\mathrm{\Phi }}_{\mathit{\text{IF}}}}^b$$ E23

$$P_{\mathit{\text{IF}}}^r={\rho }^r+\mathit{cdt}+{\mathrm{ISB}}_{r,g}+M^r{d}_{\mathit{zwd}}+{\epsilon }_{P_{\mathit{\text{IF}}}}^r$$ E24

$${\mathrm{\Phi }}_{\mathit{\text{IF}}}^r={\rho }^r+\mathit{cdt}+{\mathrm{ISB}}_{r,g}+M^r{d}_{\mathit{zwd}}+N_{\mathit{\text{IF}}}^r+{\epsilon }_{{\mathrm{\Phi }}_{\mathit{\text{IF}}}}^r$$ E25

$$P_{\mathit{\text{IF}}}^e={\rho }^e+\mathit{cdt}+{\mathrm{ISB}}_{e,g}+M^e{d}_{\mathit{zwd}}+{\epsilon }_{P_{\mathit{\text{IF}}}}^e$$ E26

$${\mathrm{\Phi }}_{\mathit{\text{IF}}}^e={\rho }^e+\mathit{cdt}+{\mathrm{ISB}}_{e,g}+M^e{d}_{\mathit{zwd}}+N_{\mathit{\text{IF}}}^e+{\epsilon }_{{\mathrm{\Phi }}_{\mathit{\text{IF}}}}^e$$ E27

where ISB_b,g, ISB_r,g, and ISB_e,g are the GPS-BDS, GPS-GLONASS, and GPS-Galileo inter-system biases, respectively, d_zwd is the tropospheric zenith wet delay (ZWD), and M is the tropospheric mapping function.

For the QIPPP processing, we employ a Kalman filter approach. Actually, the geometric range ρ in Eqs. (20)–(27) is a function of receiver coordinates and satellite coordinates. After the precise ephemeris data is used to determine the satellite coordinates and the ρ is linearized, the unknown parameters include phase ambiguity parameters equal to the number of the observed GNSS satellites, one ZWD, three ISB, one receiver clock offset, and three receiver coordinates. We should also provide appropriate dynamic models for the state vector and proper stochastic models for the measurements in the Kalman filter. The IF combined phase and code observables are obtained using the raw phase and code measurements on two different frequencies in a form of linear combination. It is assumed that the measurements on different frequencies are independent. Based on the law of random error propagation [23], we can derive the initial variances of the IF combined observables. The satellite elevation angles can be further introduced to acquire the actual variances [24]. As for the dynamic models, the phase ambiguity parameters can be modeled as constants, whereas the ZWD, ISB, receiver clock offset, and kinematic receiver coordinates may be modeled as a random walk (RW) process [25, 26, 27].

4.2 Performance analysis of QIPPP solutions

The datasets from stations SEYG, JFNG, and MAR7 on 1 March 2019 are used for numerical analysis in this section. The three stations are located at low-, middle-, and high-latitude regions, respectively, and all of them are able to provide multi-constellation observations. The in-house MIPS-PPP software capable of processing quad-system observation data as well as single-system measurements of Galileo, BDS, GLONASS, and GPS, which is developed at Central South University, China, is employed for the quad-constellation GNSS PPP processing. Regarding the specific PPP position determination, the Galileo E1/E5A, BDS B1/B2, GLONASS G1/G2, and GPS L1/L2 dual-frequency observations are used. The cutoff satellite angle is set to 10°, while the observations are recorded at a sampling rate of 30 s. The IGS Analysis Center, German Research Centre for Geosciences (GFZ), has been generating and releasing the mixed multi-GNSS final precise satellite clock offset and orbit products with a sampling interval of 30 s and 5 min, respectively, and the mitigation of the satellite clock and orbit errors is carried out using them in this study. For the ISB, receiver clock offset, and ZWD parameters, the spectral density values are set to 10⁻⁷, 10⁵, and 10⁻⁹ m²/s, respectively [28]. The initial STD values for GPS and GLONASS code observations are set to 0.3 and 0.6 m, while for phase observations they are both set to 2 mm [28]. Since there is a relatively lower accuracy for the satellite orbits and clocks of Galileo and BDS [29, 30], their observations are down-weighted with a factor of four. That is, for both Galileo and BDS, the code and phase observation accuracies are set to 0.6 m and 4 mm, respectively. Kinematic processing is made on an epoch-by-epoch basis using the static data. No constraints between epochs are imposed so as to simulate kinematic situations. In the Kalman filtering, the coordinates of the dynamic receiver are modeled as a RW process, and the spectral density is set as 10² m²/s.

Figure 4 shows the positioning errors of GPS, GPS/BDS, GPS/BDS/GLONASS, and four-system PPP in three directions of east, north, and up in the simulated kinematic test. It can be seen from Figure 4 that, compared with the PPP solutions of GPS-only system, the error curve of the PPP of multi-constellation combinations converges to the stable value faster in the east, north, and up directions. For all processing schemes, the positioning errors in the vertical direction are larger than those of horizontal directions.

Figure 4.

PPP kinematic positioning errors at stations SEYG, JFNG, and MAR7 for four different processing cases.

In order to assess the kinematic positioning accuracy, Table 2 provides the RMS statistical values using the position errors in the last 1 h in which the position solutions in all three components have reached stable values. The results show that, taking JFNG station as an example, the positioning accuracy of the GPS-only PPP in three directions is 0.160, 0.097, and 0.192 m, respectively. After the combination of GPS and BDS, compared with the single GPS system, the positioning accuracy is improved by 11, 55, and 18%, respectively. Compared with the GPS/BDS PPP, the positioning accuracy of the combination of three systems is significantly improved. After further adding Galileo observations, the three-dimensional (3D) accuracy of PPP of the four-system combination is slightly improved by 1.2 cm. Table 3 shows the convergence time in three directions. Taking the JFNG station as an example, the PPP solutions of the four-system combination requires 22.5, 28.5, and 79.0 min to converge to the accuracy level of 1 dm. Compared with the single- and dual-system cases, the convergence time in the horizontal directions of the triple- and quad-system integrated PPP is dramatically shortened. In view of the fact that the positioning accuracy of kinematic PPP in the vertical direction is worse than that of the horizontal directions, the convergence standard may be too strict for the vertical direction, which even leads to the failure of effective convergence in a short period of time sometimes.

Stations	Directions	GPS	GPS/BDS	GPS/BDS/GLO	GPS/BDS/GLO/GAL
SEYG	East	0.090	0.061	0.058	0.043
	North	0.062	0.057	0.028	0.023
	Up	0.158	0.126	0.082	0.071
	3D	0.192	0.151	0.104	0.086
JFNG	East	0.160	0.142	0.052	0.050
	North	0.097	0.044	0.020	0.020
	Up	0.192	0.158	0.119	0.106
	3D	0.268	0.217	0.131	0.119
MAR7	East	0.066	0.054	0.040	0.032
	North	0.045	0.035	0.028	0.021
	Up	0.081	0.068	0.034	0.030
	3D	0.114	0.094	0.060	0.049

Table 2.

Convergence accuracy of kinematic PPP (unit: m).

Stations	Directions	GPS	GPS/BDS	GPS/BDS/GLO	GPS/BDS/GLO/GAL
SEYG	East	70.5	64.0	51.5	51.5
	North	37.5	33.0	15.0	14.5
	Up	90.0	88.0	77.0	76.5
JFNG	East	94.0	85.0	23.0	22.5
	North	86.0	79.5	28.5	28.5
	Up	88.0	82.5	79.0	79.0
MAR7	East	79.0	78.0	68.5	68.5
	North	43.0	41.0	23.5	23.5
	Up	86.0	84.5	61.0	60.0

Table 3.

Convergence time of kinematic PPP (unit: min).

5. Conclusions

The GNSS-based absolute positioning technologies can provide reliable kinematic position services anywhere, in all weather, and anytime using a single receiver. With single-frequency code measurements and broadcast satellite ephemeris, the SPP technology can provide meter-level positioning accuracy. With dual-frequency code and carrier phase measurements as well as precise satellite orbit and clock products, the PPP technology can offer centimeter-level positioning accuracy. In recent years, the satellite systems have been booming. In view that both SPP and PPP belong to the satellite-based kinematic absolute positioning technologies, the multi-constellation combination provides new prospects for their performance improvement, due to more visible satellites, increased measurement redundancy, and enhanced satellite sky distribution. The quad-constellation integrated SPP and PPP models with GPS, GLONASS, BDS, and Galileo measurements are developed, respectively. The results indicate significantly improved positioning performance of the multi-GNSS integration, which will further promote the applications of SPP and PPP technologies.

Acknowledgments

The contribution of data and products from IGS is appreciated.

Conflict of interest

The authors declare no conflict of interest.

List of abbreviations

BDS

BeiDou Navigation Satellite System

BDS-1

BeiDou navigation demonstration system

BDS-2

regional BDS

BDS-3

global BDS

BDS-3S

BDS-3 demonstration system

BDT

BDS time system

FDMA

frequency-division multiple access

FOC

full operational capability

GEO

geostationary earth orbit

GFZ

German Research Centre for Geosciences

GLONASS

Global Navigation Satellite System

GNSS

Global Navigation Satellite System

GST

Galileo System Time

IF

ionosphere-free

IGS

International GNSS Service

IGSO

inclined geosynchronous orbit

IOV

in-orbit validation

ISB

inter-system bias

MEO

medium earth orbit

MGEX

multi-GNSS experiment

PDOP

position dilution of precision

PNT

positioning, navigation, and timing

PPP

precise point positioning

QIPPP

quad-constellation integrated PPP

QISPP

quad-constellation integrated SPP

RMS

root mean square

RW

random walk

SPP

single point positioning

STD

standard deviation

SU

Soviet Union

USNO

United States Naval Observatory

UTC

Coordinated Universal Time

ZWD

zenith wet delay

3D

three-dimensional