Method, apparatus and system for estimating user position with a satellite positioning system in poor signal conditions

[0002] The present invention relates to time-of-arrival ranging receivers, and more specifically to global positioning system (GPS) receivers.

[0003] A global positioning system (GPS) receiver measures the time of arrival of a satellite ranging code provided by a GPS satellite. Based on the measured time of arrival of the ranging code (which includes orbit parameters and a timestamp indicating when the spacecraft transmitted the ranging code), the receiver determines the time required for the signal to propagate from the satellite to the receiver, and then multiplies this time by the speed of light to obtain a so-called pseudorange; the term pseudorange is used because there are a number of errors involved in the above measurement, the primary error being the error in the time of arrival of the signal because of an offset in the receiver clock from true GPS system time. The true range is just the difference between the time of reception of a specific code phase (of the ranging code) by the GPS receiver and the time of transmission of the specific code phase by the satellite multiplied by the speed of light, i.e.

[0004] r=(T_u−T_s)c,  (1)

[0005] where r is the (true) range, T_uis the GPS system time at which a specific code phase is received by the GPS receiver (i.e. the user of the GPS), T_sis the GPS system time at which the specific code phase was transmitted, and c is the speed of light. Because of an offset t_uof the receiver clock from GPS system time, what is measured is not the true range r but a pseudorange ρ, related to the true range according to,

[0006] ρ=c[(T_u+t_u)−T_s]=r+ct_u.  (2)

[0007] The true range r depends on the sought-after user position {right arrow over (u)} (at the time of receipt of the specific code phase) as well as on the position {right arrow over (s)} of the satellite transmitting the specific code phase (i.e. the position of the satellite at the time of transmission of the specific code phase of the ranging code), according to simply,

[0008] r=∥{right arrow over (s)}(t_tx)−{right arrow over (u)}(t_rx)∥,  (3)

[0009] where the time of transmission t_txand the time of reception t_rxare both according to the same clock, the GPS system clock. (FIG. 1 illustrates the notation of equation (3).) Therefore, from equations (2) and (3), the (measured) pseudorange depends on the sought-after user position and receiver clock offset according to,

[0010] ρ=∥{right arrow over (s)}(t_tx)−{right arrow over (u)}(t_rx)∥+ct_u,  (4)

[0011] in which the satellite position {right arrow over (s)}(t_tx) is known from the navigation data modulating the ranging code, since the navigation data provide the orbit parameters and a timestamp indicating the time at which the spacecraft transmitted the ranging code, and so indicating the location of the spacecraft along its orbit at the time of transmitting the (specific code phase of the) ranging code.

[0012] The problem of determining the offset of the receiver clock is called the time-recovery problem. The prior art teaches solving for the receiver clock offset, as part of the solution for the user position, using various techniques, including linearizing a system of equations for pseudoranges from four or more satellites (using least squares estimation in case of more than four satellites) as well as other techniques such as using a Kalman filter. In a typical solution, at least four satellites are needed because there are four unknowns being solved for all three coordinates of the user position and the receiver clock offset. The approach based on linearizing a system of four pseudorange equations begins with the pseudorange equation, i.e. equation (4), for each of four satellites, $\begin{array}{cc}\begin{array}{c}{\rho }_i=\Vert {\stackrel{\rightharpoonup }{s}}_i\left(t_{\mathrm{tx}}\right)-\stackrel{\rightharpoonup }{u}\left(t_{\mathrm{rx}}\right)\Vert +{\mathrm{ct}}_u\\ =\sqrt{{\left(x_i-x_u\right)}^2+{\left(y_i-y_u\right)}^2+{\left(z_i-z_u\right)}^2}+{\mathrm{ct}}_u\\ =f_i\left(x_u,y_u,z_u,t_u\right),\mathrm{for}i=1,2\dots ,4.\end{array}& \left(5\right)\end{array}$

[0013] The method the Taylor-series expands each of the equations (5) about an estimated position {right arrow over ({circumflex over (u)})} (related to the sought-after user position {right arrow over (u)} according to {right arrow over (u)}={right arrow over ({circumflex over (u)})}+Δ{right arrow over (u)}) and an estimated offset {circumflex over (t)}_urelated to the sought-after receiver offset t_uaccording to t_u={circumflex over (t)}_u+Δt_u), and neglects second-order and higher terms in the expansion, leaving, $\begin{array}{cc}{\rho }_i=f_i\left(\widehat{\stackrel{\rightharpoonup }{u}},{\hat{t}}_u\right)+{▽}_u{f}_i\left(\widehat{\stackrel{\rightharpoonup }{u}},{\hat{t}}_u\right)·\Delta \stackrel{\rightharpoonup }{u}+\frac{\partial f_i\left(\widehat{\stackrel{\rightharpoonup }{u}},{\hat{t}}_u\right)}{\partial t_u}\Delta t_u,\mathrm{for}i=1,2\dots ,4& \left(6\right)\end{array}$

[0014] where the notations ∇_u𝑓_i({right arrow over ({circumflex over (u)})},{circumflex over (t)}_u) and ∂𝑓_i({right arrow over ({circumflex over (u)})},{circumflex over (t)}_u)/∂t_uindicate that the function 𝑓 is to be treated as a function of {right arrow over (u)} and t_ufor the indicated differentiation, and the result evaluated at {right arrow over ({circumflex over (u)})} and {circumflex over (t)}_u. Using equations (5), the equations (6) become, $\begin{array}{cc}{\rho }_i={\hat{\rho }}_i-\frac{{\stackrel{\rightharpoonup }{s}}_i-\widehat{\stackrel{\rightharpoonup }{u}}}{{\hat{r}}_i}\Delta \stackrel{\rightharpoonup }{u}+c\Delta t_u,\mathrm{for}i=1,2\dots ,4,& \left(7\right)\end{array}$

