Estimating Position from Millisecond Samples of GPS Signals (the “FastFix” Algorithm)

Abstract

A new approach to GPS positioning is described in which the post-processing of ultra-short sequences of captured GPS signal data can produce an estimate of receiver location. The algorithm, called ‘FastFix’, needs only 2–4 ms of stored L1-band data sampled at ∼16 MHz. The algorithm uses a least-squares optimization to estimate receiver position and GPS time from measurements of the relative codephase, and Doppler-shift of GNSS satellite signals. A practical application of this algorithm is demonstrated in a small, lightweight, low-power tracking tag that periodically wakes-up, records and stores 4 ms of GPS L1-band signal and returns to a low-power state—reducing power requirements by a factor of ∼10,000 compared to typical GPS devices. Stationary device testing shows a median error of 27.7 m with a small patch antenna. Results from deployment of this tag on adult Royal Albatross show excellent performance, demonstrating lightweight, solar-powered, long-term tracking of these remarkable birds. This work was performed on the GPS system; however, the algorithm is applicable to other GNSS systems.

1. Introduction

Tracking devices based on GNSS (Global Navigation Satellite System) are now commonplace, their many uses include wildlife tracking [1] and vehicle tracking [2]. Wildlife applications, in particular, are constrained by size and weight, and, until recently, GNSS devices have been limited to tracking animals such as moose [3] and marine species [4] capable of carrying a large receiver, limiting the range of species that can be tagged [5]. Improvements in semiconductor fabrication have steadily reduced power consumption [6], and improvement in battery technology [7] have reduced the weight of such tracking devices to the point that GNSS tracking devices can be deployed to track small mammals and medium-sized birds [8,9]. With each technical development allowing smaller and lower-power tracking devices, the range of species that can be tracked, and the duration of the deployment is increased.

The most widely used GNSS is the GPS (global Positioning System) network [10]. GNSS tracking devices operate by using codephase and navigation message data transmitted from the satellites to determine the time of transmission of a signal. This can be used to determine the distance (ranges) from the current position to each satellite as well as the precise local time (see for example Kaplan [10]).

Navigation messages transmitted from the receiver contain sufficient information to determine orbital parameters (ephemerides) of the satellites and the time of propagation of the signals from the satellite, and hence the ranges. Using the measured ranges from four or more satellites, the position of the receiver is determined—commonly using a least-squares procedure [11]. The time to first fix typically exceeds 30 s, as this is the time required for all of the satellite ephemerides to be transmitted. As a result, traditional GPS receivers are required to operate from cold-start for more than 30 s—typically consuming tens of milliamps—in order to measure their position. GNSS position accuracy varies, and in the case of GPS, the distribution of positioning errors is not Gaussian [12] and depends sensitively on the antenna characteristics [10]; typically, the first fix from a single-frequency commercial GPS has a horizontal positioning error of ∼10 m [13] in open skies.

Alternative Algorithms

Advances in algorithms provide another technique to reduce power consumption, and several alternative approaches to the traditional positioning algorithm have been developed. Of these, the most widely used is Assisted GPS (A-GPS). This employs satellite ephemeris data, and estimates of position, time to eliminate the need to receive this data via the GNSS signal [14]. This has limitations in tracking applications as a separate telemetry channel is needed for the augmentation data. Time-to-first-fix from a cold start situation, signal sensitivity, and position-accuracy are all improved with respect to a stand-alone GPS receiver [15].

Another algorithm family estimate position using only millisecond samples of data, with the satellite ephemeris provided externally. This technique is referred to as time-free, or snapshot positioning and was first proposed by Peterson et al. in 1995 [16]. This was refined in 2009 [17] to show that positioning was possible using only codephase measurements. Van Dierendonck et al. [18] was proposed in 2015 that uses an initial approximation of receiver time and position, to reconstruct the transmit time relative to the received time such that pseudorange can be calculated, and position estimated. In 2012, Othieno and Gleason [19] introduced a positioning method that use a combination of Doppler and code-phase measurements to provide a position estimate unknown initial position. Jing et al. [20] in 2017 introduced an integer ambiguity searching method based on the relationship between the maximum pseudorange difference and the inter-satellite distance. The accuracy possible with these snapshot algorithms is typically ∼20 m [21].

