Author: rtklibexplorer

Loosely coupled sensor fusion: Part II

[Update 12/18/25: I found a couple more small errors in the GNSS_IMU code since publishing this and so am updating this post with the results from the latest code. The results between the two code sets are now nearly identical for this mid-range IMU dataset.]

In my previous post, I introduced GNSS_IMU, a Python tool I’ve recently written for demonstrating loosely coupled GNSS/IMU sensor fusion based on Paul Groves’ textbook and sample Matlab code. Shortly after publishing that post, I discovered a related project from Wuhan University. It is also a demonstration of loosely coupled GNSS/IMU sensor fusion shared on GitHub. The project is called KF-GINS and has two repositories, one written in C and one in Matlab. They recently published a paper in the GPS Solutions journal introducing the project, for which they have also shared a publicly accessible copy in their GitHub repo.

It appears to be well-written code and includes a detailed document describing the theory behind the implementation, along with a variety of other supporting documents, some in English and some in Chinese. It is intended for use with IMUs ranging from consumer-grade all the way to navigation-grade, and hence has significantly greater mathematical complexity than the GNSS_IMU code, which was intended to be more lightweight and focused on consumer-grade IMUs. I definitely recommend their libraries as a great resource for anyone interested in this topic.

In addition to the code and documentation, the Matlab library includes three sample datasets, two of which use the ADIS16465 IMU from Analog Devices, a mid-range MEMS IMU. What is especially useful about these datasets is that they also include ground truth from a high-end INS system. For comparison, the ADIS16465 is available on DigiKey for about $621, over 100 times the cost of the TDK ICM45686 IMU I used during development of my code.

While the GNSS_IMU code wasn’t necessarily intended for use with this class of IMU, I thought it would be an interesting exercise to compare how it performed versus the more capable and more complex KF-GINS code, and that using one of these datasets would be a good way to make this comparison. In particular, their Dataset 2 seemed like a good choice, since it includes not just missing GNSS data, but also a mix of poor-quality GNSS data in the gaps. This is more realistic than simply creating pseudo-gaps by deleting intervals of good GNSS data, but it is also more difficult to evaluate without a high-end INS system to generate ground truth.

To make the comparison more useful, I added support files to the GNSS_IMU repo so that anyone can run it themselves, either in its original form or after making any desired changes to the input parameters or code. If you are more interested in the final results than in the details of how it was done and how to do it yourself, you can skip to the end of the post.

Before jumping into the comparison, it’s probably worth a very brief discussion of what is the same and what is different between the two codebases. Both use an extended error-state Kalman filter with states for position, velocity, attitude, gyro biases, and accelerometer biases. KF-GINS also includes states for gyro and accelerometer scale factors; in GNSS_IMU, these states are optional. In GNSS_IMU, the position and velocity states are in the ECEF frame, while in KF-GINS they are in the NED frame. For attitude calculations, GNSS_IMU uses coordinate transformation matrices, while KF-GINS primarily uses quaternions. KF-GINS includes corrections for sculling and coning and uses trapezoidal integration in the INS mechanization; GNSS_IMU does not include these corrections. Also worth mentioning is that KF_GINS has a GPL3 license which has some restrictions on commercial use, whereas the GNSS_IMU uses the less restrictive BSD3 license.

Let’s start with some details about the dataset. Below is a plot from RTKPLOT of the RTK GNSS data from this dataset. Note that the data file does not include fix status, so to better illustrate the quality of the data, I created pseudo fix statuses based on the 3-axis position sigmas. I arbitrarily set points with sigmas less than 0.05 m to fixed, less than 0.25 m to float, and everything else to single. In the plot, green indicates fix, yellow indicates float, and red indicates single. The sigmas are shown with error bars, which are only large enough to be visible on the up/down axis. Note that there is a roughly 75-second GNSS outage starting at 14:24, with only a few low-quality positions in the gap.

The GNSS_IMU code uses RTKLIB position files for GNSS input and CSV files for IMU input with a slightly different format than the KF-GINS code, so I first created a Python script to convert the KF-GINS data files to GNSS_IMU-compatible files. This script is in the src folder and is named convert_KF_GINS_files.py. Before running it, I downloaded the dataset2 folder from the KF-GINS-Matlab repo and renamed it KF_GINS. This folder contains one file with GNSS position and velocity data, one file with IMU gyro and accelerometer data, and one file with ground-truth position, velocity, and attitude data. Running the script generates GNSS_IMU-compatible input files and saves them to the same folder.

Next, I created a config file for this data called ADIS16465_config.py. I used configuration parameters identical to those used by KF-GINS in order to ensure an apples-to-apples comparison. This was made a little more difficult because the KF-GINS code uses random-walk–based IMU parameters, while the GNSS_IMU code uses power spectral density (PSD)–based parameters. Rather than trying to directly convert between the two, I added an option to GNSS_IMU to allow the parameters to be specified in either format.

My original code only output positions and velocities in the GNSS frame, since I was using GNSS as my ground truth. Because the KF-GINS ground truth data is referenced to the IMU frame, I added a configuration option to output position and velocity in either the IMU frame or the GNSS frame. These will differ due to the lever arm between them. I also turned off the zero-velocity update, velocity matching, and yaw-alignment options in GNSS_IMU, since KF-GINS does not include these. KF-GINS uses the ground-truth position, velocity, and attitude to initialize the Kalman filter. GNSS_IMU uses the initial GNSS measurement for position and velocity initialization, but it does have the option to specify ground-truth attitude, which I did. KF-GINS always includes accelerometer and gyro scale factor states in the Kalman filter. GNSS_IMU also supports these states, but I left them disabled since I generally prefer to run without the additional complexity, and enabling them had a negligible effect on the results.

Next, I added a few lines to the header of GNSS_IMU.py, the main Python file, to point to the KF-GINS data folder, and then ran it. This generated a solution file containing positions, velocities, and attitudes.

The next step was to download the full KF-GINS Matlab library and run the solution on the dataset2 data. I don’t have a recent version of Matlab, but with a few minor modifications to the input/output code, I was able to run it using Octave, the open-source Matlab alternative. I found that it ran very slowly under Windows, but significantly faster when run under Linux using WSL (Windows Subsystem for Linux). Because of its additional complexity, it was still noticeably slower than GNSS_IMU, although it’s difficult to compare directly since one is running in Octave under Linux and the other in Python under Windows.

At this point, I had a GNSS_IMU solution file, a KF-GINS solution file, and a ground-truth file, all in slightly different formats. The next step was to generate a Python script to plot the errors of the two solutions relative to the ground truth. This script is named KF_GINS_compare.py and is also located in the src folder.

The initial results showed larger errors for GNSS_IMU than for KF-GINS, which was not entirely unexpected given that the code is relatively new and had not been subject to this level of scrutiny before. After digging through the results and the code, I did not find any major issues, but I did find several small errors that, when combined, were large enough to explain most of the discrepancies. Having both the ground truth and the KF-GINS results was extremely valuable during this debugging process. I tried to be careful to only fix actual errors, not to tune parameters or make changes that would artificially improve performance on this specific dataset, and I believe I was successful in doing so. All of these fixes have been added to the repo in recent commits.

OK, so finally we get to the actual data comparison. Below are the results from the first run after fixing the bugs in the code. In this case, I made every effort to ensure that the conditions were identical for both solutions. In each figure, the GNSS_IMU solution errors are shown on top, the KF-GINS solution errors are shown below, and for position and velocity, the bottom plot also shows the GNSS measurement errors. The orientation plots show roll, pitch, and yaw errors, while the position and velocity plots show north, east, and down errors. Overall, the GNSS_IMU errors are nearly identical to the KF_GINS errors. This is quite reassuring, given that the actual code in each code set is quite different.

Now, let’s do a little experimenting with some of the additional options in GNSS_IMU. One of the additional options is to allow float and single measurements to be weighted differently than from fix measurements since I have found in the past that the reported sigmas for fixed measurements are not only smaller, but also more reliable than those for float and single measurements. Specifically, the options allow different multipliers to be applied to the measurement sigmas for float and single measurements. Since KF_GINS does not have this option, I initially set both multipliers to 1. However, in my other sample configuration files, I typically use a float sigma multiplier of 2 and a single sigma multiplier of 5, and I consider these to be reasonable default values based on my limited experience with this code.

Below are the results after changing these two values. Note that in this case, as I mentioned previously, this data set does not include fix-statuses, so I am using the pseudo fix-statuses I derived from the position sigmas. The improvement in results is quite noticeable. The position and velocity peak errors are now less than half of the previous results. I was surprised that such a seemingly minor change could make such a large effect on the results.

Let’s try one more change. In my previous post, I described the vel_match option, which uses the first GNSS position after a GNSS outage to immediately adjust the position and velocity trajectory during the outage. This is obviously not useful for true real-time applications, but in my current application—plotting mapping data on an instrument display in near real time—it is useful to have this correction applied as soon as the outage ends. Setting the float and single sigma multipliers mentioned above to zero will ignore all float and single measurements and correct the outage between points with fixed status. Below are the results with this configuration. This change has only a very small effect on orientation, so I only show the position and velocity errors. Peak position errors are reduced from the previous result by roughly another factor of two.