[0015] with {circumflex over (𝑓)}_ibeing used in place of 𝑓_i({right arrow over ({circumflex over (u)})},{circumflex over (t)}_u) in eq. (6). Using ${\stackrel{\rightharpoonup }{a}}_i=\frac{{\widehat{\stackrel{\rightharpoonup }{s}}}_i-\widehat{\stackrel{\rightharpoonup }{u}}}{{\hat{r}}_i}$

[0016] and Δ𝑓_i={circumflex over (𝑓)}_i−𝑓_i, equations (7) can be written as,

[0017] Δ𝑓_i={right arrow over (a)}_i·Δ{right arrow over (u)}+cΔt_u, for i=1, 2 . . . , 4,  (8)

[0018] and by constructing the vectors and matrix, $\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{\rho }\equiv \left(\begin{array}{c}\Delta {\rho }_1\\ \Delta {\rho }_2\\ \Delta {\rho }_3\\ \Delta {\rho }_4\end{array}\right)H=\left(\begin{array}{cccc}{a}_{\mathrm{x1}}& a_{\mathrm{y1}}& a_{\mathrm{z1}}& 1\\ a_{\mathrm{x2}}& a_{\mathrm{y2}}& a_{\mathrm{z2}}& 1\\ a_{\mathrm{x3}}& a_{\mathrm{y3}}& a_{\mathrm{z3}}& 1\\ a_{\mathrm{x4}}& a_{\mathrm{y4}}& a_{\mathrm{z4}}& 1\end{array}\right)\mathrm{and}\Delta x=\left(\begin{array}{c}\Delta x_u\\ \Delta y_u\\ \Delta z_u\\ -c{t}_u\end{array}\right),& \left(9\right)\end{array}$

[0019] equations (8) can be written as,

[0020] Δ{right arrow over (ρ)}=HΔ{right arrow over (x)}  (10)

[0021] which has the solution,

[0022] Δ{right arrow over (x)}=H⁻¹Δ{right arrow over (ρ)}.  (11)

[0023] When navigating by GPS in weak signal conditions such that the navigation component of the navigation signal cannot be decoded, the time when the signal is received is not precisely known. That time, when expressed according to the GPS clock, is called here simply GPS time at arrival and is indicated by the notation τ. It differs from the time at which the signal is transmitted by a satellite (i.e. t_txin equation (5)) by the time of flight t_𝑓 of the signal from the satellite to the GPS receiver. In terms of the GPS time at arrival τ and the time of flight t_𝑓, the pseudorange given by equation (2) can be written as,

[0024] ρ=c└(T_u+t_u)−(τ−t_𝑓)┘=c(T_u+t_u)−c(τ−t_𝑓)=c(t_u+t_𝑓),  (2′)

[0025] since, by definition, T_u=τ.

[0026] In normal circumstances (i.e. strong signal conditions), the time at which the ranging code was sent from the spacecraft, τ−t_𝑓 (i.e., in the notation of equation (2), simply T_s), is determined (to within the offset of the satellite clock δt, here taken to be negligible) by decoding the received navigation message. However, when the navigation data (containing the time of transmission/timestamp and the satellite orbit parameters) included in the ranging code sent by the satellite cannot be demodulated (decoded), but the baseband measurements (i.e. aligning the replica code with the transmitted code) can still be made, it is necessary to determine not only the receiver offset t_u, but also the time t_txat which the signal was sent by the satellite, which is the GPS time at arrival τ (the time according to the GPS clock at which the signal was received) minus the unknown time of flight, t_𝑓.

[0027] In the case of signal conditions that are so poor that the navigation component of the ranging code cannot be demodulated, the prior art teaches solving for the user position based on a five-step procedure. The procedure is based on the fact that the pseudorange ρ_ito the i^thspacecraft, which is a measured quantity in the case of strong signal conditions (as opposed to a calculated quantity like the user position), is the sum of two components, and accordingly can be written as,

[0028] ρ_i=ρ_i⁽¹⁾+ρ_i⁽²⁾  (12)

[0029] where ρ_i⁽¹⁾is obtained from the navigation data included with the ranging code, and is given by $\begin{array}{cc}{\rho }_i^{\left(1\right)}=c\lceil \frac{t_f}{1ms}\rceil & \left(13\right)\end{array}$

[0030] where ┌ . . . ┐ indicates the ceiling function (next greatest integer), and ρ_i⁽²⁾is obtained from baseband processing the ranging code (i.e. aligning the replica with the transmitted code), giving the submillisecond part of the pseudorange. The part ρ_i⁽¹⁾is here sometimes called the millisecond part, since it gives the pseudorange only to within a millisecond, i.e. to within the speed of light times a millisecond, or nearly 300 km. When the first component ρ_i⁽¹⁾cannot be determined (measured) because the navigation data cannot be demodulated, only the second component ρ_i⁽²⁾can be determined (measured). However, with cellular network assistance or some other means of assistance giving us the spacecraft orbit parameters (but not their precise position in the orbits) and the location of a nearby base station, (determined by methods known in the art) the prior art teaches solving for the user position in six steps as follows.

[0031] 1. Make an estimate of the GPS time at arrival. The estimate can be based on either the local clock or a clock at the reference position.

[0032] 2. Make estimates of the times of flight for each satellite, using e.g. satellite positions at the estimated GPS time at arrival. If no estimates are otherwise available, use an average of around 70 ms.

