Improving RTKLIB solution: (Receiver Dynamics and Error Ratio)

In the previous posts, I demonstrated reasonably clean looking solutions using the default RTKLIB configuration options for rover data from a roving M8N receiver and base data from a stationary COREX base station.

I also showed a second solution for the slightly more challenging problem where roving M8N receivers are used for both the base and rover inputs. I mentioned that run was not done with the default settings but did not go into the details. Let us now go back and try that run using the default settings and see what happens. Remember what we hope to see is a circle with radius equal to the separation between the two receivers, ideally plotted in green, indicating RTKLIB was able to resolve the integer ambiguities and provide a fixed solution.

For reference, here is the solution we previously calculated with RNX2RTKP using non-default settings and code. The plot on the left is from the ground track option in the RTKPLOT menu which plots the x axis vs the y axis and the second plot is the position option which plots all 3 directions separately vs time.

Below is what we get running with the default settings with the two basic modifications I mentioned earlier. The red indicates the kinematic solution is unable to converge and the plotted result is done in single mode, which is not differential, and therefore is indicating absolute position, and not distance between the two receivers.

Clearly, the default settings are not adequate when both receivers are lower quality and both are moving. Let’s see if we can improve this.

First, let’s make a few small changes that won’t improve the solution but will make the input assumptions better match our problem.

Edit the input configuration file for rnx2rtkp to make the following changes. These are in addition to the two changes we made earlier (pos1-posmode and ant2-postype)

  1. pos1-frequency: L1+L2 → L1, both receivers are L1 only, so no need to look for L2 data
  2. pos1-navsys: 1 → 5, we have GLONASS data, so let’s use it
  3. pos2-gloarmode: on->off, not ready for this yet, in current state, it prevents fixed solution
  4. pos2-bdsarmode: on->off, no Bediou data so disable
  5. out-solformat: llh → enu, save solution to output file with equal units in x and y axis
  6. out-outstat: off → residual, save residuals for later analysis

Rerunning RNX2RTKP gives us a solution very similar to the previous one with little improvement, even with the additional data from the GLONASS satellites.  I won’t bother to plot it here.

Next, let’s turn on receiver dynamics with the “pos1-dynamics” input option. This adds nine states to the kalman filter, 3 for receiver position, 3 for velocity, and 3 for acceleration. This enables the solution to take into account the previous position of the receiver when calculating the new solution and allows us to better define likely changes in position, velocity, and acceleration. For now, we will use the default filter settings. With this change, the solution looks like this:

Not quite there yet, but much better. We can now see the expected circle but there are fairly large errors when the solution is not fixed. We do sometimes resolve the integer ambiguities but only 26.6% of the time.

Next, we’ll take a very brief look at the input parameters defining the error characteristics of the input data. For the kalman filter to effectively combine the data from multiple inputs, it requires information about the error distribution of each input.  RTKLIB calculates the variances of each code and carrier phase input based on a set of input parameters and a simple model based on satellite elevation, the assumption being that lower elevation satellites are noisier. For now, we will look at only the first input parameter. This is “stats-eratio1” and it sets the ratio of standard deviations of the pseudorange errors to the standard deviations of the carrier-phase errors. I believe it is not unreasonable to expect, with the lower quality receiver and antenna we are using, that the pseudorange errors will degrade more quickly than the carrier-phase errors. If this is true, then it would make sense to increase this ratio from the default of 100. Setting this to 300 give the following solution.

This looks a lot better. We converge much more quickly to a fixed solution and when it loses lock it also re-converges more quickly.

We still need to make a few more changes to get to the quality of the solution at the top of this post, but most of those changes require making some code changes so I will leave them to another post.

The physical justification for using 300 for the error ratio is a bit weak at this point but we’ll just go with it for now and I hope to come back and do a more complete analysis of the way RTKLIB generates the variances and covariances for input to the kalman filter later.

Anyone else come up with a good way to optimize the  RTKLIB input parameters for the kalman filter convariances?  If so, please leave a comment, I’m interested in how other people handle this problem.

Update 6/12/16:  Since writing this post, I have found that increasing eratio1 can have a tendency to increase the number of false fixes if “pos2-rejionno” is not also increased at the same time.  I now recommend that if you increase eratio1 that you also increase rejionno from 30 to 1000.  This should eliminate false rejection of outliers.  For more details see this post.

 

Advertisements

Building RTKLIB 2.4.3 code with Visual C++ 2010

One of my goals in this project is to make some changes to the RTKLIB solution algorithms and evaluate them specifically for my set of inputs. RTKLIB is intended to work in many different environments. I would like to customize it to be optimal for calculating differential solutions for two single frequency receivers assuming low velocities and open skies. Some of the changes I would like to make are based on what other people have already tried but for which code is not available. Other changes are based on looking at cases in my data where RTKLIB did not provide a good solution, understanding why, and changing the algorithm appropriately.

To make these changes, I need to be able to modify the RTKLIB code and rebuild it. Since I am doing this initial evaluation work on a Windows PC I will need to build the code in that environment. The GUI versions of the RTKLIB programs require the Embarcadero C++ compiler which I do not have and which is fairly expensive to purchase. The CUI versions can be built with Visual C++, which is available for free, so I have chosen to go that path. As mentioned before, I also find the CUI interface with a matlab wrapper better for tracking configuration information for multiple runs.

The project files for Visual C++ are in the rtklib\app\”app name”\msc folders. They are configured for Visual C++ 2008. Since I already have Visual C++ 2010 loaded on my machine I had to make a few changes to make the project files work. Some of the changes would also be required with VS 2008 since it does not look like the project files have been kept up to date with recent RTKLIB changes.

  1. Convert solution files to VS 2010: Click on msc.sln file in rtklib\app\rnx2rtkp\msc folder. This will bring up the visual studio conversion wizard which will convert the VC 2008 project files to VC 2010.
  2. Add new src files to project: Add tides.c and ppp_corr.c to the list of source files in the “Solution Explorer” window by right-clicking on the source folder and selecting “Add”. This prevents build errors from unfound symbols.
  3. Modify Include Directories: Select “Properties” from “Project” tab. Select “General” under “C/C++” under “Configuration Properties” in the menu on the left hand side. Replace the existing entry with “..\..\..\src” to make it independent of code path. This enables the compiler to find the rtklib.h file.
  4. Modify Target Name: Select “General” under “Configuration Properties” in the the menu on the left hand side and change “Target Name” to “rnx2rtkp” to avoid linking errors.
  5. Add winmm library: Select “Input” under “Linker” under “Configuration Properties” in the menu on the left hand side, and add “winmm.lib” to “Additional Dependencies”. This avoids linker errors for an unresolved TimeGetTime symbol.
  6. Fix LPCWSTR conversion warning: Select “General” under “Configuration Properties” in the the menu on the left hand side and change “Character Set” from “Use Unicode” to “Not Set”

To run these apps from the data folders you will need to add the paths for the executables to your Windows path variable. The executables are in app/”app name”/Release/”app name”.exe or app/”app name”/Debug/”app name”.exe depending on whether you built in Debug or Release mode. I always build in Debug mode and point my path variable to that folder just to make it easier to switch between debugging with VS and running stand-alone, but either will work.

Evaluating the quality of moving position solutions

Stationary position solutions are easy to evaluate, an ideal solution is a single point, any deviation from this point must be error. It is possible that the point itself is in error in an absolute sense but since I am only interested in the distance between receivers I am not concerned with absolute error.

Evaluating a solution for a moving receiver is more difficult since it is not easy to know what is the correct solution to compare to. If I had an expensive high precision receiver/antenna I could mount it on the same platform as the test receiver and use that as a reference, but I don’t have that, so I need to come up with another solution.

Mounting two low cost receivers on the moving platform gives us a couple of options for evaluation. First of all, since I have data from two nearby base stations, I can calculate two solutions, each using one of the low cost moving receivers as rover, and one of the COREX base stations for base data. Since the two low cost receivers remain a fixed distance apart, the difference between the two solutions should be constant in magnitude. The direction of the difference however is varying as the orientation changes. This means the difference between the two solutions should be a circle with radius equal to the distance. Since each solution is based on different base stations, and different receivers they should be fairly independent. The fact that one of the base stations is reporting GPS only, and the other is reporting GPS and GLONASS should also increase the independence of the solutions, since they are also using different sets of satellites,at least if we include GLONASS satellites in the solution. This not the case at the moment since the default config has GLONASS turned off, but we will turn it on soon.

RTKPLOT has a nice feature which allows us to plot the difference between two solutions. Click on the “1-2” button on the menu bar to see this, after opening both solutions from the file menu.

Below are the two solutions plotted on top of each other. Since the scale is 20 meters per division, it is impossible to see any difference between them.

sol3

Here is a plot of the difference between the two solutions while the car was stationary and driving in low velocity circles. Now the scale is 10 cm per division. We can see the expected circle, and the error from that circle is generally less than 5 cm once the solution converges to a fixed solution. I’ve jumped ahead a bit, since the solution in the plots includes some input parameter changes as well as some code changes, so the default solution we just calculated won’t be as clean, but we’ll get to those changes soon. The good news is that it looks like we are getting errors of only a few centimeters after the solution has converged, even when both receivers are moving.

sol4

There is also another, simpler way to evaluate the quality of the solutions if, as in my case, we are interested only in relative error and not absolute error. Instead of generating absolute solutions using the COREX base stations, we can calculate a differential solution directly, using one of the low-cost receivers as the base and the other as the rover. The result will then be the position difference between the two receivers. Since the distance between the two receivers is fixed but the orientation between them is varying as the platform they are on is moving, the solution should be a circle, again with radius equal to the distance between the two receivers. In this case, unlike the previous one, we do not get any absolute position information, only the distance between the two receivers. Note that even though the base is moving in this case, we do not use the “Moving-Base” solution option from the positioning mode choices, but stick with “Kinematic”. The solution algorithm does not need exact position of the base, only approximate location. We use the approximate starting location of the base receiver which is in the header of the RINEX observation file and ignore the fact that the base is moving. As long as we are concerned only with relative position between the two receivers and the base does not move a large distance from its starting point, this is a valid assumption to make.

Here is a plot of the result from using both base and rover on the moving platform. Note that it is very similar to the previous plot, generated from the difference between two absolute solutions, and again the error is only a few centimeters once the solution has converged.  As in the last plot, it does include some changes from the default input parameters and code that we will get to later.

sol5

This is the form of solution I will use for evaluating the effectiveness of various input parameters and algorithm changes.

Kinematic solution with RNX2RTKP

 

In the last post we used the RTKPOST GUI to generate a kinematic position solution using “ebay.obs” for the rover data and “zdv13480.15o” for the base station data. To do the equivalent run from RNX2RTKP we use the following command line, run from the folder containing the data files.

rnx2rtkp -k test1.conf -o out.pos ebay.obs zdv13480.15o ebay.nav

The -k option is to specify the configuration file, the -o option is the output file. In this case we use the configuration file we saved from RTKPOST in the previous post.

To plot the calculated solution use the following command line:

rtkplot out.pos

This calls up the same plot GUI we previously accessed by clicking on the “Plot” button in the RTKPOST GUI.

The plots should be identical between this solution and the previous solution generated with RTKPOST.

To keep track of changes between various runs with the same data, I use a matlab wrapper to create a sub-folder, copy the config file, trace file, stats file, and result file to that sub-folder, input a short description of the run, and also save that to the sub-folder.

Kinematic solution with RTKPOST

Kinematic solution with RTKPOST

RTKPOST is the GUI tool in RTKLIB to calculate position solutions. Most of the time I find the CUI version (RNX2RTKP) better fits my needs, but just to check everything is working, it is probably easier to use RTKPOST the first time.

For a demonstration of using RTKPOST to find a kinematic position solution, I will use the ZDV1 (COREX station data) for base station data and the “EBAY” (Ublox M8N receiver data) for rover data. Zipped versions of this data is available here.  Since the exact location of ZDV1 is known and is in the observation file header, the kinematic solution will give us an absolute position by solving for the relative distance between the rover and the base, and then adding that to the base location. I set up the GUI inputs as shown below to point the program to the correct observation and navigation files. If you use my data, be sure to change the paths to match where you saved the data to.

rtkpost

For this first run, to keep things as simple as possible, we will make just two changes from the default setup. Use the “Options” button to get to the options menu. Under the “Setting 1” tab, change “Positioning mode” from “Single” to “Kinematic”. This will give us a differential solution using carrier phase info instead of an absolute solution using only pseudorange. Next, under the “Positions” tab, change the first field under “Base Station” from “Lat/Lon/Height” to RINEX Header Position. This will tell RTKPOST to get the base station location from the header of the observation file.

While you are in the options menu, click the “Save” button, and save the options setup to a location you will remember later. We will use this file as the configuration input file for the CUI version. Then click OK to exit the Options menu.

Click the “Execute” button to calculate position and then “Plot” to see the solution. Select “Gnd Trk” and zoom in and it should look like this. The two rectangles are parking lots. The yellow represents a float solution, the green a fixed solution. The fact that we were able to get a fixed solution at least part of the time is a sign things are working reasonably well.

sol1

 

Zooming into the initial time period when the car was stationary we see the plot below. Since the receiver is not moving during this time, any movement in the solution represents error. During the initial convergence of the kalman filter we see quite a lot of error, but once it does converge, we do get a fixed solution for 5 of the 20 minutes which appears as green in the plot below. During this time you can see the error is roughly +/- 1cm in the xy direction which again is a good sign things are working.

sol2

Note about restoring RTKPOST default options:  There is no button in RTKPOST to reset the options to defaults and it remembers the options from the previous session when restarted, so there is no obvious way to put it back to defaults.  The best way I have found to do this is to delete the “rtkpost.ini” file saved in the rtklib\bin folder before starting RTKPOST.

 

My reference raw GPS data set

In many of the posts to follow I hope to compare the effectiveness of different input options and algorithm modifications to the position solution. To do this it will be useful to have a single representative set of raw data to use for comparisons as well as a set of metrics to compare. I will discuss the metrics in a future post, here I put together a data set that I believe will be a good test, at least for my goals.

My interest is in relatively low velocity scenarios with open skies (I’ll talk more about the details in another post). To generate my reference data sets, I mounted two Ublox M8N receivers with basic patch antennas on top of a car, both antennas sitting on the metal roof for a decent ground plane. I then collected about twenty minutes of data with the car stationary, followed by driving around in circles in a open parking lot (no nearby trees or buildings) for roughly another twenty minutes, followed by another twenty minutes of driving on open roads with some nearby trees at higher speeds. The low velocity circles are what I am most interested in, but collected the other data to allow additional comparisons under different conditions.

I am fortunate enough to have two CORS base stations within 10 km of where I collected my data and so I downloaded the data for both stations (ZDV1 and TMGO) for the same time period. ZDV1 is GPS data only, TMGO includes both GPS and GLONASS data.  That gives me four sets of data. Two are from Ublox M8N receivers with patch antennas, both moving but staying a constant distance from each other, both containing only single frequency (L1) data. The other two are dual-frequency (L1 and L2), fixed position. For consistency, I will use this set of data when comparing results in most of my future posts. Here is a link to a zipped copy of all four data sets in case they may be of use to anybody else. In the file names, “cga” refers to data from the M8N receiver bought from CGA, and “ebay” refers to data from the M8N receiver bought on Ebay.

Base station data for RTKLIB

RTKLIB has a number of different algorithms it can use to calculate position. The first two in the list are “Single” and “DGPS”. Both of these methods use only the pseudorange data and not the carrier phase information. Without the carrier phase information, the precision of these methods is quite low, and probably worse than what the GPS receiver would provide without RTKLIB. In general they are not very interesting other than possibly for some initial debug during setup. The rest of the algorithms can be divided into two groups, differential and PPP. The differential algorithms determine position relative to a known nearby location while the PPP (precise point positioning) algorithms determine absolute position. In general, the data quality of the low cost hardware we are using is only good enough to use with the differential methods and not the PPP methods, so we will focus on those. RTKLIB supports four differential modes: Static, Kinematic, Moving-Base, and Fixed. The kinematic mode is designed to calculate the relative position between a fixed base and a moving rover and that is what we will use. It can also be used with a moving base if we are concerned only with finding the distance between the two and not trying to translate that to a fixed position, and we will use it this way as well.

In my particular application, I am interested only in the distance between two receivers and not in any absolute locations. In this case I will collect data from two low-cost receivers and use one as the base and the other as the rover. In many cases, though, it is useful to use an existing ground station as the base and the low cost receiver as the rover. I use this method for testing and verification, even though I don’t plan to use it in my final solution.

In the US, base data is available online for free from many GPS stations in the CORS network.  Here’s a map of CORS station locations.

cors_map.PNG

For the ground stations I have used, it is generally available online less than an hour after it is taken. It is fairy easy to pull it manually using a user-friendly form from the CORS webpage, or you can use the RTKGET utility in RTKLIB. Be aware that some stations have only GPS data, and some have both GPS and GLONASS data.