This algorithm described in this paper uses Doppler information to provide an a priori position estimate, in a manner similar to Othieno [19], and combines a relative-codephase measure, in a manner analogous to the work of Jing et al. [20] who instead use an inter-satellite distance. In practice, the algorithm is able to determine position anywhere on earth from as little as 2 ms of GPS L1-band signal using measurements of Doppler shift and relative codephase recovered from the data. Section 2 describes the well-known process of acquisition, and shows that small errors in the knowledge of the local time at the receiver will affect the measurement of codephase; however, this effect is the same for each satellite signal, and introduces relative codephase which is not affected by these errors. Section 3 describes the FastFix algorithm for positioning, and Section 4 presents results for stationary device performance including typical accuracy. Finally, Section 5 presents a real-world application of this algorithm in tracking albatross.

2. Satellite Signal Acquisition

The GPS L1-band civilian signal is transmitted in a 2.5 MHz bandwidth radio channel centered at 1575.42 MHz. The satellite signals are modulated using Code Division Multiple Access (CDMA) [10]. A GPS receiver digitizes this channel, and signals from all visible satellites are contained within the received data. The individual satellite signals are recovered by correlating the received data against a pseudorandom sequence unique to each satellite called a C/A (Coarse Acquisition) code. In the case of GPS, the C/A code repeats every millisecond and contains 1023 binary numbers in a pseudorandom sequence. The motion of each satellite relative to the observer causes the signal from that satellite to be Doppler-shifted, and this Doppler-shift must be compensated for in the receiver. The process of signal acquisition recovers the signal, and its Doppler-shift from each visible satellite.

There are several well-established methods for signal acquisition; see, for example, Kaplan [10] or Borre et al. [22]. In this work, parallel codephase recovery [22] using Fourier transforms was used. For the ith satellite, the acquisition process results in three measurements:

• $X_{max}/{\sigma }_X$. The maximum correlation divided by the standard-deviation of the correlation—a measure of the ratio of correlation amplitude to the noise. This is similar to the carrier-to-noise density ratio, or the correlation amplitude relative to the noise level in the data. A satellite is considered visible if its $X_{max}/{\sigma }_X$ exceeds a threshold (for the data shown in Figure 1 and

  Table 1. Typical acquisition results for a 4 ms sample of signal (shown in Figure 1). There are six satellite signals that exceed the threshold set for signal detection in this data.

  i	SV	${\mathit{X}}_{\mathit{max}}/{\mathit{\sigma }}_{\mathit{X}}$	Doppler-Shift $\mathsf{\Delta }{\mathit{f}}^{\mathit{i}}$	Codephase ${\mathit{\delta }}_{\mathit{\varphi }}^{\mathit{i}}$
  1	16	30.01	+2185.15461	0.503051794
  2	23	21.08	+1632.13669	0.576269412
  3	32	16.84	−2624.85478	0.209007485
  4	13	14.03	+2335.25122	0.197458456
  5	20	13.56	−818.94480	0.344973267
  6	4	10.70	−155.43040	0.439863498

  , a threshold of 7 is used).
• ${\delta }_{\varphi }^i$, the code-phase—a number between zero and 1 that indicates the phase of the C/A code. The starting time of the C/A code is given by $t_0+{\delta }_{\varphi }^i\tau$, where $t_0$ is the time of the start of the data, and $\tau$ is the C/A code sequence repeat time (one millisecond in the case of GPS).
• $\Delta f_m^i$, the Doppler-shift of the signal from the satellite.

