Spread spectrum receiver kalman filter residual estimator method

[0001] The present application is related to applicant's co-owned U.S. Pat. No. 6,516,021 entitled Global Positioning System and Inertial Measuring Unit Ultratight Coupling Method, now on the military critical technology list, for official use only.

[0002] The invention was made with Government support under contract No. F04701-00-C-0009 by the Department of the Air Force. The Government has certain rights in the invention.

[0003] The invention relates to the field of communications and navigation systems. More particularly, the present invention relates to communication systems coupled to inertial navigation systems having robust Kalman prefiltering generating residual estimates used for improved signal tracking.

[0004] Two fundamental components of a GPS signal are a pseudorandom noise code and a carrier. Navigation data are derived from this signal by correlating it with code and carrier replicas generated by the GPS receiver. A code tracking loop attempts to drive the code correlation error, the difference between the receiver generated replica and the received signal, to zero by advancing or retarding the replica code rate. Receiver to satellite range, or pseudorange, is computed from the code phase because the code phase is determined by the time required for the signal to travel from the satellite to the receiver. When the code correlation error is small the receiver is said to be in code lock.

[0005] The carrier loop tracks either the carrier frequency or phase by beating a replicated carrier against the received signal to generate an error harmonic. The measured frequency of the error is fed through a filter to form a correction to the carrier replica oscillator. The receiver is said to be in carrier lock when the error harmonic frequency is small, typically less than ten Hertz. Receiver to satellite range rate, or pseudorange rate, is computed from the replicated carrier frequency offset, or Doppler shift, from the nominal carrier frequency. A 50.0 Hertz navigation message superimposed on all GPS signals admits a 180° degree phase ambiguity in the carrier. A frequency locked loop or a phase locked loop implemented as a Costas loop allows tracking with no loss of lock across data bit changes. Both code and carrier tracking loops adjust the phase and frequency of the generated replicas to account for many factors, which include receiver and satellite motion, receiver clock drift, and ionospheric and tropospheric signal delays.

[0006] Nonzero tracking errors appear as unmodeled range and range rate errors that corrupt the navigation solution. With navigation updates computed nominally once per second, only the most recent data taken prior to the measurement time are included in measurements provided to the navigation filter. Most of the correlation data taken over the one-second interval are used to keep the tracking loop in lock. This is determined by the tracking loop bandwidth, which, if set too low, will not track the signal through receiver acceleration or turns.

[0007] In traditional tracking loops, there is the possibility that the errors in the loops are time correlated because there is a delay in the loops as the loops drive the error signal to zero. Depending on the order of the tracking loop, there may actually be a steady state tracking error when the order of the dynamic motion exceeds the order of the tracking loop. This is commonly called the dynamic stress error. Even with inertial aiding to the tracking loop, the tracking error can be time correlated because of time correlated errors in the aiding data. This possibility has limited tracking performance because the tracking loop measurement errors can manifest themselves as errors in the state vector after processing, which in turn can affect future aiding data and possibly cause the IMU navigation and GPS tracking system to generate a wrong but consistent navigation solution. With normal tracking loops, the GPS measurements used by the navigation filter are deliberately spaced out in time to avoid the unwanted temporal correlations of the tracking loops.

[0008] In conventional tracking loops, the measurement noise variance used by the navigation filter is usually estimated from the tracking state and the carrier to noise ratio. There is no guarantee that these values are consistent with the actually realized measurement errors. Hence, a large integration navigation filter may process measurements with an erroneous uncertainty variance of the measurement noise. The navigation filter can be implemented as a Kalman filter.

[0009] In conventional tracking loops, the measurement of pseudorange acceleration is not considered because there is no apparent way to obtain a low noise estimate of the pseudorange acceleration from the tracking loops. Conventional tracking loops do not have an ability to determine the pseudorange acceleration. The use of a pseudorange acceleration measurement has many potential benefits in the observability of user clock and IMU error instability. With a pseudorange acceleration measurement, there is a more direct measurement of higher order errors with increased observability. Normally, these errors are observed through dynamic coupling into pseudorange and pseudorange rate measurements. Hence, without the pseudorange acceleration measurements, the observability of clock and IMU instabilities is limited. This limitation can be reduced using highly accurate user clocks and inertial measuring units.

[0010] Advanced tracking loop designs have been proposed to alleviate some of the loss of lock and reacquisition problems by tightly coupling IMU data with the carrier and code tracking loops and centralizing the navigation and correlation functions. With the navigation filter essentially forming part of the tracking loop, the computational burden becomes very high because the navigation filter that now has many states must also be updated at tens of Hertz rather than the conventional one Hertz.

