**Guorui Xiao** 

# **Multi-frequency and multi-GNSS PPP**

**phase bias estimation and ambiguity resolution**

**Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution**

**Guorui Xiao** 

Guorui Xiao

Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution

## Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution

by Guorui Xiao

Karlsruher Institut für Technologie Geodätisches Institut

Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution

Zur Erlangung des akademischen Grades eines Doktor-Ingenieurs von der KIT-Fakultät für Bauingenieur-, Geo- und Umweltwissenschaften des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Guorui Xiao, M. Eng.

Tag der mündlichen Prüfung: 4. Juli 2019 Referent: Prof. Dr.-Ing. Dr. h.c. Bernhard Heck, Karlsruhe Institute of Technology Korreferent: Prof. Dr.-Ing. habil. Lambert Wanninger, Technical University of Dresden

#### **Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2022 – Gedruckt auf FSC-zertifiziertem Papier

ISBN 978-3-7315-1055-0 DOI 10.5445/KSP/1000124935

## **Abstract**

GPS PPP has found many scientific and industrial applications due to its costeffectiveness, global coverage, and high accuracy for decades. However, it suffers from a few drawbacks which limites further applications, e.g., slow convergence and unresolved ambiguity. The rapid development of multi-GNSS and multi-frequency signals offers exciting prospects for further improvements. In this thesis, multifrequency and multi-GNSS measurements are properly integrated for PPP with ambiguity resolution to achieve the state-of-the-art fast and accurate positioning, which provides an important contribution to GNSS precise positioning and applications.

Integer ambiguity resolution is the key issue for improving PPP accuracy and convergence. The core of PPP ambiguity resolution is carrier phase fractional cycle biases (FCB) estimation. In this thesis, the characteristics of phase FCB are analyzed with the recent and long temporal series. Results indicate that the temporal stabilities have been significantly improved in recent years. Taking advantage of this property, an improved FCB estimation method based on Kalman filter is proposed. The designed Kalman filter significantly reduces the dimension of the involved matrix and accelerates FCB computation, which outperforms the commonly used least square method in terms of efficiency. The method is especially useful for real-time application.

With the estimated FCB products, the PPP ambiguity resolution with the current Galileo and BDS constellation is verified. The satellite FCB is thoroughly assessed by a comparison with that of GPS in terms of data usage rate, residual distribution, as well as standard deviation of daily estimates. Results indicate that the quality of Galileo widelane (WL) FCB is better than GPS and BDS, while the quality of Galileo narrow-lane (NL) FCB is slightly worse than that of GPS but better than that of BDS. Within the Galileo constellation, the performance of In-Orbit Validation (IOV) satellites WL FCB is worse than that of Full Operational Capability (FOC) satellites as a result of a reduction in the power of the transmitted signal. The WL FCB performance of the two highly eccentric satellites is comparable to other FOC satellites. The NL FCB quality of FOC, IOV (except E19), as well as the two eccentric satellites, shows no significant difference in terms of data usage rates and residuals. Solution of PPP ambiguity resolution demonstrates that the Galileo and BDS observations can bring an obvious benefit to GPS only PPP ambiguity resolution.

PPP ambiguity resolution is also extended to multi-frequency observations in this thesis. In order to properly integrate the multi-frequency observations, a unified uncombined PPP ambiguity resolution model based on raw observations is proposed. Based on the unified model, the FCBs generated from multi-frequency observations can be flexibly used, such as for dual- or triple- frequency PPP AR. Its efficiency is verified with Galileo and BeiDou triple-frequency observations collected from globally distributed MGEX stations. The estimated FCB are assessed with respect to residual distributions and standard deviations, which shows a good consistency with the input float ambiguities. The standard deviation of linear combined FCB is much smaller than that of raw FCB, which indicates linear combined FCB is more efficient for broadcasting in real time applications. Triple-frequency PPP ambiguity resolution indicates that the positional biases are significantly reduced compared with the float solutions. The improvements are 49.2%, 38.3%, and 29.6%, respectively, in east/north/up components for positioning with BDS, while the corresponding improvements are 60.0%, 29.0%, and 21.1% for positioning with Galileo.

The uncombined signal processing brings new challenges for cycle slip detection. The most significant feature is carrier frequency identification of cycle slips. Since all carrier frequency observations are processed separately, it is essential to identify the carrier frequency of the cycle slips to avoid contaminating other observations. To provide continuous carrier phase measurements for the above research, an improved approach based on a time-differenced model for cycle slip detection and repair is proposed, which reduce the false alarms and increase the success rate of cycle slip estimation.

In summary, the GPS only PPP has been extended to multi-frequency and multi-GNSS PPP ambiguity resolution with improved accuracy and fast convergence. This is accomplished by a unified model based on the uncombined PPP. The proposed model has been carefully studied and enriched with improved cycle slip detection and repair, fast FCB estimation method and ambiguity resolution.

This article-based (cumulative) thesis consists of a detailed introductory chapter and four chapters associated with the following peer-reviewed publications:


## **Zusammenfassung**

Aufgrund seiner Wirtschaftlichkeit, globalen Abdeckung und hohen Genauigkeit hat GPS PPP seit Jahrzehnten. PPP weist jedoch einige Nachteile auf, die weitere Anwendungen einschränken, z. B. langsame Konvergenzzeit und ungelöste Mehrdeutigkeit. Die schnelle Entwicklung von GNSS und seinen Mehrfrequenz-Signalen bietet spannende Perspektiven für weitere Verbesserungen. In dieser Dissertation sind Multi-GNSS- und Mehrfrequenz-Beobachtungen für PPP-Lösung mit Mehrdeutigkeitauflösung integriert, um eine schnelle und genauige Positionsbestimmung auf dem neuesten Stand der Technik zu erreichen. Dies ist ein wichtiger Beitrag zur präzisen Positionsbestimmung mittels GNSS und seinen Anwendungen.

Festsetzung der Phasenmehrdeutigkeit auf ihre ganzzahlige Werte ist von großer Bedeutung für die Verbesserung der PPP-Genauigkeit und -Kovergenzzeit. Der Schwerpunkt in diesem Zusammenhang ist die Schätzung der (Fractional Cycle Biases) FCB der Trägerphasen. In dieser Dissertation-werden die Eigenschaften der Phasen-FCB mit kürzlichen und langzeitlichen Serien analysiert. Die Ergebnisse zeigen, dass die zeitliche Stabilität in den letzten Jahren wesentlich verbessert wird. Bei Nutzung dieser Eigenschaft wird ein verbessertes FCB-Schätzverfahren mit Kalman-Filter entwickelt. Dieser entwicklende Kalman-Filter reduziert die Dimension der betroffenen Matrix und beschleunigt damit die Berechung des FCB. Allerdings verliert die traditionale LS Methode hinsichtlich der Effizienz im Vergleich dazu. diese Methdode eignet sich besondersfür die Echtzeitiganwendungen.

Mit den geschätzten FCB-Produkten wird die PPP–Mehrdeutigkeitauflösung mit der aktuellen Konstellation von Galileo und BDS überprüft. DerStatelliten-FCB wird im Vergleich mit GPS gewonnen Wertenhinsichtlich der Datennutzungsrate, der Verteilung der Residuen und der täglichen Abschätzungen der Standardabweichung bewertet. Den Ergebnisse zeigen, dass die Qualität der Galileo Wide-Lane (WL) FCB besser als GPS und BDS ist. Mittlerweile ist die Qualität der Galileo Narrow-Lane (NL) FCB ein bisschen schlechter als GPS aber besser als BDS. innerhalb der Galileo Konstellation sind die Leistung des In-Orbit-Validation (IOV)-Satelliten bezüglich WL FCB schlechter als die Full Operational Capability (FOC) Stalliten wegen der Dämpfung der Energie von den austrahlenden Signalen. Die WL FCB Leistung der Satelliten mit hohem Exzentrizität ist vergleichbar mit den anderen FOC-Satelliten. Wenn die Datennutzungsrate und Residuen betrachtet werden, weist die NL FCB Qualität von FOC, IOV (Außer E19), und den anderen Statelliten mit hohem Exzentrizität keine wesentliche Differenz auf. Die PPP Mehrdeutigkeitauflösung zeigt, dass die Galileo und BDS Beobachtungen einen deutlichen Vorteil für GPS-only PPP Mehrdeutigkeitauflösung darstellen können.

In dieser Dissertation wird die PPP-Mehrdeutigkeitsauflösung auch auf Mehrfrequenzbeobachtungen erweitert. Um die Mehrfrequenzbeobachtungen richtig zu integrieren, wird basierend auf Rohbeobachtungen ein einheitliches unkombiniertes PPP-Mehrdeutigkeitsauflösungsmodell entwickelt. Mit solchem einheitlichen Modell können die von Mehrfrequenz-Beobachtungen erzeugten FCBs flexibel verwendet werden, beispielsweise für zwei- oder drei Frequenzen PPP Mehrdeutigkeitauflösung. Die Effizienz dieses Modells wird mit Galileo und Beidou drei-Frequenz-Beobachtungen geprüft, die von globalen verteilten MGEX Stationen gesammelt werden. Die geschätzten FCBs werden hinsichtlich der Residuenausverteilung und Standardabweichungen bewertet, was eine gute Übereinstimmung mit der float-Lösung der Mehrdeutigkeit zeigt. Die Standardabweichung von linear kombinierten FCBs ist viel geringer als die von Roh-FCBs. Diese bedeutet, dass ein linear kombinierte FCBs effizienter für das Broadcasting in Echtzeitanwendungen sind. Die drei-Frequenz-PPP-Mehrdeutigkeitauflösung zeigt, dass die Positionsabweichung im Vergleich zur float-Lösung stark reduziert sind. Die Verbesserungen für die Ost, Nord und Höhenkompenenten mit BDS Positionierung sind 49.2%, 38.3% und 29.6%, in der Reheinfolge. Die entsprechenden Verbesserungen mit Galileo Positionierung sind 60.0%, 29.0% und 21.1%.

Neue Herausforderung für Cycle-slip Erkennung entsteht durch die Verarbeitung des unkombinierten Signals. Da alle Tragerfrequenzbeobachtungen separat verarbeitet werden, ist es sehr wichtig, die Trägerfrequenz des Cycle slips zu identifizieren, damit die anderen Beobachtungen nicht davon beeinflusst werden können. Um kontinuierliche Trägerphasenmessungen für die vorher genannten Forschungen bereitzustellen, wird ein verbesserter Ansatz entwickelt, der auf einem zeitlich differenzierten Modell für die Erkennung und Reparatur des Cycle-slips basiert. Dieses Modell reduziert die Gefahr für falsche Entscheidung und erhöht die Chance auf eine erfolgreiche Abschätzung des Cycle-slips.

Zusammenfassend ist festzuhalten, dass die PPP mittles only-GPS auf Mehrfrequenz und Mehr-GNSS-PPP-Mehrdeutigkeitauflösung mit verbesserter Genauigkeit und schneller Konvergenz erweitert wurde. Dies wurde durch ein einheitlichtes Modell erreicht, das auf nicht kombinierte PPP basiert. Dieses Model wurde sorgfältig untersucht und mit verbesserter Erkennung und Reparatur von cycle-slip, schneller FCB-Schätzverfahren und Mehrdeutigkeitsauflösung angereichert.

Diese artikel-basierte (kumulative) Dissertation besteht aus einem ausführlichen Einführungskapitel und vier Kapiteln, die mit den folgenden Publikationen verbunden sind:


## **Acronyms**






## **Chapter 1 Introductory chapter**

## **1 Introduction**

### **1.1 Motivation and background**

Satellite-based positioning has become an integral part of our modern-day society and is used by people all over the world (Teunissen and Montenbruck 2017). For positioning with satellite systems, pseudorange and carrier phase are the basic observables. Due to the low precision of pseudorange measurements and broadcast ephemeris, the resulting accuracy of standard point positioning is at several meters level (Hofmann-Wellenhof et al. 2007). For high-accuracy positioning at centimeters level, the carrier phase measurements must be adopted. The carrier phase-based positioning techniques can be roughly divided as differential positioning, such as real time kinematic (RTK), and absolute positioning, i.e., precise point positioning (PPP). Differential positioning requires reference receivers to mitigate the common errors, while PPP is characterized by single receiver with sophisticated error models and precise satellite information. PPP is able to achieve centimeter to sub centimeters accuracies with dual-frequency carrier phase and pseudorange measurements. Due to its cost-effectiveness, global coverage, and high accuracy, PPP has found increased applications. However, PPP shows a few drawbacks (Bisnath and Gao 2008), e.g., slow convergence and lower accuracy compared to differential positioning. The integer ambiguity resolution (AR) technique is expected to further enhance the accuracy and shorten the convergence time. Although the theory was proposed and elaborated in (Gabor and Nerem 1999; Gao and Shen 2002), PPP AR was not practically realized and implemented until (Collins et al. 2008; Ge et al. 2008; Laurichesse et al. 2009). It is proved that AR is able to improve the PPP accuracy, especially the east component, to a comparable level as that of differential positioning. Many applications including ionospheric modelling (Banville and Langley 2011), tropospheric modelling (Li et al. 2015c), time transfer (Petit et al. 2015), seismic displacement monitoring (Li et al. 2013c), low earth satellite orbit (Montenbruck et al. 2018) and navigation satellite orbit determination (Katsigianni et al. 2018), have further confirmed its efficiency. After ten years development, the technique is so impressive that the global navigation satellite system (GNSS) community proposed to establish a specific working group for PPP AR at the IGS Workshop 2018.

While remarkable progress has been made with legacy GPS observations, the ongoing modernization and the build-up of new GNSS constellations (shown in Table 1) offers exciting prospects for further improvements. The GNSS has evolved from GPS to four major global systems including the Russian GLONASS, Chinese BDS, and European Galileo. Specifically, GLONASS has restored its full orbital constellation of 24 satellites since 2011 (Revnivykh 2012). For BDS, the construction is planned for three steps (Yang et al. 2018). The formal establishment of BDS-1 was marked by the launch of the third BeiDou satellite in 2003. The second step is the BDS-2 regional system, which was accomplished by a constellation of 14 satellites in 2012. With the successful launch of BDS third generation experimental satellites in 2015, the construction of the BDS global system (BDS-3) starts to speed up, and is expected to be completed in 2020. Similarly, Galileo also witnesses remarkable achievements in recent years (Steigenberger and Montenbruck 2017). On 15 December 2016, Galileo started offering initial services. As of March 2019, there are 22 satellites in orbit and 4 satellites under commissioning. The complete constellation with 30 satellite could be expected by 2020.

As can be seen from Table 1, there are 110 GNSS satellites in orbit with 10 satellites under commissioning at present. When all four global systems reach their full constellations, there will be more than 120 usable satellites. Compared with 32 GPS satellites, the fusion of multiple GNSS can significantly increase the number of visible satellites by a factor of nearly 4 and improve the spatial geometry. Fast convergence, improved accuracy and enhanced reliability can be expected. Under this circumstance, the contribution of multi-GNSS integration to GPS standalone PPP, as well as ambiguity resolution, should be investigated.


**Table 1** Status of multi-GNSS as of Mar. 2019

In recent years, PPP has found new opportunity in multi-frequency observations. The additional signals can be used to form combinations with excellent properties, such as low noise and reduced ionospheric delay (Feng 2008). The additional signals are expected to further enhance the performance of PPP (Geng and Bock 2013). As of March 2019, triple- or multi-frequency signals have become available from more and more GNSS satellites, as shown in Table 2. In the effort to modernize GPS, new L5 signals were added, and are available now on all satellites from GPS IIF (Tegedor and Øvstedal 2014). The L5 signals provide secure and robust transmissions for life critical applications. However, there exists additional inter-frequency clock bias (Montenbruck et al. 2012), which should be accounted for (Li et al. 2013a; Pan et al. 2018). The legacy GLONASS satellites transmit on a different frequency using a 15-channel frequency division multiple access (FDMA) technique (Glonass 2008). The 24-satellite constellation is accommodated with only 15 channels, as identical frequency channels can be used for antipodal satellite pairs which are never both in view at the same time for earth-based users. The FDMA technique introduces the receiver inter-frequency bias (IFB) between satellites, which is a problem for high-precision positioning (Reußner and Wanninger 2012; Wanninger 2012). To avoid the deficiency, new code division multiple access (CDMA) signals were researched for GLONASS since 2008. The GLONASS-K1 satellites launched in 2011 introduced L3 CDMA signals, followed by the enhanced GLONASS-K1 and GLONASS-K2 satellites in 2018 which introduced L1, L2 and L3 CDMA signals (IAC 2016). These CDMA signals may replace the FDMA signals in the future. Galileo provides an open signal in the E1 band and a wideband signal covering the E5 a&b band (Zaminpardaz and Teunissen 2017). The Alternating Binary Offset Carrier (AltBOC) modulation can either be tracked as a composite signal or as distinct signals in the E5a and E5b sub-bands. For BDS-1 and BDS-2, three signals, e.g., B1, B2, and B3, were broadcasted, while three new carrier frequencies (B1C, B2a, Bs) were added in BDS-3 (Yang et al. 2018).


**Table 2** Status of carrier frequency of GNSS as of Mar. 2019 [Unit: MHz]

The multi-frequency signals are anticipated to enable rapid ambiguity resolution in differential positioning. Results indicate that instantaneous triple-frequency ambiguity resolution can be achieved with optimal combinations (Jung et al. 2000; Vollath et al. 1999). The underlying foundation is that combinations with reduced ionospheric delays, low noise and longer wavelength can be formed from multi-frequency observations (Geng and Bock 2013). Considering the proved efficiency in differential positioning, multi-frequency observations are expected to improve the performance of PPP.

In conclusion, multi-frequency and multi-constellation GNSS observations show a clear advantage over standard dual-frequency GPS observations. For PPP, improved accuracy and fast convergence time can be expected. However, proper integration of measurements from multi-constellation and multi-frequency GNSS into PPP requires detail investigations, which motivate this research.

#### **1.2 Outline**

The outline of the rest of this introductory chapter is presented as follows: Section 2 introduces the traditional dual-frequency PPP and ambiguity resolution models. Section 3 highlights the limitations of the current model, followed by the research objectives of this thesis. While Section 4 summarizes the main contribution of this thesis, Section 5 provides an overview of the related publications as presented in Chapters 2-5. This introductory chapter closes with Section 6, in which an overall conclusion of the thesis including an outlook for the future research is presented.

## **2 Foundations**

For a satellite s observed by receiver r, the corresponding raw pseudorange and carrier phase observation equations can be expressed as (Hofmann-Wellenhof et al. 2007; Leick et al. 2015; Xu 2007)

$$\begin{cases} P\_{r,f}^{s} = \rho\_r^s + dt\_r - dt^s + dT + a\_f \cdot dl\_{r,1}^s + D\_{r,f} - D\_f^s + \varepsilon\_{P\_f} \\ \Phi\_{r,f}^{s} = \rho\_r^s + dt\_r - dt^s + dT - a\_f \cdot dl\_{r,1}^s + \lambda\_f \{N\_{r,f}^s + B\_{r,f} - B\_f^s\} + \varepsilon\_{\Phi\_f} \end{cases} \tag{1}$$

where the subscript = (1,2,3, ⋯ ) refers to a specific carrier frequency, superscript s refers to a specific satellite; indicates the geometric distance between the satellite and receiver; and are the clock errors of receiver and satellite; is the slant tropospheric delay; ,1 is the slant ionospheric delay on the first carrier frequency and = 2 1 <sup>2</sup> ⁄ is the carrier frequency-dependent factor; , and are the receiver and satellite specific code hardware delays; and , are the wavelength in meter and integer ambiguity in cycle; , and are the receiver-dependent and satellite-dependent uncalibrated phase delays; and are the pseudorange and carrier phase measurement noise, respectively. Note that the higher-order ionospheric effects are neglected, as they have limited influence on the performance of ambiguity resolution (Hadas et al. 2017).

In the classic GPS dual-frequency PPP, the ionospheric-free (IF) combination is routinely employed to eliminate the effect of the first-order ionospheric delay(Zumberge et al. 1997). After multiplied by the IF coefficients, the corresponding pseudorange and carrier phase observation equations can be expressed as

$$\begin{cases} P\_{r,IF}^{s} = \rho\_r^s + dt\_r - dt^s + dT + D\_{r,IF} - D\_{IF}^s + \varepsilon\_{P\_{IF}}\\ \Phi\_{r,IF}^{s} = \rho\_r^s + dt\_r - dt^s + dT + \lambda\_{IF} \left( N\_{r,IF}^{s} + B\_{r,IF} - B\_{IF}^s \right) + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{2}$$

Conventionally, precise orbit and clock products from the international GNSS service (IGS) analysis center are used to remove satellite orbit and clock errors. The pseudorange ionospheric-free hardware delay bias is assimilated into the clock offset in accordance with the IGS analysis convention. Due to the fact that pseudorange measurements provide the reference to clock parameters, the actual receiver clock estimate would absorb the ionospheric-free combination of the receiver pseudorange hardware delay ,. After applying the GNSS precise satellite clock products, equation (2) can be rewritten as

$$\begin{cases} P\_{r,IF}^{s} = \rho\_r^s + dT + \overline{d}t\_r + \varepsilon\_{P\_{IF}} \\ \Phi\_{r,IF}^{s} = \rho\_r^s + dT + \overline{d}t\_r + \lambda\_{IF} \overline{N}\_{r,IF}^s + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{3}$$

where the receiver clock and ambiguity can be re-parameterized as

$$\overline{dt}\_{\mathbf{r}} = \left(dt\_{\mathbf{r}} + D\_{\mathbf{r}, \mathrm{IF}}\right) \tag{4}$$

$$
\overline{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + b\_{r,IF} - b\_{IF}^{s} \tag{5}
$$

$$b\_{r,IF} = B\_{r,IF} - \frac{D\_{r,IF}}{\lambda\_{IF}} \tag{6}$$

$$b\_{IF}^{\rm s} = B\_{IF}^{\rm s} - \frac{D\_{IF}^{\rm s}}{\lambda\_{IF}} \tag{7}$$

With at least four satellites simultaneously observed, the equations can be integrated and resolved with least squares algorithm. The estimable parameters are receiver coordinates, tropospheric delay, receiver clock offset and ambiguities. Besides this estimation strategy, the state-of-the-art corrections (Kouba 2009; Zumberge et al. 1997) including the satellite and receiver antenna phase center offset and variation (Schmid et al. 2005), the relativistic clock effects (Kouba 2004), the phase wind-up (Wu et al. 1991), the site displacement effects due to solid earth tide and ocean loading (Petit and Luzum 2010), and the differential code biases (Dach et al. 2015) should be considered. Fortunately, these biases are well researched and can be corrected with models or external products.

The above discussed model is referred to as standard dual-frequency PPP model. Note that the estimated ambiguity parameter is a combination of the integer ambiguity, the corresponding code hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite ends. In this case, the integer property is lost. The fact that doubledifferenced ambiguities in global or regional networks can be resolved to integer values (Dong and Bock 1989) lays the foundation of integer ambiguity resolution for PPP. The resolved double-differenced ambiguity implies that the fractional parts of two singledifferenced ambiguities (across satellites) must agree well with each other (Ge et al. 2008). By estimating the fractional parts, denoted as fractional cycle bias (FCB) in the sequel, at the server end and applying them to single differencing PPP at the user end, PPP integer ambiguity resolution can be achieved.

For PPP ambiguity resolution, the term , is usually decomposed into the following combination

$$\overline{N}\_{r,IF}^{s} = \frac{\frac{cf\_2}{f\_1^2 - f\_2^2} N\_{r,WL}^{s} + \frac{c}{f\_1 + f\_2} \overline{N}\_{r,NL}^{s}}{\lambda\_{IF}} \tag{8}$$

where , is the integer wide-lane (WL) ambiguity and , is the float narrow-lane (NL) ambiguity. Usually, the WL ambiguity is resolved by the Hatch-Melbourne-Wübbena (HMW) combination observable (Hatch 1982; Melbourne 1985; Wübbena 1985). With the fixed WL ambiguity, the float NL ambiguity can be derived based on (8), and tested whether it is also fixable. An ionospheric-free ambiguity is fixed only when both its WL and NL ambiguities are fixed. Now, the fixing of IF ambiguity has been transformed to fixing of WL and NL ambiguities. Similarly, the estimation of IF FCB has been transformed to estimation of WL and NL FCBs. The method of WL and NL FCB determinations is presented in the following.

The float WL ambiguity can be derived from

$$HMW = \frac{f\_1 \Phi\_{r,1}^s - f\_2 \Phi\_{r,2}^s}{f\_1 - f\_2} - \frac{f\_1 P\_{r,1}^s + f\_2 P\_{r,2}^s}{f\_1 + f\_2} = \lambda\_{r,WL}^s \overline{N}\_{r,WL}^s \tag{9}$$

where , = 1−<sup>2</sup> is the wavelength of WL ambiguity. The float WL ambiguity can be further decomposed as

$$
\overline{N}\_{r,WL}^{s} = N\_{r,WL}^{s} + d\_{r,WL} - d\_{WL}^{s} \tag{10}
$$

$$d\_{r,WL} = B\_{r,1} - B\_{r,2} - \frac{f\_1 D\_{r,1} + f\_2 D\_{r,2}}{\lambda\_{r,WL}^s (f\_1 + f\_2)} \tag{11}$$

$$d\_{WL}^{\rm s} = B\_1^{\rm s} - B\_2^{\rm s} - \frac{f\_1 D\_1^{\rm s} + f\_2 D\_2^{\rm s}}{\lambda\_{r,WL}^{\rm s} (f\_1 + f\_2)} \tag{12}$$

HMW combinations eliminate the geometric distance, satellite and receiver clock errors, as well as atmospheric delays. Averaging the HMW for a continuous arc results in accurate estimations of WL ambiguities, which can be easily fixed to integers by rounding.

When the WL ambiguity , is correctly resolved to an integer value and the float IF ambiguity , is obtained from PPP solution, the float NL ambiguity observable , can be derived based on (8)

$$
\overline{N}\_{r,NL}^{s} = \frac{f\_1 + f\_2}{c} \lambda\_{IF} \overline{N}\_{r,IF}^{s} - \frac{f\_2}{f\_1 - f\_2} N\_{r,WL}^{s} = N\_{r,1}^{s} + d\_{r,NL} - d\_{NL}^{s} \tag{13}
$$

where

$$d\_{r,NL} = \frac{f\_1 + f\_2}{\mathfrak{c}} \left(\lambda\_{IF} B\_{r,IF} - D\_{r,IF}\right) \tag{14}$$

$$d\_{NL}^{\rm s} = \frac{f\_1 + f\_2}{c} (\lambda\_{IF} B\_{IF}^{\rm s} - D\_{IF}^{\rm s}) \tag{15}$$

Equations (10) and (13) serve as the basic model for estimating FCBs. Since they have the same structure, a general expression can be formulated as

$$R^s\_r = \overline{N}^s\_r - \mathcal{N}^s\_r = d\_r - d^s \tag{16}$$

for WL and NL linear combinations, respectively. represents the FCB measurements; denotes the float undifferenced ambiguities; ̂ denotes the integer part of , which is the sum of the original integer ambiguity and the integer part of the combined code hardware delays and uncalibrated phase delays from both receiver *r* and satellite *s*; and denote the receiver and satellite FCBs. Note that ̂ can be determined by rounding ̅ assuming the float ambiguities ̅ are precisely estimated.

A set of equations in the form of (16) can be generated based on a network of reference stations. Ge et al. (2008) proposed to resolve the system by averaging based on all involved float single-differenced WL and NL ambiguity estimates. The single-differencing across satellites within single station eliminates the receiver FCB and single-differenced satellite FCB can be obtained. Since the same single-differencing approach will be implemented at the user end when using the FCB product, it is not necessary to obtain the absolute FCB and the single-differenced FCB would be sufficient. Based on the same principle but instead of the averaging process, a least squares method in an integrated adjustment was adopted to enhance the estimates (Li and Zhang 2012; Zhang and Li 2013). Suppose that there are *n* satellites tracked at *m* reference stations, the system of equations can be expressed as

$$
\begin{bmatrix} R\_{r,1}^{s,1} \\ \vdots \\ R\_{r,1}^{s,n} \\ \vdots \\ R\_{r,m}^{s,1} \\ \vdots \\ R\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} \overline{N}\_{r,1}^{s,1} - \delta \underline{J}\_{r,1}^{s,1} \\ \vdots \\ \overline{N}\_{r,1}^{s,n} - \delta \underline{J}\_{r,1}^{s,n} \\ \vdots \\ \overline{N}\_{r,m}^{s,1} - \delta \underline{J}\_{r,m}^{s,1} \\ \vdots \\ \overline{N}\_{r,m}^{s,n} - \delta \underline{J}\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} 1 & \cdots & 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cdots & 0 & 0 & \cdots & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & \cdots & -1 \end{bmatrix} \begin{bmatrix} d\_{r,1} \\ \vdots \\ d\_{r,m} \\ d\_{r,1}^{s,n} \\ \vdots \\ d\_{r,n}^{s,n} \end{bmatrix} \tag{17}$$

