# Patent

## Title
Method for refining state estimates of a global positioning system (GPS) receiver based upon sequential GPS carrier phase measurements

## Abstract

Sequential GPS Doppler carrier phase count measurements are used for precision sequential determination of position and velocity of a GPS receiver, such as for orbit determination and geolocation, with minimum throughput time. Real-time orbit determination and geolocation performance is enabled with an optimal sequential filter, and near-real-time performance is enabled with an optimal fixed-lag smoother. Many problems associated with prior art orbit determination are eliminated by addressing the “cycle slip” problem, the unknown initial range problem with RANGECP measurement representations, the problem of serial correlation in the measurements due to reprocessing of overlapping thermal noise. Also, the present invention significantly attenuates the carrier signal phase variation due to rotation of receiver antenna relative to transmitter antenna because the sequential phase count time intervals are sufficiently short.

## Background
CROSS-REFERENCE TO RELATED APPLICATIONS
This application is a continuation of prior application Ser. No. 11/140,296, filed May 27, 2005 now U.S. Pat. No. 7,050,002, which is hereby incorporated by reference.


BACKGROUND OF THE INVENTION
The GPS (Global Positioning System) consists, in part, of a constellation of NAVSTAR (Navigation System with Timing And Ranging) spacecraft with transmitters that continuously broadcast L-BAND radio carrier signals on two frequencies. NAVSTAR spacecraft all have near-circular orbits with radial distances that are 4.16 Earth radii from Earth center, with heights of 20182 km above the Earth's surface. With orbit period of 717.957 minutes and inclination to Earth's equator of 55 degrees, each NAVSTAR completes two orbits about the Earth in one sidereal day (a sidereal day is defined by one complete rotation of the Earth relative to the mean locations of a particular ensemble of stars, or to a particular ensemble of quasars. A solar day has 1440 minutes). NAVSTAR radio transmissions are directed towards the Earth.
USER spacecraft employ GPS receivers, and for most effective receipt of NAVSTAR radio transmissions, have orbits with radial distances that are significantly less than 4.16 Earth radii from Earth center. The large existing class of spacecraft in Low Earth Orbit (LEO) have radial distances that are in the neighborhood of 1.1 Earth radii from Earth center, with heights of roughly 600 km above the Earth's surface (a height of 637.81 km goes with radial distance of 1.1000 er). The orbits of USER spacecraft are estimated by processing measurements derived from NAVSTAR radio signals with the aid of GPS receivers. That is, GPS measurements are used for orbit determination of USER spacecraft equipped with GPS receivers.
Also, Earth-fixed GPS receivers are placed at unknown locations and are activated to collect GPS measurements. These are used to estimate the Earth-fixed receiver locations. We refer to this capability as geolocation.
There are two very different types of measurements derived from the NAVSTAR L-BAND radio carrier signals, namely: (i) pseudo-range measurements and (ii) Doppler carrier phase count measurements. Doppler carrier phase count measurements provide Doppler information on the relative velocity of USER spacecraft with respect to NAVSTAR spacecraft. This enables much more precise orbit estimates and geolocation estimates than those that process pseudo-range measurements only. The present invention refers to USER orbit determination and/or geolocation by processing Doppler carrier phase count measurements, not by processing pseudo-range measurements.
Measurement Representations
All methods of precise orbit determination and geolocation map measurement residuals linearly to state estimate corrections. A measurement residual is the difference between a measurement produced by receiver hardware and the estimated representation of the measurement derived from a state estimate—calculated on a computer. The measurement representation is critical to precise orbit determination and geolocation. In one aspect, the present invention refers most significantly to an improvement of the estimated representation of Doppler carrier phase count measurements.
Geodesy and Geophysics
Orbit determination of USER spacecraft, and geolocation of ground-fixed receivers, has been performed for several decades by processing GPS Doppler carrier phase count measurements. This has been particularly relevant for geodesy and geophysics, and very significant progress has been thereby achieved. But the measurement representations currently used in orbit determination and geolocation capabilities for geodesy and geophysics are fraught with serious problems, needlessly so. The Doppler carrier phase measurement is represented as a range measurement with unknown initial range. How can this be?
Every rigorous measurement representation equation for Doppler carrier phase is composed of several additive terms on the right, with the measurement representation on the left. Some terms on the right are characterized by signal, and some are characterized by noise. One of the terms for signal is called range-difference, or delta-range. This difference is: range at the end of the Doppler carrier phase count interval less range at the beginning of the Doppler carrier phase count interval, where range is a measure of distance between NAVSTAR transmitter and USER GPS receiver. One can rearrange the measurement representation equation so that the range-difference appears on the left, and is a function of the Doppler carrier phase measurement and other terms on the right. Algebraically one can then add the unknown range, at the beginning of the Doppler carrier phase count interval, to both sides of the equation. Now the sum of the unknown range on the left, at the beginning of the Doppler carrier phase count interval, and range-difference on the left defines the range at the end of the Doppler carrier phase count interval. That is, the modified equation presents the range at the end of the Doppler carrier phase count interval on the left as a function of the Doppler carrier phase measurement and the unknown range at the beginning of the Doppler carrier phase count interval, together with other terms, on the right. So orbit determination and geolocation for geodesy and geophysics treats the range at the end of the Doppler carrier phase count interval as the measurement, and requires an associated range measurement representation to form a measurement residual. This is how the Doppler carrier phase measurement is represented as a range measurement with unknown initial range. We shall call it the RANGECP measurement representation, with the superscript CP meaning that RANGECP is derived, in part, from a carrier phase measurement.
Further, the carrier phase measurements used for these range representations are overlapping measurements, not sequential measurements. This has a significant detrimental impact for orbit determination and geolocation, and is described below.
Presentation of Estimation Errors
The prior art geodesist presents small constant error magnitudes (i.e., those not derived directly from his state estimate error covariance matrix function) for his estimates of orbits and geolocations using overlapped RANGECP measurements, and has apparently convinced some technical people of their validity. These small error magnitudes usually refer to position components only, a limited subset of any complete state estimate. The argument has been made that if the geodesist is doing so well, then his estimation technique must surely be correct.
An extensive survey of the relevant literature has failed to identify a single analysis of the state estimate error, including a realistic error covariance. In fact, these analyses do not exist because of limitations in the prior art, including but not limited to the use of overlapped carrier phase measurements.
Each and every technique (including prior art methods of batch least squares, sequential filtering, batch filtering, and sequential smoothing) used by the geodesist to calculate precision estimates also calculates an associated state estimate error covariance matrix, or its inverse, that is used intrinsically to calculate the estimates. If these estimates are optimal and appropriately validated estimates, then their associated state estimate error covariance matrices are realistic. Realistic state estimate error covariance matrices would be invaluable for interpretation and use of orbit and geolocation estimates, but they are never presented. They are never presented because they are always significantly unrealistic. This is frequently admitted and is a problem with the prior art. Credible validation techniques in support of realistic state estimate error covariance presentations cannot be found. The present inventors have looked for them and have asked for them in appropriate places, but they apparently do not exist.
GPS Doppler Carrier Phase Count Measurements
NAVSTAR L-BAND radio carrier signals with specified frequencies are generated by atomic clocks on board each NAVSTAR. When received by any GPS receiver, the frequencies of these signals are Doppler shifted according to the deterministic theory of special and general relativity. The Doppler shift is proportional to the velocity of the USER receiver with respect to the velocity of the NAVSTAR transmitter. Each GPS receiver has its own local clock, and it is designed to create a radio signal with approximately the same frequency as that transmitted by each NAVSTAR. The incoming NAVSTAR Doppler shifted signal to each GPS receiver is differenced with the local receiver signal to isolate the Doppler shift and thereby generate a differenced Doppler phase with Doppler frequency. The Doppler phase is counted by the GPS receiver across a time interval [tm; tn] defined by the receiver clock, where tm<tn, and where m and n are time indices.
[00 sec,10 sec], [10 sec; 20 sec], [20 sec; 30 sec], [30 sec; 40 sec], . . . sequential
[00 sec,10 sec], [00 sec; 20 sec], [00 sec; 30 sec], [00 sec; 40 sec], . . . . overlapping
Define τmn=tn−tm. Two illustrative sequences of four time intervals are displayed here, where the first sequence has a constant τmn=10 seconds (the constant sequential value for τmn will be orbit dependant and relatively small as compared to orbit period) and is sequential, but the second sequence is overlapped with non-constant τmn expanding without bound with each measurement. Each overlapped time interval is the concatenation of a sequence of adjacent sequential time intervals.
FIGS. 1A and 1B have general applicability and distinguish the sequential and overlapped Doppler measurement time intervals graphically.
The phase count measurement Nmn across each time interval [tm; tn] is the sum of an integer number of cycles plus a partial cycle. Thus the GPS Doppler carrier phase count measurement Nmn is a rational number (but typically not an integer), and has units of cycles. If the phase count time intervals are sequential, then the Doppler phase count measurements are sequential, non-overlapping, and τmn remains constant and small. But if the time intervals are overlapped, then the Doppler phase count measurements are overlapped, and τmn becomes indefinitely large. Each overlapped measurement is the sum of a sequence of adjacent sequential measurements.
Table 1 (sequential measurements) and Table 2 (overlapping measurements) present related examples (not that many other distinct and useful examples could be constructed and simplifications have been made for purposes of this example; ordinarily all calculations are performed in double precision to 15+decimals and time is represented in the Gregorian Calendar form or Julian Date) of simulated GPS Doppler carrier phase count measurement values for a GPS receiver on a spacecraft in LEO (Low Earth Orbit), referred to the same NAVSTAR transmitter. These tables are coordinated with the time intervals displayed above and with FIGS. 1A-1B. Columns 1, 2, 4, and 7 are the same for each table. Column 1 is time index, Column 2 is time t, Column 4 is range ρ, and Column 7 is the receiver phase count N0nR. Columns 3, 5, and 6 are different for each table due to the distinction in measurement time interval differences τmn. Column 3 is time interval difference τmn, Column 5 is delta-range Δρ, and Column 6 is the measurement: Nmn for sequential measurement and N0n for overlapping measurement.







