2. The snapshot-positioning problem: a GNSS model with whole-millisecond ambiguities

We begin by showing that when departure times are not known, the observation equations that relate the arrival times of GNSS codes to the unknown position of the receiver contain integer ambiguities.

We denote by $t_{D,i}$ the time of departure of a code from satellite $i$, and we assume that $t_{D,i}$ represents a whole millisecond (in the time base of the GPS system). We denote by $\mathring{t}_{i}$ the time of arrival of that code at the antenna of the receiver. We assume that the receiver estimates the arrival time of that code as $t_{i}=\mathring{t}_{i}+\mathring {b}+\epsilon _{i}$, where $\mathring {b}$ represents the bias that is caused by the inaccurate clock of the receiver and by delays that the signal experiences in the path from the antenna to the analogue-to-digital converter, and $\epsilon _{i}$ is the arrival-time estimation error. The bias $\mathring {b}$ is time dependent, because of drift in the receiver's clock, but over short observation periods this dependence is negligible, so we ignore it.

The time of arrival is governed by the equation

\[ \mathring{t}_{i}-t_{D,i}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\delta_{i}, \]

where $c$ is the speed of light, $\mathring {\ell }$ is the location of the receiver, $\rho _{i}$ is the location of the satellite (which is a function of time, since the satellites are not stationary relative to Earth observers), and $\delta _{i}$ represents the inaccuracy of the satellite's clock and atmospheric delays. We assume that $\delta _{i}$ can be modelled, for example using models of ionospheric and tropospheric delays (dual frequency receivers can estimate the ionospheric delay, but we assume a single-frequency receiver). Setting $\delta _{i}$ to the satellite's clock error correction from the ephemeris induces a location error of approximately 30 m due to the atmospheric delays (Borre and Strang, Reference Borre and Strang2012).

A receiver that decodes the timestamp embedded in the GPS data stream can determine $t_{D,i}$, which leads to the following equation:

\[ t_{i}-t_{D,i}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\mathring{b}+\delta_{i}+\epsilon_{i}, \]

the conventional GNSS code-observation equation, in which the four unknown parameters are $\mathring {b}$ and the coordinates of $\mathring {\ell }$ (we assume that $\delta _{i}$ is modelled, possibly trivially $\delta _{i}=0$, but not estimated). To simplify the notation, we ignore $\delta _{i}$ for now and write

\[ t_{i}-t_{D,i}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\mathring{b}+\epsilon_{i}. \]

We use observations from all the satellites such that all the $t_{i}$ data lie between two consecutive whole multiples of $t_{\text {code}}$ (in GPS, two round milliseconds in the local clock). This allows us to express

\[ t_{i}=(N+\varphi_{i})t_{\text{code}} \]

with a common and easily computable $N=\lfloor t_{i}/t_{\text {code}}\rfloor$ and for $\varphi _{i}\in [0,1)$. We denote $N_{i}=t_{D,i}/t_{\text {code}}$ and write

\[ (N-N_{i}+\varphi_{i})t_{\text{code}}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\mathring{b}+\epsilon_{i}. \]

Since GNSS codes are aligned with $t_{\text {code}},$ $N_{i}\in \mathbb {Z}$. We denote $n_{i}=N-N_{i}\in \mathbb {Z}$, so

\[ (n_{i}+\varphi_{i})t_{\text{code}}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\mathring{b}+\epsilon_{i} \]

or

\[ \varphi_{i}t_{\text{code}}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}. \]

We now face two challenges. One is that we have $4+m$ unknown parameters: three location coordinates, $\mathring {b}$ and $n_{i}$, but only $m$ constraints. We clearly need more constraints so that we can resolve $n_{i}$. The other is that we have a set of nonlinear constraints with continuous real unknowns, the location and $\mathring {b}$, and with integer unknowns, the $n_{i}$. The strategy, as in other cases with this structure, is to first linearise the nonlinear term, then to resolve the integer parameters, and to then substitute them and to solve the continuous least-squares problem (either the linearised system or the original nonlinear system). We cannot linearise the nonlinear term $\|\mathring {\ell }-\rho _{i}(t_{D,i})\|_{2}/c$ using a Taylor series because it is a function of both real unknowns and of the integer unknowns $n_{i}$. We cannot differentiate this term with respect to the integer $n_{i}$.

To address this difficulty, we approximate $t_{D,i}$ by approximating the range (distance) term in the equation

\[ t_{D,i}=t_{i}-\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}-\mathring{b}-\epsilon_{i}. \]

For now, we denote the approximation of the propagation delay by

\[ d_{i}\approx\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\epsilon_{i} \]

so

\[ \hat{t}_{D,i}=t_{i}-d_{i}-\mathring{b}. \]

There are several ways to set $d_{i}$, depending on our prior knowledge of $\mathring {\ell }$ and $\mathring {b}$. One option in the GPS system is to set it to approximately 76$\cdot$5 ms; this limits the error in $\hat {t}_{D,i}$ to approximately 12$\cdot$5 ms for any Earth observer, and the error

\[ \Vert \mathring{\ell}-\rho_{i}(\hat{t}_{D,i})\Vert_{2}-\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2} \]