Table 1 shows a typical set of acquisition results from the data, which are shown in Figure 1a. Some of the correlation plots from this data are shown in Figure 1b–d.

2.1. System and Local Time Frames

The acquisition process uses data that starts at time $t_0$, measured by the receiver clock that may be inaccurate. Times measured by the receiver are identified by lower case t. The GNSS system clock times are identified by upper-case time values and are typically accurate to better than 1 $\mathsf{\mu }$s [10].

All GNSS satellites synchronously send their signals at intervals of $\tau$. The start of each code epoch is synchronous across all satellites in the constellation. The system time that the ith satellite signal is sent, $T_s^i$, is given by $T_s^i=T_s+M\tau$, where $T_s$ is an arbitrary start time for the GPS constellation of satellites, and M is the integer code epoch. We choose, quite arbitrarily, that $T_s=0$. Therefore, the time of transmission for the start of the code $T_s^i$ becomes

$$T_s^i=M^i\tau ,$$

(1)

where the $M^i$ are integers that are now different for each satellite, and depend on the relative distance between receiver and each satellite.

We assume that there is a fixed offset, $t_u$, between the local receiver clock and the system time, i.e.,

$$T=t+t_u.$$

(2)

A more sophisticated model of the receiver clock would have both an offset and a different rate, but as the data are sampled for only a few milliseconds, uncertainty in the local clock rate will not significantly affect this algorithm.

2.2. Signal Propagation Time

Let $T_r^i$ be the system time that the start of the ith C/A code reaches the receiver. This is given by

$$T_r^i=t^i+t_u,$$

where $t^i\equiv t_0+{\delta }_{\varphi }^i\tau$ is the measured time for the start of the sequence at the receiver, and $t_u$ is the receiver clock offset. The system time at the receiver for the start of the C/A code from satellite i is given in terms of $t_0$ and the code-phase as

$$T_r^i=t_0+{\delta }_{\varphi }^i\tau +t_u.$$

(3)

At system time $T_r^i$, the signal must have originated from the satellite at an earlier time $T_s^i=T_r^i-\frac{{\Delta }_r^i}{c}$ due to the propagation velocity c of the signal over the distance, ${\Delta }_r^i$, between receiver and satellite i. We know that $T_s^i$ is a time when an epoch starts, so it is an integer multiple of $\tau$, i.e., $M^i\tau$, so, in this case, we write

$$\begin{array}{ccc}\hfill M^i\tau & =& T_r^i-\frac{{\Delta }_r^i}{c},\hfill \\ & =& t_0+{\delta }_{\varphi }^i\tau +t_u-\frac{{\Delta }_r^i}{c},\hfill \\ \hfill {\delta }_{\varphi }^i\tau & =& M^i\tau +\frac{{\Delta }_r^i}{c}-(t_0+t_u).\hfill \end{array}$$

(4)

The codephase ${\delta }_{\varphi }^i$ changes rapidly with time due to satellite motion, i.e, $\frac{d{\delta }_{\varphi }^i}{dt}=\frac{{\Delta }_r^i}{c\tau }\sim 67$ cycles per second. In addition, uncertainty in the system time offset, of more than a fraction of a millisecond will make prediction of position from the codephase unfeasible. For example, an uncertainty in $t_u$ of 100 μs would lead to an uncertainty in ${\delta }_{\varphi }^i$ of 0.1, which would translate to a pseudorange error of ∼30 km.

2.3. Relative Fractional Codephase

The relative fractional codephase, ${\delta }_{\varphi }^{ij}$ is the difference (modulo $\tau$) between the fractional codephases ${\delta }_{\varphi }^i$ and ${\delta }_{\varphi }^j$ from satellite i and satellite j. From Equation (4),