TABLE 1







Sequential Doppler Carrier Phase Measurements















τmn






index
t (sec)
(sec)
ρ (m)
Δρ (m)
Nmn (cy)
NOnR (cy)
















0
00

18632055.19


0.00




10

−955.17
−5019.47


1
10

18631100.02


−5019.47




10

−274.49
−1442.44


2
20

18630825.53


−6461.91




10

406.08
2133.93


3
30

18631231.61


−4327.98




10

1086.39
5709.03


4
40

18632318.00


1381.05









The time-tag for each measurement is defined here to be the time at end of the Doppler count interval. Example for Table 1: For time interval [tm; tn]=[20 seconds; 30 seconds], and sequential measurement N23=2133.93 cycles, the time-tag t3=30 seconds.







TABLE 2







Overlapped Doppler Carrier Phase Measurements















τmn






index
t (sec)
(sec)
ρ (m)
Δρ (m)
NOn (cy)
NOnR (cy)
















0
00

18632055.19


0.00




10

−955.17
−5019.47


1
10

18631100.02


−5019.47




20

−1229.66
−6461.91


2
20

18630825.53


−6461.91




30

−823.58
−4327.98


3
30

18631231.61


−4327.98




40

262.81
1381.05


4
40

18632318.00


1381.05









BRIEF SUMMARY OF THE INVENTION
The present invention provides precision sequential orbit determination with minimum throughput time using sequential GPS Doppler carrier phase count measurements. Real-time orbit determination performance is enabled with an optimal sequential filter, and near-real-time performance is enabled with an optimal fixed-lag smoother. The present invention eliminates many problems associated with prior art orbit determination by addressing the “cycle slip” problem, the unknown initial range problem with RANGECP measurement representations, the problem of serial correlation in the measurements (due to reprocessing of overlapping thermal noise), and significantly attenuates the carrier signal phase variation (due to rotation of receiver antenna relative to transmitter antenna) because the sequential phase count time intervals are sufficiently short.
Construction of Sequential Doppler Phase Count Measurement Representations
Sequential Doppler phase count measurement representations are constructed:

    1. to process sequential Doppler carrier phase count measurements according to the new representation given by Eq. 10 (natural units: cycles);
    2. to rigorously present each Doppler carrier phase count measurement as a delta-range representation according to the new expression defined by Eq. 21 (derived units: meters);
    3. to remove ionospheric effects with the two GPS L-BAND frequencies according to Eq. 22;
    4. by choosing Eq. 6 rather than Eq. 8 for modeling USER clock phenomenology;
    5. to remove USER spacecraft clock phenomenology with Doppler measurement differences referred to distinct NAVSTARs according to Eq. 30;
    6. to remove effects of USER spacecraft clock frequency steering, when it exists, with Doppler measurement differences referred to distinct NAVSTARs according to Eq. 30; and
    7. for processing in an optimal sequential filter, or optimal sequential fixed-lag smoother, with autonomous measurement editor.


