In this post I am going to look at some of the raw satellite data directly without using RTKLIB to calculate a solution. Because of the complexity of the full solution, it is easy to start treating RTKLIB as a black box, adjusting input parameters to see how they affect the output without really understanding what is going on inside. It is often possible, however, to simplify the problem down to something we can solve, or at least analyze, with some very simple calculations. Doing this can provide some real insight into what is going on inside the more complex full solution.

In this case, I have a data set from my two M8N receivers, both stationary and within a couple meters of each other with skies partially obscured by trees. There are enough cycle slips in the data to cause RTKLIB some trouble with the solution. The eventual goal here is to improve RTKLIB’s response to cycle slips, but to start with I just want to understand more about what the signals look like during the cycle slips.

Since the two receivers are very close together, double differencing should be very effective at cancelling errors from the ionosphere and troposphere as well as the receiver and satellite clock biases. Since the receivers are stationary, the position component of the double differences will be constant and can be ignored. Double differencing is a very powerful error cancellation technique and is the basis for all differential GPS solutions. I’ll skip a detailed tutorial here, but if you aren’t familiar with how and why it works, a Google search will quickly bring up many explanations better than I could give.

To get the double differences, we first calculate the single difference value for each satellite by subtracting the raw carrier phase measurement from one receiver from the raw carrier phase measurement from the other receiver. Since the satellite clock bias and atmospheric error terms are in both measurements, they will be cancelled or at least minimized by this subtraction.

Next, we subtract the single difference value from one satellite from the single difference value from a second reference satellite. The result is the double difference. This second subtraction will cancel the receiver clock biases, further reducing the errors. Since the receivers are stationary (relative to each other) the position component of the double differences will be constant and so the dominant remaining component will be the phase bias. The phase bias should also be constant unless there is a cycle slip. So, what we expect to see in the double differences is a constant value if there is no cycle slip, and a jump equal to the size of the slip if there is a slip. For this exercise we don’t care about the actual values of the phase biases, only any changes in the values.

The cycle slips are reported by the receiver in the raw data. One thing to be aware of it that the receiver is reporting possible cycle slips, not confirmed cycle slips, so not all the reported slips are real. How many of the reported slips are real is something I hope to be able to tell by looking directly at the double differences.

I wrote a simple matlab script to parse the raw observation files (*.obs). This is quite a bit easier if the observation files are in Rinex 3 format rather than the older Rinex 2 format. This is a configuration option in RTKCONV and CONVBIN to select this. Unfortunately it does not seem to work in CONVBIN so I had to use RTKCONV to create the observation files. Here is a snippet of the Rinex 3 observation file for a typical cycle slip on a single satellite.

The top line is the time stamp, each line below is for one satellite. In this example, the first column is the satellites number, the third column is the carrier phase. A cycle slip shows up as a fourth digit after the decimal point in the carrier phase measurement. In this example there is a cycle slip for satellite R8.

Here is another example, but in this case, the cycle slip is reported for all satellites.

This is very difficult for RTKLIB to deal with since the cycle slips, whether they are real or not, immediately cause it to reset the phase bias estimates. In this case, RTKLIB resets all the bias bias estimates and has to start over completely. We’ll use this example to look more closely at what is going on inside the solution.

Let’s look at the G13 satellite first. We’ll use G12 as the reference satellite since that is what RTKLIB used. Here are the carrier phase measurements for both satellites from the two receivers, pulled from the observation file.

Note the y-axis has a scale of 10^6. Even though the receivers are not moving relative to the earth, the distance between the satellites and the receivers is changing quite quickly. This is because the satellites are moving and because the earth is rotating.

Below is the single difference measurement for G13, calculated by subtracting the rover measurement (yellow line above) from the base measurement (purple line above). The mean is also subtracted to center the plot around zero. The red circles are where cycle slips occurred.

The cancellations that occur from the subtraction have reduced the y-axis scale down from 10^6 to just tens of cycles but the remaining variation is still too large to see a cycle slip. Again, this calculation is done directly on the raw data from the observation file, RTKLIB is not used at all to process the data.

Next we get the double difference by subtracting the single difference measurement for the reference satellite (G12 in this case) from the single difference for G13 plotted above. Again the mean of the result is removed to center the result around zero. This results in the plot on the left below. Again, the red circles are where cycle slips occurred.

We have now cancelled some additional errors and the variation on the y-axis is very small, in fact small enough to see a cycle slip. It is quite encouraging how this very simple mathematical operation on the raw data can give us such meaningful results.

The above plot on the right comes from data in the RTKLIB state output file and shows the kalman filter estimate of the phase-bias for satellite G13. Again the mean is removed so only the variation is relevant. The large excursions occur when the kalman filter states are reset because of the reported cycle slips. The filter states are re-initialized using the pseudorange measurements which are much less accurate than the carrier-phase measurements but give an absolute value rather than just a relative number. As you can see, in this case it takes several minutes for RTKLIB to fully recover the phase bias estimate but when it does, it converges back to the same value as before the reported slip.

Looking back at the plot on the left, we can see for the first cycle slip there is no jump in the measurement. If there was a cycle slip we would expect a jump in this measurement equal to the size of the cycle slip. For the second cycle slip there is a brief half cycle discontinuity which I don’t understand, but since it quickly returns to zero we’ll ignore it, at least for the time being.

What our very simple calculation is telling us is in this case is that neither potential cycle slip reported by the receiver was actually real and RTKLIB would have done much better to not have reset the phase bias estimates. This is also confirmed by RTKLIB since it eventually converges back to the same phase-bias estimate as it had previous to the potential cycle slip.

Let’s look at another satellite from the same reported cycle slip. On the left below is the double difference for G18 calculated directly from the raw data just as we did for G13 above. On the right again is the kalman filter estimate of the phase bias from RTKLIB and in both plots, the reported cycle slips occur at the red circles.

This time the first cycle slip is real and you can see the two cycle jump in the double difference value on the left. On the right you can see that RTKLIB eventually settles on the same two cycle slip, but takes several minutes to figure it out. The second cycle slip again causes a brief half cycle deviation in the double difference but then returns to zero, and again we will ignore it.

Based on the plots on the left which are generated with only three basic subtractions using the raw carrier phase data, it appears that RTKLIB could be doing a significantly better job of handling cycle slips. Not only in determining if potential reported cycle slips are real or not, but also in correcting the cycle slip more quickly when it is real. It is true that the baseline in this example is very short and with longer baselines the errors in the double differences will be greater, making things more difficult, but I still have to believe there is opportunity for improvement here.

In the example above, of the thirteen satellites processed by RTKLIB, all of which were reported in the raw data as having cycle slips, only four of them actually did. Even if RTKLIB reset the phase-bias estimates for the four satellites with slips, it would still handle the slip much better using the other nine satellites than it does now where it resets all thirteen phase-bias estimates. If it could correct the slip as well as correctly identify it, the robustness to slips would be even greater.

I am actually quite surprised how clear cut the cycle slip examples I looked at are. I expect this is not always the case or cycle slip correction would probably have been added to RTKLIB a long time ago. Still, I hope to use this approach to try and add some amount of cycle slip detection and correction to RTKLIB, even if it can only be applied in some subset of cases.

Anyone else have experience trying something like this?

Can i get the matlab script to parse the raw observation files (*.obs) ?

LikeLike

Hi Hakim. I have found in the past that sharing matlab scripts generates more requests for support than I can handle, especially given that mine are a bit of a hacked together mess. I would prefer to stay focused on RTKLIB and leave this for someone else. Writing your own is not a difficult task and having the ability to parse text files, whether it is with Matlab or Python or some other tool, is a skill that you will find useful over and over again.

LikeLike

Hello,

Thanks for the site!I find this very useful!

I have a question…is there a way to use RTKLIB to save a file with the cycle slips time occurrence?

It seems to me that I can only visualize them..am I wrong?

Many thanks!

LikeLike

Hi Susan. Cycle slips are listed in the solution status file (*.pos.stat) under the residuals records ($SAT) with time stamps. This is an all text file and fairly easily parsed with Python or Matlab.

LikeLike

many thanks!

LikeLike

Hi I have a quick question out of curiosity. If there were two receivers both using network RTK positioning techniques with 1m between them (both are rovers), and both of them experienced a cycle slip, would it be unlikely that they experience the same compromise in accuracy (both in direction and quantity)? I.e. would it be unlikely that both experience a cycle slip and when analysing the “measured distance” between them for that epoch, it would still be 1m?

LikeLike

Hi Omair. Cycle slips are normally flagged by the receiver and cause RTKLIB to re-initialize the phase-bias estimate for that satellite based on a single sample. The re-initialized phase-bias values for the two receivers are very unlikely to be the same since they are coming from independent sets of electronics and are somewhat noisy. However, any effect on the RTK solution will usually be non-existent or very small since the solution will be derived primarily from the other satellites that did not have a slip. So, in that sense, the difference between the two solutions are likely to remain at one meter (and both be correct). If the slip does cause an error, it is unlikely that the error would be the same in both solutions. Not sure if that answers your question or not?

LikeLike

Let me begin by saying that for the purposes of visually looking for cycle slips, I think your method is fine. If you are planning to use this information to do slip detection, you might need to address the timing detail, as my experience is the same as Michele’s.

Since your receivers aren’t sampling at the same time, the receiver sampling some dt later than the other will accumulate more phase (Doppler*dt). Even though the dt is the same for all of that receiver’s satellite measurements, the Doppler is different, and so this difference doesn’t cancel out. With the zero-baseline case (single antenna, 2-way RF splitter), the DD phase residuals should be integers, so it becomes clear that something is not canceling out (even more so on long data sets).

In your case, even though the Doppler is relatively large (~2kHz on G13), the error term per-satellite (Doppler*dt) may be nearly constant since:

a) with two TCXO-based receivers, the sampling time offset (dt) is relatively stable

b) you are looking at a relatively short time span, and the satellite velocity along your line of sight is not changing greatly over that interval (ie, low acceleration along your line of sight)

A nearly constant error term per satellite will disappear when you subtract off the mean. On top of this bulk error you are eliminating, your drift of half a cycle per hour is in line with what I would expect with an expected change over time of the satellite Doppler measurement and sampling offset.

One thing to watch for is when one receiver makes a 1ms clock adjustment. If the clock adjustment happens during your data collection, it will look a bit like a cycle slip with your current method.

Note also that since your antennas have some separation, you should not expect the DD phase residuals to remain constant as the satellite geometry changes. However, with such a short baseline, this also will only cause slow changes in the phase differences.

I enjoy reading your posts here, and look forward to your future investigations.

LikeLike

Yes the receiver clock bias goes away in inter-satellite differences (your second difference).

However, measurements taken as different times (as receivers have free running clocks) will necessarily differ. So your first difference may contain errors due to not propagating the base observables from the base time to the rover time (or vice-versa).

At least, this is what I need to do in zero baseline to achieve integer cycle DDs.. perhaps the operation of removing the mean twice leaves only the phase trend and that appears to be enough to analyse cycle slips in the short term.

It would be good to share a couple of datasets to see if we are on the same wavelength 🙂

LikeLike

I have done very similar analysis in zero baseline. Didn’t you have to interpolate measurements to a common time before differencing? I could not get anything meaningful otherwise, since the speed of the satellite is not negligible compared to a wavelenght in presence of typical clock biases.

The half-cycle ambiguity sounds like the receiver waiting to receive the last two bits of the GPS subframe?

LikeLike

Hi Michele. I did not do any interpolation of the raw data, just aligned it so the receiver time stamps matched. If I understand your concern correctly, it is that the two data streams will still be misaligned by the sum of the clock biases. This should be on the order of 1 usec. To see how this much error would affect my plots I interpolated the data to shift one receiver’s data by 1 usec and re-plotted the data. I do see a small shift in slope of the double difference but only about a half cycle per hour, not nearly enough to obscure the cycle slips. I am only looking for relative changes over short time intervals, is it possible you were trying to do more than that in your analysis?

LikeLike

A clock bias of 1 us is very small: perhaps you meant 1 ms? With a TCXO at 0.5 ppm accuracy (as in M8N), the receiver clock slips 1 us in a couple of seconds only which would lead to an awful lot of steering. If the Doppler shift on the carrier is say 2 kHz the phase will cover a couple of circles in 1ms which is enough to invalidate any DD analysis you may be doing.

Maybe RTKLIB allows for adding the clock bias term to your Rinex files at conversion time, which would ease the process of interpolation in case you change your mind.

References:

https://forum.u-blox.com/index.php?qa=1529&qa_1=is-the-receiver-local-time-steered-to-gps-time

LikeLike

Hi Michele. Sorry, you are right, I did use a clock bias of 1 msec, not 1 usec in the data shift experiment I described in the last comment.

I tried the same DD analysis on another M8N data set and two M8T data sets to see if I just got lucky somehow the first time, but I got clean results on all four data sets so it does look like what I am doing works. I also verified that the receiver estimates of clock bias are not getting added to the raw carrier phase measurements during the conversion process.

I think I am missing something in your argument. It seems like you are telling me that I need to account for the receiver clock bias terms before I do the double difference, but my understanding is that the double difference operation should eliminate the clock bias terms through cancellation. Is that not right?

LikeLike