[0011] With unlimited navigation filter processor throughput, the navigation filter could be operated at extremely high rates with no loss of optimality due to time correlated tracking loop errors. This high rate operation is similar to what has been called a vector delay lock loop. A vector delay lock loop is a method to track all satellites in view simultaneously with one common algorithm. The vector delay lock loop broadens the normal aided and unaided tracking loop design approach to both code and carrier tracking on all in view satellites. The entire algorithm must run at very high processing rates because there is no provision for federated processing. Unfortunately, current processor throughput cannot support Kalman filter rates of several tens of Hertz with the large state vectors required for GPS inertial navigation. Large Kalman filters may be decomposed into one or more federated Kalman filters within a Kalman filter processing architecture. For example, a large Kalman filter can be decomposed into two partitions including a large integration Kalman filter and a high rate optimal prefilter that are more compatible with modern processing speed requirements. The fundamental principle is to decompose the complete formulation into suitable partitions such that the important bandwidths and models are appropriate for each partition. In IMU aided GPS sets designed to date, large navigation filters have not been decomposed into suitable partitions to allow optimal processing of the raw GPS samples with the IMU samples.

[0012] Another approach to aiding the code and carrier correlation process is to design a tracking filter that provides estimates of carrier phase tracking error, rate and acceleration, satellite to user range error, and signal amplitude. The filtered estimates, which are processed at a high rate, are provided to the navigation filter where they are processed at a much lower rate to form the navigation solution. An extended Kalman filter of this design has fifth order or higher order dynamics and an associated matrix Riccati equation that must be integrated with each correlator datum. The computational burden is extraordinarily high.

[0013] Tightly coupled navigation solutions require estimation of the residual error in real time for computation of carrier phase error and code phase error during carrier demodulation and autocorrelation. Computational requirements of Kalman filter residual estimation suffer from computation complexity. These and other disadvantages are solved or reduced using the invention.

[0014] An object of the invention is to provide a method for improved code phase tracking in spread spectrum communication systems.

[0015] Another object of the invention is to provide a method for improved code phase tracking and carrier phase tracking in spread spectrum communication systems.

[0016] Yet another object of the invention is to provide a method for improved code phase tracking and carrier phase tracking in spread spectrum communication systems using a tightly coupled residual estimating Kalman filter.

[0017] Still another object of the invention is to provide a method for improved code phase tracking and carrier phase tracking in GPS spread spectrum communications using a tightly coupled residual estimating Kalman filter in a GPS navigation system.

[0018] A further object of the invention is to provide a method for generating code phase tracking errors and carrier phase tracking errors respectively for code phase tracking and carrier phase tracking of GPS spread spectrum communication signals using a residual estimating Kalman prefilter operating on Ricatti matrices for improved autocorrelation locking in a GPS navigation system.

[0019] The present invention is a communication method directed to carrier phase and code phase tracking using Kalman filter residual estimation well suited for spread spectrum communication systems such as GPS navigation systems. Inphase and quadrature correlation data are processed to provide estimates of carrier phase tracking error and rate and acceleration, as well as code phase tracking error, and signal amplitude. The filter state is a tracking residual applicable to navigation solution correction in ultratight GPS coupling with inertial measurements. The tracking residual estimation drives the code and carrier replica oscillators in tightly coupled correlation loops providing adjusted early and late code replicas and adjusted demodulation carriers for closed-loop code and carrier tracking. The close-loop code and carrier tracking can also be applied to a weak lock navigation systems also using the tracking residual estimation to the code and carrier replica oscillators.

[0020] The method improves interference robustness and navigation accuracy. Improved interference robustness is realized first, through the use of a code tracking and carrier tracking Kalman filter with a Riccati matrix computation process wherein tracking is achieved without the need for consistent signal lock in the usual sense. Second, all correlator early and late inphase and quadrature input to the Kalman filter contribute to form the measurement residual. Third, since correlator early and late inphase and quadrature do not need to be processed in real time to keep tracking loops in lock, application of more advanced algorithms in low signal to noise conditions becomes feasible.

[0021] Navigation accuracy is improved by real time computation of a residual variance using the Ricatti matrix. The real time computation allows a navigation processor to rapidly compute navigation solutions in the presence of noisy signals. Range acceleration availability enables rapid navigation solution convergence with improved observability of signal instabilities. Using correlation data taken over a computational time interval provides accurate tracking error estimates.

