Patent Publication Number: US-8525727-B2

Title: Position and velocity uncertainty metrics in GNSS receivers

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     Not applicable. 
     STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT 
     Not applicable. 
     BACKGROUND OF THE INVENTION 
     This invention is in the field of satellite navigation. Embodiments of this invention are directed to Global Navigation Satellite System (GNSS) receivers, such as of the Global Positioning System (GPS) type, and to improved estimation of position and velocity obtained by such receivers comprehending measurement uncertainties. 
     In recent years, satellite navigation equipment has become widespread. Sophisticated navigation systems for large-scale vehicles such as aircraft and ships, and positioning systems for map and survey functions, were the first implementations of GNSS technology. In recent years, however, GPS receivers have been implemented in relatively modest applications, including many consumer automobiles, add-on navigation systems for automobiles, and handheld navigation systems for hikers. Indeed, GPS systems are now included in golf carts to provide golfers with accurate distance information, and are now commonplace in mobile telephone handsets. 
     GNSS technology, currently implemented by way of the GPS system, has many advantages over previous navigation systems, particularly land-based systems such as LOng-RAnge Navigation (LORAN) systems. GPS signals can be used at all times of the day and night, and are not greatly affected by weather or atmospheric conditions. Line-of-sight transmission and worldwide coverage are facilitated by the number and positions of transmitting satellites in orbit. The GPS system also provides both position and velocity information to the user. Varying degrees of position and velocity accuracy are available, enabling a large number of users of inexpensive receivers, while still allowing sophisticated users to obtain a higher level of accuracy. 
     In a general sense, by way of background, GPS navigation is based on triangulation of the receiver location using signals received from multiple satellites of known location, in which the received signals are effectively time-stamped with the time of transmission. More specifically, GPS satellites periodically transmit a pseudo-random number (PRN) sequence via a spread-spectrum signal, in which the transmitted value of the frame boundaries of the data modulating the PRN sequence, as well as the PRN sequence itself, has a deterministic relationship to a known time (i.e., the beginning of the GPS “week”); as such, the received data from the satellites effectively include a time stamp. From the PRN value, GPS receivers determine the signal propagation time from the satellite to the receiver and, multiplying that propagation time by the speed of light, determine a measured distance (the “pseudorange”) of the receiver from each of the multiple satellites. Typically, the pseudoranges from four or more satellites of known position, in an earth-centered earth-fixed coordinate system, are then triangulated by solving a system of position equations, to determine the position of the receiver in that coordinate system. 
       FIG. 1  illustrates, in an idealized 2-D representation, the fundamental concept of GPS triangulation. In this illustration, each of satellites SAT 1 , SAT 2 , SAT 3  have transmitted a time-stamped signal that is received by a synchronized receiver at an unknown location. By subtracting the time of transmission indicated in the signal time stamp from the time of receipt, and multiplying that difference by the speed of light, the receiver can estimate its distance r 1 , r 2 , r 3  from respective satellites SAT 1 , SAT 2 , SAT 3  of known position within the coordinate system being used. Computational circuitry within the receiver can then identify point RLOC in that coordinate system at which all of the estimated distances r 1 , r 2 , r 3  coincide, to a best approximation. In the 3-D sense, each radial distance r 1 , r 2 , r 3  from satellites SAT 1 , SAT 2 , SAT 3  will of course define a sphere, requiring either four or more satellites or extrinsic information (e.g., knowing that location RLOC is on the surface of the earth) to resolve a unique receiver location RLOC. 
     As known in the art and as mentioned above, the geometric “range”, or distance, r from a GPS satellite to a GPS receiver corresponds to the propagation time between the transmission time T s  at the satellite and the receipt time T u  at the user, multiplied by the speed of light c:
 
 r=c ( T   u   −T   s )= cΔt  
 
The “pseudorange” that can be determined by the GPS receiver necessarily involves consideration of the time offset t u  between the time at the receiver and the true reference time (“system time”), and the time offset δt between the time at the satellite and system time. As such, the pseudorange ρ considering these time offsets can be expressed as:
 
                   ρ   =       ⁢     c   ⁡     [       (       T   u     +     t   u       )     -     (       T   s     +     δ   ⁢           ⁢   t       )       ]                   =       ⁢       c   ⁡     (       T   u     -     T   s       )       +     c   ⁡     (       t   u     -     δ   ⁢           ⁢   t       )                     =       ⁢     r   +     c   ⁡     (       t   u     -     δ   ⁢           ⁢   t       )                     
Or, in coordinate space, defining s as the vector from the origin (i.e., center of the earth, in an earth-centered coordinate system) to the satellite, u as the vector from the origin to the receiver location, and r as the vector from the receiver location to the satellite corresponding to the vector difference s−u:
 
ρ− c ( t   u   +δt )=∥ s−u∥ 
 
Pseudorange ρ and the time offset values are measurable or otherwise knowable, thus allowing solution of the actual range. And considering that the GPS network itself now determines corrections for the satellite clock offset δt from system time, and transmits those corrections to the satellites for rebroadcast to receivers, satellite clock offset δt is effectively compensated for at the receiver, such that:
 
ρ− ct   u   =∥s−u∥ 
 
     Conventional GPS systems determine the user position (i.e., receiver location RLOC in  FIG. 1 ) based on pseudorange measurements to four or more satellites, giving rise to a system of equations (one for each satellite):
 
ρ j   =∥s   j   −u∥+ct   u  
 
where j is the satellite index. For a system of four satellites (j=1, 2, 3, 4) in which satellite j is at position (x j , y j , z j ) in the coordinate system, the pseudoranges ρ j  can be expressed into a system of equations in which the unknowns include the user (receiver) position (x u , y u , z u ) and the user time offset t u :
 
ρ 1 =√{square root over (( x   1   −x   u ) 2 +( y   1   −y   u ) 2 +( z   1   −z   u ) 2 )}{square root over (( x   1   −x   u ) 2 +( y   1   −y   u ) 2 +( z   1   −z   u ) 2 )}{square root over (( x   1   −x   u ) 2 +( y   1   −y   u ) 2 +( z   1   −z   u ) 2 )}+ ct   u  
 
ρ 2 =√{square root over (( x   2   −x   u ) 2 +( y   2   −y   u ) 2 +( z   2   −z   u ) 2 )}{square root over (( x   2   −x   u ) 2 +( y   2   −y   u ) 2 +( z   2   −z   u ) 2 )}{square root over (( x   2   −x   u ) 2 +( y   2   −y   u ) 2 +( z   2   −z   u ) 2 )}+ ct   u  
 
ρ 3 =√{square root over (( x   3   −x   u ) 2 +( y   3   −y   u ) 2 +( z   3   −z   u ) 2 )}{square root over (( x   3   −x   u ) 2 +( y   3   −y   u ) 2 +( z   3   −z   u ) 2 )}{square root over (( x   3   −x   u ) 2 +( y   3   −y   u ) 2 +( z   3   −z   u ) 2 )}+ ct   u  
 
ρ 4 =√{square root over (( x   4   −x   u ) 2 +( y   4   −y   u ) 2 +( z   4   −z   u ) 2 )}{square root over (( x   4   −x   u ) 2 +( y   4   −y   u ) 2 +( z   4   −z   u ) 2 )}{square root over (( x   4   −x   u ) 2 +( y   4   −y   u ) 2 +( z   4   −z   u ) 2 )}+ ct   u  
 
Theoretically, of course, this system of four equations and four unknowns has a unique solution. In practice, however, each pseudorange measurement includes some amount of noise or error. In the presence of such noise, according to conventional techniques, this system of nonlinear equations can be solved, to derive user position (x u , y u , z u ) and user time offset t u , by way of iterative techniques such as least squares minimization applied to a linearization of these non-linear equations, or through the use of Kalman filtering over a time sequence of the pseudorange measurements.
 
     It is useful, for purposes of comprehension, to describe the operation of an example of the linearization approach to the solution of this system of equations. As known in the art, linearization can be carried out with reasonable accuracy from an approximate position ({circumflex over (x)} u , ŷ u , {circumflex over (z)} u ), for example as may be determined by cellular positioning or other known methods. In the absence of assistance (i.e., via “assisted GPS” as known in the art), the initial approximate position can be taken as the center of the earth; as such, the approximate position need not be very accurate. The offset of the true user position (x u , y u , z u ) from this approximate position can then be denoted by a displacement (Δx u , Δy u , Δz u ). Taylor expansion of each of the set of pseudorange ρ j  equations allows this displacement or position offset (Δx u , Δy u , Δz u ) to be expressed as linear functions of the known coordinates and the pseudorange measurements. With a pseudorange ρ j  represented by: 
                       ρ   j     =       ⁢             (       x   j     -     x   u       )     2     +       (       y   j     -     y   u       )     2     +       (       z   j     -     z   u       )     2         +     ct   u                   =       ⁢     f   ⁡     (       x   u     ,     y   u     ,     z   u     ,     t   u       )               ⁢                 
one can calculate an approximate pseudorange {circumflex over (ρ)} u  using the approximate position ({circumflex over (x)} u , ŷ u , {circumflex over (z)} u ) and a time bias estimate {circumflex over (t)} u :
 
                       ρ   ^     j     =       ⁢             (       x   j     -       x   ^     u       )     2     +       (       y   j     -       y   ^     u       )     2     +       (       z   j     -       z   ^     u       )     2         +     c   ⁢       t   ^     u                     =       ⁢     f   ⁡     (         x   ^     u     ,       y   ^     u     ,       z   ^     u     ,       t   ^     u       )                   
Each component of the true user position (x u , y u , z u ) and receiver clock offset t u  can be expressed as the sum of an approximate component and an incremental component:
 
 x   u   ={circumflex over (x)}   u   +Δx   u  
 
 y   u   =ŷ   u   +Δy   u  
 
 z   u   ={circumflex over (z)}   u   +Δz   u  
 
 t   u   ={circumflex over (t)}   u   +Δt   u  
 
in which case:
 
ƒ( x   u   ,y   u   ,z   u   ,t   u )=ƒ( {circumflex over (x)}   u   +Δx   u   ,ŷ   u   +Δy   u   ,{circumflex over (z)}   u   +Δz   u   ,{circumflex over (t)}   u   +Δt   u )
 
This function ƒ(x u , y u , z u , t u ) can be expanded about the approximate position and time bias estimate ({circumflex over (x)} u , ŷ u , {circumflex over (z)} u , {circumflex over (t)} u ) using a Taylor series and truncating the partial derivative terms after the first-order expressions, so that the pseudoranges ρ j  can be expressed as:
 
               ρ   j     =         ρ   ^     j     -           x   j     -       x   ^     u           r   ^     j       ⁢   Δ   ⁢           ⁢     x   u       -           y   j     -       y   ^     u           r   ^     j       ⁢   Δ   ⁢           ⁢     y   u       -           z   j     -       z   ^     u           r   ^     j       ⁢   Δ   ⁢           ⁢     z   u       +     ct   u             
where {circumflex over (r)} j  is defined as:
 
 {circumflex over (r)}   j =√{square root over (( x   j   −{circumflex over (x)}   u ) 2 +( y   j   −ŷ   u ) 2 +( z   j   −{circumflex over (z)}   u ) 2 )}
 
Introducing shorthand variables a xj , a yj , a zj  as follows:
 
               a   xj     =         x   j     -       x   ^     u           r   ^     j                     a   yj     =         y   j     -       y   ^     u           r   ^     j                     a   zj     =         z   j     -       z   ^     u           r   ^     j             
and defining pseudorange difference Δρ=ρ j −{circumflex over (ρ)} j , one can express the set of linear equations to be solved, based on range measurements to four satellites, as:
 
Δρ 1   =a   x1   Δx   u   +a   y1   Δy   u   +a   z1   Δz   u   −ct   u  
 
Δρ 2   =a   x2   Δx   u   +a   y2   Δy   u   +a   z3   Δz   u   −ct   u  
 
Δρ 3   =a   x3   Δx   u   +a   y3   Δy   u   +a   z4   Δz   u   −ct   u  
 