The obtained system of equations is singular. One arbitrarily selected FCB should be set to zero. In this way, the system of equations can be resolved with least squares algorithm (Li et al. 2015a). Note that the above discussions are based on single epoch data. Considering the temporal stabilities of code hardware delays and uncalibrated carrier phase delays, the satellite FCBs can be estimated as constant values for specific time intervals, e.g., one day for WL FCB and 15 minutes for NL FCB (Ge et al. 2008). This can be accomplished by accumulating all the equations within the specific time intervals. In this way, the final FCB products are produced and broadcasted to users. Since January 1, 2015, the School of Geodesy and Geomatics at Wuhan University (WHU-SGG) routinely releases GPS WL and NL FCB products with open access (Li et al. 2015a).

With the above-estimated FCBs, users are able to correct the satellite FCBs, and the single-differenced PPP ambiguities are aimed to be fixed

$$N\_r^{s,m} = \overline{N}\_r^s - \overline{N}\_r^m - d^{s,m} \tag{18}$$

For WL float ambiguities, they can be directly fixed by the rounding approach (Ge et al. 2008). For NL float ambiguities, they are fed into the Least-squares AMBiguity Decorrelation Adjustment (LAMBDA) algorithm (Teunissen et al. 1997). If not all the float NL ambiguities can be fixed by the LAMBDA method, partial ambiguity resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999). Once the WL and NL integer ambiguities are resolved and validated, a tight constraint can be reconstructed and imposed on the ionospheric-free float ambiguities

$$0 = \overline{N}\_r^s - \overline{N}\_r^m - N\_r^{s,m} - d^{s,m},\ \sigma\_d \tag{19}$$

In this way, a fixed solution can be obtained at the user end.

The above discussed model is referred as FCB based PPP AR model. Collins et al. (2008) developed a decoupled clock model by separating satellite clocks for code and phase observations. Similarly, Laurichesse et al. (2009) developed an integer phase clock model in which the NL FCBs were assimilated into receiver and satellite clock estimates of a network solution. This model has been employed to generate the precise satellite clock products by Groupe de Recherche de Géodésie Spatiale of the Centre National d'Etudes Spatiales (CNES-GRG) (Loyer et al. 2012). The only difference between the decoupled clock model and the integer phase clock model is the approach for determining the WL ambiguity. The integer phase clock model utilizes WL FCB corrections and a satelliteaveraging process to fix the integer WL ambiguity, whereas the decoupled clock model directly estimates the integer WL ambiguity along with other unknowns through the functional model. The difference between the WL/NL FCB model and the integer phase clock model is the strategy of separating NL FCBs from integer ambiguities. NL FCBs are directly estimated in the WL/NL FCB model, whereas they are assimilated into the clock estimates in the integer phase clock model. These PPP AR techniques are proven to be equivalent in theory (Shi and Gao 2014). The systematic biases between position estimates have been demonstrated to be minimal and negligible with data from a global network of reference stations (Geng et al. 2010). The FCB based model can conveniently supplement current network solutions as an additional software module, as the FCB determination is compatible with current official clock-generation methods. For this reason, the FCB based model is adopted throughout the thesis.

## **3 Problem statements and research objectives**

The PPP AR model and method described above have been proposed for decades. Due to the rapid development, and huge benefit, of multi-frequency and multi-GNSS measurements, it is required to extend the current model and method to integrate these measurements for fast and high precision positioning.

#### **3.1 Problem statements**

#### (1) Fast FCB estimation

Based on the discussion in Section 2, we find that the least squares method is routinely utilized for FCB estimation (Li et al. 2015a; Li et al. 2017b; Liu et al. 2017; Wang et al. 2017). However, the least squares method can be extremely time-consuming concerning the large number of observations from hundreds of reference stations and thousands of epochs during the FCB estimation. The dimension of the involved matrix could be further enlarged by multi-frequency and multi-GNSS observations, which also requires a lot of computation resource. This is unfavorable for the more and more popular realtime applications (Geng et al. 2018; Li et al. 2015b; Li et al. 2013c). In addition, iterations are required to deal with the one cycle inconsistency among FCB measurements. Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU-SGG and CNES-GRG products respectively, the one cycle inconsistency arises whenever the superposition of receiver FCB and satellite FCB exceeds the boundary (Xiao et al. 2018b). The additional iterations of the least squares method and computation of a large matrix demand a long time to finish. Therefore, a fast and efficient estimation method is desirable.

#### (2) Extend PPP AR to multi-GNSS

Compared with GPS, the study of other GNSS PPP AR is limited and requires more attentions considering the rapid developments in recent years. For GLONASS, special care has to be taken to deal with the receiver inter-frequency bias between satellites originating from the frequency division multiple access strategies. Usually, it is achieved by using a network of homogeneous receivers (Geng and Bock 2016; Geng and Shi 2016; Liu et al. 2017). Also, BDS-2 code measurements suffer from the satellite induced pseudorange variations which are elevation- and frequency-dependent. Correction models based on multipath combination have been proven effective (Wanninger and Beer 2015). Unfortunately, this approach does not apply to geostationary orbit satellites due to the almost constant satellite elevation angle. Currently, only regional PPP AR or PPP AR with inclined geosynchronous orbit and medium earth orbit satellites is achievable for BDS (Gu et al. 2015a; Li et al. 2017a; Li et al. 2017b). Tegedor et al. (2016) initially assessed the Galileo PPP AR with four In-Orbit Validation satellites using regional stations around Europe. Li et al. (2017b) estimate the FCBs for the four global systems and assess the performance of combined PPP AR. The studies of Galileo and BDS PPP AR have been limited by the number of available satellites. The situation has been changed, with a lot of new Galileo and BDS satellites in the last two years (Xiao et al. 2019b). It is of interest to extend PPP AR to multi-GNSS, especially with the relatively new Galileo and BDS.

#### (3) Extend PPP AR to multi-frequency

The traditional PPP based on dual-frequency ionospheric-free combination cannot properly handle the additional signals, especially the functional and integer ambiguity resolution model. As a result, the PPP model based on uncombined measurements (Schönemann et al. 2011; Zhang et al. 2012), in which the individual signal of each carrier frequency is treated as an independent observable, has drawn increasing interest in the GNSS community (Gu et al. 2015b). Its efficiency has already been confirmed in terms of convergence time and precision for single-frequency PPP (Lou et al. 2016), multi-GNSS PPP (Chen et al. 2015), as well as PPP-RTK (Feng et al. 2013; Odijk et al. 2016). Moreover, this approach has been tested effectively for ionospheric modelling (Tu et al. 2013; Zhang et al. 2012), differential code bias (DCB) estimation (Liu et al. 2018; Shi et al. 2015), and satellite orbit determination (Strasser et al. 2018; Zehentner and Mayer-Gürr 2016). Compared to the well-developed IF PPP model, the uncombined PPP model, called uncombined PPP in the sequel, requires more investigations, especially in the case of multifrequency processing and ambiguity resolution. First, how to model and constrain the ionospheric delay has a crucial impact on performance using the uncombined PPP. For example, it has been demonstrated that a white noise model is not adequate to capture the characteristics of the ionospheric delay. The external constraints developed from the ionospheric products, such as the IGS global ionosphere maps, are also not accurate enough to completely separate the ionospheric effects from the ambiguity parameters (Gu et al. 2015b). The influence of the ionospheric effects on the ambiguity fixing therefore must be reduced. Second, the method to deal with the DCB errors is more problematic with the uncombined PPP (Guo et al. 2015) than the dual-frequency IF PPP, since the latter can cancel out the DCB biases (Kouba and Héroux 2001). The problem of partial assimilation of the code bias (DCB) into phase bias (FCB) should also be carefully considered. Third, the uncombined PPP approach was proposed to deal with multi-GNSS and multi-frequency signals, so a generalized FCB estimation and AR method (Li et al. 2018), which is extendable to dual-, triple-, and multi- frequency, should be proposed.

Li et al. (2013b) verified the feasibility of the uncombined PPP AR with refined ionospheric models. The ionospheric delay was constrained from a priori spatial-temporal information and ionospheric products. The GPS dual-frequency ambiguities were fixed sequentially in the forms of WL/NL, which followed the convention of IF PPP AR. Gu et al. (2015a) further testified the uncombined PPP AR with BDS triple-frequency observations. The extra-wide-lane (EWL) and WL ambiguities were successfully fixed, whereas the B1 ambiguities were kept as float values. In addition, the performance was further limited by the satellite-induced multipath effects (Wanninger and Beer 2015). Li et al. (2018) proposed a unified FCB estimation and PPP AR method, which is extendable to multi-frequency uncombined PPP. The FCBs on each frequency were directly estimated from the raw float ambiguities derived from triple frequency observables. The model showed a great potential for multi-frequency uncombined PPP AR, although its DCB strategy may not be optimal (Xiao et al. 2019a). The satellite DCBs, together with the receiver DCB, were estimated as unknowns, and as a result the number of unknown parameters was increased. Given that the satellite GNSS DCB product is currently available on a routine basis (Wang et al. 2016), it would be beneficial to make use of these products.

In conclusion, it is required to develop a unified modeling strategy for multi-frequency uncombined PPP ambiguity resolution.

#### (4) Data preprocessing in uncombined PPP

The uncombined signal processing brings new challenges for cycle slip detection. The most significant feature is carrier frequency identification of cycle slips. Since all carrier frequency observations are processed separately, it is essential to identify the carrier frequency of the cycle slips to avoid contaminating other observations (Xiao et al. 2018a). In conventional methods, such as TurboEdit algorithm (Cai et al. 2013; Hatch 1982; Liu 2011; Melbourne 1985; Wübbena 1985) and combinations of multi-frequency observations (Dai et al. 2009; de Lacy et al. 2012; Wu et al. 2010a; Zhao et al. 2015), a common characteristic is the forming of optimal combinations to mitigate the presence of geometric and ionospheric errors, which makes it impossible to identify the carrier frequency affected by a cycle slip. All carrier frequencies are flagged as containing cycle slips even if there is only one carrier frequency suffering from a cycle slip. Therefore, a cycle slip detection method, which is capable of carrier frequency identification of cycle slip and applicable to all GNSS should be proposed.

## **3.2 Research objectives**

The general research objective of this thesis is to improve the performance of GPS-only PPP AR by multi-GNSS integration and multi-frequency observations. In detail, the following questions and aspects have been addressed within the thesis:


## **4 Contributions**

### **4.1 Fast FCB estimation for PPP AR**

Estimating satellite FCBs based on the Kalman filter is proposed and demonstrated. Since the Kalman filter is based upon least squares methods (LSM), it is theoretically possible to calculate the same solution as for the commonly used LSM. In the proposed Kalman filter, the large number of observations is handled epoch-by-epoch, which significantly reduces the dimension of the involved matrix and accelerates the computation. Hence it outperforms the commonly used LSM in terms of efficiency. In order to avoid iterations caused by the one cycle inconsistency among FCB measurements, a pre-elimination method is developed based on the temporal stability of satellite FCBs. The preelimination method shows a clear advantage over post-residual adjustment, which further improves the efficiency (Xiao et al. 2018b).

(1) Pre-elimination of one cycle inconsistency

Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU-SGG and CNES-GRGS products respectively, the one cycle inconsistency arises whenever the superposition of receiver FCB and satellite FCB exceeds the boundary. A simple example is presented to depict the situation. Assuming there is one satellite ( = 0.2) tracked at two stations (,1 = 0.6 and ,2 = 0.8), the superposition of FCBs should be (,1 = −0.4 and ,2 = −0.6) . However, the ,2 would be 0.4 due to the rounding approach. As a consequence, the corresponding satellite FCB derived from the two stations differs by one cycle. The one cycle inconsistency can be detected by examining the posterior residuals, as employed in former works. The posterior adjustments are inefficient as iterations of the whole process are required. In contrast, we propose to eliminate the inconsistency in advance (shown in Figure 1), which shows a clear advantage over posterior adjustment. The proposed method consists of the following steps:


In addition, stations can also be handled individually to obtain clean FCB measurements. After the pre-elimination of the one cycle inconsistency, the clean FCB measurements are fed into the Kalman filter.

(2) Design of Kalman filter for FCB estimation

The Kalman filter addresses the general problem of state estimates of discrete time controlled processes that are governed by a linear stochastic difference equation (Kalman 1960). The theory has been well studied and widely applied (Yang 2006; Yang et al. 2010). Since the Kalman filter is based upon the theory of least squares, it is theoretically possible to calculate the same solution as LSM and versa vice. However, the design and normal equation matrix will be huge in LSM, considering the large number of observations from hundreds of reference stations and thousands of epochs. Computation of the large matrix is time-consuming, which is inefficient and unfavorable for real-time applications. In contrast, the large number of observations are handled epoch-by-epoch in a Kalman filter, which significantly reduces the dimension of the involved matrix and accelerates the computation. One additional advantage of Kalman filter is its real-time capability. This is of particular interest for real-time PPP AR and its applications. Therefore, a Kalman filter is employed in this work. The design and flowchart of a Kalman filter are depicted in Figure 2. In our case, is the input FCB measurement for the Kalman filter, while and are the output estimates.

**Figure 1** Pre-elimination of the one cycle inconsistency for station IQQE on 2016/02/26. The top panels represent receiver WL FCB adjustment, while the bottom ones represent that of receiver NL FCB. The numbers denote GPS satellite PRNs. The left and right panels represent the raw and adjusted measurements, respectively.

**Figure 2** Flow chart of the proposed Kalman filter. The dashed blocks represent additional steps adopted in this thesis.

Two additional steps are added to the standard Kalman filter (Xiao et al. 2018b). The first step aims at establishing the dynamic model and determining the state transition matrix by analyzing the temporal stability of existing FCB products. Satellite WL FCBs are stable over several days, which can be characterized as a constant parameter on a daily basis. Satellite NL FCBs are considered as constant within 15 minutes but exhibit small variations over 15-minute intervals. Therefore, NL FCBs are characterized as random walk processes. The second step is introduced to eliminate the one cycle inconsistency of FCB measurements. Only clean measurements are sent to the Kalman filter to avoid iterations. Note that a constraint is imposed on the Kalman filter to eliminate the rank deficiency. A satellite FCB is selected and set to zero, which is accomplished by adding a pseudo observational error equation with zero variance (Yang et al. 2010).

#### **4.2 Multi-GNSS PPP AR**

In the context of multi-GNSS, the PPP ambiguity resolution is extended to Galileo and BDS. In the view of the rapid developments of these systems in recent years, the feasibility of Galileo and BDS PPP ambiguity resolution with the current constellation is demonstrated. The satellite WL/NL FCBs are estimated from globally distributed MultiGNSS EXperiment (MGEX) stations and assessed by a comparison with those of GPS (Xiao et al. 2019b).

#### (1) WL FCB comparison

Results of 60 days indicate that the quality of Galileo WL FCB is better than for GPS and BDS in terms of data usage rate, residual distribution, as well as standard deviation (STD) of daily estimates (shown in Figure 3). We attribute the good quality of Galileo WL FCB to its advanced signal modulation, AltBOC, which significantly compresses the multipath effect for code measurements as shown Figure 4.

**Figure 3** STDs of daily satellite WL FCBs with 60 days. Averaged STDs are calculated for GPS and BDS, while an individual STD is calculated for each Galileo satellite. E30 is selected as the reference satellite for Galileo, while G01 and C09 are selected for GPS and BDS, separately.

**Figure 4** HMW linear combinations of GPS, Galileo, and BDS dual frequency code measurements at station YEL2 on DoY 160, 2017.

Within the Galileo constellation, the quality of Full Operational Capability (FOC) satellite WL FCB is much better than for In-Orbit Validation (IOV) satellites. The root mean square (RMS) for each satellite in Galileo are presented in Figure 5. It is found that the RMS of IOV satellites are significantly larger than other FOC satellites, while there is no significant difference between FOC and two highly eccentric (ECC) satellites. These results indicate that the quality of FOC WL float ambiguities is better than that of IOV satellites. The worse performance of IOV satellites probably stems from the reduction of signal transmission power for IOV satellites. ESA imposed a reduction of 1.5 dB in the signal power of all four IOV satellites following a payload power problem of the fourth IOV (E20) in 2014. The 1.5 dB power decrease results in an increase of approximately 15- 20% of the thermal noise of the receiver, which roughly matches the observed increase in IOV WL FCB residuals.

**Figure 5** RMS of WL residuals after WL FCB estimations. Averaged RMSs are calculated for GPS and BDS, while an individual RMS is calculated for each Galileo satellite.

#### (2) NL FCB comparison

Figure 6 depicts the usage rate of float ambiguities for each Galileo satellite (average 94.0%), as well as the average usage rates for GPS (97.6%) and BDS (80.4%). The usage rate of GPS NL float ambiguities is highest while that of BDS is worst. Since the weights of carrier phase measurements are 10000 times larger than those of code measurements, the quality of the float ambiguities is dominated by the unmodeled errors in PPP. Therefore, the highest usage of GPS is reasonable since its precise product and models are currently the state of the art. All Galileo satellites except E19 show similar usage rates ranging from 91.8 to 96.9%, indicating a similar quality level of NL ambiguities. The usage rate of E19 is only 83.8%, which is much smaller than that for other Galileo satellites.

Inspired by Zaminpardaz and Teunissen (2017) who reported that for SEPTENTRIO receivers the carrier-to-noise density ratio for E19 lies below the value of the other two IOV satellites for elevations higher than 60 degrees, we further calculated the carrier-to-noise density ratio, data usage rates for three receiver groups, namely LEICA (37 stations), SEPTENTRIO (41 stations) and TRIMBLE (83 stations). From the results in Figure 7, it can be seen that the carrier-to-noise density ratio of E19 lies below the other satellites for all the receiver groups. In addition, the data usage rates of E19 are lower than for the other satellites for all the three groups and no clear relationship with receiver type is found. Furthermore, we notice that the STD of satellite laser ranging (SLR) residuals for E19 (0.130 m) is larger than for E11 (0.092 m) and E12 (0.086 m) when performing SLR validation (Guo et al. 2017). This may indicate a poor quality of precise satellite products for E19. Based on the above analysis, we suspect that the lower data usage rate of E19 is caused by the joint effect of poorer satellite products and larger carrier phase noise.

**Figure 6** Usage rates of float NL ambiguities for satellite NL FCB estimation. Averaged usage rates are calculated for GPS and BDS, while an individual usage rate is calculated for each Galileo satellite.

**Figure 7** Carrier-to-noise density ratio C/N0 of Galileo signals for three receiver groups, TRIMBLE (left), SEPTENTRIO (middle) and LEICA (right).

In conclusion, the quality of Galileo NL FCB is slightly worse than that of GPS but better than that of BDS. Within Galileo, the NL FCB quality of FOC and IOV (except E19), as well as the two eccentric satellites, shows no significant difference in terms of data usage rates and residuals. The reason for the worse performance of E19 is still not clear. On the one hand, it cannot, or at least not fully, be ascribed to the signal transmission power as the power difference between FOC and IOV E11/E12 is larger than the difference between E11/E12 and E19. On the other hand, the worse quality of E19 satellite orbit, indicated by SLR residuals, could also degrade NL FCB estimations. This issue remains an open question and deserves further investigation.

For static Galileo PPP AR with three-hour sessions, the positional biases can be reduced by 67, 45 and 22% for east, north and up components, respectively. In addition, PPP AR also improves the kinematic solutions by 59, 45 and 28% for east, north and vertical components, respectively.

#### **4.3 Multi-frequency PPP AR**

A unified model for multi-frequency PPP AR based on raw uncombined observations is proposed (Xiao et al. 2019a), which simplifies the concept of phase biases for AR. No assumption is made on the method used to determine FCB on the server end, which implies that the generated FCB from multi-frequency observations could be flexibly used, such as for dual- or triple-frequency ones. It is demonstrated that the model is extendable to dual- and triple-frequency observations.

#### (1) Multi-frequency uncombined PPP AR model

To deal with multi-frequency observations, the uncombined PPP model is adopted. In the uncombined PPP model, the ionospheric delay is directly estimated. Another important difference between the IF PPP and the uncombined PPP models is the strategy to deal with the DCB. The DCB is not of concern in IF PPP as the IF combination is also used for precise clock generation, which implies that the DCB could be fully absorbed by other parameters (Zhang et al. 2012) or simply cancelled out in the IF PPP (Dach et al. 2015). But this is not the case in the uncombined PPP, especially with multi-frequency observations. We propose estimation of the receiver DCB and correction of the satellite DCB with existing multi-GNSS DCB products (Wang et al. 2016). Taking triple-frequency observations as an example, the correction equation for the raw pseudorange observations can be deduced as (Guo et al. 2015)

$$\begin{cases} P\_{r,1}^{s} = \rho\_r^s + \bar{d}t\_r - dt\_{pre}^s + dT + a\_1 \cdot dl\_{r,1}^s - \frac{1}{a\_2 - 1} DCB\_{12}^s + \varepsilon\_{\mathcal{P}\_f} \\ P\_{r,2}^{s} = \rho\_r^s + \bar{d}t\_r - dt\_{pre}^s + dT + a\_2 \cdot dl\_{r,1}^s - \frac{a\_2}{a\_2 - 1} DCB\_{12}^s + DCB\_{r,12} + \varepsilon\_{\mathcal{P}\_f} \\ P\_{r,3}^{s} = \rho\_r^s + \bar{d}t\_r - dt\_{pre}^s + dT + a\_3 \cdot dl\_{r,1}^s - DCB\_{13}^s - \frac{1}{a\_2 - 1} DCB\_{12}^s + DCB\_{r,13} + \varepsilon\_{\mathcal{P}\_f} \end{cases} \tag{20}$$

where <sup>12</sup> = <sup>2</sup> − <sup>1</sup> , ,12 = ,2 − ,1 , and ̅ = + ,1 . The <sup>12</sup> and <sup>13</sup> can be obtained from multi-GNSS DCB products, while ,12 and ,13 are estimated as daily constant parameters. Similarly, the phase equations can be rewritten as

$$
\Phi\_{r,f}^{s} = \rho\_r^s + \vec{d}t\_r - dt\_{pre}^s + dT - a\_f \cdot dI\_{r,1}^s + \lambda\_f \overline{N}\_{r,f}^s + \varepsilon\_{\Phi\_f} \tag{21}
$$

where the ambiguity can be re-parameterized as

$$\begin{cases} \overline{N}\_{r,f}^{s} = N\_{r,f}^{s} + b\_{r,f} - b\_{f}^{s} \\ b\_{r,f} = B\_{r,f} - D\_{r,1}/\lambda\_{f} \\ \quad b\_{f}^{s} = B\_{f}^{s} - \frac{D\_{IF}^{s}}{\lambda\_{f}} \end{cases} \tag{22}$$

and the estimable parameters are

$$X = \begin{bmatrix} x & y & z & \overline{d}t\_r & dT & l\_{r,1}^s & D\_{r,12} & D\_{r,13} & \overline{N}\_{r,1}^s & \overline{N}\_{r,2}^s & \overline{N}\_{r,3}^s \end{bmatrix} \tag{23}$$

The estimated ambiguity parameter is a combination of the integer ambiguity, the corresponding code hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite ends. Similar to that of IF PPP, the integer property of the ambiguity parameter can be recovered provided FCB is corrected. In dual-frequency IF PPP, the float ambiguity is usually decomposed into WL/NL forms in order to recover the integer property (Ge et al. 2008). This is partly because the IF combination of L1/L2 ambiguities is, in essence, not an integer. Another reason is that the WL ambiguities possess a relatively longer wavelength and are less correlated, therefore can be easily fixed. For uncombined PPP AR, it is also important to form combinations of raw ambiguities. On the one hand, the estimated raw float ambiguities are strongly correlated (Li et al. 2018). On the other hand, the raw float ambiguities are quite sensitive to unmodeled ionospheric errors (Gu et al. 2015b). Therefore, the combinations with longer wavelengths and lower ionospheric delays are preferred. From the systematic investigation of triple-frequency combinations, it is found that the following combinations are a good compromise between ionospheric reduction and noise amplification. These combinations are used for both BDS and Galileo triple-frequency observations for simplicity

$$
\begin{bmatrix}
\overline{N}^s\_{r,LC1} \\
\overline{N}^s\_{r,LC2} \\
\overline{N}^s\_{r,LC3}
\end{bmatrix} = \begin{bmatrix}
4 & -3 & 0 \\
1 & -1 & 0 \\
1 & 0 & -1
\end{bmatrix} \begin{bmatrix}
\overline{N}^s\_{r,1} \\
\overline{N}^s\_{r,2} \\
\overline{N}^s\_{r,3}
\end{bmatrix} \tag{24}
$$

Substituting (22) into the above system produces the basic model for estimating FCBs. Now the equations are quite similar to (16), which implies that the same method can be used to estimate FCB for these combinations, i.e., (17).

With the obtained combined FCB, we are able to calculate the FCB of the raw L1/L2/L3 carrier frequency

$$
\begin{bmatrix} d\_1^s \\ d\_2^s \\ d\_3^s \end{bmatrix} = \begin{bmatrix} 4 & -3 & 0 \\ 1 & -1 & 0 \\ 1 & 0 & -1 \end{bmatrix}^{-1} \begin{bmatrix} d\_{LC1}^s \\ d\_{LC2}^s \\ d\_{LC3}^s \end{bmatrix} \tag{25}
$$

The transformation from combined FCBs to raw FCBs is important, as it provides more flexibility to the users. With the raw FCBs, users are able to choose their own linear combinations of observations, formulate the corresponding combined FCB, and conduct PPP AR. This representation allows interoperability if the server and user sides implement different AR methods. In addition, the raw FCB is suitable for the State Space Representation (SSR) of Radio Technical Commission for Maritime services (RTCM) (Weber et al. 2005), where one phase bias per phase observable is broadcasted instead of making specific combinations.

Similar to dual-frequency IF PPP AR, single differencing across satellites must be firstly performed in order to remove receiver FCBs. Then the single-differenced ambiguities from different carrier frequencies are combined, as has been done during FCB estimation

$$
\begin{bmatrix}
\overline{N}\_{r,ljk1}^{m,n} \\
\overline{N}\_{r,ljk2}^{m,n} \\
\overline{N}\_{r,ljk3}^{m,n}
\end{bmatrix} = \begin{bmatrix}
i\_1 & j\_1 & k\_1 \\
i\_2 & j\_2 & k\_2 \\
i\_3 & j\_3 & k\_3 \\
\end{bmatrix} \begin{bmatrix}
\overline{N}\_{r,1}^{m,n} \\
\overline{N}\_{r,2}^{m,n} \\
\overline{N}\_{r,3}^{m,n}
\end{bmatrix} \tag{26}
$$