[0022] A tracking filter for the code and carrier correlation process of a GPS receiver provides estimates of carrier phase tracking error, rate and acceleration, range error, and signal amplitude. Both inphase and quadrature measurements are processed by the Kalman filter residual estimator. The filter dynamics and associated Riccati equations that are independent of the phase error enable implementation of steady state point designs. The Kalman filter state is a tracking residual applicable as a navigation solution correction well suited for ultratight GPS and IMU coupling. The tracking residual is used to drive the code and carrier oscillators for adjusting early and late code replicas and demodulation carrier signals for improved interference robustness and navigation accuracy. The filter state is an observation measurement residual that can be directly used by the navigation filter as a navigation solution correction. The error dynamics are independent of carrier phase tracking error. Harmonic terms in the phase error are not present in the filter Riccati equation nor in the error dynamics for enabling implementation of a steady state design. The real time Kalman filter processing uses filter state estimates to directly command both code and carrier correlation processes. These and other advantages will become more apparent from the following detailed description of the preferred embodiment.

[0027] An embodiment of the invention is described with reference to the figures using reference designations as shown in the figures. Referring to FIG. 1, a global positioning system (GPS) inertial navigation system includes a navigation system 10 receiving GPS signals and inertial measurement unit samples from a sensor assembly 12 and provides position and velocity data to a control and display unit 14. The navigation system 10 functions as a GPS inertial navigation system tracking GPS signals from a plurality of in view satellites, not shown. The sensor assembly 12 includes an antenna 16 for receiving and providing received GPS signals and includes an inertial measurement unit (IMU) 18 providing the IMU sample signals, both of which are communicated to the navigation system 10. The navigation system 10 includes a downconverter 20 for frequency downconversion of the received GPS signals using a reference oscillator 24 providing an f_o internal frequency reference, and includes an analog to digital (A/D) converter 22 communicating digitized GPS samples to a navigation processor 26. The navigation processor 26 also receives the IMU samples and provides the position and velocity data to the control and display unit 14.

[0028] Referring to FIGS. 1, 2A and 2B, and more particularly to FIG. 2A, the navigation system 10 may be a weak lock navigation system that includes a navigation processor 26 receiving the digitized GPS samples 30 that are communicated to mixers 32 and 34 for providing inphase (I) and quadrature (Q) signals respectively using cos and sin demodulation signals from a carrier numerically controlled oscillator (NCO) 36 using the reference oscillator f_o 24. The cos and sin demodulation signals are replica carrier signals for carrier demodulating the received GPS samples 30. The I and Q sample signals are received by a correlator 40 receiving early (E) and late (L) code replica signals from a chipping code generator 42 receiving a chipping code clocking signal from a code clock NCO 44. The code generator 42 may further include a prompt replica, not shown, for improved sampling. The chipping code may be a direct sequence spread spectrum code. The correlator 40 may operate, for example, at 50 Hz. The correlator 40 provides inphase and quadrature, early and late, GPS correlated sample signals I_E, Q_E, I_L, and Q_L, communicated, for example, at 1K Hertz into a measurement residual Kalman prefilter 46 having, for example, a one Hertz output. The prefilter 46 can receive I and Q sampled outputs from one or more correlators 40 to accommodate dual f₁ and f₂ frequency integration using, for example, two respective NCOs 36 and 44 and two respective code generators 42, for providing an integrated error vector to the single integration filter 48. The Kalman prefilter 46 computes the residual carrier phase errors φ_e, carrier frequency error {dot over (φ)}_e, carrier frequency rate error {umlaut over (φ)}_e, residual code phase error r_e and a matrix P. The carrier phase errors φ_e, {dot over (φ)}_e, and {umlaut over (φ)}_e may, in the context of inertial navigation, respectively correspond to a pseudorange, pseudorange rate and pseudorange acceleration. The carrier frequency error {dot over (φ)}_e is the first time rate of change derivative of the carrier phase error φ_e. The carrier frequency rate error {umlaut over (φ)}_e is the second time rate of change derivative of the carrier phase error φ_e. The Kalman prefilter 46 also computes the amplitude a, code phase error r_e, and a measurement covariance matrix P as part of a measurement vector. The Kalman prefilter 46 computes the measurement state covariance matrix P indicating the uncertainties in the measurement vector of errors for a number of samples m over a major cycle epoch time. The Kalman prefilter 46 can be used as part of a weak lock method using local carrier and code tracking loops for improved code and carrier tracking. The Kalman prefilter 46 can further be used in an ultratight method for GPS pseudorange computation tightly coupled with the correlation process to improve the ability to maintain tracking lock upon the received GPS signals.

[0029] In the exemplar weak lock GPS navigation processor, the correlator outputs I_E, Q_E, I_L, and Q_L are also sampled by a GPS navigation calculator 49. The Kalman prefilter 46 and GPS navigation calculator 49 provide system outputs 51 including the receiver position and velocity output P/V. The code phase error r_e is communicated to the code generator 42 in a local code phase tracking loop. Concurrently, the carrier frequency error {dot over (φ)}_e is fed to the carrier NCO in a carrier tracking loop. The navigation calculator 49 computes the best estimate of the navigation state vector solution defined by the position and velocity data P/V. The GPS navigation calculator 49 computes user clock error estimates (CEE) communicated to a pseudorange and rate calculator 56.