$$\begin{array}{ccc}\hfill {\delta }_{\varphi }^i\tau -{\delta }_{\varphi }^j\tau & =& (M^i-M^j)\tau +\left(\frac{{\Delta }_r^i}{c}-\frac{{\Delta }_r^j}{c}\right)\hfill \\ \hfill {\delta }_{\varphi }^i-{\delta }_{\varphi }^j& =& \left[(M^i-M^j)+\left(\frac{{\Delta }_r^i-{\Delta }_r^j}{c\tau }\right)\right]\hfill \\ \hfill {\delta }_{\varphi }^i-{\delta }_{\varphi }^j& =& N^{ij}+{\delta }_{\varphi }^{ij},\hfill \end{array}$$

where $N^{ij}=M^i-M^j$ is the integer part of the expression, and the relative codephase ${\delta }_{\varphi }^{ij}$ is the fractional part of the expression. For example,

$$\begin{array}{ccc}\hfill {\delta }_{\varphi }^{ij}(t,{\mathbf{r}}_0)& =& \mathrm{frac}\left(\frac{{\Delta }_r^i-{\Delta }_r^j}{c\tau }\right).\hfill \end{array}$$

(5)

This relative fractional-codephase does not depend sensitively on accurate knowledge of the local time offset $t_u$; however, range errors due to incorrect time estimates will have an effect. Errors of 50 km or more will still lead to relatively small changes in ${\delta }_{\varphi }^{ij}$. At typical satellite velocities of ∼3 km s${}^{-1}$, an uncertainty of ten seconds in $t_u$ would lead to a typical ${\delta }_{\varphi }^{ij}(t,{\mathbf{r}}_0)$ uncertainty of approximately 0.1. Contrast this with the codephase ${\delta }_{\varphi }^i$, for which this level of uncertainty would be caused by a ∼100 μs uncertainty in $t_u$.

3. The FastFix Positioning Algorithm

The FastFix algorithm uses both the Doppler-shift and the codephase measurements to solve the inverse problem of estimating the position, ${\mathbf{r}}_0$, and time-of-fix, $t_0$ from stored GPS $L_1$-band data. The forward problem assumes that the position $\mathbf{r}$, local time $t_0$ and local clock offset $t_u$ are known, and satellite ephemerides are available. This is sufficient to determine from Equation (5) the relative codephase, ${\delta }_{\varphi }^{ij}$, and the Doppler-shift of the signal from the ith satellite,

$$\Delta f^i(t_0,\mathbf{r})=-\frac{1}{c}\frac{d{\Delta }_r^i}{dt}{f}_0,$$

(6)

where $f_0\sim 1575.42$ MHz, c is the speed of light, and ${\Delta }_r^i$ is the relative distance between the ith satellite and the receiver position $\mathbf{r}$. If ${\Delta }_r^i$ is decreasing, then satellite is moving towards the receiver, and the Doppler-shift is positive.

This method is rather like the method of exact fractions from interferometry (see, for example, Baird [23]), the measurements comprise only the fractional part of the codephase for the signal from each satellite and its Doppler shift. At any given moment, there are only a small number of positions near the earth’s surface that will have any specific set of relative-codephases.

3.1. Least-Squares Estimators

The algorithm uses two least-squares estimators: one based on Doppler measurements, the other based on the relative fractional codephase. The least-squares estimator for Doppler shift, $L_{\nu }$, is

$$L_{\nu }(t,{\mathbf{r}}_0)=\sum _i{(\Delta f_m^i-\Delta f^i(t,{\mathbf{r}}_0))}^2,$$

(7)

where $\Delta f_m^i$ is the measured Doppler-shift for the ith satellite, and $\Delta f^i(t,{\dot{\mathbf{r}}}_0)$ is the predicted Doppler-shift for the ith satellite from Equation (6). The least-squares measure for relative fractional codephase, $L_{\delta \varphi }$, is

$$L_{\delta \varphi }(t,{\mathbf{r}}_0)=\sum _{i,j}{({\delta }_{\varphi m}^{ij}-{\delta }_{\varphi }^{ij}(t,{\mathbf{r}}_0))}^2,$$