where ̅ , = − is the single-differenced ambiguity between satellites m and n. Based on the coefficients (,, ), the FCB for the specific combined ambiguity can also be formed. Note that the linear combinations are not necessarily to be the same as those in FCB estimation, although the three combinations mentioned above are strongly recommended. In our experiments, we have used the same combinations as in FCB generation for uncombined PPP AR. Usually, the EWL/WL float ambiguities can be directly fixed by the rounding approach after the correction of FCB (Ge et al. 2008), and the NL float ambiguities are fed into the LAMBDA algorithm to search for correct integers (Teunissen et al. 1997). However, in our study, the LAMBDA is used for each combination, regardless of its property, which simplifies the design of the algorithm. In addition, if not all the float ambiguities can be fixed by the LAMBDA method, partial ambiguity resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999). It is found that the searching and fixing of ambiguities for the combination with longer wavelengths (e.g., EWL/WL) is quite fast. When the integer ambiguities for one combination are resolved and validated, a tight constraint can be reconstructed. The number of constraints accumulate as the process repeats for all linear combinations. Afterwards, the constraints are imposed on the raw ambiguities, and yields the AR solution. Note that the ambiguities in IF PPP AR must be sequentially fixed in the order of WL/NL. An IF ambiguity is constrained only when both its WL and NL ambiguities are fixed, while the linear combined ambiguities in our approach can be fixed and constrained independently.

#### (2) Test results with Galileo and BDS triple-frequency observations

To verify the efficiency of the proposed approach, we processed 51 days of Galileo and BDS triple-frequency observations collected from globally distributed MGEX stations (Xiao et al. 2019a). The estimated FCB shows a good consistency with the input float ambiguities. The RMS of Galileo FCB residuals is 0.05 cycles, while that of BDS is 0.08 cycles. It is also observed that the residuals are smaller for the combinations with larger wavelengths. The results indicate that there may exist ionospheric errors, and combinations are required to reduce their influence. The average STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 0.10 cycles, as shown in Figure 8 and Figure 9. To reduce the communication with servers for real time applications, it would be more efficient to broadcast linear combined FCB. When comparing the results from multi-GNSS, it can be seen that the STDs of Galileo EWL/WL FCBs are smaller than those of GPS and BDS, while that of the Galileo NL FCB is worse than GPS and BDS. The better quality of Galileo EWL/WL FCBs is likely attributed to the multipath suppression of Galileo signals, while the worse quality of Galileo NL FCB is due to the poor precision of satellite orbit and clock product. It is found that the STDs of Galileo raw FCBs are smaller than that of GPS and BDS, regardless of combinations for most of the days.

**Figure 8** Mean STD of the combined FCB series for all 51 days. Daily STD is calculated for each satellite FCB series. For each day, the mean STD of all satellite daily STDs is presented.

**Figure 9** Mean STD of the raw FCB series for all 51 days. Daily STD is calculated for each satellite FCB series. For each day, the mean STD of all satellite daily STDs is presented.

The performance of triple-frequency PPP AR is assessed with 11 days of data in threehour sessions. Compared with the float solutions, the positional biases of AR solutions are significantly improved, as shown in Figure 10 and Figure 11. The improvements of east, north, up components for positioning with BDS are 49.2%, 38.3%, and 29.6%, while those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These results demonstrate the efficiency of the proposed FCB estimation approach, and that the triple-frequency PPP AR can bring an obvious benefit to the float solution.

**Figure 10** Convergence performance of BDS triple-frequency PPP float and AR solutions based on 804 3-h sessions under 68% confidence level.

**Figure 11** Convergence performance of Galileo triple-frequency PPP float and AR solutions based on 5805 3-h sessions under 68% confidence level.

When comparing triple-frequency PPP AR with that of dual-frequency, it is found that the contribution of the third frequency observations is not remarkable. The insignificant improvement of the third frequency observable may be due to its limited contribution to satellite geometry and the narrow deployment with respect to the second carrier frequency. Nevertheless, adding the third frequency increases the reliability since it is observed that the number of successful sessions is increased. It is noted that all the experiments were conducted with excellent observational condition, e.g., the state-of-the-art receiver and antenna, open environment and static mode. Under this circumstance, the dual-frequency observations are of high quality and sufficient for high-precision positioning, which also limits the contribution of the additional signals. More obvious improvements could be expected for low cost receiver and antenna, restricted environment or kinematic mode (Li 2018).

#### **4.4 Uncombined cycle slip detection and repair**

As has been proved in the above research, the uncombined PPP is a good choice for the integration of multi-frequency and multi-GNSS observations. However, the uncombined signal processing brings new challenges for cycle slip detection, which is the carrier frequency identification of cycle slip (Xiao et al. 2018a). We presented an improved approach based on a time-differenced model for cycle slip detection and repair in uncombined PPP.

For cycle slip detection, the proposed approach significantly reduces false alarms as shown in Figure 12 and Figure 13. Having access to a reliable cycle slip detection method greatly reduces the number of ambiguity parameters to be estimated for processing uncombined GNSS measurements.

**Figure 12** Results of cycle slip detection for rover (left) and reference stations (right) using geometry-free combinations (1, 1, -2) and (1, -2, 1) (Wu et al. 2010b; Zhen et al. 2008). Red triangle, green square, and blue inverted triangle represent B1, B2 and B3 BDS carrier frequency observations, respectively.

**Figure 13** Results of cycle slip detection for rover (left) and reference (right) stations using the improved algorithm.

Compared with the routine time-differenced methods (Banville and Langley 2013; Zhang and Li 2016), the proposed approach can not only identify the carrier frequency in which the cycle slip occurs, but also makes it possible to separate the observation at other frequency of the same satellite without cycle slip, which greatly contributes to cycle slip estimations. Simulation results show that the theoretical success rates increase to 99.99% for both dual and triple frequency observations, as shown in Figure 14. Results from a real dataset also indicate that much stronger model strength of cycle slip estimation can be achieved.

**Figure 14** The ambiguity dilution of precision (ADOP) based success rates of dual-frequency (left) and triple-frequency (right) observations. The yellow surface represents the routine time-differenced model, while the green one represents the improved algorithm.

Although the results shown above are obtained with BDS triple-frequency observations, the method has been implemented in our software and used for all the researches conducted in this thesis. It is proved that the method is easily applicable to dual-frequency and multi-frequency observations for all GNSS.

#### **4.5 Program development**

Program development is also an important contribution of this thesis. Two software suites, i.e., PPP ambiguity resolution engine and FCB estimation server, have been created. The PPP ambiguity resolution engine is developed based on RTKLIB (http://www.rtklib.com/), which is an open source program package for standard and precise positioning with GNSS written in C programming language. The RTKLIB package includes a portable program library and several application programs. The original RTKLIB provides excellent support for GNSS data processing, e.g., standard data format and basic functions for PPP. During my research, I have significantly extended the original RTKLIB software. These improvements mainly include:


In addition, an FCB estimation software is developed with Python programming language, accompanied by a FCB analysis tool based on MATLAB programming language. The software supports FCB generation for both dual-frequency ionospheric-free PPP and multi-frequency uncombined PPP models. In order to investigate the characteristics of the FCB product, a detail report, which includes quality control information and residuals statistics, is also generated. Besides the two software suites, various scripts have been developed for the efficient processing of huge GNSS data sets, such as automatic data preparing, batch processing and analysis of results.

## **5 Publication overview**

This thesis includes the following four peer-reviewed publications (Xiao et al. 2019a; Xiao et al. 2019b; Xiao et al. 2018a; Xiao et al. 2018b):

5) Xiao G, Li P, Gao Y, Heck B (2019a) A Unified Model for Multi-Frequency PPP Ambiguity Resolution and Test Results with Galileo and BeiDou Triple-Frequency Observations. Remote Sensing 11(2):116 doi:10.3390/rs11020116

**Author's contributions:** The first author derived the unified model for multi-frequency PPP AR. The first and second author conceived and designed the experiments. The first author collected the data, performed the experiments, analyzed the results, and wrote the paper. All authors reviewed the paper.

6) Xiao G, Li P, Sui L, Heck B, Schuh H (2019b) Estimating and assessing Galileo satellite fractional cycle bias for PPP ambiguity resolution. GPS Solutions 23:3 doi:10.1007/s10291-018-0793-z

**Author's contributions:** The first author conceived the presented idea. The first and second author designed the experiments, worked out the technical details. The first author performed the numerical calculations for the experiments, analyzed the data, and wrote the paper. All authors reviewed the paper.

7) Xiao G, Sui L, Heck B, Zeng T, Tian Y (2018b) Estimating satellite phase fractional cycle biases based on Kalman filter. GPS Solutions 22:82 doi:10.1007/s10291-018- 0749-3

**Author's contributions:** The first author devised the main conceptual ideas and designed the proposed Kalman filter for FCB estimation. The first author performed the experiments, analyzed the data, and wrote the paper. All authors reviewed the paper and helped shape the research.

8) Xiao G, Mayer M, Heck B, Sui L, Zeng T, Zhao D (2018a) Improved time-differenced cycle slip detect and repair for GNSS undifferenced observations. GPS Solutions 22:6 doi:10.1007/s10291-017-0677-7

**Author's contributions:** The first author designed the model and the experiments. The first author performed the experiments, analyzed the data, and wrote the paper. All authors discussed the results, reviewed and contributed to the final manuscript.

## **6 Conclusions and outlook**

### **6.1 Conclusions**

For decades, GPS PPP has found many scientific and industrial applications due to its cost-effectiveness, global coverage, and high accuracy. However, it suffers from a few drawbacks which prevents more applications, e.g., slow convergence, availability and reliability. The rapid development of multi-GNSS and multi-frequency signals provides an excellent opportunity for improvements. In addition, the integer ambiguity resolution technique is expected to further enhance the accuracy and shorten the convergence time. In the framework of this thesis, multi-frequency and multi-GNSS measurements are properly integrated for PPP with ambiguity resolution to achieve state-of-the-art fast and accurate positioning, which provides an important contribution to GNSS precise positioning and applications.

Integer ambiguity resolution is the key issue to improve the PPP positioning quality. The core of PPP ambiguity resolution is carrier phase biases estimation. In this thesis, the characteristics of phase fractional cycle biases are analyzed with the updated and long time series. Results indicate that the stability of phase fractional cycle biases has been improved in recent years. Taking advantage of this property, an improved FCB estimation method based on Kalman filter is proposed. In the proposed Kalman filter, the large number of observations is handled epoch-by-epoch, which significantly reduces the dimension of the involved matrix and accelerates the computation. Hence it outperforms the commonly used least squares method in terms of efficiency. In order to avoid iterations caused by the one cycle inconsistency among FCB measurements, a pre-elimination method is developed based on the temporal stabilities of satellite FCBs. The pre-elimination method shows a clear advantage over post-residual adjustment, which further improves the efficiency. The fast estimation of FCB is especially useful for real-time application. Results from about 200 IGS stations indicate that the determined GPS satellite FCBs show a good consistency with existing FCB products, e.g., from IGS analysis center CNES-GRG and Wuhan University.

With the estimated FCB products, the PPP ambiguity resolution technique with the current Galileo and BDS constellation is verified. The satellite fractional cycle biases are thoroughly assessed by a comparison with those of GPS in terms of data usage rate, residual distribution, as well as standard deviation of daily estimates. Results indicate that the quality of Galileo WL FCB is better than for GPS and BDS. We attribute the good quality of Galileo WL FCB to its advanced signal modulation, AltBOC, which significantly compresses the multipath effect for code measurements. Within the Galileo constellation, the quality of FOC WL FCB is much better than for IOV satellites. The poorer performance of IOV satellites WL FCB is a result of a reduction in the satellite transmission signal power. The performance of the two satellites with highly eccentric orbits is comparable to other FOC satellites but having a smaller number of observations. As for NL FCB, the quality of Galileo NL FCB is slightly worse than that of GPS but better than that of BDS. Since the accuracy of NL FCB estimation is dominated by unmodeled errors in float PPP, the worse quality of Galileo NL FCB is likely caused by the imperfect antenna models. Within Galileo, the NL FCB quality of FOC and IOV (except E19), as well as the two eccentric satellites, shows no significant difference in terms of data usage rates and residuals. The reason for the worse performance of E19 is still not clear. On the one hand, it cannot, or at least not fully, be ascribed to the signal transmission power as the power difference between FOC and IOV E11/E12 is larger than the difference between E11/E12 and E19. On the other hand, the worse quality of the E19 satellite orbit, indicated by SLR residuals, could also degrade NL FCB estimations. This issue remains an open question and deserves further investigation. Galileo PPP AR solutions are conducted at 20 MGEX stations with three-hour sessions for ten days. The positional biases of AR solutions are 0.7, 0.6 and 2.1 cm for east, north and up components respectively, while those for float solutions are 2.1, 1.1 and 2.7 cm, corresponding to improvements of 67, 45 and 22% respectively. These results demonstrate that the Galileo observations can bring an obvious benefit to ambiguity-fixed PPP.

PPP ambiguity resolution is also extended to multi-frequency observations in this thesis. In order to properly integrate the multi-frequency observations, a unified uncombined PPP ambiguity resolution model based on raw observations is proposed. To verify its efficiency, we processed 51 days of Galileo and BDS triple-frequency observations collected from globally distributed MGEX stations. The estimated FCB show a good consistency with the input float ambiguities. The RMS of Galileo FCB residuals is 0.05 cycles, while that of BDS is 0.08 cycles. It is also observed that the residuals are smaller for the combinations with larger wavelengths. The results indicate that there may exist ionospheric errors, and linear combinations are required to reduce their influence. The average STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 0.10 cycles. To reduce the communication with servers for real time applications, it would be more efficient to broadcast linear combined FCB. The performance of triple-frequency PPP AR is assessed with 11 days of data in three-hour sessions. Compared with the float solutions, the positional biases of AR solutions are significantly improved. The improvements of ENU components for positioning with BDS are 49.2%, 38.3%, and 29.6%, while those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These results demonstrate the efficiency of the proposed FCB estimation approach, and that the triple-frequency PPP AR can bring an obvious benefit to the float solution. When comparing triple-frequency PPP AR with that of dual-frequency, it is found that the contribution of the third frequency observations is minimal. The insignificant improvement of the third frequency observable may be due to its limited contribution to satellite geometry and the narrow deployment with respect to the second carrier frequency. Nevertheless, adding the third frequency increases the reliability since it is observed that the number of successful sessions is increased.

The uncombined signal processing brings new challenges for cycle slip detection. The most significant feature is carrier frequency identification of cycle slips. Since all carrier frequency observations are processed separately, it is essential to identify the carrier frequency of the cycle slips to avoid contaminating other observations. To provide continuous carrier phase measurements for the above research, an improved approach based on a time-differenced model for cycle slip detection and repair is proposed, which reduces the false alarms and increases the success rate of cycle slip estimation. The simulation results show that the theoretical success rates increase to 99.99% for both dual and triple frequency observations. Results from a real dataset also indicate that much stronger model strength of cycle slip estimation can be achieved.

In summary, the GPS only PPP AR has been extended to multi-frequency and multi-GNSS PPP ambiguity resolution with improved accuracy and fast convergence in this thesis. This is accomplished by a unified model based on the uncombined PPP. The proposed model has been carefully studied and enriched with improved cycle slip detection and repair, fast FCB estimation and ambiguity resolution methods.

### **6.2 Outlook**

The outlook and perspectives of this thesis are presented as follows:


## **References**


Zumberge JF, Heflin MB, Jefferson DC, Watkins MM, Webb FH (1997) Precise point positioning for the efficient and robust analysis of GPS data from large networks. Journal of Geophysical Research: Solid Earth 102(B3):5005-

# **Chapter 2 Estimating satellite phase fractional cycle biases based on Kalman filter**

#### **Guorui Xiao1,2 · Lifen Sui 2 · Bernhard Heck<sup>1</sup> · Tian Zeng<sup>2</sup> · Yuan Tian<sup>2</sup>**

<sup>1</sup> Geodetic Institute, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany <sup>2</sup> Zhengzhou Institute of Surveying and Mapping, Zhengzhou 450052, China

Submitted: 21 August 2017; Accepted: 18 June 2018; Published: 23 June 2018

GPS Solutions 22:82, 2018. DOI:10.1007/s10291-018-0749-3

This is an author-created version of the article with permission from Springer. The final publication is available at link.springer.com

**Abstract** Phase fractional cycle biases (FCBs) originating from satellites and receivers destroy the integer nature of PPP carrier phase ambiguities. In order to achieve integer ambiguity resolution of PPP, FCBs of satellites are required. In former work, least squares methods (LSM) are commonly adopted to isolate FCBs from a network of reference stations. However, it can be extremely time-consuming concerning the large number of observations from hundreds of stations and thousands of epochs. In addition, iterations are required to deal with the one cycle inconsistency among FCB measurements. We propose to estimate the FCB based on a Kalman filter. The large number of observations are handled epoch by epoch, which significantly reduces the dimension of the involved matrix and accelerates the computation. In addition, it is also suitable for realtime applications. As for the one cycle inconsistency, a pre-elimination method is developed to avoid iterations and posterior adjustments. A globally distributed network consisting of about 200 IGS stations is selected to determine the GPS satellite FCBs. Observations recorded from DoY 52 to 61 in 2016 are processed to verify the proposed approach. The RMS of wide lane (WL) posterior residuals is 0.09 cycles while that of the narrow lane (NL) is about 0.05 cycles, which indicates a good internal accuracy. The estimated WL FCBs also have a good consistency with existing WL FCB products (e.g., CNES-GRG, WHU-SGG). The RMS of differences with respect to GRG and SGG products are 0.03 and 0.05 cycles. For satellite NL FCB estimates, 97.9% of the differences with respect to SGG products are within ± 0.1 cycles. The RMS of the difference is 0.05 cycles. These results prove the efficiency of the proposed approach.

**Keywords** Precise point positioning · Integer ambiguity resolution · Fractional cycle bias · Kalman filter

## **1 Introduction**

Integer ambiguity resolution (AR) can significantly shorten the convergence time and improve the accuracy of Precise Point Positioning (PPP). Phase fractional cycle biases

(FCBs) originating from satellites and receivers destroy the integer nature of PPP carrier phase ambiguities. The receiver FCB can be eliminated by single differencing across satellites or assimilated into the receiver clock parameter by forcing one zero difference ambiguity to its nearest integer, while the satellite FCBs must be estimated from a network of reference stations (Gabor and Nerem 1999). With the additional satellite FCB estimates, PPP users are able to remove satellite FCBs and recover the integer nature of ambiguities.

The fact that double-differenced ambiguities in global or regional networks can be resolved to integer values lays the foundation of integer ambiguity resolution for PPP. The resolved double-differenced ambiguity implies that the fractional parts of two singledifferenced ambiguities (across satellites) must agree well with each other. By estimating the fractional parts at the server end and applying them to single differencing PPP at the user end, PPP integer ambiguity resolution can be achieved. In this sense, PPP integer ambiguities are in fact equivalent to double-differenced ambiguities (Teunissen and Khodabandeh 2015). Ge et al. (2008) proposed to estimate the common fractional parts of the single-differenced ambiguities in wide lane (WL) and narrow lane (NL) form. The fractional parts, denoted as single-differenced FCB (SD-FCB), were estimated by averaging the fractional parts of all involved float single-differenced WL and NL ambiguity estimates. Based on the same principle but instead of the averaging process, a least squares method in an integrated adjustment was adopted to enhance the estimates (Li and Zhang 2012; Zhang and Li 2013). Since January 1, 2015, the School of Geodesy and Geomatics at Wuhan University (WHU-SGG) routinely releases GPS WL and NL FCB products with open access (Li et al. 2015a). Similar approaches have also been applied to BDS (Li et al. 2017a; Wang et al. 2017), Galileo (Tegedor et al. 2016) and GLONASS (Geng and Shi 2016). In order to exploit the ionosphere characteristics, the model has been extended to deal with GPS L1 and L2 raw observables (Gu et al. 2015b; Li et al. 2013a) and BDS triple-frequency observables (Gu et al. 2015a).

Different from the above approaches, Collins et al. (2008) developed a decoupled clock model by separating satellite clocks for code and phase observations. Similarly, Laurichesse et al. (2009) developed an integer phase clock model in which the NL FCBs were assimilated into receiver and satellite clock estimates of a network solution. This model has been employed to generate the precise satellite clock products by Groupe de Recherche de Géodésie Spatiale of the Centre National d'Etudes Spatiales (CNES-GRG) (Loyer et al. 2012). The only difference between the decoupled clock model and the integer phase clock model is the approach for determining the WL ambiguity. The integer phase clock model utilizes WL FCB corrections and a satellite-averaging process to fix the integer WL ambiguity, whereas the decoupled clock model directly estimates the integer WL ambiguity along with other unknowns through the functional model. The difference between the WL/NL FCB model and the integer phase clock model is the strategy of separating NL FCBs from integer ambiguities. NL FCBs are directly estimated in the WL/NL FCB model, whereas they are assimilated into the clock estimates in the integer phase clock model.

These PPP AR techniques are proven to be equivalent in theory (Shi and Gao 2014). The systematic biases between position estimates have been demonstrated to be minimal and negligible with data from a global network of reference stations (Geng et al. 2010). The WL/NL FCB model can conveniently supplement current network solutions as an additional software module, as the FCB determination is compatible with current official clock-generation methods. In contrast, the integer phase clock products are incompatible with current clock products although they may perform slightly better in practice.

Based on the review of existing work, we find that the least squares method (LSM) is routinely utilized for FCB estimation. However, the LSM can be extremely time-consuming concerning the large number of observations from hundreds of reference stations and thousands of epochs during the FCB estimation. This is unfavorable for the more and more popular real-time applications. In addition, iterations are required to deal with the one cycle inconsistency among FCB measurements. Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU-SGG and CNES-GRG products respectively, the one cycle inconsistency arises whenever the superposition of receiver FCB and satellite FCB exceeds the boundary. The additional iterations of LSM and computation of a large matrix demand a long time to finish. Therefore, a fast and efficient estimation method is desirable.

In this contribution, a Kalman filter is employed to speed up the computation. The large number of observations are handled epoch by epoch, which significantly reduces the dimension of the involved matrix and accelerates the computation. The recursive computation of the Kalman filter allows real-time applications. As for the one cycle inconsistency, a pre-elimination method is developed to avoid the posterior adjustments and iterations of the whole process. The following section describes the theoretical background of FCB estimations and Kalman filter. Then a pre-elimination method of one cycle inconsistency is proposed following the analysis of temporal stabilities of satellite FCBs. A set of FCB products is determined and evaluated by comparing with those of CNES and WHU. With the estimated FCBs, the improvements from PPP AR are assessed. Finally, the methodology is summarized, and an outlook for future research is presented.

## **2 Methodology**

We start with the basic PPP float mode followed by a description of FCB estimation strategy. The pre-elimination method of one cycle inconsistency is elaborated as well as an introduction to Kalman filter.

#### **2.1 PPP float mode**

In GNSS dual-frequency PPP, the ionospheric-free (IF) combination is routinely employed to eliminate the effect of the first-order ionospheric delay. For a satellite *s* observed by receiver *r*, the corresponding pseudorange and carrier phase observation equation can be expressed as

$$\begin{cases} P\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT + D\_{r,IF} - D\_{IF}^{\mathcal{S}} + \varepsilon\_{P\_{IF}}\\ \Phi\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT + \lambda\_{IF} \left( N\_{r,IF}^{\mathcal{S}} + B\_{r,IF} - B\_{IF}^{\mathcal{S}} \right) + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{1}$$

where indicates the geometric distance between satellite and receiver; and are the clock errors of receiver and satellite; is the slant tropospheric delay; , and are the receiver and satellite specific code hardware delays; and , are the wavelength in meters and ambiguity in cycles; , and are the receiver-dependent and satellite-dependent uncalibrated phase delays; , are the pseudorange and carrier phase measurement noise.

Conventionally, precise orbit and clock products from the IGS analysis center are used to remove satellite orbit and clock errors. During the generation of IGS precise products, the pseudorange ionospheric-free hardware delay bias is assimilated into the clock offset in accordance with the IGS analysis convention. Due to the fact that pseudorange measurements provide the reference to clock parameters, the actual receiver clock estimate would absorb the ionospheric-free combination of the receiver pseudorange hardware delay ,. After applying the GNSS precise satellite clock products, equation (1) can be rewritten as

$$\begin{cases} P\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dT + \overline{d}t\_r + \varepsilon\_{P\_{IF}} \\ \Phi\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dT + \overline{d}t\_r + \lambda\_{IF} \overline{N}\_{r,IF}^{\mathcal{S}} + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{2}$$

where the receiver clock and ambiguity can be re-parameterized as

$$\overline{dt}\_r = \left(dt\_r + D\_{r,IF}\right) \tag{3}$$

$$
\overline{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + b\_{r,IF} - b\_{IF}^{s} \tag{4}
$$

$$b\_{r,IF} = B\_{r,IF} - \frac{D\_{r,IF}}{\lambda\_{IF}} \tag{5}$$

$$h\_{IF}^{\mathcal{S}} = B\_{IF}^{\mathcal{S}} - \frac{D\_{IF}^{\mathcal{S}}}{\lambda\_{IF}} \tag{6}$$

For ambiguity-float PPP solutions, the ionospheric-free ambiguity parameter , is estimated as a real-value constant. Since the estimated ambiguity parameter is a combination of the integer ambiguity, the corresponding code hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite ends, the integer property is lost.

#### **2.2 FCB estimation**

For an ambiguity-fixed PPP solution, , is usually decomposed into the following combination of integer WL and float NL ambiguities for ambiguity fixing

$$
\overline{N}\_{r,IF}^{s} = \frac{\frac{cf\_2}{f\_1^2 - f\_2^2} N\_{r,WL}^{s} + \frac{c}{f\_1 + f\_2} \overline{N}\_{r,NL}^{s}}{\lambda\_{IF}} \tag{7}
$$

Note that , is an integer WL ambiguity, which implies that WL FCB will not directly contribute to PPP. The WL ambiguity can be resolved by the Hatch–Melbourne–Wübbena (HMW) combination observable (Hatch 1982; Melbourne 1985b; Wübbena 1985a),

$$HMW = \frac{f\_1 \Phi\_{r,1}^s - f\_2 \Phi\_{r,2}^s}{f\_1 - f\_2} - \frac{f\_1 P\_{r,1}^s + f\_2 P\_{r,2}^s}{f\_1 + f\_2} = \lambda\_{r,WL}^s \left\{ N\_{r,WL}^s + d\_{r,WL} - d\_{WL}^s \right\} \tag{8}$$

where

$$d\_{r,WL} = B\_{r,1} - B\_{r,2} - \frac{f\_1 D\_{r,1} + f\_2 D\_{r,2}}{\lambda\_{r,WL}^s (f\_1 + f\_2)} \tag{9}$$

$$d\_{WL}^{\rm s} = B\_1^{\rm s} - B\_2^{\rm s} - \frac{f\_1 D\_1^{\rm s} + f\_2 D\_2^{\rm s}}{\lambda\_{r,WL}^{\rm s} (f\_1 + f\_2)} \tag{10}$$

$$\overline{N}\_{r,WL}^{s} = N\_{r,WL}^{s} + d\_{r,WL} - d\_{WL}^{s} \tag{11}$$

If the WL ambiguity is correctly resolved to an integer value based on (8), the NL ambiguity observable can be derived,

$$
\overline{N}\_{r,NL}^{s} = N\_{r,1}^{s} + d\_{r,NL} - d\_{NL}^{s} \tag{12}
$$

where

$$d\_{r,NL} = \frac{f\_1 + f\_2}{c} \left(\lambda\_{IF} B\_{r,IF} - D\_{r,IF}\right) \tag{13}$$

$$d\_{NL}^{\rm s} = \frac{f\_1 + f\_2}{c} (\lambda\_{IF} B\_{IF}^{\rm s} - D\_{IF}^{\rm s}) \tag{14}$$

Equations (11) and (12) serve as the basic model for estimating FCBs. Since they have the same structure, a general expression can be formulated as

$$R\_r^s = \overline{N}\_r^s - N\_r^s = d\_r - d^s \tag{15}$$

for WL and NL linear combinations, respectively. represents the FCB measurements; denotes the float undifferenced ambiguities; denotes the integer part of , which is the sum of the original integer ambiguity and the integer part of the combined code hardware delays and uncalibrated phase delays from both receiver *r* and satellite *s*; and denote the receiver and satellite FCBs.

A set of equations in the form of (15) can be integrated based on a network of reference stations. Suppose that there are *n* satellites tracked at *m* reference stations, the system of equations can be expressed as

$$
\begin{bmatrix} R\_{r,1}^{s,1} \\ \vdots \\ R\_{r,1}^{s,n} \\ \vdots \\ R\_{r,m}^{s,1} \\ \vdots \\ R\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} \overline{N}\_{r,1}^{s,1} - N\_{r,1}^{s,1} \\ \vdots \\ \overline{N}\_{r,1}^{s,n} - N\_{r,1}^{s,n} \\ \vdots \\ \overline{N}\_{r,m}^{s,1} - N\_{r,m}^{s,1} \\ \vdots \\ \overline{N}\_{r,m}^{s,n} - N\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} 1 & \cdots & 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cdots & 0 & 0 & \cdots & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & \cdots & -1 \end{bmatrix} \begin{bmatrix} d\_{r,1} \\ \vdots \\ d\_{r,m} \\ d\_{r,1}^{s,n} \\ \vdots \\ d\_{r,n}^{s,n} \end{bmatrix} \tag{16}
$$

