Project logo

Introduction

GNSS Compare now supports the Dual-Frequency Xiaomi Mi 8 smartphone!

Welcome to GNSS Compare’s documentation!

Please stay tuned, this is still a beta version, and we’re making updates on a daily basis…

Check out the app on Play Store.

Get it on Play Store

Please note that you need Android 7.0+ to run the application. Also, please note that not all Android 7.0+ phones support the Galileo Satellite System or crucial for this project raw GNSS measurements. The list of compatible phones can be found here.

So far we’ve been testing the application on:

  • Samsung Galaxy S8
  • Samsung Galaxy S8+
  • Samsung Galaxy Note 8
  • Xiaomi Mi 8

And we are very happy with the results!

If you want to be kept up to date with our updates, you can sign up to our newsletter!

Description

The purpose of GNSS Compare is to make the life of developers and reasearchers easier. It’s an easy to use and easy to extend Android-based framework for calculating the PVT based on the raw GNSS measurements.

Project history

GNSS Compare is the winning Android application developed by the Galfins as part of the internal European Space Agency (ESA) competition called Galileo Smartphone App Challenge introduced in November 2017. The Challenge was about creating a smartphone application, which will allow the user to choose which satellite constellation to use for Position, Velocity and Time estimation. The aim was to increase the awareness about ESA’s Galileo programme and also to allow users from common public to compare the performance of Galileo signals with the performance from other global satellite navigation constellations.

Download

GNSS Compare is always available on the Google Play Store:

Get it on Play Store

Please note that you need Android 7.0+ to run the application. Also, please note that not all Android 7.0+ phones support the Galileo Satellite System. The list of Galileo compatible phones can be found here.

So far we’ve been testing the application on:

  • Samsung Galaxy S8
  • Samsung Galaxy S8+
  • Samsung Galaxy Note 8
  • Xiaomi Mi 8

And we are very happy with the results!

Hall of Fame

In this section we present you the amazing people that were fueled by enormous amounts of tea [*] in order to develop GNSS Compare!

[*]Typically it’s coffee, but we like to approach things differently. 😊

And here they are! From the left to right we have: Mareike Burba, Sebastian Ciuban, Dominika Perz and Mateusz Krainski.

TheGalfins

They all have an impressive set of skills that brought the 1st prize for GNSS Compare at the European Space Agency Galileo Smartphone App Challenge. Their knowledge and experience related to Software Design and GNSS Signal Processing are presented in short biographies so you can have an idea who is behind GNSS Compare!

Software Design

Mateusz Krainski

Mat

Polish Young Graduate Trainee in the Directorate of Human Spaceflight and Robotic Exploration, at ESTEC, where he supports the European Robotic Arm (ERA) team. His main duties regard designing, developing and validating a robotic testbed for testing of ERA’s on-board smart cameras. During his studies, Mateusz was one of the key board members of a robotic student society, where he managed numerous projects ranging from small teams for quick projects (this includes a Space Startup Weekend, an Android app hackathon and few duringstudies assignments), organizing robotic tournaments (with a team of over 15 people), up to technical projects counting over 30 people. Thanks to the Toastmasters International community, Mateusz has developed highly his public speaking skills. He not only helped start the first English speaking club in the area, but also received awards in presenting competitions on a semi-national level.

Dominika Perz

Dominika

Polish Young Graduate Trainee in the Technology, Engineering and Quality Directorate at ESTEC, ESA. Currently working as a Project Manager for the Lunar Exploration Mission - internal project to determine a preliminary GNC (Guidance, Navigation and Control) design for the ascent, rendezvous and docking with the Deep Space Gateway station orbiting a Moon in the highly elliptical orbit. Her background is mainly in robotics and programming. She completed Control Engineering and Robotics master studies in Poland, during which she spent one semester in Madrid, Spain as an Erasmus exchange student. As a first international carrier experience, Dominika did a 6 weeks internship in R&D team at Venderlande (Eindhoven, Netherlands), where she worked on optimisation of the search algorithm. During holidays in 2016 she participated in the Aerospace Information Technology Summer School in Würzburg, Germany. Before coming to ESA, Dominika worked for a year at a software company GlobalLogic as a Junior Software Developer for embedded systems.

GNSS Signal Processing

Sebastian Ciuban

Sebastian

Romanian Young Graduate Trainee in the Directorate of Navigation at the European Space Research and Technology Centre (ESTEC) working on the Galileo project. He holds a Master of Science degree in Aerospace Systems – Navigation and Telecommunications granted by the French Civil Aviation University (ENAC) from Toulouse, France. After his first year of master studies he worked a summer (2016) at Acorde Technologies in Spain in the area of GPS/INS integration. During the second year of studies he developed his master thesis at the German Aerospace Center (DLR) in Oberpfaffenhofen. While being at DLR he was responsible with designing and implementing algorithms that fused the Precise Point Positioning (PPP) technique with Inertial Navigation Systems (INS) in a tight coupling architecture. Moreover, he has also developed suitable integrity monitoring algorithms in order to measure the reliability of the designed fused systems in terms of fault detection sensitivity and protection level stability. His research interests are related to precise positioning, sensor fusion, integrity monitoring and GNSS receiver signal processing.

Mareike Burba

Mareike

