Document [original]

Front page

Online Estimation of Inter-Frequency/System

Phase Biases in Precise Positioning

vorgelegt von

M.Sc.

YUMIAO TIAN

geb. in Henan, China

von der Fakultät VI – Planen Bauen Umwelt

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktor der Ingenieurwissenschaften

– (Dr.-Ing.) –

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender: Prof. Dr. Jürgen Oberst, Technische Universität Berlin

Gutachter: Prof. Dr. Frank Neitzel, Technische Universität Berlin

Gutachter: Prof. Dr. Lambert Wanninger, Technische Universität Dresden

Gutachter: Dr. Maorong Ge, Deutsches GeoForschungsZentrum Potsdam

Tag der wissenschaftlichen Aussprache: 26. August 2016

Berlin 2016

Technische Universität Berlin

Fakultät VI- Institut für Geodäsie und Geoinformationstechnik

Eidesstattliche Erklärung

Hiermit versichere ich, dass ich die vorliegende Arbeit selbstständig verfasst und keine anderen als die

angegebenen Quellen und Hilfsmittel benutzt habe. Alle Ausführungen, die anderen veröffentlichten oder nicht

veröffentlichten Schriften wörtlich oder sinngemäß entnommen wurden, habe ich kenntlich gemacht.

Die Arbeit hat in gleicher oder ähnlicher Fassung noch keiner anderen Prüfungsbehörde vorgelegen.

Yumiao Tian Berlin, den 16.06.2016

Abstract

Global Navigation Satellite Systems (GNSS) play an important role in precise positioning for geodesy and

surveying engineering. The key to the real-time GNSS precise positioning is the instantaneous integer ambiguity

resolution. However, some of the biases in carrier phase observations cannot be removed by differencing

between either stations or satellites, so the integer nature of the double-differenced ambiguities is destroyed and

thus the ambiguities cannot be fixed to integers. Two typical biases are the inter-frequency bias (IFB) in GLObal

NAvigation Satellite System (GLONASS) data processing and the inter-system bias (ISB) in multi-GNSS

integration. Hence, the main objective of this thesis is the investigation, estimation and correction of these biases

in carrier phase observations to achieve better positioning accuracy, reliability and availability through the

improvement of its ambiguity resolution.

The estimated parameters of the carrier phase IFB and ISB are usually the IFB rate and the fractional ISB

(F-ISB), respectively. Most of the current methods estimate IFB rate or F-ISB together with the float ambiguities

and usually need observations of relatively long time due to their high correlation. Theoretically, the

performance of the ambiguity resolution depends on the quality of the given IFB rate/F-ISB value if the

observations are precisely modelled. In other words, the closer the given IFB rate/F-ISB value to the truth value

is, the better the resolution will be. Therefore, the RATIO in the ambiguity fixing can be applied as the

qualification factor of the IFB rate/F-ISB value. Based on this fact, a new methodology based on particle filter is

developed to estimate these biases in both post-processing and real-time mode in this study. In the proposed

method, the IFB/ISB is represented by its samples (i.e. particles) with the weights determined by the designed

likelihood function of the related RATIO given the sample values, so that the true bias value can be estimated

successfully by the particle filter approach. The integer nature of the ambiguities in the models with IFB/ISB

parameters is well utilised in the ambiguity resolution with the given IFB rate/F-ISB values. Thus, the new

method can significantly reduce the convergence time and increase the reliability of the estimation without a

priori values. Besides, when more than one bias parameter is included in the model, the multi-dimensional

particle filter approach is developed to estimate more than one bias parameter simultaneously in GNSS precise

positioning. In this case, the aforementioned benefits of the method are obviously enlarged.

In the GLONASS data processing with a nonzero IFB rate, the method can estimate the IFB rate from

observations of a few epochs. With the estimated IFB rate, the GLONASS fixed solutions are as accurate as the

GPS fixed solutions in the experiments with short baselines. In addition, the bias in the estimated IFB rate when

the state noise is set to a very small value or even zero is significant, but this bias can be removed by utilising the

regularized particle filter (RPF) and the precision of the estimated IFB rate is continuously improved by new

observations. An approach for adapting the number of particles in the estimation of the IFB rate is also proposed

to reduce the calculation burden by relating the number of particles to the standard deviation of the weighted

particles. In the estimation of the F-ISB in multi-GNSS integration, the new method based on particle filter

largely reduces the convergence time and improves the reliability of F-ISB estimation when satellites from each

system are not sufficient for independent positioning. Due to the periodic characteristics of ISB, the F-ISB

particles can be separated into different groups leading to the divergence of the filtering. This problem is solved

successfully by introducing the cluster analysis method which can detect the groups automatically so that they

can be shifted together into one group in the filtering.

The estimation of the phase IFB rate with the new method enables the usage of GLONASS in real-time

kinematic positioning even when the IFB between receivers is large. The estimation of the phase F-ISB with the

new method allows the precise positioning to be carried out with fewer satellites from each system than the

number of satellites required by the current methods. Therefore, the IFB rate/F-ISB estimation significantly

extends the application of real-time kinematic GNSS positioning. It also proves that the developed new method

is capable of estimating biases quickly and accurately, which initiates a new way of bias estimation in GNSS

precise positioning.

Keywords: GNSS carrier phase ∙ Integer ambiguity resolution ∙ RATIO ∙ Inter-frequency bias ∙ Inter-system bias

∙ Multi-GNSS integration ∙ Inter-system models ∙ Particle filter ∙ Regularized particle filter ∙ Adaptive number of

particles ∙ Multi-dimensional particle filter approach

Zusammenfassung

Global Navigation Satellite Systems (GNSS) spielen eine wichtige Rolle bei der präzisen Positionierung für

Geodäsie und Vermessungstechnik. Der Schlüssel für die präzise GNSS Echtzeitpositionierung ist die sofortige

Auflösung der ganzzahligen Mehrdeutigkeit. Jedoch können einige der Bias in Trägerphasenbeobachtungen

nicht durch Differenzbildung entweder zwischen Stationen oder Satelliten entfernt werden, so dass die

ganzzahlige Natur der Mehrdeutigkeit durch doppelte Differenzbildung zerstört werden kann und somit können

die Mehrdeutigkeit nicht als ganze Zahlen festgelegt werden. Zwei typische Biasarten sind die Inter-Frequency-

Bias (IFB) in der Prozessierung von GLObal NAvigation Satellite System (GLONASS) Daten und die Inter-

System-Bias (ISB) für Integration von mehreren GNSS. Daher ist das Hauptziel dieser Arbeit die Untersuchung,

Schätzung und Korrektur dieser Bias in Trägerphasenbeobachtungen um bessere Positionierungsgenauigkeit,

Zuverlässigkeit und Verfügbarkeit durch die Verbesserung der Auflösung von Mehrdeutigkeiten zu erreichen.

Die geschätzten Parameter der IFB und ISB von Trägerphasen sind normalerweise die IFB Rate und der

bruchzahlige Teil von ISB (F-ISB). Die meisten der aktuellen Methoden schätzen die IFB Rate oder F-ISB

zusammen mit den nicht-ganzzahligen Mehrdeutigkeiten und brauchen aufgrund ihrer hohen Korrelation

meistens relativ lange Beobachtungszeitintervalle. Theoretisch hängt die Leistung der Auflösung von

Mehrdeutigkeiten von der Qualität des gegebenen Wertes von IFB Rate/F-ISB ab, wenn die Beobachtungen

präzise modelliert werden. Mit anderen Worten, je näher der gegebene Wert von IFB Rate/F-ISB an dem wahren

Wert liegt, desto besser ist die Auflösung. Daher kann RATIO in der Festlegung von Mehrdeutigkeiten als

Qualifizierungsfaktor des Wertes von IFB-Rate/F-ISB angewendet werden. Aufgrund dieser Tatsache wurde in

dieser Arbeit eine neue auf dem Partikelfilter basierende Methode entwickelt, um diese Bias sowohl in Post-

Prozessierung als auch im Echtzeit-Modus zu schätzen. Bei dem vorgeschlagenen Verfahren wird die IFB/ISB

durch deren Stichproben (d.h. Partikel) repräsentiert, mit den Gewichten die durch die konstruierte

Wahrscheinlichkeitsverteilung (Likelihood-Funktion) von dem dazugehörigen RATIO für die gegebenen

Stichprobenwerte so festgelegt werden, dass der wahre Biaswert mit dem Partikelfilterverfahren erfolgreich

geschätzt werden kann. Die ganzzahlige Natur der Mehrdeutigkeiten in den Modellen mit IFB/ISB-Parametern

wird in der Auflösung von Mehrdeutigkeit mit dem gegebenen Wert von IFB Rate/F-ISB vorteilhaft verwendet.

Somit kann das neue Verfahren die Konvergenzzeit erheblich verringern und die Zuverlässigkeit der Schätzung

ohne a priori Werte erhöhen. Außerdem, für den Fall wenn mehr als ein Bias-Parameter in dem Modell enthalten

ist, wurde der mehrdimensionale Partikelfilter-Ansatz entwickelt, um mehr als einen Bias-Parameter gleichzeitig

innerhalb der präzisen GNSS-Positionierung abzuschätzen. In diesem Fall sind die oben genannten Vorteile des

Verfahrens noch offensichtlicher.

In der GLONASS-Datenverarbeitung mit einer Nicht-Null IFB Rate kann das Verfahren die IFB Rate aus den

Beobachtungen von einigen wenigen Epochen abschätzen. Mit der geschätzten IFB Rate sind in den

Experimenten mit kurzen Basislinien die GLONASS Lösungen mit festgesetzten Mehrdeutigkeiten so genau wie

die dazugehörigen GPS Lösungen. Zusätzlich ist das Bias in der geschätzten IFB Rate, wenn das

Zustandsrauschen als ein sehr kleiner Wert oder sogar Null festgelegt wird, signifikant, kann aber mit dem

regularisierten Partikelfilter (RPF) entfernt werden und die Präzision des geschätzten IFB wird mit neuen

Beobachtungen kontinuierlich verbessert. Ein Ansatz für die Anpassung der Anzahl der Teilchen in der

Schätzung der IFB Rate wurde auch vorgeschlagen, um die Berechnungslast zu reduzieren, indem die Anzahl

der Partikel mit der Standardabweichung der gewichteten Teilchen in Beziehung gesetzt wurde. In der Schätzung

der F-ISB in der Integration von mehreren GNSS reduziert das neue auf dem Partikelfilter basierende Verfahren

weitgehend die Konvergenzzeit und verbessert die Zuverlässigkeit der F-ISB Schätzung, wenn Satelliten von

jedem einzelnen System zu wenige für eine unabhängige Positionierung sind. Aufgrund der periodischen

Eigenschaften von ISB, können die F-ISB Partikel in verschiedene Gruppen getrennt werden, was zur Divergenz

der Filterung führen kann. Dieses Problem wird durch die Einführung der Cluster-Analyse, die die Gruppen

automatisch erkennen kann, so dass sie in der Filterung in eine Gruppe zusammengeführt werden können,

erfolgreich gelöst.

Die Schätzung der IFB Rate der Phase mit dem neuen Verfahren ermöglicht die Nutzung von GLONASS in

Echtzeit für kinematische Positionierung, auch wenn das Bias zwischen den Empfängern groß ist. Die Schätzung

des F-ISB der Phase mit dem neuen Verfahren erlaubt, dass die präzise Positionierung mit weniger Satelliten von

jedem System durchgeführt wird, als erforderlich für gängige Methoden. Daher erweitert die Schätzung der IFB

Rate/F-ISB bedeutend die Anwendung von kinematischen GNSS Echtzeitpositionierung. Es beweist auch, dass

die entwickelte neue Methode Bias schnell und genau schätzen kann, und eine neue Art der Biasschätzung in der

präzisen GNSS Positionierung einführt.

Stichwort: GNSS-Trägerphasen ∙ Ganzzahlige Mehrdeutigkeit Auflösung ∙ RATIO ∙ Inter-Frequency-Bias ∙

Inter-System-Bias ∙ Integration mehrerer GNSS ∙ Inter-System-Modelle ∙ Partikelfilter ∙ Regularisierte

Partikelfilter ∙ Adaptive Anzahl von Partikeln ∙ Mehrdimensionale Partikelfilter

Table of Contents

1 Introduction ................................................................................................................................................... 1

1.1 Research Background ............................................................................................................................... 1

1.1.1 Global Navigation Satellite Systems ............................................................................................. 1

1.1.2 Integer Ambiguity resolution ......................................................................................................... 2

1.1.3 GLONASS Inter-Frequency Biases ............................................................................................... 2

1.1.4 Multi-GNSS Inter-System Biases .................................................................................................. 3

1.1.5 Particle Filter Method .................................................................................................................... 4

1.2 Objective and Methodology ...................................................................................................................... 5

1.3 Main Contributions ................................................................................................................................... 5

1.4 Outline ...................................................................................................................................................... 6

2 Multi-GNSS Data Processing .................................................................................................................... 7

2.1 Satellite Constellations and Signals .......................................................................................................... 7

2.2 GNSS Observables and Error Sources .................................................................................................... 11

2.2.1 GNSS Observables ...................................................................................................................... 11

2.2.2 GNSS Error Sources .................................................................................................................... 12

2.3 Relative Positioning Models ................................................................................................................... 13

2.3.1 General Relative Positioning Models .......................................................................................... 13

2.3.2 Intra-System Models .................................................................................................................... 15

2.3.3 Inter-System Models .................................................................................................................... 16

2.3.4 Simplified General Relative Positioning Models for Short Baselines ......................................... 17

2.4 Parameter Estimation .............................................................................................................................. 17

3 Particle Filter ............................................................................................................................................... 20

3.1 Discrete-Time State-Space Model .......................................................................................................... 20

3.2 Bayesian filtering .................................................................................................................................... 20

3.3 Kalman Filter .......................................................................................................................................... 21

3.4 Sequential Importance Sampling ............................................................................................................ 22

3.5 Resampling ............................................................................................................................................. 24

3.6 Bootstrap Filter ....................................................................................................................................... 25

3.7 Regularized Particle Filter ...................................................................................................................... 26

3.8 Other Particle Filter Methods .................................................................................................................. 27

4 Online GLONASS Ambiguity Resolution Based on Online Phase IFB Estimation ................... 29

4.1 Between-Receiver Phase IFB Characteristics and Existing Methods for GLONASS IFB Estimation ... 29

4.2 Relationship between RATIO and Phase IFB Rate................................................................................. 30

4.3 Procedure for Phase IFB Rate Online Estimation ................................................................................... 33

4.4 Results and Analysis ............................................................................................................................... 34

4.4.1 Phase IFB Rate Estimation .......................................................................................................... 34

4.4.2 Computational Efficiency ............................................................................................................ 37

4.5 Regularized Approach............................................................................................................................. 37

4.5.1 Problem in State Noise Setting .................................................................................................... 38

4.5.2 Experiment with Regularized Approach ...................................................................................... 38

4.6 Adaptive Method for Setting the Number of Particles ............................................................................ 39

4.6.1 Proposed Adaptive Method ......................................................................................................... 41

4.6.2 Experiment with the Adaptive Method ........................................................................................ 41

4.7 Estimated Phase IFB Rates and their Characteristics .............................................................................. 42

4.8 Summary ................................................................................................................................................. 44

5 Online Inter-System Ambiguity Resolution Based on Online F-ISB Estimation ....................... 46

5.1 Existing Phase ISB Estimation Methods ................................................................................................. 46

5.2 Multi-GNSS Mathematic Models for the New Approach ....................................................................... 47

5.3 Relationship between RATIO and F-ISB ................................................................................................ 48

5.3.1 RATIO versus ISB of GPS L1 and Galileo E1 ............................................................................ 48

5.3.2 RATIO versus ISB of GPS L1 and GLONASS L1 ..................................................................... 51

5.3.3 RATIO versus ISB of GPS L1 and BDS B1 ................................................................................ 56

5.3.4 Half-Cycle Problem and Cluster Analysis Method ...................................................................... 59

5.3.5 Conclusions ................................................................................................................................. 59

5.4 Procedure for F-ISB Online Estimation .................................................................................................. 60

5.5 Results and Analysis ............................................................................................................................... 60

5.5.1 F-ISB Estimate Results ................................................................................................................ 60

5.5.2 Performance of the Solution for the Half-Cycle Problem ............................................................ 62

5.5.3 Computational Efficiency ............................................................................................................ 63

5.6 Analysis of F-ISB Characteristics in Multi-GNSS Integration ............................................................... 63

5.7 Summary ................................................................................................................................................. 65

6 Two-Dimensional Approach .................................................................................................................... 68

6.1 Motivation for Multi-dimensional Approach .......................................................................................... 68

6.2 Relationship between RATIO and Two F-ISB parameters ..................................................................... 68

6.3 Two-dimensional Particle Filter .............................................................................................................. 69

6.4 Experiments with Two Dimensional Approach ...................................................................................... 70

6.5 Summary ................................................................................................................................................. 72

7 Application of the Phase IFB Rate and F-ISB Estimation for Precise Positioning .................... 73

7.1 GLONASS Data Processing with Estimated Phase IFB Rate ................................................................. 73

7.1.1 Single-Epoch Processing ............................................................................................................. 73

7.1.2 Kinematic Positioning with Continuous Ambiguity .................................................................... 73

7.2 Multi-GNSS Data Processing with Estimated F-ISB .............................................................................. 75

7.2.1 Single-Epoch processing.............................................................................................................. 75

7.2.2 Kinematic Positioning with Continuous Ambiguity .................................................................... 77

7.3 Summary ................................................................................................................................................. 82

8 Conclusions and outlook ........................................................................................................................... 83

8.1 Contributions and Conclusions ............................................................................................................... 83

8.2 Outlook ................................................................................................................................................... 84

Bibliography ...................................................................................................................................................... 86

List of Tables ..................................................................................................................................................... 92

List of Figures ................................................................................................................................................... 93

List of Abbrevations ........................................................................................................................................ 96

Acknowledgments ............................................................................................................................................ 97

1 Introduction

The Global Navigation Satellite Systems (GNSS) precise positioning plays an important role in geodesy and

surveying engineering. The precise positioning can be realised by a single system, or by multi-GNSS integration

which can significantly enhance the positioning performance. The key to fast and precise positioning is the

carrier phase integer ambiguity resolution. However, the biases in carrier phase observations destroy the integer

nature of the ambiguities and hence lay obstacles on fast and precise positioning. Two typical biases are the

inter-frequency bias (IFB) in GLObal NAvigation Satellite System (GLONASS) data processing and the

inter-system bias (ISB) in multi-GNSS integration. In this thesis, the characteristics and estimation of these

biases will be investigated. The background and motivations of the research, the main contributions, as well as

the content of this thesis, are presented in this chapter.

1.1 Research Background

1.1.1 Global Navigation Satellite Systems

At present there are two developed GNSS, the United States (US) Global Positioning System (GPS) and the

Russian GLONASS, as well as two developing GNSS, the European Galileo system and the Chinese BeiDou

navigation satellite system (BDS). Besides, two other regional navigation satellite systems are also under

development, the Japanese Quasi-Zenith Satellite System (QZSS) and the Indian Regional Navigation Satellite

System (IRNSS). The systems GPS, GLONASS and BDS are providing fully or partially positioning, navigation

and timing (PNT) services, while the others are still in their test phase.

On one hand, these systems have the compatibility, which indicates that the PNT services of these systems can

be used separately or together without interfering with each other. On the other hand, GNSS have also

interoperability which means that the PNT services of all systems can be employed jointly to provide better

capabilities at the user level than the services of each single system as described in (Hein 2006, Li et al. 2015).

Therefore, the systems can be seen as one global system and the observation data can be processed unitedly,

which is referred to as multi-GNSS integration in this study. The signals of the same frequency from different

systems are even considered to be interchangeable and hence can be used together to estimate the parameters for

the services just as signals from one system, such as the signals of GPS L1 and Galileo E1 (Melgard et al. 2013,

Hein 2006). However, the multi-GNSS data processing encounters different reference frames which can be

unified in the orbit determination. Also, the carrier phase integer ambiguity resolution of the integration

potentially faces different biases in devices.

The available satellites from each system at one location at one moment are limited. For each complete single

system of GNSS, there are around ten satellites but not all of them are available even under good observation

conditions due to the station environment. The integration of multi-GNSS is able to provide a much larger

number of visible satellites, bringing benefits in three aspects (Ge et al. 2012, Force and Miller 2013, Odolinski

et al. 2014, Li et al. 2015). Firstly, multi-GNSS integration improves the availability of GNSS by increasing the

chances of observing enough satellites to estimate the parameters for the PNT services, which is important in the

scenarios such as in cities and mountain areas, where tall buildings and high mountains block the satellite signals

easily. Secondly, the integration improves the reliability which means that the results are more confident due to

additional observations from other providers. Thirdly, the integration improves the accuracy of the solutions,

especially in kinematic cases. Representatively, the performance of the integer ambiguity resolution can be

improved (Pratt et al. 1998, Odolinski et al. 2014, Ge et al. 2012).

Therefore, the multi-GNSS integration has become a hot topic, where GPS and GLONASS integration is firstly

focused (Wang et al 2001), and then Galileo, BDS and QZSS are also included (Odijk and Teunissen 2013a). A

main event in GNSS community was the launch of Multi-GNSS Experiment (MGEX) project by the

International GNSS Service (IGS) in 2012 to provide the multi-GNSS infrastructure for scientific and

engineering community (Rizos et al. 2013, Steigenberger et al. 2015). At present, the researches of multi-GNSS

integration cover all GNSS applications, such as GNSS positioning (Odolinski et al. 2014), time transfer (Dach

et al. 2006), atmosphere monitoring (Mayer et al. 2011) and so on.

1.1.2 Integer Ambiguity resolution

In precise GNSS positioning, usually both code pseudorange and carrier phase observations are employed. The

accuracy of the code pseudorange measurements is of decimetre to meter level, whereas the carrier phase

measurements have the accuracy of millimetre level. However, the carrier phase observations encounter

unknown cycle ambiguities which have to be estimated correctly.

GNSS receivers can only measure the fractional part of the carrier phase and record the accumulated cycle

numbers with sub-centimetre precision (Blewitt 1989), but the total number of cycles between the satellite and

the receiver is unknown, leading to the problem of ambiguity. Assuming other error sources, such as atmospheric

delays, have been accurately corrected, this ambiguity is still a float value in the non-difference (ND) model due

to the existence of uncalibrated phase delays (UPD), but is an integer after the UPD correction or by forming

double-difference (DD) model (Ge et al. 2008). In this thesis UPD are referred to as hardware delays which also

include the delays in the digital signal processing (DSP) in chips of devices and the initial phases. If the

hardware delays are properly handled, the ambiguity is an integer number which can be resolved. Besides, if the

GNSS signals are lost temporarily due to factors such as trees and buildings sheltering, the integer number of

cycles will encounter discontinuity known as the cycle slip which needs to be corrected or the integer ambiguity

has to be resolved again.

In the carrier phase observation model, the ambiguities are usually parameterised and estimated together with

other unknown variables, such as coordinate parameters. Firstly they are estimated as float values with the

ordinary least square (OLS) method and then the integer ambiguities are resolved. As the ambiguities are integer

numbers after hardware delays are removed or eliminated, constraining the ambiguities to integers can

significantly improve the accuracy of positioning results (Ge et al. 2008, Blewitt 1989, Dong and Bock 1989)

and shortens the convergence time in kinematic positioning (Li et al. 2013).

The integer ambiguity resolution methods are to get the integer ambiguity values based on the float solutions and

their corresponding variance-covariance (VC) matrix. The existing ambiguity resolution methods can be divided

into three classes (Kim and Langley 2000). The first class includes the methods in the measurement domain,

such as the ambiguity determination from the code pseudorange measurement or the code pseudorange

measurement with smoothing technique (Cocard and Geiger 1992). The second class includes the resolution

techniques in the coordinate domain, such as the ambiguity function method (AFM) (Counselman 1981, Han

1996). The third class includes the methods in the ambiguity domain, such as the Least-squares AMBiguity

Decorrelation Adjustment (LAMBDA) method (Teunissen 1995, Chang 2005). The LAMBDA is based on the

theory of integer least square (ILS), which provides optimal solution maximising the success rate of ambiguities

(Teunissen 1999). Due to the high efficiency and reliability, LAMBDA method has been widely used.

The resolved integer ambiguities are not always true in deterministic sense and wrong integer ambiguities can

seriously bias the fixed solution (Verhagen et al. 2012). Therefore, ambiguity validation methods have been

proposed to test whether the determined integer ambiguities should be accepted or refused. These methods

include the F-ratio test (Frei and Beutler 1990), R-ratio test (Euler and Schaffrin 1990), W-ratio test (Wang

1998), and difference test (Tiberius and Jonge 1995).

Among the ambiguity validation methods, the R-ratio test has been used widely. The R-ratio value test performs

very well, even though there is no theoretical criterion to select the threshold value because R-ratio value does

not strictly obey to any known regular distribution (Teunissen 1996, Teunissen 2003, Verhagen 2005, Verhagen

and Teunissen 2013). The traditional way of selecting the threshold is to set it to a fixed value by experiences.

Because the unknown parameters are estimated by comparing the RATIO values corresponding to different bias

samples in this study, the threshold value is not important. Thus, the R-ratio is employed and the threshold value

is set to a fixed value if it is needed.

1.1.3 GLONASS Inter-Frequency Biases

Among the existing satellite navigation systems, GLONASS employs the frequency division multiple access

(FDMA) technique which identifies satellites by different frequencies, while other GNSS employ code division

multiple access (CDMA) technique which identifies satellites by different pseudo-random noise (PRN) codes.

Consequently, for each satellite GLONASS uses different frequencies within the same frequency band. The

signals for each frequency pass through different paths inside the devices and therefore are biased with different

hardware delays and the bias is referred to as IFB.

GLONASS IFB exists on both code pseudorange and carrier phase observations. The IFB on code pseudorange

is not as well-regulated as that on carrier phase and the code pseudorange observations can be down-weighted in

precise positioning, so only the IFB on carrier phase is focused in this study. IFB values lump with the

ambiguities and are not integer multiple of the wavelength, resulting in the loss of the integer nature of the

ambiguities. If the IFB value is large, then the ambiguity resolution will fail and thus the accurate fixed solution

cannot be estimated.

To fix the integer DD-ambiguities in GLONASS relative positioning with non-zero IFB, a number of studies

have been carried out to investigate the IFB characteristics. Although it is confirmed that receivers of the same

manufacturer have in principle similar bias, there are anyway outliers (Wanninger 2012). Besides, in practice

employing devices from the same manufacturer cannot be always promised, especially when the number of

receiver manufacturers is increasing. Also, the antenna and cable, as well as the restart of receivers can also

contribute to the IFB (Wanninger and Wallstab-Freitag 2007). Therefore, we should not assume that IFB values

are the same and can always be eliminated in the differenced observations/ambiguities. It is also confirmed that

the GLONASS receiver IFB is nearly linearly correlated to the frequency number, and can be, therefore,

represented by a constant offset and the IFB rate with respect to the frequency number (Povaliaev 1997, Pratt et

al. 1998, Wanninger and Wallstab-Freitag 2007). Furthermore, it is also presented that the IFBs in L1 and L2 are

similar in distance as well. These characteristics can be utilized in the IFB modelling and estimating (Wanninger

2012, Al-Shaery et al. 2013).

Based on the linear relationship between IFB and the signal frequency number, several approaches have been

developed to estimate the IFB rate. Wanninger and Wallstab-Freitag (2007) and Wanninger (2012) employed

GPS and GLOANSS SD-observations between two stations to determine the GLONASS IFBs. However, this

method needs an a priori value of the IFB rate with certain accuracy, so that at least one of the ambiguities can

be fixed. Afterwards, the remaining ambiguities are estimated along with the IFB rate parameter. It was not

clearly addressed how to obtain such an initial IFB rate value. Zhang et al. (2011) ignored the differences in

wavelengths and estimated DD-ambiguities using one-day GPS and GLONASS data. The ambiguities are fixed

by simply rounding the float estimate to the nearest integer and then the IFB rate is calculated. Al-Shaery et al.

(2013) presented a method which estimates both the code pseudorange IFB rate and the carrier phase IFB rate

along with the float ambiguity solution. After the integer DD-ambiguities are fixed as integers, the IFB rates are

refined. The estimation method was applied to a zero-baseline with 23-hour GPS and GLONASS data in their

research.

Sleewagen et al. (2012) demonstrated that the main part of the linear correlation between IFB and the channel

number is caused by the code-phase bias in the DSP of the receivers. The code and phase measurements are

usually considered to share the same clock, but the time difference between the two measurements exists due to

the adjustment for code measurements in DSP, and the different paths from generator to correlator for code and

carrier phase. This code-phase bias induced by DSP is considered to compose the main part of the IFB. Based on

this conclusion, Banville et al. (2013) also proposed an approach to fix GLONASS ambiguity without any

external IFB calibration. However, this approach requires that two GLONASS satellites with adjacent frequency

numbers are observed simultaneously. Furthermore, in the demonstration the fixing rate is only 70% in the case

of using GLNOASS alone.

In general, almost all current approaches try to estimate the ambiguities and the IFB rate simultaneously.

However, due to the high correlation between the two sets of parameters, the estimation needs a long data set and

including the simultaneous GPS observations. Consequently, none of these methods can provide a fast or

real-time solution of IFB rate for GLONASS integer ambiguity resolution without an accurate a priori IFB value.

Besides the IFB, the wavelength difference of GLONASS satellites also affects the accuracy of DD-ambiguities

in the float solution, as the SD-ambiguities cannot merge together directly after being scaled into distances and

therefore initial SD-ambiguity values are needed. In this case, the bias in initial SD-ambiguity values affects the

performance of the ambiguity resolution. This problem can be solved by employing more accurate initial SD-

ambiguity values and has been thoroughly investigated and analysed (Leick 1998, Wang 2000). The

investigation in this thesis is based on short baseline, the initial SD-ambiguity values calculated directly from

code pseudorange observations can have satisfactory performance.

1.1.4 Multi-GNSS Inter-System Biases

The differences between GNSS have to be properly handled in multi-GNSS integration. These differences

include different geodetic references and timing systems, as well as different hardware delays in satellites and

multi-GNSS receivers. The differences of geodetic references can be removed by known transformation

parameters or by employing the orbit products which are free of these differences as the geodetic references have

been adjusted in the integrated satellite orbit determination, while the hardware delays vary with specific device

and are usually unknown. The hardware delays lump with the float ambiguities and do not affect the float

solution but have to be considered in integer ambiguity resolution.

In relative positioning, the multi-GNSS integration includes mainly two strategies in the ambiguity resolution.

The first strategy is to fix only intra-system DD-ambiguities of each system. The second strategy is to fix both

intra- and inter-system DD-ambiguities of the systems. The later strategy has more DD-ambiguities with integer

nature and is supposed to have better performance. However, the between-receiver ISB has to be removed so that

the DD-ambiguities in inter-system models can be fixed as integers (Odijk and Teunissen 2013a, Odijk and

Teunissen 2013b, Odolinski et al. 2014). In this thesis, the integrations between GPS L1 and Galileo E1 with the

same frequency, GPS L1 and GLONASS L1 with different frequencies, as well as GPS L1 and BDS B1 with

different frequencies but both employing CDMA technique, are taken as typical examples.

GPS L1 and Galileo E1 have the same frequency, thus the wavelengths of Galileo E1 and GPS L1 are the same.

The SD-ambiguity parameters can be differenced directly to form the DD-ambiguities even after being scaled

into distances. However, the ISB in the inter-system DD-model is necessary to be considered. This is the same as

the other GNSS integrations with systems of the same frequencies, such as the integration of GPS L5 and Galileo

E5a, Galileo E5b and BDS B2 (Odijk and Teunissen 2013b, Odolinski et al. 2014). Until now, the ISB value in

the integrations with the same frequency is considered to be the same for each pair receiver types (Odijk and

Teunissen 2013b).

GPS and GLONASS are still the only two GNSS systems with complete satellite constellations at present. The

integration with the strategy fixing only intra-system DD-ambiguities has been investigated in (Dai et al. 2001,

Han et al. 1999, He et al. 2016). The strategy fixing additional inter-system DD-ambiguities leads to two

problems. Besides the problem of ISB, the other problem is that the wavelength difference between the two

systems leads to the fact that the a priori values of SD-ambiguities which are usually calculated with code

pseudorange with relatively large errors, are included in the carrier phase DD-model (Wang et al. 2001, Meindl

2011). This problem is actually similar to the case of GLONASS only data processing (Leick 1998). The GPS

L1 and BDS B1 have also different frequencies and hence their integration is in a similar situation.

In the multi-GNSS integration with systems of the same frequency, the part of ISB, which is integer multiples of

the wavelength, lumps with the integer DD-ambiguities and does not affect the ambiguity resolution, but the

remaining fractional part of ISB (F-ISB), which is smaller than one wavelength, destroys the integer nature and

hence affects the ambiguity resolution. Therefore, some F-ISB estimation methods have been proposed. Odijk

and Teunissen (2013a) added the ISB parameter into the inter-system DD-models to preserve the integer nature

of the inter-system DD-ambiguities. One of the inter-system DD-ambiguity parameters, which include the

DD-ambiguity and ISB, is used to correct the ISB of other inter-system models so that the rank-deficiency

caused by the ISB parameter can be removed. This method is essentially the same as the method that fixes only

the intra-system DD-ambiguities of both systems and then refines the DD-ambiguity in the inter-system model.

In this method, the integer ambiguity resolution cannot benefit from the inter-system model before the F-ISB is

known (Odijk and Teunissen 2013a). Paziewski and Wielgosz (2015) tried to separate the DD-ambiguities and

ISB in the inter-system model by introducing a constraint condition which sets the F-ISB parameter to zero value

with a STD value equals to a half cycles. It is obvious that if the actual F-ISB is not zero but a half cycle, such

constraint condition is not helpful in ambiguity fixing in the estimation of F-ISB.

1.1.5 Particle Filter Method

Particle filter, also known as the sequential Monte Carlo method, is a Bayesian filtering method which is

implemented using Monte Carlo method. Monte Carlo methods are to approximate the solution of a problem

with computational mathematics, by a random process which determines the evolution of a sequence of states by

random events. As the number of random events approaches to infinitely large, the error of the approximation

can be infinitely small (Dimov and McKee 2008).

Even though the Bayesian filtering is optimal in principle, it is almost impossible to obtain its general optimal

analytical expressions of the posterior probability density function (PDF) except for some special cases. For

example, with the linear Gauss-Markov assumption, the optimal analytical expression exists as Kalman filter. In

most cases, only suboptimal models are available. One of the ways to get the suboptimal models is via the Monte

Carlo method, which expresses the PDF by a number of samples in a simulation way. If the number of the

samples is large enough, these samples with their associated weights can approximate the true PDF with errors

smaller than a given threshold. The particle filter is mainly composed of three steps, the prediction, resampling

and the measurement update (Gustafsson et al. 2002).

Since the first practical particle filter was proposed (Gordon et al. 1993), several algorithms of particle filter have

been developed for solving different problems. The auxiliary particle filter is developed to improve the

distribution of the samples according to the temporary measurements before update step (Pitt and Shephard

1999). The regularized particle filter (RPF) is also proposed to resolve the problem of diversity loss which is

usually caused by resampling step (Doucet et al. 2001). If the prediction model is linear, the transition of the

state vector can be completed with approach which is the same as in Kalman filter to lower the complexity. This

kind of particle filter is called Gaussian particle filter (Kotecha and Djuric 2003). Besides, adaptive particle filter

to resist large noise in measurements and distribution particle filter to decentralize the calculation have also been

proposed (Doucet et al 2001, Haug 2012).

Particle filter is able to solve the non-Gaussian and non-linear state space problems, which enables it to be

widely used within a short time. Its applications include target tracking, digital data processing, terrestrial

navigation, indoor navigation and others (Doucet et al. 2001). Even though particle filter has been utilised in

positioning, it is not widely used for geodetic GNSS precise positioning yet. The problems of IFB and F-ISB

estimation have attracted the attention of many researchers, which will be solved by employing particle filter in

this investigation.

1.2 Objective and Methodology

The main objective of this thesis is the estimation and applications of the phase IFB in GLONASS data

processing and the phase F-ISB in multi-GNSS integration for precise positioning, which can recover the integer

nature of the ambiguities so that the ambiguity fixing can succeed. The estimation procedures should be able to

estimate the carrier phase IFB rate and the F-ISB with short convergence time without an a priori value and to

track the biases online. Consequently, the availability, reliability and accuracy of GNSS can be improved

especially in severe environments.

It is obvious that a more accurate bias value can better remove the bias in carrier phase observations and recover

the integer nature of the ambiguities, leading to a relatively larger RATIO value. Therefore, if a number of bias

values are given, RATIO values can be employed to judge their qualities and hence to estimate the bias. This

methodology can be implemented by particle filter via designed likelihood function of RATIO in the update step.

This method considering the integer nature of the ambiguities in the estimation and hence the biases can be

estimated precisely with short convergence time.

Both precise relative positioning and precise point positioning (PPP) are based on the carrier phase observations

and encounter the biases mentioned above when integer ambiguity resolution is demanded. Although the

principle is also applicable to data processing in PPP, we focus on the problem in precise relative positioning.

1.3 Main Contributions

The main contributions of this thesis are as follows:

 A new method based on particle filter is developed to estimate the IFB rate according to the RATIO

distribution that has relatively larger values corresponding to the pre-defined IFB rate samples which

are closer to the true IFB rate value. This approach can estimate the IFB rate accurately with short

convergence time and can track it online reliably even without known station coordinates. It can be

applied for precise IFB rate calibration with long data sets, or for fast and real-time calibration with