[0033] 3. From the estimates of the GPS time at arrival and times of flight, compute the times of transmission for each satellite and then the satellite positions at those times using orbital parameters for each of the satellites acquired via cellular network assistance.

[0034] 4. With the satellite positions approximately so determined, use the position of a nearby base station (called here the reference position) to calculate the millisecond part {tilde over (ρ)}_i⁽¹⁾of the pseudoranges for each satellite. The full pseudoranges to the receiver (not to base station) are computed by combining the millisecond part {tilde over (ρ)}_i⁽¹⁾and the submillisecond part {tilde over (ρ)}_i⁽²⁾, where the submillisecond part {tilde over (ρ)}_i⁽²⁾is taken (determined) from the receiver baseband.

[0035] 5. With the pseudoranges and satellite positions approximately so determined, use the iterative least squares method to determine an estimate of the receiver position. The iterative least squares method gives not only the position, but also an estimate of the a uncertainty in the solution.

[0036] 6. Repeat steps 1-5 using different estimates of the GPS time at arrival selected from a grid of values around the value provided for example by the cellular network assistance. Use as the user position the value for which the uncertainty is the smallest.

[0037] The third step above is quite apparent from the description above. If the time of transmission is known, then the satellite positions are calculated from this time and from orbital parameters using formulas well known in the field of GPS.

[0038] In the fifth step, using the spacecraft positions calculated in step three, the user position {right arrow over (u)}(t_rx) and the receiver clock offset t_uare solved using

[0039] {tilde over (ρ)}_i⁽¹⁾+ρ_i⁽²⁾=∥{right arrow over (s)}_i(τ−{tilde over (t)}_i^(𝑓))−{right arrow over (u)}(t_rx)∥+ct_u  (15)

[0040] At least four such equations are needed, so as to have a system of at least four equations in four unknowns (user position coordinates and receiver clock offset). If we have more than four equations, we can solve the system of equations using least squares, which will give us an uncertainty (error) as part of the solution, and we can iterate on the initial guesses (of first component of the pseudoranges, the times of flight, and the user position) to make the uncertainty smaller, in what is called iterated least squares.

[0041] As is evident from the above description, the procedure for solving for user position in case of weak signal conditions is time consuming. It would be advantageous to have a more efficient (faster) way of determining a GPS position-velocity-time (PVT) solution in weak signal conditions, conditions such that the navigation data cannot be demodulated (decoded) but the baseband measurement can still be performed.

[0042] Accordingly, the present invention provides a method and corresponding apparatus and system for solving for the position of a satellite positioning system receiver, the receiver having a receiver clock with a receiver clock offset and responsive to ranging codes from some number of satellites, the ranging code from each satellite including a navigation component conveying information about the satellite and having a time of flight from the satellite, the method of use in signal conditions where the ranging codes can be baseband processed so that at least a repeating part of the pseudorange for each satellite can be determined, the method including the steps of: including the system time at arrival as an independent and additional variable of the pseudorange function for each of the ranging codes, an independent variable in addition to the receiver position; providing estimates of the times of flight of the ranging codes and of the repeating part of the pseudoranges for each of the ranging codes; receiving information about the satellite sufficient to determine the orbits of the satellites, although not necessarily the positions of the satellites in their orbits; and solving simultaneously for the receiver position and the system time at arrival.

[0043] In a further aspect of the invention, the method also includes the receiver clock offset as an independent and additional variable. In some such applications, the method also includes the step of performing a Taylor-series expansion of the pseudorange function about an initial guess or estimate of each of the independent variables for as many of the satellites as there are independent variables thereby providing a system of pseudorange equations. The pseudorange equations are sometimes solved using a determinate single-point solution that is either an iterative least squares solution or an exact algebraic solution, depending on the number of pseudorange equations, and sometimes solved using a Kalman filter.

[0044] In another further aspect of the invention, the information about each satellite includes the orbit parameters for each satellite.

[0045] In yet another further aspect of the invention, the satellite positioning system is the global satellite positioning system (GPS) and the system time at arrival is the GPS time at arrival.

[0046] In yet still another further aspect of the invention, the repeating part of the pseudorange is the submillisecond part of the pseudorange.

[0047] In yet still even another further aspect of the invention, the method is of use in signal conditions where the ranging codes can be baseband processed providing at least the submillisecond part of the pseudorange.

[0048] Compared to the prior art procedure explained above, the present invention solves for the user position using the following steps, of which the first four are the same as in the prior art.

[0049] 1. Make an estimate of the GPS time at arrival. The estimate can be based on either the local clock or a clock at the reference position.

[0050] 2. Make estimates of the times of flight for each satellite, using e.g. satellite positions at the estimated GPS time at arrival. If no estimates are otherwise available, use an average of around 70 ms.

[0051] 3. From the estimates of the GPS time at arrival and the times of flight, compute the times of transmission for each satellite and then the satellite positions at those times using orbital parameters for each of the satellites acquired via cellular network assistance.

[0052] 4. With the satellite positions approximately so determined, use the position of a nearby base station (called here the reference position) to calculate the millisecond part {tilde over (ρ)}_i⁽¹⁾of the pseudoranges for each satellite. The full pseudoranges to the receiver (not to base station) are computed by combining the millisecond part {tilde over (ρ)}_i⁽¹⁾and the submillisecond part {tilde over (ρ)}_i⁽²⁾where the submillisecond part {tilde over (ρ)}_i⁽²⁾is taken (determined) from the receiver baseband.