to approximately 10 m (Van Diggelen, Reference Van Diggelen2009). We substitute $\rho _{i}(\hat{t}_{D,i})=\rho _{i}(t_{i}-d_{i}-\mathring {b})$ for $\rho _{i}(t_{D,i})$,

(2.1)\begin{equation} \varphi_{i}t_{\text{code}}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{i}-d_{i}-\mathring{b})\Vert_{2}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}^{(D)}.\end{equation}

The superscript $(D)$ on the error term indicates that the error term now represents not only the arrival-time estimation error, but also the error induced by the inexact departure time.

We linearise around an a priori solution $\bar {\ell }$ and $\bar {b}$ (usually $\bar {b}=0$, otherwise we can simply shift the $t_{i}$s),

(2.2)\begin{align} \varphi_{i}t_{\text{code}} & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}^{(D,L)} \nonumber\\ & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix}\mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-n_{i}t_{\text{code}}+(\mathring{b}-\bar{b})+\bar{b}+\epsilon_{i}^{(D,L)}, \end{align}

where $\mathrm {J}$ is the Jacobian of the Euclidean distances with respect to both the location of the receiver and to the bias, with the derivatives evaluated at $\bar {\ell }$ and at $t_{i}-d_{i}-\bar {b}$. The superscript $(D,L)$ on the error term indicates that it now includes also the linearisation error.

There are now several ways to resolve the $n_{i}$.

3. Shadowing

Peterson et al. (Reference Peterson, Hartnett and Ottman1995) introduced a somewhat surprising modelling technique, which we refer to as shadowing. The idea is to replace the unknown $b$ by two separate unknowns that represent essentially the same quantity, the original $b$ and a shadow $s$. In principle, they should obey the equation $b=s$, but the model treats $s$ as a free parameter; the constraint $b=s$ is dropped. In the literature, $s$ is called the coarse-time parameter (and is often represented by $tc$ or $t_{c}$). We express this technique by splitting $b$ and $s$:

\begin{align*} \varphi_{i}t_{\text{code}} & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}^{(D,L)}\\ & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ s-\bar{b} \end{bmatrix}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}^{(D,L)}. \end{align*}

We now have five unknowns, not four.

As far as we can tell, there is no clear explanation in the literature as to the benefits of shadowing. One way to justify the technique is to observe that Equation (2.2) is very sensitive to small (nanosecond scale) perturbations in the additive $\mathring {b}$, but it is not highly sensitive to the $\mathring {b}$ (now $s$) by which we multiply,

(3.1)\begin{equation} \mathrm{J}_{i,4}=\frac{\partial}{\partial t_{D,i}}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}.\end{equation}

For example, in GPS, the derivative is bounded by approximately $800$ m/s for any $\bar {\ell }$ on Earth (Van Diggelen, Reference Van Diggelen2009), so $\mathrm {J}_{i,4}/c<3\times 10^{-6}$ (versus $1$ for the additive $\mathring {b}$). Therefore, the dependence of the residual (the vector of $\epsilon _{i}^{(D,L)}$ terms for a given setting of the unknown parameters) on $\mathring {b}$ in Equation (2.2) is highly non-convex. There are many different values of $\mathring {b}$ that are almost equally good, a millisecond apart, with each of these nearly optimal hypotheses being locally well defined; if we increase $\mathring {b}$ by one millisecond and also add $1$ to each $n_{i}$, the residual changes very little, because $\mathrm {J}{}_{i,4}$ is so small.

Shadowing turns this non-convexity into explicit rank deficiency, with which it is easier to deal. With one instance of $\mathring {b}$ replaced by the shadow $s$, the constraints no longer uniquely define $\mathring {b}$, only up to a multiple of $t_{\text {code}}$. For any hypothetical solution $\ell ,s,b,n$, the solution $\ell ,s,b+kt_{\text {code}},n+k$ gives exactly the same residual. We perform a change of variables, replacing the partial sum $-n_{i}t_{\text {code}}+\mathring{b}$ by $-\nu _{i}t_{\text {code}}+\beta$, where $-\nu _{i}=-n_{i}+\lfloor \mathring{b}/t_{\text {code}}\rfloor$ and $\beta =\mathring {b}-\lfloor b/t_{\text {code}}\rfloor t_{\text {code}}$, so $\beta \in [0,t_{\text {code}})$:

(3.2)\begin{equation} \begin{aligned} \varphi_{i}t_{\text{code}} & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ s-\bar{b} \end{bmatrix}-\nu_{i}t_{\text{code}}+\beta+\epsilon_{i}^{(D,L)}.\\ \beta & \in [0,t_{\text{code}}). \end{aligned} \end{equation}

3.1 Resolving the integer ambiguities: Van Diggelen's method

Van Diggelen's method exploits the fact that the $\nu _{i}$ values are very insensitive to $\mathring {\ell }$ and to $s$. It therefore sets $\mathring {\ell }=\bar {\ell }$ and $s=\bar {b}$, truncating the Jacobian term from Equation (3.2):