Δρ 4   =a   x4   Δx   u   +a   y4   Δy   u   +a   z4   Δz   u   −ct   u  
 
By expressing the coefficients and variables in matrix form:
 
                 Δ   ⁢           ⁢   ρ     =     [           Δ   ⁢           ⁢     ρ   1                 Δ   ⁢           ⁢     ρ   2                 Δ   ⁢           ⁢     ρ     3   ⁢                           Δ   ⁢           ⁢     ρ   4             ]       ;     H   =     [           a     x   ⁢           ⁢   1             a     y   ⁢           ⁢   1             a     z   ⁢           ⁢   1           1             a     x   ⁢           ⁢   2             a     y   ⁢           ⁢   2             a     z   ⁢           ⁢   2           1             a     x   ⁢           ⁢   3             a     y   ⁢           ⁢   3             a     z   ⁢           ⁢   3           1             a     x   ⁢           ⁢   4             a     y   ⁢           ⁢   4             a     z   ⁢           ⁢   4           1         ]       ;       Δ   ⁢           ⁢   x     =     [           Δ   ⁢           ⁢     x   u                 Δ   ⁢           ⁢     y   u                 Δ   ⁢           ⁢     z   u                   -   c     ⁢           ⁢   Δ   ⁢           ⁢     t   u             ]             
the system of equations can be expressed as:
 
Δρ= HΔx  
 
and, in theory, can be solved for the vector of position and time displacement Δx from the estimated position and time offset via:
 
Δ x=H   +   Δp  
 
where H +  represents a pseudoinverse of linearized derivative matrix H, calculated according to the particular solution technique (e.g., least-squares, weighted least-squares, etc.).
 
     In similar fashion, GPS navigation provides the capability of determining the three-dimensional user velocity {dot over (u)} as an approximation of the time rate of change of the calculated user position u over a time interval: 
               u   .     =         ⅆ   u       ⅆ   t       =         u   ⁡     (     t   2     )       -     u   ⁡     (     t   1     )             t   2     -     t   1                 
Doppler shifts in the frequencies of the received satellite signals are analyzed to solve the velocity problem, in many modern GPS receivers. In a linearization approach similar to that described above for determining the user position, a variable d j  is used for the j th  satellite signal:
 
 d   j   ={dot over (x)}   u   a   xj   +{dot over (y)}   u   a   yj   +ż   u   a   zj   −c{dot over (t)}   u  
 
For signals received from four satellites, a system of four equations in the four unknowns (i.e., the derivative terms {dot over (x)} u , {dot over (y)} u , ż u , {dot over (t)} u ) can be expressed in matrix form:
 
 d=Hg  
 
where g T =[{dot over (x)} u  {dot over (y)} u  ż u  {dot over (t)} u ], permitting solution of vector g from this system via:
 
 g=H   +   d  
 
     As mentioned above, in practice, these calculations by way of which GPS systems can determine a user position and a user velocity are corrupted by noise in the signal. Referring to user position, the measured pseudoranges ρ amount to a sum of the computed pseudoranges h(x), based on the true user position and time offset vector x, with a composite measurement noise vector v:
 
ρ= h ( x )+ v  
 
For the least-squares solution method, the linearized system of equations including the effects of noise is more accurately expressed as:
 
Δρ= HΔx+v  
 
and, in theory, can be solved for the vector of position and receiver clock bias Δx from the estimated position and time offset via:
 
 Δ{circumflex over (x)}=H   + Δρ
 
where
 
Δ {circumflex over (x)}=Δx+e =( x−{circumflex over (x)} )+ e  
 