(8)

where ${\delta }_{\varphi m}^{ij}$ is the measured relative fractional codephase between the ith and jth satellite, and ${\delta }_{\varphi }^{ij}(t,{\mathbf{r}}_0)$ is the predicted relative fractional codephase given by Equation (5).

3.2. Algorithm Outline

The algorithm begins with an estimate of the local time $t_0$, and, by acquiring the codephase, Doppler-shift and $X_{max}/{\sigma }_X$ for each visible satellite from the previously stored GPS $L_1$-data. These measurements are then used for three successive least-squares estimations of position. Starting with a global minimization based on Doppler-shift measurement, the successive stages each further refine the position from the preceding stage. A flow-chart for the algorithm is shown in Figure 2. Each of the minimization stages requires knowledge to the satellite orbital elements, and it is assumed that these are downloaded separately from a suitable catalog.

The sensitivity of the $\Delta f^i$ to the solution is quite low, so, from an initial measurement of $t_0$, an initial first pass is made searching for a position that minimizes $L_{\nu }$ from Equation (7). This results in a rough position estimate ${\mathbf{r}}_1$.

This Doppler-estimate is refined by searching in the region near ${\mathbf{r}}_1$ for a position ${\mathbf{r}}_2$ that optimizes the relative fractional codephase estimator, $L_{\delta \varphi }$ at time $t_0$.

The final step is a full four-dimensional space-time optimisation using a Nelder–Meade [24] optimisation algorithm. This yields an estimate for the receiver position, ${\mathbf{r}}_3$, as well as the local time $t_1$.

4. Results

The FastFix algorithm was initially tested using stationary receivers in a known locations (Dunedin, New Zealand). In this test, 118 four-millisecond samples of GPS L1-band data were digitized and stored. From these, the FastFix algorithm (described in the previous section) was applied.

4.1. Doppler Search

Figure 3 shows contours of the Doppler least-squares estimator $L_{\nu }$ plotted over the globe centered on the receiver location. This shows a global and local minimum at the receiver location. This estimate is not accurate, however it is sufficient to locate the receiver within a small enough region for the next step which is the minimisation of the relative codephase estimator $L_{\delta \varphi }$. Figure 4 shows a histogram of the Doppler-only location error, showing a median error of 175 km and a 95th percentile error of 537 km.

4.2. Space Search

Figure 5 shows the relative codephase estimator $L_{\delta \varphi }$ plotted for six visible satellites over a region ten degrees in width and height centred at the actual receiver location. This figure shows that the search for the global minimum is a complex nonlinear minimisation problem. However, the Doppler search is sufficiently accurate to place the receiver within the central local minimum region and the relative codephase least-squares estimator can be easily minimised once this is known.

4.3. Measured Accuracy

Figure 6 shows a histogram of the fix error for a known location. This shows that the relative codephase least-squares minimisation produces accurate fixes most of the time. The median error is 27.7 m, and all results are based on 4 ms of sampled data from a GPS front end (Maxim MAX2769B Universal GPS Receiver).

There are outliers that are caused by convergence on to a local, but not global minimum. Work is continuing on improving the performance of this; however, the results are already good enough to produce a device based on this technique. Section 4.4 outlines some techniques that may improve the robustness and accuracy of the results. This positioning error is similar to other snapshot positioning algorithms [21].

4.4. Discussion and Potential Improvements

One avenue for improving robustness is establishing better criteria for rejecting poor position estimates. Cook’s distance is a commonly used estimate of the influence of a data point when doing least-squares regression. Cook’s distance measures the effect of deleting a given observation. Another technique is to use Huber’s M-estimation [25] to handle outliers due to bad measurements.

Time-Separated Measurements

