RTKLIB solution accuracy

OK, so let’s say you’ve just gone out and collected some decent raw measurement data, everything looks like it is working properly and you get a nice-looking solution from RTKLIB, maybe something like this.

stdev1

The next thing you might want to know is how accurate is that solution?  Fortunately, RTKLIB does include accuracy estimates in the solution file.  These are expressed as standard deviations, one for each axis.  Here’s an example for a solution calculated in ENU (east/north/up) coordinates.  The sde, sdn, and sdu columns are the standard deviation estimates for each solution point and are given in meters.  If the solution were calculated in XYZ coordinates, the columns would be labelled sdx, sdy, and sdz.

stdev5

You can display these on the solution plot with RTKPLOT by setting “Error Bar/Circle” in the options menu to “Dots” or “Bar/Circle”.  You can see the gray lines in the plot below that represent the solution value plus and minus one standard deviation.

stdev2

Visually, I prefer to see three standard deviations plotted, rather than just one.  Three standard deviations represents 99.7% of the distribution, so the correct value should nearly always appear between the two gray lines.  One standard deviation, on the other hand,  represents only 68% of the distribution which is hard for me to attach any physical meaning to.  This can be easily changed in RTKLPLOT by replacing SQRT(s[j]) with 3*SQRT(s[j]) all three places it appears in the plotdraw.cpp file and recompiling RTKPLOT which gives this:

stdev3

Of course these are just estimates of the solution uncertainty made by the kalman filter in RTKLIB and we can not assume they are accurate without some analysis.  They are a function of the input parameters that set the input measurement uncertainties and process noises for the kalman filter.  These are all in the “stats” section of the input configuration file.  Even if these are all set to correctly match your actual configuration, there will be additional errors in these estimates caused by non-linearities in the system, non-gaussian distributions, and time-correlated measurements among other things.

So, let’s do a little analysis to get a feel for how useful these estimates are.  We’ll start with the above data set which includes eight hours of measurements from a stationary u-blox M8T receiver and a CORS reference station 8 km away.  I enabled fix-and-hold for ambiguity resolution with relatively low tracking gain (pos2-varholdamb=0.1).  I ran both static and kinematic solutions on the data, then turned off ambiguity resolution and re-ran the solutions to get examples of both fixed and float solutions for static and kinematic modes.

After running the solutions, I pulled the solution files into matlab for plotting and analysis.  Since the rover was stationary and I knew it’s precise location I was able to easily calculate the true errors.  I then plotted for each solution point, the absolute value of the errors and three times the standard deviation estimate for each point, similar to the above RTKPLOT plot.  The raw data is in E/N/U coordinates but I combined the East and North values into a single horizontal error and error estimate.  After finding that in some cases the error estimates were too low, I also plotted a line using double the error estimate, to see if that might be a usable value .  In the title of each plot I’ve included the percent of samples that were inside the error lines, the first number is for the undoubled estimate, the second is for the doubled estimate.  I also made one more modification to the error estimates by limiting them to a minimum value.   In the case of the static solutions, the RTKLIB error estimate will continue to decrease almost all the way to zero but this is not realistic since there are errors that are not visible to the solution.  I arbitrarily limited each horizontal axis standard deviation to no less than 5 mm, and the vertical axis to 1 cm.   I don’t have good justification for these particular numbers but do believe there needs to be some limits set.

For a perfect Gaussian distribution, plus or minus three standard deviations should include 99.7% of the data.  Given that these are not perfect distributions and will also contain outliers, I would consider the estimates useful if more than 95% of the data fell within them.

Here are the plots for the above data set run with static solutions, float on the left, and only the fixed values of the solution with ambiguity resolution enabled on the right.  The blue line is the actual error, the red line is the RTKLIB estimate, and the yellow line is double that.  The top plots are for the horizontal axes and the bottom plots for the vertical axis.  In the static cases, all the error lines dropped below the minimum estimate values I set above fairly quickly, so only the initial convergence of the float solution was really of any interest.   Still the error estimates in that region look somewhat reasonable, although the undoubled estimates are clearly optimistic. The undoubled estimates (red) are not visible on the plots where both estimates have dropped below the minimum and they get overwritten by the doubled estimate.

stdev8

Next are the same plots for the kinematic solutions.  Just to confuse you, I have accidentally put the float solution on the right this time.  Again, the solution with ambiguity resolution enabled (on the left) only includes points that were fixed.  These are more interesting than the previous plots since they are not dominated by the arbitrary minimum error estimate limits I added.  Looking at the float solution first, the undoubled estimates again look too optimistic with only 92.2% and 75.2% of the data within that limit.  At 100% and 99.6% the doubled values (yellow lines) appear to be better estimates.  For the fixed values on the left however the undoubled estimate percentages and plot lines appear to be fairly reasonable.

stdev7

In one last example, I used a different data set, this one is one of the sample data sets on my website which was taken with two M8T receivers, one on each end of a kayak out on the ocean.  In this case we know that the distance between the two receivers will remain constant so we can use that information to do a similar analysis as above.  Because the kayak is moving in all three axes with the waves we can not separate the horizontal and vertical components but can do an actual error to estimated error comparison for the three axes combined.

The float solution is back on the left again, just the fixed values from the ambiguity resolution enabled solution on the right.  Based on the calculated percentages and the shapes of the curves, it seems fairly reasonable to again choose the doubled estimates for the float solution and the undoubled estimates for the fixed values.  Even with the doubling, the percent of solution values within three standard deviations is still only 88.5% which is rather low but most of the outliers were in the initial convergence where we might expect more issues.

stdev6

So, based on these few examples I do see correlation between the RTKLIB error estimates and the actual errors, not as much as I would have liked, but there is some.   Maybe these adjustments I made would hold up for other cases as well, but I would not have a lot of confidence without trying more examples first.

Clearly these error estimates can not be taken at face value.  Do they have any value?  It probably depends on what you are trying to do with them.  If you are just trying to get a rough idea of the accuracy of the measurements, especially if it is in a statistical sense, and you have done some previous analysis with similar data sets, then you can most likely derive some value from them.  If the absolute accuracy of a single measurement is important to you then you will probably need to find another solution.

What about trying to adjust the input configuration parameters (measurement and process noises) to make the accuracy estimates more accurate?  You may be able to make small improvements but I suspect they will not be significant enough to avoid post-solution adjustments.  You will also find that adjusting these parameters will affect the quality of the solution and it will be difficult to optimize for both.

I am interested, though, if anyone has any ideas for simple (or at least only moderately complicated) code changes that could be made to RTKLIB to improve these estimates or any other ideas on how to improve them.

 

 

 

Advertisements

4 thoughts on “RTKLIB solution accuracy”

  1. I think the easiest thing to start with is to calculate a scaling factor in Matlab as follows (I’ll show for the vertical component only as 1D is easier to handle):

    s = sz;  % RTKLIB output uncertainty, in meters, vector
    r = z - mean(z);  % vertical residuals, in meters, vector
    
    s0 = sqrt(mean(r.^2./s.^2))  % reduced chi-square stat, unitless, scalar
    %s0 = 1.4826*median(abs(r./s))  % outlier robust version of s0.
    
    sms = s0.*s;  % scaled uncertainty of the mean
    sps = sqrt(s0.^2 + ssm.^2);  % scaled uncertainty of instantaneous positions
    

    Then expand the uncertainty for a given probability:

    cl = 99/100;  % confidence level
    sl = 1-(1-cl)/2  % two-tailed significance level
    smse = sms*norminv(sl, 0, 1)  % expanded uncertainty of the mean
    
    sl2 = sl/length(r);  % Bonferroni correction
    dof = 10;  % degree of freedom, vector -- ideally RTKLIB would output this
    spse = sps*tinv(sl2, dof)  % expanded uncertainty of instantaneous positions
    

    A second issue is what type of uncertainty you need: uncertainty of the mean or uncertainty of an instantaneous position. For the distinction between the two, see page 39, Section 3.5, “Confidence intervals for predictions”, in Faraway (2002), “Practical Regression and Anova using R” (available online).

    The expanded uncertainty of the mean has the following interpretation: if you run multiple sessions of the same duration (eight hours) under similar conditions (same base station and time of day), about 1% of the session means will fall outside the expanded and scaled uncertainty of the mean (at the end of the session). The 1% will be approached over many sessions (asymptotically).

    It’s only the instantaneous uncertainty that is supposed to bound the epoch-by-epoch position residuals. And even so, it is more difficult to get it right, as it depends on more assumptions, including the number of raw measurements per epoch (which AFAIK is not output by RTKLIB), and the correction for multiple comparisons (Bonferroni, Sidak, etc.).

    I think a practical recommendation to get realistic uncertainty estimates is to break up the session into sessions as small as possible (say, soon after convergence) and use the repeatability over individual sessions to infer the noise level. For example, having eight one-hour sessions in this case.

    Finally, errors with respect to truth may include hidden systematic errors which will remain undetectable, such as a blunder in the base station coordinates.

    Like

  2. Hi Tim,

    There is a paper by Han and Rizos (1995) titled “Standardization of the variance-covariance matrix for GPS rapid static positioning” that I find quite useful for scaling the covariance matrix in static mode. If you Google it, you will find a link to download the paper.

    The paper presents a scaling factor that takes into account the data sampling interval, the session length and an approximate correlation time for measurement errors. It is easy to compute and provides decent results under normal conditions.

    I don’t think that the scaling factor fully applies to kinematic positioning. It could provide more realistic standard deviations during the convergence period, as long as you don’t have discontinuities in your solution.

    I hope you find it useful!

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s