{circumflex over (x)} being the estimated user position and receiver clock bias, and e representing the position-domain error caused by the measurement noise v in the pseudorange measurements, as applied to the derivative matrix H:
 
 e=H   +   v  
 
     Measurement noise vector v of the signals as received from the satellites depends on such factors as pure measurement noise, errors included within the applicable model considered by the receiver, and any unmodeled effects. In conventional GPS systems, to the extent that random variables are included within vector v, those random variables are typically modeled as uncorrelated with one another, and as having equal variances over satellite signals. 
     It is contemplated, however, according to this invention, that the equal variance assumption is invalid in many practical situations. For example, a GPS receiver located within an “urban canyon” will receive GPS signals from some satellites in highly attenuated form, if not completely blocked by the surrounding buildings. In this case, some of the satellite signals received by that GPS receiver may be reflected signals, or may be received over multiple paths because of the reflections (i.e., multipath signal transmission). According to this invention, it is contemplated that these suboptimal received satellite signals will have a measurement noise vector v that is time-varying due to the changing channel environment as the user moves through the physical environment. These signals will have unequal variances among one another, and will generally exhibit time-varying non-zero biases. As such, it has been observed, in connection with this invention, that treatment of the random variables in the measurement noise vector v as uncorrelated and having equal variances will lead to inaccurate position and velocity determinations. 
     BRIEF SUMMARY OF THE INVENTION 
     Embodiments of this invention provide a system and method of determining realistic uncertainty measurements for signals received from each of multiple navigation satellites, such as in the GPS system. 
     Embodiments of this invention provide such a system and method in which the uncertainty measurements are useful for calculating either or both the position and the velocity of the receiver. 
     Embodiments of this invention provide such a system and method in which multipath effects due to satellite signal reflection are accurately estimated and weighted in determining position and velocity. 
     Embodiments of this invention provide such a system and method that is compatible for use with Kalman filter-based solution techniques. 
     Embodiments of this invention provide such a system and method in which the uncertainty measurements for each satellite are useful in balancing the weight assigned to GPS versus other position-determining technologies in determining position and velocity. 
     Other objects and advantages of this invention provided by embodiments of this invention will be apparent to those of ordinary skill in the art having reference to the following specification together with its drawings. 
     Embodiments of this invention may be implemented into receiver circuitry for a mobile device, such as a mobile telephone handset, portable GPS receiver, automobile or other vehicle navigation system, and the like. The receiver circuitry includes a measurement engine that processes received GPS or GNSS signals from multiple satellites, and that provides separate independent estimates of pseudorange and “delta range” (i.e., time derivative of pseudorange) measurement noise variances for each of the received satellite signals. A position engine receives the pseudorange and delta range measurements, and the estimated measurement noise variances for each satellite, from the measurement engine, and processes those measurements and variance estimates to derive position or velocity results, or both, of improved accuracy, and uncertainty values corresponding to those position and velocity determinations that accurately reflect current signal conditions. 
     According to one aspect of the invention, the measurement engine determines the pseudorange measurement noise variance for each satellite signal based on a discriminator output signal in a delay-locked loop. The discriminator output determines the relative time position of the received demodulated satellite signal based on one or more correlations of that signal with a replica code. The correlations in the delay-locked loop are adjusted by the discriminator output, and a pseudorange measurement noise variance is determined based on the variance of that discriminator output signal over time. A delta range measurement noise variance for each satellite signal can also be determined from the variance in a discriminator output signal in a frequency-locked loop (or phase-locked loop), in which one or more correlations of the received signal with a replica code are performed. The pseudorange and delta range measurements, and the measurement noise variance for each of those measurements, are then used by the position engine in determining position, velocity, and uncertainties in those position and velocity determinations. 
     According to another aspect of the invention, the pseudorange measurement noise variance and/or delta range measurement noise variance are determined, by the measurement engine, in a manner that contemplates the effects of correlation bias due to multipath effects in one or more of the received satellite signals. 
     According to another aspect of the invention, the measurements of pseudorange, delta range, pseudorange measurement noise variance, and delta range measurement noise variance, for the signals from each of the satellites are applied by the position engine according to a Kalman filter approach, such that the determination of position and velocity, and their corresponding uncertainties, are obtained in a time-filtered manner. 
     According to another aspect of the invention, the position engine incorporates dynamic estimates of process noise, for example as based on user dynamics, in determining the position and velocity uncertainties. 
    
    
     
       BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING 
         FIG. 1  is a conceptual diagram illustrating the operation of conventional GPS navigation in a general sense. 
         FIG. 2  is an electrical diagram, in block form, of a navigation system constructed according to embodiments of this invention. 
         FIG. 3  is an electrical diagram, in block form, of an applications processor in the navigation system of  FIG. 2  according to embodiments of this invention. 
         FIGS. 4   a  through  4   c  are electrical diagrams, in block form, of a measurement engine and a channel within the measurement engine, in the system of  FIG. 3  according to embodiments of this invention. 
         FIGS. 5   a  through  5   e  are plots illustrating correlation plots of received satellite signals and various error-inducing effects in those correlation plots. 
         FIG. 6  is a flow diagram illustrating the operation of the measurement engine in the navigation system of  FIG. 2  according to embodiments of this invention. 
         FIGS. 7   a  through  7   c  are flow diagrams illustrating the operation of a position engine in the navigation system of  FIG. 2  according to embodiments of this invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention will be described in connection with its preferred embodiment, namely as implemented into a mobile Global Positioning System (GPS) receiver such as a personal navigation device, because it is contemplated that this invention will be especially beneficial in such an application, and in other similar handheld devices such as mobile telephone handsets with GPS capability, and the like. However, those skilled in the art having reference to this specification will readily comprehend that embodiments of this invention may be beneficially implemented in a wide range of Global Navigation Satellite System (GNSS) receivers other than in such handheld devices, including surveying instrumentation, navigation systems deployed in larger-scale vehicles and vessels, and the like. Accordingly, it is to be understood that the following description is provided by way of example only, and is not intended to limit the true scope of this invention as claimed. 
     Referring now to  FIG. 2 , the construction of an electronic system according to an embodiment of this invention will now be described, by way of navigation system  20 . Navigation system  20  may correspond to a stand-alone handheld or automotive navigation device (i.e., a “personal navigation device”), to a larger-scale navigation system in a large vehicle, and in any of these applications may be combined with other functions in a multipurpose device. For example, modern mobile telephone handsets referred to as “smartphones” often include GPS capability, and can thus serve as navigation system  20 . For purposes of this specification, it is contemplated that navigation system  20 , shown in  FIG. 2  by way of example, can correspond to any such realizations; of course, if navigation system  20  is a mobile telephone handset, it will also include such functionality necessary to carry out that and any other primary or ancillary functions. As such, the description of navigation system  20  in this specification is not intended to be interpreted in any limited sense. 
     As shown in  FIG. 2 , the digital functionality controlling the operation of navigation system  20  is centered on applications processor  22 , which is a conventional single or multiple processor core integrated circuit device, such as an OMAP4xxx applications processor or an AM35x ARM-based processor, both types of processors available from Texas Instruments Incorporated. As known in the art, such applications processors are capable of managing functions and applications including GPS position and velocity determinations, and can also manage other functions and operations as desired, such as mobile telephony, digital camera functions, data and multimedia wireless communications, audio and video streaming, storage, and playback, and the like. As such, applications processor  22  in navigation system  20  provides substantial computational power. As typical in conventional navigation systems, applications processor  22  includes or has access to program memory that stores computer-executable instructions and software routines that, when executed by applications processor  22 , carries out certain functions executed by navigation system  20  according to embodiments of this invention. Applications processor  22  also includes or has access to data memory that stores the results of its various processing routines and functions. Of course, other architectures and capabilities (lesser or greater) may be realized by applications processor  22  within navigation system  20  in this embodiment of the invention. 
     Applications processor  22  cooperates with various interface functions within navigation system  20 , some of which are shown in  FIG. 2 . Audio codec  21  manages the acquisition of audio input from microphone M and the communication of audio signals to applications processor  22 , and also the output of signals from application processor  22  as audio output via speaker S. Touchscreen interface  29  provides the appropriate graphics interface between applications processor  22  and touchscreen D, as well as receiving, processing, and forwarding user inputs entered via touchscreen D. 
     Navigation system  20  is powered from battery B, under the control and management of power management unit (PMU)  28 . PMU  28  may be realized as a separate integrated circuit (from applications processor  22 ), to control and apply the appropriate system power from battery B to applications processor  22  and the other functions within navigation system  20 , including any necessary regulated internal voltages. PMU  28 , as typical in the art, also can include power savings functionality, an interface to an on/off switch, a display of a remaining-battery-power indication, and the like. 
     Navigation system  20  according to this example includes inertial navigation functionality, to determine one or both of the user position and velocity according to measurements obtained by inertial measurement unit  30 . Inertial measurement unit  30  includes one or more transducers capable of sensing the position or movement of navigation system  20 . As shown in  FIG. 2 , inertial measurement unit  30  includes accelerometer  32   a  and magnetometer  32   b . As known in the art, accelerometer  32   a  detects physical movement of navigation system  20 , while magnetometer  32   b  can detect the earth&#39;s magnetic field and movement therewithin, in compass fashion. Other transducers such those useful to provide the function of an altimeter, gyroscope navigation, a humidity sensor, and the like can also be included in inertial measurement unit  30 , as known in the art. Conventional navigation functions and operations carried out by applications processor  22  based on such measurements can deduce the position of navigation system  20 . 
     Any or all of the various peripheral functions of navigation system  20  may be realized in one or more separate integrated circuits from the integrated circuit realizing applications processor  22 . Alternatively, some or more of these functions may be implemented into a single integrated circuit along with applications processor  22 . 
     As shown in  FIG. 2 , navigation system  20  according to embodiments of this invention includes Global Positioning System (GPS) receiver  33  that is coupled to antenna A to receive GPS signals from satellites in orbit, and that is coupled to applications processor  22  to forward those signals to applications processor  22  for navigation processing. GPS receiver  33  refers to one or more integrated circuits having interface circuitry to receive the incoming satellite signals, and having substantial computation capability for carrying out various operations involved in determining position and velocity based on those satellite signals. In this example, GPS receiver  33  includes multiple computational “engines” as realized by general-purpose processor capability such as that provided by dual-core ARM CORTEX A9MPCORE processors, such engines carrying out specific functions involved in the overall operation of GPS receiver  33 . An example of an integrated circuit known in the art for providing GPS and assisted GPS navigation and map functions is the NL5350 NAVILINK GPS device available from Texas Instruments Incorporated. It is contemplated that GPS navigation function  35  according to embodiments of this invention can be realized by way of a single integrated circuit with computational capability similar to that provided by that NL5350 device, and constructed and programmed to perform the operations described herein. The functions of certain ones of these computational engines directed to navigation functions, according to embodiments of this invention, will be described in further detail below. As shown in  FIG. 2 , GPS receiver  33  includes program memory  25   p , which stores computer-executable instructions and software (or “firmware”) routines that, when executed by GPS receiver  33 , carries out the functions executed by navigation system  20  according to embodiments of this invention. It is contemplated that program memory  25   p  will be realized as some form of embedded non-volatile memory such as electrically erasable programmable read-only memory (EEPROM), considering that the applications program instructions stored therein should persist after power-down of navigation system  20 . Of course, program memory  25   p  may alternatively be realized in other ways, for example by way of a memory resource external to the integrated circuit that realizes GPS receiver  33 , such as external non-volatile memory  31  in communication via memory interface  26 . GPS receiver  33  also includes data memory  22   d , either contained within the same integrated circuit as the processing circuitry of GPS receiver  33  or external thereto, for storing results of its various processing routines and functions. Data memory  22   d  may be realized as conventional volatile random access memory (RAM), or as non-volatile read/write memory such as EEPROM, or as some combination of the two. 
     In these embodiments of the invention, GPS receiver  33  executes program instructions, stored in program memory  25   p  or elsewhere, to calculate estimates of the geographical position and/or velocity of navigation system  20 . In addition, as known in the art, GPS receiver  33  can also or instead provide an accurate estimate of the current time from the received GPS satellite signals. The architecture and operation of GPS receiver  33  in determining these position and velocity measurements from the GPS system, according to embodiments of this invention, will be described in further detail below. 
     Referring now to  FIG. 3 , functions within navigation system  20  directed to navigation (i.e., determining one or both of position and velocity), according to embodiments of this invention, will be described in further detail. Not all functions performed by navigation system  20  are shown in  FIG. 3 , but rather only those functions concerned with the navigation function of embodiments of this invention are illustrated. As shown in  FIG. 3 , inertial measurement unit  30  communicates with applications processor  22 , specifically with inertial navigation system  44 , to provide inertial navigation system  44  with measurements from one or more of its inertial transducers (e.g., accelerometer, magnetometer, gyro, etc.). Inertial navigation system  44  evaluates these inertial transducer measurements in the conventional manner, to derive estimates of the current position of navigation system  20  and, perhaps also its current velocity. According to embodiments of this invention, inertial navigation system  44  forwards this position and velocity information to blending Extended Kalman Filter (EKF) function  45 , for combining with GPS position and velocity information to determine an integrated position and velocity result on output POS_INT. 
     The cooperation between applications processor  22  and GPS receiver  33  is also shown in  FIG. 3 . While the GPS navigation functions illustrated in  FIG. 3  are shown as being carried out within GPS receiver  33 , and thus typically realized in GPS receiver  33  by way of a separate integrated circuit from other functions of applications processor  22 , those skilled in the art will realize that the partitioning of the GPS functions between GPS receiver  33  and applications processor  22  may vary, depending on the particular realization. For example, the particular navigation functions shown as within GPS receiver  33  in  FIG. 3  may be alternatively realized by way of one or more specific processor cores within the integrated circuit embodying applications processor  22 , as may be the case. 
     As shown in  FIG. 3 , GPS receiver  33  includes analog front end (AFE)  38 , which receives and processes signals from GPS satellites, specifically from multiple GPS satellites in this implementation. AFE  38  thus implements such physical layer circuitry that interfaces with antenna A, including first stage amplification and the like. AFE  38  performs conventional downconversion of the received GPS signals, analog filtering, automatic gain control (AGC), analog to digital conversion, and other conventional analog domain processing. AFE  38  may also include some digital functionality, as known in the art. The received and processed signals from AFE  38  are coupled to measurement engine  40 , according to embodiments of this invention, for processing and computations in the digital domain on the processed signals received at antenna A from multiple GPS satellites. 
     Measurement engine  40  may be realized as customized hardware or logic, or as programmable logic executing computer program instructions stored in program memory  25   p  ( FIG. 2 ) or elsewhere, or some combination thereof, for carrying out the computations and processing of received digitized GPS satellite signals in the manner described in this specification. In particular, in the embodiment of the invention illustrated in  FIG. 3 , measurement engine  40  is operable to receive one or more digital data streams from AFE  38  including digital versions of the signals received from the GPS satellites, and to process those data streams into estimates of “pseudorange” and “delta range”, and estimates of measurement noise variances of the “pseudorange” and “delta range” values, for each of the individual received satellite signals. It is contemplated that those skilled in the art having reference to this specification, and particularly those portions of this specification that describe the functions and operations carried out by measurement engine  40 , will be readily able to realize measurement engine  40  in a manner suitable for their particular applications, without undue experimentation. 
     As shown in  FIG. 3 , measurement engine  40  produces signal groups  46   0  through  46   k , each associated with a corresponding received satellite signal, and each forwarded to position engine  42 . Typically, signals from four or more satellites are processed by measurement engine  40 , with the number of signal groups  46  varying accordingly. Each signal group  46  includes separate signals for the various measurements and metrics, according to embodiments of this invention. For example, as shown in  FIG. 3  relative to signal group  46   0  for a corresponding individual satellite and received signal, measurement engine  40  produces pseudorange signal ρ, delta range signal δρ (commonly referred to in the art by the symbol {dot over (ρ)}), pseudorange variance signal σ 2 (ρ), and delta range variance signal σ 2 (δρ). In this embodiment of the invention, each of signal groups  46   1 ,  46   2 , etc. include these same signals for a different received satellite signal (i.e., received from a different GPS satellite). While  FIG. 3  illustrates the signals within each signal group  46  to be communicated in parallel, and illustrates the signal groups  46  themselves as communicated in parallel, it is contemplated that the signals communicated by measurement engine  40  may be modulated, multiplexed, or otherwise combined with one another, or alternatively stored in data memory  25   d  or elsewhere, or communicated in some other known manner to position engine  42 . The representation of this signal communication in  FIG. 3  is presented in a functional manner only; the particular physical signal communication technique from measurement engine  40  to position engine  42  is not of particular importance. 
     Position engine  42  is another function within GPS receiver  33 , and receives signal groups  46  from measurement engine  40  at its inputs. As in the case of measurement engine  40 , position engine  42  may be realized as customized hardware or logic, or as programmable logic executing computer program instructions stored in program memory  25   p  ( FIG. 2 ) or elsewhere, or some combination thereof, for carrying out the computations and processing of the contents of signal groups  46  from measurement engine  40  in the manner described in this specification. In the embodiment of the invention illustrated in  FIG. 3 , position engine  42  is operable to receive pseudorange and delta range measurements, and measurement noise variances of these pseudorange and delta range measurements, for each satellite on an individual basis, and to compute one or both of the position and velocity of navigation system  20 , and uncertainty metrics for either or both of those position and velocity determinations. In some embodiments of this invention, position engine  42  may optionally provide measurement engine  40  with feedback, such as position bias signals {b i } indicative of a difference between measured and predicted positions of navigation system  20 . Position engine  42  also provides inputs to inertial navigation system  44 , to facilitate calibration of sensors in inertial measurement unit  30 ; in some cases, inertial navigation system  44  may provide inputs to position engine  42  for its calibration or for other purposes, as known in the art. The manner in which position engine  42  performs its operations, according to embodiments of this invention, will be described in further detail below. In any event, it is contemplated that those skilled in the art having reference to this specification will be readily able to realize position engine  42  in a manner suitable for carrying out the descried operations for their particular applications, without undue experimentation. 
     In the overall architecture of navigation system  20  shown in  FIG. 3 , the position and velocity determinations generated by position engine  42 , and the corresponding uncertainty metrics for those measurements, are forwarded to blending EKF  45  for combination with position (and optionally velocity) measurements based on the inertial transducer measurements described above. It is contemplated that, according to embodiments of this invention, the uncertainty metrics for position and velocity generated by GPS receiver  33  will assist in properly blending the measurements, for example by weighting the GPS position and velocity determinations according to the detected uncertainty level in those determinations. Blending EKF  45  produces the resulting integrated position and velocity results on line POS_INT, for forwarding to other functions within applications processor  22  or elsewhere in navigation system  20 . 
     This combination of the GPS and inertial position and velocity determinations is provided by way of example only. It is contemplated that navigation system  20  or such other navigation system incorporating embodiments of this invention may rely solely on GPS navigation for its position and velocity determinations, and as such inertial measurement unit  30  and inertial navigation system  44  are entirely optional insofar as this invention is concerned. In such a case, the uncertainty metrics for GPS position and velocity are still useful, for example to advise the user of the accuracy of the GPS results. However, it is contemplated that the uncertainty metrics for GPS position and velocity as derived according to embodiments of this invention are especially useful in the blended application shown in  FIG. 3 . 
     The functional architecture of measurement engine  40  according to embodiments of this invention will now be described in connection with the generalized block diagram of  FIG. 4   a . For purposes of this description, it will be assumed that GPS receiver  33  also includes the appropriate functionality (not shown) for acquiring the incoming satellite signals; as such, the description of the operation of measurement engine  40  will be based on that assumption. The received satellite signals, after processing and filtering by AFE  38 , are applied to inputs of multiple digital measurement engine channels  65 , which are arranged in parallel. Each digital measurement engine channel  65  is associated with one of the k satellite signals and, among other functions, will strip the satellite-specific spreading code from the combined signal stream to recover its corresponding satellite signal contents and convert those contents to baseband. 
     The generalized functional construction of one of digital measurement engine channels  65 , is also illustrated in  FIG. 4   a . It is to be understood that the other digital measurement engine channels  65  within measurement engine  40  will be similarly arranged. The construction illustrated in  FIG. 4   a  is a functional architecture. As such, each digital measurement engine channel  65  may be physically constructed as custom logic circuitry within an integrated circuit, or as programmable logic (e.g., programmable general purpose logic such as a microprocessor, application-specific programmable logic such as a digital signal processor, etc.) executing program instructions stored in program memory, or as a combination of the two approaches (i.e., some portions realized in dedicated logic circuitry and other portions as program instructions executable by programmable logic circuitry). It is contemplated that those skilled in the art can readily realize the functions of these multiple digital measurement engine channels  65  in a manner best-suited for a particular application. Alternatively, it is contemplated that a single digital measurement engine channel  65  could be provided, with that single channel  65  sequentially processing the received satellite signals from each of multiple satellites; the following description will refer to the case in which multiple digital measurement engine channels  65  are used, by way of example. 
     Digital measurement engine channel  65 , receives the incoming satellite signals at its downconverter  63  of measurement engine  40 , which down-converts the frequencies of these signals to an intermediate frequency (i.e., at an IF below the received RF frequency, but still above baseband), with separate in-phase and quadrature-phase components resulting from the IF demodulation. Each digital measurement engine channel  65   i  includes its own IF downconverter  63 , because the signals received from each satellite can exhibit very different Doppler frequencies. Each digital measurement engine channel  65   i  operates on the IF signal from its downconverter  63  to produce a pseudorange signal ρ, delta range signal δρ, pseudorange variance signal σ 2 (ρ), and delta range variance signal σ 2 (δρ), based on that signal from satellite i. Digital measurement engine channel  65   i  includes pseudorange computation function  86 , delta range computation function  88 , and variance computation function  90 , each of which compute components of the signals output by digital measurement engine channel  65   i ; alternatively, multiple instances of measurement engine channel  65 , may share pseudorange computation function  86 , delta range computation function  88 , and variance computation function  90 , depending on the architecture of measurement engine  40 . 
     In the arrangement shown in  FIG. 4   a , digital measurement engine channel  65 , includes delay-locked loop  70 , which receives the IF in-phase and quadrature-phase data streams I,Q, and demodulate these data to baseband by way of orthogonal pseudo-noise (PN) spreading codes, and from the demodulated data determine a correlation peak time τ PK  indicative of the time-of-arrival of the satellite signal, from which pseudorange computation function  86  determines pseudorange signal ρ for satellite i. Similarly, in this arrangement, frequency-locked loop  71  determines a Doppler phase shift measurement ƒ PK  from IF signal I,Q, from which delta range computation function  88  determines delta range signal δρ for satellite i. Frequency-locked loop  71  also produces a feedback signal on line CARR_PH_ADJ to IF downconverter  63 , by way of which the carrier phase used in the IF downconversion can be incremented in response to sensed Doppler frequency shifts. 
     In a general sense, digital measurement engine channel  65   i  performs a correlation between the received satellite signal and the known spreading code for that signal, according to both a code-phase hypothesis and also a carrier-phase hypothesis. In the arrangement of  FIG. 4   a , delay-locked loop  70  determines a correlation value using the code-phase hypothesis of an estimated phase relationship between the known spreading code and the received signal; by way of an feedback loop, this code-phase hypothesis is adjusted according to the correlation result, to minimize an error value between the code-phase estimate of the hypothesis and the code-phase that produces the maximum correlation. That code-phase of the correlation peak indicates the time of receipt of the satellite signal at navigation system  20 , from which the user position can be determined. Similarly, frequency-locked loop  71  determines a correlation using the carrier-phase hypothesis of an estimate of the carrier frequency (or phase) of the received satellite signal (the received signal possibly including Doppler shift due to motion of navigation system  20  relative to the satellite), as applied by downconverter  63 . Feedback generated by frequency-locked loop  71  adjusts the carrier frequency estimate based on the correlation, to minimize an error value between the measured correlation and a correlation peak frequency; the correlation peak frequency is indicative of the instantaneous frequency of the received satellite, from which the user velocity can be determined. 
     According to some embodiments of this invention, digital measurement engine channel  65   i  also includes variance computation function  90 , which calculates pseudorange variance signal σ 2 (ρ), and delta range variance signal σ 2 (δρ) accordingly, for the signal received from satellite i. In some embodiments of this invention, variance computation function  90  determines these variances from a measure of the signal-to-noise ratio (or carrier-to-noise ratio, or some other noise measurement), as performed by SNR measurement function  72 . For example, SNR measurement function  72  can determine the signal-to-noise ratio (or carrier-to-noise ratio) in the conventional manner from the inputs and outputs of the discriminator functions within delay-locked loop  70  and frequency-locked loop  71 . In other embodiments of the invention, variance computation function  90  bases those variances on discriminator outputs D i , D i,f  from delay-locked loop  70  and frequency-locked loop  71 , respectively. In any event, these output signals representing pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) are forwarded on to position engine  42 , as described above. 
       FIG. 4   b  illustrates the generalized construction of an example of delay-locked loop  70   i , within representative digital measurement engine channel  65   i  for the single “channel” associated with the signal from satellite i; additional channels are provided for signals from others of the k satellites, as discussed above. In-phase and quadrature signals I, Q produced by IF downconverter  63  are applied to inputs of correlators  74 E,  74 L,  74 P, which correlate signals I,Q with a replica of the spreading code at baseband, as generated by code generator  82 . Correlators  74 E,  74 L,  74 P correlate the incoming signal with respective spreading code replicas C E , C L , C P , which are at different delay times relative to one another. Prompt correlator  74 P applies the replica code at a delay time (or phase) that corresponds to a recent prediction of the correlation time τ PK  at which the correlation result R(τ) is a maximum, for example as predicted by a prior iteration. That predicted correlation time τ PK  corresponds to the predicted time of arrival of the satellite signal. “Early” correlator  74 E applies the spreading code replica at a delay time that is in advance of the predicted peak correlation time (e.g., time τ PK −d/2, or ½ “chip” early), while “late” correlator  74 L applies the spreading code replica at a delay time after the predicted peak correlation time (e.g., time τ PK +d/2). Correlators  74 E,  74 L,  74 P generate correlation results E, L, P, respectively, which correspond to the correlation between the replica code and the early, late, and prompt replica code versions C E , C L , C P , respectively. These correlation results are forwarded to discriminator function  75 . 
     As known in the art, multiple correlation results at different delay times are useful in improving the accuracy of the estimate the time of arrival of the satellite signal, particularly in providing feedback control of the delay times of the correlations. In this example, discriminator function  75  receives these multiple correlation results, and from these results effectively determines an error between the correlation peak as measured by those results, and a predicted correlation peak, for example as based on prior correlations in earlier time intervals. In this example, discriminator function  75  has inputs receiving early and late correlation results E, L from correlators  74 E,  74 L, respectively (and, optionally, result P from correlator  74 P), upon which its error measurement is determined and presented at its output DIFF. Time-alignment of the replica code applied in correlation can then be adjusted in response to that error, so that the measured correlation peak tracks the actual peak time, allowing accurate derivation of the time of receipt of the satellite signal, and thus accuracy in the measured pseudorange ρ for that signal. 
     In the arrangement of  FIG. 4   b , filter function  80  receives output DIFF from discriminator  75 , and generates filter output D i  for application to variance computation function  90 . Filter function  80 , or alternatively some other function within delay-locked loop  70   i , also produces a signal τ PK  that indicates the current estimate of the location of the correlation peak, which approximates the time-of-arrival of the satellite signal. Code generator  82  produces replica code versions C E , C L , C P  for application to correlators  74 E,  74 L,  74 P, at delay times adjusted in response to current correlation peak estimate signal τ PK . As mentioned above, this feedback control allows these correlations to lock onto and track the time of arrival. 
     Other approaches or modifications to determining the time of arrival can alternatively be used within digital measurement engine channel  65   i  in embodiments of this invention. For example, the correlation results can be obtained by coherently adding multiple correlations within the time interval, effectively correlating to a longer PN sequence consisting of concatenated copies of the actual PN sequence. Non-coherent combining of correlations may also be used, to improve the SNR, in the conventional manner. 
     In the operation of the discriminator-based arrangement of  FIG. 4   b , the correlation used in determining the time of arrival is generally based on the power of the product of the received signal with the replica signal (i.e., spreading code), at varying delay times τ of the replica signal.  FIG. 5   a  illustrates an example of the correlation characteristic R(τ) over correlation delay time τ in an ideal case, in which the maximum correlation value R(τ) appears at correlation delay time τ=0 (i.e., τ PK =0). The peak correlation time provides a best estimate the time of arrival of the signal. The values of correlation characteristic R(τ) at early delay times (τ&lt;0) and late delay times (τ&gt;0) are below the peak value at the “prompt” delay time τ=0 in this ideal case. Because of the ideal alignment as shown in  FIG. 5   a , the early correlation value E=R(τ=−d/2) and the late correlation value L=R(τ=+d/2) equal one another. 
     As mentioned above, discriminator-based delay-locked loop  70  adjusts the code phase at code generator  82  in response to the feedback of discriminator output D i  to move its prompt correlation P to the actual peak correlation delay time.  FIG. 5   b  illustrates an example of a non-ideal case in which the current estimate of the peak correlation time is in error, in this case earlier than the true correlation peak. In this example, correlation time error e t  between the assumed prompt point P and the true correlation peak is about −d/4 (i.e., early by ¼ chip); delay-locked loop  70   i  is thus operating under the erroneous assumption that the correlation characteristic is following plot  81  of  FIG. 5   b , while the actual correlation characteristic of the received signal follows plot  83 . Because of the error in the current estimate, the prompt correlation value P=R(τ=−d/4) on plot  83  is significantly lower than that which would result if the current estimate were correct (i.e., at expected correlation value P′ on curve  81 ). Discriminator function  75  thus has the task of shifting the delay times of replica code versions C E , C L , C p  to properly align with the actual correlation characteristic of curve  83 . 
     Discriminator function  75  in delay-locked loop  70 , of  FIG. 4   b  accomplishes this by deriving a measure of the error in correlation time, and applying that measure, via filter function  80 , as feedback to code generator  82 . Various discriminator techniques are known in the art, and may be used in discriminator function  75 . In this example, discriminator function  75  compares the early and late correlation values E, L, respectively, to generate difference signal DIFF. In this example, difference signal DIFF corresponds to the late correlation value R(τ−d/2) subtracted from early correlation value R(τ−d/2), where τ is the current assumed correlation peak location. Filter function  80  smoothes variations in the difference signal DIFF over time in producing discriminator output D i , and may use prompt correlation value P determined by prompt correlator  74 P in that filtering. For example, filter function  80  may operate so that discriminator output D i  averages: 
               D   i     =       E   -   L       2   ⁢   P             