[0053] 5. With the pseudoranges and satellite positions known, instead of using the iterative least squares method as in the prior art, include in the pseudorange equations the unknown GPS time at arrival τ as an additional, independent variable. Then use the iterative least squares technique for the modified equations, the iteration here being different than in the prior art. In the conventional solution of the prior art, the satellite positions are fixed while iterations refine the estimate of the user position. Here the satellite positions and the user position are refined at the same time. The satellite positions are refined by iterating on the GPS time at arrival τ, which amounts to adjusting the positions of the satellites in their orbits, with the orbits themselves known via cellular network assistance or some other source. (New estimates of the GPS time at arrival and previous estimates of the times of flight are used to compute new times of transmission. The new times of transmission are used to compute new satellite positions. Then the times of flight and pseudoranges are updated.) Each iteration gives a new estimate of the user position, the GPS time at arrival, and the receiver clock offset. The iterations here converge, i.e. the solutions stop changing, and so the iteration here is fundamentally different than that of the prior art.

[0054] In essence, the invention collapses the last steps of the solution of the prior art into one step by introducing the GPS time at arrival τ (the time at which the user receives the ranging code) as a fifth unknown in the pseudorange equation,

[0055] {tilde over (ρ)}_i=∥{right arrow over (s)}_i(τ−{tilde over (t)}_i^(𝑓))−{right arrow over (u)}(t_rx)∥+ct_u  (16)

[0056] where the pseudoranges {tilde over (ρ)}_iand the times of flight {tilde over (t)}_i^(𝑓)are (approximately) known by some means. Thus, the invention provides five unknowns (τ, u_x, u_y, u_z, t_u) which are solvable by at least five pseudorange equations (i.e. equation (16) for five satellites). In case of five such pseudorange equations, the invention provides either a determinate solution or a Kalman filter solution. In case of more than five pseudorange equations, the invention provides a solution using either the method of least squares or a Kalman filter.

[0057] The method provided by the invention is faster than methods known in the art (from five to ten times faster), and it produces clock times (receiver offsets) that are as accurate as at least some of the prior art methods (in particular, the so-called SSE method). In the method using a filter such as a Kalman filter, the time error can be determined more precisely and the extra computation load is minimized when the time error solution is integrated into the Kalman filter.

[0058] The above and other objects, features and advantages of the invention will become apparent from a consideration of the subsequent detailed description presented in connection with accompanying drawings, in which:

[0059] FIG. 1 is an illustration of the coordinate system used in the method of the invention;

[0060] FIG. 2 is a schematic of a GPS receiver executing a program for solving for the position of the GPS receiver (the user) according to the invention; and

[0061] FIG. 3 is a flowchart of the method of the invention.

[0062] The present invention will now be described in the particular context of the so-called global positioning system (GPS). There are, however, other satellite positioning systems in use today, and it is possible that there will be new satellite positioning systems developed in the future. It is to be understood that although the description here is made specific to the GPS context, the same description could be made in the specific context of any other analogous satellite positioning system by changing the terms of art used here to the corresponding terms in the other satellite positioning system. For example, the GPS term pseudorange used here would be changed to the corresponding term (which might also be pseudorange) used in the context of the other satellite positioning system. As another example, the term GPS time at arrival should be understood to be used in the GPS context, and in the more general context of any other analogous satellite positioning system, the term system time at arrival would be used. As yet another example, the term submillisecond part of the ranging code tends to be GPS specific, and that term when used here should be understood to correspond more generally to a repeating part (the less significant part) of the ranging code.

[0063] The present invention extends the methods of the prior art for determining the position of a global positioning system (GPS) receiver when the ranging code cannot be demodulated (decoded) because of poor signal conditions by Taylor-series expanding the pseudorange function 𝑓({right arrow over (u)},_u,τ), given by the equation

[0064] {tilde over (ρ)}_i=∥{right arrow over (s)}_i(τ−{tilde over (t)}_i^(𝑓))−{right arrow over (u)}(t_rx)∥+ct_u=𝑓({right arrow over (u)},t_u,τ),  (16)