The obtained system of equations is singular. The integer ambiguities need to be determined for each equation on the left side, while one arbitrarily selected FCB should be set to zero on the right side. Assuming the float ambiguities ̅ are precisely estimated, can be determined by rounding ̅ . In this way, the system of equations can be resolved. However, the rounding approach may introduce one cycle inconsistencies.

#### **2.3 One cycle inconsistency**

Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU-SGG and CNES-GRGS products respectively, the one cycle inconsistency arises whenever the superposition of receiver FCB and satellite FCB exceeds the boundary. A simple example is presented to depict the situation. Assuming there is one satellite ( = 0.2) tracked at two stations (,1 = 0.6 and ,2 = 0.8), the superposition of FCBs should be (,1 = −0.4 and ,2 = −0.6) . However, the ,2 would be 0.4 due to the rounding approach. As a consequence, the corresponding satellite FCB derived from the two stations differs one cycle. The one cycle inconsistency can be detected by examining the posterior residuals, as employed in former works. The posterior adjustments are inefficient as iterations of the whole process are required. In contrast, we propose to eliminate the inconsistency in advance, which shows a clear advantage over posterior adjustment. The proposed method consists of the following steps:

(1) For all the FCB measurements at a single station, the satellite FCBs are eliminated using previous estimates. The underlying assumption is that satellite FCBs are stable over successive epochs and can be eliminated to a large extent by previous estimates, which will be proved in the next section.

(2) The residual of each FCB measurement after the subtraction of the satellite FCB yields an initial estimate of the receiver FCB. A set of initial receiver FCB estimates could be obtained from all simultaneously observed satellites.

(3) In theory, the receiver FCBs obtained from all simultaneously observed satellites are expected to be consistent. Therefore, a one cycle inconsistency can be detected by examining the group of initial receiver FCBs. If an inconsistency exists, the corresponding receiver FCB will differ by ±1 cycle.

(4) To compensate the one cycle inconsistency in receiver FCBs, the corresponding integer ambiguity is adjusted by ±1 cycle.

(5) The procedures described above can be performed in an iterative way until a consistent receiver FCB is obtained from all satellite measurements.

In addition, stations can also be handled individually to obtain clean FCB measurements. After the pre-elimination of the one cycle inconsistency, the clean FCB measurements are fed to the Kalman filter.

### **2.4 Kalman filter**

The Kalman filter addresses the general problem of state estimates of discrete time controlled processes that are governed by a linear stochastic difference equation (Kalman 1960). The theory has been well studied and widely applied (Yang 2006; Yang et al. 2010). Since the Kalman filter is based upon the theory of least squares, it is theoretically possible to calculate the same solution as LSM and versa vice. However, the design and normal equation matrix will be huge in LSM, considering the large number of observations from hundreds of reference stations and thousands of epochs. Computation of the large matrix is time-consuming, which is inefficient and unfavorable for real-time applications. In contrast, the large number of observations are handled epoch-by-epoch in a Kalman filter, which significantly reduces the dimension of the involved matrix and accelerates the computation. One additional advantage of Kalman filter is its real-time capability. This is of particular interest for real-time PPP AR and its applications. Therefore, a Kalman filter is employed in this work. The design and flowchart of a Kalman filter are depicted in Figure 1. In our case, is the input FCB measurement for the Kalman filter, while and are the output estimates.

**Figure 1** Flow chart of the proposed Kalman filter. The dashed blocks represent additional steps adopted in this study

Two additional steps are added to the standard Kalman filter. The first step aims at establishing the dynamic model and determining the state transition matrix by analyzing the temporal stability of existing FCB products. Satellite WL FCBs are stable over several days, which can be characterized as a constant parameter on a daily basis. Satellite NL FCBs are considered as constant within 15 minutes but exhibit small variations over 15 minute intervals. Therefore, NL FCBs are characterized as random walk processes. The second step is introduced to eliminate the one cycle inconsistency of FCB measurements. Only clean measurements are sent to the Kalman filter to avoid iterations. Note that a constraint is imposed on the Kalman filter to eliminate the rank deficiency. A satellite FCB is selected and set to zero, which is accomplished by adding a pseudo observational error equation with zero variance (Yang et al. 2010).

## **3 Results, comparison, and analysis**

In this section, the temporal stability of FCB is first analyzed based on existing products. Then, the proposed pre-elimination method of one cycle inconsistency is demonstrated. Ten sets of FCBs are computed and evaluated with respect to existing products. At last, with the estimated FCBs, the improvements from PPP AR are assessed.

#### **3.1 Temporal stability of FCB**

In order to study the characteristics of satellite FCBs, products from CNES GRG and WHU SGG are employed. Both organizations provide WL FCBs derived from the same strategy, while NL FCBs are only available for SGG products since they are absorbed by satellite phase clock offsets in GRG products. The products for entire 2016 are downloaded and analyzed in this study. Due to the rank deficiency of FCB equations, a receiver FCB for GRG products and a satellite FCB for SGG products are selected as references. Therefore, a single differencing process across satellites is performed to remove the datum before comparison. The single-differenced WL FCBs are presented in Figure 2. Note that the two products have opposite signs.

**Figure 2** Single-differenced WL FCBs of CNES GRG (top) and WHU SGG (bottom) products in 2016. The lines represent individual satellite WL FCBs with respect to that of G01

The WL FCBs of satellites are remarkably stable over days with standard deviations (STDs) of less than 0.08 cycles for both products. An exception is observed for G32 satellite. The STD of G32 is significantly larger than the others with 0.2 cycles. The reason is that a new Block IIF satellite (G32) was launched on 5 February 2016 and activated since 10 February. After the replacement, the WL FCB exhibits a similar temporal stability as the other satellites. The statistics further confirm the conclusion that the satellite WL FCB is stable over at least several days. Therefore, WL FCBs are estimated as constant on a daily basis in the proposed Kalman filter.

NL FCB products are estimated with respect to specific IGS precise clock products in the WL/NL FCB based method. Due to the daily boundaries of precise clock products, NL FCB estimates are only continuous within one day. The single-differenced SGG NL FCBs on DoY 001, 2016, are presented in Figure 3. The top panel shows the raw SD NL FCBs for each satellite while G01 is selected as reference. It can be seen that the datum changes frequently, which introduces a difference of one cycle between two consecutive epochs. After adjustments, the NL FCB series are continuous, as presented in the middle panel.

One can easily discern that satellite NL FCBs are also stable. The STDs of all single-differenced satellite NL FCBs are within ± 0.05 cycles for DoY 001. In addition, daily STDs of all satellite NL FCBs for the whole year 2016 are calculated. The bottom panel depicts the distribution of daily STDs. 98.37% of the daily STDs are below 0.15 cycles while 97.02% are below 0.1 cycles. These results show that satellite NL FCBs may be more stable than reported in former research. Note that a small number of daily STDs may reach 0.5 cycles. These abnormal STDs may be caused by maneuvers during satellite eclipse or unmodeled errors. And we could not exclude the possibility that there may be blunders. Nevertheless, satellite NL FCBs are modeled as random walk processes in the proposed Kalman filter. Also note that the characteristics of receiver NL FCB depend on the receiver type, environment, and other factors. It is hard to characterize receiver NL FCBs with a general model. Therefore, they are modeled as white noise. The temporal stability of satellite FCBs also implies that it can be utilized for pre-elimination of the one cycle inconsistency.

**Figure 3** Raw (top) and adjusted (middle) Single-differenced NL FCBs of SGG products on DoY 001. The bottom panel presents the histogram of all daily STDs in 2016.

#### **3.2 Pre-elimination of the one cycle inconsistency**

For a successful pre-elimination of the one cycle inconsistency, the satellite FCBs should be removed in advance. The key issue is if or at what extent the satellite FCBs of the current epoch can be counteracted by the previous estimation. In the proposed pre-elimination method, the differences of satellite FCBs between two successive epochs will be absorbed by receiver FCB. Therefore, epoch differences of satellite FCBs should be small enough to avoid impacts on the detection of the one cycle inconsistency. In order to validate the hypothesis, an additional epoch differencing process is performed on the single-differenced FCBs. The histograms of epoch differences for WL FCBs in 2016 are presented in Figure 4.

**Figure 4** Histograms of epoch differenced WL FCBs of GRG (top) and SGG (bottom) products in 2016

99.78% and 95.81% of the differences are within ± 0.05 cycles for GRG and SGG products respectively. The max differences, with a magnitude of 0.386 cycles, are found in SGG products. Since these values are much smaller than one cycle, it will scarcely affect the detection and elimination of the one cycle inconsistency. Note that the distribution of epoch differences for GRG products is more concentrated around zero than that of SGG, which may indicate better quality.

The epoch differencing process is also conducted for NL FCB products, as shown in Figure 5. 99.73% of the differences are within ± 0.05 cycles for SGG products. The result is even better than that of WL FCB. This can be attributed to the fact that the intervals (15 minutes) of NL FCB are much shorter than those of WL FCB (24 hours).

**Figure 5** Histogram of epoch differenced NL FCBs of SGG products

**Figure 6** Pre-elimination of the one cycle inconsistency for station IQQE on 2016/02/26. The top panels represent receiver WL FCB adjustment, while the bottom ones represent that of receiver NL FCB. The left and right panels represent the raw and adjusted measurements respectively

According to the above analysis of satellite FCB products, it can be concluded that the major parts of both WL and NL satellite FCBs can be mitigated by previous estimations. After the removal of satellite FCBs, the one cycle inconsistency among FCB measurements can be detected by examining the residuals. After the detection of inconsistency,

it is eliminated by adjusting the corresponding integer ambiguity. An example of the process is shown in Figure 6. One can easily observe that the one cycle inconsistencies occur on several satellite measurements and can be effectively eliminated by the proposed method. It should be noted that elimination of the one cycle inconsistency will not improve the precision of FCB estimation as the fractional part remains the same. In this sense, the magnitude of the epoch difference is less of concern as long as it is sufficient to detect potential inconsistencies. Also, note that this can be accomplished by posterior residuals adjustments. However, pre-elimination shows a clear advantage over posterior adjustment in terms of efficiency.

In this way, all epochs and all station data can be processed individually to obtain clean FCB measurements. For post-processing, the above procedure can be simplified. Daily WL FCB measurements can be adjusted simultaneously as the receiver WL FCBs are also constant within one day. For the initialization of the Kalman filter, a preliminary set of satellite FCBs should be provided for the first epoch. The initial satellite WL FCBs are taken from the estimations of the last day, while the initial satellite NL FCBs are determined with the route method (Li et al. 2015a).

### **3.3 Comparison of FCB products**

A globally distributed network consisting of about 200 stations is selected from IGS network, as shown in Figure 7. Observations from Day of Year (DoY) 52 to 61 (2016/02/21 - 2016/03/01) are processed to determine GPS satellite FCBs. GFZ final precise products are applied to remove satellite orbit and clock errors. The other errors are corrected according to the IGS standard error models (Kouba 2009). FCB estimations are conducted in three sequential steps. First, float ambiguities containing FCBs are obtained by HMW combinations and PPP processing from the network of reference stations. Secondly, FCB measurements are generated from the float ambiguities. Thirdly, the proposed Kalman filter is adopted to estimate FCBs from all FCB measurements.

The quality of FCB estimation can be indicated by the posterior residuals of FCB measurements. In general, a highly consistent FCB estimation is expected if the residuals are close to zero. Figure 8 shows the distribution of the residuals of FCB estimations for the ten days. It can be seen that both histograms are symmetric and nearly centered at zero. The Root Mean Square (RMS) of WL residuals is 0.09 cycles while the RMS of NL residuals is about 0.05 cycles, which indicates a good consistency between the estimated FCBs and the input float ambiguities. The total number of input float ambiguities is 934908. Figure 9 shows the usage rates of WL and NL float ambiguities for each satellite. The averaged usage rate for WL is 96.66%, while that for NL is 98.17%. We find that the quality of NL FCBs is better than that of WL. The possible reason is that, WL FCBs are estimated as daily invariants and easily affected by pseudorange noise, while the NL FCBs are updated every 15 minutes and derived from much more precise carrier phase measurements.

For better visualization, the derived satellite FCBs are shifted with integer cycles, as presented in Figure 10. It can be seen that the derived satellite FCBs show similar temporal stability as in the above analysis. Satellite WL FCBs exhibit extremely small variations over the ten days period, while NL FCB shows a larger variation over tens of minutes.

**Figure 7** Geographical distribution of the selected reference stations

**Figure 8** Distributions of posterior residuals after WL (top) and NL (bottom) FCB estimation

**Figure 9** Usage rates of the float WL and NL ambiguities

**Figure 10** Time series of satellite WL (top) and NL (bottom) FCBs for all the other 31 satellites while G01 is selected as reference

In order to validate the external accuracy of our estimation, the derived satellite FCBs are also compared with those of SGG and GRG. As discussed before, WL FCB measurements are obtained from HMW combinations and are relatively independent of PPP processing. The determined WL FCBs should be consistent across all products. Figure 11 presents the comparison of the derived satellite WL FCBs with those of GRG and SGG products.

**Figure 11** Histograms of the differences between our WL FCB estimation and SGG (top) and GRG products (bottom)

It can be seen that our WL FCB estimates show good consistency with those of SGG and GRG. All of the differences are within ± 0.1 cycles. Compared with SGG products, 76.7% of the differences are within ± 0.05 cycles with an RMS of 0.05 cycles. Compared with GRG products, 91.0% of the differences are within ± 0.05 cycles with an RMS of differences of 0.03 cycles. Based on the above analysis, we can safely conclude that there is no systematic bias between our WL FCB estimates and those of SGG and GRG. The differences are actually minimal and negligible. However, one can realize that our WL FCB

estimates agree better with GRG products than those of SGG. Note that the RMS of differences between SGG and GRG is 0.03 cycles. We suspect that the GRG products may perform slightly better in practice.

**Figure 12** Histogram of the differences between our NL FCB estimation and SGG products

NL FCBs are only available for SGG products. The derived satellite NL FCBs are compared with SGG products, as shown in Figure 12. It can be seen that our NL FCBs agree well with SGG NL FCB products. 97.9 % of the differences are within ± 0.1 cycles while 70.4% of the differences are within ± 0.05 cycles. The RMS of the differences is 0.05 cycles, corresponding to 5.1 mm. The discrepancy between the two results may be ascribed to different PPP processing strategies. Since NL FCB measurements are directly derived from PPP float ambiguities, discrepancies between error correction models, such as tropospheric models, may introduce small systematic biases. Another possible reason could be the different distributions of reference stations. Since FCB estimates are contaminated by unknown temporally and spatially correlated errors, these unknown errors may change with the distribution of reference stations. GRG products should be employed for comparison in further investigations. Note that the discrepancy will not affect PPP AR at user end if the same error model as at the server end is used.

#### **3.4 PPP AR solution**

In order to validate our FCB estimates and to assess the improvement by ambiguity resolution, all IGS network stations are processed in PPP float and AR mode with the obtained satellite FCB estimates. The 24 hours observations from over 500 stations in the ten days are divided into three-hour sessions. The solutions with data integrity less than 80 percent, wrongly fixed, and incomplete convergence are removed. After the pre-processing, there are 22953 sessions. The performance in terms of convergence time and positional accuracy is evaluated under different confidence levels for the sake of reliability (Lou et al. 2016).

**Figure 13** Convergence performance of GPS PPP float and AR solutions based on 22953 three-hour sessions under 95% confidence levels

**Table 1** Accuracy comparison of GPS PPP float and AR solutions based on 22953 three-hour sessions under different confidence level [Unit: cm]


Figure 13 presents the positional errors of GPS PPP float and AR solutions based on the statistics over all the sessions. It can be seen that the convergence time is significantly shortened by ambiguity resolutions, especially for the east and up components. It takes 64 minutes for float solutions to converge to three-dimensional 10 cm accuracy, while that for AR solutions is only 48.5 minutes, corresponding to an improvement of 24%. In order to assess the improvements with respect to time length, the positional errors are calculated for one-hour, two-hour, and three-hour solutions, as presented in Table 1. It can be seen that PPP ambiguity resolutions are able to enhance the accuracy for all the schemes. The most significant improvement is found for the east component. Due to the design of satellite navigation systems, the accuracy of the east component is worse than that of north component in mid-latitude because of weaker model strength. Ambiguity resolutions can improve this situation by imposing tight constraints on ambiguity parameters. However, the improvements decrease as the time length increases. The improvements of ENU components for one-hour solutions are (70%, 45%, 28%), while that for two-hour solutions and three-hour solutions are (62%, 32%, 21%) and (51%, 22%, 20%), respectively. The reason is that the model strength of the float solution is improved with more observations. Also note that the improvements differ with respect to the confidence level. They are less significant under 95% level, which implies the PPP AR is not effective for some stations. Since FCB estimates are contaminated by temporally and spatially correlated errors, a small and dense network is preferred for better performance.

## **4 Conclusions**

In this contribution, estimating satellite FCBs based on the Kalman filter is proposed and demonstrated. Since the Kalman filter is based upon least squares methods (LSM), it is theoretically possible to calculate the same solution as for the commonly used LSM. In the proposed Kalman filter, the large number of observations is handled epoch-byepoch, which significantly reduces the dimension of the involved matrix and accelerates the computation. Hence it outperforms the commonly used LSM in terms of efficiency. In order to avoid iterations caused by the one cycle inconsistency among FCB measurements, a pre-elimination method is developed based on the temporal stability of satellite FCBs. The pre-elimination method shows a clear advantage over post-residual adjustment, which further improves the efficiency.

A globally distributed network consisting of about 200 IGS stations has been selected to determine GPS satellite FCBs. The estimated WL FCBs have a good consistency with existing WL FCB products (e.g., CNES-GRG, WHU-SGG). The RMS of differences with respect to GRG and SGG products are 0.03 and 0.05 cycles, which indicates the consistency of the proposed approach. For satellite NL FCB estimates, 97.9% of the differences with respect to SGG products are within ± 0.1 cycles. The RMS of the differences is 0.05 cycles. These results prove the efficiency of the proposed approach.

The state-based approach of the Kalman filter also allows for more realistic modeling of stochastic parameters, which will be investigated in future research.

**Acknowledgements** We thank the IGS, GFZ and CNES for providing GNSS data and products. Thanks also go to Prof. ZHANG Xiaohong and Dr. LI Pan from WHU-SGG for FCB products and valuable discussions. Guorui XIAO is supported by the China Scholarship Council. This work was supported in part by the National Natural Science Foundation of China (grants nos. 41674016 and 41274016).

## **References**


## **Author Biographies**

**Guorui Xiao** is currently a Ph.D. candidate at Geodetic Institute, Karlsruhe Institute of Technology (KIT). He received his B.Sc. degree from the School of Geodesy and Geomatics in Wuhan University in 2011. His current research focuses on multi-frequency and multi-constellation GNSS precise positioning and applications.

**Lifen Sui** is a professor at Zhengzhou Institute of Surveying and Mapping. She obtained her Ph.D. degree from Zhengzhou Institute of Surveying and Mapping in 2001. Her main research interests include surveying data processing and GNSS positioning.

**Bernhard Heck** received the Ph.D. degree and the Habilitation in Geodesy from the University of Karlsruhe (now: KIT), Germany, in 1979 and 1984, respectively. Since 1991 he fills the position of Professor and Chair of Physical and Satellite Geodesy at KIT. His research interests include geometrical geodesy, deformation analysis, gravity field modeling, and GNSS positioning.

**Tian Zeng** is a Ph.D. student at Zhengzhou Institute of Surveying and Mapping. He obtained his B.Sc. degree in 2014. His research interests specialize in GNSS data processing and precise orbit determination.

**Yuan Tian** is a master student at Zhengzhou Institute of Surveying and Mapping. He obtained his B.Sc. degree in 2015. His research focuses on GNSS/INS integrated navigation.

# **Chapter 3 Estimating and assessing Galileo satellite fractional cycle bias for PPP ambiguity resolution**

**Guorui Xiao1,2 · Pan Li3 · Lifen Sui2 · Bernhard Heck1 · Harald Schuh<sup>3</sup>**

<sup>1</sup> Geodetic Institute, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany

<sup>2</sup> Zhengzhou Institute of Surveying and Mapping, 450052 Zhengzhou, China

<sup>3</sup> The German Research Centre for Geosciences (GFZ), D-14473 Potsdam, Germany

Submitted: 6 February 2018; Accepted: 10 October 2018; Published: 19 October 2018

GPS Solutions 23:3, 2018. DOI: 10.1007/s10291-018-0793-z

This is an author-created version of the article with permission from Springer. The final publication is available at link.springer.com

**Abstract:** Due to the rapid deployment of the Galileo constellation, Galileo is now able to contribute to GNSS precise point positioning (PPP) ambiguity resolution (AR) with 17 operational satellites as of December 2017. We estimate the satellite fractional cycle bias (FCB) based on globally distributed MGEX stations and assess the Galileo FCB quality by a comparison with that of GPS and BDS. Results of 60 days indicate that the quality of Galileo wide-lane (WL) FCB is better than GPS and BDS in terms of data usage rate, residual distribution, as well as standard deviation of daily estimates. The RMS of Galileo WL FCB residuals is 0.071 cycles, while that of GPS and BDS are 0.089 and 0.117 cycles respectively. The standard deviation of Galileo daily WL FCB is 0.010 cycles, while that of GPS and BDS are 0.018 and 0.043 cycles. We attribute the better quality of Galileo WL FCB to its signal modulation, AltBOC, which significantly compresses the multipath effect for pseudorange measurement. Within the Galileo constellation, the performance of In-Orbit Validation (IOV) satellites WL FCB is worse than that of Full Operational Capability (FOC) satellites as a result of a reduction in the power of the transmitted signal. The performance of the two highly eccentric satellites is comparable to other FOC satellites. The overall quality of Galileo narrow-lane (NL) FCB is slightly worse than that of GPS but better than that of BDS. The RMS of Galileo NL FCB residuals is 0.062 cycles, while that for GPS and BDS are 0.050 and 0.086 cycles respectively. In addition, the NL FCB quality of FOC, IOV (except E19), as well as the two eccentric satellites, shows no significant difference in terms of data usage rates and residuals. Galileo PPP AR solutions are conducted at 20 MGEX stations with three-hour sessions for ten days. The positional biases of AR solutions are 0.7, 0.6 and 2.1 cm for east, north and up components respectively, while those for float solutions are 2.1, 1.1 and 2.7 cm, corresponding to improvements of 67, 45 and 22% respectively. These results demonstrate that currently Galileo FCB can be estimated with an accuracy comparable with GPS and BDS, and the Galileo observations can bring an obvious benefit to ambiguity-fixed PPP.

**Keywords:** Galileo Precise point positioning · Integer ambiguity resolution · Fractional cycle bias Full Operational Capability

## **1 Introduction**

The fact that double differenced ambiguities in network solutions can be resolved to integer values lays the foundation of integer ambiguity resolution for precise point positioning (PPP). By estimating the fractional cycle biases (FCB) at the server and applying them to single differenced PPP at the user, the PPP integer ambiguity resolution can be achieved, convergence time significantly shortened, and the accuracy improved (Teunissen and Khodabandeh 2015). Ge et al. (2008) proposed to estimate the FCBs of the single differenced ambiguities in wide-lane (WL) and narrow-lane (NL) form by averaging the fractional parts of all involved float single-differenced WL and NL ambiguity estimates. Instead of an averaging process, a least squares method in an integrated adjustment was adopted to enhance the FCB estimates in Li et al. (2015a) and Li and Zhang (2012). Different from the above WL/NL FCB approaches, Collins et al. (2008) developed a decoupled clock model by separating satellite clocks for code and phase observations. Similarly, an integer phase clock model in which the NL FCBs were assimilated into receiver and satellite clock estimates of a network solution was developed (Laurichesse et al. 2009; Loyer et al. 2012). These PPP ambiguity resolution (AR) techniques have been proven equivalent in theory (Shi and Gao 2014), and the positional biases have been demonstrated to be minimal (Geng et al. 2010).

In addition to the development of functional models, PPP AR has been extended to multiple GNSSs. For GLONASS, special care has to be taken to deal with the receiver interfrequency bias (IFB) between satellites originating from the frequency division multiple access strategies. Usually, it is achieved by using a network of homogeneous receivers (Geng and Shi 2016; Liu et al. 2017b; Yi et al. 2017). Also, BDS-2 code measurements suffer from the satellite induced pseudorange variations which are elevation- and frequency-dependent. Correction models based on multipath combination have been proven effective (Wanninger and Beer 2015). Currently, only regional or IGSO/MEO PPP AR is achievable for BDS (Gu et al. 2015a). Tegedor et al. (2016) initially assessed the Galileo PPP AR with four In-Orbit Validation (IOV) satellites using regional stations around Europe. Li et al. (2017b) estimate the FCBs for the four global systems and assess the performance of combined PPP AR.

Compared with GPS and BDS, the study of Galileo PPP AR has been limited by the number of available satellites. Galileo is currently at the initial service phase. After the decommission of Galileo In-Orbit Validation Element-A (GIOVE-A) and GIOVE-B, four IOV spacecraft were launched between 2011 and 2012. Two Full Operational Capability (FOC) Galileo satellites were launched, however, into wrong, highly eccentric orbits in 2014 (Sonica et al. 2017). While the two satellites are unlikely to ever become a part of the operational constellation, they offer proper navigation signals and broadcast navigation messages, which allows for real-time navigation and PPP (Montenbruck et al. 2017). In 2015 another six, in May 2016 another two, and in Nov 2016 another four FOC Galileo satellites were successfully launched. The full constellation with 30 satellites is expected to be realized by 2020. As of December 2017, there are 17 active Galileo satellites as shown in Table 1. Note that E20 is excluded as it transmits only E1 signal and is not available for dual-frequency positioning. These satellites are separated into three groups: four IOV, two highly eccentric (ECC) and eleven FOC satellites. Galileo provides an open signal in the E1 band and a wideband signal covering the E5 a&b band. The Alternating Binary Offset Carrier (AltBOC) modulation can either be tracked as a composite signal or as distinct signals in the E5a and E5b sub-bands. Zaminpardaz and Teunissen (2017) compared the signal power, multipath performance, code and phase noise between IOV and FOC satellites. Results show that their characteristics differ significantly. It is necessary to investigate the impact on FCB estimation. In addition, it was found that the orbit quality of the two highly eccentric satellites is comparable with that of IOV and FOC satellites (Sonica et al. 2017; Steigenberger and Montenbruck 2017). Thus, it is also worthwhile to evaluate the performance of FCB estimation of the two eccentric satellites.