German National Trainee with EOP-SMS, currently working on the atmospheric correction for the Fluorescence Explorer (FLEX). She holds a M.Sc. in Meteorology with a minor in Scientific Computing from the University of Hamburg. Her Master’s thesis was about optimizing the numerical computation of atmospheric radiative transfer for the Infrared Atmospheric Sounding Interferometer (IASI) on MetOp. Mareike previously joined the World Climate Research Program Joint Planning Staff in the World Meteorological Organization for a 6 months internship in the framework of the Carlo Schmid program. Thanks to her studies and jobs as student research assistent, she speaks fluently Python, Matlab, Bash and Fortran, is communicative in C, Java and IDL, and has some experience with parallel computing.

Contact

If you would like to reach any of the GNSS Compare developers, here are their e-mails:

Glossary

The aim of this section is to provide a small glossary on the subject of satellite navigation, Android and software engineering in general. You can find below descriptions of all used technical terms in this documentation.

Android Glossary

Software Engineering Glossary

Polymorphism

According to Wikipedia, Polymorphism is the provision of a single interface to entities of different types. In Java this is achieved due to class inheritance.

Acknowledgements

We would like to mention the libraries that we are using in the GNSS Compare software framework:

Also we would like to thank the European Space Agency for the continuous support throughout the development of this application. Especially we would like to mention the ESA employees that shared their experience and knowledge with us: Paolo Crosta, Nityaporn Sirikan, Gaetano Galluzzo and Paolo Zoccarato.

Newsletter

Dual-Frequency capabailities for GNSS Compare are currently under development and testing. Stay tuned!

GNSS Compare is currently under heavy development. If you want to be up to date with our changes, sign up to the newsletter here.

Introduction

Let’s start with the basics. GNSS compare allows the user to calculate the phone’s location, based on the phone’s GNSS pseudorange measurements. This is normally done by the phone automatically, however with the recent release of Android API 24 and the GnssMeasurement class, developers gained access to unprocessed pseudorange measurements. This can be used by GNSS scientists and researchers to come up with new, more precise or less resource intense methods for precise positioning. GNSS Compare is basically a tool for such scientists to compare their algorithms. And if you’re not an GNSS expert, it can be a tool for you to learn more on this subject. And believe me – there’s a lot of interesting things to learn.

In order to fully understand what’s under the hood for GNSS Compare, we must introduce a few terms we are using through the application. It’s going to be a very general explanation – if something will be not clear, please take a look at our Glossary. We hope you’ll find answers to more detailed questions there. This section aims to provide a bridge between general GNSS knowledge, software practices, and GNSS Compare.

Getting started with the User Interface

Here we show what you can see at the User Interface (UI) level and we also describe how to set up different processing schemes/calculation modules, effect of which can be studied in real-time.

Application’s Views

We refer to a view as the current content displayed by the application. The user can change these views by swiping on the phone’s screen. Curently we have the following views: Main View, Satellite signal strength, Positioning error plot and Google Maps view.

Main View

When you launch the application this is the first view.

_images/MainView.gif

On top there is a blue “stripe” with the name of the application, a “+” and a “gearbox” icon. What are the functionalities of those icons we will see in the Processing schemes part.

Next is the Constellation status header and the information below that shows what GNSS constellations and how many satellites are used to compute PVT. In the GIF above we can see that a combination of Galileo+GPS, GPS only and Galileo only are considered in the algorithms. Moreover, you can notice that not all Visible satellites are being Used in the calculations. The reason behind this is explained in a dedicated chapter of this documentation called Android GNSS raw measurements. Shortly, it is because not all obtained pseudoranges pass a criteria that would allow them to be used in the PVT estimation.

Below the header called Calculation results are the results of the PVT estimators (EKF for this particular example) in terms of: latitude ( Lat ), longitude ( Lon ), altitude ( Alt ) and the receiver’s clock bias ( C.bias ). The UI allows it’s user to make some interesting analysis and to gain some intuitions about the importance of the number of the used satellites in PVT. As example, because there are only 3 Galileo satellites used in the EKF we do expect the estimations of the unknowns to be degraded, which is the case.

And lastly, the START RAW LOG allows the logging of the Android GNSS raw measurements in the exact same format as the Google’s Application GNSS Logger. This feature allows you also to analyze your data in post-processing!

To get to the next view swipe from right to left.

Satellite signal strength

This view is quite straight forward. Here you can monitor the signal strenght of the satellites that are Used in the calculations.

_images/SatelliteSignalStrenght.gif

To get to the next view swipe from right to left or to return to the previous one, from left to right.

Positioning error plot

To have an idea of how well the position is estimated, we provide this view that contains a plot with the horizontal position errors using as reference the Android FINE location (i.e., the best location output by the phone). The errors are expressed in meters in the north and east direction (local frame).

_images/PosErrorPhone.gif

Below the plot there is the legend with the specific colors for the chosen processing schemes.

To get to the next view swipe from right to left or to return to the previous one, from left to right.

Google Maps view

In the last view there is the Google Maps on which the position estimations are displayed to be monitored. This can be useful especially when you are testing new PVT algorithms or change the settings of the existing ones (e.g., tunning the EKF). In the GIF below are presented the position estimations by the EKF while the user was in a bus. In this way you can study if your algorithms and their tunning are able to output estimations that follow your dynamics in real-time.

_images/BusEkf.gif

Another useful study that can be made in this view is the comparison of different PVT algorithms. In the example below, one can gain insights about the difference between WLS and EKF. It is interesting to see the performance of an estimator that relies only on measurements relative to an estimator that uses a dynamic model in addition.