Considered individually, elements (1), (2), and (4) are new and unique. Considered individually, different versions of elements (3), (5), and (6) have been used in previous orbit determination capabilities. The optimal sequential filter referred to in element (7) has been produced by Analytical Graphics Inc., of Malvern Pa. (AGI), the assignees of the present invention, and is in current use by AGI customers.
The optimal fixed-lag smoother referred to in element (7) has been designed by AGI, but has not yet been implemented. Since the new element (1) is embedded in elements (3), (5), and (6), they are new in an embedded sense. Taken together, this combination of seven elements defines a new capability to perform precision sequential orbit determination with minimum throughput time using sequential GPS Doppler carrier phase count measurements. Real-time orbit determination performance is enabled with an optimal sequential filter, and near-real-time performance is enabled with an optimal fixed-lag smoother.
As mentioned above, the “cycle slip” problem is eliminated, the unknown initial range problem with RANGECP measurement representations is eliminated, serial correlation in the measurements (due to reprocessing of overlapping thermal noise) is eliminated, and the carrier signal phase variation (due to rotation of receiver antenna relative to transmitter antenna) is significantly attenuated because the sequential phase count time intervals are sufficiently short.
Elimination of the Cycle Slip Problem
The following example illustrates the cycle slip problem for overlapping measurements. Vast experience with overlapping measurements has shown the cycle slip problem to be a very significant problem in general. The present inventors submit, generally, that the processing of overlapping measurements is significantly inferior to the processing of sequential measurements by an optimal sequential filter.
EXAMPLE
From Table 1: At time-tag t=30 seconds, the sequential phase count measurement is N23=2133.93 cycles, across time interval [20; 30] seconds.
From Table 2: At time-tag t=30 seconds, the overlapping phase count measurement is N03=−4327.98 cycles, across time interval [00; 30] seconds.
Notice that: (−4327.98 cycles)=(−6461.91 cycles+2133.93 cycles), where (−4327.98 cycles) and (−6461.91 cycles) are overlapping measurements and (2133.93 cycles) is simultaneously a sequential measurement and a difference of overlapping measurements. These numbers represent valid measurements without cycle slip in this example; i.e., without loss of lock in the PLL (Phase Lock Loop).
Suppose now, as frequently happens with real data, that there is a loss of lock for sequential measurement N23. That is, the receiver erroneously reports say N23=10000.00 cycles, instead of the correct value N23=2133.93 cycles for the sequential measurement. Then an erroneous overlapping measurement N03 is reported as N03=(−6461.91 cycles+10000.00 cycles)=(−3538.10 cycles), instead of the correct value N03=(−4327.98 cycles).
Suppose further, as frequently happens, that the sequential measurement N34=5709.03 cycles is reported correctly by the receiver. Then the overlapping incorrect measurement value N04=(−3538.10 cycles+5709.03 cycles)=(2170.93 cycles) is reported by the receiver instead of the correct value of N04=(1381.05 cycles). Sequential measurement N34 is correct, but overlapping measurement N04 is incorrect because overlapping measurement N03 was incorrect.
Suppose further, as frequently happens, that sequential measurements N45, N56, N67, . . . following N34 in this tracking pass, are reported correctly by the receiver. Then all overlapping measurements following N04 will be reported incorrectly because they are additively chained to the incorrect overlapping measurement N04. Overlapping measurements N04, N05, N06, N07, . . . are useless.
Overlapping Measurements
There are two cases to consider for the response of any estimator to the overlapping measurements. (i) If the estimator has an effective measurement editor, then overlapping measurements N04, N05, N06, N07, . . . will be appropriately discarded (note that this takes significant time in any operations) together with their associated measurement information, but (ii) If the estimator does not have an effective measurement editor, then some subset of the overlapping measurements N04, N05, N06, N07, . . . will be accepted and processed for orbit determination. This will destroy the orbit estimate in the sense that the orbit estimate will incur very significant error.
In the first case the orbit estimate suffers a significant loss of information and is thereby degraded, and in the second case, the orbit estimate is destroyed. This is due to a cycle slip (loss of lock in the PLL) in a single sequential measurement N23, and is referred to here as the cycle slip problem.
Sequential Measurements
Consider now the response of our optimal sequential filter to the sequential measurements: Erroneous measurement (N23=10000.00 cycles) is rejected by the optimal sequential autonomous editor, but all other sequential measurements N34, N45, N56, N67, . . . are accepted by the optimal sequential autonomous editor and appropriately used for sequential orbit determination. The only measurement information necessarily lost is that associated with measurement N23. We refer to this as the elimination of the cycle slip problem.
The autonomous identification and rejection of any sequential carrier phase measurement that suffers loss of phase lock is accomplished in the following manner. Suppose a good sequential carrier phase measurement yk at time tk is used to correct the propagated a priori state estimate Xk|k−1 at time tk to calculate an improved state estimate Xk|k at time tk. Suppose the next measurement yk+1 at time tk+1 had a loss of phase lock, and is thus known to have a cycle slip. However, there is no way to determine that yk+1 is a bad cycle-slipped measurement by looking solely at the carrier phase measurement value. Rather the state estimate Xk|k and its realistic error covariance matrix Pk|k is propagated to time tk+1 to get Xk+1|k and Pk+1|k, and Xk+1|k is used to calculate a measurement representation y(Xk+1|k). The measurement residual Δyk+1=yk+1−y(Xk+1|k) is formed and the absolute value |Δyk+1| is compared to a threshold Tk+1 derived from Pk+1|k and from the measurement error variance matrix Rk+1. If |Δyk+1|<Tk+1, then measurement yk+1 is accepted for processing by the optimal filter. But in the case of a cycle slip, we will find that |Δyk+1|≧Tk+1 and the measurement yk+1 is therefore autonomously rejected in real time.
Suppose the next sequential measurement yk+2 at time tk+2 is a good measurement (no cycle slip). Then we will find |Δyk+2|<Tk+2, and the measurement yk+2 is accepted for processing by the filter.
The autonomous real-time identification and rejection of cycle slips in sequential carrier phase measurements is enabled with the autonomous sequential editor, sequential estimator, and realistic covariance function.
This autonomous identification and rejection of flawed data is impossible with overlapped carrier phase measurements. A batch of overlapped carrier phase measurements are processed simultaneously with a batch estimator that prohibits autonomous editing, prohibits use of a realistic state error covariance function, and requires significant time to repair the damage done by a single sequential cycle slip.
Elimination of the Measurement Serial Correlation Problem
The processing of sequential measurements, in preference to overlapping measurements, eliminates the repeated exposure to the same thermal noise in the overlapping measurements that generates the serial correlation.
Elimination of the Unknown Initial Range Problem
Given only the GPS Doppler carrier phase count measurements Nmn and N0n and their time-tags tn, all values of range ρn are unknown (the values of GPS pseudo-range have error magnitudes significantly larger than those of Δρmn, derived from Doppler phase count). However, the range-difference Δρmn, or delta-range, can be calculated from Nmn according to Eq. 21, and Δρ0n can be calculated from N0n, also according to Eq. 21. From the definition of Δρ0n:

Δρ0n=ρn−ρ0 

That is:

ρn=ρ0+Δρ0n  (1)

for n∈{1, 2, . . . }. We shall refer to ρ0 as the unknown initial range. Given the calculation of Δρ0n from measurements N0n for nε{1, 2, . . . }, each unknown range ρn is treated by the geodesist as though it were a measurement, for n∈{1, 2, . . . }, with unknown initial range ρ0. We refer to ρn above as ρn=RANGECP.

The Unknown Initial Range Problem has three components: Communication, Estimation, and Throughput.
Communication
By definition, a measurement is a known value of an object or quantity determined by measuring it, but ρn=RANGECP is not measured and is not known. We have multiple contradictions (n∈{1, 2, . . . }) because the ρn=RANGECP are called measurements by the geodesist when they are not measurements. This is a problem in communication.
Estimation
From Eq. 1, we see that ρ0 and Δρ0n have distinct time-tags t0 and tn>t0. The unknown initial range ρ0 is estimated, together with orbit parameters and other state parameters, from measurements N0n. But the execution of a state error transition function between times t0 and tn is required for each Δρ0n to perform the estimation. The state error matrix includes errors on orbit parameters and errors on ρ0. Significant cross correlations between orbit errors and errors on ρ0 alias orbit errors and other state estimate errors into ρ0 and thereby into ρn via Eq. 1. Conversely, errors on ρ0 are aliased into orbit errors and other state estimate errors.
Orbit errors are due, in significant part, to force modeling errors. Thus ρ0, ρn, orbit errors, force modeling errors, and other state estimate errors are mutually contaminated. The geodesist wants to call ρn=RANGECP a measurement. Then his measurement modeling errors include orbit errors, force modeling errors, and other state estimate errors. This is a serious estimation problem where measurement errors are confused with, and not distinguished from, force modeling errors.
The estimation of ρ0 must occur at its epoch to, but the geodesist prefers to ignore this fact and refers to ρ0, implicitly and incorrectly, as a bias without epoch, to be estimated as an arbitrary constant. These estimation problems are serious.
Throughput
If one appropriately refers ρ0 to its epoch to for estimation, then Δρ0n=tn−t0 becomes large, and the propagation of state estimate and estimate error covariance matrix between times t0 and tn for each ρn consumes extremely excessive CPU time. This is a throughput problem that destroys, or degrades, any attempt to provide real-time state estimation.
Elimination
All three components of the Unknown Initial Range Problem are eliminated by processing sequential Doppler carrier phase count measurements Nmn with an optimal sequential filter, in preference to estimation of ρ0 and calculation of ρn via Eq. 1 for use as pseudo-measurements.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 illustrates graphically examples of sequential and overlapping time intervals for GPS Doppler phase count measurements;
FIG. 2 illustrates a flow chart of an embodiment of the process of the present invention.
FIG. 3 illustrates a flow chart for autonomous real-time identification and rejection of cycle slips in accordance with the present invention.