over time. According to this example, and in general, the feedback to code generator  82  in the form of discriminator output D i  is proportional to the correlation time error e t  ( FIG. 5   b ), i.e., D i =G*e t , where G is the loop gain.
 
     Referring back to  FIG. 4   a , frequency-locked loop  71  operates in a similar fashion as delay-locked loop  70  to synchronize the correlation of a replica code with changes in frequency, specifically to identify and track Doppler shifts in frequency caused by navigation system  20  moving at some velocity relative to the position of satellite i, and accounting for movement of the satellite itself. In the manner known in the art, frequency-locked loop  71  modulates correlator frequencies to maximize a correlation value R(ƒ) as a function of frequency (or, in the instantaneous sense, phase). The frequency at which this correlation peak occurs can then be used by position engine  42  to determine a current user velocity. Those skilled in the art having reference to this specification will recognize that various approaches are available for realization of frequency-locked loop  71   i , including by way of a phase-locked loop (PLL) approach rather than a frequency-locked loop (FLL). 
     In summary,  FIG. 4   c  illustrates a generalized construction of frequency-locked loop  71   i , according to one conventional approach for this function. In this example, correlator  76  receives in-phase and quadrature phase components I,Q of the received satellite signal associated with digital measurement engine channel  65 , from downconverter  63 , along with prompt code version C P  from code generator  82  in delay-locked loop  70   i  (to eliminate the need for a redundant code generator in frequency-locked loop  71   i ). Alternatively, correlator(s)  76  may receive early and late code versions C E , C L  as inputs. In the example of  FIG. 4   c , correlator  76  performs a correlation of prompt code version C p  with the downconverted received signal in the conventional manner, and produces in-phase and quadrature-phase correlation results IPS, QPS that are forwarded to discriminator function  77 . Discriminator function  77  operates in the conventional manner to determine a correlation value R(ƒ) based on the correlator results IPS, QPS. Alternatively, discriminator function  77  may receive one or more correlation results E, L, P from delay-locked loop  70   i , so that its operation is based on correlation hypotheses determined by delay-locked loop  70   i  rather than its own. Correlation value R(ƒ) is filtered by filter function  78 , in similar manner as filter function  80  in delay-locked loop  70   i , producing in this case frequency peak estimate ƒ PK  indicative of an estimate of the code frequency at which the peak correlation by correlator  76  is obtained. Filter  78  also produces discriminator output D i,f  that is forwarded to variance computation function  90  ( FIG. 4   a ). In this example, the frequency peak estimate ƒ PK  is feedback to downconverter  73  to adjust the carrier frequency it uses in downconversion, via signal CARR_PH_ADJ, and is also forwarded to delta range computation function  88  for use in determining a current delta range value δρ from satellite i, from which a velocity can be derived by position engine  42 . As a result of this construction, frequency-locked loop  70   i  operates to adjust the carrier frequency in order to attain maximum correlation with the received satellite signal, to sense variations in the peak carrier frequency due to the Doppler effect, such variations indicative of the velocity and direction of navigation system  20 . 
     As known in the art, navigation system  20  is sometimes in a location in which multipath signals from a given satellite i are received, often due to reflections from buildings or other manmade or natural structures. As fundamental in the art, these reflected versions of the satellite signal arrive at different times from the line-of-sight (LOS) signal; indeed, in some situations, navigation system  20  may not receive any LOS signal from a given satellite i, but may only receive reflected versions of the signal. These multipath reflections render quite difficult the determination of a time of arrival of the satellite signal, particularly by causing a “bias” error in the received signal and its correlation. More specifically, error in the calculated pseudorange ρ, in this multipath situation can be quite high, even if the signal-to-noise ratio of the satellite signal is itself very high. 
       FIG. 5   c  illustrates an example of the multipath problem, in a situation in which both a line-of-sight (LOS) signal and a reflected version of the signal are received. Plot  87  represents the correlation value C(τ) of a received satellite signal including a line-of-sight signal and a single reflection. Plot  88  is the correlation plot that would result from the LOS signal alone, and plot  89  is the correlation plot that would result from the reflected signal alone. In this situation, one could express correlation value C(τ) as:
 
 C (τ)= R (τ)+α R (τ+ B )
 
