UBX Parser UBX Offline Corrections SP3, CLK Parser Precise Products DOP calculation EKF Estimated Position SP3, CLK interpolation APO correction Adjust reception time DCB correction Adjust Earth rotation Cycle slip detection Az and El calculation WL ambiguity Solid tides Relativistic effects Plotting  Tool True Position ATX Parser Troposphere calculation Position estimation Ambiguity resolution Precise Point Positioning Multi-Constellation Ionospheric-Free Positioning Using L1/L5 Frequencies on a Low-Cost Receiver Master’s Thesis in Systems, Control and Mechatronics Marielle Melander Philip Pettersson DEPARTMENT OF MECHANICS AND MARITIME SCIENCES CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2024 www.chalmers.se www.chalmers.se Master’s Thesis in Systems, Control and Mechatronics Precise Point Positioning Multi-Constellation Ionospheric-Free Positioning Using L1/L5 Frequencies on a Low-Cost Receiver Marielle Melander Philip Pettersson Department of Mechanics and Maritime Sciences Chalmers University of Technology Gothenburg, Sweden 2024 Precise Point Positioning Multi-Constellation Ionospheric-Free Positioning Using L1/L5 Frequencies on a Low-Cost Receiver Marielle Melander Philip Pettersson © Marielle Melander, Philip Pettersson, 2024. Supervisor: Lucas Rochard, CPAC Systems AB Examiner: Peter Forsberg, Department of Mechanics and Maritime Sciences Master’s Thesis 2024 Department of Mechanics and Maritime Sciences Chalmers University of Technology SE-412 96 Gothenburg Sweden Telephone +46 31 772 1000 Cover: Overview of implemented Precise Point Positioning solution Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Gothenburg, Sweden 2024 iv Precise Point Positioning Multi-Constellation Ionospheric-Free Positioning Using L1/L5 Frequencies on a Low-Cost Receiver Marielle Melander Philip Pettersson Department of Mechanics and Maritime Sciences Chalmers University of Technology Abstract This thesis revolves around the area of satellite-based positioning using Global Naviga- tion Satellite System (GNSS) signals from the Global Positioning System (GPS) and Galileo constellations. GNSS positioning plays a vital role in daily life, particularly in applications like navigation and timing. This project aims to develop a high-accuracy, Precise Point Positioning (PPP) solution using data from a low-cost receiver. To achieve this an ionosphere-free measurement combination based on L1 and L5 frequencies is ex- plored instead of the conventional choice of L1 and L2. The research investigates the impact of different error sources, such as biases and atmospheric disturbances, on the accuracy of PPP. An Extended Kalman Filter (EKF)-based solution that incorporates offline corrections, cycle slip detection, and ambiguity resolution is proposed. The hard- ware used for data collection is the u-blox evaluation kit EKV-F9P which includes the GNSS receiver module ZED-F9P and an antenna. This low-cost dual-frequency hardware is proven to provide accurate data to implement a high-accuracy PPP solution. When enough satellites are visible, the results demonstrate a decimetre-level accuracy for sta- tionary positioning estimation. However, due to the limited availability of L5 frequency signals, continuous precise positioning is not always possible. Despite this limitation, the work contributes to advancing GNSS positioning technology and lays the groundwork for future enhancements in PPP-based applications. Keywords: PPP, GNSS, EKF, dual-frequency, Ionosphere-free, low-cost, EVK-F9P, u- blox v Acknowledgements We would like to extend our gratitude to our examiner, Peter Forsberg, and our super- visor, Lucas Rochard, for their support and feedback during our thesis project. We also want to thank CPAC Systems AB for providing the opportunity and facilities to perform this master’s thesis. Marielle Melander, Philip Pettersson, Gothenburg, June 2024 Supervisor: Lucas Rochard, CPAC Systems AB Examiner: Peter Forsberg, Department of Mechanics and Maritime Sciences vii List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alpha- betical order: APC Antenna Phase Centre APO Antenna Phase Offset BDS BeiDou Navigation Satellite System CDDIS Crustal Dynamics Data Information System CPA Carrier Phase Ambiguity DCB Differential Code Biases DGNSS Differential GNSS DOP Dilution Of Precision DOY Day Of Year ECEF Earth-centred Earth-Fixed EKF Extended Kalman Filter ENU East, North, Up GF Geometry-Free GLONASS Global Navigation Satellite System GNSS Global Navigation Satellite Systems GPS Global Positioning System HDOP Horizontal Dilution Of Precision IEEE Institute of Electrical and Electronics Engineers IF Ionosphere-Free IERS International Earth Rotation and Reference System Service ITRF International Terrestrial Reference Frame LAMBDA Least-Squares Ambiguity Decorrelation LOS Line of Sight MLAMBDA Modified LAMBDA MW Melbourne-Wübbena NavIC Navigation with Indian Constalletion NL Narrow-Lane NTRIP The Networked Transport of RTCM via Internet Protocol PPP Precise Point Positioning PPP-AR Ambiguity-Resolution PPP PRN Pseudo-Random Noise QZSS Quasi-Zenith Satellite System RF Radio Frequency RTK Real-Time Kinematic ix SD Single Difference SNR Signal-To-Noise Ratio SSR State Space Representation TEC Total Electron Content TRS Terrestrial Reference System VDOP Vertical Dilution Of Precision WARKT Wide Area RTK WL Wide-Lane ZHD Zenith Hydrostatic Delay ZWD Zenith Wet Delay x xii Nomenclature Below is the nomenclature of indices, parameters, and variables that have been used throughout this thesis. Indices s Index for satellite r Index for receiver j Index for band Constants c: 299, 792, 458 Speed of light in vacuum (light celerity) [m/s] h2: 0.6078 Nominal degree 2 Love number l2: 0.0847 Nominal degree 2 Shida number GME: 3.986004418 · 1014 Earth’s gravitational parameter [m3/s2] GM2: 4.9048695 · 1012 Moon’s gravitational parameter [m3/s2] GM3: 1.32712440018 · 1020 Sun’s gravitational parameter [m3/s2] Variables f Frequency [Hz] λ Wavelength [m] ρ Geometric distance [m] P Pseudorange [m] I Ionospheric delay [m] T Tropospheric delay [m] L Carrier phase measurement [Cycles] Φ Carrier phase measurement [m] xiii N Integer Carrier phase ambiguity [Cycles] δb Code hardware bias [m] δB Phase hardware bias [m] w Carrier phase windup effect [Cycles] dt Clock errors [m] r Position in ECEF [m] v Velocity in ECEF [m] ∆ttp Signal travel time using pseudorange [s] ∆ttg Signal travel time using geometric distance [s] az Azimuth angle to the satellite from the receiver [deg] el Elevation angle to the satellite from the receiver [deg] ∆rel Special relativistic effect [m] ωe Earth’s rotation speed [m/s] DCBs i,j Differential Code Bias between band i and j for satellite s DPBs i,j Differential Phase Bias between band i and j for satellite s xiv Contents List of Acronyms ix Nomenclature xii List of Figures xix List of Tables xxi 1 Introduction 1 1.1 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Limitations / Demarcations . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 3 2.1 GNSS signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Radio Frequency (RF) Carrier . . . . . . . . . . . . . . . . . . . . 3 2.1.2 Pseudo-Random Noise (PRN) ranging codes . . . . . . . . . . . . 4 2.1.3 Navigation data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 GNSS Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.1 Pseudorange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.2 Carrier Phase Measurement . . . . . . . . . . . . . . . . . . . . . 6 2.2.3 Combination of signals . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.3.1 Ionosphere-free (IF) . . . . . . . . . . . . . . . . . . . . 7 2.2.3.2 Geometry-free (GF) . . . . . . . . . . . . . . . . . . . . 7 2.2.3.3 Narrow-lane (NL) . . . . . . . . . . . . . . . . . . . . . . 8 2.2.3.4 Wide-lane (WL) . . . . . . . . . . . . . . . . . . . . . . 8 2.2.3.5 Melbourne-Wübbena (MW) . . . . . . . . . . . . . . . . 8 2.3 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.1 Conventional Terrestrial Reference System (TRS) . . . . . . . . . 9 2.3.1.1 International Terrestrial Reference Frame (ITRF) . . . . 9 2.3.1.2 East, North, Up (ENU) . . . . . . . . . . . . . . . . . . 9 2.3.1.3 Geodetic . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 GNSS Error Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.1 Emission and Reception Time . . . . . . . . . . . . . . . . . . . . 10 2.4.2 Earth Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.3 Relativistic Effects . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.4 Cycle Slip . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.5 Antenna Phase Centre (APC) . . . . . . . . . . . . . . . . . . . . 12 2.4.6 Atmospheric Effects . . . . . . . . . . . . . . . . . . . . . . . . . 13 xv Contents 2.4.6.1 Ionosphere . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.6.2 Troposphere . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.7 Multipath . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.8 Receiver Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.9 Clock Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.10 Orbit Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.11 Hardware-induced biases . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.11.1 Differential Code Bias (DCB) . . . . . . . . . . . . . . . 16 2.4.12 Tidal Deformations . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.12.1 Solid Tides . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.5 Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.6 Satellite-Based Positioning . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.6.1 Relative Positioning . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.6.1.1 Differential GNSS (DGNSS) . . . . . . . . . . . . . . . . 19 2.6.1.2 Real-Time Kinematic (RTK) . . . . . . . . . . . . . . . 19 2.6.1.3 Network RTK . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6.1.4 Wide Area RTK (WARTK) . . . . . . . . . . . . . . . . 20 2.6.2 Precise Point Positioning (PPP) . . . . . . . . . . . . . . . . . . . 20 2.6.2.1 Standard PPP . . . . . . . . . . . . . . . . . . . . . . . 21 2.6.2.2 Ambiguity-fixed PPP . . . . . . . . . . . . . . . . . . . 21 2.6.2.3 PPP-RTK . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Evaluation of GNSS Positioning . . . . . . . . . . . . . . . . . . . . . . . 22 2.7.1 Vertical and horizontal accuracy . . . . . . . . . . . . . . . . . . . 22 2.7.2 Dilution Of Precision (DOP) . . . . . . . . . . . . . . . . . . . . . 22 3 Methods 25 3.1 Hardware . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Positioning Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 GNSS Frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4 GNSS Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4.1 Offline Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4.1.1 Cycle Slip . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.1.2 Ionosphere . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.1.3 Clock And Orbit Errors . . . . . . . . . . . . . . . . . . 28 3.4.2 Online Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.2.1 Troposphere . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Data Collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6 Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.7 Ambiguity Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.7.1 Wide-Lane Ambiguity Fixing . . . . . . . . . . . . . . . . . . . . 34 3.7.2 Narrow-Lane Ambiguity Fixing . . . . . . . . . . . . . . . . . . . 35 3.7.3 IF Ambiguity Fixing . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.8 Implementation Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.9 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4 Results 41 4.1 Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Dilution Of Precision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 xvi Contents 4.3 Number of Satellites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.4 Vertical Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.5 Horizontal Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.6 RTKLIB Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5 Discussion 51 5.1 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.2 Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.3 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.4 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6 Conclusion 55 A Appendix 1 I A.1 Transformation matrix from ECEF to ENU . . . . . . . . . . . . . . . . I A.2 Coefficients for hydrostatic and wet mapping functions . . . . . . . . . . I A.3 Setup of EVK-F9P to send raw data . . . . . . . . . . . . . . . . . . . . II xvii Contents xviii List of Figures 2.1 Ranging code, carrier wave and the modulated signal. . . . . . . . . . . . 4 2.2 ECEF and ENU frames [1]. . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Simplified depiction of a receiver antenna with the corresponding antenna phase centre, antenna phase offset and antenna reference point. . . . . . 13 2.4 Atmospheric layers with their height, temperature and electron density [2]. 14 2.5 Simplification of tropospheric delay. . . . . . . . . . . . . . . . . . . . . . 15 2.6 Illustration of multipath effect on GNSS signal. . . . . . . . . . . . . . . 15 2.7 Simple depiction of relative positioning. . . . . . . . . . . . . . . . . . . . 19 2.8 Simple depiction of Network RTK. . . . . . . . . . . . . . . . . . . . . . 20 2.9 Precise positioning system. . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Receiver hardware. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Illustration of the ionospheric effect during 04-03-2024 where the bottom row indicates the hour of the day in UTC+1 [42]. . . . . . . . . . . . . . 30 3.3 Overview of software structure. . . . . . . . . . . . . . . . . . . . . . . . 37 3.4 Overview of EKF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 RTKLIB options that was used. . . . . . . . . . . . . . . . . . . . . . . . 39 4.1 Relativistic clock effect in meters for all Galileo and GPS satellites in each time step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Close-up of relativistic clock effect in meters for all satellites in each time step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.3 Solid Tides effect in meters for each time step. . . . . . . . . . . . . . . . 42 4.4 Earth Rotation effect in meters for 15,000 time steps for one satellite. . . 42 4.5 Clock Correction in meters for each time step. . . . . . . . . . . . . . . . 43 4.6 Troposphere correction for 19,000 time steps for one satellite. . . . . . . . 43 4.7 Detected cycle slips for 70 time steps for one satellite. . . . . . . . . . . . 44 4.8 Antenna Phase Offset Correction in meters for the time steps available in the precise products for one satellite. . . . . . . . . . . . . . . . . . . . . 44 4.9 Horizontal DOP value for each time step using satellites with elevation angle larger than 10◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.10 Vertical DOP value for each time step using satellites with elevation angle larger than 10◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.11 Number of satellites for the different frequency bands L1 and L5. . . . . 46 4.12 Number of satellites broadcasting L1 and L5 above 30° elevation angle. . 46 4.13 Number of usable satellites for the ionosphere-free combination of L1 and L5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.14 The offset in the height plane from the true position for 19,000 time steps. 47 xix List of Figures 4.15 The offset in the East direction from the true position for 19,000 time steps. 48 4.16 The offset in the North direction from the true position for 19,000 time steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.17 The offset in the East direction from the true position for 11,000 time steps. 49 4.18 The offset in the North direction from the true position for 11,000 time steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.19 RTKLIB result in East, North and up directions, indicated by E-W, N-S and U-D respectively. The time in UTC is on the x-axis and the offset from the true position in each direction is shown on the y-axis. . . . . . . 50 4.20 RTKLIB number of satellites. . . . . . . . . . . . . . . . . . . . . . . . . 50 xx List of Tables 2.1 Radar-frequency letter bands. . . . . . . . . . . . . . . . . . . . . . . . . 3 A.1 Rounded coefficients of the hydrostatic mapping function [50]. . . . . . . I A.2 Rounded coefficients of the wet mapping function [50]. . . . . . . . . . . II xxi List of Tables xxii 1 Introduction Global Navigation Satellite Systems (GNSS) have become an indispensable part of mod- ern life, providing accurate location information for a myriad of applications ranging from navigation and mapping to timing and synchronization. The backbone of GNSS sys- tems includes satellite constellations such as the U.S. Global Positioning System (GPS), the Russian Federation Global Navigation Satellite System (GLONASS), the European Galileo system, and the Chinese BeiDou Navigation Satellite System (BDS) [3]. These systems are mostly global, meaning they provide coverage on large parts of the Earth, in contrast to the Japanese system Quasi-Zenith Satellite System (QZSS) and India’s Navigation with Indian Constellation (NavIC), which are regional systems. At least four satellites are required to determine the receiver position, in X, Y and Z direc- tions, and the receiver clock bias. If more variables are to be determined, more satellites are needed to accurately estimate the position and by using multiple constellations, the number of available satellites increases [4]. This process is called multi-GNSS. The satel- lites transmit signals containing various data to receivers on the ground, which can be used to determine the receiver’s location. Calculations are performed to extract the dis- tance to each satellite and, through trilateration, the position of the receiver. During the transmission from satellite to receiver, the signals are affected by many disturbances, which affect the accuracy of the estimated receiver position. Some of these error sources and corresponding solutions will be discussed in the following sections. There are relative and precise positioning methods, with the difference being how many receivers are needed in order to estimate the position accurately. The relative positioning depends on two or more receivers, one or more stationary base stations and one mobile rover. The rover position is estimated using complementary data transmitted by the base station. This method can produce centimetre-level positioning accuracy but is limited by the range between the two receivers [5]. On the other hand, precise positioning only uses one receiver but has a longer convergence time, i.e., the time it takes until an estimate is calculated. This can be problematic for demanding systems such as automated vehicles [4]. This work will aim to develop and implement a precise point positioning (PPP) method, where a deep understanding of multiple error sources is needed. The next Chapter, 2, contains theory regarding the GNSS signals and their contents, relevant error sources and how to evaluate GNSS positioning. In Chapter 3, relevant literature is presented along with techniques to mitigate error sources and how to solve the carrier phase ambiguities. The methodology of the project is also described. Chapter 4 includes the results, which are then discussed in Chapter 5 and the conclusions are 1 1. Introduction presented in Chapter 6. 1.1 Research questions • How can the hardware components of a multi-constellation GNSS receiver be opti- mized for cost-effectiveness without compromising accuracy? • What strategies can be implemented to mitigate common sources of error in multi- constellation GNSS positioning, particularly in low-cost solutions, to enhance over- all accuracy? • How does the developed GNSS solution compare in terms of accuracy, reliability, and performance against both proprietary and open-source alternatives under var- ious environmental conditions? 1.2 Methodology As the research questions describe, this project aims to develop a low-cost GNSS solution by combining techniques for mitigating error sources such as biases and atmospheric disturbances. To achieve this goal, a comprehensive literature study will be performed to gather in- formation regarding GNSS theory. The focus will be on deepening the knowledge of the possible sources of errors for absolute positioning systems and what corrections must be performed to achieve high accuracy. Time will be spent delving into the area of PPP by going through similar implementations and finding areas of development that can be improved in this thesis. Research on available hardware will be conducted to find a receiver that can collect raw data, which will be used to test the implementation. The receiver also has to be able to receive two frequencies while still being a low-cost device. 1.3 Limitations / Demarcations The limitations of this thesis include how the data is collected and handled and what constellations will be used. To validate and test the implementation, GNSS data have to be used. Only data collected by the chosen hardware will be used. A stationary receiver will acquire this data, meaning the implementation will be a static positioning solution. Additionally, all the calculations for the corrections and the position estimation will be performed in post-processing, not in real-time. The project will utilise a dual-constellation implementation, using data from GPS and Galileo satellites. 2 2 Theory This section initially provides a general overview of GNSS positioning and explains the steps required to determine a receiver’s position. Then, it delves into the GNSS signals and observables, along with various coordinate systems relevant to GNSS positioning. The section also discusses common sources of error and how to manage them. Finally, it presents different methods of positioning (relative and absolute) and explains how to evaluate the accuracy of the estimations provided by these methods. 2.1 GNSS signals This section describes the structure of a GNSS signal and its three main components, the carrier wave, the ranging code and the navigation data. 2.1.1 Radio Frequency (RF) Carrier The satellite transmitter produces a sinusoidal wave with a constant frequency called the radio frequency carrier or carrier wave [3]. The carrier wave frequency ranges from 1 to 2GHz, depending on the constellation. The Institute of Electrical and Electronics Engineers (IEEE) has defined a letter band designation to categorize ranges of frequencies for specific purposes [6], as shown in table 2.1. The GNSS signals’ carrier wave frequencies fall within the L-band range. Table 2.1: Radar-frequency letter bands. Band designation Nominal frequency range HF 3 MHz to 30 MHz VHF 30 MHz to 300 MHz UHF 300 MHz to 1000 MHz L 1 GHz to 2 GHz S 2 GHz to 4 GHz C 4 GHz to 8 GHz X 8 GHz to 12 GHz Ku 12 GHz to 18 GHz K 18 GHz to 27 GHz Ka 27 GHz to 40 GHz V 40 GHz to 75 GHz W 75 GHz to 110 GHz 3 2. Theory The L-band is divided into three main frequencies for transmitting GPS signals: L1 (1575 MHz), L2 (1227 MHz), and L5 (1176 MHz). For the Galileo constellation, the main frequencies are E1 (1575 MHz), E5a (1176 MHz) and E5b (1207 MHz). The frequencies of E1 and E5a correspond to those of L1 and L5 of the GPS signal and will henceforth be referred to as L1 and L5. Each of these signals has unique characteristics. For instance, [7] states that L1 has the lowest ionospheric refraction error, L2 performs best in cross-correlation, and L5 has the highest power. Further, the L5 signal is described as ideal for high-precision positioning and is designed for safety-of-life applications due to its high power. It is also mentioned that L5 has a lower standard deviation when estimating position. The radio frequency carrier alone does not provide any information, but adding data with modulation is possible. Modulation is the process of adding information to a carrier wave by converting the data, in the form of symbols, to a carrier waveform [8]. An example of this operation is depicted in Figure 2.1. The figure includes the carrier wave, a code signal, or ranging code, and the resulting signal when modulating the two. The ranging code is explained in the next section. Figure 2.1: Ranging code, carrier wave and the modulated signal. 2.1.2 Pseudo-Random Noise (PRN) ranging codes Ranging codes consist of binary sequences that enable the receiver to calculate the signal’s travel time from the satellite to the receiver [4]. GNSS ranging codes, that are modulated onto the carrier, are referred to as pseudo-random as they do not have an apparent pattern. However, they repeat after a certain duration [9]. 4 2. Theory 2.1.3 Navigation data The navigation data contains information such as the satellite ephemeris, atomic clock data, corrections and the almanac data [4]. The satellite ephemeris consists of the nec- essary parameters to calculate the position of the satellite at a given point in time, and the almanac consists of information about all the satellites in a constellation including health status and rough orbit information [4]. This message is crucial for determining the receiver’s location by providing vital information about the satellites. 2.2 GNSS Observables GNSS observables are the measurements a receiver computes from the signals it receives from a satellite. There are two main types of observables, the pseudorange and the carrier phase. 2.2.1 Pseudorange The pseudorange, also referred to as code phase ranging or code measurement, is the distance between a satellite and receiver calculated using the signal travel time retrieved from the ranging code signal. If no disturbances or errors are present, the pseudorange is equal to the geometric distance ρ (actual distance) described in Equation (2.1) for satellite s and receiver r [10]. ρs r = c(tr − ts) = c∆t = √ (xs − xr)2 + (ys − yr)2 + (zs − zr)2 (2.1) Where • c is the speed of light in vacuum • tr is the time when the signal is received at the receiver • ts is the time when the signal is transmitted from the satellite • ∆t is the transmission time of the signal. Since there are errors and disturbances within practical scenarios, the pseudorange is not equal to the geometric distance and can be written as in Equation (2.2) for the band j. P s j = c ( t̄r − t̄s ) t̄r = tr + dtr t̄s = ts + dts (2.2) 5 2. Theory Where • P s j is the pseudorange observable • dtr is the receiver clock bias • dts is the satellite clock bias. Equation (2.2) can be further developed by introducing the geometric distance and the most prominent disturbances as described in Equation (2.3) [4]. P s j = ρs r + c (dtr − dts) + ∆rel + br − bs j + T s − gjI + ε(P s j ) (2.3) Where • ρs r is the geometric (true) distance from the antenna phase centre (APC) of the satellite to the APC of the receiver • ∆rel is the special relativistic effect • br is the code receiver instrumental delay • bs j is the code satellite instrumental delay • T s is the tropospheric delay • I is the ionospheric delay • gj is the frequency-specific factor that affects the ionospheric delay • ε(P s j ) is the residual term that denotes the unmodeled terms. 2.2.2 Carrier Phase Measurement The carrier phase measurement measures the beat frequency, which is the difference in frequency between the signal generated by the receiver’s oscillator and the signal transmitted from the satellite [11]. The beat frequency technique is utilized to produce an ambiguous measurement of the distance to the satellite in cycles, offset by an ambiguous number of whole wavelengths. In contrast to the pseudorange method, the carrier phase measurement is less influenced by noise and can theoretically provide accuracy within a couple of millimetres [11]. The carrier phase measurement is written in Equation (2.5) in the same manner as the pseudorange measurement. λjL s j = Φs j (2.4) Φs j = ρs r + c(dtr − dts) + ∆rel +Br −Bs j + T s − gjI + λj(N s j + ws) + ε(Φs j) (2.5) 6 2. Theory Where • Ls j is the carrier phase measurement in cycles • λj is the wavelength of the signal • Φs j is the carrier phase measurement in meters • Br is the code receiver instrumental delay • Bs j is the code satellite instrumental delay • w is the carrier phase windup effect • N s j is the integer ambiguity. When the receiver locks on to the signal, it can only determine the signal’s phase shift and frequency change and not the number of whole cycles between the satellite and receiver, resulting in a carrier phase ambiguity (CPA) that needs to be estimated in order to use the carrier phase measurement [12]. 2.2.3 Combination of signals When combining GNSS observables from different frequencies of the same satellite, new measurements with updated characteristics can be computed. The following combinations use pairs of signals, i.e., dual-frequency measurements. 2.2.3.1 Ionosphere-free (IF) The IF measurements are calculated according to Equations (2.6) and (2.7) [4]. ΦIF = f 2 1Φ1 − f 2 5Φ5 f 2 1 − f 2 5 (2.6) PIF = f 2 1P1 − f 2 5P5 f 2 1 − f 2 5 (2.7) This combination is free from up to 99.9% of the ionospheric disturbance [4]. The mea- surement obtained from the IF combination has a higher level of noise compared to the individual signals [13]. The noise will be amplified with a scale factor that is dependent on the squared wavelength ratios µj = λ2 j/λ 2 1, for band j [14]. The scale factor is calculated using Equation (2.8). √√√√ µ2 1 + µ2 5 (µ5 − µ1)2 ≈ 2.59 (2.8) 2.2.3.2 Geometry-free (GF) As the name suggests, the GF combination is free from the geometric parts of the mea- surement, such as the clocks and non-frequency dependent effects [4]. The combined 7 2. Theory measurements are described by Equations (2.9) and (2.10). ΦGF = Φ1 − Φ5 (2.9) PGF = P5 − P1 (2.10) 2.2.3.3 Narrow-lane (NL) Using Equations (2.11) and (2.12) NL measurements are calculated that have a particu- larly short wavelength and less noise [4]. ΦNL = f1Φ1 + f5Φ5 f1 + f5 (2.11) PNL = f1P1 + f5P5 f1 + f5 (2.12) 2.2.3.4 Wide-lane (WL) Unlike the NL combination, the WL measurement has an especially large wavelength, which is useful for ambiguity resolution [4]. The calculation of the WL combination is shown in Equations (2.13) and (2.14). ΦW L = f1Φ1 − f5Φ5 f1 − f5 (2.13) PW L = f1P1 − f5P5 f1 − f5 (2.14) 2.2.3.5 Melbourne-Wübbena (MW) The ionospheric dependence can be removed from the ΦW L and PNL measurements by combining them according to the MW combination, which is calculated with Equation (2.15). Further benefits with this combination include a longer wavelength than the original measurements due to the WL measurement and reduced measurement noise thanks to the NL measurement [4]. MW = ΦW L − PNL (2.15) 2.3 Coordinate systems Coordinate systems are fundamental for understanding the position of points in space. They provide a standardized way of describing a point’s location in a given space. Some of the relevant coordinate systems for this thesis are described below. 8 2. Theory 2.3.1 Conventional Terrestrial Reference System (TRS) The Conventional TRS, also called Earth-centred Earth-fixed (ECEF), is a spatial refer- ence system, meaning that different realisations and versions exist. All versions, however, have the origin at the centre of an estimated ellipsoid (representing the Earth), with the Z-axis pointing from the centre to the conventional terrestrial pole, the X-axis pointing through the cross-section of the prime meridian and the equator, and the Y-axis orthogo- nal to the Z and X axes forming a right-handed system [4]. The ECEF frame is depicted in Figure 2.2. Z Y X North East Up ecef ecef ecef φ λ Figure 2.2: ECEF and ENU frames [1]. 2.3.1.1 International Terrestrial Reference Frame (ITRF) The ITRF is a realization of the ECEF produced by the International Earth Rotation and Reference Systems Service (IERS). It is a highly accurate and stable coordinate system in geodesy and Earth science for global reference [4]. 2.3.1.2 East, North, Up (ENU) The ENU coordinates are horizontal coordinates based on geographical location and are thus local. The ’East’ and ’North’ represent directions on the Earth’s surface, while ’Up’ points towards the sky, perpendicular to the Earth at the location. This coordinate system is widely used to describe horizontal errors in GNSS positioning [4]. The ENU frame is depicted together with the TRS frame in Figure 2.2. 2.3.1.3 Geodetic Coordinates in an ECEF frame can be expressed as longitude, latitude and height above the reference ellipsoid. In Figure 2.2, the latitude and longitude components can be seen as φ and λ, respectively. 9 2. Theory 2.4 GNSS Error Sources GNSS positioning suffers from several error sources, which will be covered in this section. 2.4.1 Emission and Reception Time In order to accurately determine a receiver’s position using GNSS algorithms, it is nec- essary to know the precise coordinates of the satellites at the time of signal emission. During the signal transmission time, i.e. the time it takes for the signal to travel from the satellite to the receiver, the satellites travel about 300 meters [4]. Without consider- ing the transmission time, disturbances of up to 50 meters in the horizontal and vertical direction can occur [4]. To handle this error, the emission time has to be determined to calculate the satellite coordinates. The emission time can be calculated according to Equation (2.17) using the pseudorange measurement. ∆ttp = P c (2.16) ts[emission] = tr[reception] − ∆ttp (2.17) Where • ∆ttp is the travel time of the signal calculated using the pseudorange • ts[emission] is the time of signal emission from the satellite • tr[reception] is the time of signal reception in the receiver. 2.4.2 Earth Rotation As the signal travels from the satellite, the earth will rotate, meaning there will be a difference between the receiver position at the time of emission and reception if an ECEF frame is used. This will affect the estimated position, adding an offset of about 25 meters to the east. Equation (2.21) corrects the satellite position regarding this effect [4]. R(θ) =  cos(θ) sin(θ) 0 − sin(θ) cos(θ) 0 0 0 1  (2.18) ∆ttg = ∥rs − rr∥ c (2.19) δϑ = ωe∆ttg (2.20) r̃s = R(δϑ) · rs (2.21) 10 2. Theory Where • R is the rotation matrix of the earth around the z-axis • ωe is the rotation speed of the earth • δϑ is the rotation angle of the earth during the travel time of the signal • ∆ttg is the travel time of the signal using the geometric distance to the satellite • rs is the position in the ECEF frame of the satellite • r̃s is the corrected position of the satellite adjusted by the earth’s rotation. 2.4.3 Relativistic Effects According to Einstein’s theory of special relativity, two clocks travelling at different speeds will experience time dilation, meaning that the clocks will advance at different speeds. The faster an object moves relative to an observer, the slower it will appear to tick for the observer. Since the satellites travel at about 5km/s, this effect will be noticeable [4]. Another relativistic effect that affects GNSS positioning is from Einstein’s theory of general relativity, which predicts the time dilation of two clocks placed at different distances from a gravitational mass. The closer a clock is to a gravitational mass relative to an observer, the slower the clock will appear to tick. These effects must be accounted for since they can affect the position by tens of metres [4]. The special relativistic effect, which has the largest influence, can be cancelled with Equation (2.22), which corrects the GNSS observables in Equations (2.3) and (2.5). ∆rel = −2rs · vs c (2.22) Where • vs is the velocity in the ECEF frame of the satellite. Applying Equation (2.23) corrects the general relativistic effect. However, this only affects the positioning by a couple of centimetres. ∆ρrel = 2 ·GME c2 ln ( rs + rr + rs r rs + rr − rs r ) (2.23) 11 2. Theory Where • GME is Earth’s gravitational constant • rs is the geocentric distance to the satellite • rr is the geocentric distance to the receiver • rs r is the distance from the receiver to the satellite. 2.4.4 Cycle Slip Cycle slips are sudden changes in the estimated number of cycles related to the carrier phase measurement due to a temporary signal loss between the receiver and satellite [4]. This loss of lock can be caused by a weak connection or an obstacle, such as a tall building. Cycle slips can be discovered by either finding a data gap, i.e. a loss of data from a satellite, or an unexpected behaviour in the signal. There are multiple ways of detecting cycle slips, depending on how many frequencies are accessed. For a dual-frequency system, the geometry-free combination of the carrier phase measurements, shown in Equation (2.24) can be used as input data. Φs I(k) = Φs 1(k) − Φs 5(k) (2.24) Where • Φs I is the geometry-free carrier phase measurement • k is the epoch. Unexpected behaviour in the signal can be detected by fitting a second-degree polynomial to the previous data samples. A cycle slip can be declared if a specified threshold is exceeded by comparing the predicted value of the next epoch with the actual value from the input data. Further, a cycle slip is declared when there is a data gap. 2.4.5 Antenna Phase Centre (APC) The APC represents the radiation source, which does not have to align with the physical centre of the antenna [4]. Instead, the APC depends on the frequency and angle of the signal, as depicted in Figure 2.3. The difference between the physical centre, or antenna reference point, and the APC is called the Antenna Phase Offset (APO) and is also shown in the figure. This phenomenon is present in both the satellite and receiver antennas. 12 2. Theory Antenna Antenna Phase Center Antenna Phase Offset Antenna Reference Point Figure 2.3: Simplified depiction of a receiver antenna with the corresponding antenna phase centre, antenna phase offset and antenna reference point. Some corrections for the related error sources, such as clock and orbit errors, refer to the satellite mass centre, not the APC. The offset between the two points, ∆AP C , is in a satellite-fixed coordinate frame, which has to be recalculated to the ECEF reference frame. The recalculation is based on the unit vectors (̂i, ĵ, k̂) that define the satellite-fixed coordinates. These unit vectors are defined in Equations (2.25) to (2.27). k̂ = − rsMC ∥rsMC ∥ (2.25) ĵ = k̂ × rsun − rsMC ∥rsun − rsMC ∥ (2.26) î = ĵ × k̂ (2.27) Where • k̂ is the vector pointing to the centre of Earth from the mass centre • ĵ is the cross-product between k̂ and the unit vector from the satellite to the sun • î is the cross-product between ĵ and k̂. The satellite APC in ECEF coordinates is calculated with Equation (2.28). rsAP C = rsMC + R · ∆AP C (2.28) Where R = [̂i ĵ k̂]. 2.4.6 Atmospheric Effects When a satellite broadcasts a GNSS signal to a receiver on Earth, the signal passes through the Earth’s atmosphere. During this process, the signal is refracted and delayed, which causes a slight change in direction and impacts the signal’s speed [15]. These changes result in errors when determining the receiver position and thus must be handled. 13 2. Theory In Figure 2.4, the different layers of the atmosphere are depicted with their different heights, temperature and electron density. The primary source of error arises as the signal traverses through the ionosphere and the troposphere. However, the reasons for these errors differ and will be elaborated upon in the following subsections [15]. Figure 2.4: Atmospheric layers with their height, temperature and electron density [2]. 2.4.6.1 Ionosphere The main reason for the instabilities in the GNSS signal is the Earth’s ionosphere [16]. The ionosphere is a layer in the Earth’s atmosphere that gets ionized due to the sun’s ultraviolet radiation. This ionization process releases free electrons, causing the atmo- sphere’s density to fluctuate. This fluctuation is the primary cause of signal delay [16]. Figure 2.4 depicts the general electron density for different altitudes. The ionosphere’s Total Electron Content (TEC) varies throughout the day and year, making it a challeng- ing variable to estimate accurately and may cause errors [16]. The ionosphere is what is called a dispersive medium, meaning that RF signals of different frequencies will propagate through it at different speeds. This property makes it possible to estimate the ionosphere’s electron content, which could be used when estimating a position [4]. 2.4.6.2 Troposphere The troposphere, the closest atmospheric region to Earth, as seen in Figure 2.4, affects the signal locally, meaning it is hard to use estimations from stations far away [4]. The tropospheric delay, T (ϵ), consists of two parts and can be written as: T (ϵ) = ZHD ·mh(ϵ) + ZWD ·mw(ϵ) (2.29) 14 2. Theory Where • ZHD is the zenith value for the hydrostatic (dry) delay • ZWD is the zenith value for the wet delay • mh is the mapping function for the dry delay • mw is the mapping function for the wet delay. The mapping functions mh and mw depend on how the troposphere is modelled. Figure 2.5 depicts a simplification of the tropospheric delay where the purpose of a mapping function would be to estimate the propagation of the signal through the troposphere. Figure 2.5: Simplification of tropospheric delay. Contrary to the ionosphere, the troposphere is a non-dispersive medium for GNSS signals, meaning that the signal’s frequency does not affect its travel speed [4]. Thus, the same corrections must be applied to the code and phase at different frequencies. 2.4.7 Multipath A multipath error occurs when a signal sent from one satellite is received multiple times at a receiver, which causes interference [4]. This phenomenon can happen when the receiver is situated near large buildings or when the satellite is at a low elevation angle, causing the signal to bounce off the walls or the ground and reach the receiver from multiple directions, as illustrated in Figure 2.6. Figure 2.6: Illustration of multipath effect on GNSS signal. 15 2. Theory For the code signal, multipath can theoretically cause errors as large as 450m, but the error is usually closer to two meters [4]. For the phase measurement, the error is commonly less than 1 cm. 2.4.8 Receiver Noise Receiver noise arises from either the hardware or the software in the receiver. Generally, the error due to receiver noise is approximately 1% of the signal’s wavelength [17]. 2.4.9 Clock Error To estimate a position using satellites, a very precise clock is necessary at the satellite. Therefore, satellites have atomic clocks with a precision of ∆f/f ≃ 10−13 [4]. Even though the satellite clocks are very accurate, they tend to drift slightly. Even a small drift can affect the estimated position greatly when it is calculated by the receiver [16]. These errors can be corrected either by differencing between two receivers or by using clock corrections from compute stations [4]. 2.4.10 Orbit Error Similar to the satellite clock, it is necessary to know the position of the satellite to estimate a receiver’s position. The satellites transmit parameters such that it is possible to calculate their position, but in most cases, it is not good enough for precise positioning [4]. Therefore, just like the satellite clocks, corrections are needed to handle these errors. 2.4.11 Hardware-induced biases GNSS hardware has certain deficiencies that cause biases due to minor delays in the signal transmission from a satellite or in signal reception by a GNSS receiver [18]. As a result of these delays, events that should co-occur may become asynchronous. These biases can vary between code signals (code bias) and carrier waves (known as phase bias). If these effects are not dealt with properly, they can impair the estimation of the receiver’s position, extend the convergence time, and decrease the chance of integer ambiguity resolution. 2.4.11.1 Differential Code Bias (DCB) DCBs are the time delays between two GNSS signals transmitted from a single satellite [18]. DCBs are expressed as in Equation (2.30) for the hardware delays between band i and j for satellite s. DCBs are often used as corrections retrieved from a computation centre when working with combinations of different frequencies. DCBs i,j = bs i − bs j (2.30) 2.4.12 Tidal Deformations Due to the gravitational forces that the Earth is subject to, the crust deforms, thus moving the receiver’s position. To mitigate this issue, modelling the displacement of 16 2. Theory the Earth’s crust is necessary. This displacement can be divided into solid tides, ocean loading, and pole tides. Solid tides cause the most significant offsets, hence, only these will be described. 2.4.12.1 Solid Tides Solid tides refer to the movement of the Earth’s crust due to external gravitational forces from the Sun and the Moon [4]. The displacement can be modelled as Equation (2.31) where ∆rsol needs to be applied to the receiver position [19]. ∆rsol = 3∑ j=2 GMj R 4 E GME R3 j ( h2 r̂ (3 2(R̂j · r̂)2 − 1 2 ) + 3l2 (R̂j · r̂) ( R̂j − (R̂j · r̂)r̂ )) (2.31) Where • ∆rsol is the displacement of the receiver • GMj is the gravitational constant of the Moon and the Sun, respectively • GME is the gravitational constant of the Earth • RE is the radius of the Earth • R̂j is the unit vector from the Earth to the Moon and Sun in ECEF, respectively • Rj is the magnitude of the Moon and Sun position vectors • r is the magnitude of r • h2 is the nominal degree 2 Lovenumber • l2 is the nominal degree 2 Shilda number. 2.5 Extended Kalman Filter The Extended Kalman Filter (EKF) is often used to estimate position in GNSS appli- cations [20]. The EKF is based on the process equation in Equation (2.32) and the observation equation in Equation (2.33). The process equation describes how the states are expected to traverse from one time step to the next, with a known motion model, while the observation model describes the relationship between measurements and the state vector. Xk = fk(Xk−1) + νk (2.32) Zk = hk(Xk) + ωk (2.33) 17 2. Theory Where • fk is the motion model • hk is the measurement model • νk ∼ N (0,Qk) • ωk ∼ N (0,Rk). The EKF uses the previous equations to form prediction and update steps. The prediction step uses a known motion model of the states and the previous state estimates to form the predicted (a priori) estimate. The prediction step also calculates a predicted covariance matrix P that describes the covariance of the predicted estimates. The prediction step is shown in Equations (2.34) and (2.35). X̂k|k−1 = AX̂k−1|k−1 (2.34) Pk|k−1 = APk−1|k−1A + Q (2.35) Where • Ak = ∂fk ∂Xk−1|k−1 . The update step uses the knowledge from Equation (2.33) to form a posterior update, meaning that it uses information from all time steps, including the current step, to form an updated estimation of the state vector and covariance matrix. The EKF uses a Kalman gain calculated based on the innovation covariance in Equation (2.37), observation model and predicted covariance. The Kalman gain dictates how certain measurements and pre- dictions should be weighted based on covariance. The update step is shown in Equations (2.36)-(2.40). vk = Zk − h(X̂k|k−1) (2.36) Sk = HkPk|k−1H T + Rk (2.37) Kk = Pk|k−1H T k S −1 k (2.38) X̂k|k = X̂k|k−1 + Kkvk (2.39) Pk|k = (I − KkHk)Pk|k−1 (2.40) Where • Hk = ∂hk ∂Xk|k−1 . 18 2. Theory 2.6 Satellite-Based Positioning There are two main ways of constructing a GNSS solution, utilizing relative or absolute positioning. 2.6.1 Relative Positioning Relative positioning involves the process of determining the location of a receiver (rover) by utilizing one or more stationary receivers, known as base stations, with pre-established positions. 2.6.1.1 Differential GNSS (DGNSS) DGNSS is a relative positioning method using a stationary base station with a known location and a mobile rover, both of which act as receivers for the signals transmitted by GNSS satellites [21], as shown in Figure 2.7. The classical DGNSS method is based on calculating the pseudorange, see Section 2.2.1, for the receivers. Since the location of the base station is known, corrections can be calculated and transmitted to the rover. The pseudorange for the rover is calculated and the base station corrections are applied to improve the process of determining the position. This method can achieve a one-meter accuracy when the GPS constellation is used, provided the user is within a certain range, a few tens of kilometres from the reference station. Figure 2.7: Simple depiction of relative positioning. 2.6.1.2 Real-Time Kinematic (RTK) RTK orginated in the mid 1990s and is a type of DGNSS method that uses carrier-based ranging, meaning it can provide accurate positioning information but has to handle the ambiguities related to carrier phase measurements described in Section 2.2.2 [16], [4]. The RTK system can provide centimetre-level accuracy and includes a static base station and a mobile rover station, as shown in Figure 2.7. 19 2. Theory Correction data, including atmospheric corrections, can be calculated and sent to the rover station along with the distance between the two, improving the estimation of the rover’s position [16]. To access the corrections from a base station the rover must be in range and with minimal obstructions, as they reduce accuracy. 2.6.1.3 Network RTK For RTK users who do not have their own base station, there is an alternative to buy a service called Network RTK. Network RTK uses multiple stationary base stations which send information regarding their positioning to a control centre [16]. The rover data is also sent to the central centre, where it is corrected based on the data provided by the network of base stations. This method is dependent on two-way communication between the centre and the rover and is depicted in Figure 2.8. The communication can be performed via Networked Transport of RTCM via Internet Protocol (NTRIP) if the receiver and the central centre have Internet access. Figure 2.8: Simple depiction of Network RTK. 2.6.1.4 Wide Area RTK (WARTK) WARTK is a concept to solve the range problem in RTK by providing the user with computations of an accurate ionosphere model to correct for the ionospheric errors, thus providing an expected range of hundreds of kilometres from the nearest base station. This is, however, only a concept, and no operational WARTK centres exist at this time [22]. 2.6.2 Precise Point Positioning (PPP) PPP is a technique used as an alternative to RTK positioning, where only one receiver is used [23]. It is not possible to differentiate between receivers, as only one is available. Unlike RTK, PPP does not require a base station to receive corrections. Instead, PPP uses a network of global reference stations that provide necessary satellite clocks and 20 2. Theory orbital corrections (and biases), see Figure 2.9. The corrections from the global reference stations are then transmitted to the user either over satellite (L-band) or via the internet [16]. PPP uses both the pseudorange and carrier phase measurement to achieve high accuracy. PPP was introduced in 1997 by Zumberge, J et al. and has since continued to rise in popularity because of the lack of base station and freely available data [24],[25]. Figure 2.9: Precise positioning system. There are different versions of PPP whose differences are discussed in the sections below. 2.6.2.1 Standard PPP ”Standard PPP” is a method that uses dual-frequency observables to form an IF combi- nation, which eliminates almost all ionospheric disturbances [4]. This method does not fix the ambiguities to an integer number, instead, they will be floats. Float ambiguity is when the ambiguity is estimated together with the code biases and phase biases of the receiver and satellite, resulting in a float. This approach suffers from a long convergence time, as the ambiguities are not fixed [23]. The corrections used in this method can only be retrieved from the internet and are not available in real-time. Hence, this method is used for post-processing. 2.6.2.2 Ambiguity-fixed PPP PPP with ambiguity resolution (or PPP-AR) is when an IF combination is used, and the ambiguities are fixed to integers, which can improve convergence time and accuracy [23]. Multiple methods for computing the integer ambiguity exist, one of which is described in the Method, see Section 3. The problem with both float and fixed ambiguity is that when an IF estimate is used, the noise increases substantially, making it harder to estimate a position. This is why these methods require a longer initialization time [26]. As this method uses the same corrections as the standard PPP, it is also only applicable as a post-processing method. 21 2. Theory 2.6.2.3 PPP-RTK One disadvantage of standard PPP and PPP-AR is the long convergence time, however, recent studies have shown that a convergence time of a few seconds can deliver centimetre- level precision with the use of modern State Space Representation (SSR) corrections [5]. The modern SSR corrections do not only include the clock and orbital errors but also atmospheric effects, making it possible to not rely solely on the IF measurement, which has a substantial effect on noise levels [26]. Another advantage is that the SSR corrections are available in real-time and can be received via both satellite and the internet, making it the only solution possible to use in real-time systems. 2.7 Evaluation of GNSS Positioning GNSS positioning can be evaluated in various ways. It is possible to calculate both an estimated accuracy, which shows what accuracy should be expected, and an actual accuracy, which shows how well the algorithm performs. 2.7.1 Vertical and horizontal accuracy In GNSS, a receiver’s positional accuracy is determined by its horizontal and vertical precision. Typically, the focus is on horizontal precision, which refers to the accuracy of the positioning in latitude and longitude. Conversely, vertical positioning indicates the accuracy of the height plane. 2.7.2 Dilution Of Precision (DOP) DOP, or Dilution Of Precision, is a measure of the expected accuracy of a position calculated by a receiver [4]. The DOP value is derived solely from the estimated position of the receiver and the positions of the satellites. A high DOP value is observed when the tracked satellites are in close proximity to each other, making it challenging to estimate a position accurately. Conversely, a low DOP value is achieved when the satellites are widely dispersed, facilitating a more precise position estimation. Therefore, a lower DOP value (<5) is preferable as it yields a more accurate estimate [27]. There are several different DOP measurements, but the ones relevant to this thesis are the Horizontal DOP (HDOP) and the Vertical DOP (VDOP). To calculate the DOPs, the Line Of Sight (LOS) matrix has to be defined as in Equation (2.41), which includes the line of sight vector to every satellite n. LOS =  sin az1 cos el1 cos az1 cos el1 sin el1 1 sin az2 cos el2 cos az2 cos el2 sin el2 1 ... ... ... ... sin azn cos eln cos azn cos eln sin eln 1  (2.41) 22 2. Theory Where • az is the azimuth angle to the satellite • el is the elevation angle to the satellite. The DOP matrix is designed as the covariance matrix of the LOS matrix as in Equation (2.42). Q = Cov(LOS) = (LOST LOS)−1 =  qxx qxy qxz qxt qxy qyy qyz qyt qxz qyz qzz qzt qxt qyt qzt qtt  (2.42) QECEF = qxx qxy qxz qxy qyy qyz qxz qyz qzz  (2.43) To form the QENU matrix, which is used for VDOP and HDOP, QECEF from Equation (2.43), has to be transformed as shown in (2.44). QENU = RT QECEF R = qee qen qeu qen qnn qnu qeu qnu quu  (2.44) Where • R is the transformation matrix from ECEF to ENU at a given latitude and longi- tude. The HDOP and VDOP can then be calculated according to Equations (2.45) and (2.46). HDOP = √ qee + qnn (2.45) VDOP = √ quu (2.46) 23 2. Theory 24 3 Methods This chapter provides an overview of the hardware used for collecting GNSS data. It also discusses the selection of positioning methods and the frequencies utilized. Additionally, it details the corrections applied and how they are implemented. The chapter includes information on the collection of data used for validating the implementation. It further describes the use of the Extended Kalman Filter for estimating receiver positions and the process for resolving ambiguities. The implementation is further outlined and visualized through flowcharts, followed by a comparison with the GNSS processing tool RTKLIB. 3.1 Hardware One essential component of developing a complete GNSS solution is the receiver. For this project, the hardware was required to include an antenna for GNSS signals and a positioning module that could receive the raw GNSS data. The hardware was assessed based on the price, the frequencies it could receive, whether it could output raw GNSS data, the complexity of the setup process, and the performance. One receiver that met all the requirements and performed well in all the categories was the evaluation kit EVK-F9P from u-blox [28]. The evaluation kit included an ANN-MB antenna and a box containing a printed circuit board, both of which are shown in Figure 3.1. Mounted on the PCB was the u-blox ZED-F9P chip, a high-precision positioning module. (a) EVK-F9P box (b) ANN-MB antenna Figure 3.1: Receiver hardware. 25 3. Methods The EVK-F9P was priced at around 3000 SEK, making it a low-cost module. The inclusion of an antenna in the kit made it cost-efficient and a good option to consider when taking the price into account. As the kit was complete, no further parts had to be purchased. The ZED-F9P module is a dual-frequency receiver and can receive either L1/L2 or L1/L5 frequencies, depending on the model, EVK-F9P-01 or EVK-F9P-16 respectively. In both cases, it can output raw GNSS data, which was essential to this project. Being a dual- frequency receiver, it allows for combining the signals using the combinations described in Section 2.2.3. The model used in this thesis was the EVK-F9P-16, and the reasoning behind this decision is discussed further in Section 3.3. The setup process was fairly fast and simple, partly due to the inclusion of an antenna, a straightforward connection with a USB-C cable and the u-blox GNSS evaluation software u-centre. As the kit was assembled upon delivery, the only setup necessary was connecting the antenna to the box, which was then connected to a computer where the evaluation software was downloaded and used to analyse the data. Some configuration was required to extract only the raw data from the module, which is explained in Appendix A.3. The ZED-F9P module has been tested and proven to perform well in static PPP-RTK with L1/L2 frequencies [29], [30]. Good conditions resulted in high signal-to-noise (SNR) values and low multipath disturbances. With this module, it is possible to achieve high precision, with an accuracy of centimetres. 3.2 Positioning Method The choice of positioning method for this project heavily depended on the application and the most recent advancements in the field, as this project aimed to drive further innovation in the GNSS area. RTK, which is introduced in the Theory Section 2.6.1.2, has been well explored. It relies on differencing between two receivers, eliminating the need for advanced corrections. The major drawback of RTK is the need for close proximity to a stationary base station with known coordinates. This means that if a rover requires a precise position in a large area (generally larger than a 10km radius), the performance of traditional RTK will greatly decrease due to atmospheric differences between the receivers. Sending corrections from a base station to a rover generally also demands clear LOS if RF is used, which means there can be no obstructions between the two. The corrections could, however, be sent via NTRIP, which would require Internet access for both the rover and the base. Although PPP was first introduced in 1997, advancements are still being made. As discussed in the Theory Section, PPP was a post-processing method until recent years, when the corrections started to be readily available via NTRIP. These corrections are provided for free in real-time in the SSR format. PPP also has better scalability compared to network-RTK since no data has to be sent from the receiver to the reference stations, thus making it possible to have many users relying on the same corrections distributor. In the future, if self-driving cars become more common, this attribute will become more important since a lot of receivers demanding precise positions will be operational in tight areas, thus needing the same corrections. 26 3. Methods PPP in real-time with SSR corrections (PPP-RTK) has been shown to produce results with a high accuracy and low convergence time, making it feasible to use in kinematic scenarios like in a car [5]. In the specified source, an accuracy of 0.1m in just 2 seconds is achieved 90% of the time in a kinematic scenario. Since this thesis aimed to explore newer and promising areas in GNSS positioning, the focus lay on PPP. However, with the simplification of using precise products available for post-processing, which were used to simulate real-time corrections with a real-time algorithm. An ambiguity-fixing algorithm was researched and developed to fix ambiguities in a short time frame. Several sources describe that real-time ambiguity fixing is possible with proper corrections and can generate centimetre to decimetre accuracy [5], [23]. 3.3 GNSS Frequencies Since a two-frequency receiver was chosen, the choice of frequencies was between a combi- nation of L1/L2 and L1/L5. L1/L2 would have been the most straightforward one since it is the most common combination, most corrections even refer to the ionosphere-free combination of L1/L2 [31]. However, the L5 signal, which is the most modern and should right now only be used at its own risk, is designed to be a safety-of-life signal with greater bandwidth to improve jam and multipath resistance [32]. As the L2 signal is older than L5, more satellites are capable of transmitting it. As of mid-2023, 17 out of 31 GPS satellites transmit L5 while 24 GPS satellites transmit L2 [33], all operational Galileo satellites transmit both E5b and E5a [34]. Since this report aimed to explore newer areas, the L5 signal was chosen in combination with L1. 3.4 GNSS Corrections As described in the theory, the GNSS signals are affected by multiple sources of errors. For simplicity, the corrections of these errors were handled either offline, based on calculations of only the GNSS data or online, in relation to a filter estimation. Offline corrections were made for all steps of the observed data before the actual estimation of the position began, while online corrections were performed at each step. However, neither of the corrections relied on future data. Corrections like the multipath and receiver noise were not handled explicitly in the im- plementation. The choice of location and placement of the antenna for the data collection mitigated the multipath. The risk of multipath occurring was significantly decreased by selecting a location without large disruptions, like buildings. The receiver noise is depen- dent on the type of receiver and antenna used. This noise was estimated together with the receiver clock bias. 3.4.1 Offline Corrections The offline corrections were calculated for all time steps before the estimation of the position. Some of the error sources that were corrected offline were the emission time, earth rotation, relativistic effects and solid tides. The processes described in Sections 27 3. Methods 2.4.1-2.4.3 and 2.4.12.1 were followed to correct these errors, thus will not be covered in greater detail in this section. 3.4.1.1 Cycle Slip Cycle slip detection was based on finding either a data gap or unexpected behaviour in the signal. The signal observed for the cycle slip detection was the difference between the frequencies of the carrier phase signals, i.e. Φslips = Φ1 − Φ5. The data gap was identified when the last three data samples were unavailable. To find unexpected behaviour, a second-degree polynomial was fitted to the last 40 data samples, and the prediction of the next time step, based on this curve, was compared to the true value of the next time step. A cycle slip was declared if the difference between the two were larger than one meter. 3.4.1.2 Ionosphere One of the significant atmospheric effects, the ionosphere, can be corrected using the ionosphere-free combination described in Paragraph 2.2.3.1. Equations 2.7 and 2.6 were used to calculate new code and phase measurements using frequencies L1 and L5. 3.4.1.3 Clock And Orbit Errors The clock and orbit errors were corrected using precise orbit and clock corrections, which the Crustal Dynamics Data Information System (CDDIS) provides on a regular basis on the internet [35]. The orbit corrections are published at three intervals, ultra-rapid, rapid and final [36]. The three types of files are published 4 times per day, 17 hours after the end of the previous day and 13 days after the end of the solution week, respectively. These files differ in accuracy, with the most accurate being the final type. The clock corrections are only available in rapid and final form [37]. The final form of both the orbit and the clock corrections were utilised in this implementation. The clock file contains the satellite clock error, dts, for all satellites recorded at specific intervals. In this case, the clock file had a 30-second interval. Since data collection occurred at 1Hz in this implementation, interpolation of the clock errors was necessary. A first-degree polynomial interpolation was employed to achieve centimetre-level accuracy in clock error estimations [4]. The orbit files provide the precise positions of all satellites in the ECEF coordinate sys- tem. Similar to the clock file, orbit files are generated at various intervals. Here, an orbit solution with a five-minute interval between data points was utilised. A 10-degree polyno- mial interpolation was applied to capture the dynamics of satellite movement accurately [4]. 3.4.2 Online Corrections Contrary to offline corrections, online corrections were calculated within the filter loop at each time step. Therefore, these corrections were dependent on the latest estimated position and could not be calculated beforehand. 28 3. Methods 3.4.2.1 Troposphere There are multiple ways of determining the mapping functions mentioned in Section 2.4.6.2, mh and mw, depending on what input is available. One option is the Mapping of Niell, introduced in 1996 [38]. The hydrostatic mapping function, mh, and the wet mapping function, mw, are given by Equations (3.1) and (3.2). mh(el,H) = m(el, ah, bh, ch) + ( 1 sin(el) −m(el, aht, bht, cht) ) ·H (3.1) mw(el) = m(el, aw, bw, cw) (3.2) Where • el is the elevation angle • H is the height of the receiver. aw, bw, cw were determined through linear interpolation in Table A.2. aht, bht, cht can be found in Table A.1. ah, bh, ch were calculated based on Equation (3.4) and linear interpolation of Table A.1. m(el, a, b, c) is the mapping normalized to unity at zenith [39], shown in Equation (3.3). m(el, a, b, c) = 1 + a 1+ b 1+c sin(el) + a sin(el)+ b sin(el)+c (3.3) The average and amplitude values correspond to ξavg and ξamp respectively in Equation (3.4) and were used to determine ah, bh, ch through. ξ(ϕ, t) = ξavg(ϕ) − ξamp(ϕ) · cos ( 2π t− T0 365.25 ) (3.4) Where • ϕ is the latitude of the receiver • t is the day of year (DOY) • T0 is constant with a value of 28 for the northern hemisphere [40] • ξavg and ξamp are given by linear interpolation from Table A.1. 29 3. Methods The hydrostatic and wet delay, ZHD and ZWD, were determined using Equations (3.5) and (3.6) [40]. ZHD = 2.3 · e−0.116·10−3·H (3.5) ZWD = ZWD0 + ∆ZWD (3.6) Where • ZWD0 = 0.1 • ∆ZWD is estimated as a random walk process. 3.5 Data Collection GNSS data was needed to test and validate the implementation of the corrections as well as the entire system. The tests included analysing how the different corrections affect the accuracy, how fast the solution converges, and the overall accuracy the system could achieve. As implied by the project’s limitations, the data was only collected in a stationary scenario. The data was acquired with the receiver described in Section 3.1, the EVK-F9P evaluation kit. The data collection took place on March 4, 2024, in Mölndal, Sweden, on the roof of a building. The collection began at 17.29 UTC+1 and ended at 23.59 UTC+1. The site was expected to be a low-multipath environment due to the height of the receiver placement. The receiver was set to have a 1Hz reception frequency and to output raw data, the u-blox SFBRX and RAWX messages [41]. However, only the RAWX messages were used in this project as they contained the pseudoranges and carrier phase measurements. The ionospheric effect can be monitored through a service provided through the national network of stationary reference stations SWEPOS [42], which allows for real-time and historical ionosphere monitoring. The ionospheric effect during the day of the data col- lection session is visualized in Figure 3.2. These graphs depict the ionosphere changes in the Götaland region, a region in the south of Sweden, where Mölndal is located. Figure 3.2: Illustration of the ionospheric effect during 04-03-2024 where the bottom row indicates the hour of the day in UTC+1 [42]. 30 3. Methods The colours in the figure represent how much the ionosphere affects the uncertainty of the measurements. The top two red regions indicate high uncertainty, making the mea- surements difficult to trust. In the two middle rows, the uncertainty gradually decreases. The bottom two rows, the green region, imply a low uncertainty. During the data collection period, 17.29 to 23.59 UTC+1, the ionospheric effects were relatively low, with a slight elevation around 19 to 21 UTC+1, as seen in Figure 3.2. 3.6 Extended Kalman Filter The developed solution equipped an EKF to estimate the position, a popular filter choice according to [20]. It was not possible to use a regular Kalman filter in this application since the measurement equations in (2.3) and (2.5) are non-linear with respect to the position of the receiver, which needs to be estimated as a state. The EKF, described in Section 2.5, forms the Kalman matrices by linearising the process- and measurement- equation around an operating point at each time step. In this case, only the measurement equation had to be linearised since the process function was linear, as described further down in this section. The state vector in the developed EKF, shown in Equation (3.8), contained the positional states δx, δy and δz that described the change in position between time steps. The actual estimated position in ECEF at time step k was therefore retrieved by Equation (3.7). r̃k = [ δxk δyk δzk ]T + r̃k−1 (3.7) The state vector also included the states cdtr, which is the receiver clock bias combined with the light speed factor to produce an estimate in meters. An inter system bias, bISB, between GPS and Galileo, was also included since there could be biases between the different constellations. The second to last state to be estimated was the change in ZWD, ∆ZWD. The last states, denoted ÑIF , are the float ambiguities from the ionosphere-free combination of carrier phase measurements. X = [ δx δy δz cdtr bgal ISB ∆ZWD ÑIF ]T (3.8) The measurement vector was designed as in Equation (3.9), where some of the corrections were applied. Note that the variables with a superscript of gal or gps denote all of the satellites from either Galileo or GPS. Z =  P gal IF P gps IF Φgal IF Φgps IF  =  P gal IF,raw + dtgal − ∆gal rel − T gal p P gps IF,raw + dtgps − ∆gps rel − T gps p Φgal IF,raw + dtgal − ∆gal rel − T gal p Φgps IF,raw + dtgps − ∆gps rel − T gps p  (3.9) Where • T sys p = ZHD · msys h + ZWD0 · msys w is the partial tropospheric correction. 31 3. Methods The prediction step was executed using the process matrix, denoted as A in Equation (3.10). This matrix corresponded to a kinematic model for the positional states, implying that these states are independent of their preceding values. The process noise matrix, Q in Equation (3.11), describes the uncertainties of the process model. For the positional states, an uncertainty of 1m was chosen even if the position of the receiver was known to be stationary. Since the sampling frequency was 1Hz, this resulted in a variance of 1m/s. The receiver clock error, inter system bias and tropospheric uncertainty were chosen based on literature and results. The ambiguities were designed to equal their previous value with no uncertainty since they are constant as long as no cycle slip occurs. A = diag ( 05, 1n+1 ) (3.10) Q = diag ( 13, 1e6, 1e−2, 0.2, 0n ) (3.11) The initial uncertainties needed for the Kalman filter were as written in (3.12). The initial uncertainty was redundant for the states with a zero in the process matrix A due to the formula in Equation (2.35). The initial uncertainty for the ambiguities was set to a large value such that they quickly adapted to the measurement. P0 = diag ( 06, 1n · 1010 ) (3.12) To retrieve the observation matrix, H , the measurement equations (Equations (2.5) and (2.3)) had to be differentiated with respect to the states as in Equation (3.13). H =  ∂P gal,0 IF ∂x0 . . . ∂P gal,0 IF ∂xn ... . . . ... ∂P gal,i IF ∂x0 . . . ∂P gal,i IF ∂xn ∂P gps,0 IF ∂x0 . . . ∂P gps,0 IF ∂xn ... . . . ... ∂P gps,j IF ∂x0 . . . ∂P gps,j IF ∂xn ∂Φgal,0 IF ∂x0 . . . ∂Φgal,0 IF ∂xn ... . . . ... ∂Φgal,i IF ∂x0 . . . ∂Φgal,i IF ∂xn ∂Φgps,0 IF ∂x0 . . . ∂Φgps,0 IF ∂xn ... . . . ... ∂Φgps,j IF ∂x0 . . . ∂Φgps,j IF ∂xn  =  xr−xgal,0 ρ0 yr−ygal,0 ρ0 zr−zgal,0 ρ0 1 1 mwgal,0 0 . . . 0 ... . . . xr−xgal,i ρi yr−ygal,i ρi zr−zgal,i ρi 1 1 mwgal,i 0 . . . 0 xr−xgps,0 ρ0 yr−ygps,0 ρ0 zr−zgps,0 ρ0 1 0 mwgps,0 0 . . . 0 ... . . . xr−xgps,j ρj yr−ygps,j ρj zr−zgps,j ρj 1 0 mwgps,j 0 . . . 0 xr−xgal,0 ρ0 yr−ygal,0 ρ0 zr−zgal,0 ρ0 1 1 mwgal,0 1 . . . 0 ... . . . xr−xgal,i ρi yr−ygal,i ρi zr−zgal,i ρi 1 1 mwgal,i 0 . . . 1 xr−xgps,0 ρ0 yr−ygps,0 ρ0 zr−zgps,0 ρ0 1 0 mwgps,0 1 . . . 0 ... . . . xr−xgps,j ρj yr−ygps,j ρj zr−zgps,j ρj 1 0 mwgps,j 0 . . . 1  (3.13) 32 3. Methods The measurement noise shown in Equation (3.14) depended on the elevation angle of the satellite, el, where a greater elevation angle resulted in a lower noise level on the measurement. The first entry in R accounted for the pseudorange measurements with a higher noise and the latter for all the carrier phase measurements with lesser noise. R = diag ( 0.52 + 0.52 sin(el)2 , 0.0012 + 0.0012 sin(el)2 ) (3.14) Where • el is the vector containing the elevation angle for all satellites. In the EKF loop, it was necessary to handle cycle slips, as detailed in 2.4.4. This is due to the fact that measurements affected by cycle slips can negatively impact the accuracy of the position estimate. This negative effect occurs because the ambiguity linked to the cycle-slipped measurement will be incorrect, essentially forcing the other states to adjust in response. When a cycle slip was detected, increasing the uncertainty of the ambiguity associated with the cycle-slipped measurement was required. This allowed the ambiguity to be re-estimated and assigned a new, more accurate value. In practical terms, this was achieved by augmenting the uncertainty in the P matrix within the EKF loop specifically for the ambiguity state. This adjustment helped mitigate the effects of cycle slips and maintained the precision of the position estimate. A robust weighing function was also developed based on the innovation such that outliers were rejected, and estimations with a low innovation would be weighted more. The function is described in Equation (3.15). w(rk) =  1 |rk| ≤ co |rk| co ( c1−co c1−|rk| )2 co < |rk| ≤ c1 ∞ |rk| > c1 (3.15) Where • rk is the standardised residuals of the innovation rk = vk−E[vk] Var[vk] • c0 is a constant with a value of 1 • c1 is a constant with value 100 for pseudorange measurements and 4 for carrier phase measurements. This resulted in a weight equal to one for the measurements close to the mean of the innovation, a medium weight was applied to the measurements that were in the vicinity of the mean of the innovation most likely due to noise. An infinite weight was applied to the measurements that were far away from the mean of the innovation, as they most likely were outliers. After the prediction and update steps were performed, the ambiguities had to be fixed to get a precise result. The ambiguity resolution is described in the next section. 33 3. Methods 3.7 Ambiguity Resolution Ambiguity resolution in GNSS is a crucial process in PPP to accurately and quickly de- termine the distance to a satellite using carrier-phase measurement. The process consists of separating the receiver and satellite’s float biases with an integer ambiguity term. The process of fixing the ambiguities was based on Single Differenced (SD) measurements using a reference satellite. The reference satellite was chosen based on its elevation angle. In the first time step, the satellite with the largest elevation angle, i.e., closest to the zenith, was chosen. The same satellite was used until it reached a 20-degree elevation angle or a cycle slip was detected for the specific satellite. The GNSS constellations were handled separately; hence, one Galileo and one GPS reference satellite were chosen. One method of fixing the ambiguities handles the wide-lane and narrow-lane ambiguities consecutively [43]. The first step is to fix the wide-lane ambiguities. The float narrow-lane ambiguities are calculated based on the fixed wide-lane ambiguities and fixed using the Least-Squares Ambiguity Decorrelation (LAMBDA) method [44]. A modified LAMBDA (MLAMBDA) method can also be used. MLAMBDA is a faster and less computation- ally complex method [45]. Finally, the fixed WL and NL ambiguities are combined to determine the fixed IF ambiguities. Handling the WL and NL ambiguities consecutively eliminates the receiver phase bias, and the satellite phase biases are corrected. This described method was used and explained below. 3.7.1 Wide-Lane Ambiguity Fixing The Melbourne-Wübbena (MW) combination, see Equation 2.15, was used to fix the WL ambiguities. The SD MW measurements, ˜MW SD, were calculated as suggested by Equations (3.16)-(3.18). The SD operation removed the receiver hardware bias, leaving only the satellite hardware bias. MW gal SD = MW gal −MW gal ref (3.16) MW gps SD = MWgps −MW gps ref (3.17) ˜MW SD = [ MW gal SD MW gps SD ]T (3.18) The fixed WL ambiguities were calculated by fixing ˜MW SD. Due to the long wavelength, the ˜MW SD was fixed by rounding the average of the last 5 time steps to an integer as shown in Equation (3.19). MWfixed = WLfixed = round(⟨ ˜MW SD + bs W L,SD⟩5) (3.19) Where • bs W L,SD is the satellite hardware bias for WL combination [46]. 34 3. Methods 3.7.2 Narrow-Lane Ambiguity Fixing The next step was fixing the NL ambiguities by single differencing using the same ref- erence satellites that were used when calculating the ˜MW SD. The process shown in Equations (3.16)-(3.18) was repeated for the IF measurements, resulting in the IF SD float ambiguities ÑIF,SD. The NL float ambiguities were then calculated according to Equation (3.20), which combined the IF and WL float ambiguities [46]. ÑLSD = ÑIF,SD L1 + L5 λ1 · L1 − WLfixed L5 L1 − L5 + bs NL,SD (3.20) Where • bs NL,SD is the satellite hardware bias for NL combination. The covariance of the IF SD float ambiguities, PIF,SD, was calculated according to Equa- tion (3.21). PIF,SD = D · PIF · DT (3.21) Where PIF is the covariance of the float IF ambiguities calculated by an EKF, and D is the SD design matrix. The size and layout of matrix D were dependent on the actual time steps number of satellites and the index of the reference satellites. D is a modified identity matrix where the columns of the reference satellites are -1 for the rows corresponding to the satellite in its constellation. Additionally, the rows with the indices of the reference satellites were removed. An example of how D can look is shown in Equation (3.22) where the reference satellite for Galileo has an index of 2 and for GPS has an index of 3, where the number of Galileo satellites is 3 and GPS satellites is 4. D =  1 −1 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 1 0 −1 0 0 0 0 0 1 −1 0 0 0 0 0 0 −1 1  (3.22) The covariance of the NL float ambiguities, PNL,SD, was then determined using the covariance propagation law as in Equation (3.23). PNL,SD = 1 (λNL)2PIF,SD (3.23) Where • λNL = c L1+L5 is the wavelength of the narrow-lane combination. 35 3. Methods Since the ÑLSD ambiguities still were not fixed, the MLAMBDA method had to be used, which decorrelated and fixed the float ambiguities based on the ambiguities, ÑLSD, and their covariance, PNL,SD. 3.7.3 IF Ambiguity Fixing Based on the result of the MLAMBDA function, NLfixed, and the fixed WL ambiguities, WLfixed, the final fixed ambiguities, NIF,fixed, were calculated according to Equation (3.24). NIF,fixed = λNL · NLfixed − λ5 γ5 − 1WLfixed (3.24) Where γ5 = f 2 1 /f 2 5 . After the ionosphere-free ambiguities were fixed, in the next time- step, the fixed ambiguities were used instead of the float ones until a cycle slip occurred and the fixing process needed to restart. 3.8 Implementation Overview An overview of the implemented solution can be seen in Figure 3.3. The implementation included three parsers: one for the raw GNSS data, one for the precise clock and orbit products and one for the APO corrections. The parsed data from the input sources computed the offline corrections listed to the right in Figure 3.3. Based on these offline corrections, it was possible to calculate the DOP accuracies, including horizontal and vertical DOP, HDOP and VDOP. After the offline corrections, the estimation process was done using an EKF, which estimated the state vector in Equation (3.8) for each time step. 36 3. Methods UBX Parser UBX Offline Corrections SP3, CLK Parser Precise Products DOP calculation EKF Estimated Position SP3, CLK interpolation APO correction Adjust reception time DCB correction Adjust Earth rotation Cycle slip detection Az and El calculation WL ambiguity Solid tides Relativistic effects Plotting  Tool True Position ATX Parser Troposphere calculation Position estimation Ambiguity resolution Figure 3.3: Overview of software structure. The Extended Kalman Filter (EKF) is depicted in greater detail in Figure 3.4, beginning with the prediction and update phases, then proceeding to cycle slip management and ambiguity resolution. Upon detection of a cycle slip, the uncertainty of the associated ambiguity estimate was increased to mitigate its impact on other state estimates, as out- lined in Section 3.6. Without a cycle slip, the process assessed whether the ambiguities have been fixed since a previous step, if they were fixed, continue with the same ambi- guities. If the ambiguities were not fixed, try to fix them with ambiguity resolution and run the loop again. 37 3. Methods Ambiguity Estimation Yes Cycle slip? No Yes Fixed Ambiguities? No Yes Fixed Ambiguities? No Figure 3.4: Overview of EKF. 3.9 Comparison The implemented solution was compared to existing packages that process GNSS data. RTKLIB is a GNSS processing tool that can estimate position using several different modes, combinations and corrections [47]. The tool was configured to have the same information as the proposed solution, including precise orbit and clock files. The options for the tool are shown in Figure 3.5. 38 3. Methods Figure 3.5: RTKLIB options that was used. 39 3. Methods 40 4 Results This Chapter presents the outcomes of the applied corrections, the horizontal and vertical DOP values, and the estimated position’s accuracy in horizontal and vertical form. The results from RTKLIB are also depicted. The Figures in this section are based on the data described in Section 3.5. 4.1 Corrections The resulting effects of the implemented corrections are visualised in Figures 4.1-4.8. The amount of epochs included in each plot varies depending on how long the satellite is visible to the receiver. For the case where all satellite corrections are included in a single figure, the entire 19,000 time steps of data collection are visible. The relativistic clock effect in meters is shown in Figure 4.1. The relativistic effects on each satellite that are included in the data are visible in the figure. The relativistic correction, ∆rel, is applied to the pseudorange and carrier phase measurements, as shown in Equations (2.3) and (2.5). The effect of the relativistic clock varies in range between satellites. The smallest values are around 0m, while the largest effect for this data is around ±60m. Most values are within the ±7.5m range, depicted in Figure 4.2. The Figure presents a close-up of the same values as in Figure 4.1 and indicates that the correction behaves like a sinusoidal curve for each satellite. Figure 4.1: Relativistic clock effect in meters for all Galileo and GPS satellites in each time step. 41 4. Results Figure 4.2: Close-up of relativistic clock effect in meters for all satellites in each time step. Figure 4.3 depicts the effect of the solid tides. The Figure contains the correction in the X-, Y- and Z-directions (which correspond to the east, north and height directions) applied to the receiver position. The effect is relatively small and ranges between 0.2m and 0.7m. Figure 4.3: Solid Tides effect in meters for each time step. The earth rotation correction for the position of one satellite is shown in Figure 4.4. Since the correction is consistent at around ±200m for all Galileo and GPS satellites, the behaviour of a single satellite is representative. The effect in the Z-direction is 0 for each satellite, as the earth’s rotation only affects the position of the satellite relative to the receiver in the X- and Y-directions due to its rotation about the Z-axis. Figure 4.4: Earth Rotation effect in meters for 15,000 time steps for one satellite. 42 4. Results The satellite clock bias, dts, has an immense effect on the accuracy of the system, as visualised in Figure 4.5. The correction of the pseudorange and carrier phase measure- ments for all satellites in every time step is included. Most corrections are within the ±200,000m range and do not vary during the period of the collected data. Figure 4.5: Clock Correction in meters for each time step. Figure 4.6 depicts the correction of the tropospheric effect, T s = ZHD · ms h + (ZWD0 + ∆ZWD) · ms w, for one satellite. It is rather stable for the first 15,000 time steps, with a correction between 2.5m to 5m. During the following 2,500 epochs, the correction increases slightly. The correction then decreases and finally reaches a value of about -5m. The correction is applied to the pseudorange and carrier phase measurements. Though the correction will differ for each satellite, they are generally not larger than 20m. Figure 4.6: Troposphere correction for 19,000 time steps for one satellite. The behaviour of the cycle slip detection for a satellite is shown in Figure 4.7. A vertical black line represents the time steps where a cycle slip is detected. The Figure also includes the true measurement Φ1 − Φ5 (blue) and predicted (red) measurements. The true and predicted measurements overlap for the first 12 time steps, and no cycle clips are detected. However, at the 13th epoch, a cycle slip occurs in the true measurement, causing an offset between the prediction and true measurement and a cycle slip is declared. Cycle slips are present until the true and predicted measurements overlap once again and overlap for at least 10 epochs. 43 4. Results Figure 4.7: Detected cycle slips for 70 time steps for one satellite. Figure 4.8 displays one satellite’s antenna phase offset. This correction is applied to the satellite position before the position is interpolated to fit the data time steps. The precise product has 289 time steps with an available position to which the correction is applied. The correction has a ±1.5m magnitude in the x, y and z directions. The size, as well as the behaviour of the correction is similar for all satellites. Figure 4.8: Antenna Phase Offset Correction in meters for the time steps available in the precise products for one satellite. 4.2 Dilution Of Precision The expected accuracy of the observed data is shown as DOP graphs in Figures 4.9 and 4.10. Figure 4.9 depicts the horizontal DOP values for the first 19,000 epochs, and Figure 4.10 shows the corresponding vertical DOP values for the same data. An elevation angle cut-off is applied at 10◦ before calculating the DOP values to only include satellites that are used in the position estimation. 44 4. Results Figure 4.9: Horizontal DOP value for each time step using satellites with elevation angle larger than 10◦. As indicated by Figure 4.9, the horizontal DOP value is mostly between 0.7 and 1.2. Even the exceptions, the time steps where the HDOP exceeds a value of 1.3, are still close to 1. The outliers that do not follow the general graphs indicate a loss of connection to an important satellite. The vertical DOP values presented in Figure 4.10 are similar in magnitude to those of the HDOP. All VDOP values are close to 1, within a range of 0.8 to 1.35. The values indicate low uncertainties for the data. Figure 4.10: Vertical DOP value for each time step using satellites with elevation angle larger than 10◦. 4.3 Number of Satellites The number of visible satellites during the data capture is depicted in Figure 4.11, where the difference is shown between the amount of L1 capable satellites and L5 capable satellites. The number of L5 transmitting satellites is lower during the whole capture. 45 4. Results Figure 4.11: Number of satellites for the different frequency bands L1 and L5. The number of satellites able to broadcast both the L1 and L5 frequencies with an elevation angle above 30° is depicted in Figure 4.12. It can be seen in the figure that the number of satellites with a higher elevation angle (> 30°) is at its lowest at around time step 11, 000 and highest at the start of the data capture. Figure 4.12: Number of satellites broadcasting L1 and L5 above 30° elevation angle. Figure 4.13 displays the count of usable satellites, defined as those exceeding the elevation angle cut-off of 10°, free from cycle-slips and transmitting both L1 and L5 frequencies. Figure 4.13: Number of usable satellites for the ionosphere-free combination of L1 and L5. 46 4. Results 4.4 Vertical Accuracy The vertical error refers to the accuracy of the height plane and is depicted in Figure 4.14. The Figure includes the true and estimated position. The graph showing the true position includes blue and red sections, which are indications of the state of the ambiguity. The graph is blue for epochs with float ambiguities and red when the ambiguities are fixed. There is a large initial offset, off about 15m, before the estimation closes in on a 2m offset in the early epochs. The vertical accuracy is rather stable during epochs 4,000 to 11,000, averaging an offset of about 2.5m. After epoch 11,000, the accuracy decreases and is no longer stable. Figure 4.14: The offset in the height plane from the true position for 19,000 time steps. 4.5 Horizontal Accuracy The horizontal accuracy, i.e. the accuracy in the East/West and North/South directions, is determined by analysing the distance between the estimated and the true position. The offset in the East/West direction is shown in Figure 4.15, and the North/South offset is depicted in Figure 4.16. As with the vertical accuracy, the colours of the graph for the estimated position indicate the state of the ambiguity for each epoch. The ambiguities are fixed during four periods, all happening before time step 11,000. The offset in the East direction is stable and within decimetre-level during the first 11,000 epochs, with the exception of the large offset of the initial position. The accuracy is high until the second part of the data, the last 8,000 epochs. The estimated position starts to drift, and the accuracy decreases continuously. 47 4. Results Figure 4.15: The offset in the East direction from the true position for 19,000 time steps. Similar behaviours are visualised in the North direction in Figure 4.16, where the esti- mated position is stable until epoch 11,000. The initial uncertainty is present in this estimated position, but there is a fast recovery before converging to an offset of about 3m. As seen previously, the graph’s behaviour changes at step 11,000. Figure 4.16: The offset in the North direction from the true position for 19,000 time steps. The first 11,000 time steps for the East and North offsets are plotted in Figures 4.17 and 4.18 due to larger uncertainties after 11,000 time steps. The position is fairly stable in both these figures but not completely static. 48 4. Results Figure 4.17: The offset in the East direction from the true position for 11,000 time steps. Figure 4.18: The offset in the North direction from the true position for 11,000 time steps. 4.6 RTKLIB Results This section will show the results of using RTKLIB to process the collected data. In Figure 4.19, the East, North and altitude positions differ by tens of meters from the known position. 49 4. Results Figure 4.19: RTKLIB result in East, North and up directions, indicated by E-W, N-S and U-D respectively. The time in UTC is on the x-axis and the offset from the true position in each direction is shown on the y-axis. In Figure 4.20, the number of usable satellites is shown. As can be seen, many of the time steps consist of ∼ 5 usable satellites and some time steps, no usable satellites are detected. Figure 4.20: RTKLIB number of satellites. 50 5 Discussion This section covers the attained accuracy of the implementation and the factors that influence it. Additionally, it includes explanations of how the results can be enhanced. The effects of the corrections, as well as their limitations, are discussed. The potential uncertainties of the collected data are also presented. Finally, the proposed future work is outlined. 5.1 Accuracy The results presented in the previous chapter match the expectations of how well this solution would work at this stage. The most important result is the positional accuracy in the horizontal plane (North and East), which delivers decimetre-level accuracy when the ambiguities are fixed. As noted, it is possible to see in Figures 4.16 and 4.15 that the positional accuracy severely degrades after ∼ 11, 000 time steps. The reason for this is not completely known, but most likely, it is related to the amount of visible and not cycle-slipped satellites. It can be seen in Figure 4.13 that the amount of usable satellites decreases during the data collection, and at its lowest, there are only five usable satellites, which is not enough to estimate a position correctly. The implemented EKF uses six states plus the ambiguity states. When the number of satellites is below six, and the ambiguities are not fixed, there are too many unknown variables compared to measurements, making the system underdetermined. The displayed DOP values in Figures 4.9 and 4.10 suggest good outcomes for both vertical and horizontal positions throughout the entire sampling duration. Given the application’s focus on horizontal accuracy, the HDOP takes importance. HDOP primarily relies on satellite azimuth angles rather than elevation angles, although the latter significantly impacts signal noise. Figure 4.12 shows a significant decrease in the number of satellites above 30° during data capture, particularly at time step 11,000. This reduction in high- elevation satellite availability may contribute to the degradation of the position estimate. More constellations could be considered in the implementation to fix the issue of having too few usable satellites. Currently, only GPS and Galileo are considered for simplicity, but GLONASS and Beidou have operational satellites in the area. GLONASS however, does not transmit the L5 signal, making all of its satellites unusable for this work. Beidou transmits L5 on some of its satellites, but the amount in this area is low due to its main operation in China. 51 5. Discussion The most promising solution would be switching hardware that is able to receive and process the L1 and L2 signal combination. There are plenty of visible satellites, but not all of them transmit the L5 signal, making them unusable for this work. As can be seen in Figure 4.11, at the lowest, there are 15 visible satellites with L1 and only 10 for the L5 signal. Since all L1 frequency transmitting satellites can transmit the L2 frequency, the number of satellites would be much higher with that combination. Before 11,000 seconds, the results look promising with decimetre precision in the horizon- tal plane as seen in Figures 4.17 a