DETAILED DESCRIPTION OF THE INVENTION
One Embodiment of the Present Invention
Sequential Doppler Phase Count Measurement Representations
The transmitter on each GPS spacecraft (NAVSTAR) emits radio carrier signals at two L-BAND frequencies fnNV (NV indicates NAVSTAR and n is an integer denoting frequency, either 1 or 2). A USER spacecraft with GPS receiver detects radio carrier signals that are Doppler shifted with multiplicative Doppler shift βij, and are contaminated by random additive ionospheric effects δfnION. Thus the USER spacecraft receives radio carrier signals with frequency βijfnNV+δfnION. The USER receiver knows the nominal transmitted frequency fnNV and (using its own clock) approximates the received radio carrier signal with frequency fnU. A frequency difference is represented by the USER receiver to define (i.e., the subtraction order is defined by each actual receiver, and so our representation here may need a sign change for some receivers) the instantaneous Doppler frequency fnD:

fnD=fnU−(βijfnNV+δfnION)  (2)

with units cycles/second. Positive zero crossings (cycles) of fnD are counted and summed (integrated) by the USER receiver, inclusive of a partial cycle, across a USER defined phase count interval Δt to construct a rational phase count number Njni. Thus Njni consists of the sum of an integer number of cycles and a partial cycle. The GPS receiver phase lock loop incurs random thermal noise (Gaussian white noise) δNjni, independent for each NAVSTAR link i, and independent at each measurement time-tag tj. Thus Njni has a representation with the sum (generalized integral):





  
    
      
        
          N
          j
          ni
        
        =
        
          
            
              ∫
              
                j
                
                  
                    -
                    Δ
                  
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
              
              j
            
            ⁢
            
              
                f
                n
                D
              
              ⁢
              
                
              
              ⁢
              
                d
                t
              
            
          
          +
          
            δ
            ⁢
            
              
            
            ⁢
            
              N
              j
              ni
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




The GPS sequential carrier phase measurement Njni is defined by Eqs. 2 and 3, and has units in cycles.
Thermal Noise Serial Correlation
It is important to notice that when sequential measurements are processed by any orbit determination algorithm, there is no serial correlation due to thermal noise δNjni.
Random Clock Sampling Error
There exists a clock sampling error that generates a serial correlation coefficient value of 0.5 between each pair of sequential adjacent Doppler carrier phase count measurements. However, this correlation does not exist between measurements that are non-adjacent, and is therefore non-observable to the sequential estimator. Non-observable effects have no influence on the sequential estimator.
Calculation of Sequential Doppler Phase Count Measurement Representations
Referring now to Eq. 2, recall that NAVSTAR clock phenomenology is embedded in fnNV, and USER clock phenomenology is embedded in fnU. An explicit relation between fnU and fnNV is required, so define:

δfnU=fnU−fnNV  (4)

We wish to embed δfnU into Eq. 1 in order to identify δfnU explicitly in the Doppler frequency fnD. From Eq. 4 we can write:

fnU=fnNV+δfnU  (5)

and insert fnU into Eq. 2 to get:

fnD=fnNV(1−βij)+δfnU−δfnION  (6)

or we can write:

fnNV=fnU−δfnU  (7)

and insert fnNV into Eq. 2 to get:

fnD=fnU(1−βij)+βijδfnU−δfnION  (8)

The instantaneous Doppler frequency fnD has different representations, that are defined by Eq. 6 and that are defined by Eq. 8. We shall return to this apparent ambiguity later.
Ionosphere
The ionospheric perturbation δfnION contained in the Doppler frequency fnD has the representation:




  
    
      
        
          
            
              
                δ
                ⁢
                
                  
                
                ⁢
                
                  f
                  n
                  ION
                
              
              =
              
              ⁢
              
                -
                
                  [
                  
                    
                      
                        1.344536
                        ×
                        
                          10
                          
                            -
                            7
                          
                        
                      
                      
                        f
                        n
                      
                    
                    ⁢
                    
                      (
                      
                        
                          
                            (
                            Hz
                            )
                          
                          ⁢
                          
                            (
                            cycles
                            )
                          
                        
                        TECN
                      
                      )
                    
                  
                  ]
                
              
            
          
        
        
          
            
              
              ⁢
              
                [
                
                  
                    d
                    
                      N
                      TEC
                    
                  
                  
                    d
                    t
                  
                
                ]
              
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




where the total electron content NTEC is the integral of electron density along the path of the radio signal from transmitter to receiver, where fn is evaluated in units Hz, dNTEC/dt is evaluated in units TECN/s, and where 1 TECN=1 electron/m2. Then δfnION has units of Hz.