where B is the bias, or delay of the reflected signal relative to the LOS signal, and where attenuation α is given by:
 
α= m·e   jθ 
 
θ being the phase shift of the reflection (which is 0° in this case). In general, correlation value C(τ) is time-varying because the parameters B, m, and θ are time-varying, particularly as navigation system  20  moves through the reflective environment. As evident from  FIG. 5   c , the reflected signal causes an error in the measured pseudorange.
 
     The correlations illustrated in  FIG. 5   c  are somewhat idealized, as they assume infinite signal bandwidth. In practice, the signal bandwidth is limited, and as such the correlations are less sharp and can include additional correlation peaks.  FIG. 5   d  illustrates correlation plot  90  for the received signal, based on the sum of an LOS signal (correlation plot  92 ) and reflected signals (correlation plot  94 ) under the same conditions (same values of parameters B, m, and θ) as in  FIG. 5   c , but with a more typical signal bandwidth as encountered in practice. This limited bandwidth causes a substantial increase in the pseudorange error, from about +4 meters to about +58 meters. 
     Phase modulation in the reflections of the signal is an additional source of error in the pseudorange measurement, as evident in the correlation plots of  FIG. 5   e . In this example, in which a phase shift of π is involved in the reflected signal, but with the same parameters B, m as in  FIGS. 5   c  and  5   d , two peaks become evident in the correlation plot  97  based on the received LOS signal (correlation plot  98 ) and reflections (correlation plot  99 ). With only the change in the phase shift between the situations of  FIGS. 5   d  and  5   e , the pseudorange error changes polarity and remains substantial (from about +58 meters to about −59 meters, in other words changing by about 117 meters). 
     Those skilled in the art having reference to this specification will thus realize that many possible causes for navigation error exist. High levels of noise in the received satellite signal of course reduce fidelity in the estimated pseudorange and delta range determinations. Even with good signal-to-noise ratios, however, multipath signals can cause substantial bias in the pseudorange and delta range measurements, as described above in connection with  FIGS. 5   c  through  5   e . As mentioned above, conventional approaches to GPS and GNSS navigation often include an uncertainty measurement that is based on the SNR of the received signals. However, it has been observed, in connection with this invention, that this method of determining uncertainty often provides poor results, especially in downtown urban areas in which the multipath situation is often present. Specifically, use of the SNR in the multipath context can be especially misleading, as a strong signal (high SNR) can still be heavily biased by reflections, as described above, and thus produce pseudorange and delta range measurements that are in fact quite uncertain. 
     Discriminator-Based Measurement Noise Variance Determination 
     Measurement engine  40  according to embodiments of this invention can more precisely determine measurement noise in the received GNSS signals. The operation and functionality of measurement engine  40 , and specifically the operation of digital measurement engine channels  65  in measurement engine  40 , according to embodiments of this invention, will now be described in detail with reference to the process flow diagram of  FIG. 6 . In process  50 , measurement engine  40  receives digital data corresponding to signals received from multiple satellites, specifically those satellites from which navigation system  20  is receiving encoded GPS signals. For purposes of this description, each instance of process  50  involves the receipt of those encoded GPS satellite signals over a given time interval t. The duration of that time interval t varies with the particular GPS system, and can vary from as short as one millisecond to as long as one second or even longer. 
     Once the satellite signals are received for the current time interval t, processes  52  and  55  are performed to determine the pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) for each of the k satellites. For purposes of this description, the measured pseudorange ρ and delta range δρ are determined by measurement engine  40  in process  52 , while the pseudorange variance σ 2 (ρ) and delta range variance σ 2 (δρ) corresponding to those measurements are determined in process  55 , essentially simultaneously or in parallel with process  52 . For purposes of this description, and consistent with the parallel arrangement of digital measurement engine channels  65  ( FIG. 4 ), processes  52 ,  55  are contemplated to be performed for all of the k satellites simultaneously, and in parallel. Of course, if the architecture of measurement engine  40  may alternatively be arranged so that the pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) are sequentially determined for each of the k satellites. 
     Referring first to process  52 , each of digital measurement engine channels  65  in measurement engine  40  derives the measured pseudorange ρ and delta range δρ for its associated satellites, from the signals received in process  50  over the current interval t. The particular method used by measurement engine  40  to carry out process  52  can be any one of the available pseudorange and delta range approaches known to those skilled in the art; the approach summarized above relative to  FIGS. 4   a  and  4   b  including delay-locked loop  70  and frequency-locked loop  71 , and the corresponding computation functions  86 ,  88 , is useful in these embodiments of the invention. Process  52  thus yields a pseudorange ρ value and a delta range δρ for value from for each of the k satellites from which navigation system  20  has received a signal. 
     In parallel (or sequentially, as the case may be), pseudorange ρ and delta range δρ measurement noise variances are calculated in process  55  for each of the k satellites, by their respective digital measurement engine channels  65 . For purposes of this description, the operation of variance calculation process  55  will be described with reference to one satellite signal (from satellite i), it being understood that the parallel digital measurement engine channels  65  will be operating simultaneously on the signals from their respective satellites to also determine the corresponding variances. Of course, the particular sequence in which the pseudorange ρ and delta range δρ measurement variances are calculated for each of the k satellites can be arranged as best suited for a particular application, by those skilled in the art having reference to this description. 
     In a general sense, the result of variance calculation process  55 , over the entire set of k satellites, represents the measurement variances for a vector of measurement noise v, in the form of a covariance matrix Σ ν . As described above, for the example of a least-squares solution of the pseudorange ρ value, the position-domain error vector e caused by the measurement noise vector v in the pseudorange measurements, is expressed by:
 
 e=H   +   v  
 
The statistical characteristic of this position-domain error vector e is, as mentioned above, best represented by covariance matrix Σ e . If the solution is obtained by way of a least-squares method, for a covariance matrix Σ ν  of measurement noise v, the position-domain error covariance matrix Σ e  is given by:
 
Σ e   =H   + Σ ν ( H   + ) T  
 
The use of a covariance matrix Σ ν  of measurement noise allows a generalized approach to determining position-domain uncertainty in the case in which measurement noise variances of the various satellite signals differ from one another. Conventional uncertainty determinations often express position-domain error covariance matrix Σ e  based on the assumption that measurement noise v has equal variances over all satellites (i.e., Σ ν =σ 2 I), which results in the special case of the above equation:
 
Σ e =( H   T   H ) −1 σ 2  
 
According to embodiments of this invention, position-domain and velocity-domain uncertainties are computed separately for each satellite signal, and thus the assumption Σ ν =σ 2 I is not applied. Rather, a measurement noise covariance matrix Σ ν  in which measurement noise variances σ i   2  differ among the k satellites is used in the position and velocity determination. According to these embodiments of the invention, it is assumed that the elements of measurement noise vector v are statistically independent from one another, which renders a measurement noise covariance matrix Σ ν  in the form:
 
Σ ν =diag(σ 1   2 ,σ 2   2 , . . . , σ n   2 )
 
As will be described below, the individual satellite measurement noise variances σ i   2  are to be determined.
 
     According to embodiments of this invention, multiple methods are available for estimation of these measurement noise variances σ i   2 . These methods are illustrated in  FIG. 6  by processes  60   a ,  60   b . Of course, more than two processes  60   a ,  60   b  may be made available, if desired. In the embodiment of the invention shown in  FIG. 6 , either or both of these available processes  60   a ,  60   b  are selected for use based on attributes of the received signal from satellite i. It has been observed that a reasonable metric for determining the best method or methods for measurement noise variance is one based on a signal-to-noise ratio of the received signal. As such, measurement variance calculation process  55  begins with process  54 , in which a metric of the signal-to-noise ratio of the signal from one satellite (e.g., satellite i) is determined, for example by SNR measurement function  72  in measurement engine  40  ( FIG. 4   a ). Conventional measurements and computations used in modern GNSS receivers to determine metrics of the signal-to-noise ratio (SNR) of the received signal, or of the carrier-to-noise ratio (C/No), are suitable for process  54 . Decision  56  selects one or more of a set of available measurement noise variance processes  60   a ,  60   b , based on the results from process  54  for satellite i. For example, process  60   b  may be suitable for satellite signals with a high SNR, while process  60   a  may be more suitable for satellite signals with a low SNR. Alternatively, both of processes  60   a ,  60   b  may be executed, with the highest variance value selected to provide a conservative result. It is contemplated that those skilled in the art having reference to this specification and its description of processes  60  herein, will be readily able to derive criteria for selecting one or more measurement noise variance processes  60  for use in connection with a received signal from satellite i, or for a particular receiver or application. 
     Measurement noise variance process  60   a  is a relatively simple approach to estimating measurement noise variance σ i   2  for the received signal from satellite i. As mentioned above, process  60   a  may be well-suited for signals with low SNR or carrier-to-noise ratio; alternatively, process  60   a  may be used for all signals. According to process  60   a , measurement noise variance σ i   2  is based directly on the SNR or carrier-to-noise ratio (C/No) of the received signal, for example as determined by SNR measurement function  72  within digital measurement engine channel  65   i  for the signal from satellite i. According to one approach, in which the bias in the pseudorange estimate is assumed to be zero, measurement noise variance σ i   2  for satellite i is determined by variance calculation function  90  in process  60   a , from the signal-to-noise ratio determined in process  54 :
 
σ i   2 =(α·10 −SNR     i/20   ) 2  
 
or, if carrier-to-noise ratio (C/No) is used:
 
σ i   2 =(α·10 −(C/No)     i/20   ) 2  
 