_images/WlsVsEkf.gif

To get to the the previous view swipe from left to right.

Processing schemes

By a processing scheme or a calculation module we refer to a set of settings that are considered for the estimation of the smartphone’s position. The user can create a processing scheme in which he/she can choose the following:

  • Constellation: Galileo+GPS, GPS, Galileo
  • Correction modules: Tropospheric correction, Klobuchar Iono Correction (only for GPS), Relativistic path range correction
  • Positioning method: Weighted Least Squares, Static EKF, Dynamic EKF, Pedestrian EKF
  • Logging format: Simple format, NMEA
  • Name: The user shall specify the name of the processing scheme

We will continue to show how to create a processing scheme or modify an existing one at the UI level.

Creating a new one

From the Main View select the “+” icon on the top right corner as shown below. A new view will prompt that will allow you to select among the available options.

_images/MainViewCreate2.jpg

Let’s begin by selecting the desired Constellation and Correction modules, one at a time.

_images/Const_CorrModule.jpg

Next we would like to select the Positioning method and the Logging format.

_images/Pos_Log.jpg

Finally we have to give a Name to the processing scheme and then press Create.

_images/Name_Create.jpg

Modifying an existing one

To modify an existing processing scheme, from the Main View press the “gearbox” icon as shown below.

_images/Gear_Existing.jpg

In the right image you can find the created processing schemes. Additionally, there is a summary for each scheme with the chosen options. To change the settings of an existing scheme, press Modify and choose among the existing options. Moreover, in this view you also can enable/disable a processing scheme by pressing the Activate switch or enable/disable data logging with the Save log switch. A dedicated section will be made regarding the logging formats of GNSS Compare.

Getting started with the code

Importing the project to Android Studio

It’s not necessary to use Android Studio to make modifications and build code, there are many other IDE’s for Android development. It has been our choice so far, so we will be focusing our tutorials on that IDE at the moment.

Before importing the project, you have to clone the project repository. The directory GNSS_Compare within that repository is the Android Studio project.

To import the project into Android Studio, find the Import Project option. It should be either under File -> New -> Import Project, or somewhere on the main screen of Android Studio. Navigate to the GNSS_Compare directory. You should find the GNSS_Compare project folder with the Android Studio logo next to it, as shown on image below.

_images/path_selection.PNG

Mark it, click Ok. The project should open and gradle will start to synchronize everything and genrate its configuration files. This can take a moment… After gradle has finished, you’re free to connect your phone and build the application by pressing the Run app button in the top, right corner.

_images/run_app.PNG

Using the Google Maps Viewer

In order to use the Google Maps Viewer in GNSS Comapre, you’ll need to get your own Google Maps API key and paste it into the Android manifest (kind of).

To get the Google Maps SDK key, follow this guide.

After you have the key (it will look like a string of random characters, starting with AIza), all you need to do is to copy and paste it to the map_api_key.xml file located in the res/values directory. Find the lines containing the following:

<string name="map_api_key">YOUR_API_KEY</string>

And replace YOUR_API_KEY with your API key. It should work right away. If not – you might need to clean and rebuild your project, or manually uninstall the application on your phone.

Remember not to share the api key with anyone! If you’re using git, you can mark that file as one which should not be tracked with

git update-index --assume-unchanged GNSS_Compare/app/src/main/res/values/map_api_key.xml

This way, the file will remain in your repository, but any changes made to it will be not pushed to the remote.

GNSS basics

Constellation

There are a few Global Navigation Satellite Systems on the planet. Or actually - around the planet: the Global Positioning System (GPS), owned by the US Government, GLONASS, owned by the Russian Federation, Chineese system BeiDou, and of course Galileo, owned by the European Union, which is set to soon become fully operational.

Each of those constellations consists of a set of 20-30 satellites orbiting above our heads at an altitude of around 20000 km. So many satellites are needed, because the system operator assures that from each spot on Earth, at any time of the day, at least four satellites are visible (and usually a lot more). Those satellites are constantly broadcasting a signal, which among other parameters contains a timestamp. The receiver then receives the signal, compares the timestamp with the current time, and thanks to that is able to calculate the distance to the satellite. Knowing the time of transmission and the satellites orbital parameters (retrieved from a special server in the form of ephemeris data or extracted from the received signal iteself) it’s also possible to calcualte the satellite’s position in space. Of course, the details of how those calculations are performed vary slightly from constellation to constellation.

In the context of GNSS Compare, a Constellation, is a class, which defines those two properties:

  • is capable of converting the raw measurements extracted from the phone into pseudoranges,
  • is able to calculate satellite’s position in the time of transmission.

In this context, the constellations can be treated individually, e.g. we separate GPS from Galileo, and perform position determinating calculations for them separately, or they can be treated together, as in our Galileo+GPS example. Developers must take care as combining constellations is not always that easy!

Corrections

As the signal travels from the satellite, it’s prone to a number of sources of error (e.g., ionosphere, troposphere), which the user will have to take into account. For a more accurate positioning, we need to estimate the distance to the satellite with as good as possible, so we need to remove from the signal any disturbances we might know of.

Those disturbances are estimated using various, more or less complicated, mathematical models of the natural phenomena. Those models allow us to calculate corrections, which are later applied to the pseudoranges. The most commonly used corrections are for the ionosphere, troposphere and for those including the relativistic effects.