[0030] In the exemplar ultratight GPS navigation processor, IMU samples 52 are communicated to a 100 Hz inertial navigation calculator 54 providing the receiver position and velocity outputs P/V. The Kalman prefilter measurement vector of errors and the state covariance matrix P are communicated to the one Hertz integration Kalman filter 48 having a one second major cycle time period between Kalman filter epochs. The integration Kalman filter 48 in turn repetitively provides estimates of the position and velocity errors in an error state vector (ESV) that are used to compute the position and velocity output data 50 communicated to the control unit 14. The IMU samples 52 are communicated to a 100 Hz inertial navigation calculator 54 that computes the best estimate of the navigation state vector solution defined by GPS receiver position and velocity data P/V. The term P/V indicates the calculated position and velocity data. The position and velocity data P/V is calculated from the IMU samples 52 and the error state vector ESV. At the end of each of the one second major cycles, the Kalman filter 48 updates the ESV and then communicates the ESV to the navigation calculator 54. The IMU samples 52 are repetitively generated, for example, one hundred times a second, and include differential velocity samples ΔV and differential attitude samples Δθ that are communicated to the 100 Hz inertial navigation calculator 54 also receiving the error state vector (ESV) every second from the integration Kalman filter 48. The integration Kalman filter 48 computes the user clock error estimates (CEE) that are communicated to the pseudorange and rate calculator 56. The 100 Hz inertial navigation calculator 54 is used for strapdown inertial navigation processing.

[0031] In the exemplar weak lock navigation processor of FIG. 2A, the navigation calculator 49 provides the best estimate of the navigation state vector solution indicated by the position and velocity data (P/V) communicated to the 100 Hz delta pseudorange and delta pseudorange rate calculator 57 that in turn provides delta pseudorange rate data to the carrier NCO 36 and provides delta pseudorange data to the code NCO 44. The delta pseudorange and delta pseudorange rate calculator 57 receives calculated GPS satellite position and velocity (P-V) data from a 100 Hz GPS satellite position and velocity calculator 58 and receives the user clock error estimates (CEE). The term P-V references the GPS satellite calculated position and velocity data. The clock error estimates (CEE) include a clock phase error and clock frequency error that is the time derivative of the clock phase error. The calculated delta pseudorange is equal to the change in the geometric line of sight range to a satellite adjusted by the clock phase error estimates CEE and certain signal propagation compensation terms, such as the nominal ionospheric and tropospheric delays that are not shown. The delta pseudorange is a change in the geometric range taken within respect to the one Hertz prefilter 46 updates applied to the carrier NCO 36 and the code generator 42. The geometric range is computed from the GPS calculated position and velocity data P-V, and from the navigation calculated position and velocity data P/V. The delta pseudorange rate is equal to the change in the relative geometric velocity to the satellite from the user antenna 16 plus the clock frequency error estimate. The GPS satellite position and velocity calculator 58 receives timing signals from a timer 60 using the reference f_o 24 and receives demodulated ephemeris data from an ephemeris demodulator 62 using the GPS samples 30 to compute the GPS satellite position and velocity data P-V. The weak lock formulation accommodates larger errors between the replica signal and the received signal by observing the errors over many data samples and correcting the replica signal at a much lower rate than the sampling rate of the correlation process inputs.

[0032] In the exemplar ultratight navigation processor of FIG. 2B, the navigation calculator 54 provides the best estimate of the navigation state vector solution indicated by the position and velocity data (P/V) communicated to the 100 Hz pseudorange and pseudorange rate calculator 56 that in turn provides pseudorange rate data to the carrier code NCO 36 and provides pseudorange data to the code NCO 44. The pseudorange and pseudorange rate calculator 56 receives calculated GPS satellite position and velocity (P-V) data from a 100 Hz GPS satellite position and velocity calculator 58 and receives the user clock error estimates (CEE). The term P-V references the GPS satellite calculated position and velocity data. The clock error estimates (CEE) include a clock phase error and clock frequency error that is the time derivative of the clock phase error. The calculated pseudorange is equal to the geometric line of sight range to a satellite adjusted by the clock phase error estimates CEE and certain signal propagation compensation terms, such as the nominal ionospheric and tropospheric delays that are not shown. The geometric range is computed from the GPS calculated position and velocity data P-V, and from the IMU calculated position and velocity data P/V. The pseudorange rate is equal to the relative geometric velocity to the satellite from the user antenna 16 plus the clock frequency error estimate. The GPS satellite position and velocity calculator 58 receives timing signals from a timer 60 using the reference f_o 24 and receives demodulated ephemeris data from an ephemeris demodulator 62 using the GPS samples 30 to compute the GPS satellite position and velocity data P-V.

