EP1949126A2

EP1949126A2 - Method and apparatus for reconstructing time of transit from assisted or weak signal gps observations

Info

Publication number: EP1949126A2
Application number: EP06788960A
Authority: EP
Inventors: Eamonn Glennon; Roderick Bryant
Original assignee: SigNav Pty Ltd
Current assignee: SigNav Pty Ltd
Priority date: 2005-07-29
Filing date: 2006-07-28
Publication date: 2008-07-30
Also published as: WO2007016474A2; EP1949126A4; CA2617142A1; WO2007016474A3; WO2007016474A9

Abstract

The invention provides a GPS device that resolves timing errors as applicable to an AGPS navigation solution. A single state algorithm provides rapid convergence, typical with standard Newton-Raphson performance. Problems due to ill-conditioning of the input matrix are identified by an augmented DOP calculation procedure, which may be used to identify unreliable results or to estimate the errors for those results.

Description

Method and Apparatus for Reconstructing Time of Transmit from Assisted or Weak Signal GPS Observations

This invention claims the priority of US provisional patent application 60/703,637 filed July 29, 2005 entitled "Solution of Timing Errors for AGPS".

Field of the Invention

This invention relates to Global Positioning Satellite systems where the time of transmit of satellite signals is not available from the navigation message data. In particular it concerns the reconstruction of time-of-transmit from the course acquisition code of weak signals or where assisted-GPS is employed.

Background of the Invention

Assisted GPS (AGPS) navigation solutions differ from normal GPS navigation solutions due to the use of ambiguous code-phases rather than full time-of-transmits for each GPS satellite observation. As such, it is necessary to reconstruct the full time-of- transmit using a-priori knowledge such as a position estimate within 150 km of the true position and an estimate of the time-of-receipt. Errors arise if the initial time-of-receipt used to construct the time-of-transmit observations is in error.

The Enhanced-E911 requirement for mobile cellular communications and the subsequent use of GPS in order to fulfill this requirement has necessitated the development of new methods to perform GPS navigation solutions. Unlike standard GPS, weak-signal or AGPS is not able to extract all the information from the GPS signal due to extremely weak signal to noise ratios. As a result, a weak-signal or AGPS satellite observation generally consists of a 1-ms ambiguous code phase and a measured Doppler frequency compared to a standard GPS observation which consists of a full time-of-transmit (TOT) and a measured Doppler frequency. Nonetheless, AGPS is still able to perform a navigation solution through the use of prior information, including a rough estimate of the position of the receiver and a time tag for the time-of-receipt (TOR). The rough estimate of position of the receiver is generally based on the location of the cell-site, although it could be based on the use of previous estimates or calculated using a Doppler based solution. These parameters can then be used to estimate ranges to each satellite which together with the initial TOR estimate can be used to estimate a full TOT for each satellite. However, since the initial TOR typically contains errors, the reconstructed TOTs will all be subject to a common timing error thereby resulting in navigation position errors in the solution process.

The problem of time-recovery for AGPS is discussed in J. Syrjarinne,

"Possibilities for GPS Time Recovery with GSM Network Assistance," presented at ION GPS 2000, Salt-Lake City, UT, 2000 and M. M. Chanarkar, "Resolving time ambiguity in GPS using over-determined navigation solution." United States of America: Sirf Technology, Inc, 2003. These references outline general algorithms for solution of the timing error through addition of an additional variable and its solution using least squares techniques. However those algorithms suffer from a poor rate of convergence requiring a large number of iterations in order to reach an acceptable solution. In some cases, up to 500 iterations were required; this being a significant problem when being implemented on a small embedded processor as employed in a typical mobile phone. In some scenarios the algorithm simply failed to converge and hence no solution was obtained. Furthermore, a weighting process employed to avoid poor timing error corrections for near zero range-rates, did not fully solve the problem.

Furthermore, those systems lack any predictive mechanism for determining when their algorithms are unreliable and thus leave the user uncertain whether or not to rely on their resolution of the time-of-transmit and accordingly on their position solution. Brief Description of the Invention