In the context of GNSS Compare, a Correction is a class, which provides a method to calculate the value of the correction, based on few parameters, which include:

  • time of signal reception,
  • receiver’s approximate position,
  • the satellite’s position,
  • additional data, stored in the ephemeris data

The general rule is simple – the more corrections are applied, the more accurate the final position.

PVT Estimator

The PVT estimators are algorithms which take as input satellite positions and pseudoranges to those satellites and aim to estimate the receiver’s position, velocity and time. In some applications, it’s sufficient to estimate just the position and time. Let’s take a look at the parameters we wish to estimate. Position is quite obvious - that’s what we would want to get from this whole process. In some cases, we can use the signal characteristics to improve the estimation of the receiver’s velocity (e.g., using doppler measurements), thus increasing the accuracy of the position estimations. Additionally to the position related parameters we also need to estimate the receiver clock bias with respect to a certain GNSS time frame (e.g., Galileo System Time). This is handled by having the clock bias as one of the paramters to be estimated alongside with the position and velocity.

In the context of GNSS Compare, the PvtMethod class does exactly that. It’s supposed to calculate the receiver’s position, based on observed satellite parameters. Internally, it should be storing the calculated velocity and clock bias for enhanced processing, but from the point of view of GNSS Compare’s framework, at the moment, the only value used outside of the PvtMethod class is the calculated position. But hey – there’s of course room to improve.

Android GNSS raw measurements

(To be updated to support GNSS Dual-Frequency Android phones)

In the following sections we describe the algorithms used to compute the pseudoranges taken into account the used satellite navigation system. The following algorithms are based on the European GNSS Agency’s (GSA) White Paper on using GNSS Raw Measurements on Android devices.

At the code level, you can find the algorithms in the following Java classes:

  • GalileoConstellation
  • GpsConstellation

The variable names used in the description of the algorithms are the same as the ones in the GNSS Compare’s code. Moreover, the definition of each used variable (e.g., ReceivedSvTimeNanos) can be found on the Android Developer webpage or in the white paper mentioned above. We will keep things straight forward in this section.

As an additional informative note, the pseduroanges computed here are based on the ranging codes that modulate the L1 carrier signal.

Galileo

Roughly speaking, the pseudorange is the difference between the time of signal reception and the time of signal transmission multiplied by the speed of light (for a more detailed definition check the Glossary). Therefore, let’s see how we compute the time of signal reception with the Android GNSS raw parameters:

galileoTime = TimeNanos - (FullBiasNanos + BiasNanos);
tRxGalileoTOW = galileoTime % Constants.NUMBER_NANO_SECONDS_PER_WEEK;
tRxGalileoE1_2nd = galileoTime % Constants.NumberNanoSeconds100Milli;

It may look a bit strange that we compute two times of reception ( tRxGalileoTOW and tRxGalileoE1_2nd) however there is a reason behind this. We have to be aware of the fact that Galileo signals have more complex modulation schemes if compared with the legacy signals of GPS. In this sense, processing Galileo signals requires more effort from the GNSS receiver. Now in order to use the Galileo pseudoranges in the PVT estimation, these pseudoranges have to pass some kind of health check. One of these checks looks if the Time Of Week (TOW) parameter is decoded or determined from other sources (e.g., mobile network), and the other one checks if the smartphone’s GNSS receiver is locked on the Galileo E1 secondary code. We will see soon how this is handled. However, we will not deal with the theoretical background in order to reason the approach presented here because it does require some advanced receiver signal processing knowledge and at this point this is outside of our aims. In exchange, we can advise the curious minds to check a book on GNSS signal structures, like Engineering Satellite-Based Navigation and Timing: Global Navigation Satellite Systems, Signals and Receivers by John W. Betz.

Therefore we will use tRxGalileoTOW and tRxGalileoE1_2nd to compute two pseudoranges and we will use only one of them, the one that manages to pass the health check of course! Now let’s compute the time of signal transmission:

tTxGalileo = ReceivedSvTimeNanos + TimeOffsetNanos;

The two pseudoranges are:

pseudorangeTOW = (tRxGalileoTOW - tTxGalileo) * 1e-9 * Constants.SPEED_OF_LIGHT;
pseudorangeE1_2nd = ((galileoTime - tTxGalileo) % Constants.NumberNanoSeconds100Milli) * 1e-9 * Constants.SPEED_OF_LIGHT;

We have said that we need to test these two pseudoranges for some criteria. And the Java variable containing the health status or the states that we wish to find if they are true or not is:

int measState = measurement.getState();

With the help of the bitwise AND operation we can identify if the seeked states are true or not. Please check the Android Developer website to have a better insight of this process:

boolean towKnown = (measState & GnssMeasurement.STATE_TOW_KNOWN) > 0;
boolean towDecoded = (measState & GnssMeasurement.STATE_TOW_DECODED) > 0;
boolean codeLock = (measState & GnssMeasurement.STATE_GAL_E1C_2ND_CODE_LOCK) > 0;

Finally, we do the following check and we decide which computed pseudorange we use:

if ((towKnown || towDecoded)) {

    // use pseudorangeTOW

}else if (codeLock){

   // use pseudorangeE1_2nd

}

GPS

We follow a similar approach for GPS also by starting to compute the time of signal reception:

gpsTime = TimeNanos - (FullBiasNanos + BiasNanos);
tRxGPS  = gpsTime + TimeOffsetNanos;