[0033] The GPS satellite position and velocity calculator 58 is used for satellite GPS position and velocity calculation. The pseudorange and pseudorange rate calculator 56 is used for line of sight pseudorange and pseudorange rate predictions. The delta pseudorange and delta pseudorange rate calculator 57 is used for indicating changes in line of sight pseudorange and pseudorange rate predictions. The carrier NCO 36 is used for generating carrier replica signals for quadrature demodulation of the received GPS samples 30. The code generator 42 is used for early and late code replica signal generation. The correlator 40 is used for integrate and dump correlation of the I and Q carrier demodulated signals correlated with the early and late code replica signals.

[0034] The correlation process is based upon the early and late replica code signal generation by the code generator 42 and upon the replica carrier generation by the carrier NCO 36. As the pseudorange and pseudorange rate data from the calculator 56 are refreshed during each cycle, the replica carrier cos and sin signals from the carrier NCO 36 and the early and late replica codes signals from the code generator 42 are adjusted under close loop processing. The correlator 40 may operate at 50 Hz enabling conventional early and late correlation by providing the early and late I and Q samples to the Kalman prefilter 46. The Q and I samples contain amplitudinal informational that can be used with signal correlation functions to estimate the code and carrier errors. The carrier and code phase offsets of the replica signals as compared to the received GPS sample signal can be determined as offsets knowledge of the signal correlation function. The carrier and code phase offsets are used to estimate the residual errors.

[0035] Referring to all of the Figures, and more particularly to FIG. 3, the measurement residual Kalman prefilter 46 is based on a GPS signal residual estimation process that receives the process inputs I_E, Q_E, I_L, and Q_L that are taken as the correlator early (E) and late (L), and, inphase (I) and quadrature (Q) samples from the correlator 40. The process inputs I_E, Q_E, I_L, and Q_L may be processed at a rate of 1 KHz. The process inputs are represented collectively by the symbol set Y=[I_E,Q_E,I_L,Q_L]. The GPS signal residual prefilter 46 has four process components or algorithms, including a 1 KHz measurement residual generator 82, a 1 KHz state estimate calculator 84, a state propagator 86, and a 5 Hz state error covariance estimator 88. The GPS signal residual estimator process receives the process input y and generates the residual measurements φ_e, {dot over (φ)}_e, {umlaut over (φ)}_ee, r_e, and a.

[0036] The measurement residual generator 82 computes a measurement residual r at 1 KHz from the process input y and the propagated state {overscore (x)}. The propagated state {overscore (x)} has five components where {overscore (x)}=[{overscore (φ)}_e, {overscore({dot over (φ)})}_e, {umlaut over({overscore (φ)})}_e, r_e, a]. The propagated state {overscore (x)} includes the carrier phase error state {overscore (φ)}_e, the carrier frequency error state {dot over({overscore (φ)})}_e, the carrier frequency rate error state {overscore({umlaut over (φ)})}_e, code phase error state {overscore (r)}_e that is related to the range from the receiver to the satellite, and the signal amplitude state ā.

[0037] Each component of the measurement residual r depends on the carrier phase error state {overscore (φ)}_e, code phase error state {overscore (r)}_e, the early to late replicated code offset δ_r, the broadcast data d, and the signal amplitude state ā of the received spread spectrum signal 30. There is a nonlinear relationship between the residual components and the measurement residual r defined by an r matrix. $r=\left[\begin{array}{c}{I}_e-\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e\\ Q_e-\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e\\ I_l-\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e\\ Q_l-\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e\end{array}\right]$

[0038] In the residual matrix, the code correlation function c(•) depends on the pseudorandom noise code chip width λ_ca as defined by a c(ρ) pseudorandom noise code equation. $c\left(\rho \right)=\max \left(1-\frac{\rho }{{\lambda }_{\mathrm{ca}}},-\frac{1}{{\lambda }_{\mathrm{ca}}}\right)$

[0039] In the c(ρ) pseudorandom noise code equation for a coarse acquisition (CA) receiver, the nominal code rate ω_ca is known as the chipping rate where ω_ca=1.023×10⁶ chips/sec. Because the speed of light is C=2.99792458×10⁸ m/sec, the CA code chip width is λ_ca=293.1 m. The true code rate differs slightly from ω_ca due to user to satellite relative motion and atmospheric transmission effects. For the case where the early and late codes are offset from the correlation peak by ¼ chip, δ_r=λ_ca≅73.26 m. The broadcast message data is d=±1 and has a nominal 20.0 millisecond period.