where α is a tunable constant set to normalize the variances or based on another criterion for the application. If the zero bias assumption is not valid, an estimate b i  of the DC bias in the measurement noise can be included in the measurement noise variance σ i   2 . For pseudorange measurement noise, for example, this bias can be estimated based on the pseudorange residual from a predicted user position (i.e., the difference between a measured pseudorange value and an otherwise-predicted pseudorange value); another approach to estimating the bias can be based on combining or filtering signals from other satellites, or as the residual determined from a previous position determination generated by position engine  42  ( FIG. 3 ). Upon obtaining estimate b i  of the pseudorange bias, this bias estimate b i  can be included in the measurement noise variance σ i   2  by:
 
σ i,adj   2 =√{square root over ( b   i   2 +σ i   2 )}
 
or, alternatively,
 
σ i,adj   =|b   i |+σ i  
 
These measurement noise variances σ i   2  are generated in similar fashion, in process  60   a , for both the pseudorange ρ and delta range δρ measurements, and are communicated as pseudorange variance signal σ 2 (ρ l ) and delta range variance signal σ 2 (δρ), which position engine  40  will arrange into entries of covariance matrices Σ ν,ρ , Σ ν,δρ  corresponding to satellite i. In this approach, the discriminator error signals D i , D i,f  shown in  FIG. 4   a  are not provided to variance computation function  90  by delay-locked loop  70  and frequency-locked loop  71 , but rather are provided by SNR measurement function  72  for the associated satellite.
 
     Process  60   b  estimates the measurement noise variance σ i   2  for the received signal from satellite i utilizing the “tracking loop” function provided by delay-locked loop  70 , and frequency-locked loop  71   i , which operate as described above relative to  FIGS. 4   a  and  4   b . In this discriminator-based measurement noise variance process  60   b , discriminator error signals D i , D i,f  for each satellite are used by variance computation function  90  to generate pseudorange variance signal σ 2 (ρ l ) and delta range variance signal σ 2 (δρ) for each satellite. Direct use of the discriminator error signals D i , D i,f  for determining the variance signals σ 2 (ρ l ) and σ 2 (δρ) is best suited for non-multipath signal receipt; in addition, the mean value of discriminator error signals D i , D i,f  is zero in steady-state if the tracking loop function is operating properly. In this case, the variances of discriminator error signals D i , D i,f  are used in producing the measurement noise variances σ 2 (ρ l ) and σ 2 (δρ), respectively, as will be described in further detail below. 
     In the multipath situation, assuming proper operation of the tracking loops (e.g., delay-locked loop  70  for pseudorange measurement and frequency-locked loop  71  for delta range measurement), the correlation peak C(τ) of the composite signal will be continually tracked. But because this peak will have significant movement over time, as the multipath parameters B, m, and θ of the reflected signal change over time, the discriminator output values D i , D i,f  will exhibit a high variance as these multipath parameters change from time interval to time interval. For example, variation of only the parameter of phase shift θ between the plots of  FIGS. 5   d  and  5   e  causes the correlation peak amplitude to range from 0.73 to 1.45 (about 6 dB), but causes the location of the correlation peak to range from −59.3 meters to +57.7 meters, assuming a constant true pseudorange. Use of the variance of the discriminator output to compute measurement noise will thus result in high values of measurement noise, and thus indicate high uncertainty in the pseudorange and delta range measurements, even if the actual measurements are in fact not that uncertain. 
     In this multipath situation, therefore, it is useful to also incorporate the multipath parameters B, m, and θ into the variance determination of process  60   b , by determining the variation in correlation peak amplitude. One can use either the SNR or the discriminator error signal to detect the multipath parameter variations, whichever computation corresponds best to the rate of variation of the multipath parameter. Typically, the SNR is computed by SNR measurement function  72  over a certain signal duration, such as 200 msec or one second; the discriminator error used in its feedback control function is also averaged over a certain signal duration. It is useful to select either the SNR or the discriminator error as the appropriate parameter for estimating the variance of the change in peak amplitude, based on which value better gives visibility into the multipath variations. 
     In this multipath case, it is straightforward to determine a difference output signal (also referred to as D i  and D i,f  for pseudorange and delta range, respectively) as the difference between the current estimated value of the peak amplitude P i  for a time interval and a mean peak amplitude value  P . The estimated value of peak amplitude P i  can be determined in various ways, including simply using the most recent output P of prompt correlator  74 P (or its equivalent in frequency-locked loop  71 ), and also including the output of SNR measurement function  72 , either in the form of SNR or carrier-to-noise ratio over a set time duration. The variances of the peak amplitude from the mean in both the time and frequency domains will be used to generate the corresponding measurement noise variances σ 2 (ρ l ) and σ 2 (δρ), as will be described below. 
     In any event, variance computation function  90  completes process  60   b  by calculating variance signals σ 2 (ρ) and σ 2 (δρ) by one of a number of available approaches, which are similar in each of the cases of the discriminator output or peak amplitude cases. In general, variance computation function  90  determines a sample variance S i  for each satellite i, and then scales that sample variance S i  into the appropriate units to arrive at the variance signals σ 2 (ρ l ) and σ 2 (δρ). 
     In one approach, a sample variance S is statistically calculated from a number N F  of samples of the error values over time. For example, sample variance S i  of discriminator output values (or correlation peak differences) D i  can be calculated as: 
               S   i     =       1     N   F       ⁢         ∑     i   =   1         N   F       ⁢       (       D   i     G     )     2               
One could also calculate the sample variance S i  based on the correlation peak amplitude (i.e., the value of P(τ) for each of N F  samples, as in this case gain G is unity). Another approach to estimating sample variance S i  is by way of an infinite impulse response (IIR) digital filter realized within variance computation function  90 . For example, a first order IIR filter can essentially compute a moving average over a time period defined by the reciprocal of the filter bandwidth:
 
                 S   ~     i     =         (              D   i     G          -       S   ~       i   -   1         )     ·   K     +       S   ~       i   -   1               
where the coefficient K is selected to provide the desired filter bandwidth, and thus the desired duration of the moving average. This approximation of the variance may require a scale factor for accuracy; for example, a scale factor of 1.25 applied to the IIR sample variance {tilde over (S)} i  is useful for Gaussian error distribution. Similar sample variance values can be determined for the delta range measurements, based on the discriminator outputs (or correlation peak differences) D i,f  from frequency-locked loop  71 .
 
     Following this determination of sample variation S i  according to the desired approach, conversion into the appropriate units (meters, for pseudorange variance signal σ 2 (ρ); meters/second, for delta range variance σ 2 (δρ)) completes process  60   b  as performed by variance computation function  90 . In the pseudorange domain, because discriminator output D i  is typically in time units of “chips”, pseudorange variance signal σ 2 (ρ) is based on measurement variance σ i   2  as follows: 
               σ   i   2     =     S   ·       (     c   1023000     )     2             
where c is the speed of light (m/s) and for the case in which the pseudo-noise spreading code has 1023000 chips per second. For generating delta range variance signal σ 2 (δρ) based on measurement variance σ i   2  in the delta range domain, the discriminator output of frequency-locked loop  71  is in units of Hz, such that:
 
               σ   i   2     =     S   ·       (     c     L   1       )     2             
where L 1  is the carrier frequency (Hz). Variance determination process  60   b  according to this embodiment of the invention, in which variance in discriminator output (correlation error) is used as a measure of measurement noise variance, is then complete.
 
     Referring back to  FIG. 6 , variance computation function  90  then finalizes the current values of the variance signals σ 2 (ρ) and σ 2 (δρ) in process  62 , if both approaches of processes  60   a ,  60   b  were executed. For example, process  62  can select the higher-valued result from the two processes  60   a ,  60   b  to select the more conservative value. Other approaches (average, weighted average, etc.) to finalizing these variance signals σ 2 (ρ) and σ 2 (δρ) in process  62  are also contemplated. 
     In process  64 , measurement engine  40  communicates the values of pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ), for each of the k satellites, to position engine  42 . These values are communicated in a manner consistent with the architecture of navigation system  20 , for example by way of corresponding signal groups  46  as shown in  FIG. 3 . The measurement and variance process of  FIG. 6  is then repeated, for the received satellite signals over a next time interval t, beginning again from process  50  as shown in  FIG. 6 . In this manner, the pseudorange and delta range values indicated by the received satellite signals, and their respective uncertainties for each satellite individually, are determined on a real-time basis, for example to reflect the movement of navigation system  20  within the coordinate system. 
     According to embodiments of this invention which base the determination of pseudorange and delta range measurement noise variances on the discriminator output signal in the tracking loop, position engine  42  can determine the current position or velocity, or both, of navigation system  20 , and corresponding uncertainty measurements for position and velocity, according to any one of a number of approaches, including both conventional techniques and also in a manner corresponding to embodiments of this invention. As such, the following description of the operation of position engine  42  will summarize conventional least-squares approaches, as well as a time-filtered approach according to embodiments of this invention. 
     As will now be described in detail relative to  FIGS. 7   a  and  7   b , position engine  42  uses the values of pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) to determine the current position and velocity of navigation system  20 , as well as an overall uncertainty measurement for those position and velocity estimates. 
     Position engine  42  is contemplated to be realized as customized hardware or logic, or as programmable logic executing computer program instructions stored in program memory  25   p  ( FIG. 2 ) or elsewhere, or some combination thereof, for carrying out the computations and processing of the contents of signal groups  46  from measurement engine  40 . More specifically, given the functions of position engine  42  described below, it is contemplated that position engine  42  will generally be realized by general purpose microprocessor logic or as application-specific processor (e.g., a digital signal processor core or device) executing program instructions stored in program memory  25   p , or in program memory more locally deployed to position engine  42  within GPS receiver  33 . It is therefore contemplated that those skilled in the art having reference to this specification will be readily able to realize position engine  42  from the functional and process-based description of its operation in connection with embodiments of this invention, which will now be provided in connection with  FIGS. 7   a  and  7   b.    
       FIG. 7   a  illustrates the operation of position engine  42  to determine the current position and velocity of navigation system  20 , and uncertainties in those position and velocity determinations, by way of a weighted or unweighted least-squares solution process. In process  100 , position engine  42  receives, from measurement engine  40 , values of pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) for each of the k received satellite signals, for example by way of signal groups  46  ( FIG. 3 ). These signal values typically correspond to the measurements based on satellite signals received over a particular time interval t, for example ranging from 200 msec to one second for consumer level navigation systems  20 . 
     As shown in  FIG. 7   a , process  101  is executed by position engine  42  to determine the current position and velocity of navigation system  20 , from the pseudorange values ρ and the delta range values δρ from each of the k satellites received in process  100  for the current time interval. Process  101  can be performed in the conventional manner for least-squares navigation, with or without a weighting matrix, as described above in connection with the background of the invention and as known in the art. If a weighted least-squares (“WLS”) approach is used, the weighting matrix W is typically constructed from the measurement noise variances σ 2 (ρ) and σ 2 (δρ) generated by measurement engine  40 . For example, this weighting matrix W may be constructed as the inverse of measurement noise covariance matrix Σ ν :
 
 W=Σ   ν   −1  
 
Typically, least-squares solution process  101  involves the determination of a pseudoinverse of a linearized derivative matrix H in which the entries correspond to linearized derivatives of the x, y, z, t coordinate values represented by each pseudorange (or delta range) value, at a nearby predicted position or velocity based on previous measurements or on results from other navigation techniques (e.g., inertial, cellular signal triangulation, etc.). The result of the least-squares solution process  101  is then converted to a convenient form, for example a latitude (φ) and longitude (λ) position on the surface of the earth, for terrestrial navigation or navigation at sea. Those skilled in the art having reference to this specification will comprehend the manner in which least-squares solution process  101  is carried out by position engine  42 , and can readily implement such a process by way of conventional program instructions for the particular architecture of position engine  42 , without undue experimentation.
 
     In parallel with process  101 , position engine  42  also determines uncertainty values for the position and velocity determinations made in process  101 , according to embodiments of this invention as shown in  FIG. 7   a . Determination of the position uncertainty in connection with this least-squares approach begins with process  102   p , in which the individual pseudorange measurement noise variance values σ 2 (ρ) from each of the k satellites are arranged into a measurement noise covariance matrix Σ ν,p  of the form:
 