convergence judged by standard deviation (STD) value.

 To satisfy the requirements on high precision in special cases, the RPF is introduced and hence the

noise in the prediction model can be set to small values or even value zero to achieve more precise

results. Besides, relating the number of particles to the STD values leads to the reduction of number of

particles in the tracking process after convergence and therefore the computation time for each epoch is

reduced.

 The particle filter method is utilized to estimate the phase F-ISB in multi-GNSS data processing. The

characteristics of RATIO distribution with different F-ISB values are investigated firstly, which shows

that the correct F-ISB values lead to relatively larger RATIO values and therefore the F-ISB can be

estimated with the approach similar as the phase IFB estimation. Besides, the half-cycle problem caused

by the periodic characteristic of ISB is solved by a cluster analysis method. In addition, the inter-system

models with different frequencies are employed in the multi-GNSS integration and their integer DD-

ambiguities are investigated for the first time.

 The multi-dimensional approach is demonstrated. The two-dimensional particle filter approach, which

can estimate two F-ISB values simultaneously in the case of observations from three systems, is taken

as a case study. With even only two satellites from each system, the F-ISB values can still be

determined and the fixed solutions are estimated successfully.

 The applications of estimated IFB rate in GLONASS data processing and F-ISB in multi-GNSS data

processing in precise positioning are investigated with short baselines. With the IFB rate estimated by

the proposed approach, the DD-ambiguities can be fixed as integers in real-time and the corresponding

fixed solutions are consistent with GPS fixed solutions. With the estimated F-ISB rate, both intra- and

inter-system DD-ambiguities can be fixed as integers and therefore, the fixed solution has larger chance

to be available especially in the environments with fewer satellites from each constellation.

1.4 Outline

The eight chapters of this thesis are organized as below.

Chapter 1 introduces the background, objective and main contribution of this work, as well as the outline of this

thesis.

Chapter 2 gives a brief introduction to the existing satellite navigation systems, and presents the mathematic

models and the procedures in GNSS data processing.

Chapter 3 presents the principles of particle filter. Bayesian filtering, Kalman filter, bootstrap particle filter and

RPF are introduced.

Chapter 4 proposes a new approach based on particle filter to estimate the carrier phase IFB rate. Afterwards, the

estimation approach based on RPF and the procedure relating the number of particles to the STD are presented.

Chapter 5 utilises the approach based on particle filter to estimate the F-ISB with the half-cycle problem solved

by a cluster analysis method. Later on, the long term characteristics for F-ISBs are investigated.

Chapter 6 demonstrates the multi-dimensional particle filter approach. The two-dimensional approach which can

estimate two F-ISB values between three systems simultaneously is implemented.

Chapter 7 is dedicated to the applications of IFB rate and F-ISB which are estimated by the proposed method.

Finally, chapter 8 draws conclusions of the work and provides outlook for the future research.

2 Multi-GNSS Data Processing

The different GNSS have similar principles. In all the systems, the ranging signal emitted by a navigation

satellite usually contains the PRN code, the navigation data message and the radio frequency carrier. The PRN

code carrying the time information of the satellite clock is firstly combined with the binary navigation data

message which provides data about the satellite orbit, clock correction and is modulated on the radio frequency

carrier. The signals are then received by the receivers which generate observables such as code pseudorange and

carrier phase. For multi-GNSS receivers, signals from more than one constellation can be received. The

observations from different systems can be processed together, which can improve the accuracy, reliability and

availability of GNSS precise positioning.

This chapter aims to introduce the multi-GNSS status and the observation models in data processing. The GNSS

satellite constellations and the satellite frequencies are presented in section 2.1, followed by the observable

models and error sources, as well as the relative positioning models including the IFB encountered in GLONASS

data processing and the ISB among different systems in section 2.2 and 2.3. The parameter estimation in GNSS

data processing is described in section 2.4.

2.1 Satellite Constellations and Signals

GPS

GPS built by US department of Defence is the earliest system, which is fully operated since 1995, and now it is

operated and maintained by the US Air Force. At the beginning, GPS had the Selective Availability (SA)

strategy to reduce the accuracy for civil users, which was stopped in 2000 and the new GPS satellites have no SA

function. GPS employs CDMA technology to identify satellites (Hofmann-Wellenhof et al. 2007).

GPS is composed of 24 satellites which are scattered in six orbit planes with inclination angles of 55 degrees and

with the orbit altitude of approximate 20200 km. The GPS orbit period is a half sidereal day or about 11 hours 58

minutes. In each orbit plane there are four satellites, which are not evenly spaced but with angles of 30 105 120

and 105 degrees so that at least six satellites are visible almost everywhere from the earth’s surface. After a

constellation expansion in June 2011, three extra satellites are added to the constellation and other satellites are

repositioned. As a consequence, there are 27 satellites in the constellation and the coverage in most parts of the

world is improved. In the orbit planes, there are usually extra satellites which are not considered to be part of the

core constellation, such as on 5th June 2015, 31 operational satellites (IAC_GPS 2015) in total are in orbit. Due

to the new demands and technology, the GPS satellites are being modernized since 1999 by employing a new

type of satellites and a new operational control system. The GPS satellites in orbits have evolved from Block I to

Block IIF and they will be replaced gradually by GPS III in the future to maintain the constellation and to

improve the services.

Table 2.1 Satellites of GPS on 5th June 2015

Satellite type

# Launched

satellites

# Operational

satellites

Signal band

Code

Block I

L1 L2

C/A P(Y)1 P(Y)2

Block II

L1 L2

C/A P(Y)1 P(Y)2

Block IIA

L1 L2

C/A P(Y)1 P(Y)2

Block IIR

L1 L2

C/A P(Y)1 P(Y)2

Block IIR(M)

L1 L2

C/A P(Y)1 P(Y)2

L2C

Block IIF

L1 L2 L5

C/A P(Y)1 P(Y)2

L2C L5C

GPS III

L1 L2 L5

C/A P(Y)1 P(Y)2

L2C L5C L1C

The terrestrial reference system of GPS is the World Geodetic System 1984 (WGS-1984), its time system is

related to atomic time and referenced to coordinated universal time (UTC). The GPS time has a constant offset

with international atomic time (TAI), which is referred to as leap seconds. To reduce the major error caused by

the ionospheric refraction, GPS has two carrier frequencies including L1 (1575.42 MHz) and L2 (1227.60 MHz).

Moreover, new carrier frequency L5 (1176.45 MHz) is being introduced in the GPS modernization program

(Hofmann-Wellenhof et al. 2007, NOAA, 2015). The operational satellites on 5th June 2015 are presented in

Table 2.1 (IAC_GPS 2015).

GLONASS

The development of GLONASS system was started in 1976 by the former Union of Soviet Socialist Republics

(USSR). The GLONASS constellation was completed in 1996. Now GLONASS is operated by Russian Space

Forces. Even though the constellation kept declining in late 1990s and had only six to eight satellites at the worst

time, it recovered to full constellation again in October 2011. GLONASS stations have been widely established

until now, but the usage of GLONASS is still far from that of GPS. One of the reasons is the employed FDMA

technique, which leads to different wavelengths of satellites and therefore it is not easy to fix the integer

DD-ambiguities (Wang et al. 2001, Meindl 2011). In the GLONASS modernization, CDMA signals will be

added.

GLONASS has three orbit planes whose ascending nodes are separated by 120 degrees and their inclination

angle to the equator is 64.8 degrees. The eight satellites in each orbit plane are evenly spaced so that at least five

satellites can be observed over more than 99% of the earth’s surface (Hofmann-Wellenhof et al. 2007). The orbit

period of GLONASS satellites is 11 hours 15 minutes 44 seconds and the full constellation includes 24 satellites

with the altitude of about 19100 km. The satellite types of GLONASS have evolved from GLONASS Block I to

GLONASS-M. At present, all 24 satellites belong to GLONASS-M. In the future, GLONASS-K1 will be

launched. GLONASS will have its own CDMA carrier frequencies and even CDMA carrier frequencies

overlapping with these of GPS (Stupak 2010). The satellites of GLONASS on 5th June 2015 are presented in

Table 2.2 (IAC_GLO 2015), where frequency bands G1, G2 and G5 overlap with GPS carrier frequency bands.

The terrestrial reference system of GLONASS is the Earth 1990 (PE-90 or in Russian PZ-90). It can be

transformed to WGS-84 by a seven parameter transformation. The time system for GLONASS has a strong

relationship with UTC. Except for the three-hour difference due to the difference between Moscow Time and

Greenwich Time, the remaining difference between GLONASS time system and UTC is less than 1 millisecond.

Table 2.2 Satellites of GLONASS on 5th June 2015

Satellite version

# Launched

satellites

# Operational

satellites

Signal band

Code

FDMA

CDMA

GLONASS Block I

L1 L2

L1OF L1SF L2SF

GLONASS Block IIA

L1 L2

L1OF L1SF

L2OF L2SF

GLONASS Block IIB

L1 L2

L1OF L1SF

L2OF L2SF

GLONASS Block IIV

L1 L2

L1OF L1SF

L2OF L2SF

GLONASS-M

L1 L2

L3(launched

since 2014 )

L1OF L1SF

L2OF L2SF

L3OC(launched

since 2014 )

GLONASS-K1

L1 L2

L1OF L1SF

L2OF L2SF

L3OC

GLONASS-K2

L1 L2

L1 L2 L3

L1OF L1SF

L2OF L2SF

L1OC L1SC L2OC

L2SC L3OC

GLONASS-KM

L1 L2

L1 L2 L3

G1 G2 G5

L1OF L1SF

L2OFL2SF

L1OC L1SC L2OC

L2SC L3OC L3SC

L1OCM L3OCM

L5OCM

As mentioned before, GLONASS employs FDMA to identify satellites, resulting in different frequencies for

satellites. The frequencies for each frequency band can be expressed by

𝑓𝑗,𝑘=∆𝑓𝑗(2848+𝑘), (2.1)

where 𝑗=1,2 is the frequency band number; 𝑘 is the frequency channel number; ∆𝑓𝑗 is the frequency increment

for two adjacent channels within the same frequency band (ICD-GLONASS 2008), which are 0.5625 MHz and

0.4375 MHz, respectively. The two frequency bands can be expressed by

𝑓1,𝑘=1602.00+0.5625 𝑘, (2.2a)

𝑓2,𝑘=1246.00+0.4375 𝑘. (2.2b)

Although the frequency bands were wider in the past because the frequency numbers were from 0 to 24, the

satellites launched after 2005 have frequency numbers only from -7 to 6 to avoid interfering with radio

astronomy signals and signals of satellite communication services.

Galileo

Galileo is initiated by the European Commission (EC) and the European Space Agency (ESA), and is still under

development at present. Just as GPS and GLONASS, Galileo is expected to be another GNSS but is designed to

provide the services of highest precision and intended to be operated as a civilian GNSS. The construction of

Galileo was started in 2011 and the constellation is expected to be completed by 2019. Until July 2015, eight

satellites have been launched (Lekkerkerk 2015) and three of them were in operation.

Galileo will be composed of 27 operational and 3 spare satellites in three nearly circular medium earth orbit

(MEO) planes with inclination angle to the equator of 56 degrees. In each orbit plane, nine operational satellites

with altitude of 23222 km will be evenly distributed (Nurmi et al. 2015). The orbit period is 14 hours 4 minutes

45 seconds and the constellation configuration is repeated every ten days. The Galileo program has two phases,

the In-Orbit Validation (IOV) phase and the Full Operational Capability (FOC) phase. After the constellation is

fully developed, six to eight satellites will be always visible at most locations on Earth’s surface.

The coordinate system of Galileo is the Galileo terrestrial reference frame (GTRF), which is very close to the

International terrestrial reference frame (ITRF) and the three-dimensional differences is within 3 centimetres.

The time reference of Galileo is the Galileo system time (GST), which is a continuous atomic time scale and has

a constant offset with the TAI. The carrier frequency bands of Galileo include E1 (1575.12 MHz), E6

(1278.75 MHz), E5a (1176.450 MHz) and E5b (1207.140 MHz) (ICD-Galileo 2015). Signal bands for each

satellite type are shown in Table 2.3.

Table 2.3 Satellites of Galileo on 5th June 2015

# Launched

satellites

# Operational

satellites

Signal band

Code

Galileo IOV

E1 E6 E5a E5b

E1A E1B E1C E6A E6B

E6C E5a-I E5a-Q E5b-I

E5b-Q

Galileo FOC

E1 E6 E5a E5b

E1A E1B E1C E6A E6B

E6C E5a-I E5a-Q E5b-I

E5b-Q

BDS

The Chinese BDS began to offer services in Asia-Pacific region since December 2012 and will provide global

services from 2020. Until August 2015, 16 navigation satellites have been launched. BDS is the first

constellation that enables a systematic assessment of three frequencies (Montenbruck et al. 2013).

BDS is planned to be composed of 27 MEO satellites, 5 Geostationary Orbit (GEO) satellites and 3 Inclined

Geosynchronous Satellite Orbit (IGSO) satellites. The MEO satellites have an inclination angle of 55 degrees

which is the same as IGSO satellites, but the altitudes of the MEO and IGSO satellites are 21528 km and 35786

km, respectively. 24 MEO satellites will be evenly distributed in the orbit planes and other three satellites will be

spare. The GEO satellites have altitude of 35786 km with the inclination angle of 0 degrees (ICD-BDS 2013).

Until June 2015, the regional system composed of 5 GEO satellites, 5 IGSO satellites and 4 MEO satellites has

been completed. At present, satellites for its global system are being launched. The three carrier frequencies of

BDS satellites are B1 (1561.098 MHz), B2 (1207.14 MHz) and B3 (1268.52 MHz).

QZSS

QZSS is operated by Japan’s Cabinet Office. QZSS aims to augment the GPS services in Japan and to maximize

the interoperation ability with GPS (IS-QZSS 2014). The system will comprise 4 satellites between 2018 and

2022, and 7 satellites after 2023. Until June 2015, there was only one QZSS satellite in orbit.

In the system of four satellites, three satellites will deploy evenly in a quasi-zenith orbit with inclination angle to

the equator of 43±4 degrees so that at least one satellite will be located near the zenith at any moment; the rest

one is a GEO satellite on equator. Its coordinate system is the Japanese satellite navigation Geodetic System

(JGS) and the offset between QZSS coordinate system and WGS-84 will be less than 2 cm. QZSS time scale is

aligned to the TAI and has the same integer-second offset to TAI as for GPS. The difference between time offset

scales of QZSS and GPS is less than 2 m in distance and is emitted in the navigation message (Hofmann-

Wellenhof et al. 2007). The satellites of QZSS transmit almost the same signals as these of the GPS satellites.

Except for L1, L2 and L5 carrier frequencies which are the same as GPS, there is another frequency which is the

same as Galileo E6, modulated by an experimental signal LEX for high precision service (3 cm).

IRNSS

IRNSS is a regional satellite navigation system developed by Indian Space Research Organization (ISRO) which

is under the control of the Indian government. It provides both civilian and military services (ICD-IRNSS 2014).

The system is supposed to have 7 satellites including 3 GEO satellites and 4 IGSO satellites, and the IGSO

satellites have an inclination angle of 29 degrees. All 7 satellites will be fully operational since 2016. Until

August 2015, four satellites have been in orbit. IRNSS has two carrier frequencies, SPS-L5 and SPS-S

(ICD-IRNSS 2014).

Among these systems, some frequency bands overlap with each other or are close to each other. The frequencies

of above systems are summarized in Table 2.4, where (F) refers to FDMA signal while (C) means CDMA signal

for GLONASS frequency bands.

Table 2.4 Frequency distribution of the satellite systems

Central Frequency

(MHz)

Wavelength

(cm)

GPS

Galileo

BDS

QZSS

IRNSS

GLONASS

2492.028

12.03

SPS-S

1602.000

18.71

L1(F)

1575.420

19.03

1561.098

19.20

1278.750

23.44

LEX

1268.520

23.63

1246.000

24.06

L2(F)

1227.600

24.42

1207.140

24.83

E5b

1202.025

24.94

L3(C)

1191.795

25.15

1176.450

25.48

E5a

SPS-L5

2.2 GNSS Observables and Error Sources

GNSS receivers usually provide code pseudorange and carrier phase observables. The observations depend not

only on the geometric distance between the satellite and receiver, but also on the error sources of the instruments

and on the path of signal transmission. For example, the clocks on the satellite and at the receivers are not strictly

synchronized i.e. clock biases; the atmosphere changes the signal propagation speed and the path; the multipath

effects interfere with the signal reception; the hardware delays both on the satellite and at the receivers are

included in the measurements. In order to obtain the accurate distance values, these error sources must be

carefully handled in data processing.

2.2.1 GNSS Observables

According to (Teunissen 1996), the code pseudorange observations can be modelled as

𝑃𝑎𝑖=𝜌𝑎𝑖−𝑐(𝛿𝑡𝑎−𝛿𝑡𝑖)+𝑑𝑎

𝑖−𝑑𝑖+𝐼𝑎𝑖+𝑇𝑎𝑖+𝑅𝑎

𝑖+𝑆𝑎𝑖+𝑀𝑎𝑖+𝜀𝑎𝑖, (2.3)

where 𝑃 is the code pseudorange measurement; 𝑖 and a refer to the satellite number and the observation station

number, respectively; 𝜌 is the initial value of the distance; 𝛿𝑡𝑎 and 𝛿𝑡𝑖 are the receiver clock bias and the satellite

clock bias. 𝑑𝑎 and 𝑑𝑖 are the receiver hardware delay and the satellite hardware delay in code observations; 𝐼 is

the ionospheric delay; 𝑇 the tropospheric delay; 𝑅 is the effects of the relativity; 𝑆 is the sagnac effect i.e. earth

rotation correction; 𝑀 refers to the multipath effects on the code pseudorange measurement; 𝜀 denotes the

remaining errors which are considered as white noise.

The distance 𝜌𝑎𝑖 in (2.3) is calculated from the initial coordinates of the satellite and the station. Assuming the

antenna phase centres of satellite i and station a are (𝑥𝑖 𝑦𝑖 𝑧𝑖) and (𝑥𝑎 𝑦𝑎 𝑧𝑎) in Earth-fixed frame at observation

time, respectively, as well as considering the tide effects, the distance in (2.3) can be expressed by

𝜌𝑎𝑖=√(𝑥𝑖−𝑥𝑎)2+(𝑦𝑖−𝑦𝑎)2+(𝑧𝑖−𝑧𝑎)2 +𝑇𝑖𝑑𝑒𝑎𝑖, (2.4)

where 𝑇𝑖𝑑𝑒𝑎𝑖 refers to the effect of the solid earth tide and ocean loading. Besides, the antenna phase centre

coordinates can be converted to the mass centre coordinates of the satellite or the station reference point

coordinates by antenna eccentricity correction and antenna phase centre correction. The satellite coordinates,

either antenna phase centre coordinates or mass centre coordinates, are normally seen as known values in

positioning as they can be usually determined from either the broadcast ephemerides or the precise ephemerides.

However, the station coordinates are usually unknown in positioning, but can be estimated by an iterative

calculation after the linearization of (2.4).

The carrier phase has a similar observation equation as the code pseudorange, as they propagate at the same time

on the same path. However, the ionospheric delay has a minus sign, because the ionosphere accelerates the phase

propagation speed. Thus, the GNSS phase observation model can be expressed in distance as (Teunissen 1996)

𝜆𝑖𝛷𝑎𝑖=𝜌𝑎𝑖−𝑐(𝛿𝑡𝑎−𝛿𝑡𝑖)+𝜇𝑎

𝑖−𝜇𝑖+𝜆𝑖𝑁𝑎𝑖+𝜆𝑖𝜓𝑎𝑖−𝐼𝑎𝑖+𝑇𝑎𝑖+𝑅𝑎

𝑖+𝑆𝑎𝑖+𝑚𝑎

𝑖+𝜉𝑎𝑖, (2.5)

where 𝛷 is the carrier phase measurement; 𝜆 is the wavelength; 𝜓 is the initial phase value; 𝜉 is the noise on

carrier phase observations; 𝜇𝑎

𝑖 and 𝜇𝑖 are the receiver and satellite hardware delays in phase, respectively; 𝑁 is

the phase ambiguity which is an integer number; 𝑚 refers to the multipath effects.

The signal can usually be measured at the accuracy of better than 1% wavelength. The equivalent wavelength of

the code signal is very long, such as around 300 m for GPS, which limits the accuracy of the positioning results

with code pseudorange observations. The carrier frequency has a much smaller wavelength, such as around

20 cm for GPS. Therefore, the carrier phase observations enable much more accurate positioning.

All the terms in the observation equations must be carefully handled, especially for precise positioning with

carrier phase observations. Usually, some of the terms can be eliminated or significantly reduced by differencing

the observations such as in relative positioning; some can be corrected with physical models; some must be

estimated as unknown parameters.

2.2.2 GNSS Error Sources

Time Clock Offset

The range measurements between satellites and receivers are based on the signal propagation time which is

calculated with the emission time on the satellite and the reception time at the receiver with reference to the

satellite and receiver clocks. Since the clocks are not strictly synchronized, the clock biases are included in the

observations.

The satellite clocks are atomic clocks which are very stable, for example with stability of 10-13 to 10-14 over one

day for GPS rubidium clocks and 10-14 to 10-15 over one day for GPS hydrogen clocks. The satellite clock bias is

monitored and predicted, then broadcasted to users so that it can be corrected in the data processing. Much more

accurate satellite clock corrections can be calculated using data of a GNSS network. For example, the accuracy

of the IGS final precise ephemeris can reach 0.1 ns (Dow et al. 2009), and that of the real-time service can reach

0.28 ns for GPS and 0.82 ns for GLONASS (Hadas and Bosy 2015).

However, the receivers are usually equipped with quartz crystal clocks which are not stable and drift fast

compared to the atomic clocks. Hence, the receiver clock biases are usually estimated every epoch if they cannot

be eliminated. Almost all the modern receivers have only one receiver clock for signals of all the satellite

systems (Melgard et al. 2013). This indicates that, without consideration of the hardware delays, the clock biases

for code and carrier phase observations from the same receiver including observations from different systems

and different frequency bands can be considered to be the same. Therefore, differencing the observations

between satellites can eliminate receiver clock biases.

The Hardware Delay and Initial Phase

Signals are delayed when they travel through hardware in devices, leading to hardware delays which can vary

with signals due to different hardware paths. The differences of hardware delays may exist between different

GNSS systems, as well as between different channels of GLONASS signals due to FDMA technology. The

different hardware delays between GLONASS channels lead to IFB, while the hardware delays between

channels of different satellite systems lead to ISB (Zinoviev 2005, Odijk and Teunissen 2013a, Wanninger and

Wallstab-Freitag 2007).

Even with the same frequency, signals from different constellations may encounter different delays at DSP step

in the device as they may be processed differently in the firmware. These delays are deterministic and not likely

to be affected by the environment (Melgard et al. 2013). Although the hardware delays and the delays caused by

the DSP may have different characteristics, they actually cannot be separated and lump together resulting in ISB.

Thus, they are not distinguished in this investigation and are referred to as hardware delays.

The initial phase bias is caused by the non-synchronisation of the satellite and receiver clocks at the first epoch

of an observation session. This bias is constant for a continuous observation session and only the fractional part

is important because the part of integer multiple of wavelength is absorbed by the integer ambiguities. The

receiver designers are supposed to make sure that all the initial phase biases are the same for all tracking

channels of one satellite system (O'Driscoll 2010). O’Driscoll (2010) and Wang et al. (2001) demonstrated that

the initial phases are the same if the same heterodyne for all systems is used in a receiver, which is usually the

case for observations of the same system. Besides, the initial phases may change when there is a power off

(Kozlov and Tkachenko 1998). However, the initial phase does not change during a continuous observation

session and is correlated with the hardware delay, and hence can be considered to lump with the phase hardware

delay to avoid being parameterized, separately.

Atmospheric Delays

According to the electromagnetic structure, the atmosphere is divided into the troposphere (actually the neutral

atmosphere) and the ionosphere. Both of them affect the observations.

The troposphere layer extends to the height of around 12 km and contains about 80% of the mass of the earth’s

atmosphere. The tropospheric delay does not vary with the radio frequency, i.e. is non-dispersive with respect to

the radio frequency. The delay can be separated into two parts, the hydrostatic delay and the wet delay. The

former is caused by dry air and the latter is caused by water vapour. The hydrostatic delay contributes around 90%

of the tropospheric delay and can be modelled accurately, while the wet delay contributes only 10%. The wet

delay depends on the amount of the water vapour and is difficult to be modelled due to the high variability of the

water vapour distribution. Some models have been proposed to correct the tropospheric delay, such as Hopfield

model and Saastamoinen model (Saastamoinen 1973, Janes et al. 1991). Besides, both the hydrostatic delay and

the wet delay have smallest value for paths oriented along the zenith direction and they increase as the elevation

angle decreases. The delay along a path of arbitrary elevation can be modelled as the zenith delay multiplies a

mapping function which describes the dependence on the elevation angle. Therefore, for long baselines where

the atmospheric delays are not well known, the zenith delays can be estimated with the mapping function so that

the effects of the atmospheric delays can be removed or largely reduced (Bevis et al. 1992, Ge et al. 2000).

The ionosphere is composed of electrons and electrically charged atoms and molecules. Mainly due to the

ultraviolet radiation, the electrons and ions are energetic and separate from each other. They significantly affect

the radio wave propagation. Because the electrons contribute much more delays than the ions, the total electron

content is usually employed to model the ionospheric delay. The ionosphere has different effects to the phase

velocity and group velocity, which leads to delays of the same magnitude but different signs in code

psuedorange and carrier phase measurements. Thanks to the frequency-dependent characteristic, the ionospheric

delay can be calculated with multi-frequency observations. These observations of different frequencies can also

compose the ionosphere-free combinations. For example, the ionosphere-free combination with GPS L1 and L2

frequencies can eliminate the first-order ionospheric delay which accounts for more than 99.9% of the total value

(Hernández-Pajares et al. 2011).

Multipath

Multipath is mainly caused by the environment around a GNSS receiver. When the signal is reflected by some

surfaces near the receiver, such as a wall surface, water surface, ground surface and so on, the receiver will

receive both the direct and the reflected signals from the same satellite. In this case, the signals cannot be

separated and the reception time recorded in the receiver is biased. The effects of multipath depend on the

properties of the surfaces and their distances to the antenna, and are different for code and carrier phase. The

code is more likely to be affected by multipath and the introduced error can amount to 10-20 m, whereas the

effect in carrier phase is smaller than one quarter of the wavelength, i.e. below 5 cm. This value is usually below

1 cm under the conditions of good satellite geometry and properly long observation interval (Hofmann-

Wellenhof et al. 2007).

Multipath effects can be reduced by careful antenna design, improved receiver technology or signal and data

processing procedure. But the multipath is difficult to be eliminated once it has happened as it depends on the

surroundings, especially in the case of kinematic positioning. Because the satellite with lower elevation angle is

more likely to suffer from multipath effects, setting a larger elevation mask value in the data processing to

exclude the satellites with lower elevation angles is very helpful to reduce the multipath effects.

Others

Besides, the solid earth tide and ocean tide loading, which are mainly caused by the moon and the sun, displace

the station positions. The relativity including the special and general relativity, as well as the sagnac effects also

affects the range measurements. Fortunately, all these effects can be corrected by models to the extent which can

be neglected even for precise poisoning. Moreover, the phase wind-up effect, which is caused by the orientations

of both receiver and satellite antennas, must be considered for PPP and relative positioning of long baselines, but

this effect can be neglected in relative precise positioning with baselines shorter than 100 km.

2.3 Relative Positioning Models

2.3.1 General Relative Positioning Models

Since the inter-system difference in space and time reference frame can be aligned with known parameters or

united in integrated orbit determination, they will not be considered in the following models. Besides, the effects

of antenna eccentricity centre, antenna phase centre, solid earth tide, ocean tide loading, relativity, sagnac, as

well as the phase wind-up can be accurately modelled, so they will not be listed in the following equations. The

multipath effects at two stations are usually not the same, but they are difficult to estimate. Hence, it is assumed

that the stations have taken measures to mitigate the multipath effects and the remaining parts are seen as white

noise.

Then the SD-model between two stations 𝑎 and 𝑏 from the same satellite 𝑖 for the same carrier frequency can be

derived from (2.3) and (2.5) by differencing equations between stations (Teunissen and Kleusberg 1996,

Hofmann-Wellenhof et al. 2007), and expressed by

𝑃𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖−𝑐𝛿𝑡𝑎𝑏+𝑑𝑎𝑏

𝑖+𝐼𝑎𝑏

𝑖+𝑇𝑎𝑏

𝑖+𝜀𝑎𝑏

𝑖, (2.6a)

𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖−𝑐𝛿𝑡𝑎𝑏+𝜇𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖+𝑇𝑎𝑏

𝑖+𝜆𝑖𝑁𝑎𝑏

𝑖+𝜆𝑖𝜓𝑎𝑏

𝑖+𝜉𝑎𝑏

𝑖. (2.6b)

The satellite clock biases are eliminated in model (2.6). The hardware delay 𝑑𝑎𝑏

𝑖 includes not only the code

biases caused by receivers, but also satellite differential code biases if different code types and tracing modes are

implemented in the two receivers. For the same frequency of satellite 𝑖, the hardware delay 𝜇𝑎𝑏 includes only the

differential delays between receivers as the hardware delays in satellite are eliminated. The initial phase 𝜓𝑎𝑏

𝑖 may

be different between the two receivers and therefore it also remains.

By differencing equations between satellites, the SD-model between two satellites 𝑖 and 𝑗 at the same

observation station 𝑎 is

𝑃𝑎𝑖𝑗=𝜌𝑎𝑖𝑗+𝑐𝛿𝑡𝑖𝑗+𝑑𝑎

𝑖𝑗+𝐼𝑎𝑖𝑗+𝑇𝑎𝑖𝑗+𝜀𝑎𝑖𝑗, (2.7a)

𝜆𝑗𝛷𝑎𝑗−𝜆𝑖𝛷𝑎𝑖=𝜌𝑎𝑖𝑗+𝑐𝛿𝑡𝑖𝑗+𝜇𝑎

𝑖𝑗+𝜆𝑗𝑁𝑎𝑗−𝜆𝑖𝑁𝑎𝑖+𝜆𝑗𝜓𝑎𝑗−𝜆𝑖𝜓𝑎𝑖−𝐼𝑎𝑖𝑗+𝑇𝑎𝑖𝑗+𝜉𝑎𝑖𝑗. (2.7b)

In model (2.7), the receiver clock biases cancel each other but the satellite clock biases remain. If the two

observations are identical in frequency, the two integer SD-ambiguities can be merged together to form an

integer DD-ambiguity. If a common heterodyne is applied in the receiver, 𝜓𝑎𝑗 will equal 𝜓𝑎𝑖.

Differencing SD-model (2.6) between satellites, or differencing SD-model (2.7) between stations, the DD-model

between two stations 𝑎 and 𝑏, and two satellites 𝑖 and 𝑗 is

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝑑𝑎𝑏

𝑖𝑗 +𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑏

𝑖𝑗, (2.8a)

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏

𝑖𝑗 +𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖+𝜆𝑗𝜓𝑎𝑏

𝑗−𝜆𝑖𝜓𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.8b)

The DD-model (2.8) can be constructed with observations of the same system or different systems, with the

same carrier frequency or different carrier frequencies. As described in section 2.2.2, although the initial phases

may be different for observations of different systems, the initial phase values do not change during a continuous

observation session and cannot be separated from hardware delays. Hence, they can be considered to lump with

the phase hardware delays to avoid being parameterized, separately.

The stochastic model for (2.8) is constructed with consideration of two aspects, the variances of the raw GNSS

observations and the correlation of the observation combinations. The variances of the raw GNSS observations

are dominated by the systematic errors, such as the multipath effects, the remains of the atmospheric delays, and

have a close connection with the satellite elevation angles (Jin and Wang 2004). Thus, a function of satellite

elevation angle is usually used to describe the variances of the raw GNSS observations. Due to the complexity of

unknown factors affecting the variances of the raw observations, the function can only be expressed as an

approximation, such as via sine or cosine functions. This study employs the function with the sine of the

elevation angle (King and Bock 1999, Jin and Wang 2004). The variance of a raw observation is calculated by

𝜎2=𝑎2+𝑏2/sin 2(𝐸𝑙), (2.9)

where 𝜎 is the variance of the raw observation; 𝑎 and 𝑏 are constant values, which are different for code

pseudorange and carrier phase observations; 𝐸𝑙 is the elevation angle of the satellite. The correlation of the

observations includes the physical correlation and the mathematical correlation (Hofmann-Wellenhof et al. 2007).

The physical correlation is caused by the fact that some observations are emitted or received by the same device

and is usually not considered in practice. This indicates that the raw observations are regarded as independent

from each other. The mathematical correlations introduced by computation of differences need to be modelled

according to the rules of variance propagation.

2.3.2 Intra-System Models

Intra-System DD-Model with the Same Frequency

The intra-system DD-model with satellites of the same frequency can eliminate not only the receiver and satellite

clock biases, but also the hardware delays as described in 2.2.2. After removing these terms in model (2.8), the

code pseudorange and carrier phase observations are expressed by

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑏

𝑖𝑗, (2.10a)

𝜆𝛷𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝜆𝑁𝑎𝑏

𝑖𝑗 −𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.10b)

Intra-System DD-Model of GLONASS

Due to the FDMA technology in GLONASS, the hardware delays on both code pseudorange and carrier phase

DD-model observations cannot be eliminated at least in the case of employing receivers of different types. The

difference in frequency results in that the two unknown SD-ambiguities for the two satellites cannot merge

together directly to form an integer DD-ambiguity in (2.8). So the code pseudorange and carrier phase

DD-observations for GLONASS can be expressed by

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝑑𝑎𝑏

𝑖𝑗 +𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑏

𝑖𝑗, (2.11a)

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏

𝑖𝑗+𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.11b)

The separated SD-ambiguities lead to a rank-deficiency of the normal equations (NEQ) in data processing,

because the unknown SD-ambiguity parameters is one more than the phase DD-equations. There are three

approaches to deal with this problem (Leick 1998). The first one is to rewrite the SD-ambiguity terms 𝜆𝑗𝐵𝑎𝑏

𝑗−

𝜆𝑖𝐵𝑎𝑏

𝑖 into two terms, the term including integer DD-ambiguity and the term including SD-ambiguity, as shown

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏

𝑖𝑗 +(𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑗𝑁𝑎𝑏

𝑖)+(𝜆𝑗𝑁𝑎𝑏

𝑖−𝜆𝑖𝑁𝑎𝑏

𝑖)−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀

=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏

𝑖𝑗 +𝜆𝑗(𝑁𝑎𝑏

𝑗−𝑁𝑎𝑏

𝑖)+(𝜆𝑗−𝜆𝑖)𝑁𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀

=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏

𝑖𝑗 +𝜆𝑗𝑁𝑎𝑏

𝑗𝑖 +(𝜆𝑗−𝜆𝑖)𝑁𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.12)

Because the frequency difference in GLONASS is small within the same frequency band, 𝜆𝑗−𝜆𝑖 is seen as a

small value. The maximum values of 𝜆𝑗−𝜆𝑖 for L1 and L2 are 0.85 and 1.10 mm, respectively. If the error in

𝑁𝑎𝑏

𝑗𝑖 caused by inaccurate 𝑁𝑎𝑏

𝑖 is smaller than 0.1 cycles, according to the coefficients in (2.12), the biases in 𝑁𝑎𝑏

𝑖

should be smaller than 4.1 m and 5.3 m for L1 and L2, respectively. After the initial value of the SD-ambiguity

𝑁𝑎𝑏

𝑖 is calculated, the term (𝜆𝑗−𝜆𝑖)𝑁𝑎𝑏

𝑖 can be removed from model (2.12). Then model (2.12) has the same

form as (2.10b) except for the IFB parameters.

The second approach is to transform the SD-model (2.6b) in distance into the model in cycle as

𝛷𝑎𝑏

𝑖=1

𝜆𝑖(𝜌𝑎𝑏

𝑖−𝛿𝑡𝑎𝑏𝑐)+𝑁𝑎𝑏

𝑖+1

𝜆𝑖(𝜇𝑎𝑏−𝐼𝑎𝑖+𝑇𝑎𝑖)+1

𝜆𝑖𝜉𝑎𝑏

𝑖. (2.13)

The coefficient of the SD-ambiguity disappears in (2.13) and therefore two SD-ambiguities from two stations

can be merged together directly, but the receiver clock bias terms cannot cancel each other in this case as their

coefficients are not the same. The DD-model can be expressed by

𝛷𝑎𝑏

𝑖𝑗 =1

𝜆𝑗(𝜌𝑎𝑏

𝑗−𝜌𝑎𝑏

𝑖)−(1

𝜆𝑗−1

𝜆𝑖)𝛿𝑡𝑎𝑏𝑐+𝑁𝑎𝑏

𝑗𝑖 +1

𝜆𝑗(𝜇𝑎𝑏−𝐼𝑎𝑗+𝑇𝑎𝑗)

−1

𝜆𝑖(𝜇𝑎𝑏−𝐼𝑎𝑖+𝑇𝑎𝑖)+1

𝜆𝑗𝜉𝑎𝑏

𝑗−1

𝜆𝑖𝜉𝑎𝑏

𝑖. (2.14)

Even though the ambiguities are supposed to be integers in (2.14), the clock biases have to be known or

estimated.

The third approach is to use equation (2.11) directly. The SD-ambiguities are estimated as unknown parameters

instead of the DD-ambiguities. In this case, the unknown SD-ambiguity parameters are one more than the