[0040] A state update calculator 84 computes the state estimate x at a rate of 1 KHz. From the measurement residual r, from the propagated state {overscore (x)}, and from a state error covariance estimate P, the state estimate x is equal to {overscore (x)}+K({overscore (x)},P)r, where K is the gain. The gain K depends on components of the propagated state {overscore (x)} and depends on the Riccati matrix P. The gain K is computed to minimize the error variance as a linear Kalman filter gain K({overscore (x)}, P)=Ph^T({overscore (x)})V⁻¹δ_t, where h^T is a transposed linearized observation matrix and V⁻¹ is an inverse observation covariance. The term δ_t is a time step and may be, for example, 0.001 second. The term h({overscore (x)}) is a linearized observation matrix. The term V is the measurement noise variance as used in an observation equation of the y process input.

[0041] The filter measurement y∈R⁴, where y is the process input including observations I_E, Q_E, I_L, Q_L, is modeled as a nonlinear function of the true state x_o, which is not known, with additive noise ν as y=G(x)+ν with G(x) given by a G(x) observation matrix of the true state x_o. $G\left(x\right)=\left[\begin{array}{c}\mathrm{adc}\left(r_e+{\delta }_r\right)\cos {\varphi }_e\\ \mathrm{adc}\left(r_e+{\delta }_r\right)\sin {\varphi }_e\\ \mathrm{adc}\left(r_e-{\delta }_r\right)\cos {\varphi }_e\\ \mathrm{adc}\left(r_e-{\delta }_r\right)\sin {\varphi }_e\end{array}\right]$

[0042] In the G(x) observation matrix, the amplitude a, the code phase error r_e, and the carrier phase error φ_e components can be a true state x_o having a_o, r_o, and φ_o components. The process input y is the sum of the observation matrix G(x) and the additive noise ν, as in an observation equation of the y process input. $y=\left[\begin{array}{c}{I}_e\\ Q_e\\ I_l\\ Q_l\end{array}\right]=\left[\begin{array}{c}\mathrm{adc}\left(r_e+{\delta }_r\right)\cos {\varphi }_e\\ \mathrm{adc}\left(r_e+{\delta }_r\right)\sin {\varphi }_e\\ \mathrm{adc}\left(r_e-{\delta }_r\right)\cos {\varphi }_e\\ \mathrm{adc}\left(r_e-{\delta }_r\right)\sin {\varphi }_e\end{array}\right]+\left[\begin{array}{c}{\upsilon }_{\mathrm{ie}}\\ {\upsilon }_{\mathrm{qe}}\\ {\upsilon }_{\mathrm{il}}\\ {\upsilon }_{\mathrm{ql}}\end{array}\right]$

[0043] The additive disturbance ν∈R⁴ is taken to be a random process with a zero-mean Gaussian distribution and the observation covariance V. The observation covariance V reflects correlation with respect to the early and late, inphase and quadrature correlator output signals but with the additive disturbance ν where ν=[ν_ie, ν_qe, ν_il, ν_ql], where ν is uncorrelated both in time and with respect to the I and Q correlator output signals. With δ=½ being the offset, normalized to the code wavelength, between the early and late correlation signals, then the observation noise covariance V is given by an observation expectation covariance equation. $V=E\left[{\mathrm{\upsilon \upsilon }}^T\right]={\left(a\eta \right)}^2\left[\begin{array}{cccc}1& 0& \gamma & 0\\ 0& 1& 0& \gamma \\ \gamma & 0& 1& 0\\ 0& \gamma & 0& 1\end{array}\right]$

[0044] The observation expectation covariance equation includes a measurement noise structure matrix including the terms 0, 1 and γ that define a measurement noise structure. In the covariance expectation equation, E is a statistical expectation of the additive disturbance ν, γ=1−δ, the scalar a is the signal amplitude is a component of the state x, and 1/η is the signal-to-noise ratio. The observation covariance matrix V can be shown as an exemplar observation matrix V_o. $V_o={\left(0.3\right)}^2\left[\begin{array}{cccc}0.990& -0.004& 0.478& 0.218\\ -0.004& 1.004& -0.227& 0.488\\ 0.478& -0.227& 1.001& -0.005\\ 0.218& 0.488& -0.005& 1.005\end{array}\right]$

[0045] The state estimation error e is defined as e=x_o−{overscore (x)} as in a state estimation error equation.
ė=[A−Kh({overscore (x)})]e−Kν+ω