(3.3)\begin{equation} \begin{aligned} \varphi_{i}t_{\text{code}} & = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}-\nu_{i}t_{\text{code}}+\beta+\epsilon_{i}^{(D,L,A)}\\ \beta & \in [0,t_{\text{code}}). \end{aligned} \end{equation}

The new subscript $(D,L,A)$ indicates that the error term now compensates also for the use of the a priori estimates $\bar {b}$ and $\bar {\ell }$ for $s$ and $\mathring {\ell }$.

Van Diggelen uses these constraints to set the $\nu _{i}$ values in a particular way. The method selects one index $j$ that is used to set $\nu _{j}$ and $\beta$, and then resolves all the other $\nu _{i}$ values so they are consistent with this $\beta$. That is, he assumes that $\epsilon _{j}^{(D,L,A)}=0$ so

\begin{align*} \nu_{j} & = \left\lceil \frac{\frac{1}{c}\Vert \bar{\ell}-\rho_{j}(t_{j}-d_{j}-\bar{b})\Vert -\varphi_{j}t_{\text{code}}}{t_{\text{code}}}\right\rceil\\ \beta & = (\nu_{j}+\varphi_{j})t_{\text{code}}-\frac{1}{c}\Vert \bar{\ell}-\rho_{j}(t_{j}-d_{j}-\bar{b})\Vert_{2}. \end{align*}

The method now substitutes this $\beta$ in all the other constraints and assigns the other $\nu _{i}$ values by setting $\epsilon _{i}^{(D,L,A)}=0$ and rounding,

(3.4)\begin{equation} \nu_{i}=\left\lfloor \frac{\frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert -\varphi_{i}t_{\text{code}}+\beta}{t_{\text{code}}}\right\rceil .\end{equation}

When $\|\mathring {\ell }-\bar {\ell }\|$ and $|s-\bar {b}|$ are small enough, this gives a set of $\nu _{i}$ values that are correct in the sense that they all differ from the correct $\nu _{i}$ values by the same integer.

Van Diggelen chooses $j$ in a particular way: he chooses the $j$ that minimises the magnitude of Equation (3.1), which corresponds to the satellite closest to the zenith of $\bar {\ell }$ at $t_{j}-\bar {b}$. In our notation, Van Diggelen's justification for this choice is as follows. He searches for a $j$ for which Equation (3.3) approximates well Equation (3.2). The difference between the two is

\[ \frac{1}{c}\mathrm{J}_{j,:}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ s-\bar{b} \end{bmatrix}. \]

For each satellite, $\mathrm {J}_{i,1:3}$ is the negation of the so-called line-of-sight vector $e_{i}^{T}$, which is the normalised direction from the satellite to the receiver; element $\mathrm {J}_{i,4}$ is the range rate. Van Diggelen's choice of $j$ leads to a row of $\mathrm {J}$ in which the first three elements are almost orthogonal to $\mathring {\ell }-\bar {\ell }$ and in which the fourth element, the range rate, is small. This leads to an estimated $\beta$ that is relatively accurate, which helps resolve the correct $\nu _{i}$ values.

Van Diggelen also shows that if we resolve the $\nu _{i}$ values by setting each $\epsilon _{i}^{(D,L,A)}=0$ separately, then the resolved $\beta$ values might be close to $0$ in one equation and close to $t_{\text {code}}$ in another; this leads to inconsistent $\nu _{i}$ values and to a huge position error.

3.2 Final resolution of the receiver's location

Van Diggelen's method resolves the integer $\nu _{i}$ values in Equation (3.3). Now we need to resolve the continuous unknowns. We do so using Gauss–Newton iterations on Equation (3.2), iterating on $\delta _{\ell }=\mathring {\ell }-\bar {\ell }$, $\delta _{s}=s-\bar {b}$ and $\beta$ but keeping $\nu$ fixed. We start with $\delta _{\ell }$, $\delta _{s}$ and $\beta$ set to zero.

In every iteration, we use the current iterates to produce estimates of the location and bias,

\begin{align*} \hat{\ell} & = \bar{\ell}+\delta_{\ell}\\ \hat{b} & = \bar{b}+\delta_{s}. \end{align*}

We use them to improve the estimate of the ranges $d_{i}$, setting

\[ \hat{d}_{i}=\Vert \hat{\ell}-\rho_{i}(\hat{t}_{D,i})\Vert_{2}. \]

This allows us to reduce the errors in Equation (2.1),

\begin{align*} \varphi_{i}t_{\text{code}} & = \frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{i}-\hat{d}_{i}-\mathring{b})\Vert_{2}-n_{i}t_{\text{code}}+\mathring{b}+\epsilon_{i}^{(\hat{D})}\\ & = \frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{i}-\hat{d}_{i}-\mathring{b})\Vert_{2}-\nu_{i}t_{\text{code}}+\beta+\epsilon_{i}^{(\hat{D})} \end{align*}

(the second line holds because $n_{i}t_{\text {code}}+\mathring {b}=\nu _{i}t_{\text {code}}+\beta$, by definition). We again linearise this and solve the constraints