The invention is embodied in a GPS device having a processor that runs an algorithm for determining TOT, the time of transmit of a satellite signal. Inaccuracies in the receiver clock are the source of errors in the TOT. After adjusting the TOT by the phase shift between a locally generated C/A code and the satellite gold code, a calculation in a single stage is performed in which both the receiver clock error Δt and the TOT error ΔT. A measurement matrix is used whose row elements are the partial derivatives of the pseudorange with respect to the unknowns, in this case including an addition timing error term. The solution is obtained by applying the pseudoinverse of the measurement matrix to the difference between the predicted and measured pseudorange vector.

At the same time, the Dilution of Precision is calculated based on the measurement matrix. The DOP provides a prediction of the accuracy of the solution. Other embodiments take use the accurate time determination when time signals are required.

Brief Description of the Drawings

Figure 1 depicts the timing for a satellite signal. In the present invention. Figure 2 depicts a satellite configuration for an example used to demonstrate the invention.

Detailed Description of Preferred Embodiments

GPS observations used in AGPS Mobile Station Assisted (MS-Assisted) or MS- Based navigation solutions generally consist of a set of C/A code phases measured at a given time-instant as well as Doppler frequency measurements. Since the C/A code has a repetition frequency of 1 kHz, this means that the C/A code phases indicate the TOT modulo 1 ms and for a navigation solution to be performed, it is necessary to reconstruct the full TOT for each satellite. This process generally starts with an estimate of the TOR₁ the time instant at which the observation was made, which is generally obtained using assistance from the cell phone or real time clock (RTC). Subtracting an estimate for the satellite signal's time-of-flight (TOF) from the TOR then provides an estimate for each satellite TOT. The time of flight is essentially the phase shift in the C/A signal modulo the length of a C/A code epoch of 1 ms. Each TOT can then have it's sub- millisecond portion replaced by the measured code phase and the resulting TOTs can then be corrected by a integer number of millisecond epochs such that the times are consistent with the specified coarse receiver location. Boundary truncation or round-off problems can be avoided by adjusting the initial TOR such that after application of the measured code phases, at least one satellite has a millisecond adjustment that is exactly zero. This procedure produces a set of reconstructed TOTs that is consistent with the true location at which the measurement was made and that can be then used to perform a navigation solution.

As shown in Figure 1, a satellite 1 in an orbit 3 transmits a signal to a receiver 5 located on or near the surface of the earth. A portion of the signal, typically the start of an epoch, commences its transmission at a time TOT as measured by an accurate satellite clock, and commences its reception at a time TOR as measured by a less accurate receiver clock. When the signal arrives at the receiver it is detected by its correlation with a locally generated C/A code. The phase difference between the locally generated code and the received code translates into a time difference modulo the length of a code epoch. The epoch difference is constrained by knowledge of the approximate location 7 of the receiver to an accuracy of about 150 km, since the speed of light is approximately 300km/ms and the epoch length is 1 ms. A problem with the above procedure is that any error in the initial TOR will result in biased TOTs. Where this causes a problem is that each TOT is used to estimate the position of the corresponding satellite in its orbit 3 which means that if the TOT is incorrect then the satellite position will be incorrect as well. Since each satellite range rate lies in the range of ±1000 m/s, a timing error in the TOT of as little as 0.01 seconds can result in pseudorange predictions that are in error by 10 m, while errors of 1 second can result in 1000 m pseudorange errors. These errors then result in biases in the calculated positions.

The problem can be solved if there is redundancy in the number of satellite observations used to perform the navigation solution. In a standard GPS solution, four observations are used to solve for four different state variables, namely (x, y, z, cAt), where x, y and z are the WGS84 coordinates of the receiver, c is the speed of light and Δt is the clock error in the TOR. It should be emphasized that this clock error Ai in the TOR is different to the timing error ΔT that is present in the estimated TOTs when AGPS is carried out. In an AGPS solution, the dimension of the state vector is increased by addition of the timing error ΔTand an additional observation used to solve for this additional unknown.

The present invention is best understood by comparison with the following two stage algorithm. In the two stage algorithm, the pseudorange error for each satellite / is modeled as:

Λ, = IR₁ I = ^Ri - R, (2)

-Vfo-tf+CK-tf ⁺fo-'J¹

₌ {^χ, - ^χj^χ, -^χ)+ (y, -yjy, - ϊ)+ (^z, ^~ Φr - *)