Another possible improvement would be to combine two sets of relative codephase and Doppler measurements taken from two samples of L1-band data a short time, $\Gamma$, apart. Assuming $\Gamma$ is large enough for the satellites to have moved by ∼1 km, yet small enough for the receiver to have moved by less than the desired spatial resolution, then the following relative codephase estimator should be more robust,

$$\begin{array}{ccc}\hfill L_{\delta \varphi }(t,{\mathbf{r}}_0)& =& \sum _i{\left({\delta }_{\varphi m}^{ij}-{\delta }_{\varphi }^{ij}(t,{\mathbf{r}}_0)\right)}^2+\hfill \\ & & \sum _i{\left(({\delta }_{\varphi m}^{ij}-{\delta }_{\varphi }^{ij}(t+\Gamma ,{\mathbf{r}}_0)\right)}^2.\hfill \end{array}$$

Taking multiple samples would, however, increase the power required.

5. Application to Wildlife Tracking

There are many possible applications for devices that use this algorithm. As only very short snapshots of GPS data are required, a device only needs to operate for 4 ms to acquire the data necessary for a fix. This dramatically reduces the power required (from at least 30 s of operation for a traditional GPS fix), by a factor of approximately 10,000.

A disadvantage of tracking with this technique is that any tracking device has to be recovered in order for the position to be estimated by post-processing the stored data. This technique is suitable for tracking species that return to a known location so that the data can be recovered. Here, we describe how this was done on nesting albatross.

5.1. Tag Electronics

Tags were developed that include a standard GPS RF front-end (Maxim MAX2769B Universal GPS Receiver) connected to a high-performance low-power 32-bit STM32F103CB [26] microprocessor. The microprocessor is programmed to wake at intervals, read data from the GPS RF front-end and store this data in a file on an micro Secure Digital (microSD) Memory Card. The un-packaged electronics from one of these FastFix tags are shown in Figure 7. A packaged solar-powered version of the tag, designed to attach to seabird feathers was developed [9] for deployment on seabirds. This tracking tag uses a small Taoglas AP.12F.07.0045 active patch antenna [27].

Each fix stores 4 ms of 1-bit baseband data sampled at 16.368 MHz, which requires 64 kb of data. A 2 GB SD-card can store approximately 244,000 fixes. A solar-powered tag could operate continuously storing data for a position estimate every five minutes for more than two years, limited only by storage limits.

5.2. Wildlife Tracking Results

Twelve solar-powered FastFix tags were attached to nesting adult Royal Albatross [28] (shown packaged and attached in Figure 8). While nesting, the adults carry out feeding flights of up to two weeks of duration. Of the 12 tags attached, eight were recovered undamaged, and from these eight tags, 26,042 positions were calculated. The position tracks recovered from the tags are shown in Figure 9.

6. Conclusions

The algorithm described here allows ultra-low-power, reasonably accurate global position determination using very short sequences of GPS data. The computational effort required to measure position is a function of both the initial uncertainty in position and the local receiver clock offset. The results in stationary testing show a median error of 27.7 m which is worse than the median error of traditional GPS single-frequency receivers, but similar to existing snapshot algorithms [21]. I speculate that this is due to the lack of a code-tracking filter that can refine the pseudorange estimates during the relatively long period that a traditional GPS device spends receiving navigation information. Despite this limitation, there are several interesting avenues to explore for improving the robustness and accuracy which were discussed in Section 4.4. The performance of the algorithm in practical moderate-precision applications such as wildlife tracking is sufficient as demonstrated in Section 5, and the power consumption was extremely low, which enabling long-term solar-powered, or primary-cell tracking of wildlife.

Another avenue for further work would be to use Bayesian inference techniques to estimate a posterior distribution of receiver positions. This would also enable uncertainty quantification and, consequently, a systematic exploration of the lower limit on the amount of stored data necessary to recover position to a given accuracy. An additional advantage of a Bayesian approach would be avoiding the decision to reject or accept a satellite as ‘visible’ from its correlation strength $X_{max}/{\sigma }_X$.