(3.5)\begin{equation} \varphi_{i}t_{\text{code}}=\frac{1}{c}\Vert \hat{\ell}-\rho_{i}(t_{i}-\hat{d}_{i}-\hat{b})\Vert_{2}+\frac{1}{c}\hat{\mathrm{J}}_{i,:} \begin{bmatrix} \mathring{\ell}-\hat{\ell}\\ s-\hat{b} \end{bmatrix}-\nu_{i}t_{\text{code}}+\beta+\epsilon_{i}^{(\hat{D},L)}.\end{equation}

for $\mathring {\ell }$, $s$ and $\beta$ using in the generalised least-squares sense, where the Jacobian is evaluated at $\hat {\ell }$ and $t-\hat {d}-\hat {b}$.

We can now explain why Van Diggelen's method resolves the integers only once and iterates only on the continuous unknowns. The $\nu _{i}$ values that Van Diggelen's method resolves are not equal to the $n_{i}$ values in the nonlinear Equation (2.1). But when the linearisation error is small enough, the two integer vectors differ by a constant, $\lfloor \mathring{b}/t_{\text {code}}\rfloor$. This difference is compensated for by the integer part of the continuous variable $\beta$, which is not constrained to $[0,t_{\text {code}})$ in the Gauss–Newton iterations. This is the actual function of shadowing; to allow $\beta$ to compensate not only for the clock error, but also for the constant error in $\nu$. When the initial linearisation error is so large that $n-\nu$ is no longer a constant, the method breaks down.

4. A mixed-integer least-squares approach

A different approach, which has never been proposed for snapshot positioning, is to add regularisation constraints that will allow us to resolve all the $4+m$ unknowns in the $m$ instances of Equation (2.2)

\[ \varphi_{i}t_{\text{code}} = \frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-(\mathring{n}_{i}-\bar{n}_{i})t_{\text{code}}-\bar{n}_{i}t_{\text{code}} +(\mathring{b}-\bar{b})+\bar{b}+\epsilon_{i}^{(D,L)} \]

using mixed-integer least-squares techniques. Note that we have rewritten Equation (2.2) in a way that emphasises a change of variables that facilitate iterative improvements: the new unknowns are $\mathring {\ell }-\bar {\ell }$, $\mathring {b}-\bar {b}$ and $\mathring {n}-\bar {n}$. We initially set $\bar {b}$ and $\bar {n}$ to zero.

We denote the vector of delays by $g$,

\[ g_{i}=\frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}. \]

This section proposes two sets of regularising equations and explains how to use this approach in an iterative Gauss–Newton solver.

4.1 Resolving the ambiguities: regularisation using a priori estimates

The first set of regularising equations that we propose are

\[ \frac{1}{c}\mathrm{J}_{i,:}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}=0. \]

We do not enforce them exactly, only in a (weak) least-squares sense. They favour solutions of the mixed-integer least-squares problem that are in the vicinity of the a priori solution. This leads to the following weighted mixed-integer least-squares problem:

\begin{align*} \hat{\ell},\hat{b},\hat{n}&=\arg\min_{\ell,b,n} \left\Vert W\left(\left[\begin{array}{cc} \dfrac{1}{c}\mathrm{J}+[\mathbf{0}_{m\times 3} \quad \mathbf{1}_{m\times 1}] & -t_{\text{code}}I_{m\times m}\\[6pt] \dfrac{1}{c}\mathrm{J} & \mathbf{0}_{m\times m} \end{array}\right]\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ \mathring{n}-\bar{n} \end{bmatrix}\right.\right.\\ &\quad \left.\left.-\begin{bmatrix} t_{\text{code}}\varphi-g+\bar{n}t_{\text{code}}-\bar{b}\\ \mathbf{0} \end{bmatrix}\right)\right\Vert_{2}^{2}, \end{align*}

where $W$ is a block-diagonal weight matrix derived from the covariance matrix $C$ of the error terms $\epsilon$, $W^{T}W=C^{-1}$. Now we have $2m$ constraints, which for $m\geq 4$ should allow us to resolve the integer $n_{i}$ values.

We propose to choose a diagonal $W$ as follows. We set the first $m$ diagonal elements of $W$ to the standard deviation of the arrival-time estimator, say $W_{i,i}=1/\sigma (t_{i})\approx 1/(10\,\text {ns})$. To set the rest, we use box constraints on the a priori estimates $\bar {\ell }$ and $\bar {b}$, denoted as

\begin{align*} |x-\bar{x}| & \leq x_{\max}\\ |y-\bar{y}| & \leq y_{\max}\\ |z-\bar{z}| & \leq z_{\max}\\ |\mathring{b}-\bar{b}| & \leq b_{\max}. \end{align*}

By the triangle inequality

\[ \left|\mathrm{J}_{i,:}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}\right|\leq |\mathrm{J}_{i,1}|x_{\max}+|\mathrm{J}_{i,2}|y_{\max} +|\mathrm{J}_{i,3}|z_{\max}+|\mathrm{J}_{i,4}|b_{\max}. \]

We define

\[ r_{i}=|\mathrm{J}_{i,1}|x_{\max}+|\mathrm{J}_{i,2}|y_{\max}+|\mathrm{J}_{i,3}|z_{\max} +|\mathrm{J}_{i,4}|b_{\max}, \]