In the next step we compute in a more straight forward way the GPS pseudorange:

weekNumberNanos = Math.floor((-1. * FullBiasNanos) / Constants.NUMBER_NANO_SECONDS_PER_WEEK)*onstants.NUMBER_NANO_SECONDS_PER_WEEK;
pseudorange = (tRxGPS - weekNumberNanos - ReceivedSvTimeNanos) / 1.0E9 * Constants.SPEED_OF_LIGHT;

We have to check if the computed pseudorange is usable in PVT or not. Therefore, we get the states status:

int measState = measurement.getState();

We apply again the bitwise AND operator to see if the TOW is decoded and if the receiver is locked on the signal’s code:

boolean codeLock = (measState & GnssMeasurement.STATE_CODE_LOCK) > 0;
boolean towDecoded = (measState & GnssMeasurement.STATE_TOW_DECODED) > 0;

Additionaly we can add an extra criteria, a criteria that checks for the uncertainty in the determined TOW:

private static final int MAXTOWUNCNS = 50;    // [nanoseconds]
boolean towUncertainty = measurement.getReceivedSvTimeUncertaintyNanos() < MAXTOWUNCNS;

Finally we decide to use the GPS pseduorange if the following check is true:

if(codeLock && towDecoded && towUncertainty && pseudorange < 1e9){

   // use pseudorange

}

Implemented PVT Algorithms

In this section we provide the theoretical aspects behind the GNSS Compare’s PVT algorithms. The information here can be associated with the following Java classes:

  • StaticExtendedKalmanFilter
  • DynamicExtendedKalmanFilter
  • PedestrianStaticExtendedKalmanFilter - this one sounds a bit strange, however bear with us as explanations will be given when the filter tunning is explained
  • WeightedLeastSquares

Extended Kalman Filter

One of the estimation techniques implemented in the GNSS Compare framework is the Kalman Filter. Taking into account that the measurement model is linearized about the time predicted position, in fact the implementation is an Extended Kalman Filter (EKF).

In this section we describe the theoretical aspects of the EKF implementation such that the curious minds can understand easily what is behind GNSS Compare’s awesome algorithms. We are interested to implement the EKF for two types of users: a static user and a dynamic user.

Therefore we will describe how the state vector is defined, or in other words, the vector containing the parameters that we wish to estimate (hint: the parameters are related to the GNSS Compare’s user position!), and also what dynamic and measurements models we have considered. And as bonus we will also write about the tunning of the EKFs.

First things first! Let’s remember the Kalman Filter equations, the implemented ones, in order to make the rest of this section more enjoyable.

We have the time prediction of the state vector (x) and it’s variance-covariance matrix (P):

\[\hat{\mathbf{x}}^-_k = \mathbf{F}_k \hat{\mathbf{x}}^+_{k-1}\]
\[\mathbf{P}^-_k = \mathbf{F}_k \mathbf{P}^+_{k-1} \mathbf{F}^{\text{T}}_k + \mathbf{Q}_k .\]

In the next step we can compute the innovation vector (gamma) and it’s variance-covariance matrix (S) with the help of the obsevation vector (z), the observation matrix (H) and the measurement noise matrix (R):

\[\boldsymbol{\gamma}_k = \mathbf{z}_k - \mathbf{H}_k\hat{\mathbf{x}}^-_k\]
\[\mathbf{S}_k = \mathbf{H}_k \mathbf{P}^-_k \mathbf{H}_k^{\text{T}} + \mathbf{R}_k.\]

We are almost there, we just need to compute the famous Kalman gain (K)!

\[\mathbf{K}_k = \mathbf{P}^-_k \mathbf{H}_k^{\text{T}} \mathbf{S}^{-1}_k.\]

Finally the measurement update step is:

\[\hat{\mathbf{x}}^+_k = \hat{\mathbf{x}}^-_k + \mathbf{K}_k \boldsymbol{\gamma}_k\]
\[\mathbf{P}^+_k = \left(\mathbf{I}_k - \mathbf{K}_k \mathbf{H}_k \right) \mathbf{P}^-_k.\]

However, before explaining how the EKF for the static user and the dynamic user was implemented, we still need to talk about the measurement model based on the GNSS pseudoranges retrieved from the smartphone’s GNSS receiver. If you are familiar with this concepts, you can skip the following section.

Pseudorange measurement model

For a code-based pseudorange (PRc) we have the following (non-linear) equation taking into account the satellite clock bias (dtS), the delay caused by the ionosphere (dion), the delay caused by the troposphere (dtrop) and the receiver noise (epsilon).

\[PR_c = \rho + \delta t_R - \delta t^S + d_{\text{ion}} + d_{\text{trop}} + \mathbf{\epsilon}\]

We know, there are more effects that are perturbing the GNSS measurements, however we wish to keep things as simple as possible and the interested persons can always consider more error sources in the GNSS Compare’s code.

The above equation is non-linear because of the geometric distance (rho) between the receiver and the GNSS satellite. Luckly we can linearize it if we have knowledge about an approximated position of the receiver (X0, Y0, Z0), which we do! We do have from the time prediction step of the EKF. Taking this into account and applying a first order Taylor series expansion we obtain:

\[PR_c - \rho_0 + \delta t^S - d_{0,\text{ion}} - d_{0,\text{trop}} = -\frac{X^S-X_0}{\rho_0}\Delta X-\frac{Y^S-Y_0}{\rho_0}\Delta Y-\frac{Z^S-Z_0}{\rho_0}\Delta Z+\delta t_R.\]

