Define metrics for evaluating the quality of a solution

So far, the improvements in the quality of the solution with each change have been quite obvious in the resulting plots from RTKLIB. Going forward though, the improvements are going to be more subtle and it will be more difficult to evaluate the changes.

Let us define some metrics we can use to help evaluate these future changes.

There are many metrics we could come up with, but I will start with these:

1) Stdev of xy distances between base and rover: This is a standard deviation of the distance between receivers. The distance between receivers should remain constant regardless of the change in orientation as discussed in a previous post, so ideally the standard deviation should be zero.

2) Median AR (ambiguity ratio) between best solution and 2nd best solution: This is a measure of the confidence of the solution coming from the integer ambiguity resolution algorithm. In general higher is better, but there is a trade-off with number of satellites used. A lower AR with a large number of satellites may be preferable to a high AR with only a small number of satellites

3) Ratio of epochs with fixed solution to total epochs: A measure of how often the result of the integer ambiguity resolution algorithm exceeds the threshold for using a fixed solution. Higher is better provided there are no significant errors in any of the fixed solutions. For integrity, we rely on a fixed solution indicating a very high confidence in the solution.

4) Median number of satellites used in float solution: This is the total number of satellites used to determine position. A higher number of satellites gives us more confidence in the result provided, if all other metrics are equal. Including low quality satellites in the solution will increase this metric while decreasing other metrics.

5) Median number of satellites used in fixed solution: This is the number of satellites considered high enough quality to use for the fixed solution. Similar trade-offs apply as with the number of satellites used in the float solution.

6) Fraction of incorrect fixed solutions: Number of erroneous fixed solution epochs divided by total number of fixed solution epochs. Any fixed solutions with significant errors dramatically reduces our confidence in the answer, so this metric should always be zero or very close to zero. We need to define a threshold for this metric based on how much error we consider acceptable in the solution.

Some of these metrics will probably prove to be more valuable than others, but I will start by monitoring all of these until have a better understanding of which are most important. Some of these metrics will tend to move in opposite directions. For example, if the algorithm excludes all but the best quality satellites, the median AR and ratio of fixed epochs will increase because the quality of the data is better, but the number of satellites used will decrease. Increasing the number of satellites used will increase the confidence that the solution is correct but may reduce some of the metrics because of the inclusion of lower quality data.

All of the metrics except #5 can be derived from data output to the status file, which is generated when RTKLIB calculates a solution. Metric #5, the satellites used in the fixed solution, requires a code modification to make a small change to the format of the status data. The code currently calculates the solution status for each satellite separately, but at the end of each epoch, changes the solution status of all satellites to the status of the solution. In other words, if five of eight satellites are used to determine the fixed solution, but the fixed solution does not meet the AR threshold criteria, the solution status for all 8 satellites are set to float. The code modification comments out one line from the end of the relpos() function in rtkpos.c as shown below in green to retain the individual status for each satellite. This change does not affect the functionality of the code in any way since this these variables are not used again after they are changed.

for (i=0;i<MAXSAT;i++) for (j=0;j<nf;j++)  {
      // Don’t lose track of which sats were used to try and resolve the ambiguities
     // if (rtk->ssat[i].fix[j]==2&&stat!=SOLQ_FIX) rtk->ssat[i].fix[j]=1;
     if (rtk->ssat[i].slip[j]&1) rtk->ssat[i].slipc[j]++;
}

I will track these metrics as I make changes in input parameters and algorithm to quantify any improvement. To add further insight, I will track the three different conditions in my input data separately. Those are static, moving slowly with open skies, and moving quickly with partially obstructed skies.

Here is the plot of the metrics for the changes we have made so far in the three previous posts for the period when the receivers are moving.

metrics1

Case 1 on the x-axis represents the baseline parameters, case 2 is after increasing eratio1 from 100 to 300, case 3 is after fixing the bug in lock count with arlockcnt=20, and case 4 is after increasing arlockcnt to 120.

The metrics in the first three plots are all moving in the desired direction. The most important, stdev of the error is reduced from over 15 cm to less than 0.5 cm indicating a significant improvement in accuracy. Increases in % fixed solutions and median AR ratio both indicate increased confidence in the solution. In the last plot, the large difference between number of satellites used for the float solution and number of satellites used for the fixed solution is because we have disabled integer ambiguity resolution for the GLONASS satellites. This suggests a large amount of information is not being used and represents opportunity for future improvement. The slight drop in number of satellites used for the fixed solution in cases 3 and 4 occur because we are not using position data from satellites that have lost lock until the kalman filter has had some time to reconverge.

Plots for the other two environments, (receivers stationary, and receivers moving quickly) show similar trends.

For all three cases after the first improvement, the last metric, fraction of incorrect fixed positions, is zero for the period when the receivers were moving, which is what we want to see. I defined an erroneous position as anytime an error in a fixed solution exceeded 5 cm. Note that for the period of time when the receivers were moving quickly, this value was non-zero and hence this setup probably would not be adequate for that environment.

I’m interested in what metrics other people are using to evaluate the quality of solutions? Please leave a comment if you have any thoughts or other metrics that might be more useful.

Improving RTKLIB solution: (AR lock count and elevation mask)

In the previous post, after a couple changes to the input parameters, we ended up with the following solution for the difference between two stationary, then moving Ublox receivers:

During the time when the receivers were stationary (17:25 to 17:45), we were able to hold a fixed solution (green in plot) continuously after the initial convergence, but started to revert to float solution (yellow in plot) fairly frequently once the receivers started moving.

Let’s look at the data a little closer to see if we can discern why we are reverting back to the float solution. If we switch to the “Nsat” view in RTKPLOT plotted below, we can see a couple of useful plots.

sol9c

The top plot shows the number of satellites used for the solution and the bottom plot show the AR ratio factor. The AR ratio factor is the residuals ratio between the best solution and the second best solution while attempting to resolve the integer cycle ambiguities in the carrier-phase data. In general, the larger this number, the higher confidence we have in the fixed solution. If the AR ratio exceeds a set threshold, then the fixed solution is used, otherwise the float solution is used. I have left this threshold at its default value of 3.0.

Notice in the plot that each time we revert to the fixed solution, it happens when the number of valid satellites increases. At first this seems counter-intuitive, you would think that more information would improve the solution, not degrade it. What is happening though, is that the solution is not using the satellite data directly, but rather, the kalman filter estimate of the data. This estimate improves as the number of samples increases. This is also why the solution takes some time to converge at the beginning of the data. By using the data from the new satellite before the kalman filter has had time to converge, we are adding large errors to the solution, forcing us back to the float solution.

Fortunately, RTKLIB has a parameter in the input configuration file to address this issue. “Pos2-arlockcnt” sets the minimum lock count which has to be exceeded before using a satellite for integer ambiguity resolution. This value defaults to zero, but let’s change it to 20 and see what happens. While we are it, we will also change the related parameter “pos2-arelmask” which excludes satellites with elevations lower than this number from being used. We will change this from it’s default of 0 to 15 degrees.

Unfortunately, re-running the solution with these changes makes almost no difference as seen below. The jumps back to float solutions still occur immediately after a new satellite is included, and the additional satellites are being included at the same epochs. Something is clearly wrong.

 

Digging into the code a bit reveals the problem. First, the arlockcnt input parameter is copied to a variable called “rtk->opt.minlock” Then the relevant lines of code from rtkpos.c are:

rtk->ssat[sat[i]-1].lock[f]=-rtk->opt.minlock;

if (rtk->ssat[i-k].lock[f]>0&&!(rtk->ssat[i-k].slip[f]&2)&&
     rtk->ssat[i-k].azel[1]>=rtk->opt.elmaskar) {
     rtk->ssat[i-k].fix[f]=2; /* fix */
     break;
}
else rtk->ssat[i-k].fix[f]=1;

So far, everything looks OK, but when we look at the definition for the structure member “lock” in rtklib.h we find:

unsigned int lock [NFREQ]; /* lock counter of phase */

This won’t work. We are assigning a negative number (-rtk->opt.minlock) to an unsigned variable, then checking if the unsigned variable is greater than zero, which by definition, will always be true. Hence, satellites are always used immediately, regardless of lock count.

We fix this bug by changing the definition of lock from unsigned to signed:

int lock [NFREQ]; /* lock counter of phase */

Rerunning after rebuilding the code, gives us the following solution:

Much better! All but a couple of the jumps back to float solutions have been eliminated.

Looking at the plot of the AR ratio for this solution, we can see that the kalman filter requires more than 20 seconds (the value we chose for arlockcnt) to re-converge. Two minutes looks like a more reasonable value from the plot, so let’s try rerunning the solution with arlockcnt=120.

Even better.  That eliminates the last couple of jumps back to the float solution and also significantly reduces the number of jumps in the AR ratio plot as shown below. The green line in the bottom plot is arlockcnt=20 and the blue line is arlockcnt=120.

sol12c

It’s possible we will need to revisit this number after looking at more data, but for now we will use 120 secs.

As always, please leave a comment if you have any comments, discussion points, or questions.

Update 4/11/16:  

  • In the above discussion, my data is one epoch/sec and I equate epochs with seconds.  Arlockcnt is specified in epochs, not seconds, so if your data is at a different sample rate you will need to adjust accordingly.
  • The suggestion above to set arlock to 120 secs is probably too conservative for more challenging environments.  If there are few satellites and/or many cycle slips for example it is probably better to start using a just-locked satellite before it is fully converged.  I have started using 30 secs rather than 120 secs for my solutions.
  • You can pull this code fix from my github repository either by itself from the arlockcnt branch or as part of a larger set in the demo1 branch

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.

 

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.