**Table 1** List of Galileo satellites as of December 2017 (https://www.gsc-europa.eu/system-status/Constellation-Information)

In the context of multi-GNSS PPP, the aim of this study is to extend the PPP ambiguity resolution to Galileo based on the globally distributed MGEX observations. Additionally, how FCB estimation can benefit from the advanced technique of Galileo signal modulation is investigated. Furthermore, within the Galileo satellite constellation, we also conduct a comparison among IOV, highly eccentric and FOC satellites. With the estimated FCBs, the improvements from Galileo PPP AR are evaluated. Finally, the methodology is summarized, and an outlook for future research is presented.

## **2 Methodology**

We start with a presentation of the basic PPP float mode, followed by a description of PPP data processing for ambiguity resolution. Then our FCB estimation strategy is introduced. With the obtained FCB products, the PPP AR algorithm at the user end is elaborated.

#### **2.1 Ambiguity-float PPP mode**

In GNSS dual-frequency PPP, the ionospheric-free (IF) combination is routinely employed to eliminate the effect of the first-order ionospheric delay. For a satellite s observed by receiver r, the corresponding pseudorange and carrier phase observation equations can be expressed as

$$\begin{cases} P\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT + D\_{r,IF} - D\_{IF}^{\mathcal{S}} + \varepsilon\_{P\_{IF}}\\ \Phi\_{r,IF}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT + \lambda\_{IF} \{N\_{r,IF}^{\mathcal{S}} + B\_{r,IF} - B\_{IF}^{\mathcal{S}}\} + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{1}$$

where indicates the geometric distance between satellite and receiver; and are the clock errors of receiver and satellite; is the slant tropospheric delay; , and are the receiver and satellite specific code hardware delays; and , are the wavelength in meter and integer ambiguity in cycle; , and are the receiver-dependent and satellite-dependent uncalibrated phase delays; , are the pseudorange and carrier phase measurement noise.

Conventionally, precise orbit and clock products from the IGS analysis center are used to remove satellite orbit and clock errors. The pseudorange ionospheric-free hardware delay bias is assimilated into the clock offset in accordance with the IGS analysis convention. Due to the fact that pseudorange measurements provide the reference to clock parameters, the actual receiver clock estimate would absorb the ionospheric-free combination of the receiver pseudorange hardware delay ,. After applying the GNSS precise satellite clock products, equation (1) can be rewritten as

$$\begin{cases} P\_{r,IF}^{s} = \rho\_r^s + dT + \overline{d}t\_r + \varepsilon\_{P\_{IF}} \\ \Phi\_{r,IF}^{s} = \rho\_r^s + dT + \overline{d}t\_r + \lambda\_{IF} \overline{N}\_{r,IF}^s + \varepsilon\_{\Phi\_{IF}} \end{cases} \tag{2}$$

where the receiver clock and ambiguity can be re-parameterized as

$$\overline{\mathbf{d}}t\_r = \left(dt\_r + D\_{r,IF}\right) \tag{3}$$

$$
\overline{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + b\_{r,IF} - b\_{IF}^{s} \tag{4}
$$

$$b\_{r,IF} = B\_{r,IF} - \frac{D\_{r,IF}}{\lambda\_{IF}} \tag{5}$$

$$
\lambda b\_{IF}^{\mathbb{S}} = B\_{IF}^{\mathbb{S}} - \frac{D\_{IF}^{\mathbb{S}}}{\lambda\_{IF}} \tag{6}
$$

Since the estimated ambiguity parameter is a combination of the integer ambiguity, the corresponding code hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite ends, the integer property is lost.

#### **2.2 FCB estimation**

In order to fix PPP ambiguity, the term , is usually decomposed into the following combination

$$
\overline{N}\_{r,IF}^{s} = \frac{\frac{cf\_2}{f\_1^2 - f\_2^2} N\_{r,WL}^{s} + \frac{c}{f\_1 + f\_2} \overline{N}\_{r,NL}^{s}}{\lambda\_{IF}} \tag{7}
$$

where , is the integer WL ambiguity and , is the float NL ambiguity. Usually, the WL ambiguity is resolved by the Hatch-Melbourne-Wübbena (HMW) combination observable (Hatch 1982; Melbourne 1985b; Wübbena 1985a). With the fixed WL ambiguity, the float NL ambiguity can be derived and tested whether it is also fixable. An ionospheric-free ambiguity is fixed only when both its WL and NL ambiguities are fixed.

The float WL ambiguity can be derived from

$$HMW = \frac{f\_1 \Phi\_{r,1}^s - f\_2 \Phi\_{r,2}^s}{f\_1 - f\_2} - \frac{f\_1 P\_{r,1}^s + f\_2 P\_{r,2}^s}{f\_1 + f\_2} = \lambda\_{r,WL}^s \overline{N}\_{r,WL}^s \tag{8}$$

where , = 1−<sup>2</sup> is the wavelength of WL ambiguity. The float WL ambiguity can be further decomposed as

$$
\overline{N}\_{r,WL}^{s} = N\_{r,WL}^{s} + d\_{r,WL} - d\_{WL}^{s} \tag{9}
$$

$$d\_{r,WL} = B\_{r,1} - B\_{r,2} - \frac{f\_1 D\_{r,1} + f\_2 D\_{r,2}}{\lambda\_{r,WL}^s (f\_1 + f\_2)} \tag{10}$$

$$d\_{WL}^{\rm s} = B\_1^{\rm s} - B\_2^{\rm s} - \frac{f\_1 D\_1^{\rm s} + f\_2 D\_2^{\rm s}}{\lambda\_{r,WL}^{\rm s} (f\_1 + f\_2)} \tag{11}$$

HMW combinations eliminate the geometric distance, satellite and receiver clock errors, as well as atmospheric delays. The key factors that affect HMW accuracy would be code noise and multipath effect. The noises of code measurements will be smoothed by the more precise phase measurements when averaging the HMW for a continuous arc.

If the WL ambiguity , is correctly resolved to an integer value and the float IF ambiguity , is obtained from PPP solution, the float NL ambiguity observable , can be derived based on (7)

$$\overline{N}\_{r,NL}^{s} = \frac{f\_1 + f\_2}{c} \lambda\_{IF} \overline{N}\_{r,IF}^{s} - \frac{f\_2}{f\_1 - f\_2} N\_{r,WL}^{s} = N\_{r,1}^{s} + d\_{r,NL} - d\_{NL}^{s} \tag{12}$$

where

$$d\_{r,NL} = \frac{f\_1 + f\_2}{c} \left(\lambda\_{IF} B\_{r,IF} - D\_{r,IF}\right) \tag{13}$$

$$\lambda d\_{NL}^{\rm s} = \frac{f\_1 + f\_2}{\mathcal{C}} \left( \lambda\_{IF} B\_{IF}^{\rm s} - D\_{IF}^{\rm s} \right) \tag{14}$$

As we can see, there are no code measurements directly involved in the calculation of NL ambiguities. Since the phase measurement is very accurate at millimeter level, the unmodeled errors in PPP will be important for the quality of NL ambiguity.

Equations (9) and (12) serve as the basic model for estimating FCBs. Since they have the same structure, a general expression can be formulated as

$$R\_r^s = \overline{N}\_r^s - \hat{N}\_r^s = d\_r - d^s \tag{15}$$

for WL and NL linear combinations, respectively. represents the FCB measurements; denotes the float undifferenced ambiguities; ̂ denotes the integer part of , which is the sum of the original integer ambiguity and the integer part of the combined code hardware delays and uncalibrated phase delays from both receiver *r* and satellite *s*; and denote the receiver and satellite FCBs. Considering the temporal stabilities of code hardware delays and uncalibrated carrier phase delays, the satellite FCBs can be estimated as constant values for specific time intervals, e.g., one day for WL FCB and 15 minutes for NL FCB (Ge et al. 2008).

A set of equations in the form of (15) can be generated based on a network of reference stations. Suppose that there are *n* satellites tracked at *m* reference stations, the system of equations can be expressed as

$$
\begin{bmatrix} R\_{r,1}^{s,1} \\ \vdots \\ R\_{r,1}^{s,n} \\ \vdots \\ R\_{r,m}^{s,1} \\ \vdots \\ R\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} \overline{N}\_{r,1}^{s,1} - \delta \underline{\mathcal{I}}\_{r,1}^{s,1} \\ \vdots \\ \overline{N}\_{r,1}^{s,n} - \delta \underline{\mathcal{I}}\_{r,1}^{s,n} \\ \vdots \\ \overline{N}\_{r,m}^{s,1} - \delta \underline{\mathcal{I}}\_{r,m}^{s,1} \\ \vdots \\ \overline{N}\_{r,m}^{s,n} - \delta \underline{\mathcal{I}}\_{r,m}^{s,n} \end{bmatrix} = \begin{bmatrix} 1 & \cdots & 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cdots & 0 & 0 & \cdots & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & \cdots & -1 \end{bmatrix} \begin{bmatrix} d\_{r,1} \\ \vdots \\ d\_{r,m} \\ d\_{s,1}^{s} \\ \vdots \\ d\_{s,n}^{s} \end{bmatrix} \tag{16}$$

The obtained system of equations is singular. Assuming the float ambiguities ̅ are precisely estimated, ̂ can be determined by rounding ̅ and one arbitrarily selected FCB should be set to zero. In this way, the system of equations can be resolved.

#### **2.3 PPP AR at the user end**

With the above-estimated FCBs, the satellite FCBs can be corrected, and the single-differenced PPP ambiguities are aimed to be fixed. It is noted that, although single differencing operations could be conducted on raw observation equations and float ambiguities level, single differencing at the ambiguity level may be better in practice as it provides more flexibility of the choice of a reference satellite and the formation of satellite pairs. Therefore, we recommend performing the single differencing at ambiguity level,

$$
\overline{N}\_r^s - \overline{N}\_r^m = N\_r^{s,m} + d^{s,m} \tag{17}
$$

For WL float ambiguities, they can be directly fixed by the rounding approach (Ge et al. 2008). For NL float ambiguities, they are fed into the LAMBDA algorithm (Teunissen et al. 1997). If not all the float NL ambiguities can be fixed by the LAMBDA method, partial ambiguity resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999).

For sequential PPP AR processing realized with Kalman filter, such as real-time applications, there exists a discrepancy concerning whether the states of the fixed solution should be used for the time update of the subsequent epoch. The float and fixed solution loops are independent. Even though the fixed solution is obtained, the states of the floating loop do not change accordingly (Wang et al. 2017). We argue that the benefits of ambiguity resolution are not fully exploited in this approach. On the other hand, the fixed ambiguity of the current epoch should be tightly constrained to the succeeding epoch in a continuous arc (Takasu and Yasuda 2010). This approach utilizes all the states of the fixed solution to the next time update of the float solution and is introduced for PPP AR here. Once the WL and NL integer ambiguities are resolved and validated, a tight constraint can be reconstructed and imposed on the ionospheric-free float ambiguities

$$0 = \overline{N}\_r^s - \overline{N}\_r^m - N\_r^s + N\_r^m - d^{s,m},\ \sigma\_d \tag{18}$$

In this way, a fixed solution can be obtained at the user. It is noted that only the NL FCB is used in the reconstructed fixed ambiguity and has direct influence on the estimated parameters. In this sense, the NL FCB has to be estimated as precisely as possible.

## **3 Data and processing strategy**

The International GNSS Service (IGS) established the Multi-GNSS Experiment (MGEX) in order to prepare operational services for new and upcoming GNSS, such as the European Galileo, Chinese BeiDou, and regional systems, such as Japanese QZSS and Indian NAVIC system (Montenbruck et al. 2017). Currently, there are over 210 MGEX stations. Figure 1 shows the geographic distribution of all the MGEX stations. Among these stations, 208 stations support Galileo signal tracking, and 165 stations support BDS signal tracking. The MGEX stations, equipped with multiple brands of professional receivers, provide almost full and continuous tracking of Galileo dual-frequency signals from the 17 satellites. Therefore, all the data collected by the MGEX network from June 4 to August 2, 2017, in total 60 days, are processed.

**Figure 1** Geographic distribution of all the MGEX stations. Almost all stations support Galileo, while the number of BDS stations, represented by the red circles, is slightly smaller.

In the processing, E1/E5a are used for Galileo, while L1/L2 and B1/B2 are used for GPS and BDS respectively. These signals were chosen according to a high availability, and in accordance with the principle of orbit and clock generation (Uhlemann et al. 2016). In order to obtain accurate float ambiguities, a GPS/Galileo/BDS combined PPP mode is used at the server end. Inter-system biases are estimated per station and per system. The combination of multiple GNSSs increases the data usage rates of the incomplete Galileo and BDS systems, as the number of observable satellites for either of the two systems may be not enough for a standalone PPP solution (Li et al. 2017a). Figure 2 shows the global distribution of the number of observable Galileo satellites with a cut-off elevation of 10 degrees. It can be seen that a standalone Galileo PPP solution is not feasible in a few regions. It is noted that this phenomenon is temporary, and significant improvements can be expected when the constellation is completed. In addition, multi-GNSS PPP is robust against outliers and improves the precision of the float ambiguities. A forward and backward filter is used to avoid the convergences of ambiguity parameters. For the sake of comparison, the basic settings of data edit are the same for all three systems: the cut-off elevation angle is set to 10 degrees; float ambiguities with an elevation below 30 degrees or with standard deviation (STD) larger than 0.1m are removed. It is noted that the BDS satellite-induced code multipath effects are corrected for IGSO and MEO satellites according to Wanninger and Beer (2015), while GEO satellites are excluded from the processing. The third generation BDS satellite, which are no longer affected by such effects (Lei et al. 2017), are also excluded since the data is not publicly available.

**Figure 2** Distribution of the number of observable Galileo satellites at DoY 160, 2017 with a cut-off elevation of 10 degrees

Throughout the processing, MGEX precise products, including precise satellite orbit, clock, and earth rotation parameters, provided by GFZ (Uhlemann et al. 2016) are used. The satellite phase center offsets and variations are corrected according to the IGS antenna file (igs14\_1949.atx). As for the receiver antenna phase center offsets and variations, the correction values for GPS are employed for both Galileo and BDS in accordance with the principle of orbit and clock generation (Prange et al. 2017). The detail of our FCB estimating strategy could be found in Xiao et al. (2018).

## **4 Results, comparisons, and discussions**

In this section, the quality of Galileo WL/NL FCB estimation is thoroughly analyzed. We further compared the FCB estimations with that of GPS and BDS in order to characterize the Galileo signal properties. The generated FCBs are used for static and kinematic PPP AR solutions.

#### **4.1 Satellite WL FCB estimation**

Figure 3 presents the time series of Galileo satellite WL FCBs during the 60 days. For better visualization, the derived satellite FCBs are shifted with integer cycles. It can be seen that the derived Galileo satellite WL FCBs are quite stable over time, with a STD of 0.010 cycles. It would be sufficient to broadcast daily WL FCB, which is similar to other GNSSs (Ge et al. 2008). For comparison, we also present the STD statistics for GPS and BDS in Figure 4. Note that BDS C07 and C11 are excluded from the statistics as there exist jumps in the time series due to unknown reasons, which has been reported by Li et al. (2017a) and Li et al. (2017b). The average STDs are calculated for GPS and BDS, while an individual STD is calculated for each Galileo satellite. The STDs of WL FCB are 0.018 cycles for GPS and 0.041 cycles for BDS. This implies that Galileo WL FCB is less noisy than GPS and BDS. We attribute the better performance of Galileo WL FCBs to its signal modulation. The advanced AltBOC modulation can significantly compress the multipath effect for pseudorange measurements and decrease the fluctuations of HMW combinations, especially for low elevations (see Figure 5). Also, we notice that there is a difference between FOC and IOV satellites. IOV satellites show the largest STDs with 0.013 cycles, while that of the two eccentric satellites is at the same level as other FOC satellites with 0.010 cycles. To further investigate the possible cause, more analyses including data usage rates, residual statistics are performed.

**Figure 3** Time series of Galileo satellite WL FCBs for all the 17 satellites while E30 is selected as the reference

**Figure 4** STDs of daily satellite WL FCBs with 60 days. Averaged STDs are calculated for GPS and BDS, while an individual STD is calculated for each Galileo satellite. E30 is selected as the reference satellite for Galileo, while G01 and C09 are selected for GPS and BDS separately

**Figure 5** HMW linear combinations of GPS, Galileo, and BDS dual frequency code measurements at station YEL2 on DoY 160, 2017

**Figure 6** Usage rates of the float WL ambiguities. Averaged usage rates are calculated for GPS and BDS, while an individual usage rate is calculated for each Galileo satellite

The float ambiguity observations, with residuals larger than 0.3 cycles, are rejected in FCB estimation. Figure 6 shows the usage rates of WL ambiguities for each satellite. We can see that almost all of the WL ambiguities are used for WL FCB estimation. The averaged usage rate is 96.2% for Galileo, while that for GPS and BDS were 95.3% and 90.1%, respectively. This indicates that the number of code multipath events for Galileo is smaller than for GPS, which may be ascribed to the advanced signal modulation of E5a. The usage rates of the two ECC satellites are comparable to other FOC satellites, while that for the IOV satellites are slightly lower. This indicates that the quality of FOC WL float ambiguities is slightly better than that of IOV satellites. The averaged data usage for BDS is about 90% which is much lower than that of Galileo and GPS.

**Figure 7** Distribution of posterior residuals for Galileo (left), GPS (middle) and BDS (right) after WL FCB estimation

The quality of FCB estimation can also be indicated by the posterior residuals of FCB estimation. Figure 7 shows the distribution of the residuals of WL FCB estimations. The mean values are -0.001, -0.003 and -0.003 cycles for Galileo, GPS, and BDS, respectively. The percentages of residuals smaller than 0.1 cycles are 86.5, 77.3 and 61.5% for Galileo, GPS and BDS, respectively. The RMS of Galileo WL residuals is 0.071 cycles, while that for GPS and BDS are 0.089 and 0.117 cycles, respectively. The smaller and better distribution of Galileo residuals indicates a better quality of WL ambiguities. This further confirms that the code multipath effects of Galileo are remarkably reduced.

In order to further compare the WL FCB quality of IOV and FOC satellites, the RMS for each satellite in Galileo are presented in Figure 8. The figure also presents linear averages of RMS for each GPS and BDS satellite. It is found that the RMS of IOV satellites are significantly larger than other FOC satellites, while there is no significant difference between FOC and ECC satellites. The average RMS of IOV satellites is 0.089 cycles, while that for FOC satellites is only 0.066 cycles. The percentage of residuals smaller than 0.1 cycles is 76.4% and 88.9% for IOV and FOC satellites, respectively. These results indicate that the quality of FOC WL float ambiguities is better than that of IOV satellites. The worse performance of IOV satellites probably stems from the reduction of signal transmission power for IOV satellites. ESA imposed a reduction of 1.5 dB in the signal power of all four IOV satellites following a payload power problem of the fourth IOV (E20) in 2014. The 1.5 dB power decrease results in an increase of approximately 15-20% of the thermal noise of the receiver, which roughly matches the observed increase in IOV WL FCB residuals.

**Figure 8** RMS of WL residuals after WL FCB estimations. Averaged RMSs are calculated for GPS and BDS, while an individual RMS is calculated for each Galileo satellite

From the above results, we can safely conclude that the overall performance of WL FCB can be expressed by the following inequality:

FOC ≈ ECC > IOV ≈ GPS > BDS

The WL FCB quality of the two ECC satellites is comparable to that of FOC satellites. Despite fewer observations acquired by receivers, the two satellites can be used for WL ambiguity resolution. Concerning the lower data usage rates, larger residuals and larger STDs of WL FCB estimations for IOV satellites, the quality of IOV WL FCB is worse than that of FOC satellites. Despite the difference between the Galileo IOV and FOC satellites, the overall performance of Galileo WL FCB is better than that of GPS and BDS.

### **4.2 Satellite NL FCB estimation**

Figure 9 presents the time series of Galileo/GPS/BDS satellite NL FCBs on DoY 160 in 2017. For better visualization, the derived satellite FCBs are shifted with integer cycles. It can be seen that the derived satellite Galileo NL FCBs are only stable in tens of minutes. It is noted that the common mode shape of NL FCB variations is caused by the reference satellite as epoch-wise differencing is applied.

**Figure 9** Time series of satellite NL FCBs for Galileo (upper), GPS (middle) and BDS (bottom), while satellite PRN E30, G10, C10 are selected as the references respectively. Each line represents one satellite, while the lines for reference satellites are zero and excluded.

Figure 10 depicts the usage rate of float ambiguities for each Galileo satellite (average 94.0%), as well as the average usage rates for GPS (97.6%) and BDS (80.4%). The usage rate of GPS NL float ambiguities is highest while that of BDS is worst. Since the weights of carrier phase measurements are 10000 times larger than that of code measurements, the quality of the float ambiguities is dominated by the unmodeled errors in PPP. Therefore, the highest usage of GPS is reasonable since its precise product and models are currently state of the art. All Galileo satellites except E19 show similar usage rates ranging from 91.8 to 96.9%, indicating a similar quality level of NL ambiguities. The usage rate of E19 is only 83.8%, which is much smaller than that for other Galileo satellites. Inspired by Zaminpardaz and Teunissen (2017) who reported that for SEPTENTRIO receivers the carrier-to-noise density ratio for E19 lies below the value of the other two IOV satellites for elevations higher than 60 degrees, we further calculated the carrier-to-noise density ratio, data usage rates for three receiver groups, namely LEICA (37 stations), SEPTENTRIO (41 stations) and TRIMBLE (83 stations). From the results in Figure 11, it can be seen that the carrier-to-noise density ratio of E19 lies below the other satellites for all the receiver groups. In addition, the data usage rates of E19 are lower than for the other satellites for all the three groups and no clear relationship with receiver type is found. Furthermore, we notice that the STD of satellite laser ranging (SLR) residuals for E19 (0.130 m) is larger than for E11 (0.092 m) and E12 (0.086 m) when performing SLR validation (Guo et al. 2017). This may indicate poor quality of precise satellite products for E19. Based on the above analysis, we suspect that the lower data usage rate of E19 is caused by the joint effect of poorer satellite products and larger carrier phase noise.

**Figure 10** Usage rates of float NL ambiguities for satellite NL FCB estimation. Averaged usage rates are calculated for GPS and BDS, while an individual usage rate is calculated for each Galileo satellite

**Figure 11** Carrier-to-noise density ratio C/N0 of Galileo signals for three receiver groups, TRIMBLE (left), SEPTENTRIO (middle) and LEICA (right)

**Figure 12** Distributions of posterior residuals for Galileo (left), GPS (middle) and BDS (right) after NL FCB estimation

Figure 12 shows the distribution of the residuals of NL FCB estimations for the 60 days. The mean values are 0.000, 0.000 and -0.001 cycles for GPS, Galileo, and BDS, respectively. The percentages of residuals that are smaller than 0.1 cycles are 94.6, 91.4 and 83.6% for GPS, Galileo, and BDS, respectively. The RMS of Galileo is 0.062, while that for GPS and BDS are 0.050 and 0.086 cycles, respectively. These results indicate that the quality of Galileo NL FCB estimations is slightly worse than that of GPS but better than that of BDS. As discussed above, the quality of NL FCB estimation could be affected by unmodeled errors. The lack of precise Galileo antenna PCO and PCV models both at satellite and user end certainly has influences on the PPP float ambiguities. The worst performance of BDS can be attributed to the poor quality of precise satellite products (Montenbruck et al. 2017), as well as the satellite-induced code multipath effects.

**Figure 13** RMS of NL residuals after NL FCB estimations. Averaged RMSs are calculated for GPS and BDS, while an individual RMS is calculated for each Galileo satellite

In order to compare the NL FCB quality of IOV and FOC satellites, the RMSs for each Galileo satellite are presented in Figure 13. It can be found that there is no significant difference among FOC, IOV and ECC satellites. The average RMSs are 0.061, 0.065 and 0.058 cycles for FOC, IOV, and ECC, respectively. The percentages of residuals that are smaller than 0.1 cycles are 91.9, 90.3 and 92.7% for FOC, IOV, and ECC, respectively. These results indicate that the satellites in the incorrect orbital planes behave not worse than those in nominal FOC and IOV orbits in terms of PPP. We can further deduce that the orbit and clock quality of ECC satellites from MGEX precise product are at a similar level as that of FOC and IOV. The largest RMS is found for satellite E19 with values at about 0.069 cycles. Concerning the lower data usage of this satellite, the quality of its NL FCB is the worst. From the results of the usage rates of NL float ambiguities and distributions of the residuals, we conclude that the overall performance of Galileo NL FCB is worse than that of GPS but better than BDS. It is noted that the ranking is temporary, and promising improvements can be expected for Galileo and BDS considering their rapid developments.

### **4.3 PPP AR solution**

In order to validate our FCB estimates, as well as to assess the performance of PPP AR, a subset of MGEX network stations is processed in PPP AR mode. These test stations have been excluded from the FCB estimation. Due to the incomplete constellation of Galileo, the number of observable satellites differs with respect to geographical locations. 20 MGEX stations, as shown in Figure 14, with relatively higher percentage of Galileo observations are selected. GPS and Galileo observations from DoY 155 to 164 are processed. For static solutions, the 24 hours observations are divided into eight three-hour sessions. For each session, the epochs with less than 4 visible Galileo satellites are deleted. Then

the sessions with data integrity less than 80 percent are removed. After the pre-processing, there are 1484 sessions. For pseudo-kinematic solutions, the 24 hours observations are regarded as one session. Galileo-only solutions, GPS-only solutions, as well as GPS/Galileo-combined solutions, are calculated for each session. The positional biases of static Galileo PPP float solutions and AR solutions on DoY 160, 2017 are presented in Figure 15. The statistics of all sessions in the ten days are provided in Table 2.

**Figure 14** Geographical distribution of the selected 20 MGEX stations

**Figure 15** Absolute NEU biases of static Galileo PPP float and AR solutions with three-hour sessions on DoY 160, 2017. In total, there are 98 sessions for the 20 MGEX stations


**Table 2** Statistics of static PPP float and AR solutions with three-hour sessions for 20 stations in ten days [Unit: m]

It can be seen that the accuracy of the Galileo-only PPP float solution is several centimeters with the current constellation. After the ambiguity resolutions, the accuracy is significantly improved. The most significant improvement is found for the east component, from 2.1 to 0.7 cm, corresponding to an enhancement of 67%. For GPS, the improvement of the north component by PPP AR is insignificant as the accuracy of the north component in the float solution is relatively high. However, in Galileo, we also observe an obvious improvement for the north component, from 1.1 to 0.6 cm with an enhancement of 45%. This is because the accuracy of the north component in Galileo float solution is low due to bad geometry. The improvement for the vertical component is about 22%, which is less significant than for the horizontal components. The possible reason is that the vertical component is tightly coupled with receiver clock and zenith wet delay parameters, which are simultaneously estimated.

Compared with the GPS-only solutions, both float and AR solutions of Galileo are slightly worse, although it has been proven that the quality of Galileo WL FCB is better than GPS. This can partly be ascribed to the smaller number of available satellites. It can also be ascribed to the worse quality of Galileo NL FCB as the accuracy of NL FCB directly affects the accuracy of PPP AR.

**Figure 16** NEU biases of Galileo, GPS, GPS+Galileo PPP AR in static mode, at station YEL2 on DoY 160, 2017

Figure 16 shows a typical example for standalone single- and combined PPP AR solutions in static mode. One can see that the Galileo-only solution requires more time to achieve a reliable accuracy, while the convergence time of GPS-only solutions is much shorter. The fastest convergence can be achieved in GPS/Galileo PPP AR. Although the final accuracy of the GPS/Galileo combined solution is almost the same as that of the GPS-only solution, integrating Galileo with GPS increases the fixing rate and results in a more stable solution. This is further confirmed by the reduction of incorrectly fixed sessions (see Table 2), where the success rate increases from 95.8% to 97.2% when adding Galileo to GPS.

**Figure 17** NEU biases of kinematic Galileo PPP float and AR solutions at station YEL2 on DoY 160, 2017