[0065] not only about some user position {right arrow over (u)}(t_rx) and some receiver offset t_u, as in equations (6) of the prior art, but also about the GPS time at arrival τ, i.e. the (single) time, according to the GPS clock, at which the user (GPS receiver) receives the ranging codes from the different spacecraft. Although a precise determination of the receiver clock offset t_ushould include the same information as a precise determination of what is here called the GPS time at arrival τ (in the sense that if you know t_uprecisely, you should know τ), the GPS time at arrival as used in the invention likely does not convey the same information as the receiver clock offset. The GPS time at arrival is treated in the invention as independent of the receiver clock offset. By introducing the GPS time at arrival τ as one additional degree of freedom for the solution of the system of pseudorange equations (16), the invention approximately accounts for the uncertainties (errors) in the times of flight estimates {tilde over (t)}_i^(𝑓)for the different spacecraft at positions {right arrow over (s)}_i(τ−{tilde over (t)}_i^(𝑓).

[0066] Referring now to FIG. 2, in the preferred embodiment of the invention, a GPS receiver 10 includes a nonvolatile memory 10afor storing a program 10cthat, when running on a processor 10balso included in the GPS receiver, determines the position of the GPS receiver (user) by solving a system of pseudorange equations treating the GPS time at arrival τ as a fifth variable, in addition to the user position (three variables, one for each coordinate) and the GPS receiver clock offset from GPS system time (i.e. time according to the GPS clock). The program 10bsolves the system of pseudorange equations by linearizing the pseudorange equations as in the prior art, but instead of Taylor-series expanding only about some initial guess or estimate of the user position and about some initial guess or estimate of the receiver offset, as in equations (6) above, the invention also Taylor-series expands the pseudorange equations about the GPS time at arrival. It should be understood that although in the preferred embodiment of the invention there is a Taylor series expansion about some initial guess or estimate of the five variables of the pseudorange equation (16′), the invention is the use of the GPS time at arrival as at least one variable in solving for the position (or other aspect of the state of motion of) a user (GPS receiver), regardless of whether in the course of the solution a Taylor series expansion is performed.

[0067] In linearizing the pseudorange equation (16′) including the GPS time at arrival τ, the pseudorange equations for each of the various (at least five) satellites are written as, $\begin{array}{cc}\begin{array}{c}{\rho }_i=\Vert {\stackrel{\rightharpoonup }{s}}_i\left(t_{\mathrm{tx}}\right)-\stackrel{\rightharpoonup }{u}\left(t_{\mathrm{rx}}\right)\Vert +{\mathrm{ct}}_u\\ =\sqrt{{\left(x_i\left(\tau -t_f^{\left(i\right)}\right)-x_u\right)}^2+{\left(y_i\left(\tau -t_f^{\left(i\right)}\right)-y_u\right)}^2+{\left(z_i\left(\tau -t_f^{\left(i\right)}\right)-z_u\right)}^2}+{\mathrm{ct}}_u\\ =f_i\left(x_u,y_u,z_u,\tau ,t_u\right),\mathrm{for}i=1,2,\dots ,N.\end{array}& \left(17\right)\end{array}$

[0068] where {right arrow over (s)}_i(t_tx), with components (x_i,y_i,z_i), is the position of the i^thsatellite (out of N satellites, where N could be larger than four) at the time t_tx=τ−t_𝑓⁽ⁱ⁾at which the satellite transmits the specific code phase of the ranging code, τ is the GPS time at arrival, t_𝑓⁽ⁱ⁾is the time of flight of the ranging code, {right arrow over (u)}(t_rx) is the position of the GPS receiver (user) at the time t_rxof reception of the specific code phase of the ranging code, t_uis the offset from GPS system time of the receiver clock, and c is the speed of light.

[0069] Now instead of Taylor-series expanding the pseudorange function 𝑓 about only some estimated receiver position {right arrow over ({circumflex over (u)})} and some estimated receiver offset {circumflex over (t)}_u, i.e. about ({right arrow over ({circumflex over (u)})},{circumflex over (t)}_u), the function 𝑓 is Taylor-series expanded about the five-dimensional quasi-point ({right arrow over ({circumflex over (u)})},{circumflex over (τ)},{circumflex over (t)}_u), i.e. $\begin{array}{cc}{\rho }_i=f_i\left(\stackrel{\hat{\rightharpoonup }}{u},\hat{\tau },{\hat{t}}_u\right)+{\nabla }_u{f}_i\left(\stackrel{\hat{\rightharpoonup }}{u},\hat{\tau },{\hat{t}}_u\right)·\Delta \stackrel{\rightharpoonup }{u}+\frac{\partial f_i\left(\stackrel{\hat{\rightharpoonup }}{u},\hat{\tau },{\hat{t}}_u\right)}{\partial \tau }\Delta \tau +\frac{\partial f_i\left(\stackrel{\hat{\rightharpoonup }}{u},\hat{\tau },{\hat{t}}_u\right)}{\partial t_u}\Delta t_u& \left(18\right)\end{array}$

[0070] for i=1, . . . , N satellites, where the estimated user (GPS receiver) position {right arrow over ({circumflex over (u)})} is related to the sought-after user position {right arrow over (u)} according to {right arrow over (u)}={right arrow over ({circumflex over (u)})}+Δ{right arrow over (u)}, {circumflex over (τ)} is related to the sought-after GPS time at arrival τ according to τ={circumflex over (τ)}+Δτ, and the estimated receiver offset {circumflex over (t)}_uis related to the sought-after receiver offset t_uaccording to t_u32 {circumflex over (t)}_u+Δt_u, and, finally, where the notations ∇_u𝑓_i({right arrow over ({circumflex over (u)})},{circumflex over (τ)},{circumflex over (t)}_u), ∂𝑓_i({right arrow over ({circumflex over (u)})}, {circumflex over (τ)},{circumflex over (t)}_u)/∂τ and ∂𝑓_i({right arrow over ({circumflex over (u)})},{circumflex over (τ)},{circumflex over (t)}_u)/∂t_uindicate that the function 𝑓 is to be treated as a function of {right arrow over (u)}, τ and t_ufor the indicated differentiation, and the result evaluated at {right arrow over ({circumflex over (u)})}, {circumflex over (τ)} and {circumflex over (t)}_u. Using equations (17), the equations (18) become, $\begin{array}{cc}{\rho }_i={\hat{\rho }}_i-\frac{{\stackrel{\rightharpoonup }{s}}_i-\stackrel{\hat{\rightharpoonup }}{u}}{{\hat{r}}_i}·\Delta \stackrel{\rightharpoonup }{u}+\frac{{\stackrel{\rightharpoonup }{s}}_i-\stackrel{\hat{\rightharpoonup }}{u}}{{\hat{r}}_i}·{\stackrel{\rightharpoonup }{v}}_i+c\Delta t_u,\mathrm{for}i=1,2,\dots ,N,& \left(19\right)\end{array}$

[0071] where {right arrow over (s)}_iis the satellite position at the estimated GPS time at arrival {circumflex over (τ)}, i.e. {right arrow over (s)}_i({circumflex over (τ)}), {circumflex over (ρ)}_i=∥{right arrow over (s)}_i({circumflex over (τ)})−{right arrow over (u)}∥+c{circumflex over (t)}_u={circumflex over (r)}_i+c{circumflex over (t)}_u, and {right arrow over (ν)}_i=∂{right arrow over ({circumflex over (s)})}_i({circumflex over (τ)})/∂τ is the velocity of the i^thsatellite evaluated at the approximate GPS time at arrival {circumflex over (τ)}.

[0072] Using ${\stackrel{\rightharpoonup }{a}}_i=\frac{{\stackrel{\rightharpoonup }{s}}_i-\stackrel{\hat{\rightharpoonup }}{u}}{{\hat{r}}_i}$

[0073] (with {right arrow over (s)}_ievaluated at {circumflex over (τ)}) and Δρ_i={circumflex over (ρ)}_i−ρ_i, equations (19) can be written as,

[0074] Δρ_i={right arrow over (a)}_i·(Δ{right arrow over (u)}+{right arrow over (ν)}_iΔτ)+cΔt_u,  (20)

[0075] for i=1,2, . . . , N satellites, and by constructing the vectors and matrix (using for clarity of illustration only five satellites), $\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{\rho }\equiv \left(\begin{array}{c}\Delta {\rho }_1\\ \Delta {\rho }_2\\ \Delta {\rho }_3\\ \Delta {\rho }_4\\ \Delta {\rho }_5\end{array}\right)H=\left(\begin{array}{ccccc}{a}_{\mathrm{x1}}& a_{\mathrm{y1}}& a_{\mathrm{z1}}& {\stackrel{\rightharpoonup }{a}}_1·{\stackrel{\rightharpoonup }{v}}_1& 1\\ a_{\mathrm{x2}}& a_{\mathrm{y2}}& a_{\mathrm{z2}}& {\stackrel{\rightharpoonup }{a}}_2·{\stackrel{\rightharpoonup }{v}}_2& 1\\ a_{\mathrm{x3}}& a_{\mathrm{y3}}& a_{\mathrm{z3}}& {\stackrel{\rightharpoonup }{a}}_3·{\stackrel{\rightharpoonup }{v}}_3& 1\\ a_{\mathrm{x4}}& a_{\mathrm{y4}}& a_{\mathrm{z4}}& {\stackrel{\rightharpoonup }{a}}_4·{\stackrel{\rightharpoonup }{v}}_4& 1\\ a_{\mathrm{x5}}& a_{\mathrm{y5}}& a_{\mathrm{z5}}& {\stackrel{\rightharpoonup }{a}}_5·{\stackrel{\rightharpoonup }{v}}_5& 1\end{array}\right)\mathrm{and}\Delta \stackrel{\rightharpoonup }{x}=\left(\begin{array}{c}\Delta x_u\\ \Delta y_u\\ \Delta z_u\\ \Delta \tau \\ -{\mathrm{ct}}_u\end{array}\right),& \left(21\right)\end{array}$

[0076] equations (19) can be written as,

[0077] Δ{right arrow over (ρ)}=HΔ{right arrow over (x)}.  (22)

[0078] The use of the invention will be illustrated by application first to algorithmic solutions to the user position problem (i.e. either exact algebraic solutions or least squares solutions, both of which can be categorized as single-point solutions because of relying on only a single quasi-point as an initial guess to the user position problem), and then to Kalman filter solutions (which provide a smoothed navigation solution in which the effect of noise is reduced by taking into account two or more previous measurements as well as noise statistics).

[0079] FIG. 3 is a flowchart of the invention, showing that in the preferred embodiment the invention uses five independent variables, as explained above, namely the user position, the receiver clock offset, and the GPS time at arrival, and with outside assistance determines a system of pseudorange equations which according to the invention are then solved simultaneously by any of a number of methods, two of which are as follows.

[0080] Algorithmic (or Single-point) Solution

[0081] For five satellites, provided that H is non-singular, equation (22) has the solution,

[0082] Δ{right arrow over (x)}=H⁻¹Δ{right arrow over (ρ)}.  (23)

[0083] For more than five satellites, it has the least squares solution,

[0084] Δ{right arrow over (x)}=(H^TH)⁻¹H^TΔ{right arrow over (ρ)},  (24)

[0085] i.e. the value of Δ{right arrow over (x)} given by equation (24) minimizes the square of the residual vector {right arrow over (R)} defined as,

[0086] {right arrow over (R)}(Δ{right arrow over (x)})=HΔ{right arrow over (x)}−Δ{right arrow over (ρ)}.  (25)

[0087] The least squares solution at each iteration stage uses some guess for the receiver position {right arrow over (u)}, and then iterates to solve for the time of flight from each of the N satellites, using,

[0088] t_𝑓,i=∥{right arrow over (s)}_i(τ−t_𝑓,i)−{right arrow over (u)}(τ)∥/c  (26)

[0089] where i=1, . . . , N. If the solution results in a position/GPS time at arrival/offset estimation that is too far from the initial guess, i.e. if Δ{right arrow over (x)} is large, then another initial guess can be made and a new Δ{right arrow over (x)} can be calculated.

[0090] (Extended) Kalman Filter Solution

[0091] Background Information

[0092] A Kalman filter provides an estimate of a user's position based not only on current measurements, but also based on past estimates (more specifically, based on predictions of the current measurements). In essence, a Kalman filter provides a solution that filters out the effect of noise in current measurements. The discrete extended Kalman filter (EKF) is generally used with a stand-alone GPS receiver (as opposed to a network-assisted GPS receiver). Nevertheless, in cases where the tracking is often interrupted because of poor satellite signal quality, poor enough that the navigation component of the ranging code cannot be demodulated, cellular assistance in the form of orbit parameters and the location of a nearby cellular base station allows a Kalman filter to still provide the position of a user (GPS receiver), albeit an approximate position.

[0093] The Kalman filter formulation assumes a process model by which the state vector for the system including both the user and the GPS system (the space component, i.e. the satellites) is propagated from an earlier time to a later time. The state vector used here is taken to be,

[0094] x_L=[xyz|{dot over (x)}{dot over (y)}{dot over (z)}|t_u𝑓_clock|τ{dot over (τ)}],  (27)

[0095] where {right arrow over (u)}=(x,y,z) is the user (GPS receiver) position at the time of reception of the specific code phase, {right arrow over ({dot over (u)})}=({dot over (x)},{dot over (y)},{dot over (z)})=d{right arrow over (u)}/dt is the user velocity with respect to an internal time variable t (with the same rate of progression as the GPS time at arrival but with an unknown and irrelevant offset) and designated below by the discrete internal variable k (in e.g. equation (28)), t_uis the offset in the user (GPS receiver) clock, 𝑓_clockis the rate of progression of the user clock, τ is the GPS time at arrival according to the satellite clock (presumed to be correct), {dot over (τ)} is the drift in the satellite clock with respect to the internal time variable t. (Although the state vector of equation (27) is what is typically used in a Kalman filter solution, sometimes a state vector is defined to include fewer variables.)

[0096] A Kalman filter operates by processing all available measurements at each epoch, in what is called a measurement update step, and then propagating the state vector estimate and process noise covariance ahead to the next epoch, in what is called a time update step. Estimates of the state vector and process noise covariance after the measurement update are here indicated by a superscript “+”, while estimates and covariances propagated ahead to the next epoch are indicated here by a superscript “−”.

[0097] The process model described here (as one example) is the so-called low-dynamics model, i.e. the dynamical model usually used for GPS receivers in cars or boats, but not airplanes, for which the so-called high-dynamics model would be used. The process model, i.e. the state vector update equation, is here (in the case of the low-dynamics model) written as,

[0098] x_L(k)=Φ_L(Δt)x_L(k−1)+w_L(k−1)  (28)

[0099] where x_L(k) is the state vector at the discrete time instance k (the subscript L being used to indicate that the low-dynamics model is being used), w_L(k−1) is the process noise at instance k−1, Δt is the discretization interval (i.e. the time difference between instance k and instance k−1) , and Φ_L(Δt) is the one-step state transition matrix. (The Kalman filter described here is just an example. Once we obtain the modified pseudorange equation including the GPS time at arrival as one of the unknowns, we can use other dynamic models, besides the low-dynamics model, and other types of Kalman filters.)

[0100] The measurement model is here written as,

[0101] y_L(k)=G_L(k)x_L(k)+ε_L(k),  (29)

[0102] where y_L(k) is a vector of measured (observed) quantities, G_L(k) is the so-called connection matrix, connecting the state variables (components of x_L(k)) to the measured quantities (components of y_L(k)), and ε_L(k) is the measurement noise vector (the noise being assumed to be white, distributed normally as N(0,R_L(k)), i.e. with some predetermined measurement noise covariance matrix R_L(k)).

[0103] The state transistion matrix is here defined as $\begin{array}{cc}{\phi }_L\left(\Delta t\right)=\left[\begin{array}{cccc}I& \Delta \mathrm{tI}& 0& 0\\ 0& I& 0& 0\\ 0& 0& {\phi }_{\mathrm{clock}}\left(\Delta t\right)& 0\\ 0& 0& 0& {\phi }_{\mathrm{sat}}\left(\Delta t\right)\end{array}\right],& \left(30\right)\end{array}$

[0104] where each element is a 2×2 matrix, where I is the 2×2 identity matrix, where Δt is the discretization interval as above, where Φ_clock(Δt) is the transition matrix of a receiver clock process model, such as one of the widely used two-state models, $\begin{array}{cc}{\phi }_{\mathrm{clock}}\left(\Delta t\right)=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right],& \left(31\right)\end{array}$

