Measurement quality from affordable GNSS receivers

Measurement quality from affordable GNSS receivers#

In this section we will assess some basic differences in terms of GNSS observables between a geodetic-grade and an affordable receiver.

Using the roktools library, the RINEX data from the receivers will be stored into a pandas DataFrame (df for short) for manipulation.

# DataFrame with data from a geodetic receiver
df_geodetic = rinex.to_dataframe('../assets/SEYG00SYC_R_20140581500_05H_01S_MO.rnx')

# DataFrame with an affordable receiver
df_afford = rinex.to_dataframe('../assets/MTIC00ESP_R_20191221131_05H_01S_MO.rnx')

The DataFrame can be consider a CSV file of sorts, where each row has a time tag, satellite, observables and other fields that will be explained and used later in the notebook. A preview of the contents can be obtained with the head method

df_geodetic.head()

	epoch	constellation	sat	channel	signal	range	phase	doppler	snr
0	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	1X	E121X	2.492298e+07	1.309712e+08	2196.344	47.0
1	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	5X	E125X	2.492299e+07	9.780326e+07	NaN	47.0
2	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	7X	E127X	2.492299e+07	1.003546e+08	NaN	46.3
3	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	8X	E128X	2.492299e+07	9.907894e+07	NaN	52.2
4	2014-02-27 15:00:00	ConstellationId.GALILEO	E19	1X	E191X	2.500360e+07	1.313950e+08	-2145.535	48.7

Observable types#

Once we have loaded the RINEX files into DataFrames, we can perform some basic checks on the differences between geodetic and affordable GNSS data.

The following example gives you the channels tracked by the receiver for each GNSS constellation.

The channel corresponds to the last two characters (i.e. band and attribute) of the RINEX observation code (see Section 5.2.17). For instance 1C for GPS means the observables obtained with the C/A tracking at the L1 frequency.

To know the observables for constellation, we will use the groupby method of pandas.

# We will group the data using various criteria
columns = ['constellation', 'channel']

# Use groupby() to group by the two columns and apply unique()
unique_combinations = df_geodetic.groupby(columns).size()

# Print the unique combinations (along with the number of samples for each
# tracking channel)
print(unique_combinations)

constellation            channel
ConstellationId.GALILEO  1X         18910
                         5X         18910
                         7X         18910
                         8X         18910
ConstellationId.GPS      1C         87774
                         2W         87774
                         2X         87774
                         5X         87774
ConstellationId.GLONASS  1C         54304
                         1P         54304
                         2C         54304
                         2P         54304
ConstellationId.SBAS     1C         25134
dtype: int64

We can now perform the same for the affordable receivers

unique_combinations = df_afford.groupby(columns).size()

print(unique_combinations)

constellation            channel
ConstellationId.BEIDOU   2I          95280
ConstellationId.GALILEO  1C         106261
                         7Q         106261
ConstellationId.GPS      1C         138731
                         2L         138731
ConstellationId.GLONASS  1C         111380
                         2C         111380
dtype: int64

In these examples, the following basic differences are observed:

the geodetic receiver tracks various frequencies (GPS L1/L2/L5, Galileo E1/E5a/E5b/E5, …) whereas the affordable receiver tracks typically two frequencies (GPS L1/L2, Galileo E1/E5b, …)
the affordable receiver does not attempt to track encrypted codes (i.e. GPS P code) by means of e.g. z-tracking loops. Only civilian codes (e.g. GPS L2C) are used.

Some other strenghts of affordable receivers:

Availability of SNR and Doppler measurements (not always available in 30s or high rate CORS data)
High rate up to 0.1s (or even higher) available for affordable measurements

Code noise: Detrended code-minus-carrier#

The observable code noise of a GNSS receiver can be estimated using the code-minus-carrier combination. This section illustrates how to estimate to check some basic differences between receiver types

First we will import several methods that will be required to process data (will be illustrated later)

# Import the modules from the custom 'source' package of this book
from gnss.observables import compute_code_minus_carrier
from gnss.edit import mark_time_gap, detrend, remove_mean

The steps to be followed to compute the unbiased detrended Code minus carrier are detailed in the following steps (taking as example the Geodetic receiver and then applying it to the Affordable receiver)

Edit the data and find phase breaks such as cycle slips. In this example, since receivers already provide with Loss-of-Lock Indicator (LLI), we will only mark phase breaks due to data time gaps

