RTKLIB: GLONASS Ambiguity Resolution

Several readers now have suggested that adding inter-channel bias calibration to RTKLIB would be a worthwhile enhancement so I thought I would give it a shot, or at least explore the idea.

Inter-channel biases are what make integer ambiguity resolution more difficult on the GLONASS satellites than on the GPS satellites. In the GPS system, double-differencing the carrier-phase observations is a very powerful trick that causes the largest sources of error to all cancel each other out. These errors include the satellite and receiver clock biases, as well as the majority of the ionospheric and tropospheric delays and ephemeris errors. The remaining errors are quite small, meaning that the estimates of the unknown phase-biases we are trying to solve for will be quite close to integer values. We can then use integer ambiguity resolution to select the most likely set of integers, thus rejecting most of the remaining error.

Unfortunately, unlike the GPS system in which all satellites operate at the same frequency, in the GLONASS system, the satellites use a set of multiple frequencies. The different frequencies introduce an additional set of errors from the receiver into the solution, called the inter-channel biases. There will be a different bias for each frequency but the error between two satellites should be roughly proportional to the frequency difference between the two satellites. For a more complete discussion, see here

A quite detailed description of at least one method to solve for the inter-channel biases can be found here. Fitting this solution into the existing RTKLIB architecture, however, looks to be quite challenging and more than I really want to take on, but I suspect this might be the “right” way to do it.

Instead, I decided to try a “feedback short-circuit” similar to what is is used in the RTKLIB fix-and-hold feature. If you haven’t already done so, you might want to read my previous post in which I discuss that feature and a modification I made to it I call “fix-and-bump” which I will use here.

Fix-and-hold uses the fact that we know the correct answer must be a set of integers to push the phase-bias estimates towards integer values, using the difference between fixed and float solutions as feedback, rather than deriving it from physical principles. My premise is that just as we can blindly adjust the phase-biases estimates based on the difference between fixed and float solutions, we ought be able to do something similar for the inter-channel biases.

The adjustments made to the phase-biases during fix-and-hold are done in the holdamb() function in the rtkpos.c file in RTKLIB. This function is called only when we have a high confidence in our fixed solution, to avoid chasing an incorrect solution. The key line of code from this function is

    v[nv]=(xa[index[0]]-xa[index[i]])-(rtk->x[index[0]]-rtk->x[index[i]]);

where:

v[] is a pseudo-innovation (or pseudo measurement residual) used to update the     kalman filter
xa [] is the set of fixed solution phase-bias estimates
rtk->x[] is the set of float solution phase-bias estimates
index[i] is the satellite to update
index[0] is the reference satellite.

So what this line of code does is generate a double-differenced error input to the kalman filter for each satellite based on the difference between the fixed solution phase-biases and the float solution phase-biases. This will drive the phase-bias states towards a set of integers, with all the advantages and disadvantages as discussed in the previous post.

So far I have just reviewed what is in the existing RTKLIB code. Now lets talk about an addition to the fix-and-hold code for the GLONASS satellites. It is reasonable to drive the phase-bias estimates for the GPS satellites towards integer values because we have accounted for all the significant errors (mostly through cancellation in the double-differences) and so we expect them to be integers. On the other hand, we can not expect the phase-bias estimates for the GLONASS satellites to converge towards integer values since we have not accounted for the errors from the inter-channel biases. These are different for each satellite and hence do not cancel. So, in order to account for them, I add an additional variable for each satellite to hold the estimated inter-channel biases. These values are subtracted from the measurement residuals every epoch when calculating the kalman filter innovations, just as the other error terms are. To do things properly I should make them kalman filter states, but for now they are just stand-alone variables.

The inter-channel bias estimates are initialized to zero and then only adjusted during a fix-and-hold update. The new update is done after the existing fix-and-hold updates are made to the phase-biases for the satellites with fixed solutions. At this point, we assume that all the remaining error in the GLONASS phase-bias estimates is coming from the inter-channel biases, so we move that error from one to the other. To make the new update, the fractional components of the GLONASS phase-bias estimates are removed (subtracted) from the phase-bias estimates and moved (added) to the inter-channel bias variables. Ideally we would subtract the fixed solutions from the float solutions as we did with the GPS satellites, but since we don’t have fixed solutions yet for the GLONASS satellites, we use our best guess for the fixed solution, which is just the closest integer. In a final implementation, this process would be iterated over multiple epochs to improve the estimates but at this point I am only making a single adjustment to simply evaluation of the solution.

These are the key lines of code added to holdamb() to do this:

dd=rtk->x[IB(j+1,f,&rtk->opt)]-rtk->x[IB(i+1,f,&rtk->opt)];
dd=dd-ROUND(dd); // throw out integer component
rtk->x[IB(j+1,f,&rtk->opt)]-=dd; // remove fractional part from phase bias
rtk->ssat[j].icbias[f]+=dd; // and move to IC bias

With the inter-channel biases now accounted for, we can turn on integer ambiguity resolution for the GLONASS satellites (gloarmode=1) Since we have forced the phase-bias estimates to integers, the integer ambiguity resolution criteria will always be satisfied initially, regardless of the quality of the estimates. If we then let the system run without any more fix-and-hold updates, meaning we have removed any preference in the solution calculations for integer results, the solution should hopefully either stay converged if it correct, or drift away from integer values if it is not. To do this,I modified fix-and-hold to occur only once (“fix-and-bump”) instead of every epoch, as described in my previous post.

So how well did it work? Below are some plots taken after making the code changes described above. The first plots are for a baseline taken with GPS ambiguity resolution enabled and GLONASS disabled, as RTKLIB is typically run. I have modified the meanings of green and yellow in the RTKPLOTs of residuals. Green now represents a satellite that was successfully used for ambiguity resolution, yellow represents a satellite that was valid but not successfully used for ambiguity resolution. If no fixed solution was found, then all satellites will be yellow. The GPS satellites are on the left, the GLONASS satellites on the right. You can see the GPS satellites are yellow during the initial lock-on and occasionally after a cycle slip but otherwise stay green, while the GLONASS satellites are always yellow because ambiguity resolution is disabled for GLONASS.

Now I make a second run with ambiguity resolution disabled for GPS (left) and enabled for GLONASS (right). I do this only for evaluation purposes, normally we would run with both enabled. Note that we still require GPS ambiguity resolution to calibrate the inter-channel biases which shows up below at the start of the plots.

As you can see, the solution remains fixed when using only the GLONASS satellites for ambiguity resolution for a significant period of time and the residuals remain small. This is very encouraging.

To further verify things are working the way we think they are, lets look at the output traces for three runs, GPS AR only, GLONASS AR only, and GPS+GLONASS AR. I’ve added some trace output to list the satellite pairs used to make it easier to follow. RefSats are the reference satellites for each double difference pair, fixSats are the other satellite in each pair. N(0) is the float solution, N(1) is the fixed solution, and ratio is the AR ratio. All three outputs are for the same epoch but in three different runs. I picked a random epoch from the middle of the run.

Integer ambiguity resolution with GPS satellites (17:56:52)

3 resamb_LAMBDA : nx=98
3 ddmat: gps=1/1 glo=0/0
3 refSats= 10 10 10 10 10 10 10
3 fixSats= 13 15 18 20 21 22 29
3 N(0)= -1021543.012 -583342.043 -321520.024 -1128921.004 -351484.032 -1244853.988 -1342905.008
3 N(1)= -1021543.000 -583342.000 -321520.000 -1128921.000 -351484.000 -1244854.000 -1342905.000
3 N(2)= -1021543.000 -583342.000 -321520.000 -1128921.000 -351484.000 -1244854.000 -1342904.000
3 resamb : validation ok (nb=7 ratio=160.92 s=132.06/21251.08)

 

Integer ambiguity resolution with GLONASS satellites (17:56:52)

3 resamb_LAMBDA : nx=98
3 ddmat: gps=0/0 glo=1/1
3 refSats= 33 33 33 33 33
3 fixSats= 39 40 54 55 56
3 N(0)= -1572711.936 -785246.994 -1482008.966 -496531.985 151149.024
3 N(1)= -1572712.000 -785247.000 -1482009.000 -496532.000 151149.000
3 N(2)= -1572711.000 -785247.000 -1482008.000 -496532.000 151149.000
3 resamb : validation ok (nb=5 ratio=165.97 s=38.82/6443.06)

 

Integer ambiguity resolution with GPS and GLONASS satellites (17:56:52)

3 resamb_LAMBDA : nx=98
3 ddmat: gps=1/1 glo=1/1
3 refSats= 10 10 10 10 10 10 10 33 33 33 33 33
3 fixSats= 13 15 18 20 21 22 29 39 40 54 55 56
3 N(0)= -1021543.029 -583342.057 -321520.028 -1128921.020 -351484.040 -1244853.987 -1342905.016 -1572711.928 -785246.982 -1482008.949 -496531.972 151149.020
3 N(1)= -1021543.000 -583342.000 -321520.000 -1128921.000 -351484.000 -1244854.000 -1342905.000 -1572712.000 -785247.000 -1482009.000 -496532.000 151149.000
3 N(2)= -1021543.000 -583342.000 -321520.000 -1128921.000 -351484.000 -1244854.000 -1342905.000 -1572712.000 -785247.000 -1482008.000 -496532.000 151149.000
3 resamb : validation ok (nb=12 ratio=53.40 s=473.13/25263.69)

 

Everything looks as we would expect, so again that is encouraging.

Next I made a comparison using my four standard metrics for four solutions using different AR settings. The first is no AR, the second is GPS AR, the third is GLONASS AR, and the fourth is GPS AR+GLONASS AR. The plot is below, the x-axis is the four cases.

All looks good except that there is a small increase in the std(dist) metric for GLONASS AR which is measuring position error in the solution. This is not surprising and probably would be concerning if it were not so. The code has made a single guess at the inter-channel biases and so we would not expect to get them precisely correct. Any error in our inter-channel bias estimates will translate directly into error in the fixed solutions. If we want to improve the accuracy of the GLONASS fixed solutions we will need to improve the estimates of the inter-channel biases. This should not be difficult if we iterate the estimations over multiple epochs. Even without improved accuracy, we should still get increased robustness from having more satellites to use for the fixed solution. Note the increase from 8 to 13 satellites used on average for the fixed solution in the bottom right plot. Having more satellites available for the fixed solution should improve our ability to stay locked.  To take advantage of this, however, we may need to exclude individual satellites from being used for ambiguity resolution or break the satellites into multiple sets since with the current algorithm a single bad satellite can prevent getting a fix.

Anybody have experience with using this sort of approach to estimate the inter-channel biases?  Any obvious pitfalls I’m missing?  Any other thoughts?

This is still just an experiment, not meant to be a functional feature yet, but I’ve posted the code to my Github repository in the exp1 branch in case anyone else would like to play with these options.  The config file in the exp1\data\demo1 folder is setup to run with these two features both enabled.