[0105] and where Φ_sat(Δt) is the process model of satellite time, defined as, $\begin{array}{cc}{\phi }_{\mathrm{sat}}\left(\Delta t\right)=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right].& \left(32\right)\end{array}$

[0106] The initial condition used is

[0107] {dot over (τ)}=1.  (33)

[0108] With the initial condition {dot over (τ)}=1, the satellite time τ will be a linear function of time offset Δt, according to, $\begin{array}{cc}\left[\begin{array}{c}\tau \left(k\right)\\ \stackrel{.}{\tau }\left(k\right)\end{array}\right]=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right]\left[\begin{array}{c}\tau \left(k-1\right)\\ \stackrel{.}{\tau }\left(k-1\right)\end{array}\right].& \left(34\right)\end{array}$

[0109] In the absence of process noise, the equation (34) provides that,

[0110] τ(k)=Δt+τ(k−1).  (35)

[0111] Assuming the GPS time at arrival is accurate, the process noise covariance matrix is, $\begin{array}{cc}{Q}_L\left(\Delta t\right)=\left[\begin{array}{cccc}{Q}_p& Q_{\mathrm{pv}}& 0& 0\\ Q_{\mathrm{pv}}& Q_v& 0& 0\\ 0& 0& Q_{\mathrm{clock}}\left(\Delta t\right)& 0\\ 0& 0& 0& 0\end{array}\right],& \left(36\right)\end{array}$