where PR, is the pseudorange for satellite /, R,- is the range for satellite /, R₁ is the range-rate for satellite /^', Δt is the TOR clock error and ΔT is the timing error. R, is calculated using the satellite position at the current TOT x_; or (x_lt y,, Z₁) and the receiver position x or (x, y, z).

Given this model, the two stage algorithm used solved for (x, y, z, cΔt) using a standard GPS solution process after which equation (1 ) was used to create a set of equations permitting a value for /lTto be obtained. The calculated ΔT value would then be used to correct all of the TOTs₁ the satellite positions (x_it y,, Z₁) were then recalculated and the process repeated. There are a number of ways in which (1 ) can be used to solve for the value of ΔT.

One method is to substitute the calculated Δt value from the position solution thereby providing a set of five or more equations each with a single unknown, namely the ΔT quantity. Using these equations, five or more estimates ΔT, are obtained by re- writing (1) as:

_Aτ PR₁ -R₁ -CAt ₍₃₎

^R,

These values can then be averaged in order to obtain the final ΔT value. The problem with this approach is that the range-rate quantity can lie anyway in the range of -1000 m/s to 1000 m/s, including zero. When the range rate R₁ is close to zero, the

resulting ΔT, estimate can be poor and this can contaminate the final estimate AT . In order to fix this problem, an additional weighting process was introduced whereby rather than giving each AT₁ estimate equal weighting, as is the case for an average, the weighting depends on the value of R₁. Two different weightings were tested as given in (4) and (5).

_Δf ₌ y Δ^_ ₍₅₎

The effect of these weighted means is to give low weighting to estimates derived from range-rates close to zero and larger weights to values derived from higher range- rates. Estimates derived from values very near to zero were excluded from the solution completely.

An alternate method of performing the solution for AT is to recalculate both At and ΔT at each position solution step. This can be done by constructing a vector of Pf?;, R, and Jj₁. as PR, R and R respectively and rewriting (1) as

PR -R = [l R] ^CΔ1 (6) where 1 denotes a vector of 1's. This can then be used to solve for cΔt and AT using a standard pseudo-inverse. This method did not seem to be as susceptible to the problem of zero R₁ values as the previous technique, probably because of the presence of the cΔt variable constraining the solution in some way.

The iteration process between solving for a new ΔT value and recalculating the solution process using the calculated ΔT value continued until the magnitude of the applied corrections became reasonably small. Both of these techniques generally resulted in reasonable position solutions being obtained for a range of initial timing errors, although the methods were observed to suffer from a number of problems. One problem that was observed early in the process was a poor rate of convergence requiring a large number of iterations in order to reach an acceptable solution. In some cases, up to 500 iterations were required; this being a significant problem when being implemented on an small embedded processor as employed in a typical mobile phone. Furthermore, although the weighting process did improve the performance with near zero range-rates, it appeared as though this solution did not fully solve the problem. Finally, it was observed that occasional scenarios would be encountered in which the algorithm failed to converge and hence no solution was obtained. The combination of all these problems led to the development of the preferred embodiment of the present invention.

The single stage algorithm just described was developed as a minimum change addition to an existing process. This was achieved by simply nesting the new step within the existing process and iterating until convergence took place. Unfortunately good convergence did not always take place and it was conjectured that part of the reason for this was due to the fact that the solution process was being performed in two separated steps rather than a single step.

To overcome these defects, the entire solution process was modified so that all of the unknowns were updated at in a single operation. This required modification of the measurement matrix H so that rather than solving for (x,y,z,cAt) the process solved for (x, y, z, cΔt AT) instead. As is the case for a standard GPS solution, the rows of H are given as the partial derivatives of the PR₁ with respect to the unknowns, except in this case an additional timing error term ΔT is present. Row / of the H matrix is H_t given by

H₁ * l -A (7) where the- — '-, — '- and — '- contributions have been ignored since their magnitudes are dx dy dz small.