Σ ν,p =diag(σ 1   2 (ρ),σ 2   2 (ρ), . . . , σ k   2 (ρ))
 
where σ i   2 (ρ) is the measurement noise variance determined by measurement engine  40  for the signal from satellite i according to the processes described above in connection with  FIG. 6 , according to embodiments of the invention.
 
     Upon arranging covariance matrix Σ ν,p  for the position uncertainty in process  102   p , position engine  42  then executes process  104   p  to determine the position-domain uncertainty value. According to this embodiment of the invention, process  104   p  solves the matrix equation:
 
Σ e,p   =H   + Σ ν,p ( H   + ) T  
 
if process  101  for determining the current position estimate is carried out according to a conventional least-squares approach. The solution of this matrix equation is position-domain covariance matrix Σ e,p , which as known in the art and as described above is a statistical representation of the position-domain error vector e caused by the measurement noise vector v. For a weighted least-squares approach used in process  101 , in which inverse matrix Σ ν,p   −1  expresses the weight matrix (i.e., such that the position solution is weighted according to the measurement noise variances determined by measurement engine  40  according to this embodiment of the invention), the solution for the position-domain error covariance matrix Σ e,p  is carried out in process  104   p  by way of:
 
Σ e,p =( H   T Σ ν,p   −1   H ) −1  
 
It is contemplated that the matrix operations of process  104   p  according to this embodiment of the invention can be readily programmed or otherwise implemented within position engine  42  by those skilled in the art having reference to this specification, without undue experimentation.
 
     According to this embodiment of the invention, the position-domain error covariance matrix Σ e,p  may be determined by process  104   p  in the earth-centered, earth-fixed (ECEF) coordinate system, considering that satellite position is typically expressed according to this system. For most navigation purposes, it is more convenient to determine the uncertainty in the position measurement in the form of uncertainty values in the form of uncertainties along the north-south axis, the east-west axis, and the vertical axis (up-down), as well as an uncertainty value regarding the time offset. Accordingly, in this embodiment of the invention, process  106   p  is performed to convert the position-domain error covariance matrix Σ e,p  from its ECEF frame of reference into the NED (north-east-down) frame of reference. 
     One approach to conversion process  106   p  follows the matrix equation: 
               Σ     e   ,   p   ,   NED       =       [           C   e   n         0           0       1         ]     ⁢       Σ     e   ,   p       ⁡     [             (     C   e   n     )     T         0           0       1         ]               
where C e   n  is a 3×3 coordinate transformation matrix for transforming ECEF to the desired frame of reference. For the case of transforming ECEF to NED, as in this example, matrix C e   n  is expressed as:
 
               C     e   ,     3   ×   3       n     =     [             -     sin   ⁡     (   ϕ   )         ⁢     cos   ⁡     (   λ   )                 -     sin   ⁡     (   ϕ   )         ⁢     sin   ⁡     (   λ   )               cos   ⁡     (   ϕ   )                 -     sin   ⁡     (   λ   )               cos   ⁡     (   λ   )           0               -     cos   ⁡     (   ϕ   )         ⁢     cos   ⁡     (   λ   )                 -     cos   ⁡     (   ϕ   )         ⁢     sin   ⁡     (   λ   )               -     sin   ⁡     (   ϕ   )               ]           
where φ, λ represent the latitude and longitude, respectively, of the position determined in process  101 . The result Σ e,p,NED  of this matrix evaluation is then converted by position engine  42  in process  106   p  into the final position-domain uncertainty metrics in the desired frame of reference. For example, if the position uncertainty metrics are to be determined in the NED frame of reference, converted covariance matrix Σ e,p,NED  can be transformed in process  106   p  according to the following:
 
σ N =β√{square root over (Σ e,p,NED (1,1))}
 
σ E =β√{square root over (Σ e,p,NED (2,2))}
 
σ D =β√{square root over (Σ e,p,NED (3,3))}
 
σ T =β√{square root over (Σ e,p,NED (4,4))}/ c  
 
where the values σ N , σ E , σ D , σ T , are the position-domain uncertainty values in the north, east, down (vertical), and time dimensions, respectively. The scaling factor β is set by the programmer for a realistic value, and the value Σ e,p,NED (i, i) refers to the (i, i) th  entry of matrix Σ e,p,NED . In addition, process  106   p  can produce position-domain uncertainty values in other frames of reference using the values σ N , σ E , σ D , σ T . For example, process  106   p  can generate a “geometric” uncertainty value σ G  as:
 
σ G =√{square root over (σ N   2 +σ E   2 +σ D   2 +σ T   2 )}
 
Process  106   p  can additionally or alternatively generate a “total position” uncertainty value σ P  as:
 
σ P =√{square root over (σ N   2 +σ E   2 +σ D   2 )}
 
and additionally or alternatively generate a “horizontal” uncertainty value σ H  as:
 
σ H =√{square root over (σ N   2 +σ E   2 )}
 
Additional measures or metrics of uncertainty can also or alternatively be generated in process  106   p  by position engine  42 .
 
     Similar processes are executed by position engine  42  to generate uncertainty values based on the measurement noise variances σ 2 (δρ) generated by measurement engine  40  relative to the delta range measurements for each received satellite signal. In this embodiment of the invention, processes  102   v ,  104   v ,  106   v  are executed by position engine  42  to operate on those delta range measurement noise variances σ 2 (δρ) for the velocity domain in much the same manner as processes  102   p ,  104   p ,  106   p  operate in the position domain. As a result, processes  102   v ,  104   v ,  106   v  produce uncertainty values for the velocity of navigation system  20  determined in process  101 , in the desired frame of reference (NED, etc.). 
     Following completion of processes  101 ,  106   p ,  106   v , in either parallel or sequential fashion, position engine  42  has produced a position, a velocity, and uncertainty values for the position and velocity values, based on the GNSS satellite signals received over the current time interval t and using a least-squares (weighted according to measurement noise, if desired) approach. In the implementation of  FIG. 3 , these values are then forwarded to blending EKF function  45 , or are alternatively used by navigation system  20  in the conventional manner. And these processes are then repeated for signals received over the next time interval, as shown in  FIG. 7   a.    
     Kalman-Filter-Based Determination of Position, Velocity, and Position and Velocity Uncertainties 
     As described above in connection with the Background of the Invention, many modern GPS receivers use a Kalman filter to solve the systems of equations based on pseudorange and delta range measurements to determine position, velocity, and time offset (“PVT”) of the navigation system. An example of the operation of the Kalman filter is described in Levy, “The Kalman Filter: Navigations Integration Workhorse”,  GPS World , Vol. 8, No. 9 (September, 1997) pp. 65-71, incorporated herein by reference. As described in that article and as known in the art, the Kalman filter uses prior estimates of position and velocity, and current measurements of pseudorange, delta range, and measurement noise to generate an updated estimate of position and velocity, and uncertainties in those measurements. 
     According to these embodiments of this invention, the pseudorange, delta range, pseudorange variances, and delta range variances may be determined by measurement engine  40  either as based on SNR or C/No ratios, or by way of using the tracking loop output as described above in connection with embodiments of this invention. Position engine  42  receives these measurements and measurement noise variances, on a satellite-by-satellite basis, from measurement engine  40 . According to these embodiments of this invention, those measurements are applied to Kalman filter methods of solving the navigation equations for PVT, as will now be described in connection with  FIG. 7   b.    
     Kalman filter process  115  is illustrated in  FIG. 7   b , in the form of an iterative process in which particular operations are repeated for each time interval of interest, using the results from the previous time interval and also measurements and measurement noise variances acquired by measurement engine  40 . Kalman filter process  115  for a current time interval begins with process  110 , in which an a priori state vector estimate is predicted based on an a posteriori state vector estimate from the previous time interval. In embodiments of this invention, for the example of an 8-state EKF implementation, the state vector x corresponds to:
 
 x=[x,y,z,−ct   u   ,{dot over (x)},{dot over (y)},ż,−c{dot over (t)}   u ] T  
 
in which [x, y, z] and [{dot over (x)}, {dot over (y)}, ż] represent the current position and velocity, respectively, of navigation system  20 , and where t u  and {dot over (t)} u  represent the clock bias and drift, respectively, of navigation system  20  (c being the speed of light). In process  110 , as mentioned above, an a posteriori state vector estimate {circumflex over (x)} k-1   +  is applied to a system equation:
 
 x   k   =Ax   k-1   +w   k  
 
to derive the a priori state vector estimate {circumflex over (x)} k   −  for the current time interval k iteration. In this system equation, vector w k  represents a model of the process noise, as known in the art, and system matrix A is given by:
 
             A   =     [         1       0       0       0       T       0       0       0           0       1       0       0       0       T       0       0           0       0       1       0       0       0       T       0           0       0       0       1       0       0       0       T           0       0       0       0       1       0       0       0           0       0       0       0       0       1       0       0           0       0       0       0       0       0       1       0           0       0       0       0       0       0       0       1         ]           
where T is the sample time period (i.e., time difference between the successive values of the state vectors x k-1  and x k ). In this embodiment of the invention, the matrix operation:
 
{circumflex over ( x )} k   −   =A{circumflex over (x)}   k-1   + 
 
is carried out in process  110 .
 
     The a priori state vector estimate {circumflex over (x)} k   −  from process  110  is then applied by position engine  42  to process  112 , in which an a priori covariance matrix P k   −  is evaluated from the a posteriori covariance estimate P k-1   +  from the previous time interval, in combination with process noise vector w k , which is assumed to be a Gaussian random vector with zero mean and covariance Q k . Process noise vector w k , and thus its covariance Q k , are assumed to remain constant over time, as this noise refers to the system noise and precision of navigation system  20  itself. More specifically, process  112  is executed by position engine  42  to evaluate the matrix equation:
 
 P   k   −   =AP   k-1   +   A   T   +Q   k  
 
     The measurement equation:
 
 z   k   =h   k ( x   k )+ v   k  
 
defines the relationship among measurement vector z k , which includes up to two measurements (either or both of pseudorange ρ and delta range δρ) for each satellite, the true state vector x k , and measurement noise vector v k , which is assumed to have zero mean and covariance R k . According to embodiments of this invention involving execution according to a Kalman filter approach, this non-linear measurement equation is linearized about the current value of the state equation:
 
 z   k   =H   k   x   k   +v   k  
 
An example of this linearization of the measurement equation, illustrating the arrangement of the specific entries and values for the case in which signals are received from four satellites, can be shown by:
 
               [           ρ   1               ρ   2               ρ   3               ρ   4                   ρ   .     1     ⁢                         ρ   .     2                 ρ   .     3                 ρ   .     4           ]     =         [           a     x   ⁢           ⁢   1             a     y   ⁢           ⁢   1             a     z   ⁢           ⁢   1           1       0       0       0       0             a     x   ⁢           ⁢   2             a     y   ⁢           ⁢   2             a     z   ⁢           ⁢   2           1       0       0       0       0             a     x   ⁢           ⁢   3             a     y   ⁢           ⁢   3             a     z   ⁢           ⁢   3           1       0       0       0       0             a     x   ⁢           ⁢   4             a     y   ⁢           ⁢   4             a     z   ⁢           ⁢   4           1       0       0       0       0           0       0       0       0         a     x   ⁢           ⁢   1             a     y   ⁢           ⁢   1             a     z   ⁢           ⁢   1           1           0       0       0       0         a     x   ⁢           ⁢   2             a     y   ⁢           ⁢   2             a     z   ⁢           ⁢   2           1           0       0       0       0         a     x   ⁢           ⁢   3             a     y   ⁢           ⁢   3             a     z   ⁢           ⁢   3           1           0       0       0       0         a     x   ⁢           ⁢   4             a     y   ⁢           ⁢   4             a     z   ⁢           ⁢   4           1         ]     ⁡     [         x           y           z             -     ct   u                 x   .               y   .               z   .                 -   c     ⁢       t   .     u             ]       +     [           v   x               v   y               v   z               v     ct   u                 v     x   .                 v     y   .                 v     z   .                 v     c   ⁢       t   .     u               ]             