[0112] where the submatrices Q_p,Q_pv,Q_v are components of the covariance matrix due to the position (p), velocity (v) or both (pv). The process noise covariance matrix Q_L(Δt) can be modified to include a nonzero covariance submatrix (as the leftmost and bottommost element) for the satellite time process noise, either for the purpose of including corrections to the satellite clock or for tuning the Kalman filter. Using Q_sat(Δt) as the appropriate submatrix, the process noise covariance matrix becomes, $\begin{array}{cc}{Q}_L\left(\Delta t\right)=\left[\begin{array}{cccc}{Q}_p& Q_{\mathrm{pv}}& 0& 0\\ Q_{\mathrm{pv}}& Q_v& 0& 0\\ 0& 0& Q_{\mathrm{clock}}\left(\Delta t\right)& 0\\ 0& 0& 0& Q_{\mathrm{sat}}\left(\Delta t\right)\end{array}\right].& \left(37\right)\end{array}$

[0113] The Kalman filter also uses an error covariance matrix, here designated as P(k), defined as,

[0114] P(k)=E└(x_L(k)−{circumflex over (x)}_L(k))(x_L(k)−{circumflex over (x)}_L(k))^T┘,  (38)

[0115] where x_L(k) is the true value of the sought-after state vector, and {circumflex over (x)}_L(k) is its estimate, and E[ ] denotes the mathematical (statistical) expectation of the bracketed quantity.