number of equations, so the corresponding equation set is rank-deficient. To remove the rank-deficiency, initial

SD-ambiguity values are assigned to the parameters. After the float SD-ambiguities are estimated, the

DD-ambiguities in cycle can be estimated from the SD-ambiguity solutions (Kozlov and Tkachenko, 1997).

The three approaches are essentially the same. The information, which is usually the code pseudorange

measurments, is required by all of the three approaches in different manners. In fact, Li and Wang (2011)

compared the first and second approaches and demonstrated that they have similar performaces. Al-Shaery et al.

(2012) compared the second and third approaches and found that no one obviously surpasses the other one.

However, the third approach is more convinent to deal with the reference satelltes (Kozlov and Tkachenko 1997)

and to combine observations of different systems. Therefore, it is utilised in the data processing of this study.

Another problem caused by different wavelengths is the IFB due to the different hardware delays, at least when

receivers from different manufacturers are employed. The code pseudorange IFB is not so important as the

observations can be down-weighted. However, the carrier phase IFB has to be removed or estimated to recover

the integer nature. Fortunately, the IFB is proved to be linear with the channel number. This indicates that the

carrier phase IFB can be precisely expressed by a constant value and its rate. In the DD-model, the constant part

of IFB is eliminated and only the part with IFB rate is left (Wanninger 2012, Al-Shaery et al. 2013). Hence, the

linear model for carrier phase IFB is

𝛾𝑎𝑏𝑖𝑗=(𝑘𝑗−𝑘𝑖)∆𝛾𝑎𝑏, (2.15)

where 𝛾𝑎𝑏𝑖𝑗 and ∆𝛾𝑎𝑏 are the DD-IFB and the rate, respectively; k is the channel number. Usually, IFB rate has

similar values on both L1 and L2 frequency bands, which can also be utilised in the IFB modelling (Reussner

and Wanninger 2011, Wanninger 2012).

Including the parameter of the relative IFB rate, the carrier phase DD-model (2.11b) can be expressed by

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖+(𝑘𝑗−𝑘𝑖)∆𝛾𝑎𝑏−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.16)

In the float solution, the IFBs lump with DD-ambiguity parameters and therefore can be absorbed by ambiguity

parameters. If the IFB rate ∆𝛾𝑎𝑏 is known, the term (𝑘𝑗−𝑘𝑖)∆𝛾𝑎𝑏 can be removed directly from (2.16).

Consequently, the ambiguity resolution can be applied to estimate the integer DD-ambiguities. If it is unknown,

IFB rate must be estimated along with float ambiguities and other unknowns.

2.3.3 Inter-System Models

Inter-System DD-Model with the Same Frequency

The signals of different systems are received at the same time in a receiver as they share the same receiver

reference clock and therefore the receiver clock biases can cancel each other. Moreover, the SD-ambiguities in

(2.6b) and (2.7b) can merge together directly due to the same wavelength. However, the hardware delays cannot

be eliminated at least when the receivers of different types are employed (Odijk and Teunissen 2013a). So the

DD-model can be expressed by

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝑑𝑎𝑏+𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑖𝑗, (2.17a)

𝜆𝛷𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏+𝜆𝑁𝑎𝑏

𝑖𝑗 −𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑖𝑗. (2.17b)

Inter-System DD-Model with Different Frequencies

The inter-system model with different frequencies between different systems can also be constructed. In this case,

the hardware delays cannot be eliminated. Thus, the inter-system DD-model between systems of different

frequencies can be expressed by

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝑑𝑎𝑏+𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑏

𝑖𝑗, (2.18a)

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜇𝑎𝑏+𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖−𝐼𝑎𝑏

𝑖𝑗 +𝑇𝑎𝑏

𝑖𝑗 +𝜉𝑎𝑏

𝑖𝑗. (2.18b)

If the hardware delay in (2.18) is known, the same procedure described in GLONASS single system data

processing in section 2.3.2 can be adopted to solve (2.18). Besides, GLONASS can also be included in this

model, with IFB removed or set as additional unknown parameters.

The same as in GLONASS data processing, the number of SD-ambiguity parameters is one more than the

number of DD-equations in (2.18) and hence at least one of the initial SD-ambiguities is needed. The effects of

the biases in the initial SD-ambiguities in the DD-model can be roughly seen in

𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖=𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑗𝑁𝑎𝑏

𝑖+𝜆𝑗𝑁𝑎𝑏

𝑖−𝜆𝑖𝑁𝑎𝑏

𝑖

=𝜆𝑗𝑁𝑎𝑏

𝑖𝑗 +(𝜆𝑗−𝜆𝑖)𝑁𝑎𝑏

𝑖, (2.19)

which is similar with the ambiguity-related terms on the right side of (2.12) for GLONASS.

The initial values can be calculated from code pseudorange observations based on SD-model (2.6) and be

expressed by

𝑁𝑎𝑏

𝑖=−1

𝜆𝑖𝑃𝑎𝑏

𝑖+𝛷𝑎𝑏

𝑖+1

𝜆𝑖(𝑑𝑎𝑏−𝜇𝑎𝑏+2𝐼𝑎𝑏

𝑖). (2.20)

Considering the unknown parameters 𝑑𝑎𝑏 and 𝜇𝑎𝑏, the error 𝛿𝑁𝑎𝑏

𝑖 in (2.20) is

𝛿𝑁𝑎𝑏

𝑖=1

𝜆𝑖(𝑑𝑎𝑏−𝜇𝑎𝑏+2𝐼𝑎𝑏

𝑖−𝜀𝑎𝑏

𝑖+𝜉𝑎𝑏

𝑖), (2.21)

where the code pseudorange observation error 𝜀𝑎𝑏

𝑖 is the main part of the error. In the research of GLONASS

ambiguity resolution in precise point positioning (Banville 2016) and the GPS/GLONASS ambiguity resolution

in relative positioning (Wang 2001), with the same clock offset parameters, the hardware delays related terms

include only the carrier phase and pseudorange IFB in the non-difference model. This indicates that the hardware

delays are the same for one system after IFB have been corrected. Therefore, the value 𝑑𝑎𝑏 and 𝜇𝑎𝑏 can be seen

as equivalent to each other in data processing and therefore they cancel each other in (2.20) and (2.21).

2.3.4 Simplified General Relative Positioning Models for Short Baselines

To avoid dealing with tropospheric and ionospheric delays, the investigation is constrained to short baselines, i.e.

shorter than 10-15 km. All the above models can be summarized into a simple general model. Considering the

ISB in multi-GNSS integration and the IFB in GLONASS, the general model can be written as

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +𝑑𝑎𝑏+𝛿𝑎𝑏

𝑖𝑗 +𝜀𝑎𝑏

𝑖𝑗, (2.22a)

𝜆𝑗𝛷𝑎𝑏

𝑗−𝜆𝑖𝛷𝑎𝑏

𝑖=𝜌𝑎𝑏

𝑖𝑗 +𝜆𝑗𝑁𝑎𝑏

𝑗−𝜆𝑖𝑁𝑎𝑏

𝑖+𝜇𝑎𝑏+ (𝑘𝑗−𝑘𝑖)∆𝛾𝑎𝑏+𝜉𝑎𝑏

𝑖𝑗. (2.22b)

The two satellites can be from the same system or different systems, i.e. model (2.22) includes both intra- and

inter-system models. In the case of the same system, the ISB parameter 𝜇𝑎𝑏 will be zero. If GLONASS is not

included, the IFB rate parameter ∆𝛾𝑎𝑏 will be zero.

If only the float solution is needed, both IFB and ISB can lump with the ambiguity parameters and therefore they

do not need to be handled. However, for integer ambiguity resolution, the ISB and IFB parameters have to be

known or estimated.

2.4 Parameter Estimation

Four steps are needed to obtain the fixed solution (Teunissen and Verhagen 2007, Verhagen et al. 2012). Firstly,

the float solution is estimated with OLS method; secondly, the float ambiguities are mapped into integer values,

i.e. integer ambiguity resolution; thirdly, the integer ambiguity candidates are validated; finally, the remaining

parameters are recalculated with the integer ambiguity constraint.

Float solution

The unknown parameters in the carrier phase equation (2.22) include coordinates, ISB, IFB and SD-ambiguities,

for long baselines also atmospheric delay and others. After linearization, the models can be written in the form of

matrices as

𝒗=𝑨𝒙+𝑫𝒃+𝑪𝒚+𝒍 𝑷, (2.23)

where 𝒗 is the vector of observation residuals; 𝒃 is composed of the unknown SD-ambiguities

(𝑁𝑎𝑏

𝑖1,𝑁𝑎𝑏

𝑖2,…,𝑁𝑎𝑏

𝑖𝑛), where 𝑖 is the reference satellite and n is the number of the DD-equations; 𝒚 includes the

unknown carrier phase ISB and IFB rate; vector 𝒙 contains the unknown station coordinate and the remaining

parameters; 𝒍 is the observation vector; 𝑷 is the weight matrix of the observations; 𝑨 is the design matrix

containing the SD of the receiver-satellite unit direction vectors and coefficients of other elements in 𝒙; 𝑫 is the

design matrix transforming SD-ambiguities to DD-ambiguities with elements of zero and the corresponding

wavelength; 𝑪 is the design matrix with elements of zero and the SD of the channel numbers for phase IFB rate

parameter in 𝒚, with elements of zero and 1 for phase ISB parameter in 𝒚.

The NEQ is

[𝑨𝑻𝑷𝑨 𝑨𝑻𝑷𝑫 𝑨𝑻𝑷𝑪

𝑫𝑻𝑷𝑫 𝑫𝑻𝑷𝑪

𝒔𝒚𝒎 𝑪𝑻𝑷𝑪][𝒙

𝒃

𝒚]=[𝑨𝑻𝑷𝒍

𝑫𝑻𝑷𝒍

𝑪𝑻𝑷𝒍]. (2.24)

For simplification, the notation

[𝑵𝒙𝒙 𝑵𝒙𝒃 𝑵𝒙𝒚

𝑵𝒃𝒃 𝑵𝒃𝒚

𝒔𝒚𝒎 𝑵𝒚𝒚][𝒙

𝒃

𝒚]=[𝑾𝒙

𝑾𝒃

𝑾𝒚] (2.25)

is introduced.

Both the IFB and the ISB parameters in model (2.25) lead to rank-deficiency, so the equation cannot be solved.

If inaccurate prior information for IFB rate and ISB are employed to remove this rank-deficiency, the parameter

estimation via (2.25) needs observations of long time to converge. For example, observations of several hours up

to one day are needed for the IFB rate estimation even with the assistance of GPS data, which is not applicable

for IFB online calibration or real-time applications.

But if the elements in y are precisely known, the solution with float SD-ambiguities can be written as

[𝒙



𝒃

]=[𝑵𝒙𝒙 𝑵𝒙𝒃

𝑵𝒃𝒙 𝑵𝒃𝒃]−𝟏[𝑾𝒙−𝑵𝒙𝒚𝒚

𝑾𝒃−𝑵𝒃𝒚𝒚]=[𝑸𝒙𝒙 𝑸𝒙𝒃

𝑸𝒃𝒙 𝑸𝒃𝒃][𝑾𝒙−𝑵𝒙𝒚𝒚

𝑾𝒃−𝑵𝒃𝒚𝒚]. (2.26)

With the help of initial SD-ambiguity values, the SD-ambiguities can be successfully estimated. Afterwards, the

SD-ambiguities and their VC matrix calculated by (2.26) are transformed into DD-ambiguities and the

corresponding VC matrix, because the DD-ambiguities have integer nature and thus can be fixed to integers. The

transformation process can be expressed by

𝒃

𝐷𝐷=𝑬𝒃

, (2.27a)

𝑸𝒃

𝑫𝑫 =𝑬𝑸𝒃𝒃𝑬𝑻, (2.27b)

where 𝑬 is the transformation matrix transforming SD-ambiguities into DD-ambiguities. 𝑬 has the same form as

𝑫 in (2.23), but the elements, which equal wavelengths, are replaced by value ones.

Integer Ambiguity estimation

The elements of 𝒃

𝐷𝐷 in (2.27) are float values, but the DD-ambiguities are intrinsically integer values. If the

correct integer ambiguities are successfully resolved, they can be used to constrain the float solution 𝒃

𝐷𝐷 so that

the accuracy of coordinates can be improved to sub-centimetre level with fewer observations (Blewitt 1989,

Dong and Bock 1989, Ge et al. 2008). The aim of the ambiguity resolution can be expressed by

𝒃

=𝐹(𝒃

), (2.28)

where function 𝐹( ) maps the ambiguities from a set of floats to integers. Many methods have been proposed,

such as the ambiguity determination methods with code measurements or code measurements with smoothing

technique (Cocard and Geiger 1992), AFM (Counselman and Gourevitch 1981, Han and Rizos 1996), integer

bootstrapping, LAMBDA method (Teunissen 1995, Chang et al. 2005).

The LAMBDA method has been widely used and will be employed in this study as it efficiently mechanizes the

ILS procedure which maximizes the probability of correct integer estimation (Verhagen et al. 2012). Assuming

that the estimated float DD-ambiguity vector is 𝒃

 with associated VC matrix 𝑸𝒃

, the LAMBDA method is to

solve the ILS problem described by

min(𝒃

−𝒃

)𝑇𝑸𝒃

𝒃

−1(𝒃

−𝒃

), with 𝒃

∈𝒁𝑛, (2.29)

where 𝒃

 is the vector of integer ambiguity candidates. The objective of ILS is to determine the solution 𝒃

∈𝒃

 for

the minimization problem (2.29). The LAMBDA procedure contains mainly two steps, the reduction step and the

search step (Chang et al. 2005). The former step is to decorrelate the elements in 𝒃

 and order the diagonal entries

by so called Z-transformations, which can shrink the search space. The later step is to find the optimal ambiguity

candidates in a hyper-ellipsoidal space by searching. More details can be found in (Teunissen 1995, Chang et al.

2005).

Integer Ambiguity Validation

The obtained ambiguity combination 𝒃

 is supposed to be the best one. However, vector 𝒃

 could be a wrong

integer ambiguity vector. Thus, validation methods have been proposed to check its reliability. They can be

classified into two classes, the approaches of performance evaluation function, which try to evaluate the

performance of the integer ambiguity parameters using the probabilistic properties of the integer ambiguity

estimators, and the discrimination function approaches that try to ensure that the best ambiguity combination,

which is the optimal solution of (2.29), is statistically better than the second best one (Kim and Langley 2000).

One of the discrimination tests, the R-ratio test (Euler and Schaffrin 1990) is employed in this thesis with the

RATIO value calculated by

𝑅𝐴𝑇𝐼𝑂= 𝝈′2

𝝈2 = ‖𝒃

−𝒃

′‖2𝑸𝒃

𝒃



‖𝒃

−𝒃

‖2𝑸𝒃

𝒃

, (2.30)

where 𝒃

′ is the second best ambiguity vector according to (2.29). If RATIO value equals to or is larger than a

threshold, the integer ambiguity vector 𝒃

 will be considered true, while 𝒃

 will be refused if RATIO value is

smaller than the threshold. The threshold is usually set to a fixed value between 1.5 and 3. In this thesis, when a

threshold is needed, the value is set to 3 as reported in (Leick 2004).

Fixed Baseline Solution

If the integer ambiguity vector passes the validation test, it will be considered as the true ambiguity vector. Then

𝒃

 will be used to adjust the float solution of other parameters leading to the corresponding fixed solution. This

process is based on the correlation between these parameters and the ambiguity parameters (Teunissen 2002),

which can be expressed by

𝒙

=𝒙

−𝑸𝒙

𝒃

𝑸𝒃

𝒃

−1(𝒃

−𝒃

), (2.31a)

𝑸𝒙

𝒙

=𝑸𝒙

𝒙

−𝑸𝒙

𝒃

𝑸𝒃

𝒃

−1𝑸𝒃

𝒙

, (2.31b)

where 𝒙 is the sate vector except for ambiguity variables; 𝒙

 is the fixed solution of 𝒙; 𝑸𝒙

𝒙

 is the VC matrix of

the fixed solution 𝒙

; 𝑸𝒃

𝒙

 denotes the VC matrix of 𝒃

 and 𝒙

; 𝒃

 is the float ambiguity solution; 𝒙

 is the float

solution of 𝒙.

The fixed solution 𝒙

 is usually very accurate, reaching sub-centimetre level. If errors in the observation models

are removed, the successful ambiguity fixing will need observations of only fewer epochs, even a single epoch,

leading to very short convergence time. Therefore, the fast and reliable integer ambiguity resolution is the basis

of the real-time and high-precision GNSS positioning.

3 Particle Filter

Particle filter is a class of filters which recursively approximate variables in the filtering in a simulation way by

weighted samples. As a result, a set of samples with proper weights, which can approximately represent the

density of the variables, is achieved. The sample points are also named particles. The theoretical details can be

found in (Douct et al. 2001, Chen 2003, Dimov and McKee 2008, Candy 2009, Gustafsson 2010, Haug 2012).

This chapter provides a review of the basic estimation theory described in the above literature. The Bayesian

filtering is firstly presented (Douct et al. 2001, Dimov and McKee 2008), and then the Kalman filter is

introduced in the frame of Bayesian filtering (Chen 2003, Candy 2009). Afterwards, the sequential importance

sampling (SIS) method, resampling method, bootstrap particle filter as well as RPF are introduced (Candy 2009,

Gustafsson 2010, Haug 2012).

3.1 Discrete-Time State-Space Model

In a discrete-time state-space system, the state vector x can be expressed by its PDF 𝑝(𝒙). In a filtering process,

the state vector x is transited from one epoch to the next and then is updated by measurements. Thus, 𝑝(𝒙)

depends on the estimate of the former epoch and the measurements at the current epoch with prediction model

and the measurement model, respectively. For each epoch k, the two models can be expressed by

𝒙𝑘=𝑓𝑘(𝒙𝑘−1,𝝐𝑘), (3.1a)

𝒚𝑘=ℎ𝑘(𝒙𝑘,𝒆𝑘), (3.1b)

where 𝒚𝑘 is the measurement vector at epoch k; 𝑓𝑘( ) is the prediction function; ℎ𝑘( ) is the measurement

function; 𝝐𝑘 and 𝒆𝑘 are the state noise and the measurement noise, respectively. Model (3.1a) describes a

first-order Markov process, which indicates the estimated state vector is only related to the solution of the

previous epoch 𝒙𝑘−1, but not to other solutions 𝒙𝑘−2, …, 𝒙0.

3.2 Bayesian filtering

The posterior density of 𝒙𝑘 estimated at epoch k depends on observations 𝒚1:𝑘={𝑦1,𝑦2,…,𝑦𝑛} and is denoted

as 𝑝(𝒙𝑘|𝒚1:𝑘). The prior density 𝑝(𝒙𝑘|𝒚1:𝑘−1) is needed to derive 𝑝(𝒙𝑘|𝒚1:𝑘) and can be calculated by the

prediction model from the previous posterior density 𝑝(𝒙𝑘−1|𝒚1:𝑘−1) at epoch k-1. With the first-order Markov

process assumption, this process is realised by the Chapman-Kolmogorov equation (Arulampalam et al. 2002)

𝑝(𝒙𝑘|𝒚1:𝑘−1)=∫𝑝(𝒙𝑘|𝒙𝑘−1) 𝑝(𝒙𝑘−1|𝒚1:𝑘−1)𝑑𝒙𝑘−1. (3.2)

Then the posterior density 𝑝(𝒙𝑘|𝒚1:𝑘) can be estimated based on the Bayes’s theorem (Candy 2009, Haug 2012)

and is expressed by

𝑝(𝒙𝑘|𝒚1:𝑘)= 𝑝(𝒚1:𝑘|𝒙𝑘)𝑝(𝒙𝑘)

𝑝(𝒚1:𝑘) = 𝑝(𝒚𝑘|𝒚1:𝑘−1, 𝒙𝑘)𝑝(𝒚1:𝑘−1|𝒙𝑘)𝑝(𝒙𝑘)

𝑝(𝒚𝑘|𝒚1:𝑘−1)𝑝(𝒚1:𝑘−1)

= 𝑝(𝑦𝑘|𝑦1:𝑘−1, 𝑥𝑘)𝑝(𝑥𝑘|𝑦1:𝑘−1)𝑝(𝑦1:𝑘−1)𝑝(𝑥𝑘)

𝑝(𝑦𝑘|𝑦1:𝑘−1)𝑝(𝑦1:𝑘−1)𝑝(𝑥𝑘)

= 𝑝(𝑦𝑘|𝑦1:𝑘−1, 𝑥𝑘)𝑝(𝑥𝑘|𝑦1:𝑘−1)

𝑝(𝑦𝑘|𝑦1:𝑘−1)

= 𝑝(𝑦𝑘|𝑥𝑘)𝑝(𝑥𝑘|𝑦1:𝑘−1)

𝑝(𝑦𝑘|𝑦1:𝑘−1), (3.3)

because the measurements at epoch 𝑘 are not related to the measurements before.

The posterior density of the state vector can be calculated recursively by (3.2) and (3.3). The solution (3.3) is a

generally conceptual solution and there are no analytical forms in most cases (Arulampalam et al. 2002). For the

linear Gauss-Markov model, the optimal analytical expression can be obtained as Kalman filter. Generally, only

suboptimal solutions with approximation are available, such as the extended Kalman filter (EKF), unscented

Kalman filter (UKF), as well as particle filter.

3.3 Kalman Filter

Kalman filter can be derived from Bayesian filtering with linear Gauss-Markov assumption where both the

prediction function and the measurement function are linear, and the variables and measurements have white

noises, and first-order Markov process applies (Candy 2009, Haug 2012). In this case, the prediction model and

update model can be expressed by

𝒙𝑘=𝑭𝒙𝑘−1+𝒗𝑘, (3.4a)

𝒚𝑘=𝑯𝒙𝑘+𝒘𝑘. (3.4b)

The predicted value of the state vector at epoch k can be calculated by the density-weighted integral

(Haug, 2012)

𝒙

𝑘|𝑘−1=∫𝑭𝒙𝑘−1𝒩(𝒙𝑘−1;𝒙

𝑘−1|𝑘−1,𝜮𝑛−1|𝑛−1)𝑑𝒙𝑘−1=𝑭𝒙

𝑘−1, (3.5a)

𝜮

𝑘|𝑘−1

𝒙𝒙 =∫𝑭(𝒙𝑘−1−𝒙

𝑘−1|𝑘−1)(𝒙𝑘−1−𝒙

𝑘−1|𝑘−1)𝑇𝑭𝑻𝒩(𝒙𝑘−1;𝒙

𝑘−1|𝑘−1,𝜮𝑛−1|𝑛−1)𝑑𝒙𝑘−1+𝑸

=𝑭{∫(𝒙𝑘−1−𝒙

𝑘−1|𝑘−1)(𝒙𝑘−1−𝒙

𝑘−1|𝑘−1)𝑇𝒩(𝒙𝑘−1;𝒙

𝑘−1|𝑘−1,𝜮𝑛−1|𝑛−1)𝑑𝒙𝑘−1}𝑭𝑻+𝑸

=𝑭𝜮

𝑘−1|𝑘−1

𝒙𝒙 𝑭𝑻+𝑸, (3.5b)

where Q is the state noise covariance matrix; 𝒩(𝒙;𝒙

,𝜮) denotes the normal distribution with the expectation 𝒙



and the corresponding covariance matrix 𝜮.

The observation update can be obtained according to (3.3). The detailed description can be found in (Chen 2003,

Candy 2009), hence only the results are presented here.

With the denotations of

𝒚

𝑘|𝑘−1=𝑯𝒙

𝑘|𝑘−1, (3.6a)

𝜮

𝑘|𝑘−1

𝒚𝒚 =𝑯𝜮

𝑘|𝑘−1

𝒙𝒙 𝑯𝑇+𝑹, (3.6b)

𝜮

𝑘|𝑘−1

𝒙𝒚 =𝜮

𝑘|𝑘−1

𝒙𝒙 𝑯𝑻, (3.6c)

where R is the covariance matrix of the measurement noise, the final solution can be expressed by

𝒙

𝑘|𝑘=𝒙

𝑘|𝑘−1+𝑲𝑘(𝒚𝑘−1−𝒚

𝑘|𝑘−1), (3.7a)

𝜮

𝑘|𝑘

𝒙𝒙 =𝜮

𝑘|𝑘−1

𝒙𝒙 −𝑲𝒌𝜮

𝑘|𝑘−1

𝒚𝒚 𝑲𝑘, (3.7b)

where

𝑲𝑘=𝜮

𝑘|𝑘−1

𝒙𝒚 (𝜮

𝑘|𝑘−1

𝒚𝒚 )−1 (3.8)

is the Kalman gain.

KF is the optimal solution of Bayesian filtering under the condition of linear models with the state vector

applying to Gaussian distribution. For the cases of nonlinear problems, two suboptimal solutions EKF and UKF

which are developed based on KF can be employed. In EKF, 𝒙

𝑘|𝑘−1 is calculated by a non-linear prediction

function, 𝜮

𝑘|𝑘−1

𝒙𝒙 is obtained with linearized transition matrix, and 𝑯 is also obtained by linearization. So EKF

expresses the non-linear problem by its fist-order linearization, which can introduce large errors in the estimated

posterior density. Therefore, EKF diverges easily especially to highly non-linear problems. The UKF introduces

the idea of Particle filter and the state distribution is represented by using a minimal set of carefully chosen

sample points (Julier and Uhlmann 1997, Wan and Der Merwe 2000). This method can approximate the

non-linear problem to the third-order of Taylor series expansion and hence has a better performance for highly

non-linear problems than EKF (Wan and Der Merwe 2000).

For the problems where the models have no specific forms or the noise distributions of variables are not

Gaussian, the Kalman filter, as well as the EKF and UKF, is not suitable. Some other methods which can solve

such problems are needed, such as the approaches via the Monte Carlo method with simulation methodology.

One of these approaches, the SIS method will be introduced in the next section 3.4.

3.4 Sequential Importance Sampling

If a PDF 𝑝(𝒙) is represented by N independently and identically distributed samples {𝑥𝑖}𝑖=1

𝑁, then we have

𝑝(𝒙)≈ 1

𝑁 ∑𝛿(𝒙−𝒙𝑖)

𝑁

𝑖=1 , (3.9)

where 𝛿( ) is the Dirac delta function. When the 𝑝(𝒙) of the state vector x is known, the expectation of x can

hence be expressed by

𝒙

=∫𝒙𝑝(𝒙)𝑑𝒙≈ 1

𝑁 ∑𝒙𝑖

𝑁

𝑖=1 . (3.10)

The posterior density 𝑝(𝒙) is usually unknown at the beginning in practice. If there is an available initial PDF

𝑞(𝒙) which may be not as accurate as 𝑝(𝒙), but the particles {𝒙𝑖}𝑖=1

𝑁can be sampled from. After the sampling, all

the particles have the same weight and are distributed according to 𝑞(𝒙). Based on these particles, 𝑝(𝒙) is

estimated when the new information in the form of measurements is available. Considering that the values of the

particles are not distributed according to 𝑝(𝒙) but 𝑞(𝒙), the expectation of the unknown state vector can be

expressed by (Haug, 2012)

𝒙

=∫𝒙𝑝(𝒙)𝑑𝒙=∫𝒙 𝑝(𝒙)

𝑞(𝒙) 𝑞(𝒙)𝑑𝒙=∫𝒙 𝑝(𝒙)

𝑞(𝒙) 𝑞(𝒙)𝑑𝒙

=∫𝒙 𝑤(𝒙) 𝑞(𝒙)𝑑𝒙, (3.11)

where (𝒙)= 𝑝(𝒙)

𝑞(𝒙); 𝑞(𝒙) is called importance density function.

With samples generated according to 𝑞(𝒙), the expectation (3.11) can be written as

𝒙

≈ 1

𝑁 ∑𝑤𝑖𝒙𝑖

𝑁

𝑖=1 , (3.12)

where 𝑤𝑖= 𝑝(𝒙𝑖)

𝑞(𝒙𝑖). Eq. (3.12) is known as the importance sampling estimate of state vector 𝒙. The values 𝑤𝑖 are

named importance weights.

In a discrete state-space system, assuming that the new information is the measurements 𝒚1:𝑘 at each epoch, the

posterior density is denoted by 𝑝(𝒙𝑘|𝒚1:𝑘) and the importance function is denoted by 𝑞(𝒙𝑘|𝒚1:𝑘). According to

(3.11), the expectation of x at epoch k is

𝒙

𝑘=∫𝒙𝑘𝑝(𝒙𝑘|𝒚1:𝑘)𝑑𝒙𝑘=∫𝒙𝑘 𝑝(𝒙𝑘|𝒚1:𝑘)

𝑞(𝒙𝑘|𝒚1:𝑘) 𝑞(𝒙𝑘|𝒚1:𝑘)𝑑𝒙. (3.13)

Then the importance weight function is

𝑤(𝒙𝑘)= 𝑝(𝒙𝑘|𝒚1:𝑘)

𝑞(𝒙𝑘|𝒚1:𝑘). (3.14)

According to the Bayesian estimation (3.3), where the denominator 𝑝(𝒚𝑘|𝒚1:𝑘−1) is a normalization term and

can be removed, equation (3.14) is changed into

𝑤(𝒙𝑘)∝ 𝑝(𝒚𝑘|𝒙𝑘)𝑝(𝒙𝑘|𝒚1:𝑘−1)

𝑞(𝒙𝑘|𝒚1:𝑘) = 𝑝(𝒚𝑘|𝒙𝑘)∫𝑝(𝒙𝑘|𝒙𝑘−1)𝑝(𝒙𝑘−1|𝒚1:𝑘−1)𝑑𝒙𝑘−1

∫𝑞(𝒙𝑘|𝒙𝑘−1,𝒚1:𝑘)𝑞(𝒙𝑘−1|𝒚1:𝑘−1)𝑑𝒙𝑘−1 . (3.15)

Assume samples {𝒙𝑘−1

𝑖}𝑖=1

𝑁 are generated from 𝑞(𝒙𝑘−1|𝒚1:𝑘−1) at epoch k-1, the right side of eq. (3.15) can be

written in the discrete form as (Candy 2009, Haug 2012).

𝑤(𝒙𝑘)=∑

𝑁

𝑖=1 𝑝(𝒚𝑘|𝒙𝑘)𝑝(𝒙𝑘|𝒙𝑘−1

𝑖)𝑝(𝒙𝑘−1

𝑖|𝒚1:𝑘−1)

𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚1:𝑘)𝑞(𝒙𝑘−1

𝑖|𝒚1:𝑘)

=∑

𝑁

𝑖=1 𝑝(𝒚𝑘|𝒙𝑘)𝑝(𝒙𝑘|𝒙𝑘−1

𝑖)𝑝(𝒙𝑘−1

𝑖|𝒚1:𝑘−1)

𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚1:𝑘)𝑞(𝒙𝑘−1

𝑖|𝒚1:𝑘)

=∑

𝑁

𝑖=1 𝑝(𝒚𝑘|𝒙𝑘)𝑝(𝒙𝑘|𝒙𝑘−1

𝑖)𝑝(𝒙𝑘−1

𝑖|𝒚1:𝑘−1)

𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚𝑘)𝑞(𝒙𝑘−1

𝑖|𝒚1:𝑘−1)

=∑

𝑁

𝑖=1 𝑤𝑘−1

𝑖 𝑝(𝒚𝑘|𝒙𝑘)𝑝(𝒙𝑘|𝒙𝑘−1

𝑖)

𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚𝑘), (3.16)

where 𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚1:𝑘)= 𝑞(𝒙𝑘|𝒙𝑘−1

𝑖,𝒚𝑘) because 𝒙𝑘 does not depend on 𝒚1:𝑘−1 ;

𝑞(𝒙𝑘−1

𝑖|𝒚1:𝑘)= 𝑞(𝒙𝑘−1

𝑖|𝒚1:𝑘−1) because 𝒙𝑘−1 is not related to 𝒚𝑘. Generate {𝒙𝑘

𝑖}𝑖=1

𝑁 from 𝑞(𝒙𝑘|𝒚1:𝑘), then

𝑞(𝒙𝑘|𝒚1:𝑘)= 1

𝑁 ∑𝛿(𝒙−𝒙𝑖)

𝑁

𝑖=1 . (3.17)

According to (3.14) and (3.17), the posterior density function can be expressed as

𝑝(𝒙𝑘|𝒚1:𝑘)=𝑤(𝒙𝑘)𝑞(𝒙𝑘|𝒚1:𝑘)

∝𝑤(𝒙𝑘)𝑞(𝒙𝑘|𝒚1:𝑘)

≈∑

𝑁

𝑖=1 1

𝑁 𝑤𝑘𝑖𝛿(𝒙−𝒙𝑖). (3.18)

Normalising the right side of (3.18) by

𝑤𝑘𝑖= 𝑤

𝑘

𝑖𝑁

⁄

∑𝑤

𝑘

𝑖𝑁

⁄

𝑁

𝑖=1 , (3.19)

the posterior density is

𝑝(𝒙𝑘|𝒚1:𝑘)≈∑𝑤𝑘𝑖𝑁

𝑖=1 𝛿(𝒙−𝒙𝑖). (3.20)

Therefore, the estimate of 𝒙𝑘 is

𝒙

𝒌=∫𝒙𝑘𝑝(𝒙𝑘|𝒚1:𝑘)𝑑𝒙𝑘≈∑𝒙𝑘

𝑖𝑁

𝑖=1 𝑤𝑘𝑖. (3.21)

Equation (3.16) includes the weights of the last epoch and thus, the weights in SIS method can be updated

recursively.

In SIS algorithm, the PDF of state vector 𝒙 is represented by samples with sample number N. Even though a

larger N results in more accurate PDF, number N is usually set to a moderate number so that the computation

burden is acceptable.

The particle weights in SIS are updated every epoch. As the SIS proceeds on, the weights of most samples

decrease and become very small, while only the remaining few samples have very large weights. In this case,

only the fewer samples with large weights affect the estimation and thus the PDF cannot be well represented.

This problem is referred to as the degeneracy, which can be solved by resampling algorithm.

3.5 Resampling

The resampling algorithm eliminates the weight differences of the particles by deleting the samples with small

weights and duplicating the samples with large weights. The resampling is usually realized based on the

cumulative distribution function (CDF) of the PDF, which has values from 0 to 1. In the resampling step, the

interval of CDF’s value [0, 1] is sampled randomly instead of sampling x directly so that the corresponding x

samples in CDF will have the same weights. These x samples will be the new particles.

To numerically demonstrate the process of resampling, a variable x with errors apply to the normal distribution is

taken as an example. The Gaussian PDF can be expressed by

𝑝(𝑥)=𝒩(𝑥;𝑥,𝜎2)= 1

𝜎√2𝜋𝑒−1

2𝜎2(𝑥−𝑥)2, (3.22)

where 𝑥 is the mean value of x and 𝜎 is the STD of the distribution.

Assume 𝑥 belongs to distribution 𝒩(0,1), the probability for a random 𝑥 value being within three STDs of the

mean is 99.73% and hence the initial interval is selected as [-3, 3]. This interval is evenly sampled with the

number of samples 200. The probabilities of these samples are calculated with (3.22) and shown in the left panel

of Fig. 3.1, which are then normalised so that the sum of these probabilities equals value 1. After that, the

accumulated sum of these normalised values i.e. the simulated CDF which is denoted by 𝑃(𝒙), is calculated and

depicted in the right panel of Fig. 3.1.

To generate samples of the same weight, the CDF 𝑃(𝑥) which is between 0 and 1 is randomly sampled so that

the new samples will have equal probabilities. These new samples correspond to different x values in the right

panel of Fig. 3.1. These new x values compose the new sample collection. In this process, the new samples can

only be chosen from the original samples i.e. from the 200 samples in this example, but the number of the new

samples can be different with that of the original samples by just sampling over [0, 1] with sampling number of a

new value.

Fig. 3.1 PDF of the normal distribution 𝒩(𝑥;0,1) within [-3, 3] (left), and the corresponding CDF (right)

According to the above approach, several resampling algorithms have been developed until now, including the

methods of multinomial resampling (Gordon et al. 1993), stratified resampling (Kitagawa 1996, Doucet et al.

2001), systematic resampling (Arulampalam et al. 2002), residual resampling (Liu and Chen 1998) and so on.

The stratified resampling method and the systematic resampling method are preferred compared with the others

and these two methods have similar procedures and performances (Hol et al. 2006 ). In the following sections,

the stratified resampling method is employed.

The procedure of the stratified resampling method is as follows.

Step 1: Construct the numerical CDF {𝑊𝑘𝑖}𝑖=1

𝑁 of x

For each i=1, …, N, 𝑊𝑘𝑖=∑𝑤𝑘𝑗

𝑖𝑗=1 .

Step 2: Generate random CDF {𝑢𝑖}𝑖=1

𝑁

For each i=1, …, N, 𝑢𝑖= (𝑖−1)+𝑢



𝑁, where 𝑢 is a random value over interval [0, 1].

Step 3: Generate the new collection {𝒙

𝑘

ℎ}ℎ=1

𝑁 of x by comparing the elements of {𝑢𝑖}𝑖=1

𝑁 with that of {𝑊𝑘𝑖}𝑖=1

𝑁

For i=1, …, N, set m=1,

For each i, compare 𝑤𝑘𝑚 with 𝑢𝑖; if 𝑤𝑘𝑚<𝑢𝑖, delete 𝒙𝑘

𝑚 by setting m=m+1, else duplicate 𝒙𝑘

𝑚 by setting

𝒙

𝑘

𝑖=𝒙𝑘

𝑚.

Set the corresponding weight for each 𝒙

𝑘

𝑖 with 𝑤𝑘𝑖= 1

