A Mixed Integer Least-Squares Formulation of the GNSS Snapshot Positioning Problem
aa r X i v : . [ ee ss . SP ] J a n A Mixed Integer Least-Squares Formulation of the GNSSSnapshot Positioning Problem
Eyal Waserman and Sivan ToledoBlavatnik School of Computer Science, Tel-Aviv UniversityJanuary 12, 2021
Abstract
This paper presents a formulation of
Snapshot Positioning as a mixed-integer least-squares prob-lem. In snapshot positioning, one estimates a position from code-phase (and possibly Doppler-shift)observations of Global Navigation Satellite Systems (GNSS) signals without knowing the time of de-parture (time stamp) of the codes. Solving the problem allows a receiver to determine a fix from shortradio-frequency snapshots missing the time-stamp information embedded in the GNSS data stream.This is used to reduce the time to first fix in some receivers, and it is used in certain wildlife trackers.This paper presents two new formulations of the problem and an algorithm that solves the resultingmixed-integer least-squares problems. We also show that the new formulations can produce fixeseven with huge initial errors, much larger than permitted in Van Diggelen’s widely-cited coarse-timenavigation method.
Keywords: GNSS, Doppler, Integer Ambiguity Resolution, Satellite Navigation
The fundamental observation equation in Global Navigation Satellite Systems (GNSS) is t i − t D,i = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + δ i + ǫ i , where t i is the observed (estimated) time of arrival of a code from satellite i , t D,i is the time of departureof the signal, c is the speed of light, ˚ ℓ is the location of the receiver, ρ i ( t D,i ) is the location of the satelliteat the time of transmission, ˚ b is the offset in time-of-arrival observation caused by the inaccurate clockat the receiver and by delays in the analog RF chain (e.g., in cables), δ i represent atmospheric delaysand the satellite’s clock error, and ǫ i is an error or noise term that accounts for both physical noiseand for unmodeled effects. Normally, the GNSS solver estimates ˚ ℓ and ˚ b by minimizing the norm ofthe error vector ǫ ; it does not know ǫ i and does not attempt to estimate it. The quantities t D,i and ρ i ( t D,i ) are usually known; t D,i because the satellite time stamps its transmission, and ρ i ( t D,i ) becausethe satellite transmits the parameters that define its orbit, called the ephemeris . The ephemeris can alsobe downloaded from the Internet.Decoding t D,i takes a significant amount of time, in GPS up to 6 seconds under good SNR conditionsand longer in low-SNR conditions. GNSS receivers that need to log locations by observing the RF signalsfor short periods cannot decode t D,i . Examples for such applications include tracking marine animalslike sea turtles, which surface briefly and then submerge again. It turns out that techniques that arecollectively called snapshot positioning or coarse-time navigation can estimate ˚ ℓ and ˚ b when t D,i is notknown. These techniques can also reduce the time to a first fix when a receiver is turned on.Snapshot receivers sample the incoming GNSS RF signals for a short period, called a snapshot . Usually(but not always), the RF samples are correlated with replicas of the codes transmitted by the satellites,therefore determining t i for the subset of visible satellites. The correlation (and Doppler search) aresometimes performed on the receiver, who then stores or transmits the t i s. This appears to be the casefor a proprietary technology called Fastloc, which is used primarily to track marine animals (Dujon,Lindstrom, and Hays 2014; Tomkiewicz et al. 2010; Witt et al. 2010). In other cases (Bavaro 2012b;Bavaro 2012a; Cvikel et al. 2015; Eichelberger, Hagen, and Wattenhofer 2019; Harten et al. 2020; Liuet al. 2012; Ramos et al. 2011), the logger records the raw RF samples and correlation is performed afterthe data is uploaded to a computer. 1echniques for estimating ˚ ℓ when the t D,i s are not known date back to a 1995 paper by Peterson etal (Peterson, Hartnett, and Ottman 1995). They proposed to view t D,i as a function of both t i and acoarse clock-error unknown that they call coarse time , which in principle is identical to ˚ b , but is modeledby a separate variable. They then show that it is usually possible to estimate ˚ ℓ , ˚ b , and the coarsetime from five or more t i s. This method does not always resolve the t D,i correctly. Lannelongue andPablos (Lannelongue and Pablos 1998) and Van Diggelen (Van Diggelen 2002; Diggelen 2009) proposedmethods that appear to always resolve the t D,i correctly when the initial estimate of ˚ ℓ and ˚ b are withinsome limits (adding up to about 150 km). Muthuraman, Brown, and Chansarkar (Muthuraman, Brown,and Chansarkar 2012) showed that the two methods are equivalent in the sense that they usually producethe same estimates. However, the method of Lannelongue and Pablos is an iterative search procedure,while Van Diggelen’s is a rounding procedure that is more computationally efficient, so Van Diggelen’smethod became much more widely used and widely cited.All three methods use a system of linearized equations with five scalar unknowns (not the usual 4),which are corrections to the coordinates of ˚ ℓ , the offset ˚ b (which he refers to as the common bias), andthe coarse time unknown, usually denoted by tc (we denote it in this paper by the single letter s ).Van Diggelen’s algorithm was used and cited in numerous subsequent papers, all of which repeat hispresentation without adding explanations. Liu et al (Liu et al. 2012), Ramos et al (Ramos et al. 2011),and Wang et al (Wang, Qin, and Jin 2019) describe snapshot GPS loggers whose recordings are processedusing the algorithm. Badia-Solé and Iacobescu Ioan (Badia-Solé and Iacobescu Ioan 2010) report on theperformance of the method. Othieno and Gleason (Othieno and Gleason 2012), Chen at al (W. et al.2014), and Fernández-Hernández and Borre (Fernández-Hernández and Borre 2016) show how to useDoppler measurements to obtain an initial estimate that satisfies the requirements of Van Diggelen’salgorithm. Yoo et al (Yoo et al. 2020) propose a technique that replaces the estimation of the coarse-timeparameter by a one-dimensional grid search, which allows them to estimate ˚ ℓ using observations fromonly four satellites, but at considerable computational expense.Bissig et al (Bissig, Eichelberger, and Wattenhofer 2017) use a direct position determination (Weiss2004) approach to snapshot positioning. They quantize the four unknowns ˚ ℓ and ˚ b and maximize thelikelihood of the received snapshot over this four-dimensional integer space. To make the search efficient,they use a branch and bound approach that prunes sets of unlikely solutions. It appears that this approachallows them to estimate locations using very short snapshots, but at the cost of fairly low accuracy andrelatively long running times, compared to methods that do estimate the t i s first.In this paper we derive the observation equations that underlie the methods of Peterson et al. (Pe-terson, Hartnett, and Ottman 1995), Lannelongue and Pablos (Lannelongue and Pablos 1998) and VanDiggelen (Van Diggelen 2002; Diggelen 2009). These authors only show the correction equations, notthe observation equations whose Jacobian constitutes the correction equations. The formulation of theobservation equations, which constitute a mixed-integer least-squares problem, allows us to apply a newtype of algorithm to estimate the integer unknowns. More specifically, we regularize the mixed-integerproblem using either a priori estimates of ˚ ℓ and ˚ b or Doppler-shift observations.Our experimental results using real-world data demonstrate that our new algorithms can resolvelocations with much larger initial location and time errors than the method of Van Diggelen. VanDiggelen’s method only works when the initial estimate is up to about 150 s or 100 km (or equivalentcombinations), whereas our mixed-integer least-squares solver works with initial errors of up 150 s and200 km. When using Doppler-shift regularization, our method works even with initial errors of 1800 sand arbitrarily large initial position errors (on Earth); if the initial position error is small, the methodtolerates initial time errors of up to 7200 s.Our implementation of the new methods and the code that we used to evaluate them are publiclyavailable .The rest of this paper is organized as follows. Section 2 presents the observation equations forthe snapshot-positioning GNSS problem. Section 3 explains how to incorporate the so-called course-time parameter into the observation equations and how Van Diggelen’s method exploits it. Section 4.1presents our first regularized formulation, which uses the initial guess to regularize the mixed-integer least-squares problem. Section 4.2 presents the Doppler-regularized formulation. Our experimental results arepresented in Section 5. Section 6 discusses our conclusions from this research. https://github.com/eyalw711/snapshot-positioning 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 relatethe 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 representsa whole millisecond (in the time base of the GPS system). We denote by ˚ t i the time of arrival of thatcode at the antenna of the receiver. We assume that the receiver estimates the arrival time of a thatcode as t i = ˚ t i +˚ b + ǫ i , where ˚ b represents bias that is caused by the inaccurate clock of the receiver andby delays that the signal experiences in the path from the antenna to the ADC and ǫ i is the arrival-timeestimation error. The bias ˚ b is time dependent, because of drift in the receiver’s clock, but over shortobservation periods this dependence is negligible, so we ignore it.The time of arrival is governed by the equation ˚ t i − t D,i = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) + δ i , where c is the speed of light, ˚ ℓ is the location of the receiver, ρ i is the location of the satellite (which is afunction of time, since the satellites are not stationary relative to Earth observers), and δ i represents theinaccuracy of the satellite’s clock and atmospheric delays. We assume that δ i can be modeled, for exampleusing models of ionospheric and tropospheric delays (dual frequency receivers can estimate the ionosphericdelay, but we assume a single-frequency receiver). Setting δ i to the satellite’s clock error correction fromthe ephemeris induces a location error of about 30m due to the atmospheric delays (Borre and Strang2012).A receiver that decodes the time stamp embedded in the GPS data stream can determine t D,i , whichleads to the following equation t i − t D,i = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + δ i + ǫ i , the conventional GNSS code-observation equation, in which the unknown parameters are ˚ b and thecoordinates of ˚ ℓ (we assume that δ i is modeled, possibly trivially δ i = 0 , but not estimated). To simplifythe notation, we ignore δ i for now and write t i − t D,i = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + ǫ i . We use observations from all the satellites such that all the t i s lie between two consecutive whole multiplesof t code (in GPS, two round milliseconds in the local clock). This allows us to express t i = ( N + ϕ i ) t code with for a common and easily computable N = ⌊ t i /t code ⌋ and for ϕ i ∈ [0 , . We denote N i = t D,i /t code and write ( N − N i + ϕ i ) t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + ǫ i . Since GNSS codes are aligned with t code ,N i ∈ Z . We denote n i = N − N i ∈ Z , so ( n i + ϕ i ) t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + ǫ i or ϕ i t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) − n i t code +˚ b + ǫ i . We now face two challenges. One is that we have m unknown parameters: three location coordinates, b , and the n i s, but only m constraints. We clearly need more constraints so that we can resolve the n i s.The other is that we have a set of nonlinear constraints with continuous real unknowns, the locationand ˚ b , and with integer unknowns, the n i . The strategy, as in other cases with this structure, is tofirst linearize the non-linear term, then to resolve the integer parameters, and to then substitute themand to solve the continuous least-squares problem (either the linearized system or the original non-linearsystem). We cannot linearize the non-linear term using a Taylor series because it is a function of both3eal unknowns and the of the integer unknowns n i . We cannot differentiate with respect to the integer n i . To address this difficulty, we approximate t D,i by approximating the range (distance) term in theequation t D,i = t i − c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) − ˚ b − ǫ i . For now, we denote the approximation of the range by d i ≈ c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) + ǫ i so ˆ t D,i = t i − d i − ˚ b . There are several ways to set d i , depending on our prior knowledge of ˚ ℓ and ˚ b . One option in the GPSsystem is to set it to about . ms. This limits the error in ˆ t D,i to about . ms for any Earth observer,and the error (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:0) ˆ t D,i (cid:1)(cid:13)(cid:13)(cid:13) − (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) to about m.We substitute ρ i (ˆ t D,i ) = ρ i ( t i − d i − ˚ b ) for ρ i ( t D,i ) , ϕ i t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:16) t i − d i − ˚ b (cid:17)(cid:13)(cid:13)(cid:13) − n i t code +˚ b + ǫ ( D ) i . (1)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 linearize around an a priori solution ¯ ℓ and ¯ b (usually ¯ b = 0 , otherwise we can simply shift the t i s), ϕ i t code = 1 c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − n i t code +˚ b + ǫ ( D,L ) i c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − n i t code + (cid:16) ˚ b − ¯ b (cid:17) + ¯ b + ǫ ( D,L ) i , (2)where J is the Jacobian of the Euclidean distances with respect to both the location of the receiver and tothe time of departure of the signal, with the derivatives evaluated at ¯ ℓ and at t i − d i − ¯ b . The superscript ( D, L ) on the error term indicates that it includes now also the linearization error.There are now several ways to resolve the n i . Peterson et al. (Peterson, Hartnett, and Ottman 1995) introduced a somewhat surprising modeling tech-nique, which we refer to as shadowing . The idea is to replace the unknown b by two separate unknownsthat represent essentially the same quantity, the original b and a shadow s . In principle, they should obeythe equation b = s , but the model treats s as a free parameter; the constraint b = s is dropped. In theliterature, s is called the coarse-time parameter (and is often represented by tc , which we find confusing).We express this technique by splitting b and s : ϕ i t code = 1 c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − n i t code +˚ b + ǫ ( D,L ) i c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓs − ¯ b (cid:21) − n i t code +˚ b + ǫ ( D,L ) i . 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) is very sensitive to small (nanosecondscale) perturbations in the additive ˚ b , but it is not highly sensitive to the ˚ b (now s ) that we multiply by J i, = ∂∂t D,i (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) . (3)4or example, in GPS the derivative is bounded by about m/s for any ¯ ℓ on Earth, so J i, < × − (versus for the additive ˚ b ). Therefore, the dependence of the residual on ˚ b in Equation (2) is highlynon-convex. There are many different values of ˚ b that are almost equally good, a millisecond apart, witheach of these nearly-optimal hypotheses being locally well defined.Shadowing turns this non convexity into explicit rank deficiency, which is easier to deal with. With oneinstance of ˚ b replaced by the shadow s , the constraints no longer uniquely define ˚ b , only up to a multipleof t code . For any hypothetical solution ℓ, s, b, n , the solution ℓ, s, b + kt code , n + k gives exactly the sameresidual. We perform a change of variables, replacing the partial sum − n i t code +˚ b by − ν i t code + β , where − ν i = − n i + j ˚ b/t code k and β = ˚ b − ⌊ b/t code ⌋ t code , so β ∈ [0 , t code ) : ϕ i t code = 1 c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓs − ¯ b (cid:21) − ν i t code + β + ǫ ( D,L ) i . (4) β ∈ [0 , t code ) . Van Diggelen’s method exploits the fact that the ν i ’s are very insensitive to ˚ ℓ and to s . It therefore sets ˚ ℓ = ¯ ℓ and s = ¯ b , truncating the Jacobian term from Equation (4): ϕ i t code = 1 c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) − ν i t code + β + ǫ ( D,L,A ) i (5) β ∈ [0 , t code ) . The new subscript ( D, L, A ) indicates that the error term now compensates also for the use of the a prioriestimates ¯ b and ¯ ℓ for s and ˚ ℓ .Van Diggelen uses these constraints to set the ν i ’s in a particular way. The method selects one index j that is used to set ν j and β and then resolves all the other ν i s so they are consistent with this β . Thatis, he assumes that ǫ ( D,L,A ) j = 0 so ν j = & c (cid:13)(cid:13) ¯ ℓ − ρ j (cid:0) t j − d j − ¯ b (cid:1)(cid:13)(cid:13) − ϕ j t code t code ' β = ( ν j + ϕ j ) t code − c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) . The method now substitutes this β in all the other constraints and assigns the other ν i s by setting ǫ ( D,L,A ) i = 0 and rounding, ν i = $ c (cid:13)(cid:13) ¯ ℓ − ρ j (cid:0) t j − d j − ¯ b (cid:1)(cid:13)(cid:13) − ϕ j t code + βt code ' . (6)When k ˚ ℓ − ¯ ℓ k and | s − ¯ b | are small enough, this gives a set of ν i s that are correct in the sense that theyall differ from the correct ν i s by the same integer.Van Diggelen chooses j in a particular way: he chooses the j that minimizes the magnitude of (3),which corresponds to the satellite closest to the zenith of ¯ ℓ at t j − ¯ b . In our notation, Van Diggelen’sjustification for this choice is as follows. He searches for a j for which Equation (5) approximates wellEquation (4). The difference between the two is c J j, : (cid:20) ˚ ℓ − ¯ ℓs − ¯ b (cid:21) , and his choice leads to a row of J with small-magnitude elements that is usually almost orthogonal to ˚ ℓ − ¯ ℓ . This leads to an estimated β that is relatively accurate, which helps resolve the correct ν i ’s.Van Diggelen also shows that if we resolve the ν i ’s by setting each ǫ ( D,L,A ) i = 0 separately , then theresolved β ’s might be close to in one equation and close to t code in another; this leads to inconsistent ν i ’s and to a huge position error. 5 .2 Final Resolution of the Receiver’s Location Van Diggelen’s method resolves the integer ν i ’s in Equation (5). Now we need to resolve the continuousunknowns. We do so using Gauss-Newton iterations on Equation (4), iterating on δ ℓ = ˚ ℓ − ¯ ℓ , δ s = s − ¯ b ,and β but keeping ν fixed. We start with δ ℓ , δ s , and β set to zero.In every iteration, we use the current iterates to produce estimates of the location and bias, ˆ ℓ = ¯ ℓ + δ ℓ ˆ b = ¯ b + δ s . We use them to improve the estimate of the ranges d i , setting ˆ d i = (cid:13)(cid:13)(cid:13) ˆ ℓ − ρ i (cid:0) ˆ t D,i (cid:1)(cid:13)(cid:13)(cid:13) . This allows us to reduce the errors in Equation (1), ϕ i t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:16) t i − ˆ d i − ˚ b (cid:17)(cid:13)(cid:13)(cid:13) − n i t code +˚ b + ǫ ( ˆ D ) i = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:16) t i − ˆ d i − ˚ b (cid:17)(cid:13)(cid:13)(cid:13) − ν i t code + β + ǫ ( ˆ D ) i (the second line holds because n i t code +˚ b = ν i t code + β , by definition). We again linearize this and solvethe constraints ϕ i t code = 1 c (cid:13)(cid:13)(cid:13) ˆ ℓ − ρ i (cid:16) t i − ˆ d i − ˆ b (cid:17)(cid:13)(cid:13)(cid:13) + 1 c ˆJ i, : (cid:20) ˚ ℓ − ˆ ℓs − ˆ b (cid:21) − ν i t code + β + ǫ ( ˆ D,L ) i . (7)for ˚ ℓ , s , and β using in the generalized least-squares sense, where the Jacobian is evaluated at ˆ ℓ and t − ˆ d − ˆ b .We can now explain why Van Diggelen’s method resolves the integers only once and iterates only onthe continuous unknowns. The ν i ’s that Van Diggelen’s method resolves are not equal to the n ′ i s in thenonlinear Equation 1. But when the linearization error is small enough, the two integer vectors differ bya constant, ⌊ ˚ b/t code ⌋ . This difference is compensated for by the integer part of the continuous variable β , which is not constrained to [0 , t code ) in the Gauss-Newton iterations. This is the actual function ofshadowing; to allow β to compensate not only for the clock error, but also for the constant error in ν .When the initial linearization error is so large that n − ν is no longer a constant, the method breaksdown. A different approach, which has never been proposed for snapshot positioning, is to add regularizationconstraints that will allow us to resolve all the m unknowns in the m instances of Equation (2) ϕ i t code = 1 c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − (˚ n i − ¯ n i ) t code − ¯ n i t code + (cid:16) ˚ b − ¯ b (cid:17) + ¯ b + ǫ ( D,L ) i using mixed-integer least-squares techniques. Note that we have rewritten Equation (2) in a way empha-sizes a change of variables that facilitate iterative improvements: the new unknowns are ˚ ℓ − ¯ ℓ , ˚ b − ¯ b , and ˚ n − ¯ n . We initially set ¯ b and ¯ n to zero.This section proposes two sets of regularizing equations and explains how to use this approach in aniterative Gauss-Newton solver. The first set of regularizing equations that we propose are c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) = 0 . ˚ ℓ and ˚ b are correct in the weighted least-squares sense. This leads to the following weighted mixed-integer least-squares problem: ˆ ℓ, ˆ b, ˆ n = arg min ℓ,b,n (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) W (cid:20) c J + (cid:2) m × m × (cid:3) − t code I m × m c J m × m (cid:21) ˚ ℓ − ¯ ℓ ˚ b − ¯ b ˚ n − ¯ n − (cid:20) t code ϕ − c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + ¯ n i t code − ¯ b (cid:21)(cid:19)(cid:13)(cid:13)(cid:13)(cid:13) , where W is a weight matrix derived from the covariance matrix C of the error terms ǫ , W T W = C − .Now we have m constraints, which for m ≥ should allow us to resolve the integer n i s.We propose to choose a diagonal W as follows. We set the first m diagonal elements of W to thestandard deviation of the arrival-time estimator, say W i,i = σ ( t i ) ≈ ns . To set the rest, we use boundson the a priori estimates ¯ ℓ and ¯ b , denoted | x − ¯ x | ≤ x max | y − ¯ y | ≤ y max | z − ¯ z | ≤ z max (cid:12)(cid:12)(cid:12) ˚ b − ¯ b (cid:12)(cid:12)(cid:12) ≤ b max . By the triangle inequality (cid:12)(cid:12)(cid:12)(cid:12) J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ≤ | J i, | x max + | J i, | y max + | J i, | z max + | J i, | b max . We define r i = | J i, | x max + | J i, | y max + | J i, | z max + | J i, | b max , so (cid:12)(cid:12)(cid:12)(cid:12) J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ≤ r i . We set W m + i,m + i = 3 r i /c , say.This mixed-integer least squares minimization problem can be solved by a generic solver, such as oneof the solvers that have been developed for real-time kinematic (RTK), a method for solving combinedcode-phase and carrier-phase GNSS constraints. It turns out that Doppler shifts allow us to regularize Equation (2) in a more effective way. GNSSreceivers estimate not only the time of arrival of the signal, but also its Doppler shift. The estimatedDoppler shift is biased, because the of the inaccuracy of the receiver’s local (or master) oscillator; it isalso inexact. We now show a a novel technique to use the Doppler-shift observations to to regularizeEquation (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 about 800 m/s. This assumptioncan be easily removed, but its removal leads to an additional unknown and more complicated expressionsthat we do not present here. The other assumption is that the local oscillator and the sampling clockin the receiver are derived from a single master oscillator in a certain (very common) way. Again, thisassumption can be removed if another unknown is added.The Doppler-shift formula for velocities much lower than the speed of light is ˚ D i ≈ − c ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:13)(cid:13)(cid:13) f . The Doppler observations that the receiver makes are D i = − c ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:13)(cid:13)(cid:13) f + ˚ f + ǫ ( δ ) i (8)7here ˚ f is the frequency offset (bias) of the receiver and ǫ ( δ ) i is an error term that represents the observa-tion error and the (negligible) slow-speed approximation. Therefore, the quantities − cD i /f are biasedestimates of the range-rate. We denote the a priory estimates of the Doppler shifts by ¯ D i .We differentiate Equation (2) by time, ddt ( ϕ i t code ) = ddt (cid:18) c (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) + 1 c J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − n i t code + (cid:16) ˚ b − ¯ b (cid:17) + ¯ b (cid:19) + ǫ ( D,L,∂ ) i . We first manipulate the equation a bit, to make it easier to differentiate: ddt (cid:0) cϕ i t code − (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) − c ¯ b (cid:1) = ddt (cid:18) J i, : (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) − cn i t code + c (cid:16) ˚ b − ¯ b (cid:17)(cid:19) + ǫ ( D,L,∂ ) i . (9)We denote H = J , J , J , J , + c − ct code ... ... ... ... . . . J i, J i, J i, J i, + c − ct code ... ... ... ... . . . J m, J m, J m, J m, + c − ct code . The first three columns of H are identical to those of J , the next is the fourth column of J but shiftedby c , and the last m columns consist of a scaled identity matrix. We now express the derivative on theright-hand side of Equation (9) as ddt H ˚ ℓ − ¯ ℓ ˚ b − ¯ bn . . .n m = H ddt ˚ ℓ − ¯ ℓ ˚ b − ¯ bn . . .n m + (cid:18) ddt H (cid:19) ˚ ℓ − ¯ ℓ ˚ b − ¯ bn . . .n m . We assume that the receiver is stationary, so ˚ ℓ − ¯ ℓ is time-independent, so ddt (˚ ℓ − ¯ ℓ ) = 0 . The derivativesof the integers n , . . . , n m are also zero. The derivative of the remaining element in the vector, ddt (˚ b − ¯ b ) , isnot zero and will need to be estimated. It represents the frequency offset of the receiver, which biases theobserved Doppler shift. It is multiplied by a column whose elements are very close to c (the range-rateis tiny relative to the speed of light), allowing it to compensate for the frequency bias.To differentiate H , we exploit the known structure of J . For each satellite, J i, is the negation ofthe so-called line-of-sight vector e Ti , which is the normalized direction from the satellite to the receiver;element J i, is the range-rate. The derivatives of these quantities are shown by Van Diggelen (Diggelen2009, Equation 8.6), Fernández-Hernández and Borre (Fernández-Hernández and Borre 2016), and othersources: ddt J i, = e T ddt k ¯ ℓ − ρ i k + (cid:0) ddt (cid:0) ¯ ℓ − ρ i (cid:1)(cid:1) T k ¯ ℓ − ρ i k . where the satellite position ρ i and its velocity ddt ρ i are taken at ˆ t D . To reduce the number of unknowns,we assume that the receiver is stationary, so ddt ¯ ℓ = 0 , so ddt J i, = e T ddt k ¯ ℓ − ρ i k − (cid:0) ddt ρ i (cid:1) T k ¯ ℓ − ρ i k . Element J i, is the range-rate of satellite i , so its derivative with respect to time is the satellite’s rangeacceleration, ddt J i, = d dt k ¯ ℓ − ρ i k . We use finite differences to evaluate this second derivative. The fourth column of H is J : , + c , but thederivative of c is obviously zero. The derivative of − ct code is also zero, so (cid:18) ddt H (cid:19) ˚ ℓ − ¯ ℓ ˚ b − ¯ bn . . .n m = (cid:18) ddt H : , (cid:19) (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) .
8e now derive the left-hand side of Equation (9), ddt (cid:0) cϕ i t code − (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) − c ¯ b (cid:1) . The derivative of the a prioi range estimate (cid:13)(cid:13) ¯ ℓ − ρ i (cid:0) t i − d i − ¯ b (cid:1)(cid:13)(cid:13) is the a-priori range-rate, which wecan compute. The derivative of c ¯ b is zero.To understand the first term, recall that ( n i + ϕ i ) t code = 1 c (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) +˚ b + ǫ i . so c ddt ϕ i t code = ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i ( t D,i ) (cid:13)(cid:13)(cid:13) + c ddt ˚ b + c ddt ǫ i . We now rewrite (8) as c ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:13)(cid:13)(cid:13) f = − D i + ˚ f + ǫ ( δ ) i or ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ i (cid:13)(cid:13)(cid:13) = − cD i f + c ˚ ff + cf ǫ ( δ ) i . We now substitute in the left-hand side of Equation (9): ddt cϕ i t code = − cD i f + c ˚ ff + cf ǫ ( δ ) i + c ddt ˚ b + c ddt ǫ i . The term ˚ f /f is the relative local-oscillator error in the receiver. If the oscillator runs too fast, ˚ f isnegative. Assuming that all the clocks in the receiver are derived from a master oscillator, if it runs toofast, ˚ b grows over time. Under this assumption − c ˚ ff = c ddt ˚ b so these terms cancel each other. If our assumption on the receiver does not hold, we would need toestimate c ˚ ff + c ddt ˚ b . That’s it. We have arrived at a system of m linear equations that we use to regularize the mixed-integerequations. The equations are: − cf D − ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ (cid:13)(cid:13)(cid:13) = H : , (cid:18) ddt (cid:16) ˚ b − ¯ b (cid:17)(cid:19) + (cid:18) ddt H : , (cid:19) (cid:20) ˚ ℓ − ¯ ℓ ˚ b − ¯ b (cid:21) . (10)In this equation, D represents the vector of observed Doppler shifts, ddt (cid:13)(cid:13)(cid:13) ˚ ℓ − ρ (cid:13)(cid:13)(cid:13) is the vector of the apriori range rates, and ddt (cid:16) ˚ b − ¯ b (cid:17) is a new scalar unknown. We have explained above how to compute H : , and ddt H : , . Solving the linearized and regularized mixed-integer least-squares system improves the initial a prioriestimates of ˚ ℓ and ˚ b , but not to the extent possible given the code phases. The most important factorthat limits the accuracy of the corrections is the fact that when the a priori estimates are large, theresolved integers, the n i ’s, are inexact. Therefore, we incorporate the mixed-integer solver into a Gauss-Newton-like iteration in which we correct all the unknowns, including the integer ambiguities, more thanonce.More specifically, once we solve the mixed-integer system for ˚ n − ¯ n , ˚ ℓ − ¯ ℓ and ˚ b − ¯ b (and for ˚ f − ¯ f inthe Doppler formulation), we use the corrections to improve the estimates of the receiver’s location andof the departure times and we linearize Equation (1) again. We now solve the newly-linearized systemagain for additional corrections, and so on.It is often not necessary to estimate all the unknowns in every iteration. We explore this behaviorexperimentally below. 9 Implementation and Evaluation
We have implemented all the methods that we described above in MATLAB.We use Borre’s Easy Suite (Borre 2003; Borre 2009) to perform many routine calculations. In par-ticular, we use it to correct GPS time (check_t), to correct for Earth rotation during signal propagationtime (e_r_corr), to read an ephemeris from a RINEX file and to extract the data for a particular satel-lite (rinexe, get_eph and find_eph), to transform Julian dates to GPS time (gps_time), to representJulian dates as one number (julday), to compute the coordinate of a satellite at a given time in ECEFcoordinates (satpos), to compute the azimuth, elevation, and distance to a satellite (topocent, which callstogeod 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 FileExchange) to approximate the ionopheric delay using the Klobuchar model. We take the parameters forthe Klobuchar model from files published by the GNSS Research Center at Curtin University .During the Gauss–Newton phase of the algorithm (after the integers have been determined), if we haveonly 4 observations, we add a pseudo-measurement constraint that constrains the correction to maintainthe height of the target, in the least-squares sense (Diggelen 2009).We use Chang and Zhou’s MILES package (Chang and Zhou 2007) to solve mixed integer least-squaresproblems.We take ephemeris data from RINEX navigation files published by NASA .We filter satellites that are lower than 10 degrees above the horizon which have the lowest SNR andare 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. Weused these to test our algorithms in the initial phases of the research. These results are not shownhere.• GPS Simulations. We generated satellite positions ephemeris files and used them to compute timesof arrival and code phases. These simulations do not include ionospheric or tropospheric delays, sothey 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 re-ceiver, connected to an ANN-MB-00 u-blox antenna mounted on a steel plate on top of a roofwith 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) .The WGS84 coordinates of the antenna are 32.1121756, 34.8055775 with height above sea level of61.15 m. The code phase measurements are included in the UBX-RXM-MEASX emitted by thereceiver. The data set includes about 700 epochs, one every minute (so they span a little more than11 hours). 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 calledVesper. It was designed and produced by Alex Schwartz Developments on the basis of an earliertag called Robin designed and produced by a company called CellGuide that no longer exists. Thetag records 1-bit RF samples at a rate of 1023000 samples per second. The tag was configured torecord a 256 ms sample every 10 minutes for a few hours. It was placed next to the ANN-MB-00antenna.Figure 1 compares the probability of success achieved by our regularized mixed-integer least-squares(MILS) solver with that achieved by Van Diggelen’s method. We considered fixes that are within 1 kmof the true location to be a success in obtaining the correct integer values. Each pixel in these heatmapsrepresents 16 different runs. Each run uses a random epoch, a random initial location estimate, and arandom 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. Ineach pixel, half of the initial time errors are positive and half are negative.The results clearly show that the MILS algorithm, even with the simple a-priori time and locationregularization from Section 4.1, outperforms Van Diggelen’s. Van Diggelen’s original algorithm obtains http://saegnss2.curtin.edu/ldc/rinex/daily/ https://cddis.nasa.gov/archive/gnss/data/daily/ https://axis-gps.com an Diggelen's (Original) initial position error (km) i n i t i a l t i m e e rr o r ( s ) MILS with A Priori Regularization initial position error (km) i n i t i a l t i m e e rr o r ( s ) MILS with Doppler Regularization initial position error (x1000 km) i n i t i a l t i m e e rr o r ( H r s ) Figure 1: The probability of obtaining an integer-correct fix on the ublox data set using three differentalgorithms. 11 Final Position Absolute Error (m) CD F Initial Time Error = 20.0 (s)Initial Position Error = 20.0 (km)Num Satellites = All Satellites In View ( m ) Final Position Absolute Error (m) CD F Initial Time Error = 20.0 (s)Initial Position Error = 20.0 (km)Num Satellites = 5 ( m ) Final Position Absolute Error (m) CD F Initial Time Error = 20.0 (s)Initial Position Error = 20.0 (km)Num Satellites = 4 ( m ) Final Position Absolute Error (m) CD F Initial Time Error = 150.0 (s)Initial Position Error = 150.0 (km)Num Satellites = All Satellites In View ( m ) Final Position Absolute Error (m) CD F Initial Time Error = 150.0 (s)Initial Position Error = 150.0 (km)Num Satellites = 5 ( m ) Final Position Absolute Error (m) CD F Initial Time Error = 150.0 (s)Initial Position Error = 150.0 (km)Num Satellites = 4 ( m ) Van Diggelen'sMILS (A Priori Reg)MILS (Doppler Reg)
Absolute Error CDF - Performance Comparison
Figure 2: Cumulative distribution of the final position error on the ublox data set, when the initialerrors are 20 km and 20 s or 150 km and 150 s. The dashed vertical lines show the 90% point for theDoppler regularized solver. Runs that only used observations from 4 satellites constrained the correctionto be perpendicular to the current iterate in ECEF coordinates (thereby maintaining the height for smallcorrections).a correct fix in almost all cases (success probability close to 1) when the initial location error is smalland the initial time error is 150 s or less, when the initial location time error is small and the initiallocation error is 100 km or less, and in other equivalent combination of time and location errors. Thecorresponding limits for the ILS algorithm are about 150 s and 200 km.With Doppler regularization, the advantage over Van Diggelen’s algorithm is dramatic. The solverobtains a correct fix as long as the initial time error is at most 1800 s; this works even with great-circledistances of 20,000 km, which means that the initial position can be essentially anywhere on Earth. Ifthe initial position error is small, the method can tolerate initial time errors of up to two hours (7200 s).Figure 2 shows the distribution of the final position error on the ublox data set. When using all thehigh-elevation satellites in view, the error was at most 20 m in about 90% of the fixes. This result mainlyquantifies the quality of the corrections that we used (e.g., for ionospheric and tropospheric delays); inprinciple the errors are also influenced by the SNR, but the setup guaranteed good SNR, so this is not asignificant factor. The figure also shows that our new MILS method can indeed produce fixes with onlyfour satellites in view. The graphs show the degradation in the error when we use a random subset offour satellites in every epoch (the random subsets were truly random and were allowed to include low-elevation satellites). The errors are naturally larger than when all the satellites in view are used, bothdue to the fewer constraints and due to the use of low-elevation satellites. As expected from Figure 1,Doppler regularization has a dramatic effect when the initial errors are large and essentially no effectwhen they are small.
We have shown that Van Diggelen’s ingenious coarse-time navigation algorithm that estimates a locationfrom GNSS observations without departure times is essentially a specialized solver for a mixed-integerleast-squares problem. Even though Van Diggelen’s algorithms has been cited and used by many authors,12he actual form of the mixed-integer optimization problem has never been presented; we present it in thispaper for the first time.We also show that the integer ambiguities can be resolved by regularizing the mixed-integer least-squares problem. We proposed two regularization techniques, one that biases solutions towards an initiala priori estimate. This extends Van Diggelen’s use of the a-priori estimate to resolve the integers, but ourregularization approach can resolve the integers with larger initial errors than Van Diggelen’s. In effect,the general mixed-integer formulation uses the available information more effectively than Van Diggelen’sspecialized solver.We also proposed a regularization method based on Doppler-shift observations. This method allowsour solver to resolve the correct integers as long as the initial time error is at most 1800 s or so (halfan hour), no matter how large the initial position error is. Doppler shifts have been used in snapshotpositioning before, but they were always used to produce an initial position and time estimate that issubsequently used as an a priori estimate in Van Diggelen’s algorithm.To achieve these results, our algorithm iterates over the entire mixed-integer least-squares problemmore than once. If one resolves and fixes the integers in the first iteration and continues to iterate onlyon the continuous unknowns, the method converges, but to fixes with larger errors than if iterating onthe integers as well.In effect, by cleanly formulating the mixed-integer optimization problem that underlies snapshotpositioning, we enabled the exploration of a wide range of solvers, including the two regularized solversthat we presented here. We believe that additional solvers can be discovered for this formulation. Incontrast, all prior research treated Van Diggelen’s algorithm as a clever black box, limiting the range ofalgorithms that can be developed.Our new methods are inspired by the real-time kinematic (RTK) method, which resolves a positionfrom both code-phase and carrier-phase GNSS observations. In RTK, the position is eventually resolvedby carrier-phase constraints, which have integer ambiguities; these constraints are regularized by pseudo-range constraints, which are less precise but have no integer ambiguities. Here the position is resolvedby integer-ambiguous code-phase constraints, which are regularized by either a priori estimates or byDoppler-shift observations. RTK also requires so-called differential constraints at a fixed receiver, becausethe carrier phase of the satellites are not locked to each other. Here we do not require differentialcorrections because the code departure times are locked to a whole milliseconds in all the satellites.
Acknowledgments
This study was also supported by grant 1919/19 from the Israel Science Founda-tion. Thanks to Aya Goldshtein for assisting with collecting data from the bat-tracking GPS snapshotlogger. Thanks to Amir Beck for helpful discussions.
References
Badia-Solé, Oriol and Tudor Iacobescu Ioan (2010).
GPS snapshot techniques . Tech. rep. Danish GPSCenter, Aalborg University. url : http://kom.aau.dk/group/10gr815/Report.pdf .Bavaro, Michele (Dec. 2012a). Fruit-Bat Tracking with Robin GPS Logger . Michele’s GNSS blog. http://michelebavaro.blogspot.com ,accessed 18 May 2020. url : http://michelebavaro.blogspot.com/2012/08/the-gps-tracker-to-beat-smallest.html .— (Aug. 2012b). The GPS tracker to beat (smallest, lightest, and most durable) . Michele’s GNSS blog. http://michelebavaro.blogspot.com , accessed 18 May 2020. url : michelebavaro.blogspot.com/2012/08/the-gps-tracker-to-beat-smallest.html .Bissig, P., M. Eichelberger, and R. Wattenhofer (2017). “Fast and Robust GPS Fix Using One Millisecondof Data”. In: , pp. 223–234.Borre, Kai (2003). “The GPS Easy Suite: Matlab code for the GPS newcomer”. In: GPS Solutions url : https://doi.org/10.1007/s10291-003-0049-3 .— (Mar. 2009). “GPS Easy Suite II: a Matlab Companion”. In: Inside GNSS , pp. 48–52. url : .Borre, Kai and Gibert Strang (2012). Algorithms for Global Positioning . Wellesley-Cambridge.Chang, Xiao-Wen and Tianyang Zhou (2007). “MILES: MATLAB package for solving Mixed Integer LEastSquares problems”. In:
GPS Solutions
11, pp. 289–294. url : https://doi.org/10.1007/s10291-007-0063-y .Cvikel, Noam et al. (2015). “Bats aggregate to improve prey search but might be impaired when theirdensity becomes too high”. In: Current Biology
25, pp. 206–211.Diggelen, Frank van (2009).
A-GPS: assisted GPS, GNSS, and SBAS . Artech House.Dujon, Antoine M., R. Todd Lindstrom, and Graeme C. Hays (2014). “The accuracy of Fastloc-GPSlocations and implications for animal tracking”. In:
Methods in Ecology and Evolution doi : . 13ichelberger, Manuel, Ferdinand von Hagen, and Roger Wattenhofer (2019). “Multi-Year GPS TrackingUsing a Coin Cell”. In: Proceedings of the 20th International Workshop on Mobile Computing Sys-tems and Applications (HotMobile) . Santa Cruz, CA, USA, pp. 141–146. isbn : 9781450362733. doi : . url : https://doi.org/10.1145/3301293.3302367 .Fernández-Hernández, I. and K. Borre (2016). “Snapshot positioning without initial information”. In: GPS Solutions
20, pp. 605–616. url : https://doi.org/10.1007/s10291-016-0530-4 .Harten, Lee et al. (2020). “The ontogeny of a mammalian cognitive map in the real world”. In: Science doi : . url : https://science.sciencemag.org/content/369/6500/194 .Lannelongue, Stephane and Pedro Pablos (1998). “Fast Acquisition Techniques For GPS Receivers”. In: Proceedings of the 54th Annual Meeting of The Institute of Navigation . Denver, CO, pp. 261–269.Liu, Jie et al. (2012). “Energy Efficient GPS Sensing with Cloud Offloading”. In:
Proceedings of the 10thACM Conference on Embedded Networked Sensor Systems (SenSys) , pp. 85–98.Muthuraman, Kannan, Jim Brown, and Mangesh Chansarkar (Jan. 2012). “Coarse Time Navigation:Equivalence of Algorithms and Reliability of Time Estimates”. In:
Proceedings of International Techni-cal Meeting of The Institute of Navigation , pp. 1115–1138. url : .Othieno, Nicholas and Scott Gleason (2012). “Combined Doppler Time-free Positioning for Low DynamicsReceivers”. In: Proceedings of the IEEE/ION Position, Location and Navigation Symposium , pp. 6–65.Peterson, Benjamin, Richard Hartnett, and Geffrey Ottman (1995). “GPS receiver structures for theurban canyon”. In:
Proceedings of the 8th International Technical Meeting of the Satellite Division(ION GPS) . Palm Springs, CA: The Institute of Navigation, pp. 1323–1332.Ramos, Heitor S. et al. (2011). “LEAP: A Low Energy Assisted GPS for Trajectory-Based Services”. In:
Proceedings of the 13th International Conference on Ubiquitous Computing (UbiComp , pp. 335–344. isbn : 9781450306300. doi : . url : https://doi.org/10.1145/2030112.2030158 .Tomkiewicz, Stanley M. et al. (2010). “Global positioning system and associated technologies in ani-mal behaviour and ecological research”. In: Philosophical Transactions of the Royal Society B doi : .Van Diggelen, Frank (July 2002). Method and Apparatus for Time Free Processing of GPS Signals . USPatent 6,417,801.W., Chen H. et al. (2014). “A new coarse-time GPS positioning algorithm using combined Doppler andcode-phase measurements”. In:
GPS Solutions
18, pp. 541–551.Wang, M., H. Qin, and T. Jin (2019). “Massive terminal positioning system with snapshot positioningtechnique”. In:
GPS Solutions url : https://doi.org/10.1007/s10291-018-0821-z .Weiss, Anthony J. (2004). “Direct Position Determination of Narrowband Radio Frequency Transmitters”.In: IEEE Signal Processing Letters
11, pp. 513–516.Witt, M. J. et al. (2010). “Assessing accuracy and utility of satellite-tracking data using Argos-linkedFastloc-GPS”. In:
Animal Behaviour issn : 0003-3472. doi : https://doi.org/10.1016/j.anbehav.2010.05.022 . url : .Yoo, Won Jae et al. (2020). “A coarse-time positioning method for improved availability”. In: GPS Solu-tions url : https://doi.org/10.1007/s10291-019-0919-yhttps://doi.org/10.1007/s10291-019-0919-y