so

\[ \left|\mathrm{J}_{i,:}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}\right|\leq r_{i}. \]

We convert the hard box constraints into soft-weighted least-squares in order to allow for using a mixed-integer least-squares solver. We need to set $W_{m+i,m+i}$; if we assume that the error in the constraint is Gaussian and that an error of $r_{i}/c$ is acceptable (from the inequality above), then setting $W_{m+i,m+i}=c/r_{i}$, say, makes sense. In practice, we use $W_{m+i,m+i}=c/(100\,\text {km})$ in the experiments below.

This mixed-integer least-squares minimisation problem can be solved by a generic solver, such as one of the solvers that have been developed for RTK.

4.2 Doppler regularisation

It turns out that Doppler shifts allow us to regularise Equation (2.2) in a more effective way. GNSS receivers estimate not only the time of arrival of the signal, but also its Doppler shift. The estimated Doppler shift is biased, because of the inaccuracy of the receiver's local (or master) oscillator; it is also inexact. We now show a novel technique to use the Doppler-shift observations to to regularise Equation (2.2).

Our technique is based on two assumptions. One is that the receiver is stationary, or more precisely, that its velocity is negligible relative to the range rate, which is up to approximately 800 m/s. This assumption can be easily removed, but its removal leads to additional unknowns and more complicated expressions that we do not present here. The other assumption is that the local oscillator and the sampling clock in the receiver are derived from a single master oscillator in a certain (very common) way. Again, this assumption can be removed if another unknown is added.

The Doppler-shift formula for velocities much lower than the speed of light is

\[ \mathring{D}_{i}\approx{-}\frac{1}{c}\frac{d}{dt}\Vert \mathring{\ell}-\rho_{i}\Vert_{2}f_{0}. \]

The Doppler observations that the receiver makes are

(4.1)\begin{equation} D_{i}={-}\frac{1}{c}\frac{d}{dt}\Vert \mathring{\ell}-\rho_{i}\Vert_{2}f_{0}+\mathring{f}+\epsilon_{i}^{(\delta)},\end{equation}

where $\mathring {f}$ is the frequency offset (bias) of the receiver and $\epsilon _{i}^{(\delta )}$ is an error term that represents the observation error and the (negligible) slow-speed approximation. Therefore, the quantities $-cD_{i}/f_{0}$ are biased estimates of the range rate. We denote the a priori estimates of the Doppler shifts by $\bar {D}_{i}$.

We differentiate Equation (2.2) by time,

\[ \frac{d}{dt}(\varphi_{i}t_{\text{code}}) = \frac{d}{dt}\left(\frac{1}{c}\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}+\frac{1}{c}\mathrm{J}_{i,:} \begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-n_{i}t_{\text{code}}+(\mathring{b}-\bar{b})+\bar{b}\right)+\epsilon_{i}^{(D,L,\partial)}. \]

We first manipulate the equation a little, to make it easier to differentiate:

(4.2)\begin{align} &\frac{d}{dt}(c\varphi_{i}t_{\text{code}}-\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}-c\bar{b})\nonumber\\ &\quad = \frac{d}{dt}\left(\mathrm{J}_{i,:}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}-cn_{i}t_{\text{code}}+c(\mathring{b}-\bar{b})\right)+\epsilon_{i}^{(D,L,\partial)}. \end{align}

We denote

\[ H=\left[\begin{array}{ccccccccc} \mathrm{J}_{1,1} & \mathrm{J}_{1,2} & \mathrm{J}_{1,3} & \mathrm{J}_{1,4}+c & -ct_{\text{code}}\\ \vdots & \vdots & \vdots & \vdots & & \ddots\\ \mathrm{J}_{i,1} & \mathrm{J}_{i,2} & \mathrm{J}_{i,3} & \mathrm{J}_{i,4}+c & & & -ct_{\text{code}}\\ \vdots & \vdots & \vdots & \vdots & & & & \ddots\\ \mathrm{J}_{m,1} & \mathrm{J}_{m,2} & \mathrm{J}_{m,3} & \mathrm{J}_{m,4}+c & & & & & -ct_{\text{code}} \end{array}\right]. \]

The first three columns of $H$ are identical to those of $\mathrm {J}$, the next is the fourth column of $\mathrm {J}$ but shifted by $c$ and the last $m$ columns consist of a scaled identity matrix. We now express the derivative on the right-hand side of Equation (4.2) as

\[ \frac{d}{dt}\left(H\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ n_{1}\\ \ldots\\ n_{m} \end{bmatrix}\right)=H\left(\frac{d}{dt}\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ n_{1}\\ \ldots\\ n_{m} \end{bmatrix}\right)+\left(\frac{d}{dt}H\right)\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ n_{1}\\ \ldots\\ n_{m} \end{bmatrix}. \]