The signs of the partial derivatives have also been reversed to make the H matrix similar to the matrix typically used to calculate Dilution of Precision (DOP)₁ thereby permitting DOP to be calculated using the same matrix. (DOP is a term that describes the geometric strength of a satellite configuration. When visible satellites are close together in the sky, the geometry is said to be weak and the DOP value is high; when far apart, the geometry is strong and the DOP value is low. Factors that affect the DOP are, besides the satellite orbits, the presence of obstructions which make it impossible to use satellites in certain sectors of the local sky.) This reversal of signs requires that the solution process adjust the signs of the input quantities to compensate accordingly. It also makes sense to express the range-rate terms with units of km/s rather than m/s thereby ensuring that the magnitude of all matrix elements are similar. This results in the solved TOR error having units of milliseconds rather than seconds.

Defining y as the measured pseudorange vector, yas the predicted pseudorange vector given a particular estimate of the state vector x, an update to state Δx can be given by

where n is the number of satellites and the sign of Δy has been adjusted to match the DOP friendly H matrix. Following the calculation of each state correction Δx, the correction is applied and the predicted range vector y updated before repeating the process. The process is terminated when the magnitude of the calculated correction Δx becomes sufficiently small. This typically takes less than 10 iterations and is significantly faster than the two-stage algorithm.

The algorithm just described may be improved in a number of ways. Firstly, it is often the case that two-dimensional fixes need to be performed during periods of reduced satellite availability, so support for this requirement may be included. One way in which this can be done is to add a 'fake' satellite in the center of the earth and then specify a required pseudorange for that satellite. The corresponding row of the H matrix is then constructed as the partial derivatives of the WGS84 reference ellipsoid evaluated at the estimated position, where the ellipsoid is described as