On the left side of the equation we have moved every term that can be computed. The subscript 0 means that those parameters are estimated by using the approximate receiver position information. On the right hand side we have the unknowns (dX, dY, dZ, dtR) and their coefficients. Based on the linearized pseudorange equation one can form the observation matrix (H).

Practical advise: Take care that the unknowns from the linearized pseudorange equations are not the same as the position related unknowns that we are estimating directly in the EKF state vector. Check the GNSS Compare code ( e.g., StaticExtendedKalmanFilter class ) to understand how this is handled.

Good, now we can see how the EKF was implemented for the static user and the dynamic user!

Static user

In the case of a static user we have the following state vector at the epoch k:

\[\mathbf{x}_k = \left(X~~Y~~Z~~\delta t_R~~\dot{\delta t}_R \right)^{\text{T}}.\]

In the above expression X, Y and Z are the coordinates in Earth Centered Earth Fixed (ECEF) frame and the last two parameters are the receiver clock bias and the receiver clock drift. All the parameters are expressed in units of meters.

Now that the state vector is defined, we can move on by choosing the dynamic model. Before moving further let’s think a bit about this aspect. A static user doesn’t change his/hers position, therefore this means that over time the X, Y, Z coordinates remain the same! We only have to take care of how we model the dynamic behavior of the receiver’s clock, which is approximated to be:

\[\delta t_{R,k} = \delta t_{R,k-1} + \Delta T~\dot{\delta t}_{R,k-1},\]
\[\dot{\delta t}_{R,k} = \dot{\delta t}_{R,k-1}.\]

Having in view all of this information we can define the transition matrix (F) of the filter as:

\[\begin{split}\mathbf{F}_k = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & \Delta T \\ 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix}.\end{split}\]

We are almost done with the dynamic model elements. The only thing that we need now is the process noise matrix (Q). Because the process noise matrix contains the uncertainty we have in the dynamic model that we consider, we have to define it accordingly. In the static case we can assume that the user is not moving and that the receiver clock has some frequency and phase errors. In order to fully understand this reasoning, the interested reader is advised to check the following book: Introduction to Random Signals and Applied Kalman Filtering by Robert Grover Brown and Patrick Y. C. Hwang.Therefore, the process noise matrix is approximated to be:

\[\begin{split}\mathbf{Q}_k = \begin{pmatrix} 0~~~~& 0~~~~&0 & 0 & 0 \\ 0~~~~& 0~~~~& 0 & 0 & 0 \\ 0~~~~& 0~~~~& 0 & 0 & 0 \\ 0~~~~& 0~~~~& 0 & S_f+\frac{S_g~\Delta T^3}{3} & \frac{S_g~\Delta T^2}{2} \\ 0~~~~& 0~~~~& 0 & \frac{S_g~\Delta T^2}{2} & S_g~\Delta T \\ \end{pmatrix}.\end{split}\]

In the above expression the receiver clock related parameters are expressed as:

\[S_g \approx 2 \pi^2 h_{-2},\]
\[S_f \approx \frac{h_0}{2}.\]

The parameter h-2 and h0 are the Power Spectral Densities (PSD) of the random walk frequency noise and of the white noise, as defined in the suggested book above. Some typical values for a low quality Temperature Compensated Crystal Oscillator (TCXO) are 2e-20 and 2e-19 (in seconds). A practical advise before using this values is to take care that we are dealing with the parameters of a variance-covariance matrix and also that they have to be converted in units of meters (remember that we have expressed the receiver clock states in units of meters).

So basically we are done with the static user case. That’s great as we can move to the dynamic one!

Dynamic user

In the case of a dynamic user there are few aspects that one has to consider. Let’s start again by defining the new state vector:

\[\mathbf{x}_k = \left(X~~U~~Y~~V~~Z~~W~~\delta t_R~~\dot{\delta t}_R \right)^{\text{T}}.\]

We can already observe that we have three more parameters to estimate (U, V, W) which are the velocities on the X, Y and Z directions. If our state vector is modified (with respect to the static case) then our intuition will tell us that we need to define a new transition matrix and a new process noise matrix. Which is exactly what we are going to do next, therefore:

\[\begin{split}\mathbf{F}_k = \begin{pmatrix} 1 & \Delta T & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & \Delta T & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & \Delta T & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & \Delta T \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix}.\end{split}\]

For the process noise matrix we use the approach presented in the book of Robert Grover Brown and Patrick Y. C. Hwang ( Introduction to Random Signals and Applied Kalman Filtering ). Indeed, is the third time we refer to this book in the implemented PVT algorithms section, however you can trust us that is a very good one!