We assume that the receiver is stationary, so $\mathring {\ell }-\bar {\ell }$ is time-independent, so $({d}/{dt})(\mathring {\ell }-\bar {\ell })=0$. The derivatives of the integers $n_{1},\ldots ,n_{m}$ are also zero. The derivative of the remaining element in the vector, $({d}/{dt})(\mathring {b}-\bar {b})$, is not zero and will need to be estimated. It represents the frequency offset of the receiver, which biases the observed Doppler shift. It is multiplied by a column whose elements are very close to $c$ (the range rate is tiny relative to the speed of light), which allows it to compensate for the frequency bias.

To differentiate $H$, we exploit the known structure of $\mathrm {J}$. For each satellite, $\mathrm {J}_{i,1:3}$ is the negation of the so-called line-of-sight vector $e_{i}^{T}$, which is the normalised direction from the satellite to the receiver; element $\mathrm {J}_{i,4}$ is the range rate. The derivatives of these quantities are shown by Van Diggelen (Reference Van Diggelen2009, Equation 8.6), Fernández-Hernández and Borre (Reference Fernández-Hernández and Borre2016) and other sources:

\[ \frac{d}{dt}\mathrm{J}_{i,1:3} = \frac{e^{T}\frac{d}{dt}\|\bar{\ell}-\rho_{i}\| +\left(\frac{d}{dt}(\bar{\ell}-\rho_{i})\right)^{\textrm{T}}}{\|\bar{\ell}-\rho_{i}\|}, \]

where the satellite position $\rho _{i}$ and its velocity $({d}/{dt})\rho _{i}$ are taken at $\hat {t}_{D}$. To reduce the number of unknowns, we assume that the receiver is stationary, so $({d}/{dt})\bar {\ell }=0$, so

\[ \frac{d}{dt}\mathrm{J}_{i,1:3} = \frac{e^{T}\frac{d}{dt}\|\bar{\ell}-\rho_{i}\| -\left(\frac{d}{dt}\rho_{i}\right)^{\textrm{T}}}{\|\bar{\ell}-\rho_{i}\|}. \]

Element $\mathrm {J}_{i,4}$ is the range rate of satellite $i$, so its derivative with respect to time is the satellite's range acceleration,

\[ \frac{d}{dt}\mathrm{J}_{i,4} =\frac{d^{2}}{dt^{2}} \|\bar{\ell}-\rho_{i}\|. \]

We use finite differences to evaluate this second derivative. The fourth column of $H$ is $\mathrm {J}_{:,4}+c$, but the derivative of $c$ is obviously zero. The derivative of $-ct_{\text {code}}$ is also zero, so

\[ \left(\frac{d}{dt}H\right)\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ n_{1}\\ \ldots\\ n_{m} \end{bmatrix}=\left(\frac{d}{dt}H_{:,1:4}\right)\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}. \]

We now derive the left-hand side of Equation (4.2),

\[ \frac{d}{dt}(c\varphi_{i}t_{\text{code}}-\Vert \bar{\ell}-\rho_{i}(t_{i}-d_{i}-\bar{b})\Vert_{2}-c\bar{b}). \]

The derivative of the a prioi range estimate $\Vert \bar {\ell }-\rho _{i}(t_{i}-d_{i}-\bar {b})\Vert _{2}$ is the a priori range rate, which we can compute. The derivative of $c\bar {b}$ is zero.

To understand the first term, recall that

\[ (n_{i}+\varphi_{i})t_{\text{code}}=\frac{1}{c}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+\mathring{b}+\epsilon_{i}, \]

so

\[ c\frac{d}{dt}\varphi_{i}t_{\text{code}}=\frac{d}{dt}\Vert \mathring{\ell}-\rho_{i}(t_{D,i})\Vert_{2}+c\frac{d}{dt}\mathring{b}+c\frac{d}{dt}\epsilon_{i}. \]

We now rewrite Equation (4.1) as

\[ \frac{1}{c}\frac{d}{dt}\Vert \mathring{\ell}-\rho_{i}\Vert_{2}f_{0}={-}D_{i}+\mathring{f}+\epsilon_{i}^{(\delta)} \]

or

\[ \frac{d}{dt}\Vert \mathring{\ell}-\rho_{i}\Vert_{2}={-}\frac{cD_{i}}{f_{0}}+\frac{c\mathring{f}}{f_{0}}+\frac{c}{f_{0}}\epsilon_{i}^{(\delta)}. \]

We now substitute in the left-hand side of Equation (4.2):

\[ \frac{d}{dt}c\varphi_{i}t_{\text{code}} ={-}\frac{cD_{i}}{f_{0}}+\frac{c\mathring{f}}{f_{0}}+\frac{c}{f_{0}}\epsilon_{i}^{(\delta)}+c\frac{d}{dt}\mathring{b}+c\frac{d}{dt}\epsilon_{i}. \]

The term $\mathring {f}/f_{0}$ is the relative local-oscillator error in the receiver. If the oscillator runs too fast, $\mathring {f}$ is negative. Assuming that all the clocks in the receiver are derived from a master oscillator, if it runs too fast, $\mathring {b}$ grows over time. Under this assumption,

\[ -\frac{c\mathring{f}}{f_{0}}=c\frac{d}{dt}\mathring{b}, \]

so these terms cancel each other. If our assumption on the receiver does not hold, we would need to estimate