[0046] In the state estimation error equation, ė is the time derivative of the state error estimate. The process input measurement y is expressed as a function of the state estimation e in a y expansion equation. $y=G\left(x_o\right)+\upsilon =G\left(\stackrel{\_}{x}\right)+\frac{\partial G}{\partial x}{❘}_{\stackrel{\_}{x}}e+\cdots +\upsilon$

[0047] The the expansion equation, y expanded about the propagated state state {overscore (x)}. In the expansion equation, y is the process input from the correlator 40 and the first partial derivative of the observation matrix G(x) is the linearized observation matrix h({overscore (x)}) and is defined by an h({overscore (x)}) partial derivative equation. $h\left(\stackrel{\_}{x}\right)=\frac{\partial G}{\partial x}{❘}_{\stackrel{\_}{x}}$

[0048] The h({overscore (x)}) partial derivative equation can be written in matrix form as an h({overscore (x)}) linearized observation matrix equation. $h\left(\stackrel{\_}{x}\right)=\hspace{1em}\left[\begin{array}{ccccc}-\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e& 0& 0& \stackrel{\_}{a}d{c}^{\prime }\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e& dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e\\ \stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e& 0& 0& \stackrel{\_}{a}d{c}^{\prime }\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e& dc\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e\\ -\stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e& 0& 0& \stackrel{\_}{a}d{c}^{\prime }\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e& dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e\\ \stackrel{\_}{a}dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\cos {\stackrel{\_}{\varphi }}_e& 0& 0& \stackrel{\_}{a}d{c}^{\prime }\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e& dc\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\sin {\stackrel{\_}{\varphi }}_e\end{array}\right]$