where ρ i  and {dot over (ρ)} i  are the pseudorange and the delta range measurements, respectively, for the satellite i, and a xi , a yi , a zi  are the x, y, and z components, respectively, of the unit norm vector pointing from the position of navigation system  20  to that satellite i. It is this linearized measurement equation that is evaluated by position engine  42  by its executing of processes  114 ,  116 , as will now be described.
 
     In practice, any number of measurements from any number of satellites may be relevant to the current navigation process. In some cases, no measurements may be available over a time interval; in those cases, the predicted state vector from the corresponding time interval will simply be the propagation of the result from the previous time interval. In other cases, some satellites may not provide new measurements in that time interval, or measurement engine  40  may discard some measurements, for example if the signal-to-noise ratio for a particular satellite&#39;s measurements is too low. It is contemplated that the process of  FIG. 7   b  according to embodiments of this invention will be capable of managing the updated state vector prediction in such cases, again by perhaps propagating the previous measurements and variances from the previous interval. 
     In process  114 , position engine  42  uses the a priori covariance matrix P k   −  from process  112 , along with measurement noise variances from measurement engine  40 , to arrive at a gain matrix K k . As shown in  FIG. 7   b , the current variances are obtained by position engine  42  in process  100 , by receiving values of pseudorange ρ, delta range δρ, pseudorange variance σ 2 (ρ), and delta range variance σ 2 (δρ) from measurement engine  40 , for each of the k received satellite signals over a particular time interval t, for example by way of signal groups  46  ( FIG. 3 ). In process  102 ′, similarly as described above in processes  102   p ,  102   v  shown in  FIG. 7   a , the individual pseudorange measurement noise variance values σ 2 (ρ) and delta range measurement noise variance values σ 2 (δρ) are arranged into a measurement noise covariance matrix, specifically covariance matrix R k  for use in the measurement equation evaluation. The most recent position estimate defines measurement matrix H k , which as described above is expressed in terms of direction cosines of the unit vector pointing from that position estimate to the satellite. As a result, process  114  evaluates gain matrix K k  for this iteration as:
 
 K   k   =P   k   −   H   k   T ( H   k   P   k   −   H   k   T   +R   k ) −1  
 
And once this gain matrix K k  is evaluated for this iteration k in process  114 , then position engine  42  executes process  116  to produce an updated (i.e., a posteriori) state vector estimate {circumflex over (x)} k   +  from:
 
 {circumflex over (x)}   k   +   ={circumflex over (x)}   k   −   +K   k   [z   k   −h   k ( {circumflex over (x)}   k   − )]
 
where vector z k  is the vector of current measurement values of pseudorange ρ and delta range δρ, as acquired in process  100 . The values of this state vector estimate {circumflex over (x)} k   +  generated in process  116 , which as described above provides the current estimate of the position and velocity of navigation system  20 , can be output by position engine  42  to blending EKF  45  ( FIG. 3 ), or output elsewhere for purposes of the navigation being carried out by navigation system  20 . Conversion of these values into the particular frame of reference (e.g., NED rather than ECEF), can be carried out in the manner described above relative to  FIG. 7   a.  
 
     In process  118 , position engine generates an updated, or a posteriori, estimate of the covariance matrix P k   + , according to embodiments of this invention. In a general sense, covariance matrix P k   +  for the error in the a posteriori state vector estimate {circumflex over (x)} k   +  is calculated in this Kalman filter process  115 , and is defined as the matrix of expected value of the squares of the difference between the current estimate {circumflex over (x)} k   +  and the true position x k :
 
 P   k   +   =E [( x   k   −{circumflex over (x)}   k   + )( x   k   −{circumflex over (x)}   k   + ) T ]
 
According to embodiments of this invention, this updated, or a posteriori, estimate of the covariance matrix P k   +  is based on the a priori estimate and gain matrix K k  by position engine  42  executing process  118  to evaluate the matrix equation:
 
 P   k   + =( I−K   k   H   k ) P   k   − 
 
This a posteriori estimate of the covariance matrix P k   +  is thus essentially the filtered output of the process noise covariance Q k  (as incorporated into the a priori estimate of covariance matrix P k   −  in process  114 ) and the measurement noise covariance R k  applied to the gain calculation of process  116 ). A posteriori covariance matrix estimate P k   +  is then used in the prediction of the a priori covariance matrix P k   −  in the next iteration of Kalman filter process  115  for the next time interval (i.e., after the incrementing of index k performed in process  119 ).
 
     The uncertainties in position, velocity, time, and time drift expressed in a posteriori covariance matrix estimate P k   +  for the current time interval are output by position engine  42  via process  120 . According to this embodiment of the invention, submatrices of a posteriori covariance matrix estimate P k   +  constitute the position-domain and velocity-domain uncertainties. More specifically, in this embodiment of the invention, position uncertainty matrix P k,pos   +  is expressed in the submatrix:
 
 P   k,pos   +   =P   k   + (1:4,1:4)
 
which in other words amounts to the submatrix of the entries in the first four rows and first four columns of a posteriori covariance matrix estimate P k   + . Similarly, velocity uncertainty matrix P k,vel   +  is expressed in the submatrix:
 
 P   k,vel   +   =P   k   + (5:8,5:8)
 
namely, the submatrix of the entries in the fifth through eighth rows and fifth through eighth columns of a posteriori covariance matrix estimate P k   + . The values along the diagonals of velocity uncertainty matrices P k,pos   + , P k,vel   +  return the uncertainty values for each of the x, y, z, and t dimensions in the ECEF frame of reference. Process  120  also includes the necessary operations, carried out by position engine  42 , to convert these uncertainties from the ECEF frame of reference to one or more other desired frames of reference, such as NED, as described above relative to processes  106   p ,  106   v  of  FIG. 7   a.  
 
     It has been observed, in connection with this invention, that in some instances the process noise variances or even the measurement noise variances used by position engine  42  to derive the position and velocity uncertainties can be somewhat unrealistic, which of course results in unrealistic position and velocity uncertainties. One such instance occurs in the case in which low quality measurements are provided by measurement engine  40 , requiring the Kalman gain matrix K k  to be reduced, applying a heavy smoothing to the iterative calculations. 
     According to an alternative embodiment of the invention, a separate Kalman filter “core” is implemented within position engine  42 , to generate separate estimates of position and velocity uncertainties using the time-varying measurement noise covariance R k  from measurement engine  40  in combination with a dynamic process noise estimate Q k * in which “realistic” (i.e., time-varying) dynamics of navigation system  20 , including those based on user inputs, are reflected. This separate Kalman filter “core” can be realized in hardware by parallel processor hardware within position engine  42 , or alternatively by a parallel or pipelined computing process relative to Kalman filter process  115 . It is contemplated that those skilled in the art having reference to this specification will be readily able to realize such a separate Kalman filter core in a manner well-suited for particular applications. 
       FIG. 7   c  illustrates the operation of a parallel approach according to this alternative embodiment of the invention. In this arrangement, separate Kalman filter process  125  is shown as performed in parallel with Kalman filter process  115  described above in connection with  FIG. 7   b . Kalman filter process  125  includes processes  132 ,  134 ,  138 , which correspond to respective processes  112 ,  114 ,  118  of Kalman filter process  115  ( FIG. 7   b ), and which are similarly iteratively executed by process engine  42  for each time interval, indexed by index k as in Kalman filter process  115 . According to this embodiment of the invention, however, process  135  is also performed by process engine  42  to dynamically estimate changes in the process noise of navigation system  20 , including any actual user dynamics, such as those corresponding to the mode of travel (i.e., motion by vehicle, pedestrian movement, bicycle, etc.). Process  135  may simply select process noise covariance matrix Q k * from among a set of predetermined covariance matrices stored in data memory  25   d  or elsewhere, each covariance matrix associated with one of the navigational modes (vehicle, pedestrian, bicycle, etc.); the selected predetermined covariance matrix would then be selected in process  135 . Alternatively, position engine  42  may perform process  135  by deriving process noise covariance matrix Q k * in n a manner that accounts for the actual user dynamics. For example, signals from inertial measurement unit  30  ( FIG. 2 ) may be used by position engine  42  to derive process noise covariance matrix Q k *, in this case based on inertial measurements of position and velocity. In any case, the values of that process noise covariance matrix Q k * determined by position engine  42  in process  135  are forwarded to process  132  for the corresponding time interval. 
     In process  132 , as in the case of process  112  described above, a priori estimate of covariance matrix P k   −  is evaluated based on the a posteriori covariance estimate P k-1   +  from the previous time interval. But in process  132 , this a priori covariance matrix P k   −  estimate is also based on process noise covariance Q k * from process  135 , which as described above includes dynamic effects in the process noise vector w k . Accordingly, position engine executes process  132  to evaluate the matrix equation:
 
 P   k   −   =AP   k-1   +   A   T   +Q   k *
 
Following process  134 , processes  134 ,  138  of Kalman filter process  125  are then executed by position engine  42  in similar manner as described above in connection with processes  114 ,  118  of Kalman filter process  115 . The result of process  138 , specifically, is an a posteriori estimate of covariance matrix P k   +*  that includes the effects of dynamic process noise changes as determined by process  135 . This filtered a posteriori covariance matrix estimate P k   +*  is forwarded to process  140 , along with a posteriori covariance matrix estimate P k   +  from process  118  in Kalman filter process  125 , from which the uncertainties in position and velocity are determined and converted to the desired frame of reference, in similar manner as described above in connection with process  120 . Process  140  may select from between the parallel a posteriori covariance matrix estimates P k   + , P k   +*  in determining the final uncertainty values, or may alternatively combine the two estimates in some sort of weighted fashion according to a separate parameter in so doing. It is contemplated that the actual position and velocity output may be based solely on the result of Kalman filter process  115  as shown in  FIG. 7   c , or alternatively may include some or entirely the contribution of Kalman filter process  125 . Those skilled in the art having reference to this specification will readily derive the particular manner in which the parallel results can be best used for particular situations.
 
     According to any of these embodiments of the invention, position engine  42  provides output signals indicative of the current position and velocity of navigation system  20  in its environment, along with measures of uncertainty for each of those position and velocity calculations. In the example of  FIG. 3 , the position and velocity uncertainty measurements are used by blending EKF function  45  in its overall operation of generating a final position and velocity signal on lines POS_INT for communication by navigation system  20  to the user, or for other purposes within navigation system  20 . It is contemplated that blending EKF function  45  will produce the integrated position and velocity signal from a combination of the GPS and inertial position and velocity results, in the extreme case by selecting one or the other of the GPS and inertial results, based on the uncertainties in the GPS position and velocity as determined by position engine  42 . For example, if the uncertainties in either or both position and velocity as generated by GPS navigation function  35  according to an embodiment of this invention are relatively large in value, blending EKF  45  will favor the position and velocity determinations generated by inertial navigation system  44  in arriving at the final integrated position and velocity result. Alternatively, these uncertainty values generated by GPS navigation function  35  may be used for other purposes, for example including the adjustment of operational parameters within navigation system  20 , triggering warnings or alarms resulting from high uncertainty, or other operations and functions, which will be apparent to those skilled in the art having reference to this specification. 
     According to embodiments of this invention, therefore, a GNSS navigation system can be constructed and can operate to provide realistic and accurate uncertainty measurements for its position and velocity determinations. These embodiments of the invention provide such realism because, among other reasons, the measurement noise is considered individually and realistically for each satellite signal, and in a manner that contemplates multipath reflected signal propagation. In addition, these improved individual measurement noise estimates are determined in a manner that is compatible with modern Kalman filter-based navigation algorithms, thus obtaining the benefit of that time-dependent approach to navigation. As a result, the resulting uncertainty determination neither underestimates nor overestimates the current signal conditions and environment of the navigation system, improving the ability of the user and the system itself to provide accurate and useful information. 
     While the present invention has been described according to its embodiments, it is of course contemplated that modifications of, and alternatives to, these embodiments, such modifications and alternatives obtaining the advantages and benefits of this invention, will be apparent to those of ordinary skill in the art having reference to this specification and its drawings. It is contemplated that such modifications and alternatives are within the scope of this invention as claimed.