\[ \frac{c\mathring{f}}{f_{0}}+c\frac{d}{dt}\mathring{b}. \]

We have arrived at a system of $m$ linear equations that we use to regularise the mixed-integer equations. The equations are:

(4.3)\begin{equation} -\frac{c}{f_{0}}D-\frac{d}{dt}\Vert \mathring{\ell}-\rho\Vert_{2} =H_{:,4}(\mathring{u}-\bar{u})+\left(\frac{d}{dt}H_{:,1:4}\right) \begin{bmatrix}\mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b} \end{bmatrix}.\end{equation}

In this equation, $D$ represents the vector of observed Doppler shifts, $({d}/{dt})\Vert \mathring {\ell }-\rho \Vert _{2}$ is the vector of the a priori range rates and $\mathring {u}-\bar {u}=({d}/{d\,t})(\mathring {b}-\bar {b})$ is a new scalar unknown. We have explained above how to compute $H_{:,4}$ and $({d}/{dt})H_{:,1:4}$. The full regularised weighted least-squares that we solve is

\begin{align*} \hat{\ell},\hat{b},\hat{n},\hat{u}&=\arg\min_{\ell,b,n,u} \left\Vert W\left(\left[\begin{array}{ccc} \dfrac{1}{c}\mathrm{J}+[\mathbf{0}_{m\times 3} \quad \mathbf{1}_{m\times 1}] & -t_{\text{code}}I_{m\times m} & \mathbf{0}_{m\times 1}\\ \left(\dfrac{d}{dt}H_{:,1:4}\right) & \mathbf{0}_{m\times m} & H_{:,4} \end{array}\right]\begin{bmatrix} \mathring{\ell}-\bar{\ell}\\ \mathring{b}-\bar{b}\\ \mathring{n}-\bar{n}\\ \mathring{u}-\bar{u} \end{bmatrix}\right.\right.\\ &\quad \left.\left.-\begin{bmatrix} t_{\text{code}}\varphi-g+\bar{n}t_{\text{code}}-\bar{b}\\ -\dfrac{c}{f_{0}}D-\dfrac{d}{dt}\Vert \mathring{\ell}-\rho\Vert_{2} \end{bmatrix}\right)\right\Vert_{2}^{2}. \end{align*}

5. Implementation and evaluation

We have implemented all the methods that we described above in MATLAB.

We use Borre's Easy Suite (Borre, Reference Borre2003, Reference Borre2009) to perform many routine calculations. In particular, we use it to correct the GPS time (check_t), to correct for Earth rotation during signal propagation time (e_r_corr), to read an ephemeris from a RINEX file and to extract the data for a particular satellite (rinexe, get_eph and find_eph), to transform Julian dates to GPS time (gps_time), to represent Julian dates as one number (julday), to compute the coordinate of a satellite at a given time in ECEF coordinates (satpos), to compute the azimuth, elevation and distance to a satellite (topocent, which calls togeod to transform ECEF to WG84 coordinates), and to approximate the tropospheric delay (tropo). We also use a MATLAB function by Eric Ogier (ionophericDelay.m, available on the MathWorks File Exchange) to approximate the ionopheric delay using the Klobuchar model. We take the parameters for the Klobuchar model from files published by the GNSS Research Center at Curtin University.Footnote ²

During the Gauss–Newton phase of the algorithm (after the integers have been determined), if we have only four observations, we add a pseudo-measurement constraint that constrains the correction to maintain the height of the target, in the least-squares sense (Van Diggelen, Reference Van Diggelen2009).

We use Chang and Zhou's MILES package (Chang and Zhou, Reference Chang and Zhou2007) to solve mixed-integer least-squares problems.

We take ephemeris data from RINEX navigation files published by NASA.Footnote ³

We filter satellites that are lower than 10 degrees above the horizon which have the lowest SNR and are more likely to suffer from multipath interference.

We evaluated the code on data from several sources:

• • Publicly available observation data files in a standard format (RINEX) distributed by NASA. We used these to test our algorithms in the initial phases of the research. These results are not shown here.

• • GPS simulations. We generated satellite positions ephemeris files and used them to compute times of arrival and code phases. These simulations do not include ionospheric or tropospheric delays, so they help us separate the issues arising from these delays from other algorithmic issues.

• • Code-phase and Doppler-shift measurements collected by us using a u-blox ZED-F9P GNSS receiver, connected to an ANN-MB-00 u-blox antenna mounted on a steel plate on top of a roof with excellent sky view. We established the precise coordinates of the antenna (to compute errors) using differential carrier-phase corrections from a commercial virtual reference station (VRS).Footnote ⁴ The WGS84 coordinates of the antenna are 32$\cdot$1121756, 34$\cdot$8055775 with height above sea level of 61$\cdot$15 m. The code phase measurements are included in the UBX-RXM-MEASX emitted by the receiver. The data set includes approximately 700 epochs, one every minute (so they span a little more than 11 h). The number of satellites per epoch ranges from 8 to 13 and after filtering by elevation, between 7 and 11.