Figure 17 presents the kinematic positional biases for station YEL2 on DoY 160, 2017. Epochs with less than five valid satellites have been deleted. It can be seen that the accuracy of all three components is significantly improved by ambiguity resolution. In addition, it is observed that the convergence time is also shortened. And yet the time to first fix (TTFF) is still larger than two hours as a result of the limited number of observable satellites.

**Figure 18** NEU biases of kinematic GPS, and GPS/Galileo PPP AR at station YEL2 on DoY 160, 2017

Although Galileo-only solutions are limited by the number of observable satellites, adding Galileo to GPS can enhance the performance of kinematic PPP solutions. Figure 18 shows a typical example for GPS-only and GPS/Galileo combined PPP AR solutions. As we can see, integrating Galileo with GPS increases the accuracy of kinematic AR solutions, especially when the accuracy of the GPS-only solution is bad. The combined solution can achieve an accuracy of 1-2 cm in horizontal and 4-6 cm in vertical components.

**Table 3** Statistics of kinematic PPP float and AR solutions for 20 stations in ten days [Unit: m]


Pseudo-kinematic PPP AR solutions are conducted for each station in the ten days. The statistics of positional biases are presented in Table 3. For Galileo-only solutions, the RMSs of the east, north and up components are improved by 59%, 45%, and 28% respectively, while for GPS-only solutions the enhancements are 50%, 31% and 14%. Similar to the static results, the improvements for Galileo are more significant than for GPS. Dualsystem PPP AR can further improve the kinematic positional accuracy.

## **5 Conclusions**

The feasibility of Galileo PPP ambiguity resolution with the current constellation has been demonstrated. The satellite WL/NL FCBs are estimated from globally distributed MGEX stations. Results indicate that the quality of Galileo WL FCB is better than for GPS and BDS in terms of data usage rate, residual distribution, as well as STD of daily estimates. We attribute the good quality of Galileo WL FCB to its advanced signal modulation, AltBOC, which significantly compresses the multipath effect for code measurements. Within the Galileo constellation, the quality of FOC WL FCB is much better than for IOV satellites. The poorer performance of IOV satellites WL FCB is a result of a reduction in the satellite transmit signal power. The performance of the two satellites with highly eccentric orbits is comparable to other FOC satellites but having a smaller number of observations. As for NL FCB, the quality of Galileo NL FCB is slightly worse than that of GPS but better than that of BDS. Since the accuracy of NL FCB estimation is dominated by unmodeled errors in float PPP, the worse quality of Galileo NL FCB is likely caused by the imperfect PCO and PCV models used in this study. Within Galileo, the NL FCB quality of FOC and IOV (except E19), as well as the two eccentric satellites, shows no significant difference in terms of data usage rates and residuals. The reason for the worse performance of E19 is still not clear. On the one hand, it cannot, or at least not fully, be ascribed to the signal transmission power as the power difference between FOC and IOV E11/E12 is larger than the difference between E11/E12 and E19. On the other hand, the worse quality of E19 satellite orbit, indicated by SLR residuals, could also degrade NL FCB estimations. This issue remains an open question and deserves further investigation. For static Galileo PPP AR with three-hour sessions, the positional biases can be reduced by 67, 45 and 22% for east, north and up components, respectively. In addition, PPP AR also improves the kinematic solutions by 59, 45 and 28% for east, north and vertical components, respectively. Although the performance of Galileo PPP AR is still worse than that of GPS, a promising improvement can be expected in the near future as the Galileo FOC satellite metadata, including satellite mass, attitude law, PCO and PCV, are published on 2017-10-06 (https://www.gsc-europa.eu/support-to-developers/galileosatellite-metadata). This informa-tion is expected to improve the precise orbit determination and PPP solution.

**Acknowledgments** Thanks go to IGS-MGEX and GFZ for providing GNSS data and products. Thanks also go to Maorong GE from GFZ for valuable discussions. Guorui Xiao is supported by the China Scholarship Council. This work was supported in part by the National Natural Science Foundation of China (grants nos. 41674016, 41274016 and 41604024). Finally, we highly acknowledge the comments by two anonymous reviewers which helped to improve the submitted version of the paper.

## **References**


## **Author Biographies**

**Guorui Xiao** is currently a Ph.D. candidate at Geodetic Institute, Karlsruhe Institute of Technology (KIT). He received his B.Sc. degree from the School of Geodesy and Geomatics in Wuhan University in 2011. His current research focuses on multi-frequency and multi-constellation GNSS precise positioning and applications.

**Pan Li** is a postdoctoral researcher at GFZ, Germany. He obtained his Ph.D. degree in 2016 from School of Geodesy and Geomatics, Wuhan University, China. His current research focuses mainly on GNSS precise point positioning and ambiguity resolution.

**Lifen Sui** is a professor at Zhengzhou Institute of Surveying and Mapping. She obtained her Ph.D. degree from Zhengzhou Institute of Surveying and Mapping in 2001. Her main research interests include surveying data processing and GNSS positioning.

**Bernhard Heck** received the Ph.D. degree and the Habilitation in Geodesy from the University of Karlsruhe (now: KIT), Germany, in 1979 and 1984, respectively. Since 1991 he fills the position of Professor and Chair of Physical and Satellite Geodesy at KIT. His research interests include geometrical geodesy, deformation analysis, gravity field modelling, and GNSS positioning.

**Harald Schuh** obtained his Ph.D. from the University of Bonn, Germany and is currently the Director of the Department 1 "Geodesy", German Research Centre for Geosciences (GFZ). He is also the President of the International Association of Geodesy (IAG) since 2015. His research interests are the space geodetic techniques in particular Very Long Baseline Interferometry and GNSS and their applications.

# **Chapter 4 A Unified Model for Multi-Frequency PPP Ambiguity Resolution**

#### **Guorui Xiao1,2,\* · Pan Li3,\* · Yang Gao4 · Bernhard Heck<sup>1</sup>**

<sup>1</sup> Geodetic Institute, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany; bernhard.heck@kit.edu

<sup>2</sup> Zhengzhou Institute of Surveying and Mapping, Zhengzhou 450052, China

<sup>3</sup> The German Research Centre for Geosciences (GFZ), D-14473 Potsdam, Germany

<sup>4</sup> Department of Geomatics Engineering, University of Calgary, Calgary, AB T2N 1N4, Canada; ygao@ucalgary.ca

**\*** Correspondence: xgr@whu.edu.cn (G.X.); panli@gfz-potsdam.de (P.L.)

Submitted: 11 December 2018; Accepted: 7 January 2019; Published: 10 January 2019

Remote Sensing 11(2):116, 2019. DOI:10.3390/rs11020116

This is an author-created version of the article with permission from MDPI. The final publication is available at www.mdpi.com

**Abstract:** With the modernization of Global Navigation Satellite System (GNSS), triple- or multifrequency signals have become available from more and more GNSS satellites. The additional signals are expected to enhance the performance of precise point positioning (PPP) with ambiguity resolution (AR). To deal with the additional signals, we propose a unified modeling strategy for multi-frequency PPP AR based on raw uncombined observations. Based on the unified model, the fractional cycle biases (FCBs) generated from multi-frequency observations can be flexibly used, such as for dual- or triple- frequency PPP AR. Its efficiency is verified with Galileo and BeiDou triple-frequency observations collected from globally distributed MGEX stations. The estimated FCB are assessed with respect to residual distributions and standard deviations. The obtained results indicate good consistency between the input float ambiguities and the generated FCBs. To assess the performance of the triple-frequency PPP AR, 11 days of MGEX data are processed in three-hour sessions. The positional biases in the ambiguity-fixed solutions are significantly reduced compared with the float solutions. The improvements are 49.2%, 38.3%, and 29.6%, respectively, in east/north/up components for positioning with BDS, while the corresponding improvements are 60.0%, 29.0%, and 21.1% for positioning with Galileo. These results confirm the efficiency of the proposed approach, and that the triple-frequency PPP AR can bring an obvious benefit to the ambiguity-float PPP solution.

**Keywords:** Galileo · BeiDou · Precise point positioning · Integer ambiguity resolution · Triplefrequency · Fractional cycle bias

## **1 Introduction**

Precise point positioning (PPP) has found increased applications due to its cost-effectiveness, global coverage, and high accuracy (Kouba and Héroux 2001; Zumberge et al. 1997). Usually PPP is able to achieve a positional accuracy of 10 cm after a convergence time of 30 min (Cai et al. 2015). The integer ambiguity resolution (AR) technique is expected to further enhance the accuracy and shorten the convergence time (Gabor and Nerem 1999; Gao and Shen 2002). In addition, (Katsigianni et al. 2018) shows that the Galileo orbit determination could be improved when employing AR in multiple Global Navigation Satellite System (GNSS) data processing. However, the uncalibrated phase delays (UPDs) originating from satellites and receivers destroy the integer nature of PPP ambiguities. By determining the UPDs, to be estimated as fractional cycle biases (FCBs) at the server end and applying them at the user end, the PPP integer ambiguity resolution could become feasible (Ge et al. 2008; Li et al. 2015a; Xiao et al. 2018b). Similarly, the decoupled clock model (Collins et al. 2008) and the integer phase clock model (Laurichesse et al. 2009) were developed. These PPP AR techniques have been proven equivalent in theory (Shi and Gao 2014; Teunissen and Khodabandeh 2015), and the positional biases have been demonstrated to be minimal (Geng et al. 2010). Beside GPS, PPP AR has been extended to GLONASS (Geng and Shi 2016; Liu et al. 2017b; Yi et al. 2017), BeiDou Navigation Satellite System (BDS) (Li et al. 2017a; Liu et al. 2017a), Galileo (Tegedor et al. 2016; Xiao et al. 2019b), and multi-GNSS (Li et al. 2017b).

As to the functional models used for PPP, the dual-frequency ionospheric-free (IF) combination is routinely employed (e.g., in the above-mentioned research). However, with emerging BDS and Galileo, as well as the modernization of GPS and GLONASS, various types of multi-frequency observables become available (Montenbruck et al. 2017). The choice of optimum combinations then becomes practically difficult given the diversity of equipment (Schönemann et al. 2011). In addition, the IF combination will amplify the measurement noise level by a factor of about 3, which will degrade the performance of the position solution. As a result, the PPP model based on uncombined measurements, in which the individual signal of each frequency is treated as an independent observable, has drawn increasing interest in the GNSS community (Gu et al. 2015b). Its efficiency has already been confirmed in terms of convergence time and precision for single-frequency PPP (Lou et al. 2016), multi-GNSS PPP (Chen et al. 2015), as well as PPP-RTK (Feng et al. 2013b; Odijk et al. 2016b). Moreover, this approach has been tested effectively for ionospheric modelling (Tu et al. 2013; Zhang et al. 2012), differential code bias (DCB) estimation (Liu et al. 2018; Shi et al. 2015), and low earth satellite orbit determination (Zehentner and Mayer-Gürr 2016).

Compared to the well-developed IF PPP model, the uncombined PPP model, called uncombined PPP in the sequel, requires more investigation, especially in the case of multifrequency processing and ambiguity resolution. First, how to model and constrain the ionospheric delay has a crucial impact on performance using the uncombined PPP. For example, it has been demonstrated that a white noise model is not adequate to capture the characteristics of the ionospheric delay. The external constraints developed from the ionospheric products, such as the IGS global ionosphere maps, are also not accurate enough to completely separate the ionospheric effects from the ambiguity parameters (Gu et al. 2015b). The influence of the ionospheric effects on the ambiguity fixing therefore must be reduced. Second, the method to deal with the DCB errors is more problematic with the uncombined PPP (Guo et al. 2015) than the dual-frequency IF PPP, since the latter can cancel out the DCB biases. The problem of partial assimilation of the code bias (DCB) into phase bias (FCB) should also be carefully considered. Third, the uncombined PPP approach was proposed to deal with multi-GNSS and multi-frequency signals, so a generalized FCB estimation and AR method (Li et al. 2018), which is extendable to dual-, triple-, and multi- frequency, should be proposed.

Li et al. (2013a) verified the feasibility of the uncombined PPP AR with refined ionospheric models. The ionospheric delay was constrained from a priori spatial-temporal information and ionospheric products. The GPS dual-frequency ambiguities were fixed sequentially in the forms of wide-lane (WL)/narrow-lane (NL), which followed the convention of IF PPP AR. Gu et al. (2015a) further testified the uncombined PPP AR with BDS triple-frequency observations. The extra-wide-lane (EWL) and WL ambiguities were successfully fixed, whereas the B1 ambiguities were kept as float values. In addition, the performance was further limited by the satellite-induced multipath effects (Wanninger and Beer 2015). Li et al. (2018) proposed a unified FCB estimation and PPP AR method, which is extendable to multi-frequency uncombined PPP. The FCBs on each frequency were directly estimated from the raw float ambiguities derived from triple frequency observables. The model showed a great potential for multi-frequency uncombined PPP AR, although its DCB strategy may not be optimal. The satellite DCBs, together with the receiver DCB, were estimated as unknowns, and as a result the number of unknown parameters was increased. Given that the satellite GNSS DCB product is currently available on a routine basis (Wang et al. 2016), it would be beneficial to make use of these products. In addition, validating the method with Galileo observations is of interest considering the recent and rapid development of Galileo.

The aim of this study is to develop a unified modeling strategy for multi-frequency and multi-GNSS uncombined PPP AR. The unified model is able to generate consistent FCB products and perform PPP AR for multi-frequency PPP. The proposed approach will be first described, and its effectiveness will be verified with Galileo and BDS triple-frequency observations collected from the globally distributed Multi-GNSS EXperiment (MGEX) stations. The estimated FCB are assessed with respect to residual distributions and standard deviations, followed by an evaluation of the performance improvements in Galileo and BDS triple-frequency PPP AR. Finally, the results are summarized, and an outlook for future research is presented.

## **2 Methodology**

The proposed uncombined PPP mode will be first described in this section, followed by a description of our FCB estimation strategy. With the obtained FCB products, the uncombined PPP AR algorithm at the user end is then elaborated.

#### **2.1 Uncombined PPP Float Ambiguity Model**

In the classic GNSS dual-frequency PPP, the first-order ionospheric delay is eliminated by the formation of the IF combination (Zumberge et al. 1997). In the uncombined PPP model, the ionospheric delay is directly estimated. For a satellite s observed by receiver r, the corresponding raw pseudo-range and carrier phase observation equations can be expressed as (Leick et al. 2015)

$$\begin{cases} P\_{r,f}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT + a\_f \cdot dl\_{r,1}^{\mathcal{S}} + D\_{r,f} - D\_f^{\mathcal{S}} + \varepsilon\_{P\_f} \\ \Phi\_{r,f}^{\mathcal{S}} = \rho\_r^{\mathcal{S}} + dt\_r - dt^{\mathcal{S}} + dT - a\_f \cdot dl\_{r,1}^{\mathcal{S}} + \lambda\_f \{N\_{r,f}^{\mathcal{S}} + B\_{r,f} - B\_f^{\mathcal{S}}\} + \varepsilon\_{\Phi\_f} \end{cases} \tag{1}$$

where the subscript = (1,2,3, ⋯ ) refers to a specific carrier frequency, superscript s refers to a specific satellite; indicates the geometric distance between the satellite and receiver; and are the clock errors of receiver and satellite; is the slant tropospheric delay; ,1 is the slant ionospheric delay on the first carrier frequency and = 2 1 <sup>2</sup> ⁄ is the carrier frequency-dependent factor; , and are the receiver and satellite specific code hardware delays; and , are the wavelength in meter and integer ambiguity in cycle; , and are the receiver-dependent and satellite-dependent uncalibrated phase delays; and are the pseudo-range and carrier phase measurement noise, respectively. Note that the higher-order ionospheric effects are neglected, as they have limited influence on the performance of ambiguity resolution (Hadas et al. 2017).

Another important difference between the IF PPP and the uncombined PPP model is the strategy to deal with the DCB. The DCB is not of concern in IF PPP as the IF combination is also used for precise clock generation, which implies that the DCB could be fully absorbed by other parameters or simply cancelled out in the IF PPP (Dach et al. 2015). But this is not the case in the uncombined PPP, especially with multi-frequency observations. Conventionally, precise orbit and clock products from the IGS analysis center are used to remove satellite orbit and clock errors. During the generation of precise clock products, the pseudo-range IF hardware delay bias = 2 2−1 1 − 1 2−1 2 is assimilated into the clock offset in accordance with the IGS analysis convention. After applying the GNSS precise satellite clock products, Equation (1) can be rewritten as

$$\begin{cases} P\_{r,f}^{s} = \rho\_r^s + dt\_r - dt\_{pre}^s + dT + a\_f \cdot dl\_{r,1}^s + D\_{r,f} - D\_f^s + D\_{IF}^s + \varepsilon\_{P\_f} \\ \Phi\_{r,f}^{s} = \rho\_r^s + dt\_r - dt\_{pre}^s + dT - a\_f \cdot dl\_{r,1}^s + \lambda\_f \left( N\_{r,f}^s + B\_{r,f} - B\_f^s \right) + D\_{IF}^s + \varepsilon\_{\Phi\_f} \end{cases} \tag{2}$$

This linear system is rank-deficient due to the DCB parameters. For dual-frequency uncombined PPP processing, the singularities can be eliminated by a re-parameterization process. This is accomplished based on the fact that the DCB parameters can be separated into satellite-related, receiver-related, and frequency-related parts, and therefore can be fully absorbed, respectively, by satellite clock, receiver clock, and ionospheric parameters (Zhang et al. 2012). This method is efficient for dual-frequency processing but becomes complicated when facing multi-frequency observations. It is also possible to estimate these DCB parameters in advance. This is typically done by employing a network of receivers and imposing a zero-mean constraint. This option is complicated and not suitable for single receivers. In our study, we propose estimation of the receiver DCB and correction of the satellite DCB with existing multi-GNSS DCB products (Wang et al. 2016). Taking triple-frequency observations as an example, the correction equation can be deduced as (Guo et al. 2015)

$$\begin{cases} P\_{r,1}^{s} = \rho\_r^s + \vec{d}t\_r - dt\_{pre}^s + dT + a\_1 \cdot dl\_{r,1}^s - \frac{1}{a\_2 - 1} DCB\_{12}^s + \varepsilon\_{P\_f} \\ P\_{r,2}^{s} = \rho\_r^s + \vec{d}t\_r - dt\_{pre}^s + dT + a\_2 \cdot dl\_{r,1}^s - \frac{a\_2}{a\_2 - 1} DCB\_{12}^s + DCB\_{r,12} + \varepsilon\_{P\_f} \\ P\_{r,3}^{s} = \rho\_r^s + \vec{d}t\_r - dt\_{pre}^s + dT + a\_3 \cdot dl\_{r,1}^s - DCB\_{13}^s - \frac{1}{a\_2 - 1} DCB\_{12}^s + DCB\_{r,13} + \varepsilon\_{P\_f} \end{cases} (3)$$

where <sup>12</sup> = <sup>2</sup> − <sup>1</sup> , ,12 = ,2 − ,1 , and ̅ = + ,1 . The <sup>12</sup> and <sup>13</sup> can be obtained from multi-GNSS DCB products, while ,12 and ,13 are estimated as daily constant parameters. Similarly, the phase equations can be rewritten as

$$
\Phi\_{r,f}^{s} = \rho\_r^s + \bar{d}t\_r - dt\_{pre}^s + dT - a\_f \cdot dI\_{r,1}^s + \lambda\_f \overline{N}\_{r,f}^s + \varepsilon\_{\Phi\_f} \tag{4}
$$

where the ambiguity can be re-parameterized as

$$\begin{cases} \overline{N\_{r,f}^s} = N\_{r,f}^s + b\_{r,f} - b\_f^s \\ \quad b\_{r,f} = B\_{r,f} - \frac{D\_{r,1}}{\lambda\_f} \\ \quad b\_f^S = B\_f^S - \frac{D\_{IF}^S}{\lambda\_f} \end{cases} \tag{5}$$

and the estimable parameters are

$$X = \begin{bmatrix} x & y & z & \overline{d}t\_r & dT & I\_{r,1}^s & D\_{r,12} & D\_{r,13} & \overline{N}\_{r,1}^s & \overline{N}\_{r,2}^s & \overline{N}\_{r,3}^s \end{bmatrix} \tag{6}$$

Compared with the model in Li et al. (2018), the structure of the unknown parameters, except (, , ), is different, due to the different strategies of DCB correction. In Li et al. (2018), the dual-frequency DCBs are absorbed by other parameters, whereas the third

frequency DCBs are estimated. In this study, the satellite DCBs for all the three frequencies are corrected with existing DCB products (Wang et al. 2016) and the receiver DCBs are estimated. Consequently, our ionospheric parameters will not be biased by DCBs, which is beneficial for ionospheric modelling. In addition, for single stations with n observable satellites, the number of DCB parameters to be estimated in their model is n, while it is 2 in our model. The degree of freedom of our model is larger, which could increase the redundancy and robustness of the positioning solutions. The estimated ambiguity parameter is a combination of the integer ambiguity, the corresponding code hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite ends. In order to recover its integer property, these biases, i.e., satellite FCB and receiver FCB ,, must be accounted for. Normally, the receiver FCB is not of concern as it can be eliminated when performing single differences of observations between satellites. The satellite FCB, however, must be estimated at the server end and broadcasted to the users.

#### **2.2 FCB Estimation Strategy**

In dual-frequency IF PPP, the float ambiguity is usually decomposed into WL/NL forms in order to recover the integer property (Ge et al. 2008). This is partly because the IF combination of L1/L2 ambiguities is, in essence, not an integer. Another reason is that the WL ambiguities possess a relatively longer wavelength and are less correlated, therefore can be easily fixed. For uncombined PPP AR, it is also important to form combinations of raw ambiguities. On the one hand, the estimated raw float ambiguities are strongly correlated. On the other hand, the raw float ambiguities are quite sensitive to unmodeled ionospheric errors (Gu et al. 2015b). Therefore, the combinations with longer wavelengths and lower ionospheric delays are preferred. The coefficients must be integers in order to preserve the integer nature of ambiguities. In addition, these combinations should be independent to avoid rank-deficiency. While the Least-squares AMBiguity Decorrelation Adjustment (LAMBDA) method can be used to automatically search for the optimal linear combinations of ambiguities (Li et al. 2018; Teunissen et al. 1997), the classic extra-/wide-lane ambiguities (EWL/WL) were found to perform equally well and were used in our work to simplify the algorithm. For triple-frequency observations, it is easy to find two optimal combinations, e.g., one EWL and one WL combination or two WL combinations. The searching of the third combination, however, is much more difficult. From the systematic investigation of triple-frequency combinations, it is found that (4, −3, 0) is a good compromise between ionospheric reduction and noise amplification (Li et al. 2018). Concerning the properties of the above combinations shown in Table 1, they are denoted as NL/WL/EWL without specific explanation herein. In addition, they are used for both BDS and Galileo triple-frequency observations for simplicity

$$
\begin{bmatrix}
\overline{N}^s\_{r,LC1} \\
\overline{N}^s\_{r,LC2} \\
\overline{N}^s\_{r,LC3}
\end{bmatrix} = \begin{bmatrix}
4 & -3 & 0 \\
1 & -1 & 0 \\
1 & 0 & -1
\end{bmatrix} \begin{bmatrix}
\overline{N}^s\_{r,1} \\
\overline{N}^s\_{r,2} \\
\overline{N}^s\_{r,3}
\end{bmatrix} \tag{7}
$$

Substituting (5) into the above system produces the basic model for estimating FCBs. Since they have the same structure, a general expression can be formulated as

$$R\_{r,LC}^{s} = \overline{N}\_{r,LC}^{s} - \hat{N}\_{r,LC}^{s} = d\_{r,LC} - d\_{LC}^{s} \tag{8}$$

for all linear combinations, , denotes the float combined ambiguities; ̂ denotes the integer part of , ; , and denotes the receiver and satellite FCBs; , represents the FCB measurements. For each linear combination, a set of equations in the form of (8) can be generated, based on a network of reference stations. Suppose that there are n satellites tracked by m reference stations, the system of equations can be expressed as

$$
\begin{bmatrix} R\_{1,LC}^{1} \\ \vdots \\ R\_{1,LC}^{n} \\ \vdots \\ R\_{m,LC}^{n} \\ \vdots \\ R\_{m,LC}^{n} \end{bmatrix} = \begin{bmatrix} \overline{N}\_{1,LC}^{1} - \overline{N}\_{1,LC}^{1} \\ \vdots \\ \overline{N}\_{1,LC}^{n} - \overline{N}\_{1,LC}^{n} \\ \vdots \\ \overline{N}\_{m,LC}^{n} - \delta \underline{1}\_{m,LC}^{1} \\ \vdots \\ \overline{N}\_{m,LC}^{n} - \delta \underline{1}\_{m,LC}^{1} \end{bmatrix} = \begin{bmatrix} 1 & \cdots & 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cdots & 0 & 0 & \cdots & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & \cdots & -1 \end{bmatrix} \begin{bmatrix} d\_{1,LC} \\ \vdots \\ d\_{m,LC} \\ d\_{1C} \\ \vdots \\ d\_{LC}^{n} \end{bmatrix} \tag{9}
$$

where *LC* stand for the linear combinations (1, 2, 3, ⋯ ). The obtained system is singular on both sides of the equations. For the left side, ̂, can be determined by rounding ̅, , assuming that the float ambiguities are precisely estimated. For the right side, one arbitrarily combined FCB should be set to zero. For all the linear combinations, we always set the combined FCB of the last satellite to zero, i.e., G32/E30/C14. Note that the FCB measurements , from different stations may differ with ±1 cycle. This is due to the rounding process and can be adjusted with the strategy described in Xiao et al. (2018b). In this way, the system of equations can be solved. With the obtained combined FCB, we are able to calculate the FCB of the raw L1/L2/L3 carrier frequency

$$
\begin{bmatrix} d\_1^s \\ d\_2^s \\ d\_3^s \end{bmatrix} = \begin{bmatrix} 4 & -3 & 0 \\ 1 & -1 & 0 \\ 1 & 0 & -1 \end{bmatrix}^{-1} \begin{bmatrix} d\_{LC1}^s \\ d\_{LC2}^s \\ d\_{LC3}^s \end{bmatrix} \tag{10}
$$

The transformation from combined FCBs to raw FCBs is important, as it provides more flexibility to the users. With the raw FCBs, users are able to choose their own linear combinations of observations, formulate the corresponding combined FCB, and conduct PPP AR. This representation allows interoperability if the server and user sides implement different AR methods. In addition, the raw FCB is suitable for the State Space Representation (SSR) of Radio Technical Commission for Maritime services (RTCM) (Weber et al. 2005), where one phase bias per phase observable is broadcasted instead of making specific combinations.


**Table 1** Properties of GPS, Galileo, and BDS triple-frequency linear combinations.

#### **2.3 Uncombined PPP AR at the User End**

Similar to dual-frequency IF PPP AR, single differencing across satellites must be firstly performed in order to remove receiver FCBs. Then the single-differenced ambiguities from different carrier frequencies are combined, as has been done during FCB estimation

$$
\begin{bmatrix}
\overline{N}\_{r, \ell jk1}^{m,n} \\
\overline{N}\_{r, \ell jk2}^{m,n} \\
N\_{r, \ell jk3}^{m,n}
\end{bmatrix} = \begin{bmatrix}
i\_1 & j\_1 & k\_1 \\
i\_2 & j\_2 & k\_2 \\
i\_3 & j\_3 & k\_3
\end{bmatrix} \begin{bmatrix}
\overline{N}\_{r,1}^{m,n} \\
N\_{r,2}^{m,n} \\
\overline{N}\_{r,3}^{m,n}
\end{bmatrix} \tag{11}
$$

where ̅ , = − is the single-differenced ambiguity between satellites m and n. Based on the coefficients (,, ), the FCB for the specific combined ambiguity can also be formed. Note that the linear combinations are not necessary to be the same as that in FCB estimation, although the three combinations mentioned above are strongly recommended. In our experiments, we have used the same combinations as in FCB generation for uncombined PPP AR.

Usually, the EWL/WL float ambiguities can be directly fixed by the rounding approach after the correction of FCB (Ge et al. 2008), and the NL float ambiguities are fed into the LAMBDA algorithm to search for correct integers (Teunissen et al. 1997). However, in our study, the LAMBDA is used for each combination, regardless of its property, which simplifies the design of the algorithm. In addition, if not all the float ambiguities can be fixed by the LAMBDA method, partial ambiguity resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999). It is found that the searching and fixing of ambiguities for the combination with longer wavelengths (e.g., EWL/WL) is quite fast. When the integer ambiguities for one combination are resolved and validated, a tight constraint can be reconstructed. The number of constraints accumulate as the process repeats for all linear combinations. Afterwards, the constraints are imposed on the raw ambiguities, and yields the AR solution. Note that the ambiguities in IF PPP AR must be sequentially fixed in the order of WL/NL. An IF ambiguity is constrained only when both its WL and NL ambiguities are fixed, while the linear combined ambiguities in our study can be fixed and constrained independently.