Sequential Doppler Carrier Phase Count Measurement Representation for Njni 
To find the representation for Njni, insert Eq. 6 (this choice as opposed to Eq. 8 will be explained later) into Eq. 3 using Eq. 9:




  
    
      
        
          
            
              
                
                  N
                  j
                  ni
                
                =
                
                ⁢
                
                  
                    
                      
                        f
                        _
                      
                      n
                      NV
                    
                    ⁡
                    
                      (
                      
                        
                          ζ
                          ij
                        
                        +
                        
                          I
                          j
                          nNi
                        
                      
                      )
                    
                  
                  +
                  
                    
                      
                        f
                        _
                      
                      n
                      U
                    
                    ⁢
                    
                      I
                      j
                      nU
                    
                  
                  +
                
              
            
          
          
            
              
                
                ⁢
                
                  
                    A
                    ⁢
                    
                      
                    
                    ⁢
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        N
                        TEC
                      
                      /
                      
                        ⌊
                        
                          c
                          ⁢
                          
                            
                              f
                              _
                            
                            n
                            NV
                          
                        
                        ⌋
                      
                    
                  
                  +
                  
                    δ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      N
                      j
                      ni
                    
                  
                
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
      
    
    
      
        (
        10
        )
      
    
  
  
    
      
        
          ζ
          ij
        
        =
        
          
            
              ∫
              
                
                  t
                  j
                
                -
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
              
              
                t
                j
              
            
            ⁢
            
              
                (
                
                  1
                  -
                  β
                
                )
              
              ⁢
              
                
              
              ⁢
              
                d
                t
              
            
          
          =
          
            
              
                Δ
                
                  ρ
                  ij
                
              
              /
              c
            
            -
            
              χ
              ij
            
          
        
      
    
    
      
        (
        11
        )
      
    
  
  
    
      
        
          x
          ij
        
        =
        
          
            ∫
            
              
                t
                j
              
              -
              
                Δ
                ⁢
                
                  
                
                ⁢
                t
              
            
            
              t
              j
            
          
          ⁢
          
            
              {
              
                
                  
                    μ
                    
                      c
                      2
                    
                  
                  ⁡
                  
                    [
                    
                      
                        1
                        
                          r
                          j
                        
                      
                      -
                      
                        1
                        
                          R
                          i
                        
                      
                    
                    ]
                  
                
                +
                
                  
                    1
                    
                      2
                      ⁢
                      
                        c
                        2
                      
                    
                  
                  ⁡
                  
                    [
                    
                      
                        
                          s
                          .
                        
                        j
                        2
                      
                      -
                      
                        
                          S
                          .
                        
                        i
                        2
                      
                    
                    ]
                  
                
              
              }
            
            ⁢
            
              d
              t
            
          
        
      
    
    
      
        (
        12
        )
      
    
  
  
    
      
        
          Δρ
          ij
        
        =
        
          
            
              ∫
              
                
                  t
                  j
                
                -
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
              
              
                t
                j
              
            
            ⁢
            
              
                
                  ρ
                  .
                
                ij
              
              ⁢
              
                d
                t
              
            
          
          =
          
            
              
                ρ
                ij
              
              ⁡
              
                (
                
                  t
                  j
                
                )
              
            
            -
            
              
                ρ
                ij
              
              ⁡
              
                (
                
                  
                    t
                    j
                  
                  -
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    t
                  
                
                )
              
            
          
        
      
    
    
      
        (
        13
        )
      
    
  
  
    
      
        
          Δ
          ⁢
          
            
          
          ⁢
          
            N
            TEC
          
        
        =
        
          
            
              N
              TEC
            
            ⁡
            
              (
              
                t
                j
              
              )
            
          
          -
          
            
              N
              TEC
            
            ⁡
            
              (
              
                
                  t
                  j
                
                -
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
              
              )
            
          
        
      
    
    
      
        (
        14
        )
      
    
  
  
    
      
        A
        =
        
          40.30816
          ⁢
          
            (
            
              
                m
                3
              
              ⁢
              
                
                  Hz
                  2
                
                /
                electron
              
            
            )
          
        
      
    
    
      
        (
        15
        )
      
    
  
  
    
      
        
          β
          ij
        
        =
        
          1
          -
          
            
              1
              c
            
            ⁢
            
              
                ρ
                .
              
              ij
            
          
          +
          
            
              μ
              
                c
                2
              
            
            ⁡
            
              [
              
                
                  1
                  
                    r
                    j
                  
                
                -
                
                  1
                  
                    R
                    i
                  
                
              
              ]
            
          
          +
          
            
              1
              
                2
                ⁢
                
                  c
                  2
                
              
            
            ⁡
            
              [
              
                
                  
                    s
                    .
                  
                  j
                  2
                
                -
                
                  
                    S
                    .
                  
                  i
                  2
                
              
              ]
            
          
        
      
    
    
      
        (
        16
        )
      
    
  
  
    
      
        
          I
          j
          nNi
        
        =
        
          
            ∫
            
              
                t
                j
              
              -
              
                Δ
                ⁢
                
                  
                
                ⁢
                t
              
            
            
              t
              j
            
          
          ⁢
          
            
              [
              
                
                  δ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    f
                    n
                    N
                  
                
                
                  
                    f
                    _
                  
                  n
                  N
                
              
              ]
            
            ⁢
            
              (
              
                1
                -
                
                  β
                  ij
                
              
              )
            
            ⁢
            
              d
              t
            
          
        
      
    
    
      
        (
        17
        )
      
    
  
  
    
      
        
          I
          j
          nU
        
        =
        
          
            ∫
            
              
                t
                j
              
              -
              
                Δ
                ⁢
                
                  
                
                ⁢
                t
              
            
            
              t
              j
            
          
          ⁢
          
            
              [
              
                
                  δ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    f
                    n
                    U
                  
                
                
                  
                    f
                    _
                  
                  n
                  U
                
              
              ]
            
            ⁢
            
              d
              t
            
          
        
      
    
    
      
        (
        18
        )
      
    
  




where:

    
    
        δNjni is receiver thermal noise, independent for each NAVSTAR i
        c is speed of light in a vacuum
        μ is the geocentric two-body gravitational constant
        rj and {dot over (s)}j are magnitudes of USER spacecraft position and velocity vectors at time tj 
        Ri and {dot over (S)}i are magnitudes of NAVSTAR spacecraft position and velocity vectors at time ti 
        {dot over (ρ)}ij is the range-rate {dot over (ρ)}ij=dρij/dt at time t=tj 
        ρij is the range between NAVSTAR spacecraft at time ti and USER spacecraft at time tj 
    
    


In an alternate representation, Eq. 10 can also be written:




  
    
      
        
          
            N
            j
            ni
          
          =
          
            
              
                (
                
                  1
                  
                    λ
                    n
                  
                
                )
              
              ⁢
              
                {
                
                  
                    Δρ
                    ij
                  
                  -
                  
                    c
                    ⁢
                    
                      
                    
                    ⁢
                    
                      χ
                      ij
                    
                  
                  +
                  
                    c
                    ⁡
                    
                      [
                      
                        
                          I
                          j
                          nNi
                        
                        +
                        
                          I
                          j
                          nU
                        
                      
                      ]
                    
                  
                  +
                  
                    [
                    
                      
                        A
                        ⁢
                        
                          
                        
                        ⁢
                        Δ
                        ⁢
                        
                          
                        
                        ⁢
                        
                          N
                          TEC
                        
                      
                      
                        
                          (
                          
                            
                              f
                              _
                            
                            n
                            NV
                          
                          )
                        
                        2
                      
                    
                    ]
                  
                
                }
              
            
            +
            
              δ
              ⁢
              
                
              
              ⁢
              
                N
                j
                ni
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
      
    
    
      
        (
        19
        )
      
    
  
  
    
      
        
          
            
              λ
              n
            
            ⁢
            
              
                f
                _
              
              n
              NV
            
          
          =
          c
        
        ⁢
        
          

        
        ⁢
        
          
            λ
            1
          
          =
          
            19.029367
            ⁢
            
              cm
              cycle
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            λ
            2
          
          =
          
            24.42102
            ⁢
            
              cm
              cycle
            
          
        
      
    
    
      
        (
        20
        )
      
    
  