\[\begin{split}\mathbf{Q}_k = \begin{pmatrix} \frac{S_X~\Delta T^3}{3}& \frac{S_X~\Delta T^2}{2}& 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{S_X~\Delta T^2}{2}& S_X~\Delta T & 0 & 0 & 0 & 0 & 0 & 0 \\ 0~~~~& 0~~~~& \frac{S_Y~\Delta T^3}{3} & \frac{S_Y~\Delta T^2}{2} & 0 & 0 & 0 & 0\\ 0~~~~& 0~~~~& \frac{S_Y~\Delta T^2}{2} & S_Y~\Delta T & 0 & 0 & 0 & 0\\ 0~~~~& 0~~~~& 0 & 0 & \frac{S_Z~\Delta T^3}{3} & \frac{S_Z~\Delta T^2}{2} & 0 & 0\\ 0~~~~& 0~~~~& 0 & 0 & \frac{S_Z~\Delta T^2}{2} & S_Z~\Delta T & 0 & 0\\ 0~~~~& 0~~~~& 0 & 0 & 0 & 0 & S_f+\frac{S_g~\Delta T^3}{3} & \frac{S_g~\Delta T^2}{2} \\ 0~~~~& 0~~~~& 0 & 0 & 0 & 0 & \frac{S_g~\Delta T^2}{2} & S_g~\Delta T \\ \end{pmatrix}\end{split}\]

The parameters Sx, Sy and Sz are the spectral amplitudes that reflect the position random process. Unfortunately, setting their values is not as straigth forward as for the receiver clock states. We have to rely on what we call the tunning process which is modifying the values in Q and R experimentally (i.e., trial and error). Just as a note, this can be avoided by designing and implementing adaptive estimators. Who knows, maybe you (the reader) will decide to implement some nice ideas now that this possibility is enabled with GNSS Compare’s flexible framework.

Practical advise: When the observation matrix (H) is being built do consider that it’s size is defined in the following way: the number of rows is the number of measurements and the number of columns is the number of unknowns. Therefore when switching from the static case to the dynamic case, H changes also. We mention this just to be sure that a possible conceptual hiccup is avoided.

Filter tunning

Because at the moment we are dealing with a standard EKF and not an adaptive one this means that we have to assign values in the process noise matrix (Q) and in the measurement noise matrix (R) such that the filter is tunned to our situation.

Let’s start with the R matrix. We set R to be a diagonal matrix containing the variances of each pseudorange measurement. The measurement noise matrix being diagonal relies on the assumption that there is no cross-correlation between the measurements coming from different satellites ( an assumption that is not entirely represeting the reality, however it fits most of the applications ). Therefore, the diagonal elements of the R matrix are:

\[\mathbf{R}_{ii,k} = \sigma^2_{ii}.\]

To keep things relatively simple, we can assign the value for the sigma 10 meters ( don’t forget to square it before putting it in R ). Another assumption that is made is that the measurements received at the k-th epoch have equal variances ( ok, this assumption is not true at all ). However here is an idea for you, maybe you can try investigating some interesting measurement weigthing methods and then compare (the main keyword of the whole project) the results you get with our not so realistic assumption. Let the researcher within you thrive!

Let’s move to the Q matrix now. For this we present three tunning examples: static, pedestrian and dynamic.

Static tunning

For the static case we have already seen that we only have to take care about the process noise of the receiver clock states. So the values that we are assigning to the PSD of the random walk of the frequency noise and of the white noise are:

\[h_{-2} = 2e-20~c^2,\]
\[h_0 = 2e-19~c^2.\]

In the above we use the c notation for the speed of light.

Pedestrian tunning

Intuitively we should have used the EKF designed for a dynamic user in this situation. It would only make sense as a pedestrian changes his/hers position over time. However, one must take into account that the raw measurements delivered by the smartphone’s GNSS receiver are quite noisy and if there are no other means to detect the motion of the user (e.g., using an Inertial Measurement Unit) then estimating the velocities can make our results not soo accurate. Having this situation in view we have found a workaround: we use the EKF designed for a static user and we let some process noise for the X and Y coordinates ( unless one of our users is not Superman we are not that interested in the Z direction ). This means that we have the following Q matrix:

\[\begin{split}\mathbf{Q}_k = \begin{pmatrix} 0.2~~~~& 0~~~~&0 & 0 & 0 \\ 0~~~~& 0.2~~~~& 0 & 0 & 0 \\ 0~~~~& 0~~~~& 0 & 0 & 0 \\ 0~~~~& 0~~~~& 0 & S_f+\frac{S_g~\Delta T^3}{3} & \frac{S_g~\Delta T^2}{2} \\ 0~~~~& 0~~~~& 0 & \frac{S_g~\Delta T^2}{2} & S_g~\Delta T \\ \end{pmatrix}.\end{split}\]

The value 0.2 was chosen by trial and error and it fits a slow walking pedestrian. We hope that the name of the Java class PedestrianStaticExtendedKalmanFilter makes a little bit more sense now.

Dynamic tunning

Finally we have arrived at the final case regarding the tunning of the dynamic EKF. Again the following values were determined empirically:

\[S_X = S_Y = 0.8,\]
\[S_Z = 0.08.\]

Weighted Least Squares

GNSS Compare offers also the possibility to change the PVT estimator if the user whishes so. By not requiring knowledge about the dynamics, Weighted Least Squares (WLS) can be used to estimate the position using only the pseudorange measurements. However there are some drawbacks like: the quality of the estimations fully depends on the quality of the measurements and also the WLS requires a minimum number of measurements (typically 4 if we want to estimate the 3D position and the receiver clock bias).

Nevertheless is useful to have such an estimator as its behavior can be studied in real-time/post-processing in comparison with an EKF. And all this thanks to GNSS Compare!

Altough the pseudorange measurement model was presented in the EKF description we will do it one more time just for the sake of completion.

Pseudorange measurement model

The linearized code-based pseudorange measurement is:

\[PR_c - \rho_0 + \underbrace{(\delta t^S - d_{0,\text{ion}} - d_{0,\text{trop}})}_{Corr} = -\frac{X^S-X_0}{\rho_0}\Delta X-\frac{Y^S-Y_0}{\rho_0}\Delta Y-\frac{Z^S-Z_0}{\rho_0}\Delta Z+\delta t_R.\]

Let’s also express the unit line of sight vector and the position related unknowns as:

\[\mathbf{u} = \left[-\frac{X^S-X_0}{\rho_0},~~-\frac{Y^S-Y_0}{\rho_0},~~-\frac{Z^S-Z_0}{\rho_0} \right],\]
\[\delta\mathbf{r} = \left[\Delta X,~~\Delta Y,~~\Delta Z \right].\]

For n observed satellites we have the following measurement model:

\[\begin{split}\underbrace{\begin{pmatrix} PR^1_c - \rho^1_0 + Corr^1 \\ PR^2_c - \rho^2_0 + Corr^2\\ \vdots \\ PR^n_c - \rho^n_0 + Corr^n\\ \end{pmatrix}}_{\mathbf{z}} = \underbrace{\begin{pmatrix} \mathbf{u}^1 & 1\\ \mathbf{u}^2 & 1\\ \vdots & \vdots \\ \mathbf{u}^n & 1\\ \end{pmatrix}}_{\mathbf{H}} \underbrace{\begin{pmatrix} \delta \mathbf{r}^{\text{T}} \\ \delta t_R\\ \end{pmatrix}}_{\mathbf{x}}.\end{split}\]

To estimate the vector of unknowns (x) in the WLS sense, we proceed in the following way:

\[\hat{\mathbf{x}}_{\text{WLS}} = \left(\mathbf{H}^{\text{T}} \mathbf{W} \mathbf{H} \right)^{-1}\mathbf{H}^{\text{T}} \mathbf{W} \mathbf{z}.\]

Practical advise: In the WLS case, as the position is concerned, we are estimating the difference between the approximated position and the true position until this difference is below a certain threshold. We encourage you to check the WeightedLeastSquares class to see how this is handled.

Example of analysis

(To be updated with an example of analysis for the GNSS Dual-Frequency Android phones)

This section provides information about the scenarios in which GNSS Compare was tested and the PVT performance obtained from its estimation algorithms (e.g., Extended Kalman Filter). Furthermore, the analysis presented here serves also as an example of how GNSS Compare can be used for algorithmic performance assessment. For a preliminary PVT performance assessment the following scenarios were considered: Static user, Pedestrian user and Dynamic user.

With this section we would like to give you an idea of how GNSS Compare can be used. The application allows data logging (e.g., results of the PVT estimations) in different formats, like NMEA and a custom one. These files can be retrieved from the phone and then processed in your favourit programming environment for analysis. More details about the logging formats of GNSS Compare will be given soon.

Note 1: Please be aware that the results presented here are specific to the environment/time when they were generated and they cannot be interpreted in a general sense.

Note 2: The Extended Kalman Filters were initialized with the Android FINE location.

Static user

Let’s take a look at some details about this scenario:

  • Reference location: Latitude 52.16954469, Longitude 4.48089101, Altitude 55.48 m
  • Data collection duration: approximately 4 minutes
  • Enabled constellations: GPS, Galileo+GPS
  • Number of used satellites: 4 Galileo and 5 GPS

After the results of the PVT estimations were obtained from the logged files of GNSS Compare, they were projected in Google Earth as seen in the figure below for an initial analysis.

StaticGoogleEarth

In this scenario one can observe in the above figure that the computations based on Galileo+GPS are closer to the reference when compared with GPS only. In order to understand these aspects in a more detailed manner, the behavior of the errors with respect to the reference can be studied. The errors are computed based on the cartesian coordinates within the Earth Centered Earth Fixed (ECEF) frame.

staticEKFgps

The error evolutions for GPS only PVT are presented in the above figure and it can be directly observed that they are quite large and with a high variance. Let’s see what happens if we add Galileo in the processing.

staticGalGPS

For the case when the PVT is computed using both Galileo and GPS, the above figure shows improvements when compared with the solution based only on GPS.

Pedestrian user

This scenario is defined in the following way:

  • User dynamics: Walking pedestrian
  • Location: The European Space Research and Technology Centre (ESTEC)’s parking lot
  • Data collection duration: approximately 4 and half minutes
  • Enabled constellations: GPS, Galileo+GPS
  • PVT estimator: Extended Kalman Filter
  • Number of satellites: On average 3 Galileo and 8 GPS

As for this case there is no reference trajectory available the results are analyzed at the observed satellite level and at the projection of the estimated position in Google Earth.

pedestrianObsSV pedestrianGoogle

In the above figure the estimation of the trajectory that is based only on GPS does not follow too accurately the real pedestrian motion. However when both Galileo and GPS satellites are used together the position estimation is improved obtaining a pedestrian path closer to reality.

Dynamic user

And the last scenario has the following characteristics:

  • User dynamics: Cycling user
  • Location: ESTEC
  • Data collection duration: approximately 3 minutes
  • Enabled constellations: GPS, Galileo+GPS
  • PVT estimator: Extended Kalman Filter
  • Number of satellites: On average 4 Galileo and 8 GPS
bikeObsSV bikeGoogleEarth

Even with this rather simplistic analysis one can gain some interesting insights. We do hope that you have now a more clear idea about the possibilities that GNSS Compare can open!