• • Recordings of RF samples made by a bat-tracking GPS snapshot logger. The tag model is called Vesper. It was designed and produced by Alex Schwartz Developments on the basis of an earlier tag called Robin, designed and produced by a company called CellGuide that no longer exists. The tag records 1-bit RF samples at a rate of 1,023,000 samples per second. (The sampling rate is a multiple of $1/t_{\text {code}}$; this is known to make time-of-arrival estimate difficult (Tran et al., Reference Tran, Shivaramaiah, Nguyen, Glennon and Dempster2018) but the rate cannot be changed in this logger.) The tag was configured to record a 256 ms sample every 10 min for a few hours. It was placed next to the ANN-MB-00 antenna.

Figure 1 shows the cumulative distribution function (CDF) of four algorithms: Van Diggelen's non-iterative method, the Doppler constraints alone (as used in the first phase of Fernández-Hernández and Borre's method), and MILS with either a priori or Doppler regularisation. The data from the u-blox receiver were used to produce these graphs. We used all 694 epochs. The initial error was of 20–21 s (uniform distribution) and 20 km in a random uniform horizontal direction. In the MILS algorithms, the final position was computed with the regularisation constraints; this is why the Doppler regularisation produced less accurate results. We can see that the accuracy of MILS with a priori regularisation and of Van Diggelen's method is essentially identical.

Figure 1. Cumulative distribution function of the absolute positioning errors of four algorithms: Van Diggelen's non-iterative method, the Doppler constraints alone (the first phase of Fernández-Hernández and Borre's method), and mixed-integer least-squares (MILS) with either a priori or Doppler regularisation

Figure 2 compares the probability of success achieved by our regularised MILS solver with that achieved by Van Diggelen's non-iterative method and by the Doppler constraints alone. We considered fixes that are within 1 km of the true location to be a success in obtaining the correct integer values. Each pixel in these heat maps represents 16 different runs. Each run uses a random epoch, a random initial location estimate and a random initial time estimate. The initial location estimates have a given distance to the true location (the $x$ axis of the heat map) but a random azimuth. The initial time estimate is a slight perturbation (uniform between zero and one second) of the given time error, which is the $y$ axis of the heat map. In each pixel, half of the initial time errors are positive and half are negative. We used all the satellites in view in each epoch.

Figure 2. The probability of obtaining a fix with an error smaller than 1 km from the u-blox data set using four different algorithms

The results clearly show that the MILS algorithm, even with the simple a priori time and location regularisation from Section 4.1, outperforms Van Diggelen's non-iterative method. Van Diggelen's method obtains a correct fix in almost all cases (success probability close to 1) when the initial location error is small and the initial time error is 150 s or less, when the initial location time error is small and the initial location error is 100 km or less, and in other equivalent combinations of time and location errors. The corresponding limits for the MILS algorithm with a priori regularisation are approximately 150 s and 250 km.

Doppler-shift observations expand dramatically the region of convergence in both approaches. The MILS algorithm with Doppler regularisation obtains a correct fix as long as the initial time error is at up to approximately 180 s (3 min); this works even with great-circle distances of 20,000 km, which means that the initial position can essentially be anywhere on Earth. If the initial position error is small, the method can tolerate initial time errors of up to approximately 80 min (5000 s). The heat map of the Doppler constraints alone, together with the CDF in Figure 1, indicate that these constraints produce an estimate good enough for initialising Van Diggelen's method, but are not accurate enough on their own. Indeed, Van Diggelen writes about the Doppler constraints alone: ‘For less than 1 Hz of measurement error, we expect a position error of the order of 1 km’ (Van Diggelen, Reference Van Diggelen2009, Section 8.3); this explains why the probabilities in the top-right plot in Figure 2 are usually far from 1, even with small initial errors. In general, both approaches have similar regions of convergence and they produce similarly accurate fixes.

Figure 3 explores how the number of satellites (observations) affects the success rates of the four methods. We repeated the experiment whose results are shown in the heat maps in Figure 2, but only on the 20 epochs in which 13 satellites were in view. We selected random subsets of the satellites in view and random initial errors, within the bounds shown in Figure 2, and computed the fraction of successful experiments. We can see that when only code phases are used, Van Diggelen's method is better when using 6–8 observations, probably because the weighting of the observations in the MILS method sometimes leads to incorrect integers when Van Diggelen's method resolves the integers correctly. However, with 5 satellites in view or more than 8, the MILS method is better. MILS with Doppler regularisation is superior to all the other methods.

Figure 3. The fraction of successful positioning (error of at most 1 km) in the spaces of initial errors shown in Figure 2 as a function of the number of satellites (observations) used

While our Matlab implementation is not designed to carefully evaluate running times and computational efficiency, we did measure the running times and we can draw from them some useful conclusions. The running times of a single Gauss–Newton correction step in Van Diggelen's method and in the solution of the Doppler equations is 8–10 $\mathrm {\mu }\text {s}$, while the running time of a single Gauss–Newton correction step in the MILS formulation is approximately 2$\cdot$5 ms when using a priori regularisation and 0$\cdot$6 ms when using Doppler regularisation. While the MILS methods are clearly more expensive, they also appear to be fast enough for real-time applications.