Then solve Eq. 19 for Δρij to write:

Δρij=λn(Njni−δNjni)+c(χij−[IjnNi+IjnU])−AΔNTEC/( fnNV)2  (21)

Eq. 21 is useful for representation of Δρij from that of Njni.

Two-Frequency Ionosphere Removal
For two-frequency ionosphere removal, let NjEi be defined by:




  
    
      
        
          
            
              N
              _
            
            j
            Ei
          
          =
          
            
              
                
                  
                    
                      f
                      _
                    
                    1
                    NV
                  
                  ⁢
                  
                    N
                    j
                    
                      1
                      ⁢
                      i
                    
                  
                
                -
                
                  
                    
                      f
                      _
                    
                    2
                    NV
                  
                  ⁢
                  
                    N
                    j
                    
                      2
                      ⁢
                      i
                    
                  
                
              
              
                
                  
                    f
                    _
                  
                  1
                  NV
                
                -
                
                  
                    f
                    _
                  
                  2
                  NV
                
              
            
            =
            
              
                
                  γ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    N
                    j
                    
                      1
                      ⁢
                      i
                    
                  
                
                -
                
                  N
                  j
                  
                    2
                    ⁢
                    i
                  
                
              
              
                γ
                -
                1
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
      
    
    
      
        (
        22
        )
      
    
  
  
    
      
        γ
        =
        
          
            
              f
              _
            
            1
            NV
          
          /
          
            
              f
              _
            
            2
            NV
          
        
      
    
    
      
        (
        23
        )
      
    
  




For orbit determination, use the representation NjEi:




  
    
      
        
          
            
              N
              _
            
            j
            Ei
          
          =
          
            
              
                
                  
                    f
                    _
                  
                  1
                  NV
                
                ⁢
                
                  H
                  j
                  
                    1
                    ⁢
                    i
                  
                
              
              -
              
                
                  
                    f
                    _
                  
                  2
                  NV
                
                ⁢
                
                  H
                  j
                  
                    2
                    ⁢
                    i
                  
                
              
            
            
              
                
                  f
                  _
                
                1
                NV
              
              -
              
                
                  f
                  _
                
                2
                NV
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
      
    
    
      
        (
        24
        )
      
    
  
  
    
      
        
          H
          j
          ni
        
        =
        
          
            
              
                f
                _
              
              n
              NV
            
            ⁢
            
              ζ
              ij
            
          
          +
          
            I
            j
            nNi
          
          +
          
            I
            j
            nU
          
          +
          
            δ
            ⁢
            
              
            
            ⁢
            
              N
              j
              ni
            
          
        
      
    
    
      
        (
        25
        )
      
    
  




It is straightforward to demonstrate that ionospheric terms are eliminated due to use of Eqs. 22 and 24.

With respect to ionospheric wave length, Eq. 22 can be written:




  
    
      
        
          N
          j
          Ei
        
        =
        
          
            
              
                λ
                2
              
              ⁢
              
                N
                j
                
                  1
                  ⁢
                  i
                
              
            
            -
            
              
                λ
                1
              
              ⁢
              
                N
                j
                
                  2
                  ⁢
                  i
                
              
            
          
          
            
              λ
              2
            
            -
            
              λ
              1
            
          
        
      
    
    
      
        (
        26
        )
      
    
  




where λnfn=c. Insert Eq. 19 into Eq. 26 and define:

Δ{circumflex over (ρ)}ij=E{Δρij|Njni}  (27)

Let λION denote the two-frequency ionospheric wave length. Then:





  
    
      
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              
                ρ
                ^
              
              ij
            
          
          =
          
            
              
                λ
                ION
              
              ⁢
              
                N
                j
                Ei
              
            
            +
            
              c
              ⁢
              
                
              
              ⁢
              
                χ
                ij
              
            
            +
            
              E
              ⁢
              
                {
                
                  
                    I
                    j
                    Ni
                  
                  ❘
                  
                    N
                    j
                    ni
                  
                
                }
              
            
            +
            
              E
              ⁢
              
                {
                
                  
                    I
                    j
                    U
                  
                  ❘
                  
                    N
                    j
                    ni
                  
                
                }
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
      
    
    
      
        (
        28
        )
      
    
  
  
    
      
        
          λ
          ION
        
        =
        
          
            
              λ
              1
            
            ⁢
            
              λ
              2
            
          
          
            
              λ
              1
            
            +
            
              λ
              2
            
          
        
      
    
    
      
        (
        29
        )
      
    
  




and where approximately λION=10.69534 cm/cycle, given λ1=19.029367 cm/cycle and λ2=24.42102 cm/cycle.

First Differences on NjEi 
With respect to first differences on NE, we define:

ΔNjE=NjEp−NjEq  (30)

where p and q refer to distinct NAVSTARS. Insert Eq. 22 into Eq. 30 to find that f1NVIj1U and f2NVIj2U are differenced out. That is, USER clock phase perturbations Ij1U and Ij2U vanish entirely due to first differences NjE on ionosphere-free Doppler phase count measurements, and due to our choice of Eq. 6 rather than Eq. 8 for modeling USER clock phenomenology. NAVSTAR clock phase perturbations Ij1Np, Ij2Np, Ij1Nq, and Ij2Nq all survive, and all USER receiver PLL thermal noise perturbations δNjNi survive.

Serial Correlation Due to Overlapping Doppler Phase Count Measurements
With respect to overlapping carrier phase measurement, one can inspect Eq. 3 for time tj=ti to write:




  
    N
    1
    ni
  
  =
  
    
      
        ∫
        
          j
          
            
              -
              Δ
            
            ⁢
            
              
            
            ⁢
            t
          
        
        j
      
      ⁢
      
        
          f
          n
          D
        
        ⁢
        
          
        
        ⁢
        
          d
          t
        
      
    
    +
    
      δ
      ⁢
      
        
      
      ⁢
      
        N
        1
        ni
      
    
  




and for time tj=t2 to write:





  
    N
    2
    ni
  
  =
  
    
      
        ∫
        
          2
          
            
              -
              Δ
            
            ⁢
            
              
            
            ⁢
            t
          
        
        2
      
      ⁢
      
        
          f
          n
          D
        
        ⁢
        
          
        
        ⁢
        
          d
          t
        
      
    
    +
    
      δ
      ⁢
      
        
      
      ⁢
      
        N
        2
        ni
      
    
  




The overlapping measurement defined by the sum of the two sequential measurements can be represented with:





  
    
      N
      1
      ni
    
    +
    
      N
      2
      ni
    
  
  =
  
    
      
        ∫
        
          2
          
            
              -
              2
            
            ⁢
            Δ
            ⁢
            
              
            
            ⁢
            t
          
        
        2
      
      ⁢
      
        
          f
          n
          D
        
        ⁢
        
          
        
        ⁢
        
          d
          t
        
      
    
    +
    
      δ
      ⁢
      
        
      
      ⁢
      
        N
        
          1
          ,
          2
        
        ni
      
    
  




where:

δN1,2ni=δN1ni+δN2ni 

The overlapping measurement defined by the sum of L sequential measurements can be defined with:

N1,Lni=N1ni+N2ni+ . . . +NLni 

By induction on the positive integers:




  
    N
    
      1
      ,
      L
    
    ni
  
  =
  
    
      
        ∫
        
          2
          
            
              -
              L
            
            ⁢
            
              
            
            ⁢
            Δ
            ⁢
            
              
            
            ⁢
            t
          
        
        2
      
      ⁢
      
        
          f
          n
          D
        
        ⁢
        
          
        
        ⁢
        
          d
          t
        
      
    
    +
    
      δ
      ⁢
      
        
      
      ⁢
      
        N
        
          1
          ,
          L
        
        ni
      
    
  




where:

δN1,Lni=δN1ni+δN2ni+ . . . +δNLni 

Thermal Noise
It is important to notice that when overlapped Doppler phase count measurements N1,1ni, N1,2ni, . . . , N1,Lni are processed by any orbit determination or geolocation capability, there is significant serial correlation in these measurements due to thermal noise component δN1ni in N1,1ni, N1,2ni, . . . , N1,Lni, also due to thermal noise component δN2ni in overlapped measurements N1,2ni, N1,3ni, . . . , N1,Lni, and, by induction on the positive integers, due to thermal noise component δNkni in overlapped measurements N1,kni, N1,k+1ni, . . . , N1,Lni. The existence of the same thermal noise components in successive overlapped measurements generates serial correlations in the overlapped measurements. Optimal orbit determination requires that these serial correlations be explicitly accounted for when processing overlapped Doppler phase count measurements so that measurement residuals are properly mapped into state estimate corrections in the orbit determination and geolocation capabilities. However, these serial correlations are not accounted for by existing prior art capabilities that process overlapped Doppler phase count measurements.
Absence of Serial Correlation for Sequential Doppler Phase Count Measurements
There is no serial correlation in sequential Doppler phase count measurements due to thermal noise because thermal noise is Gaussian white noise, and white noise is independent with time. That is, white noise has no serial correlation with time.
A method for providing orbit determination and geolocation using sequential GPS Doppler phase count measurements has been described. It will be understood by those skilled in the art that the present invention may be embodied in other specific forms without departing from the scope of the invention disclosed and that the examples and embodiments described herein are in all respects illustrative and not restrictive. Those skilled in the art of the present invention will recognize that other embodiments using the concepts described herein are also possible. Further, any reference to claim elements in the singular, for example, using the articles “a,” “an,” or “the” is not to be construed as limiting the element to the singular.

## Claims

1. A method of refining state estimates of a Global Positioning System (GPS) receiver based upon sequential GPS carrier phase measurements, comprising:
the GPS receiver receives sequential GPS Doppler carrier phase signals from Navigation System with Timing And Ranging (NAVSTAR) transmitter at a defined phase count interval;
the GPS receiver internal clock generates nominal frequencies and derives received frequencies from phase count measurements;
determining a frequency differences to define instantaneous Doppler frequencies in units of cycles per second (Hz);
counting positive zero crossings of the instantaneous Doppler frequencies to measure rational Doppler carrier phase count number in cycles;
calculating an estimated representation of Doppler carrier phase count number;
differencing measured and estimated representations of Doppler carrier phase count number to obtain measurement residuals;
mapping the measurement residuals to state estimate corrections sequentially; and
applying the state estimate corrections to propagated state estimates.

2. The method of claim 1, wherein the state estimates comprise a state estimate structure that contains GPS receiver position, GPS receiver velocity, and GPS receiver clock parameters.

3. The method of claim 1, further comprising eliminating cycle slip by autonomous identification and rejection of any sequential carrier phase measurement that suffers loss of phase lock.

4. The method of claim 3, further comprising continuing processing with a next acceptable sequential carrier phase measurement.

5. The method of claim 1, further comprising processing the state estimates in a sequential filter or a sequential smoother.