[0116] In using a Kalman filter, values must be provided for the measurement noise covariance matrix R associated with the measurement noise vector ε (see equation (29) above and its discussion following), the initial estimate P₀⁻ of the error covariance matrix, and an initial estimate {circumflex over (x)}₀⁻ of the state vector.

[0117] The Invention Applied to the Extended Kalman Filter

[0118] A Kalman filter solution according to the invention proceeds, in the preferred embodiment, is performed according to the following steps. First, compute the expected (predicted) pseudorange {right arrow over ({circumflex over (ρ)})}_k=[{circumflex over (ρ)}_1,k,{circumflex over (ρ)}_2,k, . . . ,{circumflex over (ρ)}_N,k]^t(for k=0 indicating the initial instant) based on the GPS satellite position and an a priori state estimate {circumflex over (x)}_k⁻ of the state vector:

[0119] {circumflex over (ρ)}_i,k=∥{right arrow over (s)}_i({circumflex over (τ)}_k⁻−{tilde over (t)}_𝑓⁽ⁱ⁾)−{right arrow over ({circumflex over (u)})}_k⁻∥+cΔ{circumflex over (t)}_u,k⁻,  (39)

[0120] for i=1,2, . . . , N satellites.

[0121] Next, construct the measurement connection matrix, $\begin{array}{cc}{G_k\left({\hat{x}}_k^{-}\right)=\left(\begin{array}{cc}\frac{\partial f_1}{\partial x}\frac{\partial f_1}{\partial y}\frac{\partial f_1}{\partial z}\frac{\partial f_1}{\partial t_{\mathrm{sat}}}& 1\\ \dots & \\ \frac{\partial f_N}{\partial x}\frac{\partial f_N}{\partial y}\frac{\partial f_N}{\partial z}\frac{\partial f_N}{\partial t_{\mathrm{sat}}}& 1\end{array}\right)|}_{{\hat{x}}_{\stackrel{\_}{k}}}.& \left(40\right)\end{array}$

[0122] in which 𝑓_iis the pseudorange function of equation (17), i.e.

[0123] 𝑓_i(x_u,y_u,z_u,τ,t_u)=∥{right arrow over (s)}_i({circumflex over (t)}_tx={circumflex over (τ)}−{tilde over (t)}_𝑓⁽ⁱ⁾)−{right arrow over (u)}(t_rx)∥+ct_u.

[0124] Then compute the gain matrix K_k, using,

[0125] K_k=P_k⁻G_k^t(G_kP_k⁻G_k^t+R_k)⁻¹,  (41)

[0126] and update the state vector and the covariance matrix according to,

[0127] {circumflex over (x)}_k⁺={circumflex over (x)}_k⁻+K_k({right arrow over (ρ)}_k−{right arrow over ({circumflex over (ρ)})}_k)={circumflex over (x)}_k⁻−K_kΔ{right arrow over (ρ)}_k

[0128] P_k⁺=(I−K_kG_k)P_k⁻(I−K_kG_k)^t+K_kR_kK_k^t  (42)

[0129] (The measured pseudorange can be determined from the submillisecond fractions measured by the receiver and the estimated millisecond fractions obtained from previous estimates, or from the location of a nearby cellular base station, or by some other method.)

[0130] Then propagate the covariance matrix to the next measurement epoch, using,

[0131] P_k+1⁻=ΦP_k⁺Φ^t+Q_k.  (43)

[0132] Finally, propagate the state estimate to the next measurement epoch using the assumed process dynamics, i.e. according to,

[0133] {circumflex over (x)}_k+1⁻=Φ{circumflex over (x)}_k⁺.  (44)

[0134] The Kalman filter procedure as defined allows changing what are used as the components of the state vector; it allows removing those state vector components, during the course of the Kalman filter solution, that become available through some external source. When a state vector component is removed, all the matrices used in the Kalman filter operation must be modified to correspond to the new state vector. For example, if at some stage the satellite time becomes available from baseband tracking, then the satellite time state vector components can be removed from Kalman filter, and the other vectors and the matrices used in the state update equations and the measurement update equations must be changed to correspond.

[0135] Scope of the Invention

[0136] It is to be understood that the above-described arrangements are only illustrative of the application of the principles of the present invention. In particular, instead of orbit parameters (of a satellite) being used in the receiver position estimation, information in some other form about the orbit of the satellite may instead be used, and in general information about the satellite may be used, information such as ephemeris or almanac information, provided such information includes or implies information about the satellite orbit. Numerous other modifications and alternative arrangements may be devised by those skilled in the art without departing from the spirit and scope of the present invention, and the appended claims are intended to cover such modifications and arrangements.