## **3 Results and Discussion**

In this section, the data and processing strategy is described, followed by an analysis of the quality of triple-frequency FCB estimation. We further convert the combined FCB to raw FCB on each carrier frequency in order to characterize their properties. The generated FCBs are used to evaluate the performance of triple-frequency PPP AR solutions.

## **3.1 Data and Processing Strategy**

The International GNSS Service (IGS) established the MGEX in order to prepare operational services for new and upcoming GNSS (Montenbruck et al. 2017). The MGEX network comprises over 220 MGEX stations, as of October 2017. The daily observations from September 7–October 27, 2017—in total 51 days—were collected. About 200 stations were used for Galileo FCB estimation, of which 160 stations provide E1/E5a/E5b triplefrequency observations. About 150 stations were used for BDS FCB estimation, of which 60 stations provide B1/B2/B3 triple-frequency observations. Figure 1 shows the geographic distribution of the MGEX stations with Galileo and BDS triple-frequency observations. These data provide almost full and continuous tracking of Galileo and BDS signals.

**Figure 1** Geographic distribution of the selected MGEX stations with Galileo (represented by empty circles) and BDS (represented by solid triangles) triple-frequency observations.

In the processing, the E1/E5a/E5b were used for Galileo, while B1/B2/B3 were used for BDS. Data from GPS L1/L2 was also used to test the efficiency of the proposed approach in the case of dual-frequency observations, while the L3 was excluded due to inter-frequency clock bias (Pan et al. 2018). The cut-off elevation angle was set to 10°, while the float ambiguities with an elevation below 30° or with standard deviation (STD) larger than 0.1 m were removed for FCB estimation. It is noted that the BDS satellite-induced code multipath effects were corrected for inclined geosynchronous orbit (IGSO) and medium earth orbit (MEO) satellites, according to Wanninger and Beer (2015), while geostationary orbit (GEO) satellites were excluded from the processing. The third-generation BDS satellites, which were no longer affected by such effects (Lei et al. 2017), were also excluded due to no public data. Throughout the processing, MGEX precise products provided by Deutsches GeoForschungsZentrum (GFZ) (Uhlemann et al. 2016) were used. The satellite phase center offsets and variations were corrected according to the IGS antenna file. Since the antenna correction values for the third frequency, i.e., E5b and B3, were not available, we simply used that of the second frequency, i.e., E5a/B2. It is demonstrated that the satellite antenna characteristics of the third carrier frequency were quite similar to those of the second carrier frequency (Schönemann et al. 2011). However, it is a compromised strategy considering the precision of phase measurements. We have downweighed the observations of the third frequency by a factor of 4, compared with that of the first and the second carrier frequency. As for the receiver antenna phase center offsets and variations, the correction values for GPS were employed for both Galileo and BDS, in accordance with the principle of orbit and clock generation (Prange et al. 2017). For the combined GPS, Galileo, and BDS processing, the system related weighting ratio of GPS, Galileo, and BDS code observations was assumed to be 1:1:3, while the precision of the phase observations was assumed to be at the same level (Kazmierski et al. 2018). The detail of the used software and processing standards can be found in (Xiao et al. 2019b) and (Xiao et al. 2018b).

### **3.2 FCB Residual Distributions**

The performance of PPP AR depends on the quality of the FCBs, which can be indicated by the posterior residuals. In general, a highly consistent FCB estimation can be expected if the residuals are close to zero. Figures 2–4 present the distributions of the posterior residuals after FCB estimation for Galileo, BDS, and GPS, respectively. The subfigures refer to the different linear combinations.

**Figure 2** Distributions of posterior residuals of Galileo FCB for different linear combinations.

**Figure 3** Distributions of posterior residuals of BDS FCB for different linear combinations.

In general, all the histograms are symmetric and nearly centered at zero, following Gaussian distributions. These results indicate a good consistency between the input float ambiguities and the generated FCBs, which prove the efficiency of the proposed FCB estimation strategy. However, the characteristics of residuals differ with respect to the combinations and systems.

**Figure 4** Distributions of posterior residuals of GPS FCB for different linear combinations.

For all the systems, the residuals of linear combinations with larger wavelengths are smaller. For example, the residuals of the NL combination with wavelength of around 10 cm, are larger than those of the other two combinations, i.e., WL/EWL. The reason is that the combination with smaller wavelength is susceptible to errors. An exception is that the BDS WL residuals are larger than that of NL. The reason is not clear, and we suspect that it may be related to the satellite induced multipath effect (Wanninger and Beer 2015).

When comparing the residuals from multi-GNSS, it is found that the Galileo WL/EWL outper-formed those of GPS and BDS. As discussed in Xiao et al. (2019b), the signal of Galileo possesses a better performance of multipath suppression, which may explain the results. For the NL, the residuals of GPS are the smallest, which is reasonable as the accuracy of the GPS PPP float solution is the highest.

Furthermore, the results are also different from that of IF PPP, in which the performance of WL is worse than that of the NL. The residuals of WL/NL in the uncombined PPP model are almost comparable in terms of root mean square (RMS) and distributions. The possible reason is that the WL ambiguities are directly formed from raw ambiguities in uncombined PPP, while it is derived from MW combinations in IF PPP. The noise of MW combinations is larger as code measurements are employed. In addition, the sample rate of WL FCB is 15 min in uncombined PPP, while that of IF PPP is 24 h. The larger sample interval may also increase the residuals.

#### **3.3 FCB Time Series**

For real time applications, another question of interest is the temporal stability of the FCB estimates. It would be possible to predict the FCB if they are stable over time. Figures 5–7 present the time series of FCBs for Galileo, BDS, and GPS, respectively.

**Figure 5** Time series of the combined (upper) and raw (bottom) Galileo FCB in each 15 min session on DoY 255, 2017.

**Figure 6** Time series of the combined (upper) and raw (bottom) BDS FCB in each 15 min session on DoY 255, 2017.

**Figure 7** Time series of the combined (upper) and raw (bottom) GPS FCB in each 15 min session on DoY 255, 2017.

In general, all the time series of FCB are quite stable over time. The fluctuations between adjacent sessions are smaller than 0.05 cycles, which indicates that the 15 min interval is sufficient for FCB estimation. When comparing the results from different linear combinations, it is found that the time series of NL FCB are noisier than those of WL/EWL FCB. The NL FCBs, possessing a smaller wavelength of about 10 cm, are susceptible to errors. The EWL shows extremely small variations over time, with a standard deviation of 0.01 cycles, as presented in Figure 8. For all the three systems, the raw FCBs are much noisier than that of the combined FCB (Figure 9). The raw FCBs, also having smaller wavelengths around 20 cm, are susceptible to ionospheric residuals, while it is eliminated or decreased by linear combinations in combined FCB. The average STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 0.10 cycles. It is easier to predict the combined FCB, especially for the EWL/WL FCB. A lower update rate could be used to reduce the burden of communication. In this manner, it would be more efficient to broadcast combined FCB for real time applications.

**Figure 8** Mean STD of the combined FCB series for all the 51 days. Daily STD is calculated for each satellite FCB series. For each day, the mean STD of all satellite daily STDs is presented.

**Figure 9** Mean STD of the raw FCB series for all the 51 days. Daily STD is calculated for each satellite FCB series. For each day, the mean STD of all satellite daily STDs is presented.

When comparing the results from multi-GNSS, it can be seen that the STDs of Galileo EWL/WL FCBs are smaller than those of GPS and BDS, while that of the Galileo NL FCB is worse than GPS and BDS. The better quality of Galileo EWL/WL FCBs is likely attributed to the multipath suppression of Galileo signals, while the worse quality of Galileo NL FCB is due to the poor precision of satellite orbit and clock product. From Figure 9, it is found that the STDs of Galileo raw FCBs are smaller than that of GPS and BDS, regardless of combinations for most of the days. The smaller STDs facilitate the prediction of FCBs, which indicates a promising future for real time applications.

#### **3.4 Triple-Frequency PPP AR**

In order to validate our FCB estimates, as well as to assess the performance of triplefrequency PPP AR, 11 days from DoY 250 to 260 in 2017 of MGEX network stations were processed in static PPP AR mode. The 24 h observations were divided into eight threehour sessions. The positional biases of BDS-only and Galileo-only PPP float solutions and AR solutions are presented in Figure 10 and 11. The positional biases are calculated with respect to the 24-h static GPS, Galileo, and BDS combined PPP solutions. The statistics of all the sessions in the 11 days are provided in Table 2. Note that the number of sessions for BDS is smaller than that of Galileo due to the regional BDS IGSO/MEO constellation.

**Figure 10** Convergence performance of BDS triple-frequency PPP float and AR solutions based on 804 3-h sessions under 68% confidence level.

**Figure 11** Convergence performance of Galileo triple-frequency PPP float and AR solutions based on 5805 3-h sessions under 68% confidence level.

It can be seen that the convergence time is significantly shortened by ambiguity resolution, especially for the east component. For Galileo triple-frequency observations, it takes 64.5 min for float solutions to converge to three-dimensional 10 cm accuracy, while that for AR solutions is only 56.0 min, corresponding to an improvement of 13.2%. For BDS triple-frequency observations, the corresponding numbers are 121.5 min, 97.0 min, and 20.2%, respectively. In addition, it can be seen that with the current constellation, the performance of Galileo already outperforms that of BDS, both in terms of float PPP and PPP AR.

From Table 2, it can be seen that PPP ambiguity resolution was able to enhance the accuracy for all the three components. The improvements of (east, north, up) components for positioning with BDS are (49.2%, 38.3%, and 29.6%), while that for positioning with Galileo are (60.0%, 29.0%, and 21.1%). The performance of BDS is worse than that of Galileo for both float and AR solutions. The worse performance of the BDS solution is likely ascribed to the incomplete convergence, which is due to the IGSO/MEO constellation and the limited number of observable satellites.


**Table 2** Accuracy comparison of Galileo and BDS triple-frequency float and AR solutions (Unit: cm).

The additional signals are expected to further enhance the performance of PPP AR, as has been discussed in previous research (Geng and Bock 2013; Gu et al. 2015a; Li et al. 2018). Therefore, we also conduct an experiment to investigate the benefit of the third frequency observations in addition to dual-frequency ones. The strategy is that the ambiguities of the dual-frequency observations, i.e., B1/B2 and E1/E5a, are resolved to integers, while the ambiguities from the third frequency observations are kept as float values. This is accomplished by deleting the third column and row of the matrix in Equation (10). Then the results are compared to that of resolving the ambiguities of all the threecarrier frequencies. It is found that the improvements of positional error and convergence time are minimal. The positional improvements for Galileo are 1.8%, 2.3%, and 1.7%, while that of BDS are 8.0%, 5.0%, and 7.8%, for east, north, and up components, respectively. The possible reason could be: (1) the third frequency observations coming from the same satellites as the dual-frequency observations, will not improve the geometry of satellites, i.e., the DOP value. Compared with new observations from new satellites or other GNSS, its contribution to the model strength is insignificant; (2) the third carrier frequency E5b of Galileo is very close to the second one E5a, which implies the contribution of E5b is almost negligible when E1/E5a WL ambiguities are resolved. For BDS, the contribution of B3 is slightly larger than that of E5b, as its carrier frequency difference with respect to the B2 is larger; (3) the third frequency observations have been

down-weighted due to lack of antenna corrections, which may also degrade the contribution of the third frequency.

## **4 Conclusions**

A unified model for multi-frequency PPP AR based on raw uncombined observations is proposed, which simplifies the concept of phase biases for AR. No assumption is made on the method used to determine FCB on the server end, which implies that the generated FCB from multi-frequency observations could be flexibly used, such as for dual- or triple-frequency ones. It is demonstrated that the model is extendable to dual- and triplefrequency observations.

To verify its efficiency, we processed 51 days of Galileo and BDS triple-frequency observations collected from globally distributed MGEX stations. The estimated FCB shows a good consistency with the input float ambiguities. The RMS of Galileo FCB residuals is 0.05 cycles, while that of BDS is 0.08 cycles. It is also observed that the residuals are smaller for the combinations with larger wavelengths. The results indicate that there may exist ionospheric errors, and combinations are required to reduce its influence. The average STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 0.10 cycles. To reduce the communication with servers for real time applications, it would be more efficient to broadcast linear combined FCB. The performance of triplefrequency PPP AR is assessed with 11 days of data in three-hour sessions. Compared with the float solutions, the positional biases of AR solutions are significantly improved. The improvements of ENU components for positioning with BDS are 49.2%, 38.3%, and 29.6%, while those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These results demonstrate the efficiency of the proposed FCB estimation approach, and that the triple-frequency PPP AR can bring an obvious benefit to the float solution.

When comparing triple-frequency PPP AR with that of dual-frequency, it is found that the contribution of the third frequency observations is minimal. The insignificant improvement of the third frequency observable may be due to its limited contribution to satellite geometry and the narrow deployment with respect to the second carrier frequency. Nevertheless, adding the third frequency increases the reliability since it is observed that the number of successful sessions is increased. The proposed model is also applicable to GPS triple-frequency observations, provided that the inter-frequency clock biases are accounted for, which will be investigated in the future.

**Author Contributions:** G.X. and P.L. conceived and design the experiments. G.X. performed the experiments, analyzed the data, and wrote the paper. Y.G. and B.H. reviewed the paper.

**Funding:** This research and the APC was funded by National Natural Science Foundation of China grant number [41674016, 41274016, 41604024].

**Acknowledgments:** Thanks go to IGS-MGEX and GFZ for providing GNSS data and products. G.X. is supported by the China Scholarship Council. This work was supported in part by the National Natural Science Foundation of China (grants nos. 41674016, 41274016 and 41604024).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


# **Chapter 5 Improved time-differenced cycle slip detect and repair for GNSS undifferenced observations**

#### **Guorui Xiao1,2 · Michael Mayer2 · Bernhard Heck2 · Lifen Sui1 · Tian Zeng1 · Dongming Zhao1,2**

<sup>1</sup> Zhengzhou Institute of Surveying and Mapping, 450052 Zhengzhou, China <sup>2</sup> Geodetic Institute, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany

Submitted: 17 January 2017; Accepted: 25 October 2017; Published: 04 November 2017

GPS Solutions 22:6, 2018. DOI: 10.1007/s10291-017-0677-7

This is an author-created version of the article with permission from Springer. The final publication is available at link.springer.com

**Abstract** GNSS undifferenced signal processing, in which the individual signal of each frequency is treated as independent observable, has drawn increasing interest in GNSS community. However, undifferenced signal processing brings new challenges for cycle slip detection and repair. One important feature is the carrier frequency identification of cycle slips since observationws are processed separately. An analysis of real cycle slips in a BDS triple-frequency baseline dataset illustrates the deficiencies in the cycle slip detection process commonly implemented, in the case when cycle slips occur in just one specific carrier frequency. Hence, we propose an improved cycle slip detection and repair approach based on a time-differenced model. Two major advantages characterize this proposed approach. The first one is a significant reduction of false alarms due to carrier frequency identification of cycle slips. Having access to a reliable cycle slip detection method significantly reduces the number of ambiguity parameters to be estimated. The second advantage is the benefit of separating the OSS (Observation at the other frequency of Same Satellite without cycle slip) and OSE (Observation at the Same Epoch from other satellites without cycle slip) from the OCS (Observation with Cycle Slip). The simulation results indicate that separation of the OSS and OSE can significantly improve the model strength of cycle slip estimation, especially OSS. The proposed approach is validated by cycle slip estimation with a real data set. Smaller biases and larger ratio values jointly demonstrate that a much stronger model strength can be achieved. Finally, the cycle slip repair procedure is applied to triple-frequency PPP. The stable and fast convergence, as well as the reduction of standard deviations, proves the efficiency of the proposed approach.

**Keywords** Cycle slip detection and repair · BDS triple-frequency measurements · Time-differenced model · PPP

## **1 Introduction**

Carrier phase measurements are essential for high-precision GNSS positioning, such as real-time kinematic (RTK) and precise point positioning (PPP), since they are much more accurate than pseudoranges. However, carrier phase measurements often suffer from cycle slips, which are discontinuities of phase ambiguity by an integer number of cycles. Such unexpected slips should be detected and repaired, if possible, before carrier phases are used as high-precision measurements.

In the past decades, the ionospheric-free (IF) combination is routinely employed in GNSS precise positioning. However, with emerging BDS and Galileo, as well as the modernization of GPS and GLONASS, various types of multi-frequency observables are available. Hence, the choice of optimum combinations becomes practically difficult given the diversity of equipment (Schönemann et al. 2011). Furthermore, the progress in ionospheric environment studies offers a priori knowledge on the ionospheric delay, notably accurate ionosphere correction models. This information cannot be employed due to the elimination of the ionospheric delay in the IF combination. In addition, the IF combination also amplifies the measurement noise by a factor of approx. 3, which degrades the performance of precise applications. As a result, a uniform solution, in which the individual signal of each frequency is treated as independent observable, has drawn increasing interest in the GNSS community. The efficiency of undifferenced signal processing has already been confirmed in terms of convergence and precision with multi-system PPP (Chen et al. 2015), as well as PPP-RTK (Feng et al. 2013a; Odijk et al. 2016a). Moreover, this new approach has been tested effectively for ionospheric modelling (Tu et al. 2013), Differential Code Bias estimation (Shi et al. 2016), and LEO orbit determination (Zehentner and Mayer-Gürr 2015).

However, undifferenced signal processing brings new challenges for cycle slip detection. The most significant feature is carrier frequency identification of cycle slips. Since all carrier frequency observations are processed separately, it is essential to identify the carrier frequency of the cycle slips to avoid contaminating other observations. Figure 1 provides a picture of the problem at hand. In a conventional method, such as TurboEdit algorithm, all carrier frequencies are flagged as containing cycle slips even if there is only one carrier frequency suffering from a cycle slip.

**Figure 1** Ambiguities of B1, B2 and B3 carrier frequencies from PPP for rover (left) and reference station (right). The jumps indicate cycle slips, which clearly do not occur simultaneously on all carrier frequencies.

Some algorithms dedicated to cycle slip detection/correction have been developed and implemented, which can be classified into three categories.