[0049] In the h({overscore (x)}) linearized observation matrix equation, c′({overscore (r)}_e+δ_γ)and c′({overscore (r)}_e−δ_γ) are the code correlation sensitivities to range error, and can be expressed by range partial differential equations. $\begin{array}{c}{c}^{\prime }\left({\stackrel{\_}{r}}_e+{\delta }_r\right)=\frac{\partial }{\partial {\stackrel{\_}{r}}_e}c\left({\stackrel{\_}{r}}_e+{\delta }_r\right)\\ c^{\prime }\left({\stackrel{\_}{r}}_e-{\delta }_r\right)=\frac{\partial }{\partial {\stackrel{\_}{r}}_e}c\left({\stackrel{\_}{r}}_e-{\delta }_r\right)\end{array}$ $\frac{\partial c\left(\rho \right)}{\partial {\stackrel{\_}{r}}_e}=\{\begin{array}{cc}0& \rho >{\lambda }_{\mathrm{ca}}\\ \frac{1}{{\lambda }_{\mathrm{ca}}}& 0>\rho >-\left({\lambda }_{\mathrm{ca}}+1\right)\\ -\frac{1}{{\lambda }_{\mathrm{ca}}}& 0<\rho <\left({\lambda }_{\mathrm{ca}}+1\right)\end{array}$

[0050] The state propagator 86 computes at 1 KHz the propagated state {overscore (x)} where the propagated state components are {overscore (φ)}_e, {overscore({dot over (φ)})}_e, {overscore({umlaut over (φ)})}_e, {overscore (r)}_e, ā. The state propagation dynamics are linear, such that, {overscore (φ)}_e=φ_e+{dot over (φ)}_eδ_t+{umlaut over (φ)}_eδ_t², {overscore({dot over (φ)})}_e={dot over (φ)}_e+{umlaut over (φ)}_eδ_t, {overscore({dot over (φ)})}_e={umlaut over (φ)}_e, {overscore (r)}_e=r_e+(λ_L1/2π)δ_t{dot over (φ)}_e, and ā=a, wherein δ_t=0.001 sec and δ_t is the propagation period, and λ_L1≅0.1903 m, that is, λ_L1 is the carrier wavelength. The L1 band carrier frequency is ω_L1 and ω_L1=2.99792458×10⁹ Hz, so that λ_L1=C/ω_L1≅0.1903 m, where C is the speed of light.

[0051] The state error covariance estimator 88 provides a state covariance estimate P at 5 Hz preferably only using the {overscore (x)} propagated state components {overscore (r)}_e and ā. The state covariance estimator 88 may be implemented as a look-up table. Elements of a P scheduling look-up table, that could be a trivial look-up table of only one element, are computed using a state dynamic matrix A. A continuous time representation of the state dynamics is linear with {dot over (x)}=Ax+ω, with {dot over (x)} being the time derivative of the state x, and with A being the state dynamics matrix. $A=\left[\begin{array}{ccccc}0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0\\ 0& \frac{{\lambda }_{\mathrm{L1}}}{2\pi }& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$

[0052] In the state dynamics matrix A, λ_L1≅0.1903 m and is the carrier wavelength. The additive disturbance ω∈R⁵ is taken to be a random process with a zero-mean Gaussian distribution and state covariance W. The state error covariance matrix P is provided by a continuous time Riccati equation.
{dot over (P)}=AP+P^TA+W−Ph^TV⁻¹hP

[0053] In the continuous time Riccati equation, {dot over (P)} is the time derivative of the Riccati state error covariance matrix P, A is the state dynamics matrix, V⁻¹ in an inverse observation covariance, h is the linearized observation matrix, h^T is a transposed linearized observation matrix, W is the state covariance, and P is the Riccati matrix. The measurement noise structure, as in the covariance expectation equation, leads to a product h^T(x)V⁻¹h(x) that does not depend on carrier phase error φ_e nor the carrier frequency error {dot over (φ)}_e. The inverse observation covariance V⁻¹ is given by an inverse covariance expectation equation. $V^{-1}={\left(a\eta \right)}^{-2}\left[\begin{array}{cccc}1& 0& -\gamma & 0\\ 0& 1& 0& -\gamma \\ -\gamma & 0& 1& 0\\ 0& -\gamma & 0& 1\end{array}\right]\frac{1}{1-{\gamma }^2}$

[0054] In the covariance expectation equation, γ=(1−δ) and δ=½ is the normalized early to late code offset δ_r.

[0055] The product h^TV⁻¹h is then defined by a scalar covariance equation. $h^T\left(x\right)V^{-1}h\left(x\right)=\left[\begin{array}{ccccc}{q}_{\mathrm{\varphi \varphi }}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& q_{\mathrm{rr}}& q_{\mathrm{ra}}\\ 0& 0& 0& q_{\mathrm{ra}}& q_{\mathrm{aa}}\end{array}\right]$

[0056] In the scalar covariance equation, the scalar covariance terms q_φφ, q_γγ, q_γa, and q_aa are defined by scalar covariance equations.
q_φφ=a²(c_p²+2γc_pc_m+c_m²)
q_r=a²(c′_p²+2γc′_pc′_m+c′_m²)
q_aa=c_p²+2γc_pc_m+c_m²
q_ra=a[c_pc′_p+γ(c_pc′_m+c_mc′_p)+c_mc′_m]

[0057] In the scalar equations c_p=c(r_e+δ_r), c′_p=c′(r_e+δ_r), c_m=c(r_e−δ_r), c′_m=c(r_e+δ_r). The term h^TV⁻¹h does not depend on message data d=±1 because the message data d is a multiplication factor of each term in h(x) and d²=1. The product h^TV⁻¹h depends only on range error r_e and signal amplitude a. Because h^TV⁻¹h does not depend on the harmonics sin φ_e or cos φ_e, the continuous time Riccati equation has steady-state solutions with respect to constant values for r_e and a. Point design and gain trajectories can be determined with respect to the range error r_e and the amplitude a. A constant error dynamics solution is determined by setting the range error and the signal amplitude to nominal values r_e=0 and a=1.

[0058] Because the product h^TV⁻¹h depends only on the range error, that is, the code phase r_e and the signal amplitude a, a steady-state solution of the continuous time Riccati equation can be expressed as an algebraic Riccati equation.
0=AP+P^TA+W−Ph^TV⁻¹hP

[0059] Nominal closed-loop filter dynamics are found by setting the state covariance matrix W in the observation covariance matrix V to place the closed-loop filter poles at σ(A−Kh)={−335, −152, −106, −8.7±4.9i). $W=\left[\begin{array}{ccccc}10& 0& 0& 0& 0\\ 0& {10}^3& 0& 0& 0\\ 0& 0& {10}^6& -3\times {10}^5& 0\\ 0& 0& -3\times {10}^5& {10}^5& 0\\ 0& 0& 0& 0& 1\end{array}\right]$  V=10⁻⁴I₄

[0060] In the state covariance matrix W, the cross terms are found to improve the damping in the carrier frequency error rate {umlaut over (φ)}_e response. The process outputs are the state estimate x and state error covariance estimate P with the state estimate x including the φ_e, {dot over (φ)}_e, {umlaut over (φ)}_e, r_e, and a components.

[0061] The use of the algebraic Riccati equation enables state error covariance estimation in real time in the Kalman prefilter 46. The Kalman prefilter 46 provides code phase error and carrier phase error signals for adjusting the code phase and carrier phase for improved code phase tracking and carrier phase tracking. Those skilled in the art can make enhancements, improvements, and modifications to the invention, and these enhancements, improvements, and modifications may nonetheless fall within the spirit and scope of the following claims.