𝑁.

Even though the resampling method solves the degeneracy problem of SIS, the resampling may lead to the loss

of diversity. In extreme case, there are N samples with the same values and the PDF of the variable still cannot

be well represented. This phenomenon is called sample impoverishment, which can be overcome by adding

noise or by using some regularization method. Although the former one can usually solves the problem, if the

prediction model with a very low noise level or even free of noise is employed, the regularization method should

be used. This will be demonstrated in section 4.5.

3.6 Bootstrap Filter

After introducing the resampling step to the SIS, the method is so called sequential importance resampling (SIR)

(Doucet et al. 2001). The illustration of SIR is shown in Fig. 3.2.

Fig. 3.2 Illustration of the SIR method. The positions of the dots represent the samples while the sizes represent

the weights

The optimal importance function in (3.16) is impossible to be determined analytically in many cases. So the

suboptimal importance function is usually selected. One of the methods is to let the importance function equal

the transition density (Doucet et al. 2000), which can be expressed by

𝑞(𝒙𝑘|𝒙𝑘−1,𝒚𝑘)=𝑝(𝒙𝑘|𝒙𝑘−1). (3.23)

Update

{𝒙𝑘

𝑖,𝑤𝑘𝑖}𝑖=1

𝑁

×2

×3

Resample

{𝒙

𝑘

𝑖,𝑤0𝑖}𝑖=1

𝑁

Predict

{𝒙k+1

𝑖,𝑤0𝑖}𝑖=1

𝑁

Update

{𝒙𝑘+1

𝑖,𝑤𝑘+1

𝑖}𝑖=1

𝑁

×2

×3

×2

Resample

{𝒙

𝑘+1

𝑖,𝑤0𝑖}𝑖=1

𝑁

Predict

{𝒙𝑘+2

𝑖,𝑤0𝑖}𝑖=1

𝑁

{𝒙0

𝑖,𝑤0𝑖}𝑖=1

𝑁

Initialize

Then the weight update is simplified as

𝑤𝑘𝑖=𝑤𝑘−1

𝑖𝑝(𝒚𝑘|𝒙𝑘

𝑖). (3.24)

This kind of particle filter is named bootstrap particle filter. The procedure is as follows (Candy 2009,

Gustafsson 2010, Haug 2012),

Step 1: Initialize filter

Generate samples {𝒙0

𝑖}𝑖=1

𝑁, with 𝒙0

𝑖~𝑞(𝒙0).

Assign the weights {𝑤0𝑖}𝑖=1

𝑁.

Step 2: Sequential importance sampling

Draw new samples {𝒙𝑘

𝑖}𝑖=1

𝑁, by

𝒙𝑘

𝑖=𝑓(𝒙𝑘−1

𝑖)+𝒗𝑘. (3.25)

Update the weights according to likelihood function 𝑝(𝒚𝑘|𝒙𝑘

𝑖) of measurements with

𝑤𝑘𝑖=𝑤𝑘−1

𝑖𝑝(𝒚𝑘|𝒙𝑘

𝑖). (3.26)

Normalize the weights by

𝑤𝑘𝑖= 𝑤

𝑘

𝑖

∑𝑤

𝑘

𝑗

𝑁

𝑗=1 . (3.27)

Calculate the estimated value and variance by

𝒙

𝑘≈∑𝒙𝑘

𝑖𝑤𝑘𝑖𝑁

𝑖=1 , (3.28)

var(𝒙

𝑘)≈∑(𝒙𝑘

𝑖−𝒙

𝑘)(𝒙𝑘

𝑖−𝒙

𝑘)𝑇𝑤𝑘𝑖𝑁

𝑖=1 . (3.29)

Step 3: Resample if

𝑁𝑒𝑓𝑓<𝑁𝑡ℎ, (3.30)

where 𝑁𝑒𝑓𝑓 is the effective number of samples which is calculated by

𝑁𝑒𝑓𝑓= 1

∑(𝑤𝑘

𝑖)2

𝑁

𝑖=0 (3.31)

and 𝑁𝑡ℎ is a threshold which can be set to the value of 2

3 𝑁.

Step 4: Repeat steps 2 and 3 for the following epochs.

The procedure above is utilised in the remaining part of this study as normal particle filter procedure.

3.7 Regularized Particle Filter

The resampling step in particle filter is to deal with the degeneracy problem. However, it deletes the particles

with small weights leading to the loss of diversity. In the worst case, all the particles will have the same value

and the PDF cannot be well represented. Thus, the regularization method is introduced to increase the diversity

of the particles (Doucet et al. 2001).

In the RPF, each particle is jittered by a small value which is calculated according to a local individual kernel.

The kernel is constructed to minimize the distance between the true posterior density and its regularized

empirical representation. After the regularization, the particles which have the same value move a little away

from each other and the diversity of the particles increases. The regularization procedure can be implemented

after the resampling step; it can also be located before the weight update step (Doucet et al. 2001). The

regularization does not guarantee the asymptotical approximation of the posterior PDF by the particles and is

only necessary when the diversity loss is severe (Arulampalam et al. 2002).

The PDF 𝑝(𝑥𝑘) of a variable x can be approximated as the mixture of many individual PDFs, whose kernel is

supposed to be symmetric and has the following characteristics

𝐾(𝑥)≥0,∫𝐾(𝑥)𝑑𝑥=1,∫𝑥𝐾(𝑥)𝑑𝑥=1,∫‖𝑥‖2𝐾(𝑥)𝑑𝑥<∞. (3.32)

Then the PDF 𝑝(𝑥𝑘) can be expressed by

𝑝(𝑥𝑘)≈∑𝑤𝑘𝑖𝐾ℎ(𝑥𝑘−𝑥𝑘𝑖)

𝑁

𝑖=1 =𝑝(𝑥𝑘), (3.33)

where 𝐾ℎ(𝑥)= 1

ℎ𝑛𝑥 𝐾( 𝑥

ℎ ) is the rescaled kernel; ℎ>0 is the bandwidth; 𝑛𝑥 is the dimension of the state

vector.

The difference between the true PDF 𝑝(𝑥𝑘) and the probability density 𝑝(𝑥𝑘) in (3.33) is expressed by the mean

integrated square error (MISE)

𝑀𝐼𝑆𝐸(𝑝)=𝐸[∫‖𝑝(𝑥𝑘)−𝑝(𝑥𝑘)‖2

2𝑑𝑥𝑘]. (3.34)

The choice of the kernel 𝐾 and the bandwidth ℎ should be done in such a way that the MISE is minimized

(Arulampalam et al. 2002). In the case of samples with equal weights, the optimal choice for the kernel is