The first category contains methods based on time series analysis of GNSS observations. Cycle slips can be characterized as discontinuities in a smooth signal, i.e. in a signal that can be reasonably modeled by a multiple polynomial regression. A typical example for dual-frequency approaches is the TurboEdit algorithm, which makes use of the HMW linear combination (Hatch 1983; Melbourne 1985a; Wübbena 1985b) together with ionospheric residual combinations. A few modifications to the TurboEdit algorithm have been studied to further strengthen the cycle slip detection/correction under high sampling rate data (Cai et al. 2013b) and high ionospheric activity (Liu 2011b). Bayesian theory is also applied for the detection of cycle slips and outliers ((de Lacy et al. 2008). To deal with the small cycle slips in BDS GEO observations, Ju et al. (2017) propose a robust polynomial fit algorithm to provide an adaptive detection threshold. However, these methods require several minutes of continuous phase data before and after a cycle slip in order to satisfy the criteria for phase connection (Huang et al. 2016). Thus, they are not suitable for real-time applications.

The second category includes methods based on optimal combinations of multi-frequency observations. The most popular approach is to form linear combinations of observations or parameters to mitigate the presence of some error sources, such as geometric or ionospheric errors. (Dai et al. 2009) employed two geometry-free combinations of triple-frequency carrier phase measurements to detect and determine cycle slips. (Wu et al. 2010b) also used triple-frequency carrier phase combinations to detect and repair cycle slips. (de Lacy et al. 2012) defined more than five types of linear combinations of triple-frequency GNSS observations to detect and correct cycle slips in real time. (Zhao et al. 2015a) presented a real-time cycle slip detection and repair method based on independent linear combinations of undifferenced triple-frequency GNSS observations. However, a common characteristic of the above-mentioned triple-frequency methods is the forming of optimal combinations to mitigate the presence of geometric and ionospheric errors, which makes it impossible to identify the carrier frequency affected by a cycle slip.

The third category consists of methods relying on geometry-based and time-differenced models. The dual-frequency cycle slip correction method based on a time-differenced model has been investigated by (Banville and Langley 2013) and (Xiaohong and Xingxing 2012). (Zhang and Li 2015) extended this method to deal with triple-frequency observations and analyzed the benefit of a third frequency. In their model, time difference observation equations are formed between two consecutive epochs and all satellite equations are processed in an integrated adjustment showing great potential for carrier frequency identification of cycle slips. However, the routine process is to flag all carrier frequencies as containing cycle slips when a discontinuity is detected for one satellite (Banville and Langley 2013). None of the above publications have exploited the carrier frequency identification of cycle slips and analyzed the benefit.

The review of existing cycle slip detection and repair methods reveals that most approaches lack proper handling of carrier frequency identification of cycle slips. Given the characteristics of GNSS positioning based on undifferenced measurements, a sophisticated cycle slip detection and repair approach should possess the following properties: (1) Practical for both real-time and post-processing applications. (2) Applicable to all GNSS and flexible with respect to frequency (dual frequency, triple frequency, even single-frequency). (3) Capable of frequency identification of cycle slip. (4) Robust against measurement noise. Taking these aspects into account, we extend the time-differenced model, with an emphasis on carrier frequency identification of cycle slips and the benefit for cycle slip estimation.

The following section describes the mathematical model used to detect and estimate cycle slips. Then the benefit of carrier frequency identification for cycle slip correction is analyzed based on a simulation study and numerical experiments. The cycle slips in a real baseline dataset are analyzed. In addition, the performance of different methods is assessed based on triple-frequency uncombined PPP. Finally, the methodology described in this study is summarized, and an outlook for future research is presented.

## **2 Methodology**

As pointed out in the introduction, most existing cycle slip detection approaches for undifferenced observations assume that cycle slips occur on all frequencies simultaneously. When not all frequencies are affected by cycle slips simultaneously, such approaches may lead to numerous false alarms. Although the results presented here are based on BDS, this phenomenon also exists in other GNSS (e.g., GPS, GLONASS). For this reason, a general improved algorithm is described in detail.

#### **2.1 Time-differenced model for cycle slip detection**

The functional model of pseudorange and carrier phase observation can be formulated as follows (Banville and Langley 2013; Zhang and Li 2015):

$$\begin{cases} P\_l^j = \rho + a\_l I\_1^j + dt\_r - dt^j + D\_l - D^j + dT + \varepsilon\_P \\ \Phi\_l^j = \rho - a\_l I\_1^j + dt\_r - dt^j + dT + \lambda\_l \{N\_l^j + b\_l - b^j\} + \varepsilon\_\Phi \end{cases} \tag{1}$$

where the subscript = (1,2,3) refers to a specific frequency, superscript refers to a specific satellite; denotes the range between receiver antenna and the phase center of satellite, including displacements due to earth tides, ocean loading, and relativistic effects; and are the clock offsets of satellite and receiver r respectively; is a constant= 1 <sup>2</sup> <sup>2</sup> ⁄ ; <sup>1</sup> is the slant ionospheric delay on the 1 frequency of satellite ; is the receiver hardware code delay; is satellite hardware code delay; is the slant tropospheric delay; and are the wavelength of the signal and the integer ambiguity in cycles; and are the receiver and satellite uncalibrated hardware phase delays at frequency and satellite ; , are the pseudorange and carrier phase measurement noise. Among the errors in (1), the effects of phase center offsets (PCO) and phase center variations (PCV) at satellite and receiver antenna, as well as phase windup, can be calculated using existing correction models. The corrections and can be modeled using external information such as precise clock products and tropospheric models. Differencing the observations with respect to time, which eliminates the effects of constant or slowly varying parameters, i.e. code instrumental biases, uncalibrated phase delays and wet component of tropospheric delay, a geometry-based and time-differenced model is formulated as:

$$\begin{cases} \Delta \check{\theta}\_{l}^{j} = \Delta \rho + \Delta dt\_{r} + a\_{l} \Delta l\_{1}^{j} + \Delta \varepsilon\_{P} \\ \Delta \check{\Phi}\_{l}^{j} = \Delta \rho + \Delta dt\_{r} - a\_{l} \Delta l\_{1}^{j} + \lambda\_{l} \Delta N\_{l}^{j} + \Delta \varepsilon\_{\Phi} \end{cases} \tag{2}$$

where ∆̌ and ∆̌ refer to the time-differenced pseudorange and phase observations corrected for the above mentioned errors. ∆ denotes the change in distance between two adjacent epochs. ∆ is the receiver clock offset variation and ∆<sup>1</sup> is the slant ionospheric delay variation on the <sup>1</sup> frequency of satellite .

The variation of the ambiguity parameter ∆ is zero for continuous observations, while it is an integer number when a cycle slip occurs. With the application of elevation-dependent weighting methods, a least-squares adjustment using observations from multiple satellites allows an estimation of the parameters above (Zhang and Li 2015). Assuming there are observable satellites with frequency signals tracked by a GNSS receiver, the system of equations involved can be defined as:

$$E\{L\} = A\mathbf{x}\_{non\_{dis}} + B\mathbf{x}\_{lon}, \\ D\{L\} = Q\_L \tag{3}$$

where \_ = [∆ ] in the case of a static receiver and \_ = [∆ ∆ ∆ ∆ ] in case of a moving receiver, = [∆<sup>1</sup> <sup>1</sup> ⋯ ∆<sup>1</sup> ] , and = [∆̌ 1 <sup>1</sup> ∆̌<sup>1</sup> <sup>1</sup> ⋯ ∆̌ <sup>1</sup> ∆̌ <sup>1</sup> ⋯ ∆̌ ∆̌ ] . The parameters in (3) have been divided into two groups representing non-dispersive \_ and dispersive effects . The design matrix and can easily be obtained by computing the partial derivatives of (2) with respect to the estimated parameters. is a diagonal elevation-dependent weight matrix assuming that temporal correlations between related epochs have been neglected. If external information regarding the ionosphere is available, a vector of pseudo-observations can be added to the system (Banville and Langley 2013).

Cycle slips can be detected by assessing the consistency of this solution and examining the residuals of the adjustment. Testing of the residuals can be accomplished by computing the normalized residuals and verifying if they exceed a predefined threshold (Baarda 1968):

$$V = \frac{|v\_{lj}|}{\sigma\_{v\_{lj}}} > \eta \tag{4}$$

where represents the estimated residual of the ℎ measurement of satellite ; is its precision and is the threshold value from the standard normal distribution. The observation with the largest normalized residual exceeding the threshold is then discarded from the adjustment and iterations are performed until the test of (4) is passed for all observations or there is no further redundancy. Readers are kindly referred to Teunissen (2010) for details.

In contrast, the most popular approaches to cycle slip detection are to form linear combinations of observations to mitigate the presence of some error sources, such as geometric errors (e.g. \_) or ionospheric errors (e.g. ). Then cycle slip detection can be achieved by examining the continuities of new combinations. The time-differenced model constitutes two clear advantages over combination-dependent methods. First, the time-differenced model allows identification of the exact carrier frequency at which the cycle slip occurs since it treats the signals of each carrier frequency as independent observables. Identifying cycle slips with combination-dependent approaches, however, entails that all frequencies are tagged even if only one frequency suffers from cycle slip event. This property is beneficial for precise GNSS positioning based on undifferenced measurements. For example, in the "detect-reset" model, only the ambiguity parameters from a specific carrier frequency are reset instead of resetting all ambiguity parameters. The reduction of ambiguity parameters can result in a faster and more stable solution. Second, using the time-differenced model allows separating the OSS (Observation at other frequency of Same Satellite without cycle slip) and OSE (Observation of Same Epoch from other satellites without cycle slip) from the OCS (Observation with Cycle Slip), thus the benefits of OSSs and OSEs can be fully exploited for the estimation of cycle slip.

#### **2.2 Improved model for cycle slip estimation**

Compared to the "detect-reset" model, a more ambitious approach of handling cycle slips is to estimate and correct the discontinuities in carrier phase, denoted as "detectrepair" model in this study. In order to take full advantage of the precision of carrier phase observations, once cycle slips are detected, the next step consists of estimating the size of the discontinuities. In the combination-dependent method, three independent combinations are employed to resolve the cycle slip parameters. This is typically done by a transformation matrix. It can be demonstrated, however, that those transformations will not improve the cycle slip estimation. Hence, using linear combinations, such as the Extra-Wide-Lane and Wide-Lane, will have no impact on the cycle slip estimation process (Banville and Langley 2013). In our approach, the float estimates are derived by the solution of the equations:

$$E\{L\} = A\mathbf{x}\_{non\_{dis}} + B\mathbf{x}\_{lon} + C\Delta N,\\ D\{L\} = Q\_L \tag{5}$$

where ∆ refers to the cycle slip parameters. Based on the analysis in the above sections, equation (5) can also be divided into three groups:

$$\begin{cases} OCS: L\_{OCS} = A\_1 \mathfrak{x}\_{non\_{dds}} + B\_1 \mathfrak{x}\_{lon}^1 + C \Delta N\\ OSS: L\_{OSS} = A\_2 \mathfrak{x}\_{non\_{dds}} + B\_2 \mathfrak{x}\_{lon}^1 \\ OSE: L\_{OSE} = A\_3 \mathfrak{x}\_{non\_{dds}} + B\_3 \mathfrak{x}\_{lon}^2 \end{cases} \tag{6}$$

For simplification, assuming there is a cycle slip detected at B1 frequency of the first satellite at a static receiver, the improved model for cycle slip estimation can be formulated as:

$$E\{L\} = \begin{bmatrix} OCS\begin{Bmatrix} 1\\1\\1\\\vdots\\\vdots\\1\\1\\\end{Bmatrix} & a\_m \begin{Bmatrix} 1\\\vdots\\\vdots\\a\_m \end{Bmatrix} & \vdots & \vdots\\a\_m \begin{Bmatrix} 1\\\vdots\\\end{Bmatrix} & a\_m \begin{Bmatrix} 0\\\vdots\\\end{Bmatrix} & \vdots\\a\_m \begin{Bmatrix} 0\\\vdots\\\end{Bmatrix} & a\_1\\\vdots\\\vdots\\\vdots\\\vdots\\\vdots\\\end{Bmatrix} & \vdots & \vdots\\a\_m \begin{Bmatrix} 1\\\vdots\\\vdots\\\vdots\\\end{Bmatrix} & a\_m \begin{Bmatrix} \Delta H\_r\\\vdots\\\Delta H\_1\\\vdots\\\Delta H\_r\\\vdots\\\vdots\\\vdots\\\end{Bmatrix} & \vdots\\a\_m \begin{Bmatrix} 0\\\vdots\\\end{Bmatrix} & a\_m \end{Bmatrix} & D\{L\} = Q\_L\tag{7}$$

It can be seen that the unbiased float estimates of ∆ can only be achieved through careful modeling and estimation of the two parameter groups (\_ and 1 ). For estimating \_, both OSS and OSE equations are beneficial since they have the parameter \_ in common with the OCS equations. When a receiver maintains lock on most satellites between the epochs, this parameter typically can be modeled with a precision of a few millimeters for static receivers or centimeters for moving receivers (Banville and Langley 2013). This advantage is not available for combination-dependent methods because the non-dispersive effect has been removed by forming linear combinations.

While \_ is shared by all satellite observations (OSSs, OSEs), the parameter 1 in the OCS equation is only shared in the observations of the same satellite (OSSs). This characteristic makes OSSs more critical than OSEs. If the OSS equations can be precisely solved, 1 can be achieved with phase-noise level precision, then the only parameter left in the OCS equations is the cycle slip parameter. With the precise modeling of \_ and 1 , it is easy to obtain cycle slip parameters with high precision. This procedure of obtaining the precise ionospheric delay variation constitutes a clear advantage over a priori constraints (Banville and Langley 2013; Zhang and Li 2015). Due to the capability of carrier frequency identification of cycle slips, the promising benefits of separating the OSS and OSE will motive high accuracy in cycle slip estimation. The benefits of OSSs and OSEs on the model strength are demonstrated in the following section.

Once the integer cycle slips are computed, they need to be validated before the repair procedure. The ambiguity dilution of precision (ADOP) based success rate provides a good approximation to the integer least-squares (ILS) success rate (Verhagen 2005). Therefore, the following approximation is applied (Zhang and Li 2015):

$$P\_{ILS} \approx P\_{ADOP} = \left[2\Phi\left(\frac{1}{2ADOP}\right) - 1\right]^n \tag{8}$$

where = √|̂ | 2n ,() = 1 √2 ∫ {− 1 2 2 } −∞ ,|∙| is the determinant operator ,̂ is the variance-covariance matrix, and is the dimension of ̂ . For practical applications, the ratio test is frequently applied to assess the closeness of the float solution to its nearest integer vector (Frei and Beutler 1990). Therefore, the ratio test is also applied to validate the integer cycle slip estimation.

## **3 Benefits of OSS and OSE**

In this section, we investigate the impacts of OSSs and OSEs on the model strength of cycle slip estimation by a simulation study. For the purpose of simulations, the following geometry-weighted model (Banville and Langley 2013; Zhang and Li 2015) is used. A few modifications, e.g. excluding the external constraints regarding the range and ionosphere, have been made to draw a clear distinction between internal and external information:

$$E\{\mathbf{y}\} = E\begin{Bmatrix}\Delta P\_1\\ \vdots\\ \Delta P\_m\\ \Delta \Phi\_1\\ \vdots\\ \Delta \Phi\_m \end{Bmatrix} = \begin{bmatrix} 1 & a\_1 & 0 & \cdots & 0\\ \vdots & \vdots & \vdots & \cdots & \vdots\\ 1 & a\_m & 0 & \cdots & 0\\ 1 & -a\_1 & 1 & \cdots & 0\\ \vdots & \vdots & 0 & \cdots & 0\\ 1 & -a\_m & 0 & \cdots & 1 \end{bmatrix} \begin{bmatrix} \Delta \rho\\ \Delta I\_1\\ \Delta N\_1\\ \vdots\\ \Delta N\_m \end{bmatrix}, D\{\mathbf{y}\} = Q\_{\mathbf{y}} \tag{9}$$

where = 2,3 indicates dual- and triple-frequency observations. is a diagonal matrix reflecting pseudorange noise and phase noise . The parameter ∆ which absorbs ∆ is estimated rather than the receiver displacement. Denoting the coefficient matrix in (9) by A, the covariance matrix ̂ is given as:

$$Q\_{\mathcal{X}} = \left[A^T Q\_{\mathcal{Y}}^{-1} A\right]^{-1} \tag{10}$$

The variance-covariance matrix of the cycle slip estimates ̂ can be derived from ̂ . In combination with (8), the ADOP-based success rate as function of a priori standard deviation in pseudorange and phase is investigated. Different cases are taken into consideration in case of potential poor quality of measurements.

#### **3.1 Benefits of OSS**

To investigate the impact of OSSs on the model strength, we assume that a cycle slip occurs only in the first carrier frequency for single satellite processing. This assumption implies that there are one and two OSS for the dual- and triple-frequency models, respectively. Figure 2 illustrates the ADOP-based success rate of dual- and triple-frequency models as function of a priori standard deviation in pseudorange and phase. Note that the phase standard deviation varies from 1 mm to 1 dm while the code standard deviation varies from 1 dm to 3.0 m. And "routine model" denotes the timedifferenced model presented by Banville and Langley 2013; Zhang and Li 2015, while the "improved model" represents our proposed approach.

**Figure 2** ADOP-based success rates of dual-frequency (left) and triple-frequency (right) models. The yellow surface represents the routine time-differenced model, while the blue and green ones represent the improved models for dual- and triple-frequency, respectively.

Figure 2 shows that the OSS has a significant impact on the ADOP-based success rate of dual and triple-frequency cycle slip estimation. The average success rate increases from 1.4 to 36.9% for dual-frequency and from 3.8 to 58.4% for triple-frequency models. When the standard deviations of code and phase observations are set to 3 dm and 1 mm respectively, the specific success rate increases from 8.1 to 95.8% for the dual-frequency model. The enormous increase shows that the OSS makes a significant contribution to improving the model strength of cycle slip estimation. It can be ascribed to the tight constraint for the ionospheric delay variation using precise OSS equations, while the constraint in the routine model is stemming from code measurements or previous epochs. In contrast to dual-frequency observations, the improvement for the triple-frequency model is only 7.3% (from 92.6 to 99.9%). This is reasonable since adding the third frequency already significantly improves the model strength. Thus the original success rate is relatively high and leaves narrow space for improvement.

In addition, it should be mentioned that the success rate decreases rapidly as code and phase noise increases for both models. Reliable cycle slip estimation (defined here as success rate > 99%) would be quickly impossible, which indicates that both code and phase noise are crucial for the cycle slip estimation. However, a better performance of the improved algorithm can still be observed. The gradient of decrease is much smaller, which indicates that the improved method is more efficient when pseudorange and phase noise are high. Therefore, it is possible to tolerate a larger uncertainty in the measurement noise, especially the code noise.

## **3.2 Benefits of OSE**

To investigate the impact of OSEs on the model strength and to avoid the impact of OSSs, we assume that cycle slips occur in all carrier frequencies of one satellite. If (n+1) satellites are observed in one epoch, then 2n and 3n OSEs are available for dual- and triplefrequency cycle slip estimation, respectively.

**Figure 3** ADOP-based success rates of dual-frequency (left) and triple-frequency (right) observations. The yellow surface represents satellite-by-satellite processing, while the blue and red ones represent the usage of 2 and 10 satellites, respectively.

From Figure 3, we can see that the OSE also has a clear impact on the ADOP-based success rate. The average success rate in the dual-frequency case increases by 15.4% when adding an extra OSE. However, the gradient of increment would decrease rapidly when more OSEs are applied, for example an increase of only 5.2% is observed for 8 additional OSE. The same phenomenon can be observed for triple-frequency observations. The possible reason is: OSEs and OCSs share only one common parameter ∆. Only by improving the estimation of ∆, can a OSE improve the cycle slip estimation. Thus, a remarkable improvement for ∆ estimation can be achieved when adding the first OSE. However, when a highly precise estimation of ∆ is already available, the contribution of adding more OSEs would not be as significant as the first one.

Note that the benefit of OSEs is not as significant as that of OSSs. Precise OSSs could not only improve the estimation of ∆ but also of ∆. Therefore, OSSs affect the cycle slip estimation more than OSEs.

#### **3.3 Joint benefits of OSS and OSE**

The joint benefits of OSSs and OSEs are presented under the assumption that only the first carrier frequency from the first satellite suffers from a cycle slip among all 10 observable satellites.

**Figure 4** ADOP-based success rates of dual-frequency (left) and triple-frequency (right) observations. The yellow surface represents the routine time-differenced model, while the green one represents the improved algorithm with joint application of OSSs and OSEs.


**Table 1** Average ADOP-based success rates of all noise conditions and the normal noise conditions

From Figure 4 and Table 1 we can see that the joint application of OSSs and OSEs has the strongest contribution to cycle slip estimation, especially for dual-frequency data. When routine methods are applied, it is practically impossible to reliably fix cycle slips for dual-frequency observations, while the success rate increases to 99.99% when the improved method is applied, even comparable with the triple-frequency. It is worth mentioning that the success rate decreases very slowly as the code noise increases, which indicates that the improved method is capable of handling noisy code measurements.

## **4 Results, comparison, and analysis**

A static baseline dataset was collected by two dual-system (GPS/BDS) penta-frequency (GPS-L1/L2, BDS-B1/B2/B3) receivers from Beijing Unicore Communications Incorporation on November 20, 2013. The observations were collected from 02:56:47 to 06:56:51 for the rover station and from 02:49:50 to 10:20:57 for the reference station with a distance of 5 m. It should be noted that we processed the data of the two stations separately with PPP although they were collected for RTK solution. The original observations with 1 s interval were resampled to 30 s in accordance with-common practice. Totally, 480 epochs were inspected from 02:57:00 to 06:56:30. Frequent cycle slips occur in this dataset. The BDS satellite elevations are shown in Figure 5. It can be seen that the elevation of the C04 satellite is relatively low (about 25 degrees) and a quick increment of the C12 satellite elevation from 15 to 52 degrees appears in two hours. The low elevation may introduce noisy measurements, which is an unfavorable situation for cycle slip detection and correction.

**Figure 5** Satellite elevations of the baseline dataset

#### **4.1 Cycle slip detection**

Comparing Figure 6 and Figure 7, one can easily note that a significant reduction of false alarms is achieved by the improved algorithm, particularly during the low-elevation part of C12 satellite for both stations. There are only 40 and 44 epochs detected with cycle slips by the improved algorithm for the rover and reference station respectively, compared with 127 and 138 epochs in the combination-dependent method. Most of the false alarms can be ascribed to the noisy measurements, for example due to the low satellite elevation of C12 at the start of the time series. In contrast, the improved algorithm can mitigate the impact of noisy measurements by applying elevation-dependent weighting when processing all satellite observations together.

**Figure 6** Results of cycle slip detection for rover (left) and reference stations (right) using geometry-free combinations (1, 1, -2) and (1, -2, 1) (Wu et al. 2010b; Zhen et al. 2008). Red triangle, green square, and blue inverted triangle represent B1, B2 and B3 carrier frequency observations, respectively.

**Figure 7** Results of cycle slip detection for rover (left) and reference (right) stations using the improved algorithm.

Furthermore, among the 40 and 44 epochs which suffer from cycle slips, only one epoch is affected by cycle slips in two satellites simultaneously. None of the other epochs suffers from cycle slips on more than one satellite, which clearly illustrates the benefit which can be achieved by OSEs. It is suboptimal to form linear combinations of observations and processing the data on a satellite-by-satellite basis. Among all epochs with detected cycle slips, only two epochs for the rover station and one epoch for the reference station suffer from cycle slips for all three carrier frequencies, which indicates that the benefit of OSSs could be achieved for most situations.

Besides, cycle slips are detected using both methods at five consecutive epochs from 3:36:30 to 3:38:30 of C02 for both rover and reference stations. After a detailed quality check of these observations, the discontinuities were attributed to system errors. This also provides a good opportunity to test the effectiveness of the validation procedure. Due to the capability of carrier frequency identification of cycle slips, the promising results of separating the OSS and OSE will motivate the high accuracy of cycle slip estimation, which is described in the next subsection.

## **4.2 Cycle slip estimation**

As stated in the second section, instead of flagging all frequencies as containing cycle slips, only cycle slip parameters on the affected frequency are estimated to fully exploit the advantage of OSSs and OSEs. Figure 8 presents the estimated magnitude of cycle slips for rover and reference stations, respectively.

**Figure 8** Magnitude of real cycle slips in rover (left) and reference (right) station data.

It can be seen that all cycle slip estimations are close to integer values. Besides, the improved algorithm shows great potential in dealing with various kinds of cycle slips: small cycle slips (e.g. one cycle), cycle slips occurring on two satellites simultaneously, and successive cycle slips at adjacent epochs.

In order to demonstrate the improvement, the routine time-difference method denoted as scheme 1 and the improved algorithm denoted as scheme 2 are implemented. Figure 9 and Figure 10 present the estimated biases with respect to the true amplitudes and ratio values for the above two schemes.

**Figure 9** Biases of cycle slip estimations for rover (left) and reference (left) station. Blue and red lines represent Scheme 1 (routine time-difference method) and Scheme 2 (improved algorithm), while the green line denotes the epochs with false alarms.

**Figure 10** Ratio values for rover (left) and reference (left) station.

Compared with the true amplitudes, the estimated biases of the improved algorithm are usually below 0.025 cycles, while they are about 0.05 cycles for the routine time-difference method. The closeness of the float solution makes it easier to fix the correct integers. This can be further illustrated by the ratio values (Figure 10). The mean values of the ratios derived by the improved algorithm are 202935 and 926307 for rover and reference stations respectively, while the values are 903 and 1530 derived from the routine timedifferenced method. Ratio values related to the improved algorithm are much larger, which indicates a more reliable cycle slip resolution.

The epochs with false alarms can easily be distinguished since all ratio values are below the threshold of 3. These false alarms in return proved the correctness of the developed cycle slip validation procedure. An exception can be seen at 5:02:30 at the reference station where the bias of Scheme 2 is larger compared to Scheme 1. As mentioned before, a quick elevation increment can be observed at C12 and this cycle slip occurred immediately at the beginning of the observation period of C12. Noisy measurements introduced by low satellite elevations may contaminate the estimation. Nevertheless, the magnitude of the bias is too small to affect achieving correct cycle slip integer values.

It is worth mentioning that the cycle slip estimations of the two methods are the same at epochs where all three frequencies suffer from cycle slips, for example at 5:35:30 and 6:36:00 for rover station, 5:35:30 for reference station. The reason is that there is no OSS if all carrier frequencies suffer from cycle slips, under this situation the improved algorithm is equivalent to the routine time-difference method, but still shows a clear advantage over the combination-dependent method.

These results show that much stronger model strength of cycle slip estimation can be achieved by separating OCS from OSS and OSE, which is consistent with the theoretical simulation.

## **4.3 Cycle slip repair for PPP**

The dataset is processed by a modified version of RTKLIB (Takasu 2010) to test the effectiveness of the improved algorithm. The modifications comprise PPP processing based on undifferenced triple-frequency observations, the proposed algorithm, the multi-frequency combination method (Wu et al. 2010b; Zhen et al. 2008), as well as the time-difference method (Banville and Langley 2013; Zhang and Li 2015).

Errors are corrected using the IGS standard error model (Kouba 2009), except for differential code bias (DCB) and PCO/PCV. Triple-frequency DCB biases are corrected by employing the CODE MEGX-DCB products (Guo et al. 2015). Second, empirical values are used for BDS satellite PCO corrections. Since the receiver PCO and PCV are not available, they are neglected. Considering that these terms are stable over short time interval, the proposed time-difference solution cancels their effects. The estimation of cycle slips could, therefore, be carried out without loss of performance. For BDS-PPP data processing, the precise products of satellite orbits and clocks are employed (Zhao et al. 2013). The observations are weighted according to the satellite elevation, with cutoff angle set to 15 degrees. The parameters including position, tropospheric zenith wet delay, receiver clock and ambiguities are estimated in a Kalman filter. The critical ratio values and ADOP-based success rates are chosen to be 3.0 and 0.99, respectively.

Three schemes are applied to the PPP software to investigate the positioning accuracy. The "MF detect" and "New detect" schemes represent cycle slip detection methods using multi-frequency combinations and the improved algorithm described in this study respectively, where no cycle slip repair procedure is attempted to emphasize the benefit of the proposed detection procedure. In the "Repair" scheme, the cycle slip estimation and correction procedure is applied.

As shown in Figure 11, the time series of NEU coordinate biases derived by the "New Detect" approach exhibit smaller variations and faster convergence than those obtained by "MF Detect". The advantage can also be seen from Figure 12, where the standard deviations of NEU coordinate biases derived by the "New Detect" approach decreased by 43%, 35%, 40% for the rover station, and 20%, 28%, 38% for the reference station, respectively. The better performance can be ascribed to the significant reduction in the number of ambiguity parameters achieved by the new algorithm. Based on the results of cycle slip detection, 252 ambiguity parameters need to be reset in the combinationdependent method while the number is reduced by approximately 64%, to only 91 for the "New Detect" approach. Having access to a reliable cycle slip detection method greatly reduces the number of ambiguity parameters to be estimated when processing GNSS undifferenced observations.

It should be noted that the improvement in convergence for the reference station is more remarkable than for the rover station. This can be ascribed to the larger occurrence rate of cycle slips for the reference station at the start time compared with rare cycle slips for the rover station (Figure 7). As the number of cycle slips increases, more false alarms appeared in the "MF Detect" scheme, while the "New Detect" approach can robustly identify the carrier frequency of the cycle slip.

**Figure 11** Time series of NEU coordinate biases of three schemes for rover (left) and reference (right) station

**Figure 12** NEU coordinate standard deviations of three schemes for rover (left) and reference (right) stations

The best performance of the three schemes can be achieved by cycle slip correction. After repairing cycle slips, the standard deviations decreased by 48%, 49%, 54% for rover station, and by 43%, 41%, 40% for reference station respectively, compared with the "MF detect" scheme. With the precise estimation of cycle slips and a rigorous ratio test, cycle slips are carefully repaired, and the stability of the NEU coordinate bias proved the correctness of cycle slip estimation in return.

## **5 Conclusions**

The analysis of real cycle slips in observed data clearly illustrates the deficiencies in the cycle slip detection process commonly implemented. In the dataset collected by two receivers, most cycle slips occur in just one specific carrier frequency. Under this circumstance the commonly applied combination-dependent methods of cycle slip detection, which are incapable of carrier frequency identification of cycle slips, flag all frequencies as containing cycle slips. Hence, we presented an improved approach based on a timedifferenced model for cycle slip detection and repair. For cycle slip detection, the proposed approach significantly reduces false alarms. Having access to a reliable cycle slip detection method greatly reduces the number of ambiguity parameters to be estimated for processing undifferenced GNSS measurements. Compared with the routine timedifferenced methods, the proposed approach can not only identify the carrier frequency in which the cycle slip occurs, but also makes it possible to separate the OSSs and OSEs, which greatly contribute to cycle slip estimations. The simulation results show that the theoretical success rates increase to 99.99% for both dual and triple frequency observations. Results from a real dataset also indicate that much stronger model strength of cycle slip estimation can be achieved.

Although the method is mainly applied to and validated with BDS triple-frequency observations, it is easily applicable to dual-frequency and even single-frequency observations for all GNSS, as long as phase redundancy is available. It is also applicable to GNSS positioning based on both combinations and undifferenced measurements. Further analysis on the occurrence and types of cycle slips in GNSS observations in the IGS network will be investigated in future.

**Acknowledgments:** Thanks to the IGS and Wuhan University for providing GNSS products. Thanks also go to Dr. Jinlong LI for providing data and Tomoji TAKASU for the open RTKLIB source code*.* We extend much gratitude to the editor in chief Prof. Alfred Leick and two anonymous reviewers for their valuable suggestions. Guorui XIAO is supported by the China Scholarship Council. This work was supported in part by the National Natural Science Foundation of China (grants nos. 41674016 and 41274016).

## **References**


## **Author Biographies**

**Guorui XIAO** is currently a Ph.D. candidate at Geodetic Institute, Karlsruhe Institute of Technology (KIT). He received his B.Sc. degree from the School of Geodesy and Geomatics in Wuhan University in 2011. His current research focuses on multi-frequency and multi-constellation GNSS precise positioning and applications.

**Michael Mayer** received his doctoral degree in 2005 from the Karlsruhe University (TH), where he was researching the modeling of the GPS deformation network Antarctic Peninsula. He is head of the GNSS working group of the Geodetic Institute and interested in mitigation of atmospheric and site-specific GNSS effects with a special focus on continuously operated reference sites.

**Bernhard Heck** received the Ph.D. degree and the Habilitation in Geodesy from the University of Karlsruhe (now: Karlsruhe Institute of Technology), Germany, in 1979 and 1984, respectively. Since 1991 he fills the position of Professor and Chair of Physical and Satellite Geodesy at KIT. His research interests include geometrical geodesy, deformation analysis, gravity field modelling, and GNSS positioning.

**Lifen SUI** is currently a professor at Zhengzhou Institute of Surveying and Mapping. She obtained her Ph.D. degree from Zhengzhou Institute of Surveying and Mapping in 2001. Her main research interests include surveying data processing and GNSS positioning.

**Tian Zeng** is currently a Ph.D. student of Zhengzhou Institute of Surveying and Mapping. He obtained his B.Sc. degree in 2014. His research interests specialize in GNSS data processing and precise orbit determination.

**Dongming Zhao** is an associate professor at Zhengzhou Institute of Surveying and Mapping. Since he received an Engineering Doctor degree in the field geodesy in 2004, he has more than 10 years of experience in academic research and teaching in physical geodesy and space geodesy.

## **Acknowledgments**

First of all, I would like to express my deepest gratitude to my supervisor, Prof. Bernhard Heck. His valuable guidance helped me in all the time of the research and writing of this thesis. His insightful comments and encouragements have been a continue support of my Ph.D. study. In addition, he has cared my life in Germany, and my family too. I could not have imagined having a better advisor and mentor for my Ph.D. study. I would like to thank Prof. Lambert Wanninger for serving as a second reviewer. I am quite grateful for his valuable discussions by email. Prof. Hansjörg Kutterer, who has chaired the Ph.D. defense and provided many supports, is thanked. Many thanks go to all further members of the doctoral committee for their interest in this thesis.

Thanks also go to the supervisors who support my visiting researches. They are Prof. Harald Schuh and Dr. Maorong Ge from The German Research Centre for Geosciences (GFZ), and Prof. Yang Gao from the Department of Geomatics Engineering of University of Calgary. They provided me an opportunity to join their team as intern, and gave access to the laboratory and research facilities. Thanks for the opportunities and inspiring discussions. A lot of friends at GFZ and University of Calgary had provided valuable help to me. Thanks for the warm welcome, fruitful cooperations and leisure activities we participated together.

A very special thank go to my collaborators, Prof. Lifen Sui and Dr. Pan Li, who had provided valuable suggestions on our publications.

I would also like to thank my colleagues at Geodetic institute and the Institute of Photogrammetry and Remote Sensing of Karlsruhe Institute of Technology, especially the PSG group, for their support.

Thanks also go to Tomoji TAKASU for the open RTKLIB source code. I also want to extend much gratitude to the editor in chief of GPS Solutions Prof. Alfred Leick and anonymous reviewers for their valuable suggestions.

Numerous institutions, organization, and individuals contribute to the IGS data, service and products, which are acknowledged. Special thanks to KIT GRAduate school for Climate and Environment (GRACE) for providing novel lectures and financial support. The China Scholarship Council is acknowledged for supporting my study in Germany under grant number 201503170213.

Living in a foreign land and pursuing a Ph.D. degree is not easy. Fortunately, a lot of friends were always there right beside me. They are CSC Ph.D. students at KIT, GuG students at GIK and IPF, colleagues in China, and basketball fans in Karlsruhe. Thanks for spending leisure time together, which is a continuous support of my study.

Finally, I owe my sincere gratitude to my family for their support without which I could not have focused on my study for these years. Words cannot express how grateful I am to them for all of the sacrifices that they have made on my behalf. I am looking forward to our common future.

This work and the publishing were funded by National Natural Science Foundation of China under grant number [41904039, 42274045].

**GPS precise point positioning (PPP) has found many scientifi c and industrial applications due to its cost-eff ectiveness, global coverage, and high accuracy for decades. However, it suff ers from a few drawbacks which limits further applications, e.g. slow convergence and unresolved ambiguity. The rapid development of multi-GNSS and multi-frequency signals off ers exciting prospects for further improvements. In this book, multi-frequency and multi-GNSS measurements from modernized satellites are properly integrated for PPP with ambiguity resolution to achieve the state-of-the-art fast and accurate positioning, which provides an important contribution to GNSS precise positioning and applications. The multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution, which is accomplished by a unifi ed model based on the uncombined PPP, are thoroughly evaluated with special focus on Galileo and BDS systems.** 

**Guorui Xiao** 

**Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution**