In GNSS, when considering the propagation of signals from the satellites to a receiver, it is easier to work in an ECI reference frame, since (ignoring the gravitational potential of Earth), light travels in straight lines in ECI coordinates. However, it is often common to do all the calculations in an ECEF frame, as the final goal is to obtain the receiver’s position in ECEF coordinates, and the ephemerides also use ECEF coordinates to describe the satellite positions. Therefore, a non-relativistic correction needs to be applied to account for the fact that light no longer travels in straight lines when one considers ECEF coordinates. Often, the correction is done as some kind of approximation. These types of corrections are known in the GNSS literature as the Sagnac effect.
The goal of this post is to discuss where the corrections arise from, the typical approximations that can be made, and how these corrections affects the calculation of range and range-rate. I didn’t find a good source in the literature where this is described in detail and in a self-contained way, so I decided to write it myself.
Range Earth rotation correction
Before beginning, let us fix the notation. We consider an ECEF frame and an ECI frame so that the rotation matrix that takes ECEF coordinates into ECI coordinates at time \(t\) is given by \(e^{A\omega_Et}\), where \(\omega_E = 7.292115\cdot10^{-5}\,\mathrm{rad}/\mathrm{s}\) denotes Earth’s angular rate and \[A = \begin{pmatrix}0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{pmatrix}.\]Note that this ignores all considerations about the irregularities of the rotation of Earth in space. We have also arranged the origin of the time coordinate conveniently so that at \(t = 0\) the axes of the ECEF frame coincide with those of the ECI frame.
We fix some time \(t\) and denote by \(x_R\) the position of the receiver in ECEF coordinates at time \(t\). Given a satellite, the range \(\rho\) of that satellite is defined by the fact that the time that the signal of the satellite takes to arrive to the receiver is \(\tau = \rho/c\). Now, this implies that the signal was transmitted by the satellite at time \(t – \tau\), so we denote by \(x_S\) the position of the satellite in ECEF coordinates at time \(t – \tau\).
The fact that light doesn’t travel in straight lines in ECEF coordinates implies that, in general, \(\rho\) is not equal to \(\|x_R – x_S\|\). To compute \(\rho\) we need to convert to ECI coordinates. The position of the receiver at time \(t\) in ECI coordinates is \(e^{A\omega_E t}x_R\), and the position of the satellite at time \(t-\tau\) in ECI coordinates is \(e^{A\omega_E(t-\tau)}x_S\). Since light travels in straight lines in ECI coordinates, we have\[\tag{1}\rho = \|e^{A\omega_Et}x_R – e^{A\omega_E(t-\tau)}x_S\|.\]This equation is the basis of all the correction formulas and can be found as equation (E.3.8) in RTKLIB’s manual.
We put \(\sigma = \omega_E\tau\). Since \(e^{A\omega_E t}\) is an isometry, we see that\[\tag{2}\rho = \|x_R – x_S – (e^{-A\sigma} – I)x_S\|.\]We want to do an approximation so that we can write the right hand side as \(\|x_R – x_S\|\) plus some correction. There are several ways in which this could be done, but the main idea is to do an approximation of first order in the parameter \(\sigma\), which is small. Indeed, \(\tau\) is on the order of 80 ms for MEO satellites and receivers near the surface of the Earth, so \(\sigma\) is on the order of \(6\cdot 10^{-6}\). Since \(x_R\) and \(x_S\) are on the order of 6000km and 26000km respectively, terms on the order of \(\sigma x_R\) and \(\sigma x_S\) are still metric level, which justifies the need for this correction. Terms on the order of \(\sigma^2 x_R\) and \(\sigma^2 x_S\) are on the order of 1mm or less, and can be safely ignored in most applications.
We can use the first order approximation of the exponential, \[e^{-A\sigma} – I \approx -\sigma A,\] to get\[\tag{3}\rho \approx \left\|x_R – x_S + \sigma Ax_s\right\|.\]
For vectors \(u,v\) on a real Hilbert space and a small parameter \(s\), it holds that\[\|u + s v\| = \|u\| + \frac{\langle u, v\rangle}{\|u\|}s + O(s^2).\]Applying this to (3) we get\[\tag{4}\rho \approx \|x_R – x_S\| + \frac{\sigma}{\|x_R-x_S\|} \langle x_R – x_S, Ax_S\rangle.\]It is also possible (and perhaps formally more correct) to get directly from (2) to (4) by computing the derivative of the right hand side of (2) with respect to \(\sigma\).
We know from (2) that \(\|x_R – x_S\| = \rho + O(\sigma)\), so\[\frac{\sigma}{\|x_R – x_S\|} = \frac{\sigma}{\rho} + O(\sigma^2),\] and, for an approximation of first order in \(\sigma\) it is acceptable to replace \(\sigma/\|x_R – x_S\|\) by \(\sigma/\rho = \omega_E/c\) in (4). Additionally, since \(A^T = -A\), we have \(\langle x_S, A x_S \rangle = 0\). Thus, we arrive to\[\tag{5}\rho \approx \|x_R – x_S\| + \frac{\omega_E}{c} \langle x_R, A x_S \rangle.\]This is the usual form of the Earth rotation correction for the range, and indeed can be found as (E.3.8b) in RTKLIB’s manual. To see this, write \(x_R = (x^R, y^R, z^R)^T\), \(x_S = (x^S, y^S, z^S)^T\), so that the correction is\[\tag{6}\frac{\omega_E}{c}\langle x_R, A x_S \rangle = \frac{\omega_E}{c}(y^Rx^S – x^Ry^S).\]
Range-rate Earth rotation correction
For the range-rate \(\rho’\), which is the time derivative of \(\rho\), there are several ways to obtain an approximate Earth rotation correction. The most straightforward is to compute the time derivative of (5), obtaining\[\tag{6}\rho’ \approx \langle x_R’ – x_S’, E \rangle+ \frac{\omega_E}{c} \left( \langle x’_R, A x_S \rangle + \langle x_R, A x’_S\rangle\right),\]where \(E\) denotes the line of sight vector\[E = \frac{x_R-x_S}{\|x_R – x_S\|}.\]Note that this vector is not the direction in which the signal arrives to the receiver (more on this later).
The interpretation of the first summand in the right hand side of (6) is as the projection of the relative velocity vector (in ECEF coordinates) onto \(E\). The second summand in the right hand side of (6) is the approximate Earth rotation correction. Written in coordinates it is\[\tag{7}\frac{\omega_E}{c}(y’^Rx^S – x’^Ry^S + y^Rx’^S – x^Ry’^S).\]This expression appears in (F.6.29) (note that this is a typo and should be (E.6.29)) in RTKLIB’s manual, but it seems to me that the correction used by RTKLIB has the wrong sign!
There is an alternative way to obtain (6) which perhaps makes it more clear that (6) is just an approximation of \(\rho’\) of first order in \(\sigma\). In this manner, we start again with (1) and define the vector \(N\) by\[N = \frac{x_R – e^{-A\sigma}x_S}{\|x_R – e^{-A\sigma}x_S\|}.\]Note that \(N\) is indeed a unit vector in the direction in which the signal arrives to the receiver in ECEF coordinates, since \(e^{A\omega_E t}N\) is a unit vector in the direction of the straight line in which the signal travels in ECI coordinates.
We compute the derivative of (1) as\[\begin{split}\rho’ &= \langle e^{A\omega_E t} x’_R – e^{A\omega_E(t-\tau)}x’_S, e^{A\omega_E t} N\rangle \\ &\quad + \omega_E \langle A (e^{A\omega_E t} x_R – (1-\tau’)e^{A\omega_E(t-\tau)} x_S), e^{A\omega_E t}N\rangle\\ &= \langle e^{A\omega_E t} x’_R – e^{A\omega_E (t-\tau)}x’_S, e^{A\omega_E t} N\rangle \\ &\quad + \omega_E \langle A (e^{A\omega_E t} x_R – e^{A\omega_E (t-\tau)} x_S), e^{A\omega_E t}N\rangle\ \\ &\quad + \omega_E\tau’\langle Ae^{A\omega_E(t-\tau)}x_s, e^{A\omega_E t}N \rangle.\end{split} \]Now the second summand in the last equality is zero because\[e^{A\omega_E t} x_R – e^{A\omega_E (t-\tau)} x_S = \|x_R – e^{-A\sigma} x_S\|e^{A\omega_E t} N\]and \(\langle A u, u\rangle = 0\) for all vectors \(u\). Since \(e^{A\omega_E t}\) is an isometry, we get\[\tag{8}\rho’ = \langle x’_R – e^{-A\sigma} x’_S, N\rangle + \tau’\omega_E \langle A e^{-A\sigma}x_S, N\rangle,\]which is an exact expression for the range-rate \(\rho’\).
The first summand in the right hand side of (8) can be split as\[\tag{9}\langle x’_R – e^{-A\sigma} x’_S, N\rangle = \langle x’_R – x’_S, N\rangle – \langle (e^{-A\sigma}-I)x’_S, N\rangle.\]Since we want to replace \(N\) by \(E\) in these expressions, let us compute an approximate formula for \(N-E\). First note that \(E\) can be obtained by setting \(\sigma = 0\) in the expression for \(N\). Therefore,\[N – E = \sigma \left.\frac{d}{d\sigma}\right|_{\sigma=0}\frac{x_R – e^{-A\sigma}x_S}{\|x_R-e^{-A\sigma}x_S\|} + O(\sigma^2).\]After some straightforward calculations, this yields\[\tag{10}N – E = \frac{\sigma}{\|x_R-x_S\|}(Ax_S – \langle Ax_S, E\rangle E) + O(\sigma^2).\]
We can apply this result to the first summand in the right hand side of (9) to yield the approximation of first order in \(\sigma\)\[\begin{split}\langle x’_R – x’_S, N\rangle \approx &\left( 1 – \frac{\sigma \langle Ax_S, E\rangle}{\|x_R-x_S\|}\right) \langle x’_R – x’_S, E \rangle \\ &+ \frac{\sigma}{\|x_R-x_S\|}\langle x’_R – x’_S, Ax_S\rangle.\end{split}\tag{11}\]
For the second summand in the right hand side of (9) we do a first order approximation of the exponential:\[-\langle (e^{-A\sigma}-I)x’_S, N\rangle \approx \sigma \langle Ax’_S, N\rangle.\]Since this is already a term of order one in \(\sigma\), we may replace \(N\) by \(E\) directly, obtaining\[\sigma\langle Ax’_S, N\rangle \approx \frac{\sigma}{\|x_R-x_S\|} \langle Ax’_S, x_R-x_S \rangle.\]
Putting together this last result with (9) and (11), and using \(A = -A^T\) we get\[\begin{split}\langle x’_R – e^{-A\sigma}x’_S, N \rangle \approx &\left( 1 – \frac{\sigma \langle Ax_S, E\rangle}{\|x_R-x_S\|}\right)\langle x’_R-x’_S, E\rangle \\ &+ \frac{\sigma}{\|x_R-x_S\|}(\langle x’_R, Ax_S\rangle + \langle Ax’_S, x_R\rangle).\end{split}\tag{12}\]
Now we need to take care of the second summand in the right hand side of (8), which is somewhat more delicate. Since \(\tau’ = \rho’ \tau / \rho\), we can write\[\tau’ \omega_E\langle A e^{-A\sigma}x_S, N \rangle = \rho’ \sigma \left[\frac{1}{\rho}\langle A e^{-A\sigma}x_S, N \rangle\right].\]Let us write \(C_2\) for the term in square brackets in the equation above and let us write the right hand side of (12) as \(\langle x’_R-x’_S, E\rangle + \sigma C_1\), so that\[\rho’ = \langle x’_R-x’_S,E\rangle + \sigma C_1 + \rho’ \sigma C_2 + O(\sigma^2)\]according to (8). Applying this equality to the right hand side again we get\[\begin{split}\rho’ = &\langle x’_R-x’_S,E\rangle + \sigma C_1 + \sigma C_2 \langle x’_R-x’_S,E\rangle \\ &+ \sigma^2 C_1C_2 + \rho’\sigma^2C_2^2 + C_2O(\sigma^3) + O(\sigma^2).\end{split}\]
Let us examine the order of magnitude of the summands in the right hand side of this equality. The term \(C_1\) has an order of magnitude bounded by that of \(x’_R\) and \(x’_S\). The term \(C_2\) is smaller: it can be bounded by a numeric constant (probably 3 is sufficient as this constant). Therefore, the terms \(\sigma C_1\) and \(\sigma C_2 \langle x’_R-x’_S, E \rangle\) can be grouped together as our correction of first order in \(\sigma\). All the terms in the lower row, can be subsummed under the \(O(\sigma^2)\). Hence,\[\rho’ = \langle x’_R-x’_S,E\rangle + \sigma C_1 + \sigma C_2 \langle x’_R-x’_S,E\rangle + O(\sigma^2).\]
Since the term \(C_2\) is multiplied by \(\sigma\), we can replace it by an approximation of order zero in \(\sigma\). To this end, \(\rho\) can be replaced by \(\|x_R – x_S\|\), the exponential \(e^{-A\sigma}\) can be replaced by the identity matrix by setting \(\sigma = 0\), and \(N\) can be replaced by \(E\), obtaining\[C_2 = \frac{1}{\rho}\langle A e^{-A\sigma}x_S, N \rangle \approx \frac{\langle A x_S, E\rangle}{\|x_R-x_S\|},\]so that we get the approximation\[\begin{split}\tau’ \omega_E \langle A e^{-A\sigma}x_S, N \rangle &= \rho’\sigma C_2 \\ &= \sigma \langle x’_R – x’_S, E\rangle C_2 + O(\sigma^2)\\ &\approx \frac{\sigma\langle A x_S, E\rangle}{\|x_R-x_S\|}\langle x’_R – x’_S, E\rangle.\end{split}\tag{13}\]Note that this last expression also appears in (12) with the opposite sign, so they will cancel out when (12) and (13) are added.
Indeed, using (8) to add (12) and (13), we get\[\rho’ \approx \langle x’_R-x’_S, E\rangle + \frac{\sigma}{\|x_R-x_S\|}(\langle x’_R, A x_S\rangle + \langle Ax’_S, x_R\rangle).\]Replacing \(\sigma/\|x_R-x_S\|\) by \(\omega_E/c\) as we did in the derivation of the approximation for the range correction, we get equation (6), as we wanted to show.
Geometric interpretation of the corrections
It is interesting to note that all the correction terms are of the form \(\omega_E\langle A u, v\rangle/c\). This can be interpreted geometrically as\[\tag{14}\langle A u, v\rangle = 2 \operatorname{area}(P_z u, P_z v),\]where \(P_z\) denotes the orthogonal projection onto the plane \(z = 0\) and \(\operatorname{area}(u,v)\) is the signed area of the triangle which has a vertex at zero and sides given by the vectors \(u\) and \(v\). The sign of the area is chosen as the orientation of the basis \(\{u,v\}\) (i.e., as the sign of the determinant of \(u\) and \(v\)). In fact, the easiest way to check (14) is to note that\[\langle A u, v\rangle = \operatorname{det}(P_z u, P_z v).\]
This geometric interpretation can be used to estimate the order of magnitude of the corrections in different geometric configurations. For the range, the correction is\[\frac{\omega_E}{c}\langle Ax_S, x_R\rangle.\]For a receiver on the Earth’s surface, this term is largest when the receiver is on the equator and the satellite is in the equatorial plane and seen at an elevation of 0º by the receiver. In this case, the vectors \(x_R\) and \(x_R-x_S\) are orthogonal.
Since we can also write this correction as\[\frac{\omega_E}{c}\langle A(x_S-x_R), x_R\rangle\]due to the fact that \(A^T = -A\), in this case the absolute value of the correction works out to be\[\frac{\omega_E}{c}\|x_S-x_R\|\|x_R\| = \frac{\omega_E}{c}R_E\sqrt{a^2 – R_E^2},\]where \(a\) is the orbital radius of the satellite (we assume a circular orbit) and \(R_E\) is the Earth’s radius. For a GPS-like MEO orbit with a period of half a sidereal day, the correction amounts to 40 metres. Hence, we see that it is critical to perform this kind of correction when doing most GNSS calculations involving the range.
For the range-rate, we have two correction terms,\[\frac{\omega_E}{c}\langle Ax_S, x’_R\rangle,\quad\mathrm{and}\quad \frac{\omega_E}{c}\langle A x’_S, x_R\rangle.\]The first term depends on the receiver velocity, but it is usually much smaller than the second one. For a GPS-like MEO orbit and a receiver at the vicinity of the Earth’s surface, \(\|x_S\|\) is approximately 4 times larger than \(\|x_R\|\). However, \(\|x’_S\| = 3874\,\mathrm{m/s}\), so in most situations \(\|x’_R\|\) is much smaller than \(\|x’_S\|\).
Regarding the second term, for a receiver on the Earth’s surface, this correction is largest when the receiver is on the equator, the satellite has an orbital inclination of 0º (I reckon that most GNSS constellations have inclinations around 50º, but I use this example for simplicity) and the satellite is at the zenith of the receiver. In this situation \(x_R\) and \(x’_S\) form right angles, so the correction amounts to\[\frac{\omega_E}{c}\|x’_S\|\|x_E\| = \frac{\omega_E}{c}R_E\sqrt{\frac{\mu}{a}},\]where \(\mu = GM_E\) is Earth’s gravitational parameter. For a MEO orbit with a period of half a sidereal day, this works out to be 6mm/s, so it is a very small correction. Therefore, the Earth rotation correction can be ignored in most GNSS calculations involving the range-rate.
2 comments