OK, that’s probably enough—or maybe more than enough—plots for now. As usual with my posts, the goal isn’t to do a fully comprehensive comparison or exhaustive analysis, but rather to provide a relatively quick overview. For me, this exercise provided a useful validation of the GNSS_IMU code under conditions more challenging than it was originally designed for. Being able to compare results directly with both the ground-truth data and the KF-GINS solution was extremely helpful for debugging the code. It also demonstrates that a relatively lightweight GNSS/IMU implementation can perform surprisingly well when compared with a more complex codebase when not working with navigation-grade IMUs.

Loosely coupled sensor fusion with GNSS and IMU data

I’ve recently been doing some consulting work for Seescan, a company that provides equipment for locating underground utilities. The goal was to help them improve their precision positioning in challenging conditions. In particular, they tend to see rapid transitions from reasonably open sky conditions to very challenging conditions as the operator may follow a utility under a roof or overhang, or even inside a building for a short time. Their equipment already included a GNSS receiver with RTK or PPP-RTK corrections and an IMU, so the goal was to take advantage of the IMU to dead reckon through relatively short intervals where the GNSS signals were very poor or non-existent. As part of that project, I developed for them a tool written in python to help with the development and testing of their internal sensor-fusion code. Seescan has generously allowed me to share that code as an open-source project.

The result is the GNSS_IMU repository on Github. It is based on Matlab code that Paul Groves created as a supplement to his well-known text book, “Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems”. It is a loosely coupled solution, meaning that it uses the output of the GNSS solution as an input rather than the raw GNSS measurements. This makes it more compatible with their existing solution and easier to implement. In this post, I would like to give a quick introduction and demonstration of the code.

The original Matlab code from Paul Groves is a full simulator of both the GNSS solution and the IMU solution and includes the capability to generate simulated GNSS and IMU data from a set of position and orientation data. It includes both loosely and tightly coupled solutions. From this larger library, I have extracted just the loosely coupled solution but have added some features including lever-arm corrections, zero-velocity updates, initial yaw-alignment using GNSS heading, non-holonomic constraint updates, IMU scale factors, and simulated GNSS outages. Since I am primarily interested in post-processing and near-real time solutions rather than true real-time solutions, I have also added options for backwards solutions or combined forward/backward solutions, as well as first-fix backward smoothing after GNSS outages.

The code is based on a mechanization function for the IMU and an error-state Kalman filter for the sensor fusion with 15 default states and 6 optional states . The error-state Kalman filter tracks deviations from a nominal trajectory rather than the full states, because these errors remain small and approximately linear, especially for orientation. The 15 default states include three for orientation, three for position, three for velocity, three for accelerometer biases, and three for gyro biases. The six optional states include three for accelerometer scale factors and three for gyro scale factors.

Following the original MATLAB implementation, the code comments include references to relevant equations from the Groves textbook where possible.

I have also included two sample data sets. The first is from a vehicle driving on hilly residential streets followed by tight maneuvers in a parking lot. The second is from a handheld receiver while walking around in a backyard setting. Both data sets were collected with a u-blox X20P GNSS receiver in RTK mode and a TDK ICM-45686 IMU. The ICM-45686 is an inexpensive, roughly five-dollar consumer-grade IMU that still outperforms many of the lowest cost options in its class.

Like my python version of RTKLIB, this code is designed to get your hands dirty and dig directly into the code, so it does not have a slick GUI interface like RTKLIB. However, it does use RTKLIB format solution files for both the GNSS input solution and the combined GNSS/IMU output solution, so RTKPLOT can be used to explore and compare the input and output solutions and POS2KML can be used to easily generate KML files from input or output solution files to further explore with Google Earth.

The input files and configuration parameter file are specified in the header of the main python script, GNSS_IMU.py. The header is setup so that either sample data set can be run by commenting/uncommenting lines in the header. 

##########  Data selection and run configuration ########################################

# Uncomment these lines to run an example of walking with a handheld GNSS receiver
# import walk_config as cfg
# dataDir = r'..\data\walk_0827'
# fileIn = '1730_sf'

# Uncomment these lines to run an example of driving with a roof mounted GNSS receiver
import drive_config as cfg
dataDir = r'..\data\drive_0708'
fileIn = '1934_sf'

To run your own data, you would modify this code to point to your data files. I have included a couple of python scripts to help convert IMU and GNSS data collected with the RTKLIB stream-servers with time-tags enabled into the correct format but I will leave discussion of those for a future post.

Here’s a list of the configurable parameters for a given solution. This particular example is for the driving sample data.

# IMU parameters
imu_sample_rate = 100
imu_offset   = [0., 0., -0.65]  # offset from system origin: forward, right, down
gyro_noise_PSD =  0.0038  # deg/sec/sqrt(Hz)
accel_noise_PSD = 70    # ug/sqrt(Hz)
accel_bias_PSD = 7  # ug/sqrt(Hz)
gyro_bias_PSD = 3.8e-5   # deg/sec^2/sqrt(Hz)
accel_scale_noise_SD = 0.7
gyro_scale_noise_SD = 3.8e-6
imu_misalign = np.array([180.0, 0.0, 180.0])   # IMU orientation
imu_misalign += np.array([0, -6.79, 5.35]) # IMU misalign to body frame (degrees rpy)

# GNSS parameters 
gnss_offset = [0, -0.05, -0.65]  # offset from system origin: forward, right, down (m)

# Magnetometer parameters
mag_enable = False

# ratio of specs to process noise stdevs to account for unmodeled errors
imu_noise_factors = [1, 2, 4, 2, 0.25, 0.25]  # attitude, velocity, accel bias, gyro bias, accel scale, gyro scale
gnss_noise_factors = [1, 1]  # position, velocity

# Initial uncertainties
init =Init()
init.att_unc = [10, 10, 100]   # initial attitude uncertainty per axis in (deg)
init.vel_unc = [0.05, 0.05, 0.1]  # initial velocity uncertainty per axis (m/s)
init.pos_unc = [0.05, 0.05, 0.1]  # initial position uncertainty per axis (m)
init.bias_acc_unc = 0.2 # initial accel bias uncertainty (m/sec^2)
init.bias_gyro_unc = 0.2 # initial gyro bias uncertainty (deg/sec)
init.scale_acc_unc = 0.001
init.scale_gyro_unc = 0.001

# Run specific parameters
scale_factors = False  # Kalman filter scale factor states enabled
imu_t_off = -0.125  # time to add to IMU time stamps to acccount for logging delay (sec)
run_dir = [1]  # run direction ([-1]=backwards, [1]=forwards, [-1,1] = combined)
gnss_epoch_step = 1 # GNSS input sample decimation (1=use all GNSS samples)
out_step = 1  # Output decimation [1=include all samples]
float_err_gain = 2 # mulitplier on pos/vel stdevs if in float mode
single_err_gain = 5  # mulitplier on pos/vel stdevs if in single mode

# Velocity matching
vel_match = True  # do velocity matching at end of coast
vel_match_min_t = 1 # min GNSS outage to invoke vel match (seconds)

# Zero Velocity update
zupt_enable = True
zupt_epcoh_count = 50
zupt_accel_thresh = 0.25  # m/sec^2
zupt_gyro_thresh = 0.25  # deg/sec
zupt_vel_SD = 0.01     # standard dev (m/sec)
zaru_gyro_SD = 0.01   # standard dev (deg/sec)

# Initial yaw alignment
yaw_align = True # use GNSS heading to initialize yaw
yaw_align_min_vel = 1.0  # min vel (m/sec)
yaw_off = 0 # adjustment from initial yaw or mag (deg), just needs to be approximate
init_yaw_with_mag = False

# Non-holomonic constrains (NHC) update
nhc_enable = False
nhc_epcoh_count = 100
nhc_min_vel = 1.0
nhc_gyro_thesh = 20 # deg/sec
nhc_vel_SD = .25 # standard dev m/sec 
nhc_vel_SD_coast = 0.05 # standard dev m/sec

# Testing/Debug
disable_imu = False  # enable to Run GNSS only
start_coast = 40   # start of simulated GNSS outages (secs)
end_coast = 30 # end of simulated GNSS outages (secs before end)
coast_len = 15 # length of simulated  GNSS outages (secs)
num_epochs = 0  # num epochs to run (0=all)
gyro_bias_err = [0, 0, 0]  # Add constant error to gyro biases (deg/sec)
accel_bias_err = [0, 0, 0]  # Add constant error to acc biases (m/sec^2)
gyro_scale_factor = [1, 1, 1]
accel_scale_factor = [1, 1, 1]

# Plotting
plot_results = True
plot_bias_data = True    # plot accel and gyro bias states
plot_imu_data = False  # plot accel, gyro raw data
plot_unc_data = False  # plot Kalman filter uncertainties

In this post, I am not going to get in-depth into the different parameters, I am just going to demonstrate running the sample driving dataset under a few different configurations.

To start with, I will run GNSS_IMU.py without changing any parameters in the configuration file (drive_config.py). This is setup to run with 15 second long GNSS outages separated by 30 second intervals of valid GNSS data. It is what I call a near-real time solution, not a real-time solution, which means it uses the first fixed GNSS solution point after the outage to adjust the position/velocity trajectories during the outage. This is the “velocity matching” option in the configuration file. For applications, that don’t require a true real-time solution, this significantly reduces the errors during the GNSS outages. This mode provides initial position/velocity estimates in true real-time during the outage and then, immediately after the first post-outage fix, updated position/velocity for all of the outage epochs.

Running the script will generate an output solution file in RTKLIB format and a number of plots. The first two plots (below) indicate the position and velocity for the input solution without GNSS outages, and the loosely coupled sensor fusion output solution with GNSS outages. The blue/red/green traces are north/east/down for the input GNSS solution, the red/purple/brown are the sensor fusion solution with valid GNSS data and the pink/gray/yellow are the sensor fusion solution during the outages.

To estimate the solution errors during the GNSS outages, I use the input solution which doesn’t include the GNSS outages as a ground truth, and estimate the error to be the difference between the input solution without GNSS outages and the output solution with GNSS outages.

The fifth plot in the script output (below) shows these differences: The peaks in these plots indicate the maximum error during the 15 second GNSS outages. Remember that this is after the trajectory corrections done using the first fixed point after each outage, the errors would be significantly larger without these corrections.

The fourth plot in the set of output plots compares the solution orientation with the GNSS heading derived from velocity for all epochs with a velocity over a specified minimum velocity. For applications like a vehicle, which primarily moves in a single direction in the body frame, these should be similar provided the yaw is properly initialized. In this configuration, the “yaw_align” option is enabled which means that the first GNSS heading over a specified velocity is used to initialize the yaw state of the Kalman filter. Note that the velocity heading in this plot is derived from input velocity at the GNSS antenna which means it doesn’t compensate for any lever-arm between the GNSS antenna and the vehicle origin which is usually defined as the center of the rear axle, so these can differ by a small amount. In the case of a handheld receiver the two will tend to differ by much more since the orientation and direction of travel are more independent. R,P, and Y in the plot legend refer to roll, pitch, and yaw, and VEL_HDG is the heading derived from GNSS velocity.

The third plot shows the sensor fusion output in the body frame at the vehicle origin after compensating for the lever arm as defined in the configuration file. In this case, you can see it is almost entirely in the x-axis which is defined as the forward direction of the vehicle with small errors in the y-axis during the GNSS outages.

The last two plots show the estimated biases for the three accelerometer axes and the three gyro axes.

For comparison, let’s disable the velocity matching option which will give us the equivalent of a true real-time solution with IMU sensor fusion, but because there is no end-of-outage correction, the errors in position and velocity will be larger.

You can see in the plot on the left below that the position estimates during the GNSS outages still roughly follow the correct trajectory, but as you can see on the right, the peak position errors have increased during the 15 second outages from around 0.5 meters to over 10 meters. In theory this can be reduced for a vehicle if we constrain the sideways and vertical motion by enabling the non-holonomic constraint option in the configuration parameters. In reality, this feature is still a little experimental and I have not been able to achieve much benefit by turning this on. I need to look into this further and hope to make it the subject of a future post.

To compare this to performance without the IMU, I will turn set the “disable_imu” option in the configuration file which eliminates use of the IMU entirely in the solution. In this case, the trajectory will coast at constant velocity during the GNSS outages. This generates the plots below, showing position peak errors during the outages of over 200 meters. The velocity matching option is still disabled, so the result is true real-time.

As a final example, I will enable the velocity matching option with the IMU still disabled. This gives the plots below showing roughly 6-12 meter peak position errors during the 15 second GNSS outages. I’ve zoomed in on the position plot to more easily see the shape of the position errors with velocity matching.

I think that should be enough to get anyone started with the code if they are interested. If you try it out, I encourage you to experiment with different Kalman filter tuning values and with various options enabled or disabled, and to dig into the code itself. With a little effort, you should also be able to adapt your own GNSS/IMU data as input, though I’ll leave a detailed discussion of that for a future post.

So, to quickly summarize these result, we see that without the IMU, in this example, we were seeing roughly 200 meter peak position errors during the 15 second GNSS outage for real-time operation. By adding the velocity matching option, and dropping back to near-real time operation, we are able to reduce this to 6-12 meter peak errors. Adding a low-cost IMU reduced these numbers to roughly 10 meter peak errors for real-time operation and 0.5 meters for near-real time operation.

I should emphasize that sensor fusion is a relatively new topic for me and I’ve learned a lot by putting together this code set, but it is still a work in progress. I’m sure I’ve gotten a few things wrong or at least non-optimal so if you see anything I’m doing incorrectly or could do better, please leave a comment below, or better yet, generate a pull request for the Github code.

Quectel LC29HEA with improved antenna

In a recent post, I compared an L1/L2 u-blox F9P receiver paired with a u-blox ANN-MB antenna against a Quectel L1/L5 LC29HEA receiver paired with a very low cost L1/L5 antenna from Waveshare ($17 + shipping). I did this to get the combined cost of receiver and antenna down to an absolute minimum. The performance of this combination was reasonably good in challenging conditions but noticeably lower quality than simultaneous results from the u-blox receiver and antenna. At the time, I did not explore how much of this reduced performance was due to the lower cost receiver and how much was due to the lower cost antenna.

After publishing the post I was contacted by a representative from Quectel. He mentioned that they have done comparison testing between the LC29HEA and the u-blox F9P and have found the results to be more similar than my results if antennas of equal quality are used in the testing. He specifically recommended the Quectel YB0017AA L1/L5 antenna as being very similar to the u-blox ANN-MB antenna in size and performance. This antenna is available on AliExpress for $22 + shipping so is only marginally more expensive than the Waveshare antenna that I used. It’s also available from Digikey for a little higher price but faster shipping. I ordered the antenna from Digikey and after it arrived, repeated the previous comparison to the u-blox hardware using the new antenna.

Here’s a photo of the three antennas, the Waveshare L1/L5 on the left, the Quectel L1/L5 in the middle, and the u-blox ANN-MB on the right. As you can see, the Quectel antenna is more similar in size and shape to the u-blox antenna.

Waveshare, Quectel, u-blox antennas

To start, I did some signal strength (C/No) comparisons between the Quectel and Waveshare antennas using the C/N0 values reported by the receivers. I expected the larger Quectel antenna to show higher C/N0 measurements but I was unable to detect a consistently measurable difference between the two.

However, when I repeated the simultaneous RTK measurement collection with the u-blox hardware and the Quectel hardware in a moderately challenging environment in my backyard with a mixture of tree cover and obstruction from the house, the Quectel results this time were noticeably closer in performance to the u-blox baseline. More detailed testing would likely bring up differences between the two, but with my single, somewhat unrigorous test, the results were close enough that I was not able to discern that one performed better than the other.

Here are the results from the Quectel receiver and antenna :

LC29H->LC29H real-time RTK solution (green=fix, orange=float)

Here are the results from the u-blox receiver and antenna:

F9P->F9P real-time RTK solution (green=fix, orange=float, yellow=DGPS)

Fix rates were very similar between the two solutions at 84% for the u-blox and 88% for the Quectel. I would not consider this difference significant, especially since I suspect based on the previous comparison that the Quectel is a little more prone to false fixes.

Below is a plot of the distance between the two solutions when both were fixed. The two antennas were 34 cm apart, so the distance of each point from a circle of radius 34 cm indicates the combined error from both solutions. The number and magnitude of the errors are noticeably smaller than in a similar comparison for the previous experiment.

The large deviations in the vertical axis between the two solutions we saw in the previous comparison (shown on the left below) are not present in this comparison (shown on the right), again indicating fewer and smaller errors.

RTK solutions plotted in time (left=previous comparison, right=new comparison)

Overall, I consider both solutions quite qood considering the challenging conditions.

By the time you add the extra shipping for ordering the antenna separately from the receiver, the combined cost for receiver and antenna increases from $63 for the previous configuration to $81 for the new configuration, but given the improvement in performance, I suspect this is well worth it and is still significantly less expensive than the u-blox alternative.

Note that this and the previous comparison were made using the real-time internal RTK solutions from the receivers at a 5 Hz sample rate. For those more interested in PPK (post-processing) solutions using RTKLIB, this receiver may not be a good choice. So far, I have not been able to configure the receiver to output raw RTCM3 observations at greater than 1 Hz. 1 Hz is fine for static solutions but too slow in most cases for kinematic solutions.

Overall though, for anyone primarily interested in real-time RTK solutions for static or moving rovers, or post-processed solutions for static rovers, I would definitely suggest considering the Quectel LC29HEA module as a potential lower cost alternative to the u-blox F9P.

If you have worked with the LC29H and would like to share your thoughts and experiences, please leave a comment below.

Configuring the Quectel LC29HEA receiver for real-time RTK solutions

This is a follow up to my previous post in which I describe my first experience with the Quectel LC29H receiver. In this post I will go over the details on how I configured the rover and base receivers for a real-time RTK solution. These instructions are specifically for the LC29HEA receiver boards I bought on AliExpress, but, except for the references to the UART/USB slide switches, should work on any receiver board with an LC29HEA module.

The commands to configure the receiver are all described in the LC29H Protocol Specification, although it is not completely up to date with the latest firmware. Specifically, I had to read the latest firmware notes to find that 2 Hz and 5 Hz are now options for the position output rate. The commands all consist of strings of ASCII characters, but require a checksum at the end of the command. The easiest way to send and receive these commands is by using the command console window in the Windows QGNSS app from Quectel since it will generate the checksums for you with the click of a button.

To get started, download the app, open it up, connect the receiver to your computer with a USB cable, and connect an antenna to the receiver with the SMA connector. Check that the two slide switches on the receiver board are both set down towards the “USB” label.

Then use the “Set Device Information” option in the Device tab (or the “gear” icon in the toolbar) to find and configure the receiver communication settings. Set the Model to “LC29HEA” and the baudrate to 460800 as shown below.

QGNSS device configuration window

Once this is done, open up a “Text Data” window from the “View” menu and you should see a string of NMEA text messages scrolling through the window assuming your receiver is in its default configuration. It should look something like this:

NMEA messages in QGNSS text window

Next, open up a “Command console” window from the “Tools” menu which we will use to send configuration commands to the receiver.

Rover Configuration:

The next step is to enter the commands into the command console. We will first set up the rover. I’ve listed the command I used to do this below so that they can be cut and pasted into the QGNSS console. Note that the comments are for information only and should not be copied. I’ve also included the checksum for each command which can be copied or generated by using the “Checksum” button. Be aware that some commands don’t take effect until the parameters are saved to flash or the module is rebooted.

$PQTMRESTOREPAR*13 # restore PQTM params to default and reset
$PQTMCFGRCVRMODE,W,1*2A # set receiver to rover mode
$PQTMSAVEPAR*5A # save PQTM params to flash
# manually power cycle module
$PAIR062,2,0*3C # turn off GSA messages
$PAIR062,3,0*3D # turn off GSV messages
$PAIR062,5,0*3B # turn off VTG messages
$PQTMCFGNMEADP,W,3,6,3,2,3,2*37 # set decimal precision for NMEA
$PAIR050,200*21 # set pos output interval to 200 ms
$PQTMSAVEPAR*5A # save PQTM params to flash
# manually power cycle module

QGNSS Command Console Window ( commands to configure rover)

Once you’ve entered all the commands into the console, you can run them by consecutively pressing the “Send” button for each command. If you have made any modifications to any of the commands you will need to press the “Checksum” button first to recalculate all of the checksums.

These instructions assume the firmware on the module is no older than LC29HEANR11A03S_RSA which has an associated date of 10/31/23. You can use the command at the end of the screenshot above to check the version of firmware on your module. If your firmware is older than this, then I believe the only change you will need to make to these instructions is that the 5 Hz option is not supported and you will need to set the PAIR050 option to 100 instead of 200 which will set the output rate to 10 Hz instead of 5 Hz.

Base Configuration:

If you are going to use a second LC29H for a local base receiver, you will need to go through a similar process to configure that module. These are the commands I used:

$PQTMRESTOREPAR*13 # restore PQTM params to default and reset
$PQTMCFGRCVRMODE,W,2*29 # set receiver to base mode
$PQTMSAVEPAR*5A # save PQTM params to flash
# manually power cycle module
$PAIR432,1*22 # output RTCM3 MSM7 messages
$PAIR434,1*24 # output RTCM3 antenna position (1005)
$PAIR062,0,01*0F # Enable NMEA GGA message
$PQTMCFGSVIN,W,2,0,0,x,y,z*3B # set base location in XYZ coords
$PQTMSAVEPAR*5A # save PQTM params to flash
# manually power cycle module

Note that for the PQTMCFGSVIN command, you will need to use the XYZ coordinates of your base station and update the checksum before sending. Alternatively, the LC29H supports a survey-in capability to automatically generate an approximate base position.

Once you have sent this sequence of commands to the receiver, open up a “Binary data” window from the “View” menu and you should see the RTCM3 MSM observation messages and the 1005 base location message.

RTCM3 messages in QGNSS Binary data window

One minor issue I found is that the PAIR432 command to enable RTCM3 MSM7 messages does not seem to get saved to flash properly, so after cycling power, the module will switch to outputting MSM4 messages rather than MSM7. These contain slightly less information than the MSM7 messages but in most cases this should have negligible to no effect on the RTK solution. The MSM4 messages also require less bandwidth.

RTK Solution:

Once you are getting the correct outputs from both receivers, or a single receiver and external base correction stream, it is time to link them together to get an RTK solution. If you are using an LC29H as a local base station, you will want to connect it to an NTRIP caster to broadcast the correction stream. I describe one way of doing this with RTK2GO in this post. Once this is up and running, or you have access to NTRIP corrections from an external L1/L5 base station, the next step is to feed the base corrections to the rover using an NTRIP client. In my experiment I used an RTKLIB STRSVR stream server to do this, but for just checking out the receiver, it is simpler to use the NTRIP client tool built into QGNSS. Open this window from the “Tools” menu in QGNSS and configure it with your NTRIP address, username, password, port, and mountpoint. Once it is configured, set the “Connect to Host” button to on. If you are in an open sky view and have a ground plane under your antenna, you should see the “Quality indicator” field in the QGNSS output window switch to “Float RTK” and then to “Fixed RTK” as seen in the screenshot below. In this case, you can see from the RTK Float count field that it took 41 samples to converge from float to fix. Since the output rate is set to 5 Hz, this is about 8 seconds. The time to fix will vary depending on sky view, atmospheric conditions, and antenna quality.

So, that’s it, you should be up and running now and ready to collect some position data with centimeter-level accuracy. If you find any issues with these instructions or improvements to them, please leave a comment below.

Dual frequency RTK for less than $60 with the Quectel LC29HEA

The u-blox F9P dual-frequency RTK receiver has been my go-to choice for high precision RTK or PPK solutions ever since it first came out in 2018. Available initially on a receiver board with antenna for under $300, it offered a performance and price point far better than anything that preceded it. However, that was six years ago and although the price for the F9P receiver and antenna has come down a little to around $250 since then and they now offer a few more variations, it feels like u-blox hasn’t introduced anything dramatically new in a fairly long time.

Meanwhile, there have been a number of lower cost dual-frequency receivers introduced recently from other companies which look quite promising. I hope to take a look at a few of these eventually, but have decided to start with the LC29H from Quectel.

The F9P is an L1/L2 receiver, while the LC29H, like most of the newer receivers is an L1/L5 receiver. The L5 signals have some advantages over the L2 signals, but (like the L2C signals used by the F9P and other low-cost L1/L2 receivers) are not yet available on all satellites.

The LC29H module comes in several variants. The newest ones most generally available on receiver boards are the DA, BS, and EA variants. The specs for these are shown below. The EA variant is the one most recently available and is the only one that will output raw observations or RTK solution points at greater than 1 Hz. A 1 Hz sample rate is fine for static applications, but for dynamic rovers a higher sample rate is usually necessary. For this reason, I will focus on the LC29HEA.

Quectel LC29H specs for DA, EA, and BS variants

Probably because it’s relatively new, it is a little more difficult to find receiver boards built with the EA module. I ended up ordering two receivers from this link on AliExpress. If the link no longer works by the time you read this post, you can probably find something similar using their search option. The price for this one was $57.69 for the receiver and antenna plus $4.79 shipping to the U.S. In my case, the boards arrived in less than two weeks with no issues. I ordered the boards without the antennas since I missed the combined option when I made the order. For my initial evaluation I used what I’m guessing is a similar low-cost L1/L5 antenna from Waveshare for $17. Waveshare also sells LC29H receivers but at the current time, only with the AA, BS, and DA variants.

Quectel LC29HEA receiver with L1/L5 antenna from AliExpress

To start with, I chose to compare a real-time internal RTK solution using a pair of F9P receivers for base and rover to a real-time internal RTK solution using a pair of LC29HEA receivers. If I was trying to compare the two receivers directly, I would use the same antenna for both rovers. However, for this experiment, I wanted to compare a lower cost solution that included both a lower cost receiver and a lower cost antenna. This means that the results can not be used to draw any conclusions regarding differences between the receivers directly, only differences between the combinations of receiver and antenna. I hope to do another comparison later using a single antenna to more directly compare the receivers.

I connected the two base receivers through an antenna splitter to a Harxon geodetic GPS-500 antenna on my roof. Although this antenna is advertised as L1/L2 only, I have found it gives signal strengths in the L5 band equivalent to L1/L2 so felt it was fair to use in this comparison. I connected the F9P rover to a u-blox ANN-MB L1/L2 antenna and the LC29H rover to the low-cost L1/L5 antenna from Waveshare.

I have done most of my previous comparison tests with antennas on top of vehicles. For this test I chose something a little more challenging, and walked around my backyard, sometimes in the relative open, sometimes close to large trees, and sometimes close to the house, but never underneath dense tree foilage.

I configured the base receivers to output RTCM3 MSM7 observation messages at 1 Hz and broadcast them to the internet using the RTKLIB STRSVR stream server app and a couple of RTK2go NTRIP servers similar to what I describe in this post.

The LC29H can only be configured to output RTK positions at 1 Hz or 10 Hz, so I chose 10 Hz for the LC29H rover. [Note 5/1/24: the release notes for the latest firmware indicate that options for 2 Hz and 5 Hz have just been added] According to the datasheet, the F9P has a maximum RTK position output of 7 Hz when running with all four GNSS constellations but for this experiment I left it at 5 Hz which I find an adequate sample rate for most applications.

Using a laptop running four instances of STRSVR, I configured two of these to stream the NTRIP corrections coming from the two base receivers to go to the two rover receivers over USB cables. I used the other two to stream the real-time output of the two rovers to a couple of files. I describe a similar setup for just the F9P in this post. I used the u-blox u-center app to configure the F9P, and the Quectel QGNSS app to configure the LC29H. The LC29H configuration commands are somewhat cryptic and need to be typed into a command console in the Quectel app, so this was not as easy as configuring the F9P, but the details for doing this are reasonably well documented in the Quectel LC29H Protocol Spec.

I will share the details of exactly how I configured the LC29H base and rover modules in a separate post. The biggest issue I found here was that when running at 10 Hz, not all features can be enabled at the same time and the way the module handles this can sometimes be confusing. For example, you can not output RTCM observation messages and 10 Hz RTK real-time positions at the same time. This setup is not necessary for this experiment but can be useful if you want real-time position information but also want to postprocess the raw observations later with RTKLIB. Also, some commands would take effect immediately and others only after config parameters were saved to flash, or when the module was reset.

After all four receivers were configured, I walked around the backyard with the laptop connected to the two rovers, and the two rover antennas mounted together about 30 cm apart. When close to the house or the trees I was careful to keep the spacing between both antennas and the obstruction similar to each other.

RTKPLOT will plot two real-time NMEA position streams, so I also had an instance of this running on my laptop to plot the two solutions real-time. The screenshot below is from my laptop while I was running the experiment, showing the four stream server windows and the RTKPLOT window. The yellow/green plot lines are float/fix status for the F9P and the olive/blue plot lines are float/fix status for the LC29H. I took this screenshot before I started walking around. The antennas were static and I was periodically covering both antennas to force loss-of-lock and reconvergence of the two solutions. Both receivers recovered a fix status in a reasonable time but the F9P was generally two to three times faster.

Screenshot of my laptop while collecting data

The screenshot below shows the real-time RTK position output for the F9P collected while walking around. As you can see, most of the solution is fixed, but close to the house and close to the tree line on the right where the trees are tallest and thickest, there is more float solution. The points closest to the house are actually under the eaves of the roof. The float points in the lower middle of the image are from when I was covering and uncovering the antenna as I described above. I generated this image using the POS2KML app in RTKLIB to generate a KML file and then displayed it with Google Earth.

F9P->F9P real-time RTK solution

Here is the same image for the LC29H real-time RTK solution. Again, most of the open-sky area is fixed but a higher percent of the more challenging regions are float.

LC29H->LC29H real-time RTK solution

I don’t have a ground truth for this data so I can’t directly compare the accuracy of the two solutions. What I can do is plot the difference between the two solutions to get a general idea of the solution errors. Since the antennas are a fixed distance apart, the difference between the two solutions should be a circle with radius equal to the antenna separation distance, assuming the antennas were kept fairly close to level while the data was taken. Here’s the result of plotting the solution difference with RTKPLOT

A large percentage of the points are on the circle, in this case with radius 34 cm, but there are also a fair number of points not on the circle. Given the challenging nature of this experiment, this is somewhat expected. Float points (yellow) off of the circle are less concerning since we know the accuracy is lower for them. Fix points (green) off the circle however are more concerning and generally indicate false fixes. Plotting both solutions versus time and comparing the vertical components can help tell us where the false fixes are and which solution they occurred in. This is because the trajectory covers the same ground multiple times so we know what the valid range of altitude is. The region below circled in red indicates both solutions are fixed but differ by well over the expected accuracy of a fixed solution. Since the blue (LC29H) solution is outside the range of altitudes from the rest of the data we can determine that this is the incorrect solution in this case.

F9P solution (green/yellow) and LC29H solution (olive/blue)

Based on the results of this single experiment, it appears that the Quectel LC29H with a very low cost antenna performed reasonably well in challenging conditions and deserves further investigation. It did not perform as well as the u-blox F9P with a more expensive antenna in these conditions, however given that it cost roughly one quarter as much, I think it did quite well. In less challenging conditions such as onboard a drone where sky views are generally less obstructed, I would expect the differences to be smaller.

Performance and cost are important, but of course they are not the only factors to consider when selecting a receiver. I did find that overall the documentation is more complete for the F9P, it is both easier to configure and more configurable than the LC29H, and there are some features such as event logging that are available on the F9P but not the LC29H.

One last thing to consider is that the overall performance results are a combination of, among other things, the antenna, the front end of the receiver and the internal RTK solution. I hope to do further experiments using a common L1/L2/L5 antenna and post-processing with RTKLIB to measure the differences in each of these components separately.

Well, I think that is enough for now. I hope to follow up with some more posts as I gain a little more experience with this receiver and eventually take a look at some other very low cost L1/L5 receivers as well.

If you have worked with the LC29H and would like to share your thoughts and experiences, please leave a comment below.

2023 Google Smartphone Decimeter Challenge

Google has just kicked off the 2023 Smartphone Decimeter Challenge on Kaggle. Similar to the previous two years, it is a competition to see who can generate the most accurate solutions for a large number of raw observation data sets collected with Android phones on vehicles driven around the Bay Area and Los Angeles. In the previous two competitions, all the phones were higher-end models that supported dual-frequency GNSS. This year, they are also including mid-range phones with only single frequency GNSS to better represent the total population of phones. This will make the solutions more challenging.

They have also increased the total prize money from $10,000 to $15,000. The competition started just over a week ago and will run until May 23 next year. The top three winners will be invited to present papers on their solutions at the 2024 ION GNSS+ conference.

Last year, after the competition had concluded, I shared a version of my final solution in a notebook on the competition’s Kaggle page. When copied to a local computer and run, it will generate a result that can be submitted to Kaggle and will place 5th in last year’s public leaderboard.

I have just published an updated version of this notebook on this year’s Kaggle page. It is the same solution as the previous version, just updated to run with this year’s data and also modified so it will run more easily on Linux as well as Windows. For anyone interested in joining the competition, this latest version will produce a score of 1.80 meters, which at the moment is good enough for first place. I’m sharing it to help promote a stronger competition as well as to encourage the use of RTKLIB.

Below is a screenshot of the state of the public leaderboard as of September 21. You can see the most recent version of the leaderboard here.

For reference, last year’s winning score on the public leaderboard was 1.38 meters, and I expect this year’s winning score to be lower. So, this is only a starting point, but it should give anyone interested in competing an opportunity to take advantage of previous work and jump in near the front of the pack (at least for the moment).

I’m happy to answer any questions about using RTKLIB in the competition but the rules require that I do that in the competition’s Kaggle discussion group so that the information is available to all participants. There was quite a bit of collaboration between competitors in last year’s competition as well as a lot of information shared after the competition, all on the competition’s Kaggle discussion group page, so check that out if you haven’t already.

More details of the optimizations I have made to RTKLIB for smart phone observations are described in these links:

Kalman Filter Adjustments in RTKLIB (the Options-Statistics tab)

In this post, I am going to describe how the Kalman filter in RTKLIB is configured and how it can be adjusted. In most cases, the default values work fine and there is no need to adjust these parameters, but in some cases it can be useful. Also, understanding these parameters can provide useful insight into how RTKLIB generates its solutions.

If you are interested primarily in how to set these values and not in what they do, then you may want to jump down to the end of the post where I discuss a few examples of where I have found adjusting these values to be helpful.

Parameters to enable and disable various filter states are sprinkled over the different tabs in the options menus of RTKPOST and RTKNAVI but the parameters that determine the inner details of the Kalman filter calculations are all listed on the Statistics tab as shown below.

Statistics tab in RTKPOST

Note that the parameters are divided into two sections, measurement errors and process noises. In a Kalman filter, each measurement will have a measurement error estimate associated with it, and each filter state will have a process noise estimate associated with it, both defined in terms of standard deviations (or sigmas). The Kalman treats all of the measurement errors as normally distributed with zero mean. In general, increasing the error estimates of the measurements relative to the process noises will cause the filter to rely more on the model, and reducing the error estimates of the measurements will cause the filter to rely less on the model and more directly on the measurements.

The description below will apply specifically to differential solutions (RTK/PPK) but the PPP solutions use a similar method.

MEASUREMENT ERRORS:
The primary measurement inputs to the Kalman filter are code and phase observations. The amount of error in each observation will be determined by many factors which we cannot determine directly, but in general, both satellite elevation and signal strength will correlate with measurement error. Since the receiver reports a signal strength (C/N0) for each measurement and we can determine the satellite elevation from the ephemeris messages, we can use either or both of these to generate an estimate of the error. For differential solutions (RTK/PPK), errors will also increase as the baseline between rover and base increases. Another option that is sometimes available is to use an error estimate provided directly by the receiver.

The original (2.4.3) version of RTKLIB uses the sum of four squared terms to determine the error estimate as a variance for each phase observation. The first term is a constant, the second term is a constant divided by the sine of the satellite elevation, the third term is a constant multiplied by the distance between base and rover, and the fourth term is a constant multiplied by the time since the last observation. The default values for the first two constants are three mm, and the third is zero. The last term is included to account for satellite clock stability and defaults to 5e-12 sec/sec. These are listed in the options tab as “Carrier-Phase: a+b/sinEl”, “Carrier-Phase: Baseline”, and “Satellite Clock Stability”.

The demo5 version of RTKLIB defaults to the same algorithm and values as the 2.4.3 code but has three additional terms which all default to zero. The first two terms allow the error estimate to be adjusted by the signal strengths of the base and rover observations using the formula a*10^(.1*MAX(snr_max – snr), 0). The two additional terms are a and snr_max in the above equation and correspond to Carrier-Phase: SNR / SNR maxDb” in the options tab. This squared result is added to the previous terms. Using the convention of RTKLIB, snr and snr_max are actually C/N0 values as reported by the receiver.

The last term in the demo5 code is used to include the reported error estimate from the receiver and is only available when the observation files were generated from u-blox binary files with the “-RCVSTDS” receiver option enabled. This option causes RTKCONV or CONVBIN to add the receiver error estimates to unused fields in the RINEX observation files. The receiver error estimates from the RINEX files are multiplied by this last term, squared, and added to the previous terms. This term is labelled “Carrier-Phase: Rcv Errs” in the options tab.

Thus, for the demo5 code, by setting the appropriate terms to zero and non-zero values, the phase observation error estimates can be derived from either the elevation, the signal quality, the receiver error estimates, or any combination of the three.

Finally, for both the 2.4.3 and the demo5 code, the error estimate is adjusted based on which constellation the observation is from. These adjustments are defined in the rtklib.h source file and can not be modified with the config file. GLONASS and IRNSS adjustments are set to 1.5, SBAS is 3.0, and the other constellations are set to 1.0.

Everything so far describes how the error estimates are derived for the phase observations. To determine the error estimates for the code observations, the phase error estimates are simply multiplied by a constant determined by which frequency band the observation is from. These are the “Code / Carrier-Phase Error Ratio L1/L2/L5” parameters in the options tab above. These ratios all default to 100 in the 2.4.3 code which is probably a more appropriate value for higher-end receivers, and 300 in the demo5 code which seems to work better for low-cost receivers.

There is also a parameter for setting the measurement noise for the doppler measurements but these are used only very peripherally in the solution, so adjusting this parameter will have little to no effect on your results.

PROCESS NOISES:
The primary states in the Kalman filter are the single-difference phase biases and the positions. There are nine position states if the “Rec Dynamics” option is enabled (x,y,z for position, velocity, and acceleration) or just three if it’s not enabled (x,y,z for position). For long baseline solutions, additional states can be enabled for tropospheric and ionospheric delays.

Each of these states has a process noise associated with it. The process noises are estimates of the amount of unmodeled error in the state, expressed as a standard deviation. In the case of the position states, the model does not account for any external accelerations of the receiver. Therefore we need to specify non-zero process noises for the acceleration states, assuming the receiver is not static. The position and velocity changes that result from these unmodeled accelerations are accounted for by the model, and hence the position and velocity process noises are set to zero. The process noises for the acceleration states are set with the “Receiver Accel Horiz/Vertical” parameters in the options tab and have units of m/sec^2 The smaller these parameters are set, the less the position solution can jump around from noise, since the filter will constrain the motion to smaller accelerations, but the more lag will be introduced in the filter response if there are real accelerations larger than the specified process noise. If the dynamics option is not enabled, meaning there are no velocity or acceleration states, then there will be unmodeled errors in the position states. In this case, the code does not use the acceleration state process noise values , but sets the process noises for the position states to a large but unconfigurable value.

The phase bias states also have process noises assigned to them and these are set with the “Carrier-Phase Bias” parameter in the options tab. They are set in units of carrier phase cycles but in this case they don’t translate as well into unmodeled physical errors as the acceleration errors do. By definition, these states are modeling values that can only be integers (or half-integers in the case of a phase lock loop that is not fully locked) and don’t change unless there is a cycle slip. The variance of the state is reset if a cycle slip is detected or reported, so the only unmodeled errors should be undetected cycle slips. These are large, discrete, and very infrequent, about as un-gaussian as you can get, yet we are forced by the math behind the Kalman filter to treat them as normally distributed. The default value is 0.0001 cycles which is probably based more on what gives the best results rather than anything physical. Increasing this value will cause the filter to weight recent measurements more heavily relative to earlier measurements.

There are also process noise values for the tropospheric and ionospheric delay states if these states are enabled for long baseline solutions.

Choosing Config Parameter Values:
As I mentioned before, most of the time I just use the defaults for these parameters and don’t try to tune them. There are a few exceptions however, so I will describe them here.

First, I have found that setting the code/carrier-phase error ratio to lower values for more expensive receivers (~100) and to higher values for low-cost receivers (~300) tends to give better results. Lower values increase the weight of the pseudorange measurements relative to the carrier phase measurements. I suspect that the more expensive receivers are able to improve the quality of the pseudorange measurements more than they are able to improve the quality of the carrier-phase measurements relative to the low-cost receivers but that is only conjecture.

Sometimes, I find that forward-only solutions give higher fix percentages than combined (forward+backward) solutions. The reason for this is that RTKLIB validates every fixed point in the combined solution and if the forward solution differs from the backward solution point by more than four standard deviations, then the point is downgraded to float. The purpose of this step is to identify and remove false fixes, but if the observation measurement error estimates are too small, then even quite small differences between the forward and backward solutions that were not caused by false fixes can be downgraded. In this case increasing the constant and elevation weighting terms of the phase observation error estimates from 0.003,0.003 to 0.005,0.005 or 0.006,0.006 usually eliminates some of the unnecessary downgrades and improves the fix rate.

For cell phone solution I increase these values by even more to account for the very poor signal quality.

In a more challenging receiver environment with many obstructions and reflections, observation quality may correlate less strongly to elevation and more strongly to signal quality than in more typical situations, especially for short baselines where atmospheric errors are less significant. Especially with cell phones, several researchers have reported better results when weighting observations by signal quality rather than elevation.

I also suspect that incorporating the receiver quality metrics in addition to the other factors should be helpful but have not proven this in practice, despite some attempt to do this.

I have also noticed that as the interval between base observations increases, the effect of the satellite clock stability term becomes more and more dominant on the observation measurement errors to the point where all of the other terms have almost no effect. I suspect that this is not realistic or desirable but I have not investigated it closely.

If you have found other examples where adjusting these values improved your results, and are willing to share your experience, please leave a comment below.


Raspberry Pi based PPK and RTK solutions with RTKLIB

It’s been over six years now since I published my last post on how to run RTKLIB on a Raspberry PI, so it’s more than time for an update. In my previous post, I described using a Pi Zero as a data logger for a u-blox M8N for PPK solutions. In this post I will work with a Pi Model 4 and a u-blox M8T to demonstrate both logging for PPK solutions and a real-time RTK solution. The good news is that this time no soldering is required since we are going to use the USB port on the Pi to connect the receiver. These instructions will work with any u-blox receiver that supports raw observations and any model Pi that has USB ports for peripherals. With minor modifications, they can be used with any receiver that has a USB or UART port and supports raw observations.

Here’s an image of the assembled setup. The Pi in the center, and the u-blox M8T receiver is on top. We will use a wireless connection to talk to the Pi from an external computer so there is no need for a keyboard or display.

Raspberry Pi with u-blox M8T receiver


Step 1: Configure the Pi

The first step is to configure the Pi in “headless” mode so that we can talk to it from an external computer. This is quite straightforward and well-explained in this post, so I will not describe how to do it here. Only steps 1 and 2 in the post are required for this exercise. If you plan to use this for RTK solutions, be aware that the Pi will rely on the wireless connection to the internet for base station observations. This means that if you don’t want to be limited to using it within range of your home wireless router, then you will probably want to connect to a hot spot from a cell phone. If you are just interested in collecting data for PPK solutions, then it doesn’t matter.

After you have completed steps 1 and 2 above, you should have a Putty window open and have logged into your Pi. The next step is to build and install the RTKLIB code. The commands below will clone the RTKLIB code from the Github repository, compile the stream server app (str2str) and the RTK solution app (rtkrcv), and copy the executables to a location where they can be accessed from any directory.

> sudo apt update
> sudo apt install git
> mkdir rtklib
> cd rtklib
> git clone https://github.com/rtklibexplorer/RTKLIB.git
> cd RTKLIB/app/consapp/str2str/gcc
> make
> sudo cp str2str /usr/local/bin/str2str
> cd ../../rtkrcv/gcc
> make
> sudo cp rtkrcv /usr/local/bin/rtkrcv
> cd ../../../../..


Step 2: Configure the receiver

Before we connect the u-blox receiver to the Pi, we will need to configure it to output raw observation and navigation messages. The easiest way to do this is from your computer with the u-center app which can be downloaded from the u-blox website. Connect the receiver to your computer with a USB cable, start u-center, and connect to the receiver using the “Connection” option in the “Receiver” tab as shown in the image below.

u-center: Connect to receiver

Next, use the “Messages View” window from the “View” menu to enable the RAWX and SFRBX messages as seen below. While you are in the messages view, you can also disable any unnecessary NMEA messages to save communication bandwidth.

u-center: Enable raw observation and navigation messages


Next, we will switch to the “Configuration View” window to configure any other desired settings and then save them to flash . I would recommend verifying that all constellations are enabled with the “GNSS” command and that the sample rate is set to the desired value with the “RATE” command. I usually set this to 5 Hz. I would also recommend disabling both UART ports with the “PORT” command if you are not using them. If the baud rates are set too low, they will limit bandwidth on all ports including the USB port, even if nothing is connected to those ports. Finally, use the “CFG” command to save the settings to flash as shown below.

u-center: Save settings to flash

Step 3: Verify the data stream(s)

Next, we will confirm that we are receiving data from the rover receiver and if running a real-time solution, also from the base receiver. This step is not absolutely essential, but it does verify that we have the individual pieces working before we put it all together, and also gives some practice using the RTKLIB str2str command.

Disconnect the rover receiver from the computer and connect it to the Pi using a USB cable as shown in the image at the top of this post. Enter the following commands into the Putty console to create a new folder and run the stream server. This will connect to the USB port on the Pi. If you are using a UART port, you will need to use the appropriate port name.

> mkdir data
> cd data
> str2str -in serial://ttyACM0

The output of the receiver should now scroll across the Putty console screen. If you have any NMEA messages enabled, you should be able to see them mixed in with a bunch of random characters from the binary messages. Once you’ve confirmed the data stream, hit Control C to stop it.

If we want to log the receiver output for a PPK solution, we just need to add a file name to the previous command to redirect the data stream from the screen to a file. The command below will do this, using keywords in the file name to create a name that includes the current month, day, hour, and minute.

> str2str -in serial://ttyACM0 -out rover_%m%d_%h%M.ubx

The image below shows the expected output of both commands.

Verification of receiver data stream

If you are using the Pi just to log receiver data then you are done at this point unless you want to configure the Pi to make it automatically start collecting data whenever it is turned on. There are several ways to do this, all described in this post. Modifying the rc.local file is the simplest method.

For those who would prefer to run an RTK solution rather than just log data for a PPK solution, the next step is to confirm the base data stream. We will use the “str2str” command again, but this time we will specify the input to be an NTRIP stream using the format:

“ntrip://username:password@ipaddress:port/mountpoint”

In my case, the command looks like this: (with the username and password removed)

> str2str -in ntrip://username:password@rtgpsout.unavco.org:2101/P041_RTCM3: -out temp.log

If everything is working properly, you should see non-zero transfer numbers and no errors, as in example above, in which case you can use Control C again to stop.

Note that if your NTRIP provider is using a VRS (Virtual Reference Station), then things are a little more complicated. We will need to send our local position inside of a GGA message. For this to work, you must have enabled the NEMA GGA message when configuring the receiver. To route these GGA messages back to the NTRIP server we will need to connect the stream server output to the receiver and enable the relay back feature with the “-b” option. Here’s an example I used to connect to test this with a VRS NTRIP server.

str2str -in ntrip://username:password@na.l1l2.skylark.swiftnav.com:2101/CRS -b 1 -out serial://ttyACM0

Step 4: Run the RTK solution

OK, now that we’ve confirmed that we are getting data from base and rover, it’s time to generate an RTK solution. We will use the “rtkrcv” console app in RTKLIB to do this, which we installed in Step 1.

We will need a configuration file for rtkrcv. You can use the “rtknavi_example.conf” file included with the demo5 release as a starting point but you will need to edit the stream configuration settings. Below are the settings I changed as well as a few important ones worth verifying are correct for your configuration. I have it configured to write the output to a file in LLH format. If you want the output in NMEA messages you can either change output stream 1 to “nmea” format or enable output stream 2 to get both a file and a stream of NMEA messages.

pos1-posmode =kinematic  # (0:single,1:dgps,2:kin,3:static)
pos1-frequency =l1  # (1:l1,2:l1+l2,3:l1+l2+l5)
pos1-navsys =13  # (1:gps+2:sbas+4:glo+8:gal+16:qzs+32:comp)
pos2-armode =fix-and-hold # (0:off,1:cont,2:inst,3:fix-and-hold)
pos2-gloarmode =fix-and-hold # (0:off,1:on,2:autocal,3:fix-and-hold)
out-solformat =llh #    (0:llh,1:xyz,2:enu,3:nmea)
ant2-postype =rtcm # (0:llh,1:xyz,2:sing,3:file,4:rinex,5:rtcm)
inpstr1-type =serial (0:off,1:ser,2:file,3:,...,7:ntrip)
inpstr2-type =ntripcli # (0:off,1:ser,2:file,3:,...,7:ntrip)
inpstr1-path =ttyACM0
inpstr2-path =usrname:pwd@rtgpsout.unavco.org:2101/P041_RTCM3
inpstr1-format =ubx # (0:rtcm2,1:rtcm3, ...)
inpstr2-format =rtcm3 # (0:rtcm2,1:rtcm3,...)
inpstr2-nmeareq =single # (0:off,1:latlon,2:single)
outstr1-type =file # (0:off,1:serial,2:file, ...)
outstr2-type =off # (0:off,1:serial,2:file, ...)
outstr1-path =rtkrcv_%m%d_%h%M.pos
outstr2-path =
outstr1-format =llh # (0:llh,1:xyz,2:enu,3:nmea)
outstr2-format =nmea # (0:llh,1:xyz,2:enu,3:nmea)
logstr1-type =file # (0:off,1:serial,2:file, ...)
logstr2-type =file # (0:off,1:serial,2:file, ...)
logstr1-path =rover_%m%d_%h%M.ubx
logstr2-path =base_%m%d_%h%M.rtcm3

I like to use WinSCP for editing and transferring files between the Pi and external computer but there are many other ways to do this. When you are done, the edited configuration file needs to be in the current folder you will run rtkrcv from. For my example, I renamed it “rtkrcv_pi.conf”

To run rtkrcv with a configuration file named “rtkrcv_pi.conf”, use the following commands:

> rtkrcv -s -o rtkrcv_pi.conf
  >> status  1

If all is well, you should see a status screen updated every second that looks something like this:



I changed the Putty display defaults to make this a little easier to read. I’ve also highlighted in yellow some of the numbers to check to make sure they look OK. Make sure you are seeing base RTCM location messages (usually 1005). If you want to check the input streams in more detail, you can use control c to exit the status menu, then enter “?” to see some of the other rtkrcv commands. To exit rtkrcv, use the “shutdown” command.

If all of your inputs look good, your solution is not working, and it is not obvious why, you can rerun rtkrcv with a “-t 3” in the command line. This will enable trace mode which will create a trace debug file which may offer clues as to what is wrong.

This should be enough to get you started. To explore more configuration options, see the str2str and rtkrcv sections in Appendix A of the RTKLIB users manual.


Google Smartphone Decimeter Challenge 2022

In a previous post, I described my experience using RTKLIB to analyze smartphone GNSS data from last year’s Google Smartphone Decimeter Challenge. In that case, I did not get involved until after the competition was complete. After making a few modifications to RTKLIB to handle the relatively low quality smartphone data, I was able to generate a set of solutions that would have placed 5th out of 811 teams in the final standings. I shared the code to duplicate my results in this code release on Github. It includes a custom version of RTKLIB with changes specifically made for the smartphone data, as well as a set of python scripts to automatically run solutions on all of the 2021 Google test rides.

