GPS position precision

[joelmatejka]

TLDR: gps-sdr-sim -> HackRF -> ublox GPS: drifts for dozens of meters in static and also dynamic mode, sometimes looses fix

Setup

I tried multiple configurations. But always similar scenario:

• Schematic

  ┌─────────────┐    ┌─────────────────┐    ┌────────────┐  sma cable  ┌─────────┐ ┌─────────┐ ┌────────┐  ┌────────────────────┐    ┌───────────────┐
  │ gps-sdr-sim ├─┬─►│ hackrf_transfer ├───▲│ hackrf one ├────────────►│ 30db at │►│ 30db at │►│ DC blk ├─►│ ublox neo m9v/m8n  ├─┬─►│ ublox ucenter │
  └─────────────┘ │  └─────────────────┘   ◄└─────────▲──┘             └─────────┘ └─────────┘ └────────┘  └────────────────────┘ │  └───────────────┘
                  │                        │ ▲        │                                                                           │
                  │  ┌────────────────┐    │ │    ┌───┴──┐                                                                        │  ┌────────────┐  ┌────────┐
                  └─►│ hackrf android ├────┘ │    │ txc0 │                                                                        └─►│ linux gpsd ├─►│ marble │
                     └────────────────┘      │    └──────┘                                                                           └────────────┘  └────────┘
                                             │
  ┌───────────────────┐                      │
  │ multi-sdr-gps-sim ├──────────────────────┘
  └───────────────────┘
  

• HackRF one rev. 9:

  $ hackrf_info
  hackrf_info version: 2023.01.1
  libhackrf version: 2023.01.1 (0.8)
  Found HackRF
  Index: 0
  Serial number: 0000000000000000675c62dc326875cf
  Board ID Number: 4 (HackRF One)
  Firmware Version: 2023.01.1 (API:1.07)
  Part ID Number: 0xa000cb3c 0x00614353
  Hardware Revision: r9
  Hardware appears to have been manufactured by Great Scott Gadgets.
  Hardware supported by installed firmware:
      HackRF One
  $ hackrf_clock -i
  CLKIN status: clock signal detected
  

• TX gains - I tried multiple so that I see satellites with SNR around 30 to 55 dB. And tried to skip one 30 dB attenuator.

• I have TXC0 installed and HackRF should know about that and hopefully use it.

• I also tried multiple sample rates e.g.

  ./gps-sdr-sim -b 8 -e brdc2480.23n -s 13000000 -l 39.01973077309149,125.75254486374945,42.0 -o pyongyang.bin
  

  hackrf_transfer -t pyongyang.bin -f 1575420000 -s 13000000 -a 0 -x 0 -R
  

• I also tried to play with ppm correction, but when I try more than +-1 ppm, I'm not getting fix at all or loosing more often.

And now, I'm quite out of ideas what else to try.

Expected behavior

Get fix in a few minutes and simulate precise location (+- a few meters at maximum).

Issues

In any configuration I tried, I'm able to get 3D fix. However, the position is always drifting quite a lot (see pic). Usually after cold start, I get right position, which holds for a few minutes. But then it starts to walk around expected position in both static and dynamic modes. Sometimes, I also lose fix. Altitude is almost always wrong. When I repeat the same simulation, the drift is always different.

Here, scale is missing, but the square in the middle is 100m × 100m.

Questions

1. What is a precision you usually can achieve with HackRF or other devices available?
2. I'm not expert in GPS nor RF, but... Can this be caused by bad clock source (it is probably some Chinese xtal on that TXC0 board even though it was declared that it has 0.1 ppm) or some precision losses in simulation?
3. Could it be HackRF hardware problem?
4. Have you any other idea what to try?

Thank you

[KingofTown]

Honestly, with the base code, I was having similar issues as you. I've made a TON of modifications to my code to work better (but also rewritten it in c++, made it full realtime and added a full gRPC interface/web interface). Because of my mods, it's not really possible to do a pull request.

The issue was unrelated to clock source for me. Even my $800 DS Instruments external TCXO wouldn't keep a steady lock on satellites. In UBLOX U-Center, I would watch the residual error drop and drop until it resets. (NOTE - I'm talking specifically about an external clock. The one that comes with the HackRF is NOT good enough, and no amount of software changes will help that one. I'm saying even after getting a good clock, I've modified the software a bit).

With my changes, I'm able to run a stable fix at even 3MHz sample rate with a HackRF. Residual errors all remain ~= 0.5, and my PDOP ~= 1.7 and HDOP ~=1.0.

And when I say I've done "a TON of work", most of that work as learning how GPS works. This code base is still amazing. After months of digging, I've narrowed it down to a few things that could be added for better long-term tracking.

1. A simple troposphere model - nothing crazy needed. I tried with a Saastamoinen model which used relative humidity, but that's overkill if you have no way to input humidity. Instead I switched to one that just uses altitude and az/el angles to compute tropo delays.

2. I compute Range Rate from ephemeris. The current code is computing the range rate by computing the position vectors every 0.1 seconds and calculating velocity vectors from the differences of the last measurement to the current measurement. This works OK, but it is not incorporating the satellite 1st and 2nd order clock errors at all (af1, af2). Meaning, the further you get from the TOE, the further off your rates will be as the clock errors wont be accounted for (but are still sent via LNAV messages). Perhaps settings those values in the LNAV to 0.0 would make the 2 computations more similar though.

3. With my range rate enhancements, I'm also able to correct for the phenomenon the ionosphere has on the carrier phase. Carrier phase is sped up by the ionosphere while the code phase is slowed. For example, in my code I have the following:

	// Start with the inital rate (this includes satellite clock errors and velocity)
	double ini_rate = -rho1.rate/gpssim::LAMBDA_L1;

	// Determine the rate of change from iono and tropo from the last sample
	// These are in units of meters, but we're interested in the difference
	// to get the doppler change.
	double iono_rate = (rho1.iono_delay - chan->rho0.iono_delay)/dt;
	double tropo_rate = (rho1.tropo_delay - chan->rho0.tropo_delay)/dt;

	// Carrier phase is sped up by ionosphere and slowed by troposphere.
	// (subtraction is speed up, addition is slow down)
	chan->f_carr = ini_rate - iono_rate + tropo_rate;

	// Code phase is slowed down by both equally
	chan->f_code = gpssim::CODE_FREQ + (ini_rate + iono_rate + tropo_rate)*gpssim::CARR_TO_CODE;

Notice how instead of computing a rate using position, I'm instead using the computed one from the range_t (which includes the clock errors).

4. The BIGGEST contributor to Residual error I can't fully explain. I've added an entire calibration section to my code to account for things like clock bias, clock drift, RF Cable lengths, dielectric constants, etc. Looking through old comments, came across a few from Ivan about needing to add an offset (Doppler shift, need to add a "blue" offset  #192). While I can't say I agree with why he says these changes are needed, I can say I also have to add a clock drift offset to my range rate. Currently, I have a -13.95ns clock drift that gets added to ALL range rates (~= -4 m/s).

	// Compute the initial range rate
	double rangeRate = vs.dotp(los)/range;
	rho.rate = rangeRate;

	// Satellite clock errors (clk[1] s/s) reduces rate
	// Treat sim clock errors as if they are satellite clock errors
	double rate_clock_error = calibration_.sim_clock_drift_ns.value_or(0.0) / 1000000000.0 + clk[1];

	//               m/s               *     s/s    => m/s
	// Doppler is negatively affected by clock errors
	rho.rate -= gpssim::SPEED_OF_LIGHT * rate_clock_error;

My theory (which is likely very wrong) is that there are 2 parts to clocks - their stability and their aging rate. Stability is the most important part in making sure that all samples get transmitted at a stable rate. They will state that their clocks will be within 0.1 ppm over all temperature ranges which is great. However, we have to account for the clock age rate in the sim.

https://www.dsinstruments.com/product/reference-sources/ultra-stable-frequency-reference-source-10mhz-100mhz/

For example - this clock I was using states it has a daily aging rate of +/- 1 PPB, and a yearly rate of +/- 50 PPB. Note - the aging isn't a linear line....it tends to be the highest in the first few years of the clock. Thus, my correction of -13.95 ns (or about -14PPB) doesn't seem that far fetched for me. The only reason I don't fully buy into my own logic is that this same offset works well for ALL of my clocks I have.

---

As for how this can help you now, I'd try changing your ephemeris clock errors to all be 0.0 and adding an offset similar to how it was done in issue 192. I think I'd make this offset be a command line input since I think this value depends on hardware. (i.e. Ivan subtracted 0.84 m/s where I had to subtract nearly 4 m/s. The easiest way for me to know what my offset should be was to watch the residual errors in Ublox (UBX - NAV - SAT) until they stopped drifting. I think this single change alone would help you the most....but you'll have to play around to see what the actual value should be.

There is a small bug in the ionosphere calculation lines 1198 and 1206 (https://github.com/osqzss/gps-sdr-sim/blob/master/gpssim.c#L1198), its mixing radians and semicircles in calculations. Line 1198 should divide azel[0] by pi and 1206 should also divide by pi (then not divide phi_i by pi). Basically, either make them ALL semicircles, or convert the values to all Radians (i.e. the constants of 0.0137/(E + 0.11) - 0.022 could instead be in radians). But honestly, because the current code only applies ionosphere to the pseudorange, fixing this won't change a whole lot in terms of drifting.

[cjb-eng]

FWIW, I did a fair bit of work exhaustively comparing the group/phase pseudorange observables to the same observables computed using the open source GNSS toolkit GPStk. The GPStk calculation are in agreement with truth rinex data reported by commercial simulators at the mm-level. What I discovered in this analysis is that the calculation of light time by gps-sdr-sim is a simplification which is understandable if you're trying to compute quantities in real-time or perform the calculations quickly. For stationary receivers grounded to the Earth, the direct form of the light time calculation implemented by gps-sdr-sim is a reasonable approximation. But if you're simulating a receiver in an accelerated frame of reference (drone or worst case LEO satellite) the range rate accuracy is going to deteriorate. I suspect this is at least part of the range rate issue you talked about above. Generally, light time has to be computed in a higher order fashion (e.g. iteratively) rather than the direct manner applied by gps-sdr-sim.

[joelmatejka]

@KingofTown Thank you for the comprehensive response. Do you have your rewritten code available somewhere, or is it strictly proprietary now? :) I just soldered rf shield and I'm getting repeating patterns on the same data executed multiple times, so it points to the simulator code... Strange is that I'm getting the right position with a random android phone in non-conductive mode... But maybe u-blox gps receivers are too wily and do some too smart things that ... Were you playing with u-blox device as well?

[KingofTown]

@joelmatejka

I've tested with a Gen 8 ublox, a Gen 9 ublox and various other receivers I could get my hands on. UBlox was nice in that it was easy to see Residual errors.

I would just add a value into this line in the code:

https://github.com/osqzss/gps-sdr-sim/blob/master/gpssim.c#L1324

Then re-run to see if it's more stable. Maybe try starting by adding 4.5 to that value and see how it affects the output.

I don't think I can upload my code as there is a lot of other stuff that wouldn't be allowed, but I can show my functions.

My current calibration settings are:

RX Clock Drift (ns), Rx Clock Bias (ns), Sim Clock Dift (ns), Sim Clock Bias (ns), Doppler Offset (m/s), RF Cable Len (cm), RF Cable Dielectric                                                                                                     
742                           ,0                ,-13.95             ,0,                     0,                    100,                     0.87 

Then I've modified the computeRange function to instead calculate the range-rate and store it in range_t so that the first order clock errors get included. In addition, I add my own clock errors to the range-rate here as well.

range_t
GPSSim::computeRange2(channel_t *chan, xyz_pos& xyz, gpstime_t* time)
{
	// Can comute a range for a custom time instead of using receiver time.
	if (!time) {
		time = &grx_;
	}

	int sv = chan->prn - 1;

	range_t rho;

	double pos[3],vel[3],clk[2];
	double clk[2];
	double los[3];
	double tau;
	double range,rate;
	double xrot,yrot;

	double llh[3],neu[3];
	double tmat[3][3];

	// Calculate SV position at time of the observation (grx_).
	satpos(eph_[ieph_][sv], *time, rho.pos, rho.vel, clk);

	// Subtract receiver position from satellite position to get the line of sight vector.
	subVect(los, rho.pos, xyz);

	// Calculate magnitude of range vector, then divide by speed of light
	// to get the signal travel time
	tau = normVect(los)/gpssim::SPEED_OF_LIGHT;

	// Extrapolate the satellite position backwards to the transmission time.
	// Basically, take the calculated velocity vectors and move a bit further back
	// by the signal travel time
	rho.pos[0] -= rho.vel[0]*tau;
	rho.pos[1] -= rho.vel[1]*tau;
	rho.pos[2] -= rho.vel[2]*tau;

	// Earth rotation correction. When the signal was sent, the earth rotated by
	// this much. Add this rotation distance to further correct for satellite position.
	xrot = rho.pos[0] + rho.pos[1]*OMEGA_EARTH*tau;
	yrot = rho.pos[1] - rho.pos[0]*OMEGA_EARTH*tau;

	rho.pos[0] = xrot;
	rho.pos[1] = yrot;

	// Save receiver location for this measurement
	rho.xpos[0] = xyz[0];
	rho.xpos[1] = xyz[1];
	rho.xpos[2] = xyz[2];

	// Given this updated satellite position, again difference with receiver to
	// get the line of sight vector.
	subVect(los, rho.pos, xyz);

	// range (in meters)
	range = normVect(los);
	rho.d = range;

	// Pseudorange is affected by clock errors.
	// Subtract satellite clock offset from GNSS time (in seconds) (clk[0])
	//    - sim clock bias is also subtracted as it's part of the satellite clock offset
	double clock_error = calibration_.sim_clock_bias_ns.value_or(0.0)/1000000000.0 + clk[0];

	// Pseudorange takes into consideration clock errors.
	rho.range = range - gpssim::SPEED_OF_LIGHT*clock_error;

	// This value is not used in any calculations. It is simply to display the
	// receivers computed range accounting for their current clock bias. It's only
	// used for validating as it would require input to get the current receiver clock
	// bias at this given time t.
	rho.rx_corr_range = rho.range + gpssim::SPEED_OF_LIGHT*(calibration_.receiver_clock_bias_ns.value_or(0.0)/1000000000.0);

	// Save this value for TATT calculations into the range_t struct
	// (the functions that use range_t might not have access to class variables)
	rho.rf_cable_delay_m = current_cable_delay_m_;

	// Relative velocity of SV and receiver.
	rate = dotProd(rho.vel, los)/range;

	// Satellite clock errors (clk[1] s/s) reduces rate
	// Treat sim clock errors as if they are satellite clock errors
	double rate_clock_error = calibration_.sim_clock_drift_ns.value_or(0.0) / 1000000000.0 + clk[1];

	//               m/s               *     s/s    => m/s
	// Doppler is negatively affected by clock errors
	rho.rate = rate - gpssim::SPEED_OF_LIGHT * rate_clock_error;

	// Calibration offset is the same as sim clock drift, just in different units.
	// i.e. a clock drift of 10 ns/s is the same as a doppler offset of ~ 3.3m/s
	rho.rate += calibration_.doppler_offset_meters_per_second.value_or(0.0);

	// Time of application. This is the time used to calculate all above values.
	rho.g = *time;

	// Azimuth and elevation angles.
	xyz2llh(xyz, llh);
	ltcmat(llh, tmat);
	ecef2neu(los, tmat, neu);
	neu2azel(rho.azel, neu);

	if(performanceOptions_.compute_iono_delay){
		// Add ionospheric delay
		// TODO: Does this need to affect the rho.rate as well?
		rho.iono_delay = ionosphericDelay(&ionoutc_, *time, llh, rho.azel);
		rho.range += rho.iono_delay;
	}

	if(performanceOptions_.compute_tropo_delay){
		rho.tropo_delay = tropmodel2(llh, rho.azel);
		rho.range += rho.tropo_delay;
	}

	return rho;
}

Then in computeCodePhase() - you just use use the computed range-rate instead of doing the math

void computeCodePhase(channel_t *chan, range_t rho1, double dt, double sampT)
{
	(void)sampT;
	// C/A code is slowed by the ionosphere (group delay), but the carrier wave itself appreas to
	// speed up in the ionosphere (phase delay).
	// Both are slowed equally from tropo delay.

	// Start with the inital rate (this includes satellite clock errors and velocity)
	double ini_rate = -rho1.rate/gpssim::LAMBDA_L1;

	// Determine the rate of change from iono and tropo from the last sample
	// These are in units of meters, but we're interested in the difference
	// to get the doppler change.
	double iono_rate = (rho1.iono_delay - chan->rho0.iono_delay)/dt;
	double tropo_rate = (rho1.tropo_delay - chan->rho0.tropo_delay)/dt;

	// Carrier phase is sped up by ionosphere and slowed by troposphere.
	// (subtraction is speed up, addition is slow down)
	chan->f_carr = ini_rate - iono_rate + tropo_rate;

	// Code phase is slowed down by both equally
	chan->f_code = gpssim::CODE_FREQ + (ini_rate + iono_rate + tropo_rate)*gpssim::CARR_TO_CODE;

	// Calculate the 'time at the tone' in which this phase measurement will
	// releate to (in seconds).
	// chan->g0 is put in generateNavMsg and is the time corresponding to when the
	// message was generated.
	// Add 6 seconds since the HOW corresponds to the time of signal transmission
	// of the NEXT subframe (each subframe is 6 seconds long)
	auto tatt = subGpsTime(chan->rho0.g, chan->g0) + 6.0;

	// Time it would take this message to reach the receiver
	// rho.range includes iono and tropo delays already.
	// Also, correct for the RF cable length so it arrives exactly when it should
	auto ctime = (chan->rho0.range - chan->rho0.rf_cable_delay_m)/gpssim::SPEED_OF_LIGHT;

I added a few more fields to range_t to store the tropo and iono delay, but I don't think those have a huge effect outside of accuracy. The biggest one is not incorporating clk[1] into the range-rate, and also not correcting for the sim clock drift in range-rate as these affect the perceived speed of light causing receivers to drift away.

All of my other changes make it so the computed position on receivers matches to at least 5 decimal places in lat/lon or higher. Now, the biggest source of error is due to WGS-84 model (sim and receiver are both estimating the earth shape which compounds the errors). At the equator, my output is much more stable compared to near the poles.