df_geodetic_cmc = mark_time_gap(df_geodetic)

Proceed to compute the code minus carrier (CMC) observable. Each row of the DataFrame contains the range and phase, so there will be as many CMC observables as rows. A note of caution: the phase is usually expressed in cycles, therefore we will need to get the wavelength from the tracking channel in order to be consistent with the units

df_geodetic_cmc = compute_code_minus_carrier(df_geodetic_cmc)

df_geodetic_cmc.head()

	epoch	constellation	sat	channel	signal	range	phase	doppler	snr	slip	wl	cmc
0	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	1X	E121X	2.492298e+07	1.309712e+08	2196.344	47.0	False	0.190294	-13.752709
1	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	5X	E125X	2.492299e+07	9.780326e+07	NaN	47.0	False	0.254828	-23.022230
2	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	7X	E127X	2.492299e+07	1.003546e+08	NaN	46.3	False	0.248349	-20.960017
3	2014-02-27 15:00:00	ConstellationId.GALILEO	E12	8X	E128X	2.492299e+07	9.907894e+07	NaN	52.2	False	0.251547	-22.530810
4	2014-02-27 15:00:00	ConstellationId.GALILEO	E19	1X	E191X	2.500360e+07	1.313950e+08	-2145.535	48.7	False	0.190294	-36.599890

from helpers import compute_elapsed_seconds

# Plot a CMC for a channel and satellite
signal = 'E121X'
df_sample = df_geodetic_cmc[df_geodetic_cmc['signal'] == signal]
plt.close()
plt.plot(compute_elapsed_seconds(df_sample['epoch'])/3600, df_sample['cmc'], ',k')
plt.xlabel('Elapsed hours')
plt.ylabel('CMC [m]')
_ = plt.title(f'Code-minus-carrier for {signal}')

Because the CMC contains twice the ionosphere, which is a nuisance parameter for us (because we are interested in the noise). We will proceed to remove its contribution with a simple detrending, based on a rolling window

n_samples = 5
df_geodetic_cmc = detrend(df_geodetic_cmc, 'cmc', n_samples)

# Plot the detrened CMC for a channel and satellite
signal = 'E121X'
df_sample = df_geodetic_cmc[df_geodetic_cmc['signal'] == signal]
plt.close()
plt.plot(compute_elapsed_seconds(df_sample['epoch'])/3600, df_sample['cmc_detrended'], ',k')
plt.xlabel('Elapsed hours')
plt.ylabel('Detrended CMC [m]')
plt.ylim(-1, 1)
_ = plt.title(f'Detrended Code-minus-carrier for {signal}')

# Analysis for Ublox receiver
df_afford_cmc = mark_time_gap(df_afford)
df_afford_cmc = compute_code_minus_carrier(df_afford_cmc)
df_afford_cmc = detrend(df_afford_cmc, 'cmc', n_samples)

Let’s compute now an estimate of the code noise for the geodetic grade and the affordable receiver, for a specific band and constellation.

from helpers import compute_rms

# GLONASS is excluded as the slot number to know the frequency (and hence the wavelength)
# is missing
condition1 = df_geodetic_cmc['constellation'] == ConstellationId.GPS
condition2 = df_geodetic_cmc['channel'] == '2W'
#condition2 = df_geodetic_cmc['channel'] == '1C'
df_tmp_g = df_geodetic_cmc[condition1 & condition2]

condition1 = df_afford_cmc['constellation'] == ConstellationId.GPS
condition2 = df_afford_cmc['channel'] == '2L'
#condition2 = df_afford_cmc['channel'] == '1C'
df_tmp_a = df_afford_cmc[condition1 & condition2]

noise_samples_geodetic = df_tmp_g['cmc_detrended']
noise_samples_afford = df_tmp_a['cmc_detrended']

rms_geodetic = compute_rms(noise_samples_geodetic)
rms_afford = compute_rms(noise_samples_afford)

print(f'Geodetic   receiver code noise: {rms_geodetic:.2} m')
print(f'Affordable receiver code noise: {rms_afford:.2} m')

Geodetic   receiver code noise: 0.34 m
Affordable receiver code noise: 0.32 m

Observable noise: Detrended LI#

In addition, to the estimation of the code, for multifrequency receivers, the geometry-free (ionospheric) combination can be used to estimate the noise of the phase observables

from gnss.observables import compute_geometry_free

constellation = ConstellationId.GPS
channel_a = TrackingChannel.from_string('1C')
channel_b = TrackingChannel.from_string('2W')

# Let's reuse the previous dataframe that has been already edited (time gaps marked)
df_geodetic_li = compute_geometry_free(df_geodetic_cmc, constellation, channel_a, channel_b)

The ionospheric combination (LI) contains the ionospheric delay, constant terms due to the phase ambiguities of the carrier phases used to build the combination as well as the phase hardware biases and noise of the phase observable

An example of LI for a particular satellite is shown in the Figure below

sat = Satellite.from_string('G04')

condition1 = df_geodetic_li['sat'] == sat
df_tmp = df_geodetic_li[condition1]

t = compute_elapsed_seconds(df_tmp['epoch'])
li = df_tmp['li_m']

plt.close()
plt.plot(t, li, '.')
plt.xlabel('Time [seconds]')
plt.ylabel('LI [m]')
plt.title(f'LI combination for {sat}, Geoetic receiver (SEY1)')

Text(0.5, 1.0, 'LI combination for G04, Geoetic receiver (SEY1)')

import pandas as pd

def detrend_li_combination(df:pd.DataFrame):
    """
    Perform a detrending of the LI combination (if present in the input dataframe)
    """

    df['slip_li'] =  df['slip_a'] | df['slip_b']
    df['arc_id'] = df.groupby(['signal_a', 'signal_b'])['slip_li'].transform('cumsum')
    trend = df.groupby(['signal_a', 'signal_b', 'arc_id'])['li_m'].transform(lambda x: x.rolling(n_samples).mean())
    df['li_trend'] = trend
    df['li_detrended'] = df['li_m'] - df['li_trend']

detrend_li_combination(df_geodetic_li)

condition1 = df_geodetic_li['signal_a'] == 'G041C'
condition2 = df_geodetic_li['signal_b'] == 'G042W'

df_sample = df_geodetic_li[condition1 & condition2]
plt.close()
plt.plot(df_sample['epoch'], df_sample['li_detrended'], '.')
plt.ylim(-0.05, 0.05)

(-0.05, 0.05)

We can proceed with the same analysis for the affordable receiver, taking into account that the observables are slightly different

channel_a = TrackingChannel.from_string('1C')
channel_b = TrackingChannel.from_string('2L')

df_afford_li = compute_geometry_free(df_afford_cmc, constellation, channel_a, channel_b)

detrend_li_combination(df_afford_li)

condition1 = df_afford_li['signal_a'] == 'G241C'
condition2 = df_afford_li['signal_b'] == 'G242L'

df_sample = df_afford_li[condition1 & condition2]

plt.close()
plt.plot(df_sample['epoch'], df_sample['li_detrended'], '.')
plt.ylim(-0.05, 0.05)

(-0.05, 0.05)

The detrended time series can be used to compute the noise figures of the carrier phase for both types of receivers

# Remove outliers for the statistics
noise_li_afford = df_afford_li['li_detrended'][(df_afford_li['li_detrended'] > -1) & (df_afford_li['li_detrended'] < 1)]
noise_li_geodetic = df_geodetic_li['li_detrended'][(df_geodetic_li['li_detrended'] > -1) & (df_geodetic_li['li_detrended'] < 1)]

print(f"Phase Noise Affordable receiver: {np.std(noise_li_afford)*100:.0f} cm")
print(f"Phase Noise Geodetic   receiver: {np.std(noise_li_geodetic)*100:.0f} cm") 

Phase Noise Affordable receiver: 1 cm
Phase Noise Geodetic   receiver: 2 cm

As it can be seen, the differences between receivers is not substantial. This can be further confirmed with the actual distribution of the noise samples, which follow a very similar Gaussian pattern for both receivers

plt.close()
bins = np.linspace(-0.05, 0.05, 200)
_ = plt.hist(df_geodetic_li['li_detrended'], bins=bins, histtype='step', label='Geodetic')
_ = plt.hist(df_afford_li['li_detrended'], bins=bins, histtype='step', label='Affordable')
plt.legend()
plt.title('Histogram of carrier phase noise (estimated with LI)')
plt.xlabel('Code phase error [m]')
plt.ylabel('Count')

Text(0, 0.5, 'Count')

Measurement quality from affordable GNSS receivers

Contents

Measurement quality from affordable GNSS receivers#

Observable types#

Code noise: Detrended code-minus-carrier#

Observable noise: Detrended LI#