𝐾𝑜𝑝𝑡(𝑥)={ 𝑛𝑥+2

2𝑐𝑛𝑥 (1−‖𝑥‖2) 𝑖𝑓 ‖𝑥‖<1

0 otherwise, (3.35)

where 𝑐𝑛𝑥 is the volume of the unite sphere of 𝑅𝑛𝑥. This kernel is called Epanechnikov kernel. Assuming the

underlying true density is Gaussian with a unit covariance matrix, the optimal bandwidth is

ℎ𝑜𝑝𝑡=[8𝑐𝑛𝑥−1(𝑛𝑥+4)(2√𝜋)𝑛𝑥]1

𝑛𝑥+4𝑁− 1

𝑛𝑥+4. (3.36)

In a general case of an arbitrary underlying density, the underlying density is assumed to be Gaussian with

variance 𝑆 which equal to empirical variance of the samples. Then the kernel function (3.35) is changed into the

following rescaled kernel

(det𝐴𝑘)−1

ℎ𝑛𝑥 𝐾( 1

ℎ 𝐴𝑘−1𝑥), (3.37)

with 𝐴𝑘𝐴𝑘−1=𝑆𝑘.

The equation (3.36) for h can still be used directly. The new samples after regularization are

𝑥𝑘𝑖∗=𝑥𝑘𝑖+ℎ𝑜𝑝𝑡𝐴𝑘𝜖𝑖, (3.38)

where 𝜖𝑖 is generated from the kernel by sampling.

3.8 Other Particle Filter Methods

Gaussian Particle Filter

Gaussian particle filter assumes that the posterior distribution is Gaussian distribution. Therefore after the

expectation and variance are calculated in the measurement update step, the analytical expression of the posterior

distribution is known. Thus, the resampling algorithm based on the posterior density in the SIR procedure is

replaced by an algorithm sampling from the Gaussian distribution (Kotecha and Djuric 2003).

Compared with other kinds of Gaussian filter, such as EKF and UKF, this method has a much better

performance in the highly non-linear dynamic models, and compared with the normal particle filter, Gaussian

particle filter has a lower complexity due to the absence of resampling step (Kotecha and Djuric 2003). In the

case of non-Gaussian posterior distribution, it can be approximated by weighted Gaussian mixtures and the

procedure is named Gaussian sum particle filter (Kotecha and Djuric 2003).

Auxiliary Particle Filter

The transition density 𝑝(𝒙𝑘|𝒙𝑘−1) in (3.23) is usually 𝑝(𝒙𝑘|𝒙𝑘−1,𝒚1:𝑘−1), which means that the density is under

the condition of observations 𝒚1:𝑘−1. However, when the observations 𝒚𝑘 at epoch k is known, the optimal

transition density should be 𝑝(𝒙𝑘

𝑖|𝒙𝑘−1

𝑖,𝒚1:𝑘). The absence of observation 𝒚𝑘 may lead to the relatively

inaccurate 𝑝(𝒙𝑘

𝑖|𝒙𝑘−1

𝑖,𝒚1:𝑘−1). Therefore, some particles may not be at critical positions and hence have very

small weights in the update step resulting in low efficiency, which can be improved by incorporating the

observations 𝒚𝑘.

With this idea, Pitt and Shephard (1999) proposed the method of auxiliary particle filter (APF). The APF

includes two stages. The first stage is to resample the particles at epoch k-1 by using the observations of the

current epoch and the second stage is to reweight the particles at epoch k by the likelihood weights (Doucet et al.

2001). Even though this method may improve the efficiency of particles, it is actually not proper to be used in

this thesis. To incorporate observations 𝒚𝑘 at epoch k-1, the likelihood value for each particle is calculated for an

additional time. In the study, this procedure almost doubles the time consumption which has already been critical

for the estimation approach. Moreover, both IFB rate and F-ISB are stable most of the time, which indicates the

prediction model is pretty accurate and the improvement by 𝑝(𝒙𝑘

𝑖|𝒙𝑘−1

𝑖,𝒚1:𝑘) should be insignificant.

Distributed Particle Filter

If the measurements are not convenient to be processed in a centralized way, the distributed method is needed.

This method divides the measurements into several groups and processes them, separately. Hence the

measurement likelihood is factorized into several local functions. In the measurement update step of particle

filter, the weights of particles are determined locally and then they are combined together according to some

criteria to determine the final solution (Zhao and Nehorai 2007). This method enables the selective collaboration

of sensors to reduce latency and minimize bandwidth consumption in communication (Zhao et al. 2002).

Distributed particle filter may be interesting in the in-door positioning, but it is not proper for GNSS precise

positioning at present as all the measurements at each epoch are processed simultaneously.

Adaptive Particle Filter

In particle filter, the weights of the particles depend mainly on the measurements and therefore the large noise in

the measurements affects the accuracy of the results. The adaptive particle filter proposed by (Zhao 2014) aims

to decrease the effects of the large noise in measurements by lowering the weights of the measurements. In this

adaptive particle filter, a set of virtual observations is simulated using the information from prior distribution.

Then a belief factor is set to tune the weights of the practical measurements. This method employs the

methodology which is similar with the adaptive Kalman filter proposed in (Yang et al. 2001) and enables the

procedure to achieve a robust performance in the high noisy environment.

Particle filter is based on Monte Carlo method and estimates the unknown parameters by samples. After the

samples are generated at the beginning, the following calculation of the estimation is with known parameter

values which enable the employment of stronger models. This is very meaningful to bias estimation for GNSS

ambiguity fixing. In the following chapters, the new method for IFB rate / F-ISB estimation based on particle

filter will be investigated.

4 Online GLONASS Ambiguity Resolution Based on Online Phase IFB

Estimation

This chapter first introduces the existing methods for carrier phase IFB estimation in GLONASS data processing

in section 4.1. Later on, the RATIO distribution with different IFB values is analysed with practical data in

section 4.2. A method based on particle filter is proposed to estimate the IFB rate fast and online in section 4.3.

Afterwards, the procedure based on regularization particle filter to enable the noise tuning in the prediction

model, and the procedure relating the number of particles to their STD to reduce the consumed time by reducing

the number of particles are proposed in section 4.5 and section 4.6, respectively. Finally, the estimated carrier

phase IFBs of some baselines are presented in section 4.7. The RATIO distribution analysis and the IFB rate

estimation approach have been published in (Tian et al. 2015).

4.1 Between-Receiver Phase IFB Characteristics and Existing Methods for

GLONASS IFB Estimation

Although the common parts of the IFB on the satellites and at the receivers are eliminated by the DD-model, the

differential IFB values are usually too large to be neglected when receivers of different types from different

manufacturers are employed. The values are not integer multiples of the wavelength and therefore cannot be

absorbed by the integer ambiguities leading to the failure of integer ambiguity fixing if not properly handled.

However, the IFB can be estimated and removed from the DD-model by IFB modelling. The IFB model depends

on several known IFB characteristics. The between-receiver IFB values have a linear relationship with the

channel number and are similar on GLONASS frequency bands L1 and L2; the receivers from the same

manufacturer have similar IFB values.

Wanninger and Wallstab-Freitag (2007) and Wanninger (2012) estimated the IFB rate with the linear

relationship assumption. In their research, both GPS and GLONASS data are processed together by SD-models.

The carrier pahse SD-models for GPS and GLONASS can be expressed by

𝜆𝛷𝑎𝑏

𝑖,𝐺=𝜌𝑎𝑏

𝑖,𝐺−𝛿𝑡𝑎𝑏

𝐺𝑐+𝜆𝑁𝑎𝑏

𝑖,𝐺+𝜉𝑎𝑏

𝑖,𝐺, (4.1a)

𝜆𝑖𝛷𝑎𝑏

𝑖,𝑅=𝜌𝑎𝑏

𝑖−𝛿𝑡𝑎𝑏

𝑅𝑐+𝑘∆𝛾𝑎𝑏+𝜆𝑖𝑁𝑎𝑏

𝑖+𝜉𝑎𝑏

𝑖,𝑅, (4.1b)

where 𝐺 and 𝑅 refer to GPS and GLONASS, respectively. At the beginning, an initial value of IFB rate is

needed to remove a large part of the biases and then the singularity caused by the ambiguity and clock

parameters is removed by fixing one of the SD-ambiguities to an arbitrary value. After that, the other ambiguities

are estimated without IFB rate estimation. If one of the remaining SD-ambiguities is fixed to its true value, the

IFB rate can be estimated alongside the other remaining ambiguities (Wanninger 2012).

Another approach is presented by (Al-Shaery et al. 2013) where the IFB rates in both carrier phase and code

pseudorange observations in the DD-model are estimated along with the DD-ambiguities. Both kinds of IFBs are

considered to have linear relationships with the channel number. The DD-models in this method can be

expressed by

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐺+𝜀𝑎𝑏

𝑖𝑗,𝐺𝐺, (4.2a)

𝑃𝑎𝑏

𝑖𝑗,𝑅𝑅=𝜌𝑎𝑏

𝑖𝑗,𝑅𝑅+(𝑘𝑗,𝑅−𝑘𝑖,𝑅)∆𝛿+𝜀𝑎𝑏

𝑖𝑗,𝑅𝑅, (4.2b)

𝜆𝛷𝑎𝑏

𝑖𝑗,𝐺𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐺+𝜆𝑁𝑎𝑏

𝑖𝑗,𝐺𝐺+𝜉𝑎𝑏

𝑖𝑗,𝐺𝐺, (4.2c)

𝜆𝑗,𝑅𝛷𝑎𝑏

𝑗,𝑅−𝜆𝑖,𝑅𝛷𝑎𝑏

𝑖,𝑅=𝜌𝑎𝑏

𝑖𝑗,𝑅𝑅+𝜆𝑗,𝑅𝑁𝑎𝑏

𝑗,𝑅−𝜆𝑖,𝑅𝑁𝑎𝑏

𝑖,𝑅+(𝑘𝑗,𝑅−𝑘𝑖,𝑅)∆𝛾𝑎𝑏+𝜉𝑎𝑏

𝑖𝑗,𝑅𝑅, (4.2d)

where ∆𝛿 is the IFB rate for code pseudorange observations. All the unknown parameters in (4.2) are solved at

the same time. If the estimated IFB rates in the float solution are accurate enough, the DD-ambiguities can be

fixed as integers. Afterwards, the ambiguity parameters in the equations are replaced by the fixed ambiguities,

and a more accurate IFB rate can be estimated.

Even though models (4.1) and (4.2) are different, the principles are actually similar. In both approaches, the IFB

rate is estimated together with the ambiguity parameters. Once the estimated ambiguities are considered to be

accurate enough, such as they are fixed as integer ambiguities, the IFB rate is refined. Both the models and the

examples in the research by (Wanninger 2012, Al-Shaery et al. 2013) include the measurements of GPS, which

provides assistance information. In general, almost all current approaches try to estimate simultaneously the

ambiguities and the IFB rate. Unfortunately, the estimation needs a long data set even with the assistance of

simultaneous GPS observations. Consequently, none of these methods can provide a fast or real-time solution of

the IFB rate for GLONASS integer ambiguity resolution without an a priori IFB rate value.

Assuming that the IFB rate is exactly known, then it is easy to understand that for zero or short baselines

employed in the above-mentioned studies, the integer ambiguity resolution could be carried out very reliably, i.e.

with a significantly large RATIO. It is obvious that statistically the closer the IFB rate is to its true value, the

larger the RATIO will be. As IFB rate is usually within an interval of [-0.10, 0.10] in unit of metres per

frequency number (m/FN) (Wanninger 2012), a limited number of samples of IFB rate uniformly distributed can

be defined over the interval. After introducing the samples one by one into the processing for integer ambiguity

resolution, in principle the best estimate of IFB rate can be found out according to the resulting RATIO values.

The rigorous estimation can fortunately be realized via the particle filter which was developed exactly for

providing solution to such kind of estimation problems (Gordon et al. 1993, Doucet et al. 2000, Gustafsson et al.

2002).

Therefore, instead of estimating the IFB rate and ambiguities simultaneously, a new approach is developed in

this chapter to find out the IFB rate estimation which can bring the best performance of integer ambiguity

resolution for observations over all epochs. The estimation is realized by means of particle filtering with

likelihood function of RATIO. Experimental validation shows that this approach can provide a very precise IFB

rate estimation just with GLONASS data of a few epochs, which of course depends on the inter-station distance.

As soon as IFB rate has converged, GLONASS integer ambiguity resolution is usually available and the position

accuracy can be improved significantly. Hence, the new approach can be applied to real-time applications

without any a priori IFB information.

4.2 Relationship between RATIO and Phase IFB Rate

This section aims to investigate the relationship between IFB rate and RATIO by experiments. The correct IFB

rate values should enable the accurate IFB correction and lead to relatively larger RATIO values in GLONASS

data processing.

Although the IFBs for code pseudoranges are not the same with these of carries-phases, the code pseudorange

observations are significantly down-weighted with respect to phases due to their much larger noise. Hereby, the

differences between IFBs for code pseudoranges and carrier phases can be ignored without noticeable bad effects

on the solution. Then the observation model employed here is (Tian et al. 2015)

𝑃𝑎𝑏

𝑖𝑗 =𝜌𝑎𝑏

𝑖𝑗 +(𝑘𝑖−𝑘𝑗)∆𝛾𝑎𝑏+𝜀𝑎𝑏

𝑖𝑗, (4.3a)

𝜆𝑛

𝑖𝜑𝑎𝑏

𝑖−𝜆𝑛

𝑗𝜑𝑎𝑏

𝑗=𝜌𝑎𝑏

𝑖𝑗 +𝜆𝑛

𝑖𝑁𝑎𝑏

𝑖−𝜆𝑛

𝑗𝑁𝑎𝑏

𝑗+(𝑘𝑖−𝑘𝑗)∆𝛾𝑎𝑏+𝜉𝑎𝑏

𝑖𝑗. (4.3b)

The NEQ for models with unknown and known IFB rates are (2.25) and (2.26), respectively.

Data from three baselines are employed in the following numerical analysis. The first baseline is a zero-baseline

using the same type of receivers and antennas, Trimble NetR9 and TRM55971.00, respectively, with a sampling

rate of 5 seconds. The data was collected on day of year (DOY) 182 of 2014, from 9:10:35AM to 12:20:00PM.

The second baseline consists of two IGS stations KOSG and KOS1 in Holland, with a baseline length of about

814 m. The data was collected on DOY 048 of 2014 for 24 hours, with a sampling rate of 30 s, where KOSG is

equipped with LEICA GRX1200GGPRO receiver and AOAD/M_B antenna, while KOS1 is equipped with

SEPT POLARX4 receiver and LEIAR25.R3 antenna. The third one is a kinematic baseline composed of a

reference (REF6) and a rover station (AIR5) and the data of 1 Hz sampling rate were near Munich in Germany

DOY 158 of 2012, from 4:21:05AM to 05:05:27AM, with a maximum inter-station distance of about 1 km. Both

stations were equipped with JAVAD DELTA G3T receivers, but with different antennas ACCG5ANT_42AT1

and LEIAS10 respectively.

As shown in the previous research by Wanninger (2012), it is reasonable to assume that the largest ∆𝛾𝑎𝑏 is less

than 0.10 m/FN, where FN refers to Frequency Number. Therefore, ∆𝛾𝑎𝑏 is assumed to be within [-0.10, 0.10] in

units of m/FN. This interval is evenly sampled with a step-size of 1 mm, so that there are 200 samples. Then

each sample value is used as exactly known ∆𝛾𝑎𝑏 to obtain the corresponding float solution with (2.26).

Afterwards, the LAMBDA method is applied to obtain the RATIO value of the integer ambiguity resolution.

The processing was carried out epoch by epoch for all the data over the three above-mentioned baselines using

L1 and L2 as independent observations. The three-dimensional RATIO maps for all IFB rate samples over all the

epochs are shown in the three sub-plots in Fig. 4.1 for the three baselines, respectively. The RATIO results for

the first epoch corresponding to the three baselines are presented in Fig. 4.2 to show the details of a single epoch.

Theoretically, the highest RATIO value at each epoch should correspond to the true value of IFB rate. For both

the static cases (Fig. 4.1a and Fig. 4.1b) and kinematic data (Fig. 4.1c), there is a clearly detectable peak series in

a straight line with remarkable high RATIO values. The peak series line with almost the same IFB rate for all

epochs also shows its stability. The high RATIO values of the zero-baseline are centralized at the straight line

and much higher than those of another static baseline. For the kinematic baseline, although the RATIO values

are lower than that of the others, they are definitely strong enough to make the positive fixing decision except for

a few outliers. In general, if the IFB rate is given with certain accuracy, the ambiguities can usually be fixed to

integer epoch-by-epoch for both static and kinematic baselines. This can also be seen from the distribution of

RATIOs for the first epoch of the three baselines shown in Fig. 4.2. It is also noticed that there are points very

close to the straight line with rather lower RATIOs, especially for the baseline KOSG-KOS1 and the kinematic

baseline, whereas there are also points far away from the straight line with a rather high RATIO. This is probably

mainly caused by the inaccurate handling of station-specified errors.

Fig. 4.1 Three-dimensional RATIO distribution along with epochs and IFB rate samples for the zero-baseline (a),

KOSG-KOS1 (b) and REF6-AIR5 (c) (The part corresponding to RATIOs which are larger than 50 are not

depicted)

Fig. 4.2 RATIO values of the first epoch corresponding to different IFB rate samples for the zero-baseline (a),

KOSG-KOG1 (b) and REF6-AIR5 (c)

For clarity, Fig. 4.3 shows the distribution of the points with RATIO larger than 3, which can be approximately

considered as the threshold for the acceptance of the corresponding ambiguity candidates. In each sub-plot, there

is a very narrow stripe with a width smaller than ± 4 mm/FN around a line with constant IFB rate. This means

only when the given IFB rate is within the above-mentioned width, the corresponding ambiguities could be fixed

correctly to integer. For some existing methods, where a priori IFB rate value is requested (Wanninger 2012), it

is also a tough job to provide such an accurate value for reliable initial ambiguity resolution.

From the narrow stripe of the zero-baseline, it seems that the IFB rate could be obtained by taking the mean of

the time series. However, there are numerous points scattered beyond and even far away from the stripe. The

boundary of the stripe changes along with the epochs, particularly for the KOSG-KOS1 and the kinematic

baseline. This means that the possibility that the highest RATIO can lead to wrong IFB rates exists. Fig. 4.4

shows the RATIO distribution of three typical epochs: Fig. 4.4a is the worst case where all RATIOs are very low

and the highest RATIO is related to a fully biased IFB rate; in Fig. 4.4b there are two peaks both close to the true

value; and Fig. 4.4c is the usual case with a single peak corresponding to the true value.

Fig. 4.3 Distribution of the points with RATIO larger than the threshold of 3 for the zero-baseline (a), KOSG-

KOS1 (b) and REF6-AIR5 (c)

Fig. 4.4 RATIO plot at the epoch 112 (a), 323 (b) and 800 (c) for baseline REF6-AIR5

Fig. 4.5 Statistics of the relationship of the maximum RATIO and the IFB rate over all epochs for the three

baselines, the zero-baseline (a), KOSG-KOS1 (b), and REF6-AIR5 (c)

In order to obtain a statistical interpretation of their relationship, the maximum RATIO and the corresponding

IFB rate for all epochs are depicted in Fig. 4.5 for the three baselines, where each point in the figures

corresponds to one epoch. In Fig. 4.5a, IFB rates with the maximum RATIO are almost the same, while in

Fig. 4.5b and c, some maximum RATIO values correspond to wrong IFB rate values. Obviously, selecting the

IFB rate value corresponding to the maximum RATIO is not always reliable. Although there might be a number

of methods to estimate the correct IFB rate based on the relationship between RATIO and IFB rate samples

shown in Fig. 4.1, a new method based on particle filtering will be developed to obtain a more reliable solution

in real time.

4.3 Procedure for Phase IFB Rate Online Estimation

The GLONASS observations at a single epoch are processed with the observation model (2.22). The NEQ of

(2.24) is generated at each epoch, but instead of solving the NEQ with (2.25), the NEQ (2.26) is employed with

the IFB rate set to a randomly pre-defined value. The correct IFB rate is then estimated via the methodology of

particle filter described in section 3.6.

The state variable in particle filter is the IFB rate ∆𝛾𝑎𝑏 which is very stable. The prediction model is designed as

∆𝛾𝑎𝑏𝑘

𝑖=∆𝛾𝑎𝑏𝑘−1

𝑖+𝜖∆𝛾, (4.4)

where 𝜖∆𝛾 is assumed to be normally distributed noise.

The key issue here is the likelihood function 𝑝(𝒚𝑘|𝒙𝑘

𝑖) in (3.26) which is usually estimated with observations to

update the particle weights according to (3.26). However, 𝑝(𝒚𝑘|𝒙𝑘

𝑖) does not tell anything about the quality of

the IFB rate in the observation model (2.22), because IFB rate and ambiguities are correlated. In other words, for

any IFB rate we have the same observation residuals in the adjustment. According to the relationship between

the IFB rate and the corresponding RATIO in section 4.2, the RATIO can judge the quality of the pre-defined

IFB rate and be used as the probability. Approximately, the likelihood function for the ambiguities being fixed to

the correct integers under given IFB rate parameter can be expressed by (Tian et al. 2015).

𝑝(𝒃

k|𝑥𝑘𝑖)∝𝑅𝐴𝑇𝐼𝑂𝑖, (4.5)

where ∝ denotes direct proportionality. Because there is only one IFB rate parameter, therefore 𝑥 is not a vector.

As RATIO for all the particles at each epoch are usually of very much different magnitudes, the normalized

RATIO value defined by

𝑝(𝒃

𝑘|𝑥𝑘𝑖)= 𝑅𝐴𝑇𝐼𝑂𝑖

∑𝑅𝐴𝑇𝐼𝑂𝑖

𝑁

𝑖=1 (4.6)

is selected (Tian et al. 2015). It must be pointed out that (4.6) is an empirical expression, although its efficiency

is validated in the following experimental evaluation. Based on the above definition, the particle filter for

estimating the IFB rate can be carried out as follows (Tian et al. 2015):

Step 1: Process the GLONASS observations at current epoch using the observation equations (4.3) and

generate the NEQ in the form of (2.24). Of course, accumulated NEQ over several epochs can also be

used if single epoch ambiguity resolution does not perform well.

Step 2: For the first epoch, an initial set of particles with a certain number of elements must be generated, let’s

say {𝑥0𝑖,𝑤0𝑖}𝑖=1

𝑁. These particles should be uniformly distributed over the interval [-0.10, 0.10] m/FN

and the weights of all particle are 1/N. The total number of the particles N is 200 in this study. For

other epochs k = 2, 3, … the particles are already prepared in the processing of the previous epoch.

Step 3: For each of the particles, a solution of (2.26) is obtained by inserting the IFB rate value of this particle

into the NEQ in step 1. Then the integer ambiguity resolution is undertaken using the LAMBDA

method and the RATIO value is obtained. At the end of this step, we have the RATIOs for all the

particles

Step 4: Update the weight of each particle using (3.26) with the empirical PDF 𝑝(𝒃

𝑘|𝑥𝑘𝑖) of (4.6). Then

normalize the weights and calculate the estimated IFB rate by (3.28) and its STD by (3.29) as well.

Step 5: Resample the particles as described in section 3.5 and transit each particle to next epoch by (4.4), then

the particle set for the next epoch is ready.

Step 6: Repeat the steps 1 to 5 for the epoch 𝑘+1.

This algorithm can be applied for precise IFB rate calibration, for example using long data set and even without

known station coordinates. It can also be run for fast and even real-time calibration. In this case, the particle

procedure can be stopped, as soon as the estimated IFB rate converges, for example its STD is smaller than a

threshold value, and then the IFB rate value can be fixed to perform precise positioning with ambiguity-fixing.

Certainly, a procedure to monitor its possible changes should be included in the data processing as part of the

quality control. This procedure can be shown as a flow chart in Fig. 4.6.

Fig. 4.6 Flow chart of the procedure for IFB estimation, where 𝑠𝑡𝑑𝑡ℎ𝑑 is the STD threshold value

4.4 Results and Analysis

The performance of the new approach presented in Section 4.3 is investigated in this section using the same three

data sets described in Section 4.2. It must be pointed out that in all processing for IFB rate estimation or

GLONASS ambiguity resolution only dual-frequency GLONASS data alone is employed and GPS data is

processed independently for comparison. In the following, three major results of the experimental validation are

presented and analysed in details in order to confirm whether the new approach could be applied to real-time

applications. The first part is to evaluate the convergence and accuracy of the IFB rate estimation. Then, the

performance of the baseline processing with integer ambiguity resolution is investigated with IFB rate estimation

procedure or with estimated IFB rate. The last part is to analyse the time consumption of the new approach.

4.4.1 Phase IFB Rate Estimation

For the three baselines, the estimated IFB rates of all epochs are drawn in Fig. 4.7. For the zero-baseline

(Fig. 4.7a), the mean value after convergence is -0.0017 mm/FN, with STD of 0.16 mm/FN. For baseline

KOSG-KOS1 (Fig. 4.7b), the mean value after convergence is -24.9 mm/FN, with STD of 0.36 mm/FN. For the

kinematic baseline REF6-AIR5 (Fig. 4.7c) the estimated bias is -0.043 mm/FN with STD of 0.65 mm/FN.

Fig. 4.8 shows the convergence process of the IFB rate and its STD for the baselines. It is clear that the estimated

IFB rate converges quickly and it needs maximally about three minutes to become stable. By the way, it could be

even faster if a better initial value is available. However, there are obviously some fluctuations of about few

mm/FN for baseline KOSG-KOS1 and the kinematic baseline as shown in Fig. 4.7b and Fig. 4.7c, respectively.

This is probably caused by inaccurate modelling and improper quality controls, as KOSG-KOS1 are equipped

with different type of antennas and REF6-AIR5 is in kinematic mode.

𝑤𝑘𝑖=𝑤𝑘𝑖

∑𝑤𝑘𝑖

𝑁

𝑖=1

Update: 𝑤𝑘𝑖=𝑤𝑘−1

𝑖𝑝(𝒃

𝑘|𝑥𝑘𝑖)=𝑤𝑘−1

𝑖𝑅𝑎𝑡𝑖𝑜𝑘𝑖

Initialize: {𝑥0𝑖,𝑤0𝑖}𝑖=1

𝑁,with 𝑤0𝑖=1

𝑁

Output solution: 𝑥𝑘≈∑𝑤𝑘𝑖𝑥𝑘𝑖𝑁

𝑖=1

𝑃𝑘≈∑(𝑥𝑘𝑖−𝑥𝑘)2𝑤𝑘𝑖𝑁

𝑖=1

Predict: 𝑥𝑘+1

𝑖=𝑥𝑘𝑖+𝜖𝑘𝑘+1

Resampling:

{𝑥𝑘𝑖,𝑤𝑘𝑖}𝑖=1

𝑁 with 𝑤𝑘𝑖=1

𝑁

sqrt(𝑃𝑘)<𝑠𝑡𝑑𝑡ℎ𝑑

Yes

𝑥𝑘 is fixed and

IFB rate is set as

known value

Fig. 4.7 Estimated IFB rates using the particle filter for the zero-baseline (a), KOSG-KOS1 (b) and REF6-AIR5

(c)

Fig. 4.8 Convergence of the estimated IFB rate (solid line) and STD (dash lines) for zero-baseline (a),

KOSG-KOS1 (b) and baseline REF6-AIR5 (c), the sampling rate is 5s, 30s, 1s, respectively

In order to further investigate the convergence time for the IFB rate parameter, the data of the three baselines are

processed in short sessions. The first session starts at the data beginning and then moves forwards with a

step-size of 20 epochs, leading to about 110 to 150 sessions for each baseline. The processing of each session

keeps going until the IFB rate converges with STD smaller than a threshold which is set to 2 mm/FN for all

baselines in this study.

Fig. 4.9a to Fig. 4.9c show the convergence process of the IFB rate of all sessions for the three baselines,

respectively. Each line shows the IFB rate estimates for a single session and ends at the epoch meeting the

convergence criteria. For the sake of clarity, the end point is marked with a star symbol. The statistics of

convergence times are presented in Table 4.1. From the statistics, for the zero-baseline and the kinematic

baseline, the new approach could obtain a convergent IFB rate with 10 epochs of data collected within 30 s.

However, for the KOSG-KOS1 baseline with a sampling rate of 30 s four minutes of data is needed. There is a

big difference in the averaged time needed for the convergence, but the numbers of needed epochs are closer to

each other. This can be clearly seen from Fig. 4.10 showing the distribution of the number of epochs needed.

Therefore, one possibility is that a certain number of epochs are necessary for the particle filter to obtain a

convergent IFB rate. Another reason could be that KOSG-KOS1 is the only baseline with different receivers, so

that it needs longer time for convergence. Since KOSG and KOS1 do not provide 1 s high rate data, another two

collocated IGS stations STR1 and STR2 in Australia with different type of receivers are chosen for validation.

The data processing is carried out in the same way as for the previous three baselines, but using sample rates of

1 s, 5 s and 30 s. The results shown in Fig. 4.11 confirm that a certain number of epochs is needed for the new

approach, so for baselines with different receivers IFB rate could also be precisely estimated using 1 Hz data

collected within 30 s. The quick convergence and the stable estimation provide a chance for a field and real-time

calibration of the IFBs for instantaneous ambiguity resolution without any a priori information.

For a relatively long baseline, this approach still works if only other error sources can be neglected or removed

so that models in (4.3) are accurate. However, if the effects of other error sources are significant but are not

removed, they will affect the accuracy of the estimated IFB rate. An additional baseline KOSG-APEL with

length of 11km is also tested with the data collected on DOY 351 of 2013. Both stations were equipped with

LEICA GR25 receivers, but with antennas LEIAR25.R4 LEIT and AOAD/M_B NONE, respectively. The three

dimensional RATIO distribution and the estimation results are presented in Fig. 4.12. It is clear that not all

epochs can provide excellent results, due to other error sources which exist but are not considered in the model

(4.3).

Table 4.1 Statistics of the convergence time for the IFB rate estimated by the new approach for the three

baselines

Baseline

Sampling

Rate (s)

#Total

Epochs

#Sessions

/Removed

Time(s)/Epochs For Convergence

Max

Min

Mean

Zero-Baseline

2270

113/ 0

40/ 8

15/ 3

24/ 5

KOSG-KOS1

2880

144 / 1

630/ 21

120/ 4

240/ 8

REF6-AIR5

2663

133/ 6

29/ 29

6/ 6

11/ 11

Fig. 4.9 Convergence process of the estimated IFB rate versus the number of epochs for the zero-baseline (a),

KOSG-KOS1 (b) and REF6-AIR5 (c). The star symbols denote the converging point. The sampling rate is 5 s,

30 s, and 1 s, respectively

Fig. 4.10 Statistics of the epochs needed for the convergence of IFB rate for the zero-baseline (a), KOSG-KOS1

(b) and REF6-AIR5 (c). The sampling rate is 5 s, 30 s and 1 s, respectively

Fig. 4.11 The statistics of the epochs needed for the convergence of IFB rate for baseline STR1-STR2 with

different type of receivers. The sampling rate is 1 s (a), 5 s (b) and 30 s (c), respectively

Fig. 4.12 Three-dimensional RATIO distribution (left) and the estimation results (right) for baseline

KOSG-APEL

4.4.2 Computational Efficiency

From the algorithm of the new approach, at each epoch about 200 particles must be tested for ambiguity

resolution. Therefore, the computational time at each epoch is of course a critical concern, especially for

real-time applications. In order to give an estimation of the computation efficiency, the computational time at

each epoch is recorded in a personal computer (PC) with a processor of 2.8 GHz and plotted in Fig. 4.13. The

computation time is somehow correlated with the number of satellites at the epoch. Generally, it could be

completed within 1 s for most of the epochs. The computational time could be reduced significantly if a better

initial IFB rate is available, as fewer particles are needed. For example, it takes about 0.17 s if the searching is

within [-0.04, 0] m/FN.

Fig. 4.13 Computational time for a single epoch including the particle filter for baseline KOSG-KOS1. The

upper thin dash line and the thick solid line close to the bottom are for the search interval of [-0.10, 0.10] m/FN

and [-0.04, 0] m/FN, respectively. The number of satellites is also plotted

4.5 Regularized Approach

The RPF described in section 3.7 changes the discrete distribution into continuous distribution, from which new

samples are generated. This filtering method solves the diversity-loss problem and is employed in this section so

that the noise level in the prediction model can be set freely.

4.5.1 Problem in State Noise Setting

In the case of the precise calibration, we may not only be interested in the short convergence time, but also in a

certain precision that the estimated IFB rate can reach, such as 1 mm/FN or even smaller. If the observation time

is longer, the estimated value should be more and more accurate and the threshold can be set to a smaller value.

However, if the state noise level in (4.4) is high, such as 1 mm/FN, the STD of the estimated IFB rate value may

never reach the threshold. Actually, if the IFB rate is stable, the state noise level should be set to a very small or

even to zero, but the low state noise cannot solve the diversity loss problem in particle filter.

The effect of the diversity loss is shown by experiments with the data of baseline KOSG-KOS1 on DOY 048 of

2014. The number of particles is still set to 200 and these particles are randomly sampled over initial interval

[-0.1, 0.1] m/FN at the beginning. As the initial interval is wide, only fewer particles are located close to the true

IFB rate value -24.9 mm/FN. In the following epochs, the resampling step deletes all other particles and only the

particles with large weights are left. If the diversity of these particles is poor, the PDF function cannot be well

represented by the particles and the estimated values can be biased.

The state noise in model (4.4) is set to different noise levels with STDs of 0.001 m/FN, 0.0001 m/FN,

0.00001 m/FN and 0 m/FN, respectively. With each noise level the IFB rates for the whole day are estimated. As

the results are affected seriously by the initial values of particles which are randomly generated, the IFB rate is

estimated 20 times to obtain a general view. The 20 results are plotted together in Fig. 4.14.

Obviously, when the STD of the system noise is set to 0.001 m/FN, the results fluctuate due to the error sources

in GNSS observations and the STD of the estimated IFB rate can be as large as 0.004 m/FN. When the STD of

the system noise decreases to 0.0001 m/FN, the results are smoother with small STD but still with significant

fluctuations. In a further step by setting the STD of the state noise to 0.00001 m/FN, the effects of the diversity

loss can be clearly observed at the beginning of the 24 hours. It takes long time for the particles to move towards

the true value even though the STDs have been very small.

When the system noise is set to value zero, only deletion and duplication of these particles are carried out, i.e. no

new sample values are generated even though the weights of particles change due to weight update at each epoch.

In this case, the sample values do not change during the filtering and the bias in the estimated IFB rate is kept to

the end of the data except for some jumps, which are caused by the particle migration among the initial values

generated at the beginning. Obviously, when the system noise level is very low, the results cannot converge to a

stable value within acceptable convergence time.

The RPF can solve the problem of diversity loss and allows state noise to be low. The utilization of RPF in IFB

rate estimation will be investigated in section 4.5.2. The principle of the regulation method has been simply

described in section 3.7 and can be implemented by the procedure as follows. Firstly, generate 𝜖𝑖 from the kernel

(3.35), which is carried out using a procedure similar to the resampling procedure described in section 3.5.

Secondly, calculate the optimal bandwidth by (3.36), as well as the variance 𝑆 in the case of one dimension, and

then update each particle by (3.38). This procedure can be located after the resampling step in the approach

described in section 4.3.

4.5.2 Experiment with Regularized Approach

By adding the procedure of regularization, the problem of diversity loss is well solved and hence it is not

necessary to consider the particle diversity when setting the state noise level. As the IFB rate is almost a constant

value during one day, the state noise can be set to a very small or even zero value. The same data in section 4.5.1

are employed here again to test this approach.

The same four levels of state noise are still used in the experiment and the IFB rate values are calculated still 20

times for each of them. Results are shown in Fig. 4.15 where all 20 calculations converge to similar values

quickly with all four different levels. The effects of diversity loss in Fig. 4.14 cannot be observed in Fig. 4.15.

The STD values keep decreasing as the state noise becomes lower, which can satisfy the requirement of more

precise IFB rate estimate.

Fig. 4.14 Estimated IFB rate (left) and the STD (right) for twenty calculations by particle filter with the state

noise level 𝜎𝑛 set to different values

4.6 Adaptive Method for Setting the Number of Particles

There are no optimal methods to select the number of particles in particle filter until now. The number of

particles is usually decided by experiences depending on conditions such as required precision and hardware

condition. The larger the number is, the more accurately the PDF can be represented. However, as the number

increases, the calculation burden becomes heavier. Hence, reducing the number of particles without largely

degrading the accuracy is an interesting topic. Fox (2003) adapted the sample number via Kullback-Leibler

distance. The sample number will be small if the density is focused on a small part of the state space and will be

large if it is not. This adaptive approach is most advantageous in lower dimensional state spaces when the

complexity of the posterior changes drastically over time. Closas and Fernández-Prades (2011) also proposed an

adaptive method to reduce the partile number. The particles which are close to their neighboring particles are

discarded, but new particles will be generated if the innovation error is larger than a threshold.

Fig. 4.15 Estimated IFB rate (left) and the STD (right) for twenty calculations by RPF with the state noise level

𝜎𝑛 set to different values.

In the IFB rate estimation problem, the PDF determined from RATIO values is very narrow and the prediction

model is very accurate. Therefore, after convergence the PDF of the variable can be represented by particles with

even a smaller number, for example value 30 instead of value 200. The smaller number is preferred so that the

computation quantity can be reduced, but the influences of employing a smaller number should be investigated

firstly. The experiment with practical data will be conducted in this section.

With the number of particles set to a fixed value 30, the IFB rate is estimated with data for baseline

KOSG-KOS1. To reduce the effect of the randomly generated initial particle values and the state noise so that

we can have a general view, the IFB rate is estimated 20 times and results are shown in Fig. 4.16. Besides, the

results of the 20 calculations with number of particles 200 are also presented in Fig. 4.16. The precisions for both

cases are obviously similar after convergence, but the convergence times at the beginning are different, which

can be much longer for number of particles set to 30.

This is because the IFB rate is unknown at the beginning and therefore the particles are scattered over the whole

interval and jittered for a little distance every epoch to detect the true IFB rate value. If the particles are too few,

this detecting process will take more time, which is obviously not preferred. A better way is to adapt the number

of particles to their STD, which will be investigated in the following subsection.

Fig. 4.16 Estimated IFB rate (left) and the STD (right) with the number of particles N set to different values.

4.6.1 Proposed Adaptive Method

The particles are to approximate the PDF of the state variable. So when the PDF of the state variable is very

narrow compared with the accuracy we need, a relatively larger approximation error will be acceptable. In both

the adaptive methods proposed by (Fox 2003, Closas and Fernández-Prades 2011), the number of particles

becomes samll when the PDF focuses on a small part of the state space. However, once the PDF covers a large

part of the state space, more particles will be needed. With this idea, a similar adaptive method for the number of

particles with STD parameter will be proposed by using the characteristics of the IFB rate.

In the carrier phase IFB rate estimation, the state variable is stable, so is the PDF. Besides, as the IFB rate is

constrained within [-0.1, 0.1] m/FN, no maximum threshold value is needed. Assuming the number of particles

for each unit STD is a certain value so that the PDF can be represented by particles with a certain density, a

simpler function can be designed to tune the number of particles according to the STD value. The function can

be expressed by

{𝑁=int(𝑆𝑇𝐷∙𝑛) if 𝑁>𝑁0

𝑁= 𝑁0 otherwise, (4.7)

where 𝑁0 is the lower bound of the number of particles; 𝑛 is the number of particles for each STD unit. This

function is implemented in the resampling step, where the number of particles is not set to a fixed value but

controlled by (4.7).

4.6.2 Experiment with the Adaptive Method

The same data employed to draw Fig. 4.16 are used here to compare the method with fixed and controlled

number of particles. Still 200 samples are firstly generated randomly and then this value is adapted by function

(4.7) with minimum number 𝑁0 set to value 30 in the calculation. The number of particles per millimetre of STD

is set to value 6 which leads to around 200 particles when these particles are distributed randomly over

[-0.10, 0.10] m/FN. The IFB rate values are calculated 20 times with controlled numbers and results for the

24 hours are presented in Fig. 4.17. The results of the first 15 minutes in Fig. 4.17 are zoomed in and shown in

the top panels of Fig. 4.18, alongside the results with the number of particles fixed to 200 in the bottom panels of

Fig. 4.18. It is obvious that the performances with the number of particles set to controlled value and constant

value 200 are similar regarding the precision.

The number of particles along time are plotted in Fig. 4.19 for the data of 24 hours (left) and the zoom-in of the

first 15 minutes (right). It can be observed that for the approach with adaptive method, the number of particles is

large at the beginning as the variance of the particles is large. Then, the particles converge quickly and the

number of particles also decreases quickly. When the number of particles calculated from STD is smaller than

the minimum value 30 in function (4.7), it is set to the minimum value 30. If the STD becomes larger, the

number will increase again.

The corresponding computation time is investigated and presented in Fig. 4.20 for the data of 24 hours (left) and

the first 15 minutes (right), respectively. It is clear that the approach with controlled number of particles takes

much less time after convergence. The average computation time is 0.84 s for 200 particles but 0.11 s for

particles of controlled number, but the STDs of the estimated IFB rate are the same for both cases, 1.2 mm.

Fig. 4.17 Estimated IFB rate (left) and the STD (right) for 20 calculations with controlled numbers for 24 hours

Fig. 4.18 Estimated IFB rate (left) and the STD (right) for 20 calculations with number of particles N set to

controlled value and constant value 200 for the first 15 minutes

4.7 Estimated Phase IFB Rates and their Characteristics

The IFB in GLONASS data processing is stable over time and the value is dominated by receiver type. As these

characteristics of IFB rate has been investigated with a large number of baselines and receiver types in

(Wanninger 2012), only a few data are employed to confirm the characteristics.

IFB Values for Long Period of Time

In the existing research, the IFB rate in GLONASS data processing is supposed to be stable even in a long run.

Here only two baselines STR2-STR1 and KOSG-KOS1are employed to test the characteristic.

Fig. 4.19 Number of particles N along time for 20 calculations with the method of setting N to controlled value

and constant value 200 for the 24 hours (left) and for the first 15 minutes (right)

Fig. 4.20 Computation time with the number of particles set to controlled value and constant value 200, along

with the number of satellites, for the 24 hours (left) and for the first 15 minutes (right)

The first day for STR2-STR1 is DOY 001 of 2013 and the data of 20 days scattering within the following more

than 700 days are employed. The approach described in section 4.3 is used and results are shown in Fig. 4.21

(top). Even though there were some firmware updates during this time period, the IFB rate is very stable and no

changes are observed.

In the experiment with baseline KOSG-KOS1, the first day is DOY 170 of 2013 and the data of 16 days

scattering in the following 600 days are employed. The same approach is used and results are presented in

Fig. 4.21 (bottom). It can be observed that the IFB rate for KOSG-KOS1 is also very stable for the selected days.

Results for Different Receiver Brands

The factors which affect the IFB value include the receiver type, antenna type, temperature, cable length and so

on. Some of these effects are trivial and even cannot be observed, such as temperature, but the receiver type is

dominant. An experiment is designed to test it in this section.

The data of three baselines for two days on DOY 191 of 2014 and DOY 036 of 2015 are employed and the IFB

rates are listed in Table 4.2. The receiver types, which are different for each baseline but the same on the two

days, are shown in Table 4.3, along with the antenna types and radomes. Two of the four receivers are from the

same manufacturer for every two baselines, so if the IFB rate values for receivers of the same type are exactly

the same and receiver type is the dominant factor, the sum of the three IFB rate values should be zero. In the

experiment, the sums on DOY 191 of 2014 and 036 of 2015 are 1.3mm and 1.1 mm respectively, which are

pretty small. This confirms the known characteristics of IFB.

The IFB rate is stable over long time, so the IFB rate of the baseline can be estimated once and used in all the

data processing of the baseline, or even used in the other baselines equipped with receivers of the same two types.

Fig. 4.21 IFB rates in GLONASS data processing for baselines STR2-STR1 (top) and KOSG-KOS1 (bottom)

within about two years

4.8 Summary

Due to the existence of the IFB in the carrier phase observations, GLONASS integer ambiguity resolution

encounter difficulties in relative positioning when different types of receivers are employed, especially for

real-time applications. Almost all existing methods estimate the IFB rate together with the ambiguities, which

takes long time to converge and the situation cannot be improved very much by using simultaneous GPS data.

During this process, GLONASS ambiguity resolution is hardly available, which is a critical obstacle for fast and

precise positioning.

In this chapter, It is demonstrated that the integer ambiguity resolution is very reliable even using only a few

epochs of GLONASS observations when IFB rate is precisely known. Moreover, the closer the IFB rate is to the

true value, the better the fixing performance will be. Therefore, RATIO of the ambiguity fixing can be used as a

reliable index to qualify a given IFB rate value. Based on this fact, a new method is developed to estimate IFB

rate by methodology of particle filter. From the outcome of the experimental evaluation, the new method enables

for the first time the quick estimation of IFB rate using only GLONASS data of few epochs and without an a

priori value.

Afterwards, the particle filter approach is improved in two aspects to satisfy the requirement in special cases.

Firstly, the RPF is introduced to enable the prediction model to have small state noise level so that higher

precision with observations of longer time periods can be achieved. Secondly, a function is designed to tune the

number of particles according to the particles’ STD so that the computation time in the tracking period can be

largely reduced.

Table 4.2 Lengths and estimated IFB rates for the three baselines on DOY 191 of 2014 and DOY 036 of 2015

Baseline

length (m)

F-ISB (m/FN)

DOY 191 of 2014

DOY 036 of 2015

STR2-STR1

0.0317

0.0316

KOSG-KOS1

814

-0.0247

-0.0249

TLSG –TLSE

1266

-0.0057

-0.0056

Table 4.3 Receivers and antennas for the three baselines

Station ID

Receiver type

Antenna type and the radome

STR1

LEICA GRX1200GGPRO 8.71/3.822

ASH701945C_M NONE

STR2

TRIMBLE NETR9 4.70

TRM59800.00 NONE

KOS1

SEPT POLARX4 2.5.2

LEIAR25.R3 LEIT

KOSG

LEICA GRX1200GGPRO 8.20/3.019

AOAD/M_B NONE

TLSE

TRIMBLE NETR9 4.85

TRM59800.00 NONE

TLSG

SEPT POLARX4TR 2.5.2

TRM59800.00 NONE

5 Online Inter-System Ambiguity Resolution Based on Online F-ISB

Estimation

This chapter investigates the F-ISB estimation for inter-system DD-ambiguity resolution. Firstly the existing

F-ISB estimation methods are introduced in section 5.1 and the multi-GNSS models employed in this chapter are

presented in section 5.2. Then the integration of GPS L1 and Galileo E1, GPS L1 and GLONASS L1, as well as

GPS L1 and BDS B1, are taken as typical examples to investigate the RATIO distribution and the fixed solutions

with systems of the same frequency and different frequencies in section 5.3. Afterwards, in section 5.4 the new

approach based on particle filter is proposed to estimate the F-ISB precisely with the half-cycle problem solved

by the cluster analysis method. The F-ISB estimated results with the new approach are shown and analyzed in

section 5.5, while the ISB Characteristics in Multi-GNSS Integration are investigated in section 5.6.

5.1 Existing Phase ISB Estimation Methods

The inter-system ambiguity fixing is critical in the case of few satellites in view. To recover the integer nature of

the inter-system DD-ambiguities, the ISB or F-ISB has to be known or estimated. The investigations of inter-

system ambiguity fixing in recent years mainly focused on the overlapping frequencies based on model (2.17).

The equation set obtained from model (2.17) with ISB parameter is rank-deficient and the deficiency number

equals to the number of ISB parameters for carrier phase. In addition, only the F-ISB in model (2.17) affects the

results and hence these investigations aim to remove the rank-deficiency so that the F-ISB can be determined.

Odijk and Teunissen (2013a, 2013b) proposed to solve the rank-deficiency problem by selecting two reference

satellites for the two systems, respectively. The ambiguity parameter in the inter-system model between the two

reference satellites includes the integer ambiguity and the ISB. Therefore, this ambiguity parameter can be

utilised to remove the F-ISB from other inter-system models so that their DD-ambiguities have integer nature

again and can be fixed.

In this method, the model within GPS is like (2.10), while the inter-system DD-model between GPS reference

satellite and Galileo reference satellite is written as (Odijk and Teunissen 2013a, Odijk and Teunissen 2013b)

𝑃𝑎𝑏

𝑖𝑘,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑘,𝐺𝐸+𝑑𝑎𝑏

𝐺𝐸+𝜀𝑎𝑏

𝑖𝑘,𝐺𝐸, (5.1a)

𝜆𝛷𝑎𝑏

𝑖𝑘,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑘,𝐺𝐸+𝜇𝑎𝑏

𝐺𝐸+𝜉𝑎𝑏

𝑖𝑘,𝐺𝐸, (5.1b)

where 𝜇𝑎𝑏

𝐺𝐸=𝜇𝑎𝑏

𝐺𝐸+𝜆𝑁𝑎𝑏

𝑖𝑘,𝐺𝐸=𝜇𝑎𝑏

𝐺𝐸+𝜆(𝑧𝑎𝑏

𝐺𝐸+𝑁𝑎𝑏

𝑖𝑘,𝐺𝐸) contains both the ISB and the DD-ambiguity parameter

of the two reference satellites; 𝜇 is the ISB value; 𝜇 is the F-ISB; 𝑧 is the remaining part of ISB which is integer

multiple of wavelength. 𝑖 and 𝑘 are the GPS reference satellite and the Galileo reference satellite, respectively.

The other inter-system DD-models are written as

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝑑𝑎𝑏

𝐺𝐸+𝜀𝑎𝑏

𝑖𝑗,𝐺𝐸, (5.2a)

𝜆𝛷𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝜇𝑎𝑏

𝐺𝐸+𝜆𝑁𝑎𝑏

𝑘𝑗,𝐸𝐸+𝜉𝑎𝑏

𝑖𝑗,𝐺𝐸, (5.2b)

where 𝑗 refers to another satellite; 𝜇𝑎𝑏

𝐺𝐸 is the same as in (5.1b). In this way, the F-ISB in the inter-system models

except for the model of the two reference satellites can be fixed as integers. Afterwards, 𝜇𝑎𝑏

𝐺𝐸 is estimated along

with the coordinate parameters. Even though 𝜇𝑎𝑏

𝐺𝐸 may be larger than one wavelength, only the fractional part

𝜇𝑎𝑏

𝐺𝐸 matters as the part 𝑧𝑎𝑏

𝐺𝐸 lump together with the integer ambiguities.

If the F-ISB parameter is known, the ambiguity fixing can benefit from the inter-system models as they provide

one more independent equation with integer DD-ambiguity. This additional equation can be crucial under the

condition that only few satellites from multi-GNSS are observed. However, when the F-ISB is parameterized in

the models, the additional equation includes new unknown parameters and thus cannot help to improve the

solution (Odijk and Teunissen 2013a).

Paziewski and Wielgosz (2015) proposed to solve the rank-deficiency problem by assigning the F-ISB parameter

zero value with the STD set to half phase cycles to constrain the ISB. The models for this method can be

expressed by

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝑑𝑎𝑏

𝐺𝐸+𝜀, (5.3a)

𝜆𝛷𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝜇𝑎𝑏

𝐺𝐸+𝜆𝑁

𝑎𝑏

𝑖𝑗,𝐺𝐸+𝜀, (5.3b)

𝜇𝑎𝑏

𝐺𝐸=𝜇𝑎𝑏

𝐺𝐸0, (5.3c)

where 𝑁

𝑎𝑏

𝑖𝑗,𝐺𝐸=𝑧𝑎𝑏

𝐺𝐸+𝑁𝑎𝑏

𝑖𝑗,𝐺𝐸; 𝜇𝑎𝑏

𝐺𝐸0 is an a priori F-ISB value and is always zero with variance equals to the

square of half cycle (Paziewski and Wielgosz 2015). In this way, the best case is that the F-ISB is just zero.

However, once the true F-ISB value is not zero, such as half wavelength, this approach obviously cannot

outperform the strategy only fixing the intra-system DD-ambiguities.

5.2 Multi-GNSS Mathematic Models for the New Approach

Four systems are included in the experiments about F-ISB estimation in this thesis, including GPS, GLONASS,

Galileo and BDS. The functional models within each system are presented together as follows

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐺+𝜀𝑎𝑏

𝑖𝑗,𝐺𝐺, (5.4a)

𝜆𝐺𝛷𝑎𝑏

𝑖𝑗,𝐺𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐺+𝜆𝐺𝑁𝑎𝑏

𝑗,𝐺−𝜆𝐺𝑁𝑎𝑏

𝑖,𝐺+𝜉𝑎𝑏

𝑖𝑗,𝐺𝐺; (5.4b)

𝑃𝑎𝑏

𝑖𝑗,𝐸𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐸𝐸+𝜀𝑎𝑏

𝑖𝑗,𝐸𝐸, (5.5a)

𝜆𝐸𝛷𝑎𝑏

𝑖𝑗,𝐸𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐸𝐸+ 𝜆𝐸𝑁𝑎𝑏

𝑗,𝐸−𝜆𝐸𝑁𝑎𝑏

𝑖,𝐸+𝜉𝑎𝑏

𝑖𝑗,𝐸𝐸; (5.5b)

𝑃𝑎𝑏

𝑖𝑗,𝑅𝑅=𝜌𝑎𝑏

𝑖𝑗,𝑅𝑅+𝜀𝑎𝑏

𝑖𝑗,𝑅𝑅, (5.6a)

𝜆𝑗,𝑅𝛷𝑎𝑏

𝑗,𝑅−𝜆𝑖,𝑅𝛷𝑎𝑏

𝑖,𝑅=𝜌𝑎𝑏

𝑖𝑗,𝑅𝑅+𝜆𝑗,𝑅𝑁𝑎𝑏

𝑗,𝑅−𝜆𝑖,𝑅𝑁𝑎𝑏

𝑖,𝑅+𝜉𝑎𝑏

𝑖𝑗,𝑅𝑅; (5.6b)

𝑃𝑎𝑏

𝑖𝑗,𝐶𝐶=𝜌𝑎𝑏

𝑖𝑗,𝐶𝐶+𝜀𝑎𝑏

𝑖𝑗,𝐶𝐶, (5.7a)

𝜆𝐶𝛷𝑎𝑏

𝑖𝑗,𝐶𝐶=𝜌𝑎𝑏

𝑖𝑗,𝐶𝐶+𝜆𝐶𝑁𝑎𝑏

𝑗,𝐶−𝜆𝐶𝑁𝑎𝑏

𝑖,𝐶+𝜉𝑎𝑏

𝑖𝑗,𝐶𝐶, (5.7b)

where G, E, R and C refer to GPS, Galileo, GLONASS and BDS, respectively.

The inter-system models can also be constructed. Due to the easier and more accurate ambiguity resolution, the

integration of the overlapped frequencies is preferred by researchers. In this study the inter-system models with

slightly different frequencies are also included. The clock bias parameters can be eliminated in the inter-system

model, but the hardware delays in both code pseudorange and carrier phase observations stay because they are

likely to be different for different receivers. The inter-system functional models including GPS observations are

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝑑𝑎𝑏

𝐺𝐸+𝜀𝑎𝑏

𝑖𝑗,𝐺𝐸, (5.8a)

𝜆𝐺𝛷𝑎𝑏

𝑖𝑗,𝐺𝐸=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐸+𝜇𝑎𝑏

𝐺𝐸+𝜆𝐸𝑁𝑎𝑏

𝑗,𝐸−𝜆𝐺𝑁𝑎𝑏

𝑖,𝐺+𝜉𝑎𝑏

𝑖𝑗,𝐺𝐸; (5.8b)

𝑃𝑎𝑏

𝑖𝑗,𝐺𝑅=𝜌𝑎𝑏

𝑖𝑗,𝐺𝑅+𝑑𝑎𝑏

𝐺𝑅+𝜀𝑎𝑏

𝑖𝑗,𝐺𝑅, (5.9a)

𝜆𝑗,𝑅𝛷𝑎𝑏

𝑗,𝑅−𝜆𝐺𝛷𝑎𝑏

𝑖,𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝑅+𝜇𝑎𝑏

𝐺𝑅+𝜆𝑗,𝑅𝑁𝑎𝑏

𝑖,𝑅−𝜆𝐺𝑁𝑎𝑏

𝑖,𝐺+𝜉𝑎𝑏

𝑖𝑗,𝐺𝑅; (5.9b)

𝑃𝑎𝑏

𝑖𝑗,𝐺𝐶=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐶+𝑑𝑎𝑏

𝐺𝐶+𝜀𝑎𝑏

𝑖𝑗,𝐺𝐶, (5.10a)

𝜆𝐶𝛷𝑎𝑏

𝑗,𝐶−𝜆𝐺𝛷𝑎𝑏

𝑖,𝐺=𝜌𝑎𝑏

𝑖𝑗,𝐺𝐶+𝜇𝑎𝑏

𝐺𝐶+𝜆𝐶𝑁𝑎𝑏

𝑗,𝐶−𝜆𝐺𝑁𝑎𝑏

𝑖,𝐺+𝜉𝑎𝑏

𝑖𝑗,𝐺𝐶. (5.10b)

The other inter-system models, such as the model between Galileo and BDS, are correlated with the models from

(5.4) to (5.10). All of these models are included in the frame of the general model (2.22). Therefore, they can be

solved with the procedure described in section 2.4.

The stochastic models for the equations from (5.4) to (5.10) are calculated according to the rules of variance

propagation with the variances of the raw observations calculated by (2.9). The constant parameters in (2.9) are

set as the same for all the systems, while these constant parameters for pseudorange observations are 100 times

of these for carrier phase observations.

The code pseudorange ISB estimation has little effects on the carrier phase ISB, because the weights for carrier

phase observations are much larger than these of code pseudorange observations. Therefore, in order to reduce

the number of unknown parameters in each epoch, the ISB parameter for code pseudorange is estimated in

advance by only code pseudorange observations. Then the code pseudorange ISB parameter is treated as a

known value. Therefore, the unknown parameters in (2.17) and (2.18) include only the unknown coordinates, the

ambiguities and the carrier phase ISB. This method simplifies the model by removing the code pseudorange ISB

parameters and is utilized in the following experiments in this study. This ISB can also be parametrized in the

inter-system models, and experiment results show little difference.

In the inter-system phase models of different frequencies, such as inter-system phase model between GPS L1

and GLONASS L1 observations, the initial SD-ambiguities are calculated with (2.20). As described in section

2.3.3, the right-hand side of model (2.21) includes the clock error, the ISB values for both code pseudorange and

carrier phase observations, as well as the atmospheric effects which are not considered in short baselines. Thus,

both code pseudorange ISB and carrier phase ISB are required in the calculation of the initial SD-ambiguities.

The calculation of ISB for code pseudorange is not a problem as ISB parameter can be determined easily from

single point positioning (SPP), but this is not the case for carrier phase. In view of that the code pseudorange

observations are usually employed to calculate the initial values of SD-ambiguities in GLONASS only data

processing, the ISB difference between code pseudorange and carrier phase should be small compared with the

magnitude of ISB values themselves. So the two ISB terms in (2.21) are assumed to be the same and therefore

can cancel each other. This means that the initial values of the SD-ambiguities are calculated from code

pseudorange observations directly without code pseudorange ISB correction.

5.3 Relationship between RATIO and F-ISB

The RATIO distribution corresponding to different F-ISB values, as well as the fixed solutions is investigated in

this section. The employed data includes GPS L1 and Galileo E1 with the same frequency, GPS L1 and

GLONASS L1 with different frequencies and with FDMA technique, as well as the data of GPS L1 and BDS B1

with different frequencies but both with CDMA techniques.

5.3.1 RATIO versus ISB of GPS L1 and Galileo E1

GPS L1 and Galileo E1 have the same frequency, hence DD-ambiguities can be directly formed in the

DD-observation model. But for the sake of convenience, they are solved in the form of general models (2.22).

Employed Data

The baseline employed here for the study of ISB characteristics is TLSG-TLSE, whose receiver types and

antenna types have been introduced in section 4.7. The data used are collected on DOY 001 of 2015 with epoch

interval of 30 seconds. The sky plots of GPS and Galileo at station TLSE are drawn in Fig. 5.1. The numbers of

satellites for the baseline are presented in Fig. 5.2. The Galileo system is under construction and therefore the

maximum number of observed Galileo satellites is only 3 with elevation mask of 10 degrees.

Relationship between RATIO and ISB

Although the F-ISB is smaller than one wavelength, investigation of F-ISB on a much wider initial interval may

help to confirm the periodic characteristics of RATIO distribution. Firstly, the initial interval [-0.20, 0.20] m, i.e.

about [-1.0, 1.0] cycles, is evenly sampled with the sampling interval of 1 mm, which results in 400 F-ISB

candidates in total. For the epochs having Galileo observations, RATIO values corresponding to these F-ISB

samples are calculated. GPS L2 observations are also employed to enhance the single-epoch model. The RATIO

values for the first calculated epoch are presented in Fig. 5.3a, while values for all epochs are drawn in Fig. 5.3b.

Two peaks in Fig. 5.3a and two ridges in Fig. 5.3b with relatively large RATIO values can be observed.

To investigate samples over a wider interval, the initial interval for F-ISB is then expanded to [-1.0 m, 1.0 m]

and about 2000 samples are generated with sampling interval 1 mm. In this case, most of the F-ISB values

actually include the integer multiple of a wavelength. The 3D RATIO map is presented in Fig. 5.4a and the

average values of all epochs are shown in Fig. 5.4b. The periodic characteristic can be clearly observed.

From Fig. 5.3 and Fig. 5.4, there are many ISB samples which can provide a fixed solution with RATIO value

larger than 3. Actually, with any F-ISB within about [-3 cm, 3 cm] around each peak in Fig. 5.4b a fixed solution

can be achieved. Certainly, most of them are biased and thus result in a contaminated fixed solution.

Fig. 5.1 Satellite sky plots of GPS (left) and Galileo (right) for station TLSE on day 001 of 2015

Fig. 5.2 Numbers of GPS and Galileo satellites for baseline TLSG-TLSE on DOY 001 of 2015, with elevation

mask of 10 degrees

Fixed Solution for the Baseline

In order to investigate the accuracy of the fixed solutions, for all the F-ISB candidates the positions of the fixed

solutions are calculated and then compared with the GPS only fixed solution of that epoch. As a typical example,

the differences of all fixed solutions at epoch 1406 are shown in Fig. 5.5. The fixed solutions with a peak RATIO

are indicated in red, while the others in blue. For ISBs without a fixed solution, the differences are not shown. It

is clear that the position differences show a periodic characteristic with respect to F-ISB. In other words, the ISB

values with the same fractional part have the same bias in the fixed solution. The solutions corresponding to the

peaks of the RATIO values overlap very well with solutions of GPS data only, which indicates that the

corresponding F-ISB values are most likely to be the correct value.

Fig. 5.3 Relationship between ISB and fixing RATIO for the first epoch with only one Galileo satellite (a) and

the three-dimensional RATIO distribution for all epochs involving Galileo observations (b) for baseline

TLSG-TLSE

Fig. 5.4 Three-dimensional RATIO distribution of GPS and Galileo integration (a), as well as the average values

along epoch time (b) for baseline TLSG-TLSE

Fig. 5.5 Impact of ISB biases in the baseline fixed solutions of GPS L1 and Galileo E1 integration with respect

to the GPS L1 baseline fixed solution. The data are from baseline TLSG-TLSE at epoch 1406. The F-ISB

sampling interval is 0.001 m

5.3.2 RATIO versus ISB of GPS L1 and GLONASS L1

Unlike GPS L1 and Galileo E1, GPS L1 and GLONASS L1 have different frequencies. The wavelength

difference for them is larger than the wavelength difference between GPS L1 and BDS B1 and therefore the

characteristics of RATIO distribution with ISB parameter can be more clearly presented. Also, the wavelength

difference is not so large that solutions can be improved with initial SD-ambiguities determined from code

pseudorange observations for at least short baselines. Furthermore, both constellations are fully operational at

present and hence can provide more information about the characteristics of the RATIO distribution. In this

section, the employed data and the code pseudorange ISB correction method are firstly introduced, and then the

RATIO distribution with different F-ISB values is investigated.

Employed Data

The baseline employed in this section is composed of IGS stations KOSG and KOS1, which has been introduced

in section 4.2. The data are collected on DOY 048 of 2014 with sampling interval of 30 seconds. The sky plots

of GPS and GLONASS satellites at station KOS1 are given in Fig. 5.6. The baseline is about 814 m long, so the

KOSG station has similar satellite sky plots. The numbers of observed satellites for the baseline with elevation

mask of 10 degrees are presented in Fig. 5.7.

Code Pseudorange IFB Correction

The code pseudorange IFB in the GLONASS observations are corected by the look-up table of code

pseudorange SD-IFB. To generate such a table, firstly, the residuals of the code pseudorange DD-model (5.6a)

with code pseudorange IFB of zero in the post-processing are collected. Then the non-zero DD-residuals are

regarded as DD-IFB and modeled with SD-IFB parameters. The corresponding equation system is

rank-defficient with rank-deficiency one. Assuming that the sum of all the SD-IFBs equals zero, the rank-

deficiency can be removed and the SD-IFBs can be determined sucessfully. This approach for removing the

rank-deficiency follows the same methodology utilised in (Alber et al. 2000) where single-path phase delays of

atmospheric water vapor are obtained from the DD-values in GPS data processing.

The baseline solutions by code pseudorange observations without and with IFB corrections are presented in

Fig. 5.8. The assumption that the sum of the code pseudorange SD-IFB equals to zero does not affect the

GLONASS solution because only the DD-model is used and the effects are eliminated. In the inter-system

DD-models, the bias caused by the zero-sum assumption lumps with the code pseudorange ISB.

Relationship between RATIO and ISB

The RATIO distribution with different F-ISB values in GPS L1 and GLONASS L1 integration is investigated in

this part. Firstly, values within a certain initial interval are sampled as pre-defined carrier phase F-ISB. Then

(5.4), (5.6) and (5.9) are employed to estimate the float solution. The strategy described in section 2.4 is

implemented to calculate RATIO values. The data processing employs single-epoch strategy and the satellite

with highest elevation, either a GPS or a GLONASS satellite, is selected as reference satellite.

The initial interval for F-ISB is firstly selected as [-20, 20] m which is evenly sampled 2000 times with sampling

interval of 2 cm. Because the ISB for code pseudorange and carrier phase observations are considered to have

the same value in the initial SD-ambiguity calculation as described in section 5.2, the RATIO values are

calculated with these samples plus -16.2764 m, which is the approximate ISB estimated from SPP with code

pseudorange observations. The RATIO distribution with different F-ISB values is shown in Fig. 5.9.

It is clear that there is a ridge composed of the local maximum RATIO values in the 3D RATIO distribution in

Fig. 5.9. The averages of these RATIOs along epoch time for 24 hours are drawn in Fig. 5.10a, where the main

characteristic is also the ridge composed of local maximum RATIO values. To observe the details in Fig. 5.10a,

the initial interval [-1, 1] m is sampled 2000 times with the sampling interval of 1 mm. The average values of all

the epochs are presented in Fig. 5.10b, where periodic characteristic with some local maximum RATIO values

can be clearly identified. The local maximum values in Fig. 5.10b have a slope with respect to F-ISB axis. This

is because the ISB value -16.2764 m determined with code pseudorange measurements is an approximate value

and hence the top of the ridge does not exactly correspond to zero F-ISB in both Fig. 5.10a and Fig. 5.10b.

Fig. 5.6 Satellite sky plots of GPS (left) and GLONASS (right) for station KOS1 on DOY 048 of 2014

Fig. 5.7 Numbers of GPS and GLONASS satellites for baseline KOSG-KOS1 on DOY 048 of 2014, with

elevation mask of 10 degrees

Fig. 5.8 Baseline solutions from code pseudorange observations of GLONASS L1 for baseline

KOSG-KOS1without (left) and with (right) code pseudorange SD-IFB correction

Fig. 5.9 Three-dimensional RATIO distribution of GPS L1 and GLONASS L1 integration with inter-system

models for baseline KOSG-KOS1

Fig. 5.10 Average values (a) along epoch time axis of RATIO values in Fig. 5.9 and the zoom-in of the average

values (b), where the red short lines mark the width of the peaks at RATIO of value 3

The periodic characteristic with maximum RATIO values in Fig. 5.10b is easy to explain. Because the

pre-defined F-ISB samples cover a large interval, most of the values include integer multiple of wavelength

which are supposed to lump into the integer ambiguities. When biases in these samples are just the integer

multiple of the wavelength, the DD-ambiguities are supposed to be integers and the corresponding RATIO

values are large, which forms the peaks in Fig. 5.10b.

The ridges in Fig. 5.9 and Fig. 5.10a are because the RATIO values are affected by the biases in the initial

SD-ambiguities. They can be more specifically explained with the peaks in Fig. 5.10b. The distance between

peaks in Fig. 5.10b is 18.87 cm, which is approximately the average value of the wavelength of GPS L1 and

GLONASS L1. This is due to the comparable numbers of satellites for each system and the equivalent accuracy

for their observations. If one of the systems is dominant (i.e. has more satellites or much more accurate

observations), the distance between peaks will change. For example, if there is a full constellation of GPS but

only one GLONASS satellite, the float SD-ambiguity solution will change little for GPS but a lot for GLONASS

and the distances between peaks will be equal to GLONASS wavelength. If there is a full constellation of

GLONASS but only one GPS satellite, the SD-ambiguity in the float solution will change in an opposite way and

the distances between peaks will be equal to GPS wavelength. However, if both the GPS and GLONASS have

full constellations, neither of them will be dominant. In this case, the ISB error in the inter-system model will

affect the SD-ambiguity solutions more seriously, which results in that the distances between peaks is a value

between the two wavelengths and the corresponding RATIO values are smaller. However, if the ISB errors are

very small, the local maximum RATIO values will be large. This is the reason for the existence of the ridge in

Fig. 5.9 and Fig. 5.10a.

Afterwards, an experiment is conducted to verify the above explanations. Firstly, the full GPS constellation and

only one GLONASS satellite are selected to estimate the solution. The corresponding RATIO values in single

epoch strategy for the data of the first hour on DOY 048 of 2014 are presented in the top panel of Fig. 5.11. Then

only one GPS satellite but full GLONASS constellation are selected and the RATIO values are shown in the

bottom panel of Fig. 5.11. The wavelength of the L1 for GLONASS zero channel is 18.71 cm and for GPS is

19.03 cm. This leads to the fact that the plots don’t overlap with each other at the left side of the plots in

Fig. 5.11.

If the approximate ISB value is subtracted form the inter-system models, like the overlapped part at right of the

plots in Fig. 5.11, the F-ISB values estimated from GPS and GLONASS full constellations can be used in other

constellation conditions to achieve the relatively larger local maximum RATIO values.

Fig. 5.11 Comparison of RATIO values in integration of full GPS and GLONASS constellations, full GPS

constellation and only one GLONASS satellite (top), as well as only one GPS satellite and full GLONASS

constellation (bottom)

Fixed Solution for the Baseline

The F-ISB values are still evenly sampled over initial interval of [-20, 20] m with the sampling interval of

0.02 m. As typical examples, the fixed solutions for epoch 323 and epoch 2871 are plotted in Fig. 5.12 and

Fig. 5.13, respectively.

It is clear that the fixed solutions can be obtained with F-ISB values distributed over a very wide range of

interval. In Fig. 7 for epoch 323, the fixed solutions are available with pre-defined F-ISB values distributed all

over the initial interval [-20, 20] m. In Fig. 5.13, the fixed solutions can be achieved at most F-ISB samples,

except for these samples near the boundaries of the interval.

Although the pre-defined F-ISB values are over the interval as long as 40 m on F-ISB axis, the biases of the

estimated fixed solutions are small. At the two endpoints of the drawing in both Fig. 5.12 and Fig. 5.13, the ISB

values have relatively larger biases, which affect the float solutions of the SD-ambiguities and result in relatively

smaller local maximum RATIO values, leading to larger residuals after the DD-ambiguities are fixed as integer.

In this case, the fixed solutions are affected by the residuals and have relatively larger errors. Therefore, the

drawings are not parallel with F-ISB axis in Fig. 5.12 and Fig. 5.13. But even though F-ISB sample values are

over 40 meters, the biases in the fixed solutions, if available, are only several millimetres. The solutions

corresponding to F-ISB values near zero include smaller biases and are preferred.

To observe in details the biases in estimated station solutions in Fig. 5.13, the F-ISB values are sampled over

initial interval [-1, 1] m with sampling interval 0.001 m. Then the biases of the corresponding fixed solutions for

epoch 2871 with respect to GPS L1 fixed solution are presented in Fig. 5.14. It is clear that the fixed solutions

can be obtained with most F-ISB values of each period and the largest biases of these solutions are only several

millimetres.

Fig. 5.12 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration with

respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at epoch 323. The F-ISB

sampling interval is 0.02 m

Fig. 5.13 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration with

respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at epoch 2871. The F-

ISB sampling interval is 0.02 m

Fig. 5.14 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration with

respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at epoch 2871. The F-

ISB sampling interval is 0.001 m

From above, to recover the integer nature of DD-ambiguities in the inter-system model of different frequencies,

the carrier phase ISB value is regarded as the sum of an approximate ISB, which can be estimated in the SPP

with the code pseudorange observations, and an accurate F-ISB value which is smaller than one wavelength. The

approximate ISB value is not only important for achieving larger RATIO values as shown in Fig. 5.10a, but also

for more accurate positioning solutions as shown in Fig. 5.12 and Fig. 5.13.

Besides, the code pseudorange observations are employed to calculate the initial SD-ambiguities in the general

model (2.22), which affects the performance of ambiguity fixing because the frequencies are different. However,

the tolerance on the bias of the initial SD-ambiguities is large here. As marked by the red lines in Fig. 5.10b, the

width of the RATIO peaks at RATIO value of 3 reaches 10 to 12 cm, which value can be seen as the bias

tolerance in the inter-system DD-models. Setting the bias tolerance to 5 cm and considering that the maximum

wavelength difference between GPS L1 and GLONASS L1 is 3.55 mm, as well as model (2.19), the fixed

solutions can be obtained if only the bias in the initial SD-ambiguities is less than 14 m. The condition that the

biases in the SD-ambiguities is less than 14 m cannot guarantee the success of the ambiguity fixing all the time

and can lead to fixed solution only when the estimation is with the assistance of the intra-system models and

sufficient observations. Meindl (2011) described that the uncertainty of the initial SD-ambiguities are supposed

to be better than 7 cycles, i.e. around 1.3 m, so that the DD-ambiguities can be determined better than 0.1 cycles

in (2.19), which may be a little too strict when some intra-system models are included in the estimation. With

short baselines and initial SD-ambiguities calculated from code pseudorange observations, improvement on DD-

ambiguity fixing can be achieved by introducing inter-system DD-ambiguity in GPS L1 and GLONASS L1

integration and the results will be shown in section 7.2.

5.3.3 RATIO versus ISB of GPS L1 and BDS B1

The carrier phases for GPS L1 and BDS B1 have different frequencies, but both of them employ CDMA

technique which is a different situation compared with GPS L1 and GLONASS L1 integration. The frequency

difference is 14.322 MHz between GPS L1 and BDS B1 with the wavelength difference of about 1.7 mm.

Employed Data

The data for baseline TLSG-TLSE on DOY 001 of 2015 are employed here. The sky plot of BDS for station

TLSE are presented in Fig. 5.15, while the sky plot of GPS has been shown in the left panel of Fig. 5.1. The

numbers of satellites for each system along epoch time are drawn in Fig. 5.16. GPS satellites can be observed all

the time, while BDS satellites are available only for part time.

Relationship between RATIO and ISB

The same as before, an initial interval [-20, 20] m is first sampled with sampling interval of 0.02 m. Each sample

shifted by -0.3007 m, the ISB value from SPP with code pseudorange observations, is set as ISB value for the

integration. RATIO values corresponding to these samples are presented in Fig. 5.17a along accumulated epoch

time. The average values of RATIO along epochs are presented in Fig. 5.17b. The RATIO distribution in

Fig. 5.17a vary a lot from epoch to epoch, but it can be found that the ridge characteristic in Fig. 5.17 a and b is

not so obvious as in the case of GPS L1 and GLONASS L1 integration in Fig. 5.9 and Fig. 5.10a. This is because

of two reasons. First reason is that the wavelength difference between GPS L1 and BDS B1 is 1.7 mm, which is

much smaller than 3.55 mm, the maximum wavelength difference between GPS L1 and GLONASS L1. Second

reason is that the BDS satellites are much fewer than GPS satellites as shown in Fig. 5.16, while the number of

GLONASS satellites is comparable with that of GPS satellites as shown in Fig. 5.7. In this case, the approximate

ISB value is not as important as in GPS L1 and GLONASS L1 integration and all the local maximum RATIO

values are similar with F-ISB values on different periods.

Fixed Solution for the Baseline

The position differences between the fixed solutions of the integration and the fixed solutions of GPS L1 at

epoch 1703 are drawn in Fig. 5.18. The differences are similar over the whole interval [-20, 20] m.

To observe the characteristics in detail, 2000 samples of F-ISB over [-1, ..., 1] m with sampling interval of 1 mm

are tested. The calculated RATIO values, as well as the fixed solutions when RATIO values are larger than 3, are

presented in Fig. 5.19a and b, respectively. The reference values in Fig. 5.19b are the fixed solution with GPS

data of that epoch. Obviously, fixed solutions can be determined with most F-ISB samples and the biases in the

fixed solutions are only several millimeters.

From above, the GPS L1 and BDS B1 integration with inter-system model is similar to the case of GPS L1 and

GLONASS L1 integration. Due to different frequencies, both of them require an approximate ISB value and an

accurate F-ISB value to recover the integer nature of inter-system DD-ambiguities, but the approximate ISB

vlaues are less important in the GPS L1 and BDS B1 integration becasue the wavelength difference is much

smaller and also becasue the BDS satellites are much fewer than GPS satellites in the experiments.

Fig. 5.15 Satellite sky plots of BDS for station TLSE on DOY 001 of 2015

Fig. 5.16 Numbers of GPS and BDS satellites for baseline TLSG-TLSE on DOY 001 of 2015, with elevation

mask of 10 degrees

Fig. 5.17 Three-dimensional RATIO distribution for baseline TLSG-TLSE on DOY 001 of 2015 for all epochs

(a) and averages over epochs (b)

Fig. 5.18 Impact of ISB biases in the baseline fixed solutions of GPS L1 and BDS B1 integration with respect to

the GPS L1 baseline fixed solution. The data are from baseline TLSG-TLSE at epoch 1703. The F-ISB sampling

interval is 0.02 m

Fig. 5.19 Averages of RATIO values along epoch time axis (a), as well as the impact of ISB biases on the

baseline fixed solutions with GPS L1 and BDS B1 integration compared to the values from GPS L1 solutions (b)

for baseline TLSG-TLSE. Result shown in (b) employs data at epoch 1703. The F-ISB sampling interval is

0.001 m

5.3.4 Half-Cycle Problem and Cluster Analysis Method

An accurate F-ISB is employed to recover the integer nature of inter-system DD-ambiguities for integration of

both the same frequencies and different frequencies. In the ISB estimation with the particle filter approach that

will be described in section 5.4, the initial interval is set to [-0.5, 0.5] cycles and all the samples over this interval

will be introduced as known F-ISB values to get the corresponding RATIO as their quality index. In the case that

the true value of an F-ISB is very close to 0.5 cycles, particles around -0.5 cycles and 0.5 cycles are of very

similar quality. However, current particle filter cannot handle such problem properly, as the particles are split

into two groups and cannot converge during the filtering, and thus cannot provide a precise ISB estimate. This is

referred to as the F-ISB half-cycle problem in this study. Here is an example of such problem in the GPS L1 and

Galileo E1 integration with the data from baseline GOP6-GOP7. The data are collected on DOY 001 of 2015 and

the receiver types and the antenna types can be found in Table 5.2. The RATIO values at each epoch form two

peaks at -0.095 m / 0.5 cycles and 0.095 m / 0.5 cycles, as shown in Fig. 5.20. Since no a priori information

about the expected ISB value is available, it is not possible to precisely narrow the searching interval to avoid

this problem.

Fig. 5.20 Three-dimensional RATIO map for baseline GOP6-GOP7. The F-ISB is 0.095m which is very close to

0.5 cycles, and there are two RATIO peaks within the initial interval [-0.1, 0.1] m

This problem was also pointed out while estimating the uncalibrated phase delay for the integer ambiguity

resolution of PPP (Ge et al. 2008). In this study, an approach based on cluster analysis in data mining is proposed

to classify all particles into clusters. During the filtering, the distance between the centroids of clusters equals to

one cycle and therefore the clusters can be shifted together to a single cluster. The procedure of the cluster

analysis method can be found in (Tan et al. 2006). The centroids need to be calculated more than one time in the

K-means algorithm, which is one of the traditional but widely used clustering algorithm.

The computation procedure can be refined as follows. Firstly the two particles with the largest distance are

selected as the first point of each cluster. Then, all particles are sorted to the closest cluster, and the centroids of

clusters are calculated. If the distance between two centroids is close to one wavelength, the particles in one

cluster will be shifted to another cluster by shifting one cycle. This procedure is carried out just after the update

step in the approach of section 5.4.

5.3.5 Conclusions

For the integration of GPS L1 and Galileo E1, an F-ISB value corresponding to a local maximum RATIO leads

to an accurate solution. For GPS L1 and GLONASS L1 integration, an accurate F-ISB value corresponding to

local maximum RATIO and near approximate ISB is preferred because RATIO values are larger there and fixed

solutions are closer to the single system solutions. GPS L1 and Galileo E1 integration is in the same situation as

GPS L1 and GLONASS L1 integration. The F-ISB value in multi-GNSS integration can be estimated according

to the RATIO distribution.

Therefore, in order to get accurate baseline solutions and relatively large RATIO, the integration with

inter-system model of the same frequency needs only an accurate F-ISB value corresponding to local maximum

RATIO, while the integration with different frequencies needs an approximate ISB value and an accurate F-ISB

value corresponding to local maximum RATIO. In section 5.4, the approach based on particle filter will be

presented to estimate precisely the accurate F-ISB values.

Besides, due to the periodic characteristic of ISB, the F-ISB value can be obtained with more than one local

maximum RATIO value, which can lead to the problem of the convergence of the filtering and therefore should

be considered in the estimation procedure.

5.4 Procedure for F-ISB Online Estimation

The prediction model for F-ISB is set up as

∆𝜇𝑎𝑏𝑘

𝑖=∆𝜇𝑎𝑏𝑘−1

𝑖+𝜖∆𝜇𝑎𝑏 (5.11)

where 𝜖∆𝜇𝑎𝑏 is assumed to be white noise. Model (5.11) is the same as (4.4) for the IFB rate estimation in

GLONASS data processing.

Although the PDF of carrier phase measurements in GNSS data processing does not provide information about

the integer ambiguity directly, the quality of the integer ambiguity candidates can give judgements of the pre-

defined ISB values. With the accurate F-ISB, the integer nature of the inter-system DD-ambiguities is supposed

to be recovered. Therefore, the closer the pre-defined F-ISB value is to the true value, the less bias in the

corresponding inter-system ambiguities and consequently the larger RATIO values. Similarly to the IFB rate

estimation in section 4.3, the normalized RATIO by (4.6) is employed to update the particles’ weights. The

approach to estimate the F-ISB is:

Step 1: Process the phase and code pseudorange measurements based on the models from (5.4) to (5.10).

Calculate the NEQ (2.24).

Step 2: For the first epoch, initial particles are obtained by sampling over initial interval [-0.5, 0.5] cycles.

Then the particles are assigned a weight 1/N. For other epochs, the particles have been prepared in the

previous epoch.

Step 3: For each particle, the model (2.26) from (2.24) is used to estimate the float ambiguities and the

associated VC matrix. The lambda method is then employed to determine the integer ambiguity

candidates. Next, the RATIO value is calculated.

Step 4: Normalize the RATIO values by (4.6). Update the weights with the normalized RATIOs. Calculate

the estimated fractional F-ISB and variance if needed.

Step 5: Resample the particles if (3.30) is satisfied. Then predict the particles for the next epoch according to

prediction model (5.11).

Step 6: Repeat steps 1-5 for the next epoch.

5.5 Results and Analysis

The experiments are divided into three parts. The first part is to show the accuracy and the short time variation of

the estimated F-ISB with the new approach; second part is to check the effects of the cluster analysis method for

the data with half-cycle problem; the third part provides a basic survey of the computation time.

5.5.1 F-ISB Estimate Results

F-ISB between GPS L1 and Galileo E1

The data of baseline TLSG-TLSE on DOY 001 of 2015 is processed via the new approach in section 5.4 with the

initial ISB samples over interval [-0.1, 0.1] m which is around one wavelength wide. The STD of the state noise

in (5.11) is set to 1 mm. The estimated F-ISBs of all epochs are presented in Fig. 5.21 with their corresponding

STDs. As only three Galileo satellites are observed during the observing period, there is a time span when no

Fig. 5.21 Estimated F-ISB and the three times STD for integration of GPS L1 and Galileo E1 for the baseline

TLSG-TLSE on the DOY 001of 2015

Fig. 5.22 Estimated F-ISB and the three times STD for integration of GPS L1 and GLONASS L1 for baseline

KOSG-KOS1 on DOY 048 of 2014, with approximate ISB -16.2764 m.

Fig. 5.23 Estimated F-ISB and the three times STD for integration of GPS L1 and GLONASS L1 for baseline

TLSG-TLSE on DOY 001 of 2015, with approximate ISB -4.943 m.

Fig. 5.24 Estimated F-ISB and the three times STD for integration of GPS L1 and BDS B1 for baseline TLSG-

TLSE on DOY 001 of 2015, with approximate ISB -0.3007 m

Galileo satellites are available and thus no estimated F-ISB values. In view of the number of satellites shown in

Fig. 5.2, it is clear that the STD of the estimated ISB decreases along with the increased number of Galileo

satellites. The average value of the estimated F-ISB values is 0.040 m.

The convergence time is about 8 minutes, i.e. 16 epochs in this experiment, with only one Galileo satellite at the

beginning of the day. The computation time for each epoch is around one second and will be shown in section

5.5.3.

F-ISB between GPS L1 and GLONASS L1

The GPS L1 and GLONASS L1 integration with inter-system models requires the approximate ISB value and an

accurate F-ISB value as described in section 5.4. In this experiment, the F-ISB is estimated for the data of 24

hours for baseline KOSG-KOS1 on DOY 048 of 2014 with the approximate ISB of -16.2764 m calculated with

code pseudorange observations. The STD of the state noise is still set to 1 mm and the estimated results are

presented in Fig. 5.22. Obviously, the F-ISB is pretty stable within the whole day, but with some fluctuations,

which are probably caused by error sources in GNSS observations, such as the receiver offset errors and the

atmospheric delays. The average of the estimated F-ISB values is -0.0090 m, leading to the final ISB value

of -16.2854 m.

In addition, the F-ISB in the GPS and GLONASS integration with inter-system models for baseline TLSG-TLSE

on DOY 001 of 2015 is also calculated with approximate ISB set to -4.943 m. The results are drawn in Fig. 5.23.

The average of the estimated F-ISB values is -0.0191 m, leading to the final ISB value of -4.9621 m.

F-ISB between GPS L1 and BDS B1

The data used here are still from baseline TLSG-TLSE collected on DOY 001 of 2015. The state noise in (5.11)

is set to 1 mm. An approximate ISB value and an accurate F-ISB value are needed due to the frequency

difference, which is the same as GPS L1 and GLONASS L1 integration. With an approximate ISB value

of -0.3007 m calculated with code pseudorange observations, the estimated F-ISB results are presented in Fig.

5.24. The average value of the estimated F-ISB is -0.0359 m, leading to the final ISB value -0.3366 m.

5.5.2 Performance of the Solution for the Half-Cycle Problem

The data with the half-cycle problem as shown in Fig. 5.20 is processed with and without the cluster analysis

procedure described in section 5.3.4. The estimated F-ISBs and their STDs without cluster analysis are shown in

Fig. 5.25a and b, respectively. The results with cluster analysis are shown in Fig. 5.25 c and d, respectively. The

cluster-analysis procedure can detect the clusters and unify the particles, which solves the half-cycle problem

very well. It needs very little computation time which can be ignored compared to the computation time of the

estimation procedure in section 5.4.

Fig. 5.25 Estimated F-ISB (a) and the corresponding STD (b) without the cluster analysis and results (c), (d) with

cluster analysis for GPS L1 and Galileo E1 integration for baseline GOP6-GOP7 on DOY 001 of 2015

5.5.3 Computational Efficiency

For the new approach described in section 5.4, the ambiguities are fixed many times for each epoch. As

computation time is also very important for online application, it is investigated in this subsection. The employed

computer is a PC equipped with 2.8 GHz i5 CPU and 4 GB memory card. In this section, the integrations of GPS

L1 and GLONASS L1, GPS L1 and Galileo E1 are taken as examples.

For each case, 200 particles are used, which indicates that the procedure of the data processing after forming

NEQ is computed 200 times. The left panel of Fig. 5.26 shows the computation time of GPS L1 and GLONASS

L1 integration for baseline KOSG-KOS1 on DOY 048 of 2014, as well as the satellite numbers. The

computation time of GPS L1 and Galileo E1 integration for baseline TLSG-TLSE on DOY 001 of 2015 is

presented in the right panel of Fig. 5.26 together with the satellite numbers. The calculation time in the right

panel is slightly shorter than that of the left one, which is due to the fewer satellites. As shown in Fig. 5.26, most

of the epochs can be finished in around 1 second. There are three epochs that take more than 2.0 seconds, which

is probably caused by output procedure as the long computation times appear at different epochs in repeated

experiments. Shortening the computation time in Fig. 5.26 by optimizing the program is possible.

Fig. 5.26 Computation time of Particle filter for KOSG-KOS1 GPS and GLONASS L1 integration (left).

Computation time of Particle filter for TLSG-TLSE GPS L1 and Galileo E1 integration (right)

5.6 Analysis of F-ISB Characteristics in Multi-GNSS Integration

The characteristics of the biases decide how the bias estimation approach can be used. If they are the same all the

time, then the biases can be estimated once at the beginning and be utilized in the following data processing. But

if they change occasionally, the estimate approach can be employed from time to time to check whether the

change has happened. If they are different from epoch to epoch, the approach has to be used in real-time to get

accurate value for every epoch. Therefore, data over long time are employed here to show the behaviour of the

F-ISB so that the approach can be utilized properly.

Long-time F-ISB Characteristics of GPS L1 and Galileo E1

In order to investigate the temporal stability of ISBs between GPS L1 and Galileo E1, almost all short baselines

in MGEX are selected and long-term data of these baselines are processed using the new approach from section

5.4 to obtain the F-ISBs. We selected firstly three days, DOY 001, 120 and 181 of 2015, to give a snapshot of the

F-ISB values. If there is no data on the specified date, data of the nearest day within one week is taken. In case of

significant changes among the three daily ISB estimates, more data of the related baselines will be processed for

further investigation.

There are in total 18 baselines associated with 27 stations. The F-ISB results as well as the baseline lengths are

presented in Table 5.1. As the F-ISB value may depend on the receiver types and firmware (Odijk and Teunissen

2013a), they are given in Table 5.2 for all the receivers on the three days.

From Table 5.1, most of the baselines have non-zero F-ISB because of using different type of receivers but with

two exceptions. Baselines UNX2-UNX3 and WTZ3-WTZZ have zero F-ISB, though the receivers for the former

one are from different manufacturers, JAVAD and SEPT, and for the latter one the same type of receivers are

used but with different firmware.

Table 5.1 F-ISB estimation results of short baselines in MGEX

The Baseline

name

Length (m)

F-ISB (m)

DOY 001

DOY 120

DOY 181

DUND-OUS2

6888

-0.038

-0.042

GOP6-GOP7

-0.094

-0.092

HARB-HRAG

2249

-0.040

-0.087

-0.086

KIR8-KIRU

4469

-0.039

-0.040

OHI2-OHI3

0.095

0.094

RGDG-RIO2

-0.040

-0.087

SIN0-SIN1

0.040

TLSE-TLSG

1266

-0.040

-0.087

UNB3-UNBD

-0.035

UNB3-UNBN

0.055

UNBD-UNBN

0.090

-0.100

0.090

UNX2-UNX3

0.000

WTZ3-WTZR

-0.095

0.095

WTZ3-WTZZ

0.000

WTZR-WTZZ

0.095

-0.095

ZIM2-ZIM3

0.000

ZIM2-ZIMJ

-0.039

-0.036

ZIM3-ZIMJ

-0.038

-0.039

-0.036

Comparing F-ISBs of the different days, the F-ISBs change hardly along with time for all but three baselines.

The three baselines HARB-HRAG, RGDG-RIO2 and TLSE-TLSG have a jump of -47 mm from DOY 001 of

2015 to DOY 120 of 2015. However, there were no changes in either the hardware or the firmware. Because of

the periodic characteristic, if the difference between two F-ISB values is near 0.19 m which is the wavelength,

the two F-ISB values will be considered to be the same.

It should be noticed that baselines HARB-HRAG and ZIM3-ZIMJ are equipped with receivers from the same

two manufacturers and have F-ISB value difference of only 2 mm on DOY 001 of 2015, but the F-ISB value

difference is 5 cm on DOY 181 of 2015 with unchanged receivers.

For further investigation, data over longer time for the three baselines are processed. The estimated F-ISB time

series are plotted in Fig. 5.27. The results show that the F-ISB is very stable except for the jumps. Then the

estimated F-ISB values with data around the change point are presented in Fig. 5.28, which shows that all

changes happened between DOY 013 of 2015 and DOY 014 of 2015.

From the above numerical study on long-term ISB characteristics, in general, F-ISB between GPS L1 and

Galileo E1 may not be zero if different types of receivers are employed but they are very stable in time and can

be estimated. However, there are unreasonable rapid changes which cannot be explained and need further

investigation.

Long-time F-ISB Characteristics of GPS L1 and GLONASS L1

The F-ISBs between systems of different frequencies encounter large jumps more frequently for baselines

equipped with receivers of both the same type and different types. This is also true for GPS L1 and BDS B1

integration. The results of the estimated F-ISB for baseline KOSG-KOS1 are presented here.

Table 5.2 Receiver types and firmware series for each station in the short baselines in MGEX

Station

name

Receiver type

Receiver firmware

DOY 001

DOY 120

DOY 181

DUND

Trimble NetR9

4.81

GOP6

LEICA GRX1200+GNSS

8.71/6.112

GOP7

JAVAD TRE_G3TH DELTA

3.5.1

HARB

TRIMBLE NETR9

4.85

5.01

HRAG

JAVAD TRE_G2T DELTA

3.6.1

KIR8

TRIMBLE NETR9

4.85

5.01

KIRU

SEPT POLARX4

2.5.2-esa3

OHI2

JAVAD TRE_G3TH DELTA

3.5.3

3.6.1

OHI3

LEICA GR25

3.11.1639/6.403

OUS2

JAVAD TRE_G3TH DELTA

3.5.7

RGDG

TRIMBLE NETR9

4.85

5.01

RIO2

JAVAD TRE_G3TH DELTA

3.4.7

SIN0

JAVAD TRE_G3TH DELTA

3.4.7

3.6.1

SIN1

TRIMBLE NETR9

4.80

TLSE

TRIMBLE NETR9

4.85

5.01

TLSG

SEPT POLARX4TR

2.5.2

2.9.0

UNB3

TRIMBLE NETR9

4.85

5.01

UNBD

JAVAD TRE_G2T DELTA

3.6.1

UNBN

NOV OEM6

OEM060510RN0

000

OEM060510RN0

000

OEM060510RN0

000

UNX2

JAVAD TRE_G3TH DELTA

3.4.7

3.6.1

UNX3

SEPT ASTERX3

2.3.4

WTZ3

JAVAD TRE_G3TH DELTA

3.4.14

3.6.1

WTZR

LEICA GR25

3.11.1639/6.403

WTZZ

JAVAD TRE_G3TH DELTA

3.6.0

3.6.1

3.6.2

ZIM2

TRIMBLE

NETR9 4.85

NETR9 5.01

ZIM3

TRIMBLE NETR9

4.85

4.93

5.01

ZIMJ

JAVAD TRE_G3TH DELTA

3.4.9

Data of baseline KOSG-KOS1 are employed firstly. The F-ISBs for 22 days within 600 days are calculated and

presented in the top panel of Fig. 5.29 with the approximate ISB set to -16.2764 m. The reference day, which is

day zero in Fig. 5.29, is DOY 167 of 2013. Even though there is a chance that several days are with similar

values, the F-ISB is clearly not the same during long time section and the variation is large. From the F-ISB

values of continuous 40 days, which are presented in the bottom panel of Fig. 5.29, jumps are pretty large and

can be clearly observed.

5.7 Summary

Integrations with both signals of the same frequency, GPS L1 and Galileo E1, and signals of different

frequencies, GPS L1 and GLONASS L1, as well as GPS L1 and BDS B1, are investigated and analysed. The

former case needs only an accurate F-ISB value, while the latter case requires both the approximate ISB value

and the accurate F-ISB value. The approximate ISB value can be regarded as equal to the ISB value for code

pseudorange observations, and therefore can be estimated from SPP.

With a pre-defined F-ISB value, the model without unknown F-ISB parameter is employed to fix both intra- and

inter-system DD-ambiguities. Usually, the accuracy of the pre-defined F-ISB values decides the performance of

the ambiguity fixing, and so the magnitude of RATIO values as well, leading to that the pre-defined F-ISB

values can be judged by the corresponding RATIO values. Therefore, a new approach based on particle filter is

proposed to estimate the F-ISB, which can estimate and track the F-ISB values accurately without any a priori

value. This approach needs around 1 second to finish the computation with 200 particles on a PC, and it is

possible to reduce the computation time by optimizing the program.

Since the particles can achieve similarly large weights with different F-ISB values due to the periodic

characteristics of the ISB, the particles can be divided into two or more groups, which hinder the estimation of

F-ISB and is denoted as half-cycle problem in this thesis. Thus, the cluster analysis method is introduced to solve

this problem. Once the two groups of particles are detected, they are merged together again. The experiments

show that the new approach can estimate the F-ISB precisely.

Finally, the F-ISB characteristics during a long period of time are investigated in this chapter, which shows that

the F-ISB values in the experiments are stable over time, but with very large jumps. For the short baselines of

MGEX, the F-ISB values of three baselines have jumps at the same time but no other changes are observed in

the integration of GPS L1 and Galileo E1. The changed F-ISB leads to the fact that even though two baselines

are equipped with the same receivers, they may have largely different F-ISB values. For the integration of

different frequencies, the jumps are more likely to happen for receivers of both the same type and different types.

Hence, it is necessary to monitor the F-ISB values from time to time.

Fig. 5.27 F-ISB of GPS L1 and Galileo E1 integration within about one year for HARB-HRAG (a),

RGDG-RIO2 (b) and TLSE-TLSG (c)

Fig. 5.28 F-ISB of GPS L1 and Galileo E1 integration within 5 days around the jump epoch for HARB-HRAG

(a), RGDG-RIO2 (b) and TLSE-TLSG (c)

Fig. 5.29 F-ISB of GPS L1 and GLONASS L1 integration with approximate ISB set to -16.2764m within about

two years after DOY 167 of 2013 (top), as well as short-term F-ISB values within 40 days with available data

(bottom) for baseline KOSG-KOS1

6 Two-Dimensional Approach

This chapter extends the one-dimensional particle filter approaches into multi-dimensional approaches so that the

F-ISB can be estimated even with fewer satellites from each constellation in the integration of more than two

systems. The two-dimensional approach for F-ISB estimation is taken as an example. At the beginning, the

reason for using multi-dimensional method is discussed in section 6.1. Then, the RATIO distribution with two F-

ISB parameters is investigated in section 6.2 and the two-dimensional particle filter approach is described in

section 6.3. Finally, the estimated F-ISB results in the test with full constellations and also with fewer satellites

from each constellation are shown in section 6.4.

6.1 Motivation for Multi-dimensional Approach

As aforementioned, the existence of IFB/ISB lays obstacles on integer ambiguity resolution and the traditional

methods where the bias is estimated along with the float ambiguities have lower accuracies and converge slowly.

This is because the estimated bias can be precisely determined only when the ambiguities converge well, such as

that the DD-ambiguities are successfully fixed to integers when sufficient satellites from each constellation are

available. However, in many applications of multi-GNSS integration, especially in urban areas, where signals

could be easily blocked or interrupted, only a few satellites of each system can be observed. In this case, it is

obvious that the bias cannot be estimated by traditional methods.

In fact, when the bias parameter is known, all DD-ambiguities in the equations with bias parameters can be fixed,

so that the fixed performance is significantly improved. This is well utilized by the proposed particle filter

approach, where the pre-defined bias values are introduced into the data processing as known bias values and are

judged by the performance of the ambiguity fixing. An accurate bias value can recover the integer nature of the

DD-ambiguity in DD-equations with bias parameters and therefore a better performance of the ambiguity fixing

can be achieved. For the F-ISB estimation, in the case of two systems one more integer ambiguity parameter is

included in the ambiguity resolution with the pre-defined ISB values near the true value. Therefore, the benefits

from inter-system model with integer ambiguity are well utilized even in the estimation of ISB itself. However,

if the satellites from each system are much fewer, not only the traditional approaches, but also the particle filter

approach for two systems may not be able to estimate F-ISB precisely.

For example, if two satellites from each system are available in the case of three systems, only three

DD-ambiguities with integer nature can be included without known F-ISB. In this case, it is almost impossible to

estimate the two F-ISB values accurately with traditional methods. Even with the particle filter for two systems,

only one ambiguity is added, i.e. there are four integer DD-ambiguities. In this case, the ambiguity fixing is

possible but not reliable. Therefore, the method is extended to multi-dimensional case. For the above-mentioned

example, the two F-ISB are estimated parallel with five DD-ambiguities in the way of particle filter using

RATIO as the quality index of the given F-ISB.

To show the benefits of multi-dimensional particle filter approach, the example of the two-dimensional case of

F-ISB estimation is presented in the following sections with the integration of three frequency bands, GPS L1,

BDS B1 and Galileo E1.

6.2 Relationship between RATIO and Two F-ISB parameters

The relationship of fixing RATIO and the two F-ISB parameters in the integration of GPS L1, BDS B1 and

Galileo E1 is investigated in this subsection, where two independent F-ISB parameters among the three carrier

frequency bands need to be estimated.

Firstly the initial interval [-0.2, 0.2] m, whose width equals about two times the wavelength, is sampled 40 times

with the sampling interval of 0.01 m. These samples are set as F-ISB values for GPS L1 and Galileo E1

integration directly, and set as F-ISB for GPS L1 and BDS B1 integration after plus -0.3007 m which is the

approximate ISB. Consequently, there are totally 1600 combinations. For each pair of F-ISB values, the RATIO

value is calculated for each epoch.

The result of epoch 476 at 3:58:00 AM is taken as an example. Totally 7 GPS, 2 Galileo and 3 BDS satellites are

observed with elevation mask of 10 degrees. The RATIO distribution is plotted in the left panel of Fig. 6.1. Four

local maximum RATIO values corresponding to four pairs of F-ISB values can be observed. Any of the four

pairs can be used to remove the ISB effects in the inter-system models due to the periodic characteristic. A

further calculation shows that the locations of the four maximum values are stable over time, which indicates

that both the F-ISBs are stable.

When the satellites are fewer, the local maximum values in the RATIO distribution may be located at different

places due to unreliable ambiguity fixing, but these values can still be large. For example, with 3 GPS, 2 Galileo

and 3 BDS satellites at epoch 476, the RATIO distribution is presented in the right panel. There are still four

large local maximum RATIO values but obviously they are at different locations compared to the maximum

values in the left panel of Fig. 6.1. Therefore, it is not reliable to simply select one of the four pairs of F-ISB

values corresponding to the local maximum RATIOs to correct the biases. Thus, a two-dimensional particle filter

approach will be introduced to estimate the two F-ISBs in the next subsection.

Fig. 6.1 RATIO distribution with two F-ISBs of GPS L1 and Galileo E1, GPS L1 and BDS B1 for epoch 476 for

baseline TLSG-TLSE, with all observed satellites (left) and later with only 3 GPS, 2 Galileo and 2 BDS satellites

(right)

6.3 Two-dimensional Particle Filter

The prediction models of the two F-ISB variables can be expressed as

𝜇𝑎𝑏𝑘

𝑖,𝑚=𝜇𝑎𝑏𝑘−1

𝑖,𝑚 +𝜖𝜇𝑎𝑏

𝑚, (6.1)

where 𝜖𝜇𝑎𝑏

𝑚 is assumed to be white noise; 𝑖 = 0, 1, …, N is the particle number; 𝑚 refers to GE or GB, which

refer to between GPS L1 and Galileo E1, between GPS L1 and BDS B1, respectively.

After the integer ambiguity candidates are determined by LAMBDA method, the PDF for the ambiguities fixed

to the correct integers given two F-ISB parameters can be expressed by

𝑝(𝒃

𝑘|(𝜇𝑘

𝐺𝐸,𝜇𝑘

𝐺𝐵)𝑖)= 𝑅𝐴𝑇𝐼𝑂𝑖

∑𝑅𝐴𝑇𝐼𝑂𝑖

𝑁

𝑖=1 , (6.2)

The remaining part of the procedure is similar to the one described in section 5.4. The cluster analysis method

described in section 5.3.4 is also employed.

In the one-dimensional particle filter described in section 5.4, the number of particles is set to 200. This leads to

40,000 particles in the two-dimensional approach, which is too large and leads to high computation burden. Here

the number of total particles is still set to 200, although the number of particles for each dimension is much

smaller.

The procedure for the two-dimensional estimation method is:

Step 1: Process the phase and code pseudorange measurements according to the models (5.4), (5.5), (5.7), (5.8)

and (5.10) to get the NEQ (3.24).

Step 2: For the first epoch, initial particles are obtained by sampling randomly over interval with width equal

one wavelength. As each particle has two F-ISB values, four hundred samples are generated to

compose two hundred sample pairs i.e. particles. Then all the particles are assigned the weight 1/N.

For an other epoch, the particles have been prepared in the previous epoch.

Step 3: For each particle, the equation (2.26) is used to calculate the float ambiguity and the associated VC

matrix. The LAMBDA method is then employed to obtain the integer ambiguity candidates and then

the corresponding RATIO values are calculated.

Step 4: Normalize the RATIO values by (6.2). Update the weights with the normalized RATIOs. Calculate

the estimated F-ISBs and their variances if needed.

Step 5: If the STDs of the estimated F-ISB is larger than a threshold, use the cluster analysis method to judge

whether the particles are divided into more than one group, and shift them into one group if yes.

Step 6: Resample the particles if (3.30) is satisfied. Then predict the particles for the next epoch according to

the prediction model (6.1).

Step 7: Repeat steps 1-6 for the next epoch.

6.4 Experiments with Two Dimensional Approach

The data for TLSG-TLSE on DOY 001 of 2015 are employed in this section as well. The F-ISBs for the

integration of GPS L1 and BDS B1, GPS L1 and Galileo E1 are estimated simultaneously. Firstly, the estimation

is implemented with all observed satellites. Then, it is carried out with few satellites from each system to test the

performance under severe conditions.

F-ISB Estimation with All Observed Satellites

In the F-ISB estimation, the STD of the state noise in (6.1) is set to 0.003 m. Data from 4:00:00 AM to

8:00:00 AM with epoch interval of 30 s are employed. The convergence process is presented with each panel for

one epoch in Fig. 6.2 a-g, where the particles move gradually to the area with larger RATIO values and the

estimated F-ISBs move to the true values. In the background of each panel in Fig. 6.2 a-g, four relatively larger

RATIO values can be observed because the F-ISB values on each direction are over interval [-0.2, 0.2] m, i.e. is

two wavelengths long. During the filtering, the particles are not limited within [-0.2, 0.2] but are free to move.

The cluster analysis method in the procedure in section 6.3 solving the half-cycle problem is employed here to

guarantee that all particles will converge to only one of the four local maximum RATIO values. In the

convergence process shown in Fig. 6.2 a-g, the half-cycle problem is not encountered as the F-ISB is not so close

to half cycles.

The estimated results are plotted in Fig. 6.2h. The filtering converges at seventh epochs, where the STD is

smaller than the thresholds which are 2 cm for GPS L1 and BDS B1, 6 mm for GPS L1 and Galileo E1. The

thresholds are set to values which are a little larger than the STDs of the weighted particles after convergence.

The threshold for GPS L1 and Galileo E1 integration is set much smaller because GPS L1 and Galileo E1

integration leads to relative larger RATIO values and hence the F-ISB can be better distinguished by the RATIO

distribution i.e. the areas with larger RATIO values in Fig. 6.2 a-g are narrower in the direction of F-ISB for

L1-E2 axis. Besides, in this estimation, the state noise is set higher than the estimation in section 5.5, so that the

convergence time is shorter even with fewer particles because with larger noise the particles are more likely to

reach the convergence area earlier. However, a higher state noise also leads to larger STD of the weighted

particles even after convergence.

The estimated F-ISB values for the whole 4 hours are shown in Fig. 6.3. The F-ISB value for GPS L1 and BDS

B1 is -0.0377 m. Considering the approximate ISB value -0.3007 m which is estimated with code pseudorange

observations, the final ISB value is -0.3384 m, compared to -0.3366 m in section 5.5.1, while the F-ISB for GPS

L1 and Galileo E1 is 0.0398 m, compared to 0.0400 m in section 5.5.1.

Simulated Scenario with Fewer Satellites

To investigate the performance of the two-dimensional approach under severe conditions, an observation

scenario with two satellites from each system in view is defined. Thus, the total number of satellites is six. Due

to the incomplete constellations of Galileo and BDS, only data of around five hours include such six satellites on

DOY 001 of 2015. Due to the movement of the satellites, the PRN numbers of the six satellites change during

the five hours and are presented in the left panel of Fig. 6.4.

As the satellites are from three GNSS, two independent inter-system models can be formed with two F-ISB

parameters. The two F-ISB parameters are then estimated by the two-dimensional particle filter approach

described in section 6.3. The estimated results are presented in the right panel of Fig. 6.4, where it takes around

half an hour to converge. It can also be observed that the STDs of the estimated F-ISBs are different during the

five hours, which is probably due to different satellites and different satellite conditions, such as different

elevation angles.

Fig. 6.2 Convergence process with the two-dimensional particle filter approach for epochs from one to seven

(a-g) and the estimated F-ISB results for GPS L1 and BDS B1 combination (h, top), GPS L1 and Galileo E1

combination (h, bottom)

Fig. 6.3 Estimated F-ISB for GPS L1 and BDS B1 combination (top), as well as GPS L1 and Galileo E1

combination (bottom) for baseline TLSG-TLSE

After the IFBs converged, the baseline solutions are calculated. In the data processing, the SD-ambiguities

propagate from one epoch to the next if no cycle slips occur, while the DD-ambiguity fixing is carried out each

epoch to fix the DD-ambiguities. Since the estimated F-ISB are fixed as known values, both intra- and

inter-system DD-ambiguities can be fixed. For comparison, the same observations are also processed without

inter-system models where only intra-system DD-ambiguities are fixed. Thus, there are only three integer

DD-ambiguities, instead of five for fixing both intra- and inter-system DD-ambiguities.

The results of the aforementioned two strategies for data processing are shown in Fig. 6.5. With the approach of

fixing only the intra-system DD-ambiguities, almost no successful fixed solutions are available, whereas the

availability rate of the fixed solution for strategy fixing both intra- and inter-system DD-ambiguities is about

88.1% of all epochs within the five hours. The DD-ambiguities are fixed successfully for the first time at

23 minutes and for all epochs after 37.5 minutes.

Fig. 6.4 Satellite PRN numbers (left) and the estimated F-ISB results for GPS L1 and BDS B1 (right top), GPS

L1 and Galileo E1 (right bottom) with only two satellites from each system

Fig. 6.5 Positioning differences with respect to the GPS static solution for the strategy fixing only intra-system

DD-ambiguities and that fixing both intra- and inter-system DD-ambiguities, which are denoted as Only intra

DDAF (DD-Ambiguity Fixing) and Intra and inter DDAF in the figures, respectively

6.5 Summary

The multi-dimensional particle filter is investigated for the estimation of F-ISB parameters in multi-GNSS

integration. The F-ISB estimation with a two dimensional particle filter is taken as an example. The investigation

shows that the RATIO values are larger when the values of the F-ISB pair are closer to their true values.

Therefore, they can also be used to judge the pre-defined F-ISB values in the two-dimensional case. Thus, the

two-dimensional particle filter is employed to simultaneously estimate the two F-ISB parameters.

In the test, the two F-ISB values between GPS L1, Galileo E1 and BDS B1 are determined with the

two-dimensional particle filter approach. When only two satellites from each system are observed, the two F-ISB

values can still be estimated within 30 minutes. With the estimated F-ISB parameters, the approach fixing both

intra- and inter-system DD-ambiguities can be employed and the fixed solutions are available at 23 minutes for

the data of five hours, but the accurate fixed solutions cannot be determined with strategy fixing only intra-

system DD-ambiguity, which indicates the traditional methods for F-ISB estimation fail.

7 Application of the Phase IFB Rate and F-ISB Estimation for Precise

Positioning

This chapter investigates the applications of the IFB rate and F-ISB estimated by particle filter approaches for

precise positioning in terms of the integer ambiguity resolution performance and the positioning accuracy. Both

single-epoch data processing method and continuous kinematic processing method are employed. In the former

method, the current epoch is used without any information from the previous epochs, while in the later method

the ambiguities propagate from epoch to epoch if no cycle slips occur.

The empirical availability rate (EAR) defined in (Li et al. 2015) is employed to show the performance of the

integer ambiguity fixing. The EAR is calculated by

𝑓= 𝑁precise

𝑁total , (7.1)

where 𝑁precise refers to the number of epochs with fixed solution, while 𝑁total refers to the total number of

epochs.

7.1 GLONASS Data Processing with Estimated Phase IFB Rate

7.1.1 Single-Epoch Processing

The GLONASS data of the baseline KOSG-KOS1 is processed epoch by epoch independently with the IFB rate

value estimated in Section 4.4.1. The impact of the estimated IFB on positioning results is demonstrated by

comparing the single-epoch solution with integer ambiguity resolution at each epoch. As the IFB rate for this

baseline is far from zero, the ambiguity resolution cannot succeed if the IFB rate is unknown. Without fixed

integer ambiguities, the related solution is actually the single-epoch float solution.

Fig. 7.1 shows the position difference of the solution with respect to the static result in East, North and Up

directions with (blue) and without (red) IFB rate estimation procedure at the beginning. The availability rate with

estimated IFB rate is 97.9%. As soon as the ambiguities are fixed, the positioning STDs reach values of about 2

mm, 2 mm and 5 mm in East, North and Up directions, respectively. There are very few epochs with a very large

difference due to the failure of ambiguity fixing. For the single-epoch solution with unknown IFB rate, the

position accuracy is about several decimetres for each component. Obviously the new approach can fix the

integer ambiguities with single-epoch data and the positioning results are largely improved.

7.1.2 Kinematic Positioning with Continuous Ambiguity

The data of baseline KOSG-KOS1 are employed for the kinematic positioning with IFB rate online calibration

procedure. The processing started with IFB rate estimation and at the seventh epoch the estimated IFB rate

reaches the converged status with STD of 0.78 mm/FN. Therefore, the particle filter is stopped as described in

section 4.3. The remaining data are processed with the estimated IFB rate as known value and ambiguities are

fixed at each epoch. The results show that the ambiguity resolution has an availability rate of 98.6%. The

baseline components in the results are presented in Fig. 7.2, together with the float solution. The left panel of

Fig. 7.2 shows the whole time series of the position difference, while the right panel shows the result of the first

15 minutes where the impact of the integer ambiguity resolution is clearly visible. The position time series of the

GLONASS fixed solution is also compared with that of the GPS fixed solution in Fig. 7.3. The differences in

East, North and UP directions are 1.3 mm, 1.0 mm, 1.7 mm with STD of 2.3 mm, 2.8 mm, 5.5 mm, respectively.

This indicates that the GLONASS fixed solutions for real-time kinematic positioning with IFB rate estimated by

the new method has the same performance as GPS.

The integer solution largely shortens the convergence time of the baseline solution. In order to get statistics of

the convergence time, the data of the whole day for baseline KOSG-KOS1 are divided into 72 time sessions and

each session is 20 minutes long. The data of each session is then processed with ∆𝛾𝑎𝑏 estimation at the beginning.

The IFB results estimated by particle filter approach for 72 sessions are presented in Fig. 7.4, which reveals high

precision. The corresponding time series of the baseline components are presented in Fig. 7.5 to show the

convergence process of positioning results, while Fig. 7.6 presents the final horizontal position differences of the

72 sessions, where each of the positioning differences is represented by one dot. Obviously, the solutions

converge more quickly with ∆𝛾𝑎𝑏 estimated at the beginning. This indicates that with IFB rate estimated by

particle filter approach, the fixed solutions can be estimated quickly even without an a priori IFB rate value.

Fig. 7.1 Comparison of the GLONASS single-epoch solution with (blue) and without (red) IFB rate estimation

procedure at the beginning for baseline KOSG-KOS1

Fig. 7.2 Position differences with respect to the ground truth for GLONASS without (red) and with (blue) the

initial calibration of IFB rate in kinematic mode. The right panel is a snapshot of the first 30 epochs (15 minutes)

of the left one in order to show the process of the IFB rate estimation and the first ambiguity-fixing afterwards

Fig. 7.3 Comparison of GLONASS fixed solution and GPS fixed solution

Fig. 7.4 IFB rate estimates for all the 72 sessions for baseline KOSG-KOS1, as well as the corresponding STD

Fig. 7.5 Convergence processes of the GLONASS baseline solutions for all 72 sessions with and without IFB

rate estimation procedure for baseline KOSG-KOS1

7.2 Multi-GNSS Data Processing with Estimated F-ISB

7.2.1 Single-Epoch processing

The solutions from integration of only fixing the intra-system DD-ambiguities without inter-system observations,

as well as fixing both intra- and inter-system DD-ambiguities with estimated F-ISB values, are presented. Firstly,

the results of integrations between GPS L1 and Galileo E1, GPS L1 and BDS B1 are shown. Then, the three

frequency bands, GPS L1, Galileo E1 and BDS B1 are integrated together and results are discussed. Finally, the

integration between GPS L1 and GLONASS L1, which have full constellations, are computed with different

elevation masks to investigate the performance under severe conditions.

Fig. 7.6 Final horizontal positions of the GLONASS baseline solutions for all 72 sessions with and without IFB

rate estimation procedure for baseline KOSG-KOS1

Fig. 7.7 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS static solution (left),

as well as the RATIO values (right) with the strategy fixing only GPS L1 intra-system DD-ambiguities (green)

and the strategy fixing GPS L1 and Galileo E1 both intra- and inter-system DD-ambiguities (blue)

Fig. 7.8 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS static solution (left),

as well as the RATIO values (right) with the strategy fixing GPS L1 and Galileo E1 only intra-system

DD-ambiguities (red) and the strategy fixing GPS L1 and Galileo E1 both intra- and inter-system

DD-ambiguities (blue)

Baseline Solution of GPS L1 and Galileo E1 Integration

The single epoch positioning results using GPS L1 and Galileo E1 observations with estimated F-ISB of 0.040 m

in section 5.5.1 are shown in the left panel of Fig. 7.7. Even though the Galileo satellites are fewer and are not

observed all the time as shown in Fig. 5.2, the EAR of the whole day is improved from 75.5% to 81.2% by an

increase of 5.7%. The corresponding RATIOs are presented in the right panel of Fig. 7.7, from which it is clear

that larger RATIOs can be obtained when the Galileo satellites are included.

Later on, the strategy fixing only intra-system DD-ambiguities and the strategy fixing both intra- and

inter-system DD-ambiguities are implemented with both GPS L1 and Galileo E1. The results are presneted in

Fig. 7.8. The former strategy can improve the EAR to 77.08% by an increase of 1.6% compared with solutions

of GPS L1 only, but the EAR of the later strategy is even 4.1% higher than that of the former strategy.

Baseline Solution of GPS L1 and BDS B1 Integration

The single epoch positioning result for the same baseline using GPS L1 and BDS B1 observations are shown in

the left panel of Fig. 7.9, while the corresponding RATIO values are presented in the right panel of Fig. 7.9,

where again strategies of fixing only intra-system DD-ambiguities and fixing both intra- and inter-system

DD-ambiguities are employed. In the latter case, the F-ISB value is set to -0.3366 m as estimated in section 5.5.1.

The EAR of the ambiguity fixing of the former strategy is 81.5% and that for the latter one 84.8%, which is by

an increase of 3.3%; both are much higher than 75.5% for GPS only results shown in Fig. 7.7 (left).

Baseline Solution of GPS L1, Galileo E1 and BDS B1 Integration

The positioning results are calculated in single epochs using the following three processing strategies. In the first

strategy only GPS L1 observations are employed and of course only intra-system DD-ambiguities are available.

The second strategy includes observations of all three-system but only intra-system DD-ambiguities are fixed as

integer, while the third one fixes both intra-system and inter system DD-ambiguities. The results of these three

strategies are presented in Fig. 7.10. The EAR for the first strategy is 75.5%, which is improved to 82.0% by an

increase of 6.5% with the second strategy and to 86.2% by an increase of 10.7% with the third strategy.

Baseline Solutions of GPS L1 and GLONASS L1 Integration with Different Elevation Masks

Employing observations of both GPS L1 and GLONASS L1, the baseline KOSG-KOS1 is solved in

single-epoch mode and ambiguity fixing is carried out for both intra- and inter-system DD-ambiguities. The

results are shown in the left panel of Fig. 7.11, together with the results of the strategy fixing only intra-system

DD-ambiguities. It is clear that the results agree with each other very well. The elevation mask is 10 degrees and

there are usually more than 12 satellites as drawn in Fig. 5.7 which enables integer ambiguities to be fixed 100%

for both strategies.

In order to investigate how the performance changes with limited number of satellites by selecting different

elevation masks, the EAR against increasing elevation masks are shown in the right panel of Fig. 7.11. The

strategy fixing both intra- and inter-system DD-ambiguities is slightly better than that fixing only intra-system

DD-ambiguities. The GPS L1 baseline solution has a much lower EAR and the GLONASS L1 baseline solution

is even worse. When the elevation masks are higher, such as 40 degrees, some epochs can have a RATIO value

larger than 3, but the biases in the solution are more than 2 cm in horizontal direction or more than 3 cm in

vertical direction. In this case, the integer ambiguities are simply considered to be unsuccessfully resolved.

The same experiment procedure is also employed to solve baseline TLSG-TLSE with the data collected on DOY

001 of 2015. With elevation mask of 10 degrees, the EAR is improved from 99.1% to 99.8% by fixing also the

inter-system DD-ambiguities. The solutions are shown in the left panel of Fig. 7.12. After that, the EAR with

different elevation masks are also investigated and shown in the right panel of Fig. 7.12. Obviously, the EAR of

the strategy fixing both intra- and inter-system DD-ambiguities is also slightly higher than that of the strategy

fixing only intra-system DD-ambiguities, but both of them are much better than GPS only results.

7.2.2 Kinematic Positioning with Continuous Ambiguity

In this section, The first-fixing times in the integer ambiguity resolution with strategy fixing only intra-system

and strategy fixing both intra- and inter-system DD-ambiguities are compared by taking GPS L1 and Galileo E1

as an example. Later on, the scenarios with sky shelters, which are common in deformation monitoring, are

simulated to show the performance of the integration with estimated F-ISB, where the data of GPS L1 and

GLONASS L1 are selected as both systems have full constellations.

Fig. 7.9 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static solution

(left), as well as the RATIO values (right) with the strategy fixing GPS L1 and BDS B1 only intra-system DD-

ambiguities (red) and the strategy fixing GPS L1 and BDS B1 both intra- and inter-system DD-ambiguities (blue)

Fig. 7.10 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static solution,

employing three strategies, using GPS L1 only (green), using GPS L1, Galileo E1 and BDS B1 observations

fixing intra-system DD-ambiguities (red), and using the three constellations but fixing both intra- and inter-

system DD-ambiguities (blue)

First-fixing Time with/without Inter-system DD-Ambiguity Fixing

The inter-system DD-ambiguity has integer nature after applying F-ISB correction and then it can be fixed along

with intra-system DD-ambiguities. Therefore, an experiment is designed to demonstrate the advantage of fixing

inter-system DD-ambiguity in terms of the first-fixing time. If there is a large number of GPS satellites, the

fixing can be already achievable just with GPS data of a few epochs for such short baseline. In order to clearly

show the advantage of including the inter-system ambiguity, we choose a constellation with only five GPS

satellites and one Galileo satellite. Data from 1:30:00 to 8:30:00 UTC are selected during which Galileo E12 is

observed. The data are divided into 15 sessions with a length of 30 minutes. Two processing strategies, fixing

only intra-system DD-ambiguities and the strategy fixing both intra- and inter-system DD-ambiguities, are

carried out. Both of them include all six satellites in the data processing with estimated F-ISB. The first-fixing

times for all the sessions are plotted in Fig. 7.13.

The average first-fixing times without and with the inter-system ambiguities are 11.2 minutes and 5.1 minutes,

respectively. It is clear that the strategy fixing both intra- and inter-system DD-ambiguities needs only about half

of the observation time to get fixed solutions.

Fig. 7.11 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static solution

(left) and the EARs of integer ambiguity fixing with different elevation masks (right)

Fig. 7.12 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static solution

(left) and the EARs of integer ambiguity fixing with different elevation masks (right)

Fig. 7.13 First-fixing time with and without inter-system DD-ambiguity fixing for the processing scenarios of

GPS L1 and Galileo E1

Application to Simulated Environments with Few Observed Satellites

Different sheltering environments are simulated with the data of GPS L1 and GLONASS L1 for baseline

KOSG-KOS1, which are collected on DOY 048 of 2014. The sky plots of satellite constellations and numbers of

satellites have been presented in Fig. 5.6 and Fig. 5.7, respectively. In the data processing, the float

SD-ambiguities are propagated from epoch to epoch and the DD-ambiguity fixing is carried out for each epoch.

In the first simulation, for each system three satellites with elevation angles larger than 20 degrees are selected.

The satellite PRN/slot numbers of these satellites are shown in Fig. 7.14 and the biases of baseline solutions are

presented in Fig. 7.15. The strategy fixing both intra- and inter-system DD-ambiguities shows obviously much

better perfomance. The corresponding EAR is 93.6% which is 22.8% higher than that of 70.8% achieved by

fixing only intra-system DD-ambiguities.

In addition, two other obstruction simulations are simulated for the same baseline KOSG-KOS1. The second

simulation is that only the satellites in the west half sky (180 degrees < azimuth < 360 degrees) are available

while satellites in the east half sky (0 degrees < azimuth < 180 degrees) are obstructed. The third simulation is

opposite to the second one, only the satellites with azimuth belong to [0, 180] degrees are observed. The

elevation mask is set to 15 degrees. The final resutls are presented in Fig. 7.16 and Fig. 7.17, respectively. With

strategy fixing both intra- and inter-system DD-a mbiguities, the EAR is improved from 85.3% to 90% by an

increase of 4.7% for the scenario with satellites in the west half-sky and from 81.7% to 93% by an increase of

11.3% for the scenario with satellites in the east half-sky.

Fig. 7.14 Satellite PRN/slot numbers (left) and the differences of the baseline solutions with respect to the GPS

only static solution (right) for baseline KOSG-KOS1, with three GPS satellites and three GLONASS satellites.

Fig. 7.15 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static solution,

with three GPS satellites and three GLONASS satellites

Fig. 7.16 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static solution,

with satellites in the west half-sky

Fig. 7.17 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static solution,

with satellites in the east half-sky

7.3 Summary

The baseline solutions with the IFB rate estimated by particle filter approach in GLONASS data processing are

investigated. For ambiguity fixing in both single epoch strategy and continuous strategy, the IFB rate is obtained

within very short time and so are the GLONASS only fixed solution, which agree well with GPS only solutions.

Afterwards, in the experiments of multi-GNSS integration, the integrations between GPS L1 and Galileo E1,

GPS L1 and GLONASS L1, as well as GPS L1 and BDS B1 in single-epoch strategy are taken as examples to

show the improvements on ambiguity fixing by adding inter-system DD-ambiguities. In the kinematic

continuous data processing, the reduction of the first-fixing time with known ISB is demonstrated by GPS L1

and Galileo E1 integration. In the case of five GPS satellites and one Galileo satellite, the average first-fixing

time with strategy fixing also inter-system ambiguity is only half of that fixing only intra-system ambiguities.

Finally, the improvements on EAR with known F-ISB in the integration of GPS L1 and GLONASS L1 are

presented. The results show that the EAR increases reach up to 20% in the simulated sheltering cases, which are

very meaningful for practical usages where satellites can be easily blocked, such as precise positioning in urban

area, deformation monitoring of dams and side slides.

8 Conclusions and outlook

Global Navigation Satellite Systems (GNSS) has been a powerful tool of precise positioning for geodesy and

surveying engineering. The basis of the real-time GNSS precise positioning is the instantaneous integer

ambiguity resolution. However, some of the biases in carrier phase observations cannot be removed by

differencing between either stations or satellites, so the double-differenced ambiguities lose the integer nature

and cannot be fixed to integers. Therefore, these biases result in the failure of integer ambiguity resolution and

thus have to be estimated and removed in the data processing. Two of the biases are the carrier phase IFB in

GLONASS data processing and the carrier phase ISB in multi-GNSS integration.

The traditional methods estimate these biases along with the float ambiguities and then refine the estimated

values after successful ambiguity fixing. Due to the correlation between parameters and the low accuracy of the

float solutions, it is difficult to obtain the estimate of these biases with sufficient accuracy so that the ambiguities

fixing can succeed. Thus, the bias estimation requires observations over relatively long time or a priori values.

This study aimed to develop new approaches which can estimate these biases quickly and track them precisely

without a priori values.

The main contributions and conclusions of this research are presented in section 8.1 and the outlooks are given in

section 8.2.

8.1 Contributions and Conclusions

A new methodology for estimating the IFB/ISB biases in GNSS data processing is developed in this thesis.

Firstly, some bias samples (i.e. particles) are generated over the initial interval where the true bias value is

possibly located. The generated samples can remove the bias in the model when they are close enough to the bias

true value, but cannot remove the bias or can only partially remove it when they are far away from the true value.

As the remaining part of the bias in the GNSS model lowers the performance of the integer ambiguity resolution,

the RATIO value which indicates the reliability of the integer ambiguity resolution can be regarded as the quality

index of the given bias samples. Therefore, the RATIO values can be employed to judge the bias samples via

designed likelihood function, from which the true bias value can be estimated successfully. With the

IFB rate/F-ISB samples, the ambiguities in the models with IFB/ISB can be fixed to integers in the bias

estimation and thus their integer nature is well utilised. Therefore, this method can largely reduce the

convergence time and increase the reliability of the bias estimation.

This method can be successfully implemented with the particle filter approach, which is mainly composed by

three steps, update, resampling and prediction. Firstly, the update step renews the particle weights with the

likelihood function of RATIO. Secondly, the resampling step deletes the particles with too small weights and

duplicates the particles with large weights to ensure all particles contribute to the final estimated results. Finally,

the prediction step transmits particles from one epoch to the next one. The convergence of this approach can be

judged by STD of the weighted particles. The computation burden of each epoch mainly depends on the number

of particles because the RATIO value is calculated once for each particle.

The new method is first proposed to estimate the IFB rate in GLONASS data processing. The experiments with

practical data prove that RATIO can be regarded as a quality index of the given IFB rate value. In the

investigation of the relationship between RATIO and IFB rate, the maximum value of the RATIO corresponds to

the true F-IFB rate but with few exceptions when employing single epoch data. The exceptions are seen as

outliers and can be resisted by the filtering. With the particle filter approach, the IFB rate can be estimated

successfully with observations of few epochs and without the assistance of GPS or a priori values. After

convergence, the IFB rate can be precisely tracked online. The computation burden of the new approach is also

investigated. The result shows that the computation time at each epoch is less than 1 s on a PC with 200 particles.

Improvements of the developed approach for IFB rate estimation in two special cases are achieved. It is found

that setting the state noise in the prediction model of IFB rate to a low level introduces biases to the estimated

IFB rate, because the noise is needed to increase the diversity of the particles which is important to the precision

of the estimated IFB rate. Consequently, although more observations are available, the precision of the estimated

IFB rate cannot be further improved after the convergence is achieved. This problem can be solved by utilising

the regularization procedure in the approach. The precision of the estimated IFB can be continuously improved

with new observations. Besides, the adaptive method of selecting the number of particles is proposed to reduce

the computation burden in the tracking process. Although due to the large STD of the weighted particles at the

beginning, more than 200 particles are still needed to reach the convergence, after convergence the STD

becomes small and the number of particles decreases significantly, which leads to much shorter computation

time. In the test, the mean computation time is reduced from 0.84 s to 0.11 s for each epoch, but with almost the

same performance regarding the precision.

The new method is also applied to estimate the F-ISB in multi-GNSS integration. The F-ISB parameters can be

estimated online precisely by the particle filter approach via designed likelihood function of RATIO. Because

the RATIO values with F-ISB parameter have periodic characteristic, when the F-ISB value is close to half

cycles the particles are separated into multi-groups in the filtering process and hence the filtering cannot

converge. This problem can be well solved by the cluster analysis method which can automatically detect the

multi-group particles and then shift them together during the filtering. By investigating the F-ISB values of the

short baselines in MGEX project, it is found that the F-ISB values encounter large jumps for GPS L1 and Galileo

E1 integration, and the F-ISB can be significantly different even for the baselines which are equipped with

receivers of the same type. Although the particle filter approach needs relatively long computation time, the

computation for each epoch can be finished in around one second. Compared with the other F-ISB estimation

methods which fix only intra-system DD-ambiguities because the F-ISB is unknown, this approach brings the

advantage that the integer nature of the inter-system DD-ambiguities is well utilised in the estimation.

Investigation of inter-system models with different frequencies is carried out. Due to the wavelength difference,

the error of the phase ISB introduces bias into the SD-ambiguities. Hence, not only the F-ISB, but also ISB itself

affects the performance of the integer ambiguity resolution as well as the accuracy of the fixed solutions.

Therefore, an approximate ISB value and an accurate F-ISB value are needed to remove the ISB effects in the

inter-system model so that the integer nature of the inter-system DD-ambiguities can be recovered. The

approximate ISB can be regarded as equal to the ISB of the code pseudorange observations while the F-ISB can

be estimated with the particle filter approach. Although the performance of fixing the inter-system

DD-ambiguity with known ISB is affected by the biases in the initial SD-ambiguities due to the wavelength

difference, this strategy is still more likely to succeed in the test of GPS L1 and GLONASS L1 integration for

short baselines in the scenarios of few observed satellites.

Multi-dimensional estimation approach to estimate the biases is developed based on the one-dimensional particle

filter approach. With this approach, multi-biases can be estimated simultaneously. More integer ambiguities can

be included in the ambiguity fixing. Hence, the advantage of the one-dimensional approach is enlarged. In the

F-ISB estimation, the advantage is important when the observed satellites of each constellation are few. In this

case, the strategy fixing only intra-system DD-ambiguities fails and therefore the F-ISB parameters cannot be

estimated with traditional F-ISB estimation methods without accurate a priori F-ISB values. However, with the

multi-dimensional particle filter approach the F-ISB parameters can still be successfully estimated. The more

constellations the satellites are from, the more significant the advantage will be. In the experiments with

two-dimensional approach for F-ISB estimation of three constellations, the two F-ISB parameters are determined

successfully with only two satellites from each constellation.

Applications of the estimated IFB rate and the F-ISB in GNSS data processing are investigated. Because the IFB

rate can be estimated quickly without assistance of GPS or any a priori value by the new approach, the integer

ambiguity resolution in GLONASS data processing can succeed with observations of only a few epochs even

when the unknown IFB is far from zero. With the estimated IFB rate, the GLONASS fixed solutions for the short

baselines are as accurate as the GPS fixed solutions. For the applications of F-ISB in multi-GNSS integration,

compared with the strategy fixing only intra-system DD-ambiguities, the strategy fixing both intra- and inter-

system DD-ambiguities is more likely to succeed in the scenarios with few observed satellites of each system. In

the continuous kinematic data processing where ambiguities propagate from epoch to epoch, with totally six

GPS and Galileo satellites, the first-fixing time is reduced to nearly a half by fixing both intra- and inter-system

DD-ambiguities.

8.2 Outlook

The research in this thesis provides a new way for IFB rate and F-ISB estimation in GNSS data processing and

motivates further researches at least on the following aspects.

Firstly, the reason of the jumps of F-ISB values is unknown. The jumps of F-ISB, no matter for systems of the

same frequency or different frequencies, have been observed and their magnitudes are large enough to affect the

ambiguity fixing. If the reason can be found, it may be possible to remove or predict the jumps so that the

procedure of dealing with F-ISB can be simplified. Secondly, the performance of the inter-system

DD-ambiguities with different frequencies is affected by the initial SD-ambiguities which are calculated from

the code pseudorange observations, so there is a risk that the performance can be degraded if the code

pseudorange observations include large biases. Therefore, the condition to guarantee the improvements of

ambiguity fixing is important and should be investigated. Thirdly, developing the estimation approaches for

multi-GNSS with multiple frequency bands is another challenge. Six positioning satellite systems described in

section 2.1 will be fully operated in the following years. The constellations of these systems will include more

than 100 satellites with more than 20 frequencies. Many of these frequencies overlap with each other or have

slightly different wavelengths. Therefore, developing particle filter approaches to estimate the large number of

biases simultaneously is also worth being explored.

Bibliography

Alber C, Ware R, Rocken C, Braun J (2000) Obtaining single path phase delays from GPS double differences.

Geophysical Research Letters, 27(17):2661-2664

Al-Shaery A, Zhang S, Rizos C (2013) An enhanced calibration method of GLONASS inter-channel bias for

GNSS RTK. GPS solutions, 17(2):165-173

Al-Shaery A, Zhang S, Lim S, Rizos C (2012) A Comparative Study of Mathematical Modelling for

GPS/GLONASS Real-Time Kinematic (RTK). Proceedings of ION GNSS 2012, Nashville, TN,

September 2012, pp. 2231-2238

Arulampalam S, Maskell S, Gordon N, Clapp T (2002) A tutorial on particle filters for online

nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50(2):174-188

Banville S, Collins P, Lahaye F (2013) GLONASS ambiguity resolution of mixed receiver types without

external calibration. GPS solutions, 17(3):275-282

Bevis M, Businger S, Herring T, Rocken C, Anthes R, Ware R (1992) GPS meteorology: Remote sensing of

atmospheric water vapour using the Global Positioning System. Journal of Geophysical Research,

97(D14):15787-15801.

Banville S (2016) GLONASS ionosphere-free ambiguity resolution for precise point positioning. Journal of

Geodesy, 90(5):487-496

Blewitt G (1989) Carrier phase ambiguity resolution for the global positioning system applied to geodetic

baselines up to 2000 km. J Geophys Res, 94(B8):10187-10203

Candy J (2009) Bayesian signal processing: Classical, modern and particle filtering methods. John Wiley & Sons,

Inc., Hoboken, New Jersey, US

Chang X, Yang X, Zhou, T. (2005) MLAMBDA: a modified LAMBDA method for integer least-squares

estimation. Journal of Geodesy, 79(9):552-565

Chen Z (2003) Bayesian filtering: from Kalman filters to particle filters, and beyond. Working paper, McMaster

University, Canada

Closas P, Fernández-Prades C (2011) Particle filtering with adaptive number of particles. IEEE Aerospace

Conference, pp. 1-7

Cocard M, Geiger A (1992) Systematic search for all possible widelanes. Proceedings 6th International Geodetic

Symposium on Satellite Positioning, Columbus, Ohio, 1992, pp. 312-318

Counselman C, Gourevitch S (1981) Miniature interferometer terminals for earth surveying: ambiguity and

multipath with Global Positioning System. IEEE Transactions on Geoscience and Remote Sensing,

GE-19(4):244-252

Dach R, Brockmann E, Schaer S, Beutler G, Meindl M, Prange L, Ostini L (2009) GNSS processing at CODE:

status report. Journal of Geodesy, 83(3):353-365

Dach R, Schaer S, Hugentobler U, Schildknecht T, Gaede A (2006) Combined multi-system GNSS analysis for

time and frequency transfer. Proceedings of 20th European Frequency and Time Forum, Braunschweig,

Germany, March 2006, pp. 530-537

Dach R, Schaer S, Lutz S, Meindl M, Beutler G (2010) Combining the observations from different GNSS.

Proceedings of EUREF 2010 Symposium, Gävle, Sweden, June 2010, pp. 2-5

Dai L, Han S, Wang J, Rizos C (2001) A study on GPS/GLONASS multiple reference station techniques for

precise real-time carrier phase-based positioning. Proceedings of ION GPS 2001, Salt Lake City, UT,

September 2001, pp. 392-403

Dimov I, McKee S (2008) Monte Carlo methods for applied scientists. World Scientific, London

Dong D, Bock, Y (1989) Global positioning system network analysis with phase ambiguity resolution applied to

crustal deformation studies in California. J Geophys Res, 94(B4):3949-3966

Doucet A, Freitas, N, Gordon, N (2001) Sequential Monte Carlo Methods in Practice. Springer, New York

Doucet A, Godsill S, Andrieu C (2000) On sequential Monte Carlo sampling methods for Bayesian filtering.

Statistics and Computing, 10(3):197-208

Dow J, Neilan R, Rizos C (2009) The international GNSS service in a changing landscape of global navigation

satellite systems. Journal of Geodesy, 83(3-4):191-198

Euler H, Schaffrin B (1991) On a measure of discernibility between different ambiguity solutions in the

static-kinematic GPS-mode. Proceedings of the international symposium on kinematic systems in

Geodesy, surveying, and remote sensing. Springer, Berlin Heidelberg New York, pp. 285-295

Force D, Miller J (2013) Combined Global Navigation Satellite Systems in the Space Service Volume.

Proceedings of the ION International Technical Meeting, San Diego, USA, January 2013

Fox D (2003) Adapting the sample size in particle filters through KLD-sampling. The International Journal of

Robotics Research, 22(12):985-1003

Frei E, Beutler G (1990) Rapid static positioning based on the fast ambiguity resolution approach FARA: theory

and first results. Manuscripta Geodaetica, 15:325-356

Ge M, Calais E, Hasse J (2000) Reducing satellite orbit error effects in near real-time GPS zenith tropospheric

delay estimation for meteorology. Geophysical Research Letters, 27(13):1915-1918

Ge M, Gendt G, Rothacher M, Shi C, Liu J (2008) Resolution of GPS carrier-phase ambiguities in precise point

positioning (PPP) with daily observations. Journal of Geodesy, 82(7):389-399

Ge M, Zhang H, Jia X, Song S, Wickert J (2012) What is Achievable with Current COMPASS Constellations?

Proceedings of ION GNSS 2012, Nashville, TN, September 2012, pp. 331-339

Gordon J, Salmond J, Smith F (1993) Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE

Proceedings on Radar and Signal Processing, 140(2):107-113

Gustafsson F (2010) Particle filter theory and practice with positining applications. IEEE Aerospace and

Electronic Systems Magazine, 25(7):53-82

Gustafsson F, Gunnarsson F, Bergman N, Forssell U, Jansson J, Karlsson R, Nordlund P (2002) Particle filters

for Positioning, Navigation, and Tracking. IEEE transactions on Signal Processing, 50(2):425-437

Hadas T, Bosy J (2015) IGS RTS precise orbits and clocks verification and quality degradation over time. GPS

Solutions, 19:93-105

Hahn J, Powers E (2005) Implementation of the GPS to Galileo time offset (GGTO). Proceedings of 2005 joint

IEEE international frequency control symposium and precise time time interval (PPTI) systems &

applications meeting, Vancouver, Canada, pp. 33-37.

Han, S, Rizos, C (1996) Improving the computational efficiency of the ambiguity function algorithm. Journal of

Geodesy, 70(6):330-341

Han S, Dai L, Rizos C (1999) A new data processing strategy for combined GPS/GLONASS carrier phase-based

positioning. Proceedings of ION GPS 1999, Nashville, TN, September 1999, pp. 1619-1627

Haug A (2012) Bayesian Estimation and Tracking: A Practical Guide. Wiley & Sons Press, New Jersey

He K, Xu G, Xu T, Flechtner F (2016) GNSS navigation and positioning for the GEOHALO experiment in Italy.

GPS Solutions, 20(2):215-224

Hegarty C, Powers E, Fonville B (2004) Accounting for timing biases between GPS, modernized GPS, and

Galileo signals. Proceedings of 36th Annual Precise Time and Time Interval (PTTI) Meeting,

Washington, D.C., December 2004, pp. 307-317

Hein G (2006) Achieving a Global System of Systems or “Does Everything Have to Be the Same?”. Inside

GNSS, 1(1):57-60

Hein G, Godet J, Issler J, Martin J, Erhard P, Lucas-Rodriguez R, Pratt T (2002) Status of Galileo frequency and

signal design. Proceedings of ION GPS 2002, Portland, OR, September 2002, pp. 266-277

Hernández-Pajares M, Juan J, Sanz J, Aragón-Àngel A, GarcíaRigo A, Salazar D, Escudero M (2011) The

ionosphere: effects, GPS modeling and the benefits for space geodetic techniques. Journal of Geodesy,

85(12):887-907

Hofmann-Wellenhof B, Lichtenegger H, Wasle E (2007) GNSS–global navigation satellite systems: GPS,

GLONASS, Galileo, and more. Springer, Austria

Hol J, Schon T, Gustafsson F (2006 ) On resampling algorithms for particle filters. Processing of Nonlinear

Statistical Signal Processing Workshop, Cambridge, UK, pp. 79-82

IAC_GLO (2015) GLONASS constellation status [Online]. Retrieved from Information-analytical centre official

website: https://glonass-iac.ru/en/GLONASS/ [Assessed 06 May 2015]

IAC_GPS (2015) GPS constellation status [Online]. Retrieved from Information-analytical centre official

website: https://glonass-iac.ru/en/GPS/ [Assessed 06 May 2015]

ICD-BDS (2013) BeiDou navigation satellite system signal in space interface control document open service

signal, Version 2.0. China Satellite Navigation Office, Beijing

ICD-Galileo (2015) European GNSS (Galileo) open service signal in space interface control document, Issue 1.2.

European Union

ICD-GLONASS (2008) Global Navigation Satellite System GLONASS interface control document, version 5.1.

Russian Institute of Space Device Engineering, Moscow

ICD-IRNSS (2014) Indian Regional Navigational Satellite System signal in space ICD for standard positioning

service, version 1.0. Indian Space Research Organization Satellite Center, Bangalore

Ineichen D, Brockmann E, Schaer S (2008) Processing combined GPS/GLONASS data at swisstopo’s local

analysis center. Proceedings of EUREF Symposium, Brussels, Belgium, June 2008

IS-QZSS (2014) Quasi Zenith Navigation System Satellite Service: Interface Specification for QZSS, version 1.6.

Japan Aerospace Exploration Agency (JAXA)

ISRO (2014) IRNSS SIS ICD for standard positioning service. Bangalore, India

Janes H, Langley R, Newby S (1991) Analysis of tropospheric delay prediction models: comparisons with

ray-tracing and implications for GPS relative positioning. Bulletin Geodesique, 65(3):151-161

Jin S, Wang J (2004) Impacts of stochastic modeling on GPS-derived ZTD estimations. Proceedings of ION

GNSS 2004, Long Beach, CA, September 2004, pp. 941-946

Julien O, Alves P, Cannon E, Zhang W (2003) A tightly coupled GPS/GALILEO combination for improved

ambiguity resolution. Proceedings of the European Navigation Conference, pp. 1-14

Julier S, Uhlmann J (1997) A new extension of the Kalman filter to nonlinear systems. Proceedings of

AeroSense: The 11th International Symposium on Aerospace/Defense Sensing, Simulation and

Controls, Orlando, Florida, 1997, pp. 182-193

Kim D, Langley R (2000) GPS ambiguity resolution and validation: methodologies, trends and issues.

Proceedings of the 7th GNSS Workshop–International Symposium on GPS/GNSS, Seoul, Korea, 2000,

pp. 213-221

King R, Bock Y (1999) Documentation for the GAMIT GPS analysis software. Massachusetts Institute of

Technology, Cambridge MA

Kitagawa G (1996) Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. Journal of

Computational and Graphical Statistics, 5(1):1-25

Kotecha J, Djuric P (2003) Gaussian sum particle filtering. IEEE Transactions on Signal Processing,

51(10):2602-2612

Kozlov D, Tkachenko M (1997) Instant RTK cm with Low Cost GPS+ GLONASS [TM] C/A Receivers.

Proceedings of ION GPS 1997, Kansas City, MO, September 1997, pp. 1559-1570

Kozlov D, Tkachenko M (1998) Centimeter‐level, real‐time kinematic positioning with GPS+ GLONASS C/A

receivers. Navigation, 45(2):137-147

Leick A (1998) GLONASS satellite surveying. Journal of Surveying Engineering, 124(2):91-99

Leick A (2004) GPS Satellite Surveying. John Wiley & Sons, New Jersey

Lekkerkerk H (2015) 2020 and beyond. Geoinformatics, 18(5):28-29

Li T, Wang J (2011) Comparing the mathematical models for GPS and GLONASS integration. International

Global Navigation Satellite Systems Society Symposium-IGNSS, Sydney, Australia, November 2011

Li X, Ge M, Dai X, Ren X, Mathias F, Jens W (2015) Accuracy and reliability of multi-GNSS real-time precise

positioning: GPS, GLONASS, BeiDou, and Galileo. Journal of Geodesy, 89(6):607-635

Li X, Ge M, Zhang H, Wickert, J (2013) A method for improving uncalibrated phase delay estimation and

ambiguity-fixing in real-time precise point positioning. Journal of Geodesy, 87(5):405-416

Liu J, Chen R (1998) Sequential Monte Carlo methods for dynamic systems. Journal of the American Statistical

Association, 93(443):1032-1044

Mayer C, Becker C, Jakowski N, Meurer M (2011) Ionosphere monitoring and inter-frequency bias

determination using Galileo: First results and future prospects. Advances in Space Research,

47(5):859-866

Meindl M (2011) Combined analysis of observations from different Global Navigation Satellite Systems. PhD

thesis. Astronomisches Institut der Universität Bern, Geodätisch-geophysikalische Arbeiten in der

Schweiz. Vol.83

Melgard T, Tegedor J, Jong K, Lapucha D, Lachapelle G (2013) Interchangeable integration of GPS and Galileo

by using a common system clock in PPP. Proceedings of ION GNSS, Nashville, TN, September 2013,

pp. 16-20

Montenbruck O, Steigenberger P, Khachikyan R, Weber G, Langley R, Mervart L, Hugentobler U (2013)

IGS-MGEX: preparing the ground for multi-constellation GNSS science. Proceedings of the 4th

international colloquium scientific and fundamental aspects of the Galileo programme, Prague, Czech

Republic, December 2013

NOAA (2015) GPS Space Segment [Online]. Retrieved from Official US Government information about the

Global Positioning System (GPS) and related topics website: http://www.gps.gov/systems/gps/space/,

[Assessed 06 May 2015]

Nurmi J, Lohan S, Sand S, Hurskainen H (2015) GALILEO positioning technology. Springer, Netherlands

Odijk D, Teunissen P (2013a) Characterization of between-receiver GPS-Galileo inter-system biases and their

effect on mixed ambiguity resolution. GPS Solutions, 17(4):521-53

Odijk D, Teunissen P (2013b) Estimation of differential inter-system biases between the overlapping frequencies

of GPS, Galileo, BeiDou and QZSS. Proceedings of the 4th International Colloquium Scientific and

Fundamental Aspects of the Galileo Programme, Prague, Czech Republic, December 2013

Odijk D, Teunissen P, Khodabandeh A (2014) Galileo IOV RTK positioning: standalone and combined with

GPS. Survey Review, 46(337):267-277

Odolinski R, Teunissen P, Odijk D (2014) Combined BDS, Galileo, QZSS and GPS single-frequency RTK. GPS

Solutions, 19(1):151-163

O'Driscoll C (2010) GNSS Solutions: Carrier phase and its measurement for GNSS. Inside GNSS, 5(5):18-22

Paziewski J, Sieradzki R, Wielgosz P (2015) Selected properties of GPS and Galileo-IOV receiver intersystem

biases in multi-GNSS data processing. Measurement Science and Technology, 26(9):095008

Paziewski J, Wielgosz P (2015) Accounting for Galileo-GPS inter-system biases in precise satellite positioning.

Journal of Geodesy, 89(1):81-93

Pitt M, Shephard N (1999) Filtering via simulation: Auxiliary particle filters. Journal of the American Statistical

Association, 94(446):590-599

Povaliaev A (1997) Using Single Differences for Relative Positioning in GLONASS. Proceedings of ION GPS

1997, Kansas City, MO, September 1997, pp. 929-934

Pratt M, Burke B, Misra P (1998) Single-epoch integer ambiguity resolution with GPS-GLONASS L1-L2 Data.

Proceedings of ION GPS, Nashville, TN, September 1998, pp. 389-398

Reussner N, Wanninger L (2011) GLONASS inter-frequency biases and their effects on RTK and PPP

carrier-phase ambiguity resolution. Proceedings of ION GNSS 2011, Portland, OR, September 2011,

pp. 712-716

Rizos C, Janssen V, Roberts C, Grinter T (2012) Precise Point Positioning: Is the era of differential GNSS

positioning drawing to an end? Proceedings of FIG working week, Rome, Italy, May 2012

Rizos C, Montenbruck O, Weber R, Weber G, Neilan R, Hugentobler U (2013) The IGS MGEX experiment as a

milestone for a comprehensive multi-GNSS service. Proceedings of ION PNT 2013, Honolulu, Hawaii,

April, 2013

Saastamoinen J (1973) Contributions to the theory of atmospheric refraction. Bulletin Geodesique,

105(1):279-298

Sleewagen J, Simsky A, Wilde W, Boon F, Willems T (2012) Demystifying GLONASS inter-frequency carrier

phase biases. InsideGNSS, 7(3):57-61

Steigenberger P, Hugentobler U, Loyer S, Perosanz F, Prange L, Dach R, Uhlemann M, Gendt G, Montenbruck

O (2015) Galileo orbit and clock quality of the IGS Multi-GNSS experiment. Advances in Space

Research, 55(1):269-281

Stupak G (2010) GLONASS status and development plans. Proceedings of 5th Meeting of the International

Committee on GNSS, 2010, Turin, Italy, 2010

Takac F (2009) GLONASS inter-frequency biases and ambiguity resolution. Inside GNSS, 4(2):24-28

Takasu T, Yasuda A (2009) Development of the low-cost RTK-GPS receiver with an open source program

package rtklib. Proceedings of international symposium on GPS/GNSS, Jeju, Korea, 2009

Tan P, Steinbach M, Kumar V (2006) Introduction to data mining. Pearson Addison Wesley, Boston

Teunissen P (1995) The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer

ambiguity estimatio. Journal of Geodesy, 70(1-2):65-82

Teunissen P (1996) GPS carrier phase ambiguity fixing concepts. In Kleusberg A, Teunissen P, GPS for geodesy,

pp. 263-336, Springer, Berlin Heidelberg New York

Teunissen P (1999) A theorem on maximizing the probability of correct integer estimation. Artificial Satellites,

34(1):3-9

Teunissen P (2002) The parameter distributions of the integer GPS model. Journal of Geodesy, 76(1):41-48

Teunissen P (2003) Integer aperture GNSS ambiguity resolution. Artificial Satellites, 38(3):79-88

Teunissen P (2004) Penalized GNSS ambiguity resolution. Journal of Geodesy, 78(4-5):235-244

Teunissen P, Kleusberg A (1996) GPS observation equations and positioning concepts. In Kleusberg A,

Teunissen P, GPS for Geodesy, Springer-Verlag, Berlin Heidelberg New York, pp. 175-218

Teunissen P, Verhagen S (2009) The GNSS ambiguity ratio-test revisited: a better way of using it. Survey

Review, 41(312):138-151

Teunissen P, Verhagen S (2007) On GNSS ambiguity acceptance tests. Proceedings of International Global

Navigation Satellite Systems Society Symposium 2007, Sydney, Australia, December 2007

Tian Y, Ge M, Neitzel F (2015) Particle filter-based estimation of inter-frequency phase bias for real-time

GLONASS integer ambiguity resolution. Journal of Geodesy, 89(11):1145-1158

Tiberius C, de Jonge P (1995) Fast positioning using the LAMBDA method. Proceedings of 4th International

Conference on Differential Satellite Navigation Systems, Bergen, Norway, April 1995

Verhagen S (2005) On the reliability of integer ambiguity resolution. Navigation, 52(2):99-100

Verhagen S, Teunissen P (2013) The ratio test for future GNSS ambiguity resolution. GPS Solutions,

17(4):535-548

Verhagen S, Tiberius C, Li B, Teunissen P (2012) Challenges in ambiguity resolution: biases, weak models, and

dimensional curse. Proceedings of 6th ESA Satellite on Satellite Navigation Technologies and

European Workshop on GNSS Signals and Signal Processing, Noordwijk, Netherlands, December

2012, pp. 1-8

Wan E, Der Merwe R (2000) The unscented Kalman filter for nonlinear estimation. Proceedings of IEEE

Symposium on Adaptive Systems for Signal Processing Communications and Control (AS-SPCC),

Lake Louise, Canada, October 2000, pp. 153-158

Wang J (2000) An approach to GLONASS ambiguity resolution. Journal of Geodesy, 74(5):421-430.

Wang J, Rizos C, Stewart M, Leick A (2001) GPS and GLONASS integration: modeling and ambiguity

resolution issues. GPS Solutions, 5(1):55-64

Wang J, Stewart M, Tsakiri M (1998) A discrimination test procedure for ambiguity resolution on-the-fly.

Journal of Geodesy, 72(11):644-653

Wang L, Verhagen S (2015) A new ambiguity acceptance test threshold determination method with controllable

failure rate. Journal of Geodesy, 89(4):361-375

Wanninger L (2012) Carrier-phase inter-frequency biases of GLONASS receivers. Journal of Geodesy,

86(2):139-148

Wanninger L, Wallstab-Freitag S (2007) Combined processing of GPS, GLONASS, and SBAS code phase and

carrier phase measurements. Proceedings of ION GNSS 2007, Fort Worth, TX, September 2007,

pp. 866-875.

Yamada H, Takasu T, Kubo N, Yasuda A (2010) Evaluation and calibration of receiver inter-channel biases for

RTK-GPS/GLONASS. Proceedings of ION GNSS 2010, Portland, OR, September 2010, pp. 1580-1587

Yang Y, He H, Xu G (2001) Adaptively robust filtering for kinematic geodetic positioning. Journal of Geodesy,

75(2-3):109-116

Zhang S, Zhang K, Wu S, Li B (2011) Network-based RTK positioning using integrated GPS and GLONAS

observations. International Global Navigation Satellite Systems Society Symposium-IGNSS, Sydney,

Australia, 2011

Zhao F, Shin J, Reich J (2002) Information-driven dynamic sensor collaboration. IEEE Signal Processing

Magazine, 19(2):61-72

Zhao T, Nehorai A (2007) Distributed sequential Bayesian estimation of a diffusive source in wireless sensor

networks. IEEE Transactions on Signal Processing, 55(4):1511-1524

Zhao Y (2014) Adaptive Particle Filters for Wireless Indoor Target Tracking. PhD thesis, Institut für Informatik,

Freie Universität Berlin, Germany

Zinoviev A (2005) Using GLONASS in combined GNSS receivers: current status. Proceedings of ION GNSS

2005, Long Beach, CA, September 2005, pp. 1046-1057

Zinoviev A, Veitsel A, Dolgin D (2009) Renovated GLONASS: improved performances of GNSS receivers.

Proceedings of ION GNSS 2009, Savannah, GA, September 2009, pp. 3271-3277

List of Tables

Table 2.1 Satellites of GPS on 5th June 2015 ........................................................................................................... 7

Table 2.2 Satellites of GLONASS on 5th June 2015 ............................................................................................... 8

Table 2.3 Satellites of Galileo on 5th June 2015 ...................................................................................................... 9

Table 2.4 Frequency distribution of the satellite systems ..................................................................................... 10

Table 4.1 Statistics of the convergence time for the IFB rate estimated by the new approach for the

three baselines ...................................................................................................................................... 36

Table 4.2 Lengths and estimated IFB rates for the three baselines on DOY 191 of 2014 and DOY 036

of 2015 ................................................................................................................................................. 45

Table 4.3 Receivers and antennas for the three baselines ..................................................................................... 45

Table 5.1 F-ISB estimation results of short baselines in MGEX .......................................................................... 64

Table 5.2 Receiver types and firmware series for each station in the short baselines in MGEX .......................... 65

List of Figures

Fig. 3.1 PDF of the normal distribution 𝒩(𝑥;0,1) within [-3, 3] (left), and the corresponding CDF

(right) ....................................................................................................................................................... 24

Fig. 3.2 Illustration of the SIR method. The positions of the dots represent the samples while the sizes

represent the weights ................................................................................................................................ 25

Fig. 4.1 Three-dimensional RATIO distribution along with epochs and IFB rate samples for the

zero-baseline (a), KOSG-KOS1 (b) and REF6-AIR5 (c) (The part corresponding to RATIOs

which are larger than 50 are not depicted) ............................................................................................... 31

Fig. 4.2 RATIO values of the first epoch corresponding to different IFB rate samples for the

zero-baseline (a), KOSG-KOG1 (b) and REF6-AIR5 (c) ........................................................................ 31

Fig. 4.3 Distribution of the points with RATIO larger than the threshold of 3 for the zero-baseline (a),

KOSG-KOS1 (b) and REF6-AIR5 (c) ..................................................................................................... 32

Fig. 4.4 RATIO plot at the epoch 112 (a), 323 (b) and 800 (c) for baseline REF6-AIR5 ...................................... 32

Fig. 4.5 Statistics of the relationship of the maximum RATIO and the IFB rate over all epochs for the

three baselines, the zero-baseline (a), KOSG-KOS1 (b), and REF6-AIR5 (c) ........................................ 32

Fig. 4.6 Flow chart of the procedure for IFB estimation, where 𝑠𝑡𝑑𝑡ℎ𝑑 is the STD threshold value .................... 34

Fig. 4.7 Estimated IFB rates using the particle filter for the zero-baseline (a), KOSG-KOS1 (b) and

REF6-AIR5 (c)......................................................................................................................................... 35

Fig. 4.8 Convergence of the estimated IFB rate (solid line) and STD (dash lines) for zero-baseline (a),

KOSG-KOS1 (b) and baseline REF6-AIR5 (c), the sampling rate is 5s, 30s, 1s, respectively ............... 35

Fig. 4.9 Convergence process of the estimated IFB rate versus the number of epochs for the zero-

baseline (a), KOSG-KOS1 (b) and REF6-AIR5 (c). The star symbols denote the converging

point. The sampling rate is 5 s, 30 s, and 1 s, respectively....................................................................... 36

Fig. 4.10 Statistics of the epochs needed for the convergence of IFB rate for the zero-baseline (a),

KOSG-KOS1 (b) and REF6-AIR5 (c). The sampling rate is 5 s, 30 s and 1 s, respectively.................... 36

Fig. 4.11 The statistics of the epochs needed for the convergence of IFB rate for baseline STR1-STR2

with different type of receivers. The sampling rate is 1 s (a), 5 s (b) and 30 s (c), respectively .............. 36

Fig. 4.12 Three-dimensional RATIO distribution (left) and the estimation results (right) for baseline

KOSG-APEL ........................................................................................................................................... 37

Fig. 4.13 Computational time for a single epoch including the particle filter for baseline KOSG-KOS1.

The upper thin dash line and the thick solid line close to the bottom are for the search interval

of [-0.10, 0.10] m/FN and [-0.04, 0] m/FN, respectively. The number of satellites is also plotted ......... 37

Fig. 4.14 Estimated IFB rate (left) and the STD (right) for twenty calculations by particle filter with the

state noise level 𝜎𝑛 set to different values ............................................................................................... 39

Fig. 4.15 Estimated IFB rate (left) and the STD (right) for twenty calculations by RPF with the state

noise level 𝜎𝑛 set to different values. ...................................................................................................... 40

Fig. 4.16 Estimated IFB rate (left) and the STD (right) with the number of particles N set to different

values. ...................................................................................................................................................... 41

Fig. 4.17 Estimated IFB rate (left) and the STD (right) for 20 calculations with controlled numbers for

24 hours .................................................................................................................................................... 42

Fig. 4.18 Estimated IFB rate (left) and the STD (right) for 20 calculations with number of particles N set

to controlled value and constant value 200 for the first 15 minutes ......................................................... 42

Fig. 4.19 Number of particles N along time for 20 calculations with the method of setting N to

controlled value and constant value 200 for the 24 hours (left) and for the first 15 minutes

(right) ....................................................................................................................................................... 43

Fig. 4.20 Computation time with the number of particles set to controlled value and constant value 200,

along with the number of satellites, for the 24 hours (left) and for the first 15 minutes (right) ............... 43

Fig. 4.21 IFB rates in GLONASS data processing for baselines STR2-STR1 (top) and KOSG-KOS1

(bottom) within about two years .............................................................................................................. 44

Fig. 5.1 Satellite sky plots of GPS (left) and Galileo (right) for station TLSE on day 001 of 2015 ...................... 49

Fig. 5.2 Numbers of GPS and Galileo satellites for baseline TLSG-TLSE on DOY 001 of 2015, with

elevation mask of 10 degrees ................................................................................................................... 49

Fig. 5.3 Relationship between ISB and fixing RATIO for the first epoch with only one Galileo satellite

(a) and the three-dimensional RATIO distribution for all epochs involving Galileo observations

(b) for baseline TLSG-TLSE ................................................................................................................... 50

Fig. 5.4 Three-dimensional RATIO distribution of GPS and Galileo integration (a), as well as the

average values along epoch time (b) for baseline TLSG-TLSE ............................................................... 50

Fig. 5.5 Impact of ISB biases in the baseline fixed solutions of GPS L1 and Galileo E1 integration with

respect to the GPS L1 baseline fixed solution. The data are from baseline TLSG-TLSE at

epoch 1406. The F-ISB sampling interval is 0.001 m .............................................................................. 50

Fig. 5.6 Satellite sky plots of GPS (left) and GLONASS (right) for station KOS1 on DOY 048 of 2014 ............ 52

Fig. 5.7 Numbers of GPS and GLONASS satellites for baseline KOSG-KOS1 on DOY 048 of 2014,

with elevation mask of 10 degrees ........................................................................................................... 52

Fig. 5.8 Baseline solutions from code pseudorange observations of GLONASS L1 for baseline

KOSG-KOS1without (left) and with (right) code pseudorange SD-IFB correction ................................ 52

Fig. 5.9 Three-dimensional RATIO distribution of GPS L1 and GLONASS L1 integration with

inter-system models for baseline KOSG-KOS1 ....................................................................................... 53

Fig. 5.10 Average values (a) along epoch time axis of RATIO values in Fig. 5.9 and the zoom-in of the

average values (b), where the red short lines mark the width of the peaks at RATIO of value 3 ............ 53

Fig. 5.11 Comparison of RATIO values in integration of full GPS and GLONASS constellations, full

GPS constellation and only one GLONASS satellite (top), as well as only one GPS satellite and

full GLONASS constellation (bottom) .................................................................................................... 54

Fig. 5.12 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration

with respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at

epoch 323. The F-ISB sampling interval is 0.02 m .................................................................................. 55

Fig. 5.13 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration

with respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at

epoch 2871. The F-ISB sampling interval is 0.02 m ................................................................................ 55

Fig. 5.14 Impact of ISB biases in the baseline fixed solutions of GPS L1 and GLONASS L1 integration

with respect to the GPS L1 baseline fixed solution. The data are from baseline KOSG KOS1 at

epoch 2871. The F-ISB sampling interval is 0.001 m .............................................................................. 55

Fig. 5.15 Satellite sky plots of BDS for station TLSE on DOY 001 of 2015 ......................................................... 57

Fig. 5.16 Numbers of GPS and BDS satellites for baseline TLSG-TLSE on DOY 001 of 2015, with

elevation mask of 10 degrees ................................................................................................................... 57

Fig. 5.17 Three-dimensional RATIO distribution for baseline TLSG-TLSE on DOY 001 of 2015 for all

epochs (a) and averages over epochs (b) .................................................................................................. 58

Fig. 5.18 Impact of ISB biases in the baseline fixed solutions of GPS L1 and BDS B1 integration with

respect to the GPS L1 baseline fixed solution. The data are from baseline TLSG-TLSE at

epoch 1703. The F-ISB sampling interval is 0.02 m ................................................................................ 58

Fig. 5.19 Averages of RATIO values along epoch time axis (a), as well as the impact of ISB biases on

the baseline fixed solutions with GPS L1 and BDS B1 integration compared to the values from

GPS L1 solutions (b) for baseline TLSG-TLSE. Result shown in (b) employs data at epoch

1703. The F-ISB sampling interval is 0.001 m ........................................................................................ 58

Fig. 5.20 Three-dimensional RATIO map for baseline GOP6-GOP7. The F-ISB is 0.095m which is very

close to 0.5 cycles, and there are two RATIO peaks within the initial interval [-0.1, 0.1] m .................. 59

Fig. 5.21 Estimated F-ISB and the three times STD for integration of GPS L1 and Galileo E1 for the

baseline TLSG-TLSE on the DOY 001of 2015 ....................................................................................... 61

Fig. 5.22 Estimated F-ISB and the three times STD for integration of GPS L1 and GLONASS L1 for

baseline KOSG-KOS1 on DOY 048 of 2014, with approximate ISB -16.2764 m. ................................. 61

Fig. 5.23 Estimated F-ISB and the three times STD for integration of GPS L1 and GLONASS L1 for

baseline TLSG-TLSE on DOY 001 of 2015, with approximate ISB -4.943 m. ....................................... 61

Fig. 5.24 Estimated F-ISB and the three times STD for integration of GPS L1 and BDS B1 for baseline

TLSG-TLSE on DOY 001 of 2015, with approximate ISB -0.3007 m .................................................... 61

Fig. 5.25 Estimated F-ISB (a) and the corresponding STD (b) without the cluster analysis and results (c),

(d) with cluster analysis for GPS L1 and Galileo E1 integration for baseline GOP6-GOP7 on

DOY 001 of 2015 ..................................................................................................................................... 62

Fig. 5.26 Computation time of Particle filter for KOSG-KOS1 GPS and GLONASS L1 integration (left).

Computation time of Particle filter for TLSG-TLSE GPS L1 and Galileo E1 integration (right) ........... 63

Fig. 5.27 F-ISB of GPS L1 and Galileo E1 integration within about one year for HARB-HRAG (a),

RGDG-RIO2 (b) and TLSE-TLSG (c)..................................................................................................... 66

Fig. 5.28 F-ISB of GPS L1 and Galileo E1 integration within 5 days around the jump epoch for

HARB-HRAG (a), RGDG-RIO2 (b) and TLSE-TLSG (c) ...................................................................... 66

Fig. 5.29 F-ISB of GPS L1 and GLONASS L1 integration with approximate ISB set to -16.2764m

within about two years after DOY 167 of 2013 (top), as well as short-term F-ISB values within

40 days with available data (bottom) for baseline KOSG-KOS1 ............................................................. 67

Fig. 6.1 RATIO distribution with two F-ISBs of GPS L1 and Galileo E1, GPS L1 and BDS B1 for

epoch 476 for baseline TLSG-TLSE, with all observed satellites (left) and later with only 3

GPS, 2 Galileo and 2 BDS satellites (right) ............................................................................................. 69

Fig. 6.2 Convergence process with the two-dimensional particle filter approach for epochs from one to

seven (a-g) and the estimated F-ISB results for GPS L1 and BDS B1 combination (h, top), GPS

L1 and Galileo E1 combination (h, bottom) ............................................................................................ 71

Fig. 6.3 Estimated F-ISB for GPS L1 and BDS B1 combination (top), as well as GPS L1 and Galileo E1

combination (bottom) for baseline TLSG-TLSE ..................................................................................... 71

Fig. 6.4 Satellite PRN numbers (left) and the estimated F-ISB results for GPS L1 and BDS B1 (right

top), GPS L1 and Galileo E1 (right bottom) with only two satellites from each system ......................... 72

Fig. 6.5 Positioning differences with respect to the GPS static solution for the strategy fixing only intra-

system DD-ambiguities and that fixing both intra- and inter-system DD-ambiguities, which are

denoted as Only intra DDAF (DD-Ambiguity Fixing) and Intra and inter DDAF in the figures,

respectively .............................................................................................................................................. 72

Fig. 7.1 Comparison of the GLONASS single-epoch solution with (blue) and without (red) IFB rate

estimation procedure at the beginning for baseline KOSG-KOS1 ........................................................... 74

Fig. 7.2 Position differences with respect to the ground truth for GLONASS without (red) and with

(blue) the initial calibration of IFB rate in kinematic mode. The right panel is a snapshot of the

first 30 epochs (15 minutes) of the left one in order to show the process of the IFB rate

estimation and the first ambiguity-fixing afterwards ............................................................................... 74

Fig. 7.3 Comparison of GLONASS fixed solution and GPS fixed solution .......................................................... 74

Fig. 7.4 IFB rate estimates for all the 72 sessions for baseline KOSG-KOS1, as well as the

corresponding STD .................................................................................................................................. 75

Fig. 7.5 Convergence processes of the GLONASS baseline solutions for all 72 sessions with and

without IFB rate estimation procedure for baseline KOSG-KOS1 .......................................................... 75

Fig. 7.6 Final horizontal positions of the GLONASS baseline solutions for all 72 sessions with and

without IFB rate estimation procedure for baseline KOSG-KOS1 .......................................................... 76

Fig. 7.7 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS static

solution (left), as well as the RATIO values (right) with the strategy fixing only GPS L1 intra-

system DD-ambiguities (green) and the strategy fixing GPS L1 and Galileo E1 both intra- and

inter-system DD-ambiguities (blue) ......................................................................................................... 76

Fig. 7.8 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS static

solution (left), as well as the RATIO values (right) with the strategy fixing GPS L1 and Galileo

E1 only intra-system DD-ambiguities (red) and the strategy fixing GPS L1 and Galileo E1 both

intra- and inter-system DD-ambiguities (blue) ........................................................................................ 76

Fig. 7.9 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static

solution (left), as well as the RATIO values (right) with the strategy fixing GPS L1 and BDS

B1 only intra-system DD-ambiguities (red) and the strategy fixing GPS L1 and BDS B1 both

intra- and inter-system DD-ambiguities (blue) ........................................................................................ 78

Fig. 7.10 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static

solution, employing three strategies, using GPS L1 only (green), using GPS L1, Galileo E1 and

BDS B1 observations fixing intra-system DD-ambiguities (red), and using the three

constellations but fixing both intra- and inter-system DD-ambiguities (blue) ......................................... 78

Fig. 7.11 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static

solution (left) and the EARs of integer ambiguity fixing with different elevation masks (right) ............ 79

Fig. 7.12 Position differences of the TLSG-TLSE baseline solutions with respect to the GPS only static

solution (left) and the EARs of integer ambiguity fixing with different elevation masks (right) ............ 79

Fig. 7.13 First-fixing time with and without inter-system DD-ambiguity fixing for the processing

scenarios of GPS L1 and Galileo E1 ........................................................................................................ 79

Fig. 7.14 Satellite PRN/slot numbers (left) and the differences of the baseline solutions with respect to

the GPS only static solution (right) for baseline KOSG-KOS1, with three GPS satellites and

three GLONASS satellites. ...................................................................................................................... 80

Fig. 7.15 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static

solution, with three GPS satellites and three GLONASS satellites .......................................................... 81

Fig. 7.16 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static

solution, with satellites in the west half-sky ............................................................................................ 81

Fig. 7.17 Position differences of the KOSG-KOS1 baseline solutions with respect to the GPS only static

solution, with satellites in the east half-sky .............................................................................................. 81

List of Abbrevations

Three-Dimensional

AFM

Ambiguity Function Method

BDS

BeiDou Navigation Satellite System

CDF

Cumulative Distribution Function

CDMA

Code Division Multiple Access

Double Difference

DDAF

Double Difference Ambiguity Fixing

DOY

Day Of Year

DSP

Digital Signal Processing

EAR

Empirical Availability Rate

European Commission

EKF

Extended Kalman Filter

ESA

European Space Agency

FDMA

Frequency Division Multiple Access

F-ISB

Fractional Inter-System Bias

Frequency Number

FOC

Full Operational Capability

GEO

Geostationary Orbits

GLONASS

GLObal NAvigation Satellite System

GNSS

Global Navigation Satellite System

GPS

Global Positioning System

GST

Galileo System Time

GTRF

Galileo Terrestrial Reference Frame

IAC

Information Analysis Centre

IFB

Inter-Frequency Bias

ICD

Interface Control Document

IGS

International GNSS Service

IGSO

Inclined Geosynchronous Satellite Orbit

ILS

Integer Least Square

IOV

In-Orbit Validation

IRNSS

Indian Regional Navigation Satellite System

ISB

Inter-System Bias

ITRF

International Terrestrial Reference Frame

JGS

Japanese Geodetic System

LAMBDA

Least-squares AMBiguity Decorrelation Adjustment

MEO

Medium Earth Orbit

MGEX

Multi-GNSS Experiment

NEQ

Normal Equation

OLS

Ordinary Least Square

Personal Computer

PDF

Probability Density Function

PPP

Precise Point Positioning

PRN

Pseudo-Random Noise

PVT

Positioning, Navigation and Timing

QZSS

Quasi-Zenith Satellite System

RPF

Regularized Particle Filter

Selective Availability

Single Difference

SIS

Sequential Importance Sampling

SIR

Sequential Importance Resampling

SPP

Single Point Positioning

STD

STandard Deviation

TAI

International Atomic Time

UKF

Unscented Kalman filter

United States

USSR

Union of Soviet Socialist Republics

UTC

Coordinated Universal Time

Variance-Covariance

WGS

World Geodetic System

Acknowledgments

I would like to thank all the people who encouraged and supported me during my PhD study.

Firstly, I would like to express my profound respect and gratitude to Prof. Dr.-Ing. Frank Neitzel at Institute of

Geodesy and Geoformation Science of Technische Universität Berlin (TU Berlin), for giving me the opportunity

to pursue my PhD study at TU Berlin. He has given me valuable continuous support, guidance and

encouragement during my PhD study in the past years. I would also like to express my profound respect and

gratitude to Dr. Maorong Ge at German Research Centre for Geosciences (GFZ) for his professional insight and

advices on my study which inspires my research interest and enthusiasm. Without his help, this thesis would

hardly have been completed.

I am grateful to Prof. Dr.-Ing. Guochang Xu, formerly at GFZ, now at Shandong University in China, and

Professor Dr.-Ing. Jianjun Zhu at Central South University in China for all the effort for getting me the chance to

study at TU-Berlin and for their guidance in my study. I would like to thank Dr. Svetozar Petrovic for polishing

the English language of the thesis, and Dr. Dietrich Ewert for his help during my first-year study in Berlin.

I also wish to thank my colleagues and friends at TU Berlin, especially Daniel Wujanz, Sven Weisbrich,

Georgios Malissiovas, and Mathias Burger, for their help on various procedures and language improvement, as

well as Jinpeng Feng, Nan Jiang, Xunqiang Gong, Hui Tang, Mohammed Matalqah for their kindly help and

discussions.

I am also grateful to the colleagues and friends at GFZ, to Dr. Kaifei He, who has helped and encouraged me a

lot in the past several years, to Dr. Guanwen Huang, Dr. Ming Shangguan, Dr. Xingxing Li, Dr. Zhiguo Deng,

Dr. Rui Tu, Kejie Chen, Yang Liu, Yan Xu and others for their kindly help and discussions. My thanks also go

to Dr. Miao Lin at Leibniz Universität Hannover for his help and encouragement.

Gratefully acknowledged is the China Scholarship Council (CSC) for the financial support over the past four

years. Also thanks to Prof. Dr.-Ing. Frank Neitzel and Prof. Dr.-Ing. Martin Kada for providing me additional

support for the last few months.

Finally, I have to express my deepest gratitude to my parents, especially my father who still works on the farm

and on building sites in his sixties in order to support the family, also my thanks to my brother and my sister, for

their continuous support and understanding.