PVT-SPS Documentation

Documentation of PVT-SPS.
Nicolas Padron. Copyright 2024.
https://www.spacewayfinder.com
pvt-sps snapshot.

Software Description

Follow the GNSS standalone navigation layer (GPS L1 SPS) as from RINEX measurements, to satellite position calculation following ICD algorithms, and finally calculate receiver PVT.

Explore different modes that can run in parallel and finally get statistics and plots comparing:

  • Fix with and without atmospheric corrections (Collins for Troposphere and Klobuchar for Ionosphere)

  • LS vs WLS vs MMSE estimation

  • Elevation and/or CN0 masking, avoiding satellites with low elevation angles or low CN0 values.

  • For WLS: Different weights and customizable wls.py so that user can define its own weighting function, following available parameters in WeightAlgorithmsParameters class (see Documentation).

  • For MMSE tunable prior and measurment matrices by a constant to give more or less trustness to either prior or current measurements, balancing the solution smoothness.

  • For LS/WLS: satellite exclusion by choosing the solution that results in the lowest HDOP by excluding one satellite at a time.

Predefined modes in /src/config.py:

  • MODES_ATMOSPHERIC_COMPARISON: compare LS with and without atmospheric corrections.

  • MODES_ELEVATION_FILTERING: compare LS with elevation and CN0 maskings.

  • MODES_LSE_VS_FE: compare LS with and without excluding satellites.

  • MODES_LSE_VS_WLS: compare LS vs WLS with different weights.

  • MODES_LSE_VS_MMSE: compares LS vs MMSE

Usage Description

The execution is controlled by the config.py script in /src folder. Here the user can choose what algorithm from wls.py to use, as well as define its own based on the available parameters in WeightAlgorithmsParameters class.

Steps to do before running the code:

  • 1: Remember to create a virtual environment and install the libary versions in requirements.txt.

  • 2: If no navigation file is provided, is will be downloaded from CDDIS website. For this to succeed, the user should enter the CDDIS account details in /src/.netrc file.

  • 3: Explore the /src/config.py file and run the code (from root directory) as python .\src\main.py

Folder Architecture

  • /src: Python source files for RINEX reading, PVT processing and plotting.

  • /data: default folder to store RINEX observation and navigation files.

Modules

Configuration

class config.Config

Main class for Config.

Parameters:

  • modes – defined Modes

  • imgsdir – name of directory to create within imgs/ where result plot images and CSV with statistics will be stored. The directory will be created, and if it alrady exists an error message will be triggered to either remove or rename it to avoid mistakenly ovewriting previous reaults.

  • fileObs – RINEX observation file in /data/obs folder

  • fileNav – RINEX navigation file in /data/nav folder. If empty, then the program will download the corresponding navigation file for the day in the observation file from the CDDIS server.

  • use_pvt_ref_as_origin – Bool. If True, it will use the user-defined PVT reference or RINEX XYZ reference for ENU calculation. If False it uses the 1st PVT estimated.

  • pvt_reference – user-defined PVT reference in ECEF frame. If Latitude, Longitude and Height are used, they can be converted with Rotation.conversion_wgs84_2_ecef()

  • calc_el_az_from_ref – Bool. If True, decides to calculate satellite elevation and azimuth angles from PVT reference, if False it will use the last PVT estimation.

  • plotsExcludeDOP – Bool. If True then the result plots and statistics will exclude all fixes with HDOP bigger than the HDOP’s mean + half sigma for each mode.

  • maxEpochsPercentage – Number from 0 - 100. Is the percentage of total epochs that will be processed. Example: In a static 30s daily observation file with a total of 2880 epochs, we can decide to process only the first 10% by setting the value to 10.

NOTE on fileNav: when downloaded from server it is stored in /data/nav. If when reading the downloaded file it triggers an error saying that it cannot open the file, try to manually uncompress it and write the uncompressed filename in fileNav.

NOTE on use_pvt_ref_as_origin: if True, it will use the user-defined PVT reference, only if the RINEX XYZ reference is not present. If the RINEX XYZ reference is present, it will use that one instead of the user-defined. Program assumes that reference in RINEX, if present, is correct. If False, it will use the 1st calculated PVT as reference. This corresponds to the 1st calculated mode: if LS and MMSE are run in that order, the reference will be the LS result.

NOTE on calc_el_az_from_ref: On static locations, we can calculate the satellite elevation and azimuth from reference rather than from receiver estimates.

class config.Mode(_name='', _use_ls=True, _ls_use_last_fix=False, _ls_weight_algorithm=<bound method WeightAlgorithms.algo_linear of <wls.WeightAlgorithms object>>, _ls_max_iterations=15, _ls_exclusion=False, _elev_mask=5, _cn0_mask=10, _enable_iono_corr=True, _enable_tropo_corr=True, _mmse_px_mpy=1, _mmse_r_mpy=1)

Main class defining the user defined modes.

Parameters:

  • _name – Default ‘’. Defines the mode name. Plots will take this name as label.

  • _use_ls – Default True. Bool to determine if using LS (=True) or MMSE (=False).

  • _ls_use_last_fix – Default False. Bool to determine if pick LS from last iteration solution (=True), or start iterations from scratch (=False).

  • _ls_weight_algorithm – Default WeightAlgorithms().algo_linear. This is the function from class WeightAlgorithms defining the weights. LS, is the corresponding WLS, and in MMSE, these weights are applied in the measurement noise matrix R.

  • _ls_max_iterations – Default 15. Max number of iterations in LS.

  • _ls_exclusion – Default Fasle. If True, in LS it will run the LS process N times, each time using N-1 measurements (excluding one), and finally pick the solution with smallest GDOP.

  • _elev_mask – Default 5. Elevation mask so that only SVs with higher elevation will be considered as valid. Unit degrees.

  • _cn0_mask – Default is 10. CN0 mask so that only SVs with higher CN0 will be considered as valid. Units dB/Hz.

  • _enable_iono_corr – Default True. Boolean to enbale (=True) or disable (=False) ionospheric correction.

  • _enable_tropo_corr – Default True. Boolean to enbale (=True) or disable (=False) tropospheric correction.

  • _mmse_px_mpy – Default 1. If MMSE (_use_ls = False) then the prior covariance matrix Px is multiplied by this factor to ensure more lr less trustness on prior.

  • _mmse_r_mpy – Default 1. If MMSE (_use_ls = False) then the measurement covariance matrix R is multiplied by this factor to ensure more lr less trustness on mesaurements.

Weighted LS

class wls.RobustWeightFunctions

Class to define the Weight Functions for robust estimation

weight_function_huber(alpha, isabs=False)

Huber function following Table 1 in reference

Parameters:

  • res – residuals

  • alpha – alpha parameter for the weighting function.

weight_function_tuckey(alpha)

Tukey bi-weight function following Table 1 in reference

Parameters:

  • res – residuals

  • alpha – alpha parameter for the weighting function.

class wls.WeightAlgorithms(weight_function=<function RobustWeightFunctions.weight_function_huber>)

Class to define the algorithms to use for weighted LS

algo_cn0_elev(params)

Based on CN0 and Elevation

Following equation 8 in https://gssc.esa.int/navipedia/index.php?title=Best_Linear_Unbiased_Minimum-Variance_Estimator_(BLUE)

with a = 0, b = CN0 and c = 1.

Parameters:

params – input parameters as defined in WeightAlgorithmsParameters

Return weights:

Nx1 array of weights, where N is the number of available measurements.

algo_cn0_elev_tuned(params)

Tuned weighting function based on elevation and CN0.

Parameters:

params – input parameters as defined in WeightAlgorithmsParameters

Return weights:

Nx1 array of weights, where N is the number of available measurements.

algo_linear(params)

No weight at all, returns Nx1 vector of 1s

algo_m_estimation(params)

Robust estimation (simplification) based on IRLS (M-Estimation) from:

Akram, M.A.; Liu, P.; Wang, Y.; Qian, J. “GNSS Positioning Accuracy Enhancement Based on Robust Statistical MM Estimation Theory for Ground Vehicles in Challenging Environments.” Appl. Sci. 2018, 8, 876. https://doi.org/10.3390/app8060876

Parameters:

params – input parameters as defined in WeightAlgorithmsParameters

Return weights:

Nx1 array of weights, where N is the number of available measurements.

class wls.WeightAlgorithmsParameters(nIter, _numMeas, _cn0, _el, _az, _pRes)

Class defining the avalable parameters used weighting functions

Receiver Processing

class receiver.Ionosphere

Class to handle Ionosphere delay calculation following Klobuchar Ionospheric Model https://gssc.esa.int/navipedia/index.php?title=Klobuchar_Ionospheric_Model

calc_iono_delay(alphan, betan, lat, lon, elev, azim, numMeas)

Main function to calculate the ionospheric delay

Parameters:

  • t_sow – second of week

  • alphan – alpha parameter for ionospheric model

  • betan – beta parameter for ionospheric model

  • lat – receivr latitude

  • lon – receiver longitude

  • elev – array with SVs elevation angles

  • azim – array with SVs azimuth angles

  • numMeas – number of available satellites

Return ionoDelay:

ionospheric delay in meters

class receiver.RxPVTCalculation(sys_)

Main class for Rx PVT processing

_calc_atmospheric_corrections(numMeas, enable_iono_corr, enable_tropo_corr)

Private function to handle Atmospheric Corrections: both Ionospheric and Tropospheric

Parameters:

  • numMeas – number of available measurements

  • enable_tropo_corr – boolean from user-specified mode indicating if tropospheric correction should be done.

Paramm enable_iono_corr:

boolean from user-specified mode indicating if ionospheric correction should be done.

Return iono:

array with ionospheric delays in meters for each SV.

Return tropo:

array with tropospheric delays in meters for each SV.

_calculate_ls(jacobian, weights, residuals)

LS or WLS fitting calculation

Parameters:

  • jacobian – Nx4 jacobian matrix.

  • weights – weights to apply.

  • residuals – pseudorange residuals.

_check_singular(jacobian)

Check geometry to avoid linear dependent columns among the Jacobian rows, thus avoid triggering singular matrix error.

Parameters:

jacobian – Nx4 jacobiann matrix

Return validGeometry:

boolean stating if the jeometry is valid (matrix is not singular, i.e. determinant is nonzero).

_prepare_weights(ls_weight_algorithm, dataSys, pRes)

Prepare information to be used by weight functions for WLS, and calculate weights.

Params ls_weight_algorithm:

function in modes defined by the user stating what weight function in wls module use.

Params dataSys:

system data, i.e. GPS observables.

Params pRes:

pseudorange residuals

Return weights:

array containing the weights for each measurement.

_process_jacobian()

Compute Jacobian, i.e. geometry matrix.

Return jacobian:

jacobian Nx4 matrix, where N is the number of measurements.

_process_ls(dataSys, mode)

Main fuction for processing linearized LS loops

Parameters:

  • dataSys – system data, i.e. GPS, containing system observables.

  • mode – user defined mode to run which applies to LS.

Return last_iter_ok_mode:

boolean stating if the iteration is valid (successful position calculation) or not (diverging calculation, singular matrix).

_process_mmse(dataSys, mode)

Main fuction for processing linearized MMSE

Parameters:

  • dataSys – system data, i.e. GPS, containing system observables.

  • mode – user defined mode to run which applies to LS.

Return last_iter_ok_mode:

boolean stating if the iteration is valid (successful position calculation) or not (diverging calculation, singular matrix).

_process_pseudo(pseudos)

Process pseudorange, calculate pseudorange estimation and residuals

p = r - t_sv + t_rx + iono + tropo + multipath + noise * geometric range r is calculated as the euclidean distance between estiamted Rx and Tx positions. * time delays are represented in meters * multipath and noise assumed as generic error * iono and tropo corrected after the first position estimation, meanwhile they are 0 * t_rx and xRx are estimations, thus they change in the iterative algorithm until convergence.

Parameters:

pseudos – array containing available satellites’ pseudoranges

Return pseudorangeResiduals:

pseudorange residuals as the difference between estimated and meausred pseudoranges.

_run_mode(mode)

Processing for every mode

Parameters:

mode – user defined mode to run which applies to LS.

Return is_last_iter_ok_mode:

bool stating if the calculatioon is valid or not coming from _process_ls or _process_mmse .

process()

Main interface function to call for process Receiver PVT. All modes will be executed here, and the bool parameter is_last_iter_ok_mode coming from _run_modes is evaluated to ensure that all modes are valid and can be fairly comparted.

class receiver.Troposphere

Class to handle Troposphere delay following [Collins, 1999] https://gssc.esa.int/navipedia/index.php/Tropospheric_Delay “Example of Tropospheric model for Standard Point Positioning”

calc_tropo_delay(lat, hei, elev, doy, numMeas)

Main function for Tropospheric Delay calculation. :param t_sow: second of week :param lat: receiver latitude :param hei: receiver height :param elev: array with SVs elevation angles :param doy: day of year :param numMeas: number of measurements present

Return tropoDelay:

array with each SV tropospheric delay in meters

class receiver._Helpers

Helper class for the receiver module

check_valid_sv_masks(mode)

Check the SV validity depending on the chosen mode.

Remember satellite module gives all available measurements, and then receiver module choses what mode to use.

Such modes can vary, for example, elevation and cn0 masks.

Parameters:

  • dataSys – measurements data from specific system, i.e. GPS

  • mode – chosen mode following Modes class specified in config module by the user.

Return svResultHandler:

handler with satellites which are available following the chosen mode parameters.

Satellite Processing

class satellite.SvPVTCalculation(sys_)

Main class to handle SV PVT calculation: position in ECEF frame and time offset.

calc_pos_and_clkcorr(t_sow, ephSv)

Estimation of satellite position in ECEF and clock correction.

“20.3.3.4.3 Algorithms for Ephemeris Determination” on IS-GPS-200M

Basically, Table IV is followed. The Keplerian paramaters are utilized to determine the SV position and clock correction. This is done after the Tx time estimation and initial clock correction (which did not account for relativistic effect).

In fact, relativistic effect depends on Eccentric Anomaly Ek, which is determined during this algorithm previous to the SV position estimation.

Note that, the term ‘t’ in equation ‘tk = t - toe’ in Table IV on ICD document, is taken as ‘timeSvCorrected’ in this algorithm, and corresponds to the initial Tx time estimation done in function calcTxTime(), which is then corrected by the relativistic term in every loop iteration in calc_pos_and_clkcorr() since the relativistic term is estimated together with Ek.

This way of handling the relativistic term is also commented on function calcTxTime(). NOTE: This functions is a wrapper to Numba (C-compiled) function to accelerate processing.

Parameters:

  • t_sow – second of week

  • ephSv – satellite ephemerides struct

Returns:

concatenated SV postion and time error

calc_tx_time(t_sow, ephSv, pseudo, sigTypeFreq)

Calculate Tx time and clock correction except for relativistic effect.

“20.3.3.3.3.1 User Algorithm for SV Clock Correction” on IS-GPS-200M

NOTE: This functions is a wrapper to Numba (C-compiled) function to accelerate processing.

Parameters:

  • t_sow – second of week

  • ephSv – satellite ephemerides struct

  • pseudo – pseudorange

  • sigTypeFreq – signal type frequency, currently only GPS supported.

process()

Main function to process satellite raw pseudoranges, extimate satellite clock offset and position in ECEF plane.

class satellite._Helpers

Helper class for the Numba compiled functions from satellite module.

NOTE: this class depends on ephSv which is originally a DataFrame. It is converted to NumPy and as an array and the positions corresponding to every label are taken.

This is needed because Numba does not support Pandas. The indexes for every DataFrame column can be easily seen with, for example: np.column_stack([ephSv.columns.to_numpy(), np.arange(ephSv.columns.size)]) in the process() function before converting ephSv to numpy.

calc_pos_and_clkcorr(ephSv, timeSv, deltaTimeSv, deltaTimeRel, muGravConstant, forceF, omegaEDot, secondsInHalfWeek, secondsInWeek, epsilonErrMag10)

Numba compiled function to calculate satellite position in ECEF frame and clock error.

Parameters:

  • t_sow – second of week

  • ephSv – satellite ephemerides struct

  • timeSv – first estimation of satellite time offset as calculated in calc_tx_time.

  • deltaTimeSv – first estimation of satellite delta time offset as calculated in calc_tx_time.

  • deltaTimeRel – first estimation of satellite relativistic time offset as calculated in calc_tx_time.

  • muGravConstant – Earth gravitational constant, since Constants class is not accessible from Numba.

  • forceF – force constant, since Constants class is not accessible from Numba.

  • forceF – force constant, since Constants class is not accessible from Numba.

  • omegaEDot – Earth spin rate, since Constants class is not accessible from Numba.

  • secondsInHalfWeek – seconds in half week, since Constants class is not accessible from Numba.

  • secondsInWeek – seconds in week, since Constants class is not accessible from Numba.

  • epsilonErrMag10 – constant to determine when the loop estimation of Ek is enough, since Constants class is not accessible from Numba.

Return xyzECEF:

SV position in ECEF frame.

Return svClkCorr:

total SV clock correction.

calc_sv_elev_and_az(xRx, Rot)

Numba compiled function to calculate satellite elevation and azimuth angles.

Parameters:

  • xSv – SV position in ECEF frame.

  • xRx – Receiver position in ECEF frame.

  • Rot – Rotation matrix depending on receiver latitude, longitude and height.

Return elevation:

SV elevation angle in radians.

Return azimuth:

SV azimuth angle in radians.

calc_tx_time(pseudo, ephSv, speedOfLight)

Numba compiled function to calculate satellite time offset.

Parameters:

  • t_sow – second of week

  • pseudo – pseudorange

  • ephSv – satellite ephemerides struct

  • speedOfLight – speed of light constant, since Constants class is not accessible from Numba.

Return timeTxFromSv:

the time offset from the SV.

Return deltaTimeRel:

relativistic delta time, which estimation is enhanced in calc_pos_and_clkcorr.

Return deltaTimeSv:

total delta time as stated in the ICD documentation.

Manager

class manager.HelpersDataframe

Helper class to handle data dataframe indexing for receiver results

class manager.Interface

Interface class to handle progress bar

class manager.Manager

Manager class which encompasses Monitor, Processing, Interface and Time classes to control the processing flow. This class also invokes the data configuration (user defined modes) and data loading (from OBS or NAV RINEX files into respective modules).

class manager.Monitor

Monitor class to store sanity checks.

class manager.Processing

Class read and store the data to be used in processing: observation and ephemerides modules read from their respective RINEX files.

class manager.ResultHandlers(modes)

Class containing the results for satellite and receiver modules.

class manager.RxResultHandler(modes)

Class for handling the receiver module results

class manager.SvResultHandler(_isSvValid=[], _validSVsIx=[], _pos=[], _measIdx=[], _iono=[], _tropo=[], _numMeas=[])

Class for handling the satellite module measurements under usage

class manager.Time(_t, _t_ix, interval=0)

Class for to map observation file second of week, day of year and corresponding epoch index.

manager.manager = <manager.Manager object>

Manager is accessible everywhere in the program, so it is defined here in order to be common for all modules and files. Variable name is manager.

class manager._Helpers

Class to handle the JIT accelerated processing for Ephemerides module and other helper functions

class manager._EphHandlerModule(dataNav)

Private class to handle ephemerides, check the validity and pick the closest one in past time.

class manager._MeasSystemBand(idx, dataObs)

Private class to handle measurements per band

class manager._ObsHandlerModule(dataObs)

Private class to handle measurmmenets from observation file. Currently only GPS L1 supported.

Utils

class utils.Constants

General use constants

class utils.Rotation

Class to handle rotations and convertions between planes.

conversion_ecef_2_wgs84()

ECEF to WGS84 conversion.

Following Table 2.1 “Determination of Geodetic Height and Latitude in Terms of ECEF Parameters”

Parameters:

xyz – ECEF coordinates

Return llh:

WGS84 coordinates

conversion_wgs84_2_ecef()

WGS84 to ECEF conversion

Follwoing 2.2.3.2 Conversion from Geodetic Coordinates to Cartesian Coordinates in ECEF Frame

Parameters:

  • llh – WGS84 coordinates

  • xyz – ECEF coordinates

matrix_wgs84_2_enu()

Rotation matrix from ECEF to ENU, needs WGS84 coordindates.

Parameters:

llh – WGS854 coordinates

Return rotmat:

3x3 rotation matrix

class utils.SysInfo

Class storing system constants to be used in calculation. Currently only GPS L1 supported.

class utils.TimeConverters

Utils class to handle time calculations: day of year and second of week.

utils.dot_mat_vec(X, y)

Matrix-vector multiplication.

Numba functions to calculate X @ y, since the @ operator is not allowed in Numba.

Parameters:

  • X – MxL matrix

  • y – Lx1 vector

Returns:

Mx1 array containing product multiplication

utils.dot_vec_vec(x, y)

Vector-vector dot multiplication.

Numba functions to calculate x @ y, since the @ operator is not allowed in Numba.

Parameters:

  • x – Lx1 vector

  • y – Lx1 vector

Returns:

dot product

Read Rinex

class readRinex.ReadRinex(fileObs_, fileNav_)

Main class to read rinex

get_rinex_data()

Get read observation and navigation data structures.

read_rinex_file()

Main function to read RINEX files called by Manager.

Observation and navigation files are read, and if navigation file is not provided

the program will retrieve it from the CDDIS server.

Return retcode:

bool stating if the read was successful or not.