Google is is hosting a second competition this year. It started in the beginning of May and will finish at the end of July. This year I decided to join the fun and submit some results while the competition was still ongoing.

Since last year, I had already incorporated all of the changes that were previously in only the GSDC version of RTKLIB, into the main branch of the demo5 fork of RTKLIB. These are in the latest b34f release, so the special release is no longer required.

Google changed the format of some of the files for this year’s competition and so I did have to rewrite the python scripts. One of the more significant changes they made was to include only one set of phone data for each ride in the test data set. Last year it was possible to combine results from multiple phones on a single ride to improve the results but that is not an option this time.

In order to encourage participation in this year’s competition, I have shared the code and instructions to duplicate my initial attempt on this year’s data in a Kaggle notebook . If followed correctly it will generate a score of 3.135 meters when submitted to Kaggle, the competition host. At the time I first published it, it was good enough for first place. However, the competition has picked up since then, and some teams have taken advantage of this code. It will no longer get you into first place, but it will still put you into a tie for 21st place out of 234 teams. This means that anyone interested in jumping in now can still start near the front of the pack.

Since sharing the notebook, I have made a few local tweaks to the code, config files, and python scripts which improve my score to 2.152 meters. This is currently good enough for first place, but given that there are nearly two months left in the competition, I don’t think this will be good enough to win without further improvements.

To keep things interesting, I don’t plan to share my most recent changes until the competition is complete but anyone who follows some of the suggested hints at the end of my Kaggle notebook should be able to get a good part of the way there. To get all the way there will require a little more ingenuity but I also believe there is still plenty of room for further improvement on my results.

However, I suspect that winning the competition using RTKLIB will require more than just configuration changes and python script changes. I believe it will also require making changes to the RTKLIB code itself.

As anyone who has worked with the RTKLIB code is probably aware, it can be quite a challenging environment to work in. To make things easier and to encourage innovation to the code and algorithms I have recently ported a subset of RTKLIB sufficient to generate PPK solutions into Python which I described in this post. The actual code is available on Github here. I have also generated a second Kaggle notebook with instructions on duplicating the C/C++ version results on the Google data with the Python code. I have not actually submitted the results of this code to Kaggle, but based on results from this year’s training data set, and last year’s test data set, I believe this code should give slightly better results than the C/C++ code.

The python code is primarily intended for those planning to develop or modify algorithms internal to the PPK solutions and not just running the code as-is or with just configuration changes. For those users, the C code will run much faster. However, the python version provides a friendlier development platform. When development is complete, the modified python code can either be run on the complete data set on a faster PC with a little patience, or the completed changes can be fairly easily ported back into the C code since two code sets are very closely aligned. This alignment includes file names, function names, variable names, and comments. The code does not align on a line by line basis because of extensive use of Numpy in the python code, but structurally it is very similar.

Based on the discussion threads on the Kaggle forum for this competition, it appears that most competitors are more familiar with machine learning and post-solution filtering techniques than they are with GNSS theory. I suspect anyone who already has a reasonably solid background in GNSS can do quite well in the competition without an enormous amount of effort. Using some of the tools I describe here should help to get there even more quickly.

My hope is that providing these tools will encourage at least a few more people from the GNSS community to participate and help them to do well. For any of you who decide to take the challenge, I wish you good luck and hope to see you near the top of the leaderboard!

A Python version of RTKLIB

In this post I would like to introduce a new project that I have recently been working on.

RTKLIB can be a great tool for exploring the world of precision GNSS solutions. There are many configuration options (some say too many!) which allow the user to investigate the effects of different algorithm settings. It is also open source, so the user can see exactly what the code is doing behind the scene and even modify the code to explore options outside the available configuration set.

However, the code is quite complex and written in a dense style that can be quite intimidating to the casual explorer who wants to dig deeper. The compiled nature of C/C++ also makes the development platform more difficult to work in than a more interactive environment would be. These barriers mean that, despite its open source, only a small number of RTKLIB users actually ever dig into the core code. For those who want to take things a step further and experiment with or develop their own GNSS algorithms, RTKLIB may not be an optimal choice.

I have often thought that a Python version of RTKLIB would help minimize some of these barriers and make RTKLIB more useful as a learning and development tool. Recently I was made aware of a newly developed Python version of the code through a reference in the diary of Tomoji Takasu, the creator of RTKLIB. The code is called CSSRLib, and was written by Rui Hirokawa from Mitsubishi Electric. It is primarily intended to demonstrate the use of SSR correction data but includes a port into python of the pieces of RTKLIB needed to run PPK and PPP solutions. It also includes example data and wrapper scripts to run these solutions on the sample data. It was obviously a big effort and a significant accomplishment.

Although the CSSRLib package could be quite useful as-is for exploring precision GNSS solutions, the translation from C to Python is not strict enough to allow a user to jump directly back and forth between this code and RTKLIB.

I thought it would be an interesting exercise to rewrite the python code to be more closely aligned to RTKLIB, add all of the changes and enhancements in the demo5 version, and try to match the solutions of the demo5 solutions as closely as possible. This turned out to be a significantly more time-consuming proposition than I realized, but after several months of on-and-off work on the project, I finally was able to complete a code that, for the most part, meets these criteria.

It is not intended to be a substitute for the C version of RTKLIB since it only performs a small subset of the full library capability and runs noticeably slower. It currently runs only PPK solutions, although I would like to add PPP solutions later. The CSSRLib package supported the GPS and Galileo constellations. I have added support for GLONASS, but Beidou is still not supported. Since the resulting code is quite different in implementation and purpose from the original CSSRLib, I have created a new repository on Github for it, and named it rtklib-py. The code still shares many of its CSSRLib roots and so I have left the original copyright notices in the code files and have added acknowledgements to the original code.

The CSSRLib package includes an example dataset from a geodetic quality rover but since the demo5 code focuses on low-cost receivers, I have replaced the sample data with two other datasets, one from a u-blox F9P rover mounted on a vehicle roof, and a second dataset from the Google Smartphone Decimeter Challenge containing data from a smartphone mounted inside a vehicle.

A truly literal translation of the C code would run very slowly, primarily because Python is an interpreted language and C is a compiled language. To make the python code run at a reasonable speed, many of the for loops in RTKLIB have been replaced with Numpy array operations. This need to optimize, along with inherent differences between C and Python, means the codes are not identical but I have attempted to make the two codes as similar as possible in file names, function names, variable names, logic, and even comments. I also added a similar trace debug feature to the code that, when enabled, produces trace files that are very similar to the RTKLIB trace files. These will align quite closely when compared against each other with a file compare app and provides a way to confirm that the intermediate results match between the two code sets.

Here is an example of the trace file from RTKLIB on the left and rtklib-py on the right which demonstrates how similar they are. There are still small differences in the codes, which along with the iterative nature of the solution does cause the final positions to diverge by small amounts, but the intermediate results are for the most part very close.

The inputs (rinex files) and outputs are the same between the two codes, so RTKCONV and RTKPLOT in RTKLIB can be used to generate the rinex input files and plot the solution output files.

The config parameters in the new code have names and functionality very similar to the RTKLIB code although not all options are supported, particularly the ones used more frequently in PPP solutions than in PPK solutions. For the two sample data sets, the config parameters are defined in f9p_config.py and phone_config.py respectively. The top level script to run a solution is run_ppk.py. By default, it is setup to run the F9P sample data set but includes a commented out section to run a solution for the smartphone data. You can run a solution on your own data by modifying the files specified in run_ppk.py. You may also want to create a copy of one of the existing config files and adjust it for your data.

Below is a comparison of an RTKLIB forward solution to the F9P example data set on the left and an rtklib-py solution on the right, using the same configuration parameters for both.

Here is the difference between the two solutions. As you can see, they are quite close but not exactly the same.

Rtklib-py is not meant to be a replacement for RTKLIB and would not be a good choice for someone who is interested only in the final results. However, I am hoping it can be useful in several other possible ways:

  1. As a “map” to explore and learn how the inner details of RTKLIB work
  2. As a development environment to experiment with enhancements or adjustments to the RTKLIB algorithms. Due to the close alignment between the two packages, these can then be fairly easily ported back into the C/C++ version of RTKLIB when they are complete
  3. To debug issues found in RTKLIB within a more interactive environment
  4. To cut and paste pieces of the code into more custom solutions

One last goal is for this code is for it to be available as a tool for teams competing in this year’s Google Smartphone Decimeter Challenge competition.

The purpose of this post was just to introduce the new code so I won’t go into any more detail here. There are some brief instructions in the readme file in the Github repository for running the code but I do assume users are already reasonably comfortable with running Python applications.

I personally like to run Python in the Spyder IDE which provides an easy-to-work-in interactive environment and includes Numpy and other popular libraries in the installation. I’m sure, however, that there are other good development environments as well, if you prefer another option.

I’m always interested in other people’s thoughts on these topics so please leave any comments or suggestions you have in the comment section below.

Design a site like this with WordPress.com
Get started