R₀ = R(x,y,z) = τjx² +y² +(jz (9)

a 2 z

H_n+1 » , ' " 0 0

\ R[x,y,z) R{x,y,z) b² R{x,y,z)

(a/bf defines the WGS84 ellipsoid and Ro defines the altitude constraint. One problem with 2D solutions is that when the solution is over-determined, there is no hard constraint that forces the solution to the altitude specified by the user. To deal with this and also improve the solution process in general, a weighted least squares solution process is performed. Each observation is given a corresponding weighting related to the estimated RMS error that is reported with each observation. This also solves the 2D problem since the 2D observable can be given a significantly larger weighting than the other observables thereby establishing an effective altitude constraint. When weighting is added to the solution process, equation (8) is changed to

Λx = (H⁷R-¹Hy H^TR^~lAy (10) where R^"* is the weighting matrix. R is generally a diagonal matrix in which the noise- variance of each observation is the corresponding diagonal element in R, so observations with a low noise variance end up as large values in the inverted matrix R^'1.

Another feature of a preferred embodiment is to permit the use of differenced pseudoranges for solution of the (x, y, z, ΔT). This is implemented by modifying the H and R matrices once all the measured satellite observations have been inserted (but before adding the fake 2D observation). The modification involves selecting a reference satellite (usually the highest one) and then subtracting the reference satellite row from each of the other rows. Performing a differenced observation in this way eliminates the clock error cΔt as a solution state. It is also possible to disable the solution of the timing error for situations in which exact TOTs or accurate time-aiding is available omitting the final range-rate column of H thereby resorting to the standard formulation. This is also useful for permitting calculation of the receiver velocity solution since the H matrix required for a velocity solution is the same as a standard position solution. Only the input and output vectors change for this situation, which in this case are pseudorange- rates and receiver velocity respectively.

Although the single stage algorithm resulted in significantly better performance than its two-stage predecessor, it nonetheless was observed to occasionally suffer large position errors. These position errors occurred despite the use of observations obtained from observations with high signal to noise ratios and were more likely when the minimum number of observations were being used (i.e. five). To investigate the reasons for this, the measurements from such a failure were captured and the execution of the algorithm examined within a MATLAB environment. The cause of the failure was subsequently identified as being caused by the H matrix being severely ill-conditioned, i.e. small changes in matrix elements having drastic effects on the results of calculations. The unusual aspect of this situation was that although H was ill-conditioned, the DOP resulting from the matrix, although large in magnitude was not sufficiently large to explain the error. A complete end-to-end Matlab simulation including the generation of simulated measurements and solution using those measurements subsequently confirmed these captured measurement results.

The case study prompting this solution occurred with the following dataset given in Table 1 , although the effect is readily observable in other datasets. Table 2 shows the resulting satellite positions and range-rates for the constellation at this time, while Figure 1 shows provides a graphical representation of the satellite visibility.

Latitude S 35° 19' 50.84"

Longitude E 149° 11' 13.57"

Altitude (m) 600

X -4474390

Y 2668637

Z -3668214

Week/ Time-of-Week 1284/46663

Date / Time (UTC) 15 August 200412:57:43

Almanac Week/ Toa 260/319488

TABLE 2: CASE STUDY TIME AND LOCATION

TABLE 2: CASE STUDY SATELLITE POSITIONS / RANGE-RATES The problem occurs if the subset of the above constellation given by satellites 1 ,

11 , 16, 20 and 23 is used to perform the solution, with satellites 14 and 25 omitted from the solution process. The H matrix (in WGS84 XYZ coordinates) corresponding to the use of these observables is given by

If the standard GPS DOP is calculated by defining H as H with the last column

deleted and calculating (/Ftfj the result is

The GDOP can then be calculated as the square root of the trace of this matrix, giving a GDOP value of 7.17, which although higher than normal is not excessive. However, if the DOP is calculated using the H matrix as actually used in the solution process, a result with significantly larger magnitude is obtained.

This gives a GDOP of 492536, which is much larger than the standard method of 7.23, where the new augmented DOP calculation is the natural extension of the standard method. Hence the use of these observables clearly has a geometry making the solution with this set completely infeasible. These results were obtained from the simulation using almanac to define the satellite locations and velocities, although when the ephemeris from the captured log were used the resulting GDOP had a lower value of approximately 400. This difference is due to slight differences in the satellite positions between using the YUMA almanac and captured ephemeris. An examination of Figure 2 provides an intuitive reason for this failure, namely that saellitel and 20 are almost co- located which means that the geometry is effectively a 4 satellite solution even though 5 satellites are being used.

It should be noted that when reasonably good geometries are present, the numerical values for the extended DOP are generally similar to the values obtained using the conventional approach. However, in cases where the extended DOP is high, it means that the timing error and navigation solution cannot be reliably calculated.

Application to Weak Signal Timing

In addition to the application of the technique to E911 navigation solutions, the technology is also applicable to the newly developing field of GPS weak signal timing receivers. As the name suggests, this involves the use of GPS for time-transfer, except that the receiver only has access to weak signals and obtains ephemeris, coarse location and coarse time from other sources such as the Internet or wireless links. The main reason for requiring such receivers is either convenience, with many environments not being GPS antenna friendly or the requirement to provide timing in areas of heavy vegetation where GPS signals are obscured.

A number of different options are available to achieve this goal. The first option is to formulate the solution in terms of an extended Kalman filter (P. Axelrad and R. G. Brown, "GPS Navigation Algorithms," in Global Positioning Systems: Theory and Applications Volume 1, B. W. Parkinson, J. J. Spilker, P. Axelrad, and P. Enge, Eds.: American Institute of Astronautics and Aeronautics, Inc, 1996.), rather than the single shot solution algorithm just described. The procedure for doing this is straightforward and requires the inclusion of the TOR timing error /17 as an additional element in the Kalman filter state vector. The state covariance must also increase to account for the additional dimension and the new H matrix used instead of the normal H matrix. The Kalman filter can then be operated over a number of GPS measurements until the state- covariance associated with /47 has fallen significantly below one millisecond. Provided the Kalman filter state covariance represents a realistic estimate of the error it is then possible to correct the receiver TOR using the calculated estimate after which this TOR can be locked-in and a switch to a conventional solution process performed.

In addition to the use of the above process, it is also possible to take advantage of the fact that in practice a spread of satellite signal levels will generally be present. Even if the signal to noise ratio is too low to reliably extract data; those same signal levels may still permit bit-synchronization to be achieved (i.e. determine the locations of the bit boundaries). Provided that this can be done, then the ambiguity associated with each GPS observation improves from 1 ms to 20 ms. This is of benefit for a weak signal timing application where the two pieces of information may be combined to determine the TOR with no ambiguity using the following procedure.

Using the standard code phase ambiguity resolution, generate reconstructed measurements for the TOR algorithm. At this point, the bit-synchronization information for the satellite is not used.

Run the algorithm and obtain an estimate for the TOR error.

Apply the TOR error to all of the satellite observations and the incorrect TOR. If the TOR had no error then the C/A epoch calculated from the observation with bit- synchronization should match the measured epoch. If it does not match then apply an additional correction to .dTto ensure that a match occurs.

Additional Embodiments

As an alternative embodiment bit sync could be obtained by performing noncoherent correlations over many rounds (typically 50 to several hundred) of 20ms coherent correlations using 20 offsets at 1ms spacing and choosing the alignment that yields the highest correlation. In that case, the 1ms ambiguity of the codephase is replaced by a 20ms ambiguity for the bits. In addition the invention can be used to determine time to better than 10ms, by utilizing the TOR-resolving position-Time Kalman filter. By combining this better time resolution with the 20 ms bit ambiguity one can completely resolve the ambiguity leading to precise time resolution using the codephase measurements with no ambiguity.

The validity of the time resolution can be tested by performing long coherent correlations across many bit periods after stripping data known in advance. This will yield a very high correlation if the bit ambiguity has been properly resolved. If not, this will yield a very low correlation and one could re-resolve the ambiguity.

In another embodiment one may track time continuously by monitoring codephases, counting code epochs and taking account of measured Doppler offsets. The accurate time could be further utilized by outputting synchronization pulses at any desired repetition rate with sub-microsecond precision and stamping these with time via a communications port of some sort. Similarly one could discipline the receiver's reference oscillator by estimating the frequency bias of the oscillator using the Doppler measurements and the estimated time and position and steering the local oscillator to the correct frequency.

Although the invention has been discussed in terms of specific embodiments, persons of skill in this art will see its full utility. Accordingly the invention is intended to cover what is described in the following claims.

Claims

What is claimed is:

1. A GPS device having a processor that runs an algorithm for determining, the time of transmit of a satellite signal (TOT) comprising

adjusting the TOT by the phase shift between a locally generated C/A code and the satellite gold code,

calculating position solution in a single stage in which both the receiver clock error Δt and the TOT error ΔTare determined employing a measurement matrix whose row elements are the partial derivatives of a pseudorange with respect to the unknowns

extracting the solution by applying the pseudoinverse of the measurement matrix to the difference between the predicted and measured pseudorange vector.

2. The GPS device of claim 1 further including

calculating a Dilution of Precision (DOP) based on the measurement matrix,

predicting from the DOP a prediction of the accuracy of the solution.

3. The GPS device of claim 1 , wherein a pseudosatellite is hypothecated at the center of the earth and the measurement matrix elements are calculated at the partial derivatives of a reference ellipsoid evaluated at an estimated position of the device.

4. The GPS device of claim 1 wherein weighting for noise is added to the calculating process.

5. A GPS device having a processor that runs an algorithm for determining, the time of transmit of a satellite signal (TOT) comprising

adjusting the TOT by the phase shift between a locally generated C/A code and the satellite gold code, calculating position solution in a single stage in which both the receiver clock error Δt and the TOT error ΔTare determined employing a Kalman filter that includes a TOR timing error as an element in the Kalman filter state vector.

6. The GPS device of claim 1 , wherein the solution is used to provide a timing signal for steering a local oscillator to a correct frequency.

7. A method for application in a GPS device for determining, the time of transmit of a satellite signal (TOT) comprising

8. The method for application in a GPS device of claim 7 further including

calculating a Dilution of Precision (DOP) based on the measurement matrix,

predicting from the DOP a prediction of the accuracy of the solution.

9. The method for application in a GPS device of claim 7, wherein a pseudosatellite is hypothecated at the center of the earth and the measurement matrix elements are calculated at the partial derivatives of a reference ellipsoid evaluated at an estimated position of the device.

10. The method for application in a GPS device of claim 7 wherein weighting for noise is added to the calculating process.

11. A method for application in a GPS device for determining, the time of transmit of a satellite signal (TOT) comprising

calculating position solution in a single stage in which both the receiver clock error Δt and the TOT error ΔTare determined employing a Kalrnan filter that includes a TOR timing error as an element in the Kalman filter state vector.

12. The method for application in a GPS device of claim 7, wherein the solution is used to provide a timing signal for steering a local oscillator to a correct frequency.