System and method for mitigating severe multipath interference for geolocation and navigation

Abstract

A method for increasing accuracy of geolocation based on timing measurements of RF signals. The method includes: (a) receiving at a user one or more RF signals transmitted over time, (b) correlating the received signals with a plurality of replica waveforms, (c) creating a first set of estimates of the time or arrivals of the signals at the user, (d) using the first set of estimates to create a plurality of models, wherein each model assumes that some of the measurements are affected by multipath or interference, (e) testing the models to estimate the likelihood of timing measurements for each model, and (f) determining the set of LOS measurements as one associated with a model having the largest likelihood.

Description

BACKGROUND

Technologies for personal and vehicular navigation are rapidly developing due to the availability of Global Navigation Satellite Systems (GNSS) such as the United States' Global Positioning System (GPS), the Russian GLONASS system, and the European Galileo system. These systems however are designed for environments where a clear line of sight (LOS) exists between the user receiver and the GNSS satellites. Using trilateration methods, a user can convert the ranging measurement obtained from the LOS signal into an estimate of the user's position.

Indoor and urban canyon navigation presents two main challenges for GNSS: low signal strength and severe multipath interference. Assisted GPS (AGPS) and Ultra-Tight Coupling (UTC) are current methods for tracking weak signals in the GNSS receiver operating environment. Although there are several existing techniques for mitigating multipath interference, these known techniques are not effective to adequately mitigate the multipath in the indoor and urban environments. Multipath interference exists when a receiver receives signals reflected from nearby terrain or structures, such as buildings or walls. The received multipath signal is commonly referred to as a non-line of sight (NLOS) signal. Multipath signals always arrive “late” compared to the LOS signal, thus creating an error in the measured range and corrupting the user position estimate. This problem is especially acute indoors and in urban canyons, due to multiple reflection-generating objects which surround a user (building walls, furniture, cars, etc.).

There are several known methods to mitigate multipath, but these methods help only in relatively benign multipath environments. For example, if the strength of muitipath signals is not large compared to the strength of the LOS signal, and multipath delays are not too small, narrow correlators, strobe correlators and similar methods may be used to effectively isolate or remove the multipath signals. However, these techniques are effective only if the multipath delays are on the order of or larger than the inverse signal bandwidth, e.g., at least 0.1 chip length, which is 30 meters for GPS civilian signal. However, in the indoor and urban environment, there are many multipath-generating surfaces at distances smaller than 30 meters, which makes these techniques not effective for those environments.

Another known technique for multipath mitigation is called “Multipath Estimation Delay Lock Loop” (MEDLL) from NovaTel, Inc. in Calgary Alberta, Canada. A MEDLL receiver has many correlators which integrate the satellite signal at different delays (compared to typically three correlators for a traditional receiver) against the known code of the transmitted signal. The result is a profile of a GNSS satellite signal's correlation with the code replica, sampled at an array of points. LOS and all multipath components contribute to this profile, forming a complex signature. The MEDLL can discriminate the individual signal components, thus can be used to isolate LOS from multipath. This method is effective if there are few dominant signal components, and if one of the dominant components is the LOS signal. If there are many multipath components and/or if LOS is not present or weak, the MEDLL receiver does not yield a reliable LOS measurement.

SUMMARY

The present disclosure presents novel and advantageous systems and methods for discriminating between LOS and NLOS signal paths in a radio frequency (RF) receiver such as a GNSS navigation receiver. In one aspect, the present disclosure provides for a method referred to as Synthetic Aperture Line of Sight Assessment, or SALSA. The SALSA method uses the direction of arrival for a signal to discriminate between LOS and NLOS and is especially beneficial in environments where the amplitude of the multipath signals exceeds that of the LOS signal, and where the multipath delays are within the inverse of the signal bandwidth.

In another aspect, the present disclosure provides for a method referred to as Genetic Algorithm for Multipath Elimination, or GAME. The GAME method identifies and isolates signals that have been incorrectly identified as LOS signals. For a user navigating in urban or indoor environment, some LOS signals may be completely blocked by building walls, etc. Thus, there is no guarantee that a pseudorange measurement is not corrupted by the multipath, even if the most sophisticated algorithms and hardware is used for processing each signal. Therefore, one needs to have a method to identify LOS paths among received signals. If there are multiple signals received, then it is possible to check the consistency of several measurements to find if any of them are corrupted. For example, RAIM (Receiver Autonomous Integrity Monitoring) is often used to find if there is one faulty measurement in a set, e.g. due to a GPS satellite failure. In the case of urban and indoor navigation, there may be multiple faulty (biased) measurements at each epoch, and the partition between biased and non-biased measurements may change from one epoch to another due to user motion. Thus, the algorithm to identify biased measurements must be able to react to changes in the measurements at a rapid pace, and must be able to identify multiple faults at each epoch. GAME is designed to serve this purpose.

The SALSA and GAME methods are complementary and are designed to extract a LOS measurement in a multipath environment if the LOS signal is present and if the user is moving, or if a signal is received over time. SALSA and GAME may also be used to support Direction Finding algorithms for emitter geolocation.

In another aspect, the present disclosure provides for a method referred to as Weighted Average Functionality For Limiting Error Sources (“WAFFLES.”) WAFFLES uses linear combinations of TDOA measurements in such way, that it largely cancels effects of timing, calibration, and geolocation errors. An alternative and complementary navigation method to GPS means is the use of various terrestrial signals rather than GPS signals. For example, signals for digital TV, from cell phone base stations, and Wi-Fi stations can be used to navigate in the areas where such signals are available. Since these signals are not designed for navigation, they lack several important navigation features, most importantly, accurate timing. This makes it necessary to use reference stations to timestamp signal features. Information about timestamps can be transmitted to the user via a communication channel. By comparing timing of the same features received by the user equipment with that at the reference station, the user can form a TDOA measurement, which can be used for navigation. One problem facing this concept is that location of signal sources is not always known with high accuracy. For example, location of a cell base tower may be obtained from original engineering drawings, but there is always a possibility that location data will have some inaccuracies. Errors in locating signal sources directly impact user navigation accuracy. Another problem deals with calibration and timing errors at the reference stations. TDOA at the user assumes that timing measurements at the reference stations are accurate. This may be not the case if calibration and synchronization of reference stations is inaccurate. These two problems are linked. Indeed, if data on location of the signal sources is not reliable, it may be possible to geolocate these signal sources using the same reference stations as those providing timing information to the user. However, any timing and calibration errors at the reference station will introduce errors in the source geolocation. Thus, TDOA measurements at the user will be impacted by timing and calibration errors in two ways: (1) directly, and (2) due to errors in source geolocation. WAFFLES addresses this problem by using linear combinations of TDOA measurements in such way, that it largely cancels effects of timing, calibration, and geolocation errors.

In another aspect, the present disclosure provides for a method of navigation in an environment where the navigation signals may be obstructed, referred to as LEAF. RF navigation under foliage presents its own set of challenges. The same may be true for some other complex environments, where signal propagation to the user is subject to diffused scattering. The earliest arrival of the signal is due to the direct signal component, and its measurement is the goal of the pseudorange estimate by the receiver. However, this component arrives at the receiver along with all scattered signal components, which cannot be separated from the direct component and which corrupt the pseudorange measurement.

The LEAF method is designed to estimate pseudorange using statistical assumptions about the propagation medium. Individual scatterers (e.g., leaves and twigs in the foliage) produce random phases and amplitudes of the signal at different delays; however statistical properties of the received signal are predictable and can be used. The method finds the most likely signal delay from correlation measurements obtained by an array of correlators.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a simplified pictorial representation of a multipath signal environment.

FIG. 2 is a simplified pictorial representation of signal correlations in a multipath environment.

FIG. 3 is a simplified flow diagram of one embodiment of the Synthetic Aperture Line of Sight Assessment (SALSA) method.

FIG. 4 is a simplified pictorial representation of the Genetic Algorithm for Multipath Elimination (GAME) algorithm.

FIG. 5 is a simplified pictorial representation of the effects of timing and calibration errors on source geolocation.

FIG. 6 is a simplified pictorial representation of the Weighted Average Functionality For Limiting Error Sources (WAFFLES) algorithm.

FIG. 7 is a simplified pictorial representation of the LEAF algorithm.

DETAILED DESCRIPTION

SALSA

Determining the Direction of Arrival (DOA) of a signal is a useful method to discriminate LOS from multipath. The LOS signal generally comes from the direction of the transmitting satellite, while the multipath signals come from the direction of the reflector. The DOA discriminator is not constrained by assumptions of large delay and relatively weak multipath. Traditional methods for separating signals by DOA require an antenna array. Of course, it is not always practical for a user to carry an antenna array. In the mobile receiver environment, is desirable to use a small, one-element omni-directional antenna. The present disclosure describes a method for using an omni-directional antenna for DOA by utilizing an antenna array processing effect created by exploiting user motion to synthesize the array aperture.

With reference to FIG. 1, a walking user 100 moves at the speed of approximately 2 m/s beginning at time t₀. At time t₀, user 100 may receive a multipath signal 110 and a LOS signal 120 from GNSS transmitting satellite 130. At time t₁, the user position has changed by many wavelengths (e.g., wavelength of the GPS L1 frequency is 19 cm), and may receive a multipath signal 140 and a LOS signal 150 from satellite 130. In one aspect of the present disclosure, a receiver will process the received signals at times t₀ and t₁ (and at additional times t_n) to separate LOS from multipath and, in effect, create a large synthetic aperture array. In the case of a user in a vehicle traveling at an increased speed, the resulting synthetic aperture array will be larger.

In SALSA operation, a GNSS receiver can correlate an incoming satellite signal with a known code replica to determine complex (I and Q) correlations. Each correlation requires integrating over some period of time, □t. This time interval is referred to in this disclosure as a sub epoch. A bank of N correlators, each having an equally spaced code delay, produces N correlation measurements for each sub epoch. The receiver accumulates correlation measurements from multiple sub epochs over each epoch T. Thus, by the end of the epoch, there is a two-dimensional array of correlations, where one dimension is the code delay, and another dimension is time (sub epoch).

For example, consider a simple case when the user with a receiver moves with a constant velocity during an epoch (this means at a constant speed and along a straight line). Since LOS and multipath signals transmitted by the transmitter arrive at the user from different directions, they will have different rates of phase change at the receiver as illustrated in FIG. 2.

FIG. 2 is a simplified pictorial representation of signal correlations in a multipath environment. Signal correlations are measured in two dimensions: at different time instances (sub-epochs) 200 as well as at different delays of signal replica 210. The top portion of the figure illustrates a two-dimensional matrix of resulting complex correlations, each represented with an arrow 220. A spectral analysis algorithm 230 transforms the time domain into the frequency domain 240. As a result, signals reaching the receiver via different paths become separated, and each path exhibits a maximum correlation at a particular delay, shown as a tip of a triangle 250. The earliest arrival is likely to be the line-of-sight.

If a Fourier transform is performed in the time domain (across sub epochs) for each replica delay, different paths will show in different phase rate bins. In this example, signal paths can be separated by their Doppler rates. For each phase rate bin, there may still be some correlation profile in the replica delay domain that may be due to relatively few paths which have the corresponding phase rate (e.g., arrive from approximately the same corresponding direction). However, only a few paths will contribute to a particular phase rate bin, and therefore these paths can be resolved by some conventional techniques, such as MEDLL. The result of this procedure will be multiple, separately identified paths, each with its own phase rate and delay. The earliest path (smallest delay) is a candidate for the LOS.

In real world applications, a user may not move with a constant velocity and along a straight line. The disclosed methods can be applied to any trajectory if the LOS components are added coherently, and multipath components are added destructively. This can be done if the relative trajectory is approximately known during the epoch, which generally requires using an inertial measuring unit (IMU). In this case, the correlation phases for each correlation measurement may be rotated in such a way that LOS components should have the same phase. This step will be referred to as motion compensation.

After the correlation phases have been motion compensated, a simple way to isolate the LOS is to sum correlations over all sub epochs for each replica delay. In a motion-compensated frame, LOS corresponds to the zero phase rate (user “standing still” in the framework tied to the user) and the Fourier transform at the zero phase rate is the same as integration over time. This example, though not an optimal implementation, shows that SALSA has some similarities with integration.

There are several primary benefits of doing SALSA as compared to the plain long signal integration. First, SALSA accounts for imperfections in motion compensation. In an ideal case, having an accurate estimate of relative user trajectory during the epoch from an inertial measuring unit (IMU), and an accurate internal clock, perfect motion compensation would be possible, and the LOS component is guaranteed to be in the zero phase rate bin. In this idealized case, integration would achieve the same result as the Fourier transform. However, the estimate of the user trajectory is imperfect, and the user clock has a drift. In this realistic case, the true signature of the LOS may appear not in the zero phase rate bin, but somewhere in its vicinity, and a Fourier transform may be used to find it.

Second, the position of the LOS-induced maximum in the phase rate domain is an indication of the errors in the clock drift and trajectory estimate. In essence, it is the phase rate residual. When multiplied by the epoch duration, it becomes the residual in ADR (accumulated delta range), which is a valuable measurement in itself, commonly used by many navigation systems.

Third, SALSA assists in isolating the LOS component from multipath. Some multipath components may be located in phase rate bins which are poorly separated from the location of the LOS bin. Each component will show up as a maximum in the phase rate (Fourier) domain. If maxima are not well separated, a tail of one maximum may contaminate a measurement for another maximum. A plain time-integration ignores this effect, and suffers from multipath correlations “leaking into” the integration result. SALSA uses windowing (e.g., Kaiser-Bessel window) to improve path separation.

Fourth, SALSA may be used with beamforming and direction finding to improve isolation of the LOS signal. Each path contributes a spectral peak in the phase rate (Fourier) domain. To find and isolate a particular peak, such as that for LOS, there are many sophisticated methods to do so in the beamforming and direction-finding (DF) technologies, examples include MUSIC, Maximum Entropy Methods, and weighted subspace fitting. By re-casting a navigation problem in synthetic aperture terms, the entire arsenal of beamforming and DF techniques can be used with the SALSA method.

SALSA MUSIC

In another aspect of the present disclosure, SALSA uses a modified version of the MUSIC algorithm to analyze correlations across RF channels and delays, and thus does not require user motion.

In mathematical terms, A matrix C* can be formed

$\begin{array}{cc}{C}^{*}=\left(\begin{array}{cccccc}1& c_1& \dots & c_n& & 0\\ & \ddots & \ddots & & \ddots & \\ 0& & 1& c_1& \dots & c_n\end{array}\right),\left(m-n\right)⨯m& \left[1\right]\end{array}$ where m is the number of elements in the array, and n is the rank of the signal subspace (i.e., number of spectral components we are looking for). Coefficients c₁ are unknowns and will be solved for.Let S=[s₁, . . . s_n], m×n be a matrix formed by eigenvectors of the signal correlation matrix. Then it can be shown thatC*S=0.  [2]

Terms in this matrix equation can be re-arranged to rewrite it in the form:Φc=μ  [3]where the (m−n)n×n matrix Φ and the (m−n)n×1 vector μ are entirely determined from the elements of S, and where c=[c₁, . . . c_n]. If the sample version of S is used, then c can be treated as unknowns in a linear (over-determined) system of equations and can be solved for.

Next, a polynomial can be formedA(z)=1+c₁z⁻¹+ . . . +c_nz⁻ⁿ  [4]

Roots of that polynomial for z⁻¹ happen to correspond to components of the signal:

$\begin{array}{cc}A\left(z\right)=\prod _{k=1}^n\left(1-e^{{\mathrm{i\omega }}_k}{z}^{-1}\right)& \left[5\right]\end{array}$ Although it may appear convoluted, it is recognized that this method works for correlated and for coherent signals, i.e. assumption of signal independence is not necessary.

The above algorithms are known and discussed in Spectral Analysis of Signals, P. Stoica and R. Moses, 2005.

The present disclosure adapts these algorithms for use in the specific application. In the present disclosure, all signals (LOS and multipaths) come ultimately from the same source. Thus, all signals are fully coherent. This is a somewhat extreme case, and the mathematical treatment for this case follows.

First, with respect to a signal correlation matrix, for any pair of array elements, the corresponding element of the correlation matrix is defined as:

$\begin{array}{cc}{M}_{ij}=E\left\{y_i·y_j^{*}\right\}.& \left[6\right]\end{array}$ Signal samples y₁ are a sum of multiple signal components. In the case of fully coherent signal components, all components (except noise) vary by the same phase with time, and therefore

$\begin{array}{cc}E\left\{y_i·y_j^{*}\right\}=y_i·y_j^{*}.& \left[7\right]\end{array}$ In other words, averaging over time does not do any good for fully coherent signals. This has big implications. It is easy to check by direct substitution that for any vector xM x=g y  [8]where g is a constant. Thus, matrix M projects any vector on vector y, which means that M has only one non-zero eigenvector, which is y. (If there is noise, there will be noise space eigenvectors, but the present disclosure is concerned with the signal space.)

This totally invalidates the traditional MUSIC method. However, with the modified MUSIC method presented above, the l-th equation in the system of equations (2), becomes greatly simplified by the fact that there is only one non-zero eigenvector and therefore only one column in matrix S:

$\begin{array}{cc}{\left(C^{*}S\right)}_l=\sum _{q=0}^n{c}_q{y}_{q+1}=0& \left[9\right]\end{array}$ where notation c₀=1 is used for brevity. This assumes that samples y₁ comprise n signal components arriving from different directions. If there is a uniform linear array (ULA), then

$\begin{array}{cc}{y}_j=\sum _{k=0}^{n-1}{a}_k{e}^{{\mathrm{i\omega }}_{k}j}.& \left[10\right]\end{array}$ Then

$\begin{array}{cc}{\left(C^{*}S\right)}_l=\sum _{q=0}^n{c}_q\sum _{k=0}^{n-1}{a}_k{e}^{{\mathrm{i\omega }}_k\left(q+l\right)}=0.& \left[11\right]\end{array}$ Changing the order of summation to get:

$\begin{array}{cc}{\left(C^{*}S\right)}_l=\sum _{k=0}^{n-1}{a}_k{e}^{i{\omega }_{k}l}\sum _{q=0}^n{c}_q{e}^{i{\omega }_{k}q}=0& \left[12\right]\end{array}$ Equations [12] for different l form a linear system with linearly independent coefficients α_ke^iΩsl. Thus, they can only be satisfied if

$\begin{array}{cc}\sum _{q=0}^n{c}_q{e}^{{\mathrm{i\omega }}_{k}q}=0& \left[13\right]\end{array}$ for any value of k.Denoting e^iΩt=ζ_k, then [13] is a polynomial for ζ_k, which must have n roots. Since the zero-order term is 1 (see definition of c_q), equations [13] can be rewritten in the form:

$\begin{array}{cc}\prod _{k=1}^n\left(1-b_k{\zeta }_k\right)=0& \left[14\right]\end{array}$ where b_x are defined via roots of polynomial [13]. Equation [14] must hold for any k, which means that b_k are simply ζ_k⁻¹ (not necessarily in the same order).

Thus, the method can be summarized as follows:

• 1. Get samples y_i. Average them over time to decrease noise.
• 2. Form equations

$\sum _{q=0}^n{c}_q{y}_{q+1}=0$ and solve them for c_q. Assume c₀=1.

• 3. Find roots of the polynomial, which has c_q for coefficients.
• 4. Roots will be equal to e^−iωi where ω_k defines a direction of arrival for the k-th signal components in a ULA.

It is not necessary to find eigenvectors, which is good in terms of simplicity and computational efficiency.

With reference to FIG. 3, the use of one embodiment of the SALSA method is disclosed. A receiver receives and digitizes an RF signal 300. For each receiver channel, the following sub-steps are performed:

• (a) The received signal 300 is correlated with a replica waveform using a bank of correlators 310 with the following characteristics:
  • (1) each correlator is placed at a different delay;
  • (2) the correlator delays are computed using the estimated user position, user clock error, satellite clock error and satellite position (this is known as ultra-light coupling in integrated GPS-INS systems), and account for the LOS propagation;
  • (3) the delays are set in such a way that the expected LOS signal has non-zero correlations with at least some replicas; and
  • (4) the correlations are estimated at regular time intervals (sub-epochs).
• (b) With the data generated above and illustrated (“Correlations (complex)” in FIG. 2), a two-dimensional matrix of correlations 320 is formed with the following characteristics:
  • (1) one dimension is the sub-epochs, and the other dimension is delay; and
  • (2) the 2D matrix represents measurements at a single epoch.
• (c) using the navigation estimate of the user trajectory, the 2D correlation matrix is corrected for the user motion 330 with the following process:
  • (1) the phases of correlation measurements are rotated by an amount consistent with the user motion the phase of correlation for the LOS signal will be static across sub-epochs after this correction); and
  • (2) the correlation matrix is interpolated in the delay domain to compensate for changes of the delay from one epoch to another;
• (d) A spectrum estimation algorithm 340 is applied in the time domain to the correlation matrix to separate out the different signal paths by direction of arrival (shown in FIG. 2 as “FFT over sub-epochs”). Example spectrum estimation algorithms include windowing and Fourier transform, principal component analysis, singular spectrum analysis, maximum entropy method, MUSIC method, and Modified MUSIC method. An example pictorial of the resultant Replica Delay-Phase rate plot, along with the observation of the earliest arriving signal path is shown in FIG. 2. In addition to the time domain, spectral estimation algorithms may also be applied to frequency and channel domains for additional multipath mitigation.
• (e) The LOS path is selected 350 using one or more of the following criteria:
  • (1) LOS is the earliest arriving path;
  • (2) LOS is the slowest varying principal component; and
  • (3) LOS corresponds to a nearly zero phase rate in the Fourier domain.
• (f) The pseudorange (PR) delay and Doppler for the LOS component are estimated 360.

The PR and Doppler measurements are used to estimate user position using standard navigation techniques such as Kalman filtering 370. High-rate IMU data 380 can be used to augment the measurements in the navigation filter 370.

Modified MUSIC Implementation for a Multi-Channel Signal

The SALSA MUSIC method is based on two observations. First, isolating LOS and multipaths in a multi-channel signal (such as OFDM) by a moving user is mathematically equivalent to the problem of spectral estimation in two dimensions. The first dimension is the time domain, and the second dimension is the channel domain. Second, the conventional MUSIC algorithm does not work well for the problem at hand, because LOS and multipath signals are highly correlated (basically, fully correlated). However, the modified MUSIC algorithm works.

In one embodiment, the SALSA MUSIC method can be summarized as follows:

• 1. Each frequency channel is tracked and demodulated individually. The replica signal for each channel is correlated with the incoming signal to obtain I and Q measurements (complex correlations). This creates a 2D matrix of data (number of sub-epochs×number of frequency channels complex values).
• 2. A reference station supplies its own measurements of the complex amplitude for each channel. The user rotates/scales measurements to compensate for any clock drifts at the signal source and for any variations in amplitude (e.g., due to power control at the signal source).
• 3. Using the estimated user position, the algorithm corrects for the user motion, which means that phases of measurements are rotated.
• 4. The Modified MUSIC algorithm is applied to one channel (the middle one) in the 2D matrix of the data. This algorithm determines the spectrum (which can be mapped to the directions of arrival for different paths in the synthetic aperture array). The algorithm finds a fixed number of paths, regardless of the true number of paths in the data. For example, it may look for N=10 paths. If the RF propagation channel contain fewer than 10 paths, the algorithm will find spurious paths (described below); if the RF propagation channel contains more than 10 paths, some weaker ones will not be found.
• 5. Note that if the multipath environment contains two parallel walls, then reflections from these walls will have the same DOA hut different delays. This is just one example, which shows that it is not enough to resolve DOAs to isolate every path. The second dimension should be considered. Next, the algorithm loops over the N paths, which were found in the DOA domain:
  • a. For the current path in the loop, it treats that path as a signal and all other paths as interference. It computes beamforming weights, which maximize the gain in the direction of the current path, while placing nulls in the directions of the other N-1 paths. These beamforming weights are applied to data from all frequency channels.
  • b. The result of the previous step is M complex numbers, which are in the channel domain. Modified MUSIC is applied again, this time in the channel domain. Fewer paths are now sought for two reasons. First, time domain/DOA analysis has separated many paths already, and for this particular DOA there should be one or at most few paths left, which were not resolved. Second, having only M measurements, we would not be able to find many paths anyway. The algorithm currently looks for K paths.
  • c. For each of the K paths, the algorithm computes beamforming weights, places max gain on the current path and null on K-1 other paths.
• 6. The total result of all this computation is N×K path combinations (K paths found in the channel domain for each of the N paths in the time domain). Each combination has three important values associated with it:
  • a. Spectral frequency in the time domain. This value is mapped to the delta range measurement residual for that path. (Or to DOA, if desired.)
  • b. Spectral frequency in the channel domain. This value is mapped to the delay measurement residual for that path.
  • c. Amplitude of the signal on that path, as computed by the beamforming. Note that beamforming is applied twice (in the DOA domain and then in the delay domain.)
• 7. The next task is to separate the seeds from the chaff, i.e. to eliminate spurious paths. The reasoning goes like this: if a beam is formed on the spurious path and all other paths are nulled (which includes all real paths) then the resulting amplitude of the signal will be relatively small. This applies to both stages of the beamforming. Thus, 30 path combinations are sorted by their amplitudes, and some in the bottom of the list are cut off.
• 8. Finally, the LOS path is selected. The selection is based on three criteria:
  • a. Delay residual is not too large (within a configurable threshold)
  • b. Delta range residual is not too large (within a configurable threshold)
  • c. First arrival among those satisfying criteria a and b above.

Thus, in one aspect inertial measuring units (IMU) and UTC can be used to get multiple measurements for a single epoch. Multiple correlation measurements can be collected over one epoch, e.g. 100 measurements, separated by 10 ms are collected over one second. Rather than integrating or averaging these measurements (as it is done by techniques for long integration of GPS signals), measurements are processed as a set.

In another aspect, using data from the IMU and UTC, measurements are corrected for the user motion, in such way that for the LOS path the phase of the signal remains constant or changes linearly.

In another aspect array processing techniques are used to isolate paths. Array processing techniques are applied to the set of measurements to isolate individual signal paths. These techniques may include windowing and Fourier transform, principal component analysis, MUSIC and Modified MUSIC.

In yet another aspect, the paths can be cycled through and delays can be determined. With the knowledge of DOA for all paths reaching the receiver, the algorithm can apply beamforming to all or some of the paths and determines each path's delay. The beamforming can be performed using the following steps:

• • a. Compute beamforming weights to put gain on the desired path (i.e., DOA) and/or nulls on other paths' DOA
  • b. Apply weights to multiple correlators, which have different replica delays or track different frequency channels of the same signal
  • c. Determine delay using any known method, such as early-minus-late correlators, MEDLL, strobe correlator, or super-resolution in the delay or channel domain.

In another aspect, the LOS path is chosen. Delay estimates and other characteristics can be used for each of the paths to find the LOS path. Any combination of the following criteria can be used:

• • a. TOA of the LOS path is the earliest for all paths
  • b. Post-beamforming signal amplitude is above some threshold level. The threshold may be determined by the SNR and/or noise level
  • c. DOA of the LOS path is within some tolerance from the estimated direction to the signal source.

In another aspect, the LOS delay can be output to navigation software.

GAME

SALSA alone may not be sufficient to produce good measurements in certain environments. SALSA outputs measurements for first-arrival signals. These measurements may be processed by a Kalman filter; however Kalman filters are vulnerable to any biases in the input measurements. For example, if a signal from a particular GNSS satellite has no LOS component (or the LOS component is too weak to be detected) then the first arrival will correspond to a multipath component, which is necessarily delayed. This creates a bias in some measurements, and will ruin the Kalman filter performance. GAME is designed to identify and eliminate faulty (i.e., non-LOS) measurements.

One well known prior art method used in identify faulty measurements is called Receiver Autonomous Integrity Monitoring (RAIM). RAIM is a popular algorithm for airplane navigation, where data integrity is important. Current RAIM methods check consistency of measurements to identify and eliminate one faulty measurement (there also was some work to extend RAIM to multiple faults). GAME is also directly to eliminating faulty measurements, but it is a substantial improvement upon RAIM in the following ways:

• 1. GAME is able to identify multiple faults,
• 2. GAME makes use of history of measurements, and
• 3. GAME yields a quantitative measure of the integrity of measurements in the form of measurement likelihood.

GAME builds on two known algorithms: Interacting Multiple Models (IMM) and Genetic Algorithms. In one embodiment, GAME simultaneously tracks multiple models; each model characterized by a particular allocation of LOS/non-LOS flags. For each model, GAME computes a Bayesian likelihood of that model being true. GAME assumes that multipath environment is dynamic, and accounts for a possibility of flips in LOS/non-LOS flags. Each epoch, GAME performs two major steps computing a priori and posteriori likelihoods for a set of models. The most posteriori likely model is used to select LOS-only measurements, which are then passed to the Kalman filter for processing.

By way of example, at epoch n the algorithm computes likelihoods for a set of models L_n(m), where m identifies a model. At the next epoch, there is a new set of models, which are not necessarily the same as the models in the old set. There is a probability that particular LOS signals become non-LOS and vice versa. Thus a probability of a flip in each flag can be defined and the probability transition matrix can be computed from the old set to the new one. If this probability transition matrix is applied to posteriori likelihoods in the old set, we will get a priori likelihoods in the new set, which we denote {tilde over (L)}_n+1(m)

Computation of posteriori likelihoods requires measurement redundancy. A covariance matrix for all measurements can be assigned. A basic assumption is that the error variance for a non-LOS measurement is large compared to that for LOS measurement (due to multipath-induced bias). The Bayesian probability density for having a particular set of measurements and some value of the user state (i.e., position, velocity and clock) can be computed as follows:

$\begin{array}{cc}P\left(\stackrel{\rightharpoonup }{x}|\stackrel{\_}{R}\right)=\mathrm{const}·\exp \left\{-\frac{1}{2}·{\left(\stackrel{\_}{x}-\stackrel{\rightharpoonup }{\tilde{x}}\left(\tilde{R}\right)\right)}^T·{\hat{Q}}^{-i}·\left(\tilde{x}-\stackrel{\rightharpoonup }{\tilde{x}}\left(\tilde{R}\right)\right)\right\}& \left[15\right]\end{array}$ where {circumflex over (Q)} is the covariance matrix, x( R) are expected values of measurements for user state R, and x are measurements. This probability density function is conditioned on the user state, which is unknown. The goal is to compute a Bayesian likelihood irrespective of the user state, which is the same as removing the condition in the above equation. This can be achieved by integrating the above equation over R. Integration can be performed in the closed form. The resulting equation defines the Bayesian likelihood of a particular set of measurements with a particular set of LOS/non-LOS flags (the latter affects the covariance matrix and the normalization constant). This likelihood is combined with the a priori likelihood of the model to arrive at a posteriori likelihood. As was stated above, the model with the highest likelihood “wins”, and measurements identified as LOS are passed to the Kalman filter for processing.

Conventional IMM implementation is not feasible for the applications discussed above due to a large number of models to track. To illustrate, if there are 10 signals, there are 1024 models (different allocations of LOS/non-LOS flags). It is clear that a workable solution must resort to tracking but a small fraction of all possible models. In GAME, a genetic algorithm (GA) is used to track only a subset of the possible models. There is some reasonably small set of models tracked by the algorithm concurrently. In each epoch, a model may remain in the set, may be eliminated, or may spawn more models to track. The fate of a model is dependent on its likelihood, which serves as a goal function in a typical genetic algorithm. Thus, most unlikely models face higher chances of elimination, and most likely ones have the offspring. A key aspect of integrating IMM and GA is the use of IMM-computed likelihood as the goal function of GAME.

As part of computing GAME likelihoods, a starting assumption is that measurements Z and their covariances {circumflex over (Q)} are available. Also available is the measurement matrix Ĥ. The unknowns include the position and clock vector X and the set of flags N for LOS/NLOS. The multipath signals are assumed to have additional bias in the measurements, and they may have a different variance as well. Additional bias is modeled as G* N, where G is a vector, and ‘*’ is a per-component vector product of two vectors. We also recognize that for multipath measurements, the variances should be set to a higher value. This results in matrix {circumflex over (Q)} being dependent on vector N. However, this will not affect the derivation below, and this dependence is not specified explicitly to keep equations less cumbersome.

The likelihood to have vector X and the set of flags N is computed as follows:

$\begin{array}{cc}{\Lambda }_0\left(\stackrel{\rightharpoonup }{X},\stackrel{\rightharpoonup }{N}\right)=C·\exp \left[-\frac{1}{2}{\left(\stackrel{\_}{Z}-\stackrel{\_}{H}\stackrel{\rightharpoonup }{X}-\tilde{G}*\tilde{N}\right)}^T·{\hat{Q}}^{-1}·\left(\stackrel{\_}{Z}-\hat{H}\stackrel{\_}{X}-\stackrel{\_}{G}*\stackrel{\_}{N}\right)\right].& \left[16\right]\end{array}$

In addition, prior information can be used in the likelihood computation, which comes from the output of the navigation filter. If the filter navigation solution is denoted by x and its covariance matrix by {circumflex over (P)}, then the equation for the likelihood takes the form:

$\begin{array}{cc}\Lambda \left(\stackrel{\_}{X},\stackrel{\_}{N}\right)=C·\exp \left[-\frac{1}{2}{\left(\stackrel{\_}{Z}-\hat{H}\tilde{X}-\tilde{G}*\stackrel{\_}{N}\right)}^T·{\hat{Q}}^{-1}·\left(\tilde{Z}-\hat{H}\stackrel{\rightharpoonup }{X}-\stackrel{\rightharpoonup }{G}*\tilde{N}\right)-\frac{1}{2}{\left(\stackrel{\rightharpoonup }{X}-\stackrel{\rightharpoonup }{x}\right)}^T·{\hat{P}}^{-1}·\left(\stackrel{\rightharpoonup }{X}-\stackrel{\rightharpoonup }{x}\right)\right]& \left[17\right]\end{array}$

However, to use it in updating model weights, just the likelihood of a vector of flags N is needed, irrespective of the user position and clock X. This means that the equation for Λ( X, N) is integrated with respect to X. To solve the integral, we transform the argument of the exponent to a more manageable form:

$\begin{array}{cc}E={\left(\stackrel{\_}{Z}-\hat{H}\tilde{X}-\stackrel{\rightharpoonup }{G}*\stackrel{\rightharpoonup }{N}\right)}^T·{\hat{Q}}^{-1}·\left(\tilde{Z}-\hat{H}\stackrel{\_}{X}-\stackrel{\rightharpoonup }{G}*\stackrel{\_}{N}\right)+{\left(\stackrel{\_}{X}-\stackrel{\_}{x}\right)}^T·{\hat{P}}^{-1}·\left(\stackrel{\rightharpoonup }{X}-\stackrel{\rightharpoonup }{x}\right)={\stackrel{\rightharpoonup }{X}}^T·\left({\hat{H}}^T·{\hat{Q}}^{-1}·\hat{H}+{\hat{P}}^{-1}\right)·\stackrel{\_}{X}+{\stackrel{\rightharpoonup }{X}}^T·{\hat{H}}^T·{\hat{Q}}^{-1}·\stackrel{\rightharpoonup }{\xi }+{\stackrel{\rightharpoonup }{\xi }}^T·{\hat{Q}}^{-1}·\hat{H}·\stackrel{\rightharpoonup }{X}+{\stackrel{\rightharpoonup }{\xi }}^T·{\hat{Q}}^{-1}·\stackrel{\_}{\xi }+{\stackrel{\_}{X}}^T·{\hat{P}}^{-1}·\stackrel{\rightharpoonup }{X}-{\stackrel{\_}{X}}^T·{\hat{P}}^{-1}·\stackrel{\rightharpoonup }{X}+{\tilde{x}}^T·{\hat{P}}^{-1}·\stackrel{\rightharpoonup }{x}& \left[18\right]\end{array}$ where ξ= G* N− Z is denoted for brevity.

The right hand side of the last equation is a quadratic form for X. It can be converted to the canonical form by performing a variable transformation:X={circumflex over (M)}· ρ  [19]The quadratic term is now as follows:X^T·(Ĥ^T·{circumflex over (Q)}⁻¹·Ĥ+{circumflex over (P)}⁻¹)· X= ρ^T·{circumflex over (M)}^T·(Ĥ^T·{circumflex over (Q)}⁻¹·Ĥ+{circumflex over (P)}⁻¹)·{circumflex over (M)}· ρ.  [20]This term is diagonalized (i.e., converted to the canonical form) if the variable transformation matrix is given by{circumflex over (M)}=Ĵ^T−1  [21]where Ĵ is a factor in the LU-decompositionĤ^T·{circumflex over (Q)}⁻¹·Ĥ+{circumflex over (P)}⁻¹=Ĵ·Ĵ^T.  [22]In the new variables, the argument of the exponent takes the following form:E= ρ^T· ρ+2·D^T· ρ+F=( ρ+D)²+F−D²  [23]where the following notations are used:D={circumflex over (M)}^T·(Ĥ^T·{circumflex over (Q)}⁻¹· ξ−{circumflex over (P)}⁻¹· x)F= ξ^T·{circumflex over (Q)}⁻¹· ξ+ x^T·{circumflex over (P)}⁻¹· x  [24]Substituting the expression for E into that for Λ( X, N) and integrate the result with respect to X:

$\begin{array}{cc}\Lambda \left(\stackrel{\rightharpoonup }{N}\right)=\int \int \int \Lambda \left(\stackrel{\rightharpoonup }{X},\stackrel{\rightharpoonup }{N}\right)·d^3\stackrel{\rightharpoonup }{X}=\mathrm{const}·\hat{M}·\exp \left[\frac{1}{2}\left(D^2-F\right)\right].& \left[25\right]\end{array}$ The constant multiplier and ∥{circumflex over (M)}∥ will be canceled when likelihood values for all models are normalized, and therefore do not have to be computed. This is the final (albeit still not normalized) equation for the likelihood of each model. It accounts for both the new measurements and for the prior information in the form of a navigation solution.

With Reference to FIG. 4, a simplified flow chart of the operation of the GAME method is illustrated. For the first epoch, an initial set of models is created 400. Models differ by LOS/non-LOS flags for different signals, each model having its own likelihood that is set to an a priori value. Next likelihoods and likelihood gradients in the LOS flag space can be computed 410. The current list of models is stored for future reference in memory (this list is denoted List_0). Based on likelihood gradients, likelihoods for new models can be computed, each differing from a current model by flipping one LOS flag. If new models have higher likelihoods than current ones, the new models are included in the set. This is repeated some number of times. Models are removed from the set if their likelihoods are lower than a threshold. The threshold is dynamically adjusted to keep some reasonable number of models in the set. New models are spawned and old ones pruned to produce an updated set of models, denoted as List_1 420. The most likely model is selected and its LOS flags outputted 430. Measurements corresponding to non-LOS signals will be discarded from processing by the navigation filter. All models are then propagated to the next epoch 440. All models in List_0 and all models in List_1 are then looped through 450.

Next the probabilities of transition from a model in List_0 to a model in List_1 is computed. The next step is to multiply the likelihood of the model in List_0 by the transition probability and tally it in the likelihood of the model in List_1 and proceed to the next time epoch.

One aspect the present disclosure is directed to a system and method which processes pseudorange, Doppler and ADR measurements from several navigation signal sources (such as GPS satellites), where some measurements can be substantially corrupted by the effects multipath or by other biases.

In another aspect, the present disclosure processes multiple models, where each model is characterized by assuming that some signals are non-biased (e.g., not affected by multipath, or affected insignificantly), and other signals are biased (e.g., corrupted by multipath delays). The number of models processed can be less than the total number of possible models. For example, if there are S sources, the total possible number of models, which assume all possible allocations between biased and non-biased measurements, is equal to 2^S. The number of processed models can be limited by the amount of processing power available to the user.

In another aspect, each processed model has a priori and posteriori likelihood values associated with it. The a priori likelihood computation assumes that each processed model at the current epoch has originated from a model, processed at the previous time epoch. The current model may be identical to one of the previous models, in which case its a priori likelihood depends on the probability of the signals not changing from biased to not-biased category and vice versa. Alternatively, the current model may differ from any model at the previous time epoch, in which case its a priori likelihood depends on the probability of a signal switching between the biased and non-biased category.

In yet another aspect, for each processed model, the algorithm computes the posteriori Bayesian likelihood of measurements obtained by the receiver. This likelihood is computed by combining the a priori likelihood with the probability of obtaining the last epoch's set of measurements. It will be different for different models, since a priori probability distributions for measurement residuals are different for biased and non-biased measurements. Typically, biased measurements would have larger variances, and the expected value of residuals of biased measurements may have non-zero mean.

For some or all of the processed models, the algorithm may also compute posteriori likelihoods of some derivatives of these models. For example, the algorithm may compute the likelihoods of models, where one signal source is moved from the biased to the non-biased category, or vice versa. From all candidate models (processed models and their derivatives), the algorithm selects models with highest likelihood values, which are retained. Models with lower likelihood values are destroyed. This selection process forms a set of models for computing a priori likelihoods for the next epoch. In other words, a priori likelihoods for the epoch are computed from posteriori likelihoods for the previous epoch. The model with the highest likelihood is selected. Measurements, which are flagged as non-biased for this model, are passed to the navigation filter for processing. Measurements which are deemed biased are not processed by the navigation filter.

WAFFLES

Another problem in prior art navigation systems is the errors introduced by Signals of Opportunity (SoOP) geolocation errors, and reference station calibration and timing errors. A priori locations of SoOPs may not be known accurately and therefore SoOPs must be geolocated using navigation infrastructure. Assuming that SoOPs are geolocated using the very same reference stations, which provide data to the user (possibly, even at the same time), SoOP geolocation errors are largely due to errors at the reference stations, and the latter are at the root of the problem. Thus, both types of errors must be mitigated jointly. We refer to these errors as CTG errors (for Calibration, Timing, and Geolocation).

Normally, any error in the source location, reference station calibration, or reference station timing translates into an error in user TDOA ranging. For example, if a source is hypothetically geolocated with 50 m accuracy, the same order of magnitude accuracy the the user navigation is observed, which could severely degrade navigation accuracy.

However, there is a small parameter to exploit. The errors in question are normally much smaller than other distances in the problem, e.g. those from the source to the user and from the source to reference stations. The ratio of the source geolocation error to the spatial scale of the problem is small, e.g. 50 m/5000 m˜10². The basic idea is to compute a weighted average of TDOA measurements in such a way that the effects of CTG errors cancel. If this goal is achieved, then the ranging error will be non-zero in the second-order approximation only. This may reduce the ranging error from ˜50 m to about 50×10⁻² m=0.5 m.

This method described in this disclosure is referred to as WAFFLES, which is designed to cancel effects of CTG errors in the first order approximation (with respect to the small parameter defined above). Even though this method may be an important piece of the puzzle to meet performance requirements, it should not be viewed as an excuse to relax efforts on calibrating reference stations and geolocating SoOPs accurately. Indeed, the performance of WAFFLES will depend on the geometry of the problem. Even though the first-order error is canceled, the second-order error may still be a concern for some geometries. This is when calibration and synchronization of reference stations becomes important. Moreover, the WAFFLES method allows an excellent return on any improvement in the calibration and timing. For example, reducing calibration error by a factor of 2 may reduce the second-order user ranging error by a factor of 4.

As compared to well known location of GPS satellites, SoOPs introduce a number of additional error sources. One of the most important error sources is due to timing and calibration errors at the reference station, and associated error in SoOP geolocation.

The geometry of the problem is shown in FIG. 5. SoOP location is not known precisely. If there is an error in the SoOP location estimate, it directly seeps into the error for TDOA and then into the user navigation error. For example, if the source is located at point 500 instead of point 510, the time for the signal to travel 515 from the estimated position 510 to reference station 520 is greater than the time for the signal to travel 525 from the actual source position 500 to the reference station 520 which introduces a timing error. Likewise, the time for the signal to travel 535 from the estimated position 510 to user 530 is less than the time for the signal to travel 545 from the actual source position 500 to the user 530 which introduces a timing error. The timing errors will then be introduced into the TDOA measurements. Below a way to cancel effects of reference station calibration is described.

This problem turns out to be somewhat complicated. The reason for this is the lack of the number of degrees of freedom available. If there are N reference stations, then there are N weights to play with. There are also N calibration errors to cancel, which is N constraints to satisfy. However, it is required that the sum of weights is equal to one, which adds one more constraint, for the total of N+1 constraints. Thus, all constraints can not be satisfied with N weights, and something has to be sacrificed.

As described below, it can be shown that the user ranging error due to reference station calibration errors is given by (in the linear approximation):δ ρ=K·Δ t  [26]This is the un-modeled error for the user for TDOA measurements using different reference stations. The error is due to reference station calibration errors Δ t and SoOP geolocation errors (which are in turn due to reference station calibration errors). If there are N reference stations, there are N components in the calibration error Δ t and correspondingly N TDOA measurements (each with respect to a different reference station) per SoOP.

This disclosure is directed to canceling these TDOA errors by computing a weighted sum of measurements at the cost of sacrificing two SoOP measurements. The approach is similar to GPS double differencing, but regular double differencing will not work here. In the case of GPS, double differencing cancels satellite clock error (which is a bit analogous to canceling reference station calibration errors in the present case) and satellite ephemeris errors (which are analogous to SoOP geolocation errors). However, GPS has the great advantage that a satellite is very far off, and all LOS are essentially parallel. If reference stations and the user are scattered across the area, like in our case, the plain old double differencing would not work. Hence the need for an approach like WAFFLES.

Additionally, weights can be computed for any epoch by the receiver, but they depend on the estimated user position. If the user moves around quite a bit, then weights will be time dependent. Since they are used as coefficients for terms, which define the partial derivatives of the measurement with respect to the user position, this is equivalent to changing the “virtual location” of the equivalent TOA SoOP measurement. This will not be a problem for a navigation code, which is designed from scratch properly.

By way of illustration, consider a 3D case with 4 reference stations (a 2D case with 3 reference stations is completely analogous). The pseudorange equation for the j-th reference station:

$\begin{array}{cc}{\rho }_j={\stackrel{\rightharpoonup }{R}}_j-\stackrel{\rightharpoonup }{S}+\theta +\Delta t_j& \left[27\right]\end{array}$ where R_j is the position of the reference station, S is position of the source, θ is the true timing of some feature in the SoOP signal, and Δt_j is the calibration error.Assuming some initial approximation S₀ for the SoOP location, the estimated position Ŝ of the SoOP is computed as a correction to this initial approximation:Ŝ= S₀+δ S.  [28]The unknowns are δ S, θ, which can be combined into one vector δx={δ S,θ}. Then in the first order the solution for δx is obtained from a linear system of equations:

$\begin{array}{cc}\delta \stackrel{\_}{x}=N^{-1}·\left(\rho -\stackrel{\_}{R}-{\stackrel{\_}{S}}_0-\Delta t\right)& \left[29\right]\end{array}$ where components of 4-vector ρ−| R− S₀|−Δt correspond to different reference stations (i.e., j-th component is ρ_j−| R_j− S₀|−Δt_j) and where the 4×4 matrix N is defined by the following:

$\begin{array}{cc}{N}_{jk}=\{\begin{array}{ccc}\frac{\partial \stackrel{\_}{S}-{\stackrel{\_}{R}}_j}{\partial S_k}& \mathrm{if}& k\le 3\\ 1& \mathrm{if}& k=4\end{array}& \left[30\right]\end{array}$ The solution for δ x has two additive parts: the exact correction, and the error due to reference station calibration. They can be separated:δx=δ x₀+ ε  [31]where

$\begin{array}{cc}\delta {\stackrel{\_}{x}}_0=N^{-1}·\left(\rho -\stackrel{\_}{R}-{\stackrel{\_}{S}}_0\right)\tilde{\varepsilon }=N-1·\Delta t& \left[32\right]\end{array}$ Vector ε is the SoOP geolocation error. As we shall see shortly, it affects the user position error.

Proceeding to estimating the user position, the user pseudorange equation is as follows:

$\begin{array}{cc}{\rho }_u=\stackrel{\_}{Y}-\stackrel{\_}{S}+\theta +\Delta t_u& \left[33\right]\end{array}$ where Y is the true user position, and Δt_u is the user clock error. The true TDOA with the reference station pseudorange is as follows:

$\begin{array}{cc}{\rho }_u-{\rho }_j=\stackrel{\_}{Y}-\stackrel{\_}{S}+\Delta t_u-{\stackrel{\rightharpoonup }{R}}_j-\stackrel{\rightharpoonup }{S}-\Delta t_i.& \left[34\right]\end{array}$ To process this measurement, it is compared with the estimate for TDOA, i.e. compute the residual. This residual will comprise two parts: (a) The “legitimate part”, which is due to the user clock and user position error. These are unknowns we solve for, so we need that error in the residual to have something to process; and (b) The “error part”, which we do not model and do not solve for. This part is due to term Δt_j and due to SoOP geolocation error s, which is ultimately due to Δt_j as well.

The present disclosure is directed to the second part. To make the math less cumbersome, only the second error source is retained and the residual due to the user position and clock errors is not computed:δ ρ=M ε−Δ t  [35]where 4×4 matrix M is defined by:

$\begin{array}{cc}{M}_{jk}=\{\begin{array}{ccc}\frac{\partial \stackrel{\_}{Y}-\stackrel{\_}{S}-{\stackrel{\_}{R}}_j-\stackrel{\_}{S}}{\partial S_k}& \mathrm{if}& k\le 3\\ 0& \mathrm{if}& k=4\end{array}& \left[36\right]\end{array}$ Substituting ε gets:δ ρ=(MN⁻¹−1)·Δt.  [37]Next, we introduce the notationK=MN⁻¹−1  [38]to get:δ ρ=K·Δ t  [39]

Unfortunately, computation shows that matrix K does not have a full rank. Computing it explicitly for a 2D case is somewhat involved, but important. Assuming there are three reference stations.δS·Ū_j=Δ t_j+δθ  [40]where

$\begin{array}{cc}{\overrightarrow{U}}_l=\frac{\overrightarrow{S}-{\overrightarrow{R}}_l}{\overrightarrow{S}-{\overrightarrow{R}}_l}& \left[41\right]\end{array}$ is the unit vector in the SoOP to reference station direction.Differencing equations 0 and 1, and 1 and 2 in [40] to cancel δθ and solve for δ S to get:

$\begin{array}{cc}\delta {\overrightarrow{S}}_k=\sum _j{\alpha }_{kj}\Delta {\overrightarrow{t}}_j& \left[42\right]\end{array}$ where subscript k is for x,y components and

$\begin{array}{cc}{a}_{x0}=\frac{U_{21y}}{D};a_{x1}=\frac{U_{02y}}{D};a_{\mathrm{x2}}=-\frac{U_{01y}}{D};a_{y0}=\frac{U_{12x}}{D};a_{y1}=-\frac{U_{02x}}{D};a_{y2}=\frac{U_{01x}}{D}& \left[43\right]\end{array}$ and

$\begin{array}{cc}{U}_{\mathrm{pqk}}=U_{pk}-U_{\mathrm{qk}}D=U_{01y}{U}_{02x}-U_{02y}{U}_{01x}& \left[44\right]\end{array}$ User TDOA error is given by (see (8)):δρ_j=ξ_jxS_x+ξ_jyS_y+Δt_j  [45]where

$\begin{array}{cc}{\xi }_{\mathrm{jk}}=\frac{S_k-Y_k}{\overrightarrow{S}-\overrightarrow{Y}}-\frac{S-R_{\mathrm{jx}}}{\overrightarrow{S}-{\overrightarrow{R}}_j}& \left[46\right]\end{array}$ Thus, matrix K is given by:K_jm=ξ_jxα_xm+ξ_jyα_ym+δ_jm.  [47]Next we compute K_jm explicitly for some specific m, for example for m=0 and for different values of j.Case m=0; j=1

$\begin{array}{cc}{\xi }_{1x}{a}_{x0}+{\xi }_{1y}{a}_{y0}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x1}·\left(U_{y2}-U_{y1}\right)+U_{y1}·\left(U_{x1}-U_{x2}\right)\right\}}{D}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x1}·U_{y2}-U_{y1}·U_{x2}\right\}}{D}& \left[48\right]\end{array}$ Case m=0; j=2

$\begin{array}{cc}{\xi }_{1x}{a}_{x0}+{\xi }_{1y}{a}_{y0}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x2}·\left(U_{y2}-U_{y1}\right)+U_{y2}·\left(U_{x1}-U_{x2}\right)\right\}}{D}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x1}·U_{y2}-U_{y1}·U_{x2}\right\}}{D}& \left[49\right]\end{array}$ Case m=0; j=0

$\begin{array}{cc}{\xi }_{1x}{a}_{x0}+{\xi }_{1y}{a}_{y0}+1=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x0}·\left(U_{y2}-U_{y1}\right)+U_{y0}·\left(U_{x1}-U_{x2}\right)\right\}}{D}+1=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{U_{x0}·U_{y2}-U_{x0}·U_{y1}+U_{y0}·U_{x1}-U_{y0}·U_{x2}}{D}+\frac{U_{01y}{U}_{02x}-U_{02y}{U}_{01x}}{D}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{U_{x0}·U_{y2}-U_{x0}·U_{y1}+U_{y0}·U_{x1}-U_{y0}·U_{x2}}{D}+\frac{\begin{array}{c}{U}_{y0}{U}_{x0}-U_{y0}{U}_{x2}-U_{y1}{U}_{x0}+U_{y1}{U}_{x2}-\\ U_{y0}{U}_{x0}+U_{y0}{U}_{x1}+U_{y2}{U}_{x0}-U_{y2}{U}_{x1}\end{array}}{D}=\frac{S_x-Y_x}{\overrightarrow{S}-\overrightarrow{Y}}·a_{x0}+\frac{S_y-Y_y}{\overrightarrow{S}-\overrightarrow{Y}}·a_{y0}-\frac{\left\{U_{x1}·U_{y2}-U_{y1}·U_{x2}\right\}}{D}.& \left[50\right]\end{array}$ Amazingly, K₀₀=K₁₀=K₂₀. The same holds true for other values of m.

This means that errors in a particular reference station contribute the same value to the user TDOA error regardless of whether the same or different reference station is used for making that TDOA. If it is the same station, then this contribution is via SoOP geolocation error and the reference station calibration error. If a different station is used, the TDOA user error is due to SoOP geolocation error only. Yet, in both cases, the net result is the same.

This complicates the algorithm, forcing a more costly route. First, there is no absolute concept of time (as in GPS), and therefore timing and calibration errors at one station can be set to zero. This station will be treated as an absolute time reference.

In this case, there are 2 independent non-zero calibration errors for 3 stations, e.g.Δt₀=0; Δt₁,Δt₂≠0.Considering 3 SoOPs and introducing another subscript to equations to denote a particular SoOP. In addition, only TDOA measurements from one reference station (recall that errors from all stations are equal, so any station will suffice) are considered, which is denoted with as m=0. User TDOA error (again, with any reference station, e.g. with station 0) is given by:

$\begin{array}{cc}{\mathrm{\delta \rho }}_s=\sum _j{K}_{s0j}·\Delta t_j& \left[51\right]\end{array}$ where subscript s denotes different SoOPs, and not different reference stations as in (13). Setting Δt₀=0, we combine measurements for three SoOPs with weights α_s, such that coefficients for Δt₁, Δt₂ cancel, and the sum of weights is equal to 1 (so far). This produces the following system of equations for α_s:K₀₀₁α₀+K₁₀₁α₁+K₂₀₁α₂=0K₀₀₂α₀+K₁₀₂α₁+K₂₀₂α₂=0  [52]α₀+α₁+α₂=1If TDOA measurements from 3 SoOPs are linearly combined with weights α_s, the first order error cancels.

The WAFFLES method assumes some typical reference station position errors and some geometry (i.e., user, SoOP, and reference stations positions), computes the true and the estimated value of TDOAs for 3 SoOPs, computes coefficients α_s, and then combines measurements. If the algorithm works, the difference between the two weighted averages must be relatively small.

Implementation of the algorithm has two subtleties, which are described below:

• 1. Geolocation is done in 2D. Vertical positions of SoOPs are assumed to be known with some accuracy. The truth measurements and the estimated measurements are simulated using the true and the estimated vertical positions of SoOPs respectively.
• 2. User position is not known perfectly either. If the user were located exactly at the estimated position, the net results of using TDOA measurements with WAFFLES algorithm would be the measurement residual in the amount of the error, not compensated by WAFFLES. If the user is off its estimated position (i.e., if there is a position error estimate by the navigation filter), then there is an additional error source. This error is due to an error in the LOS direction to the SoOP (LOS is computed the estimated SoOP position instead of the true one). The magnitude of this error is the difference of the projection of the user position error on the estimated LOS and the true LOS. Assuming that the vector of the user position error can point at any direction equally probable means that this additional error can be negative, positive, or zero. The average magnitude of it is 1/√{square root over (2)} of the maximum error. The same is true for the SoOP geolocation error. Two factors of 1/√{square root over (2)} produce a factor of 1/2, yielding the following:

$\begin{array}{cc}\varepsilon =\frac{\delta \overrightarrow{U}·\delta \overrightarrow{S}}{2\overrightarrow{U}-\overrightarrow{S}}& \left[53\right]\end{array}$ where δŪ is the user position error, δ S is the SoOP geolocation error, and Ū− S is the range to the SoOP.There is just one more little modification to the equations to be made. The weighted sum of TDOA measurements can be used in the Kalman filter for navigating the user. Normally, a TOA or TDOA measurement equation has the position of the user in the form like |Ū− S|+const. This means that partial derivatives of the measurement equation with respect to the user position form a unit vector in the direction of the source. Thus, the RSS (root sum squared) of the partials is equal to one. Note that α₀+α₁+α₂=1 for weighted TDOA does not guarantee this condition. In order to imitate this condition for the new weighted TDOA measurement, it should be scaled. (In general, this is not necessary, but keeps the lineage to the conventional TOA/TDOA processing; scaling becomes necessary if the Kalman filter formulation is not changed from the conventional TOA/TDOA processing.) Thus, partial derivatives can be computed with respect to the user position and then the scaling can be computed:

$\begin{array}{cc}{w}^2={\left(\sum _s{\alpha }_s\frac{\overrightarrow{U}-{\overrightarrow{S}}_{\mathrm{sx}}}{\overrightarrow{U}-\overrightarrow{S}}\right)}^2={\left(\sum _s{\alpha }_s\frac{\overrightarrow{U}-{\overrightarrow{S}}_{\mathrm{sy}}}{\overrightarrow{U}-\overrightarrow{S}}\right)}^2& \left[54\right]\end{array}$ Now the desired weighted TDOA measurement is given by:

$\begin{array}{cc}\tilde{\rho }=\frac{1}{w}\sum _s{\alpha }_s{\rho }_s& \left[55\right]\end{array}$ where ρ_s are TDOA measurements for SoOP s using reference station 0.

FIG. 6 illustrates a simplified block diagram of the WAFFLES method. Reference stations measures timing features of the transmitted signals 600. The transmitted signals can be satellite signals or terrestrial signals. The reference stations transmit this information and source geolocation to the user via a communications channel 610. Locations of the reference stations are known with a high accuracy (e.g., positions are surveyed). Locations of the signal transmitters is generally known approximately, with some accuracy (e.g., from design drawings for a cell tower). Multiple reference stations are used to geolocate each signal transmitter. The effects of calibration, timing, and geolocation (CTG) errors at the reference stations are mitigated by the user 620. The user makes pseudorange and Doppler measurements 630. Navigation equations for the user are formulated to explicitly express dependence of pseudorange on the CTG errors. Navigations equations are linearized with respect to CTG errors. Navigations equations are linearly combined to cancel the linear component of CTG-induced pseudorange errors. The linear combination measurements are output to a navigation filter 640.

Estimating Location of a Static User Under Foliage

Another aspect of the present disclosure provides a method of estimating user location under foliage or for other environments dominated by scattering of GPS signals. In one embodiment, the statistical properties of the foliage are assumed to be known with some reasonable accuracy. It is not necessary to know the positions of individual twigs and leaves, rather, correlations for the channel impulse response can be used.

In one embodiment, the user receiver has a bank of correlators. The receiver can compute correlations between the replica waveform and the signal waveform at different delays. Thus, for each SoOP there will be a large number of correlation measurements, taken at different time epochs and at different delays. Measurements, which are sufficiently separated in time and/or delay will not be correlated; however, measurements taken at adjacent time epochs or by adjacent correlators will be correlated. Correlation is due to two separate effects. First, the scattering medium (foliage) is the same for different measurements; thus there will be signal components coming from same branches, twigs and leaves. Second, the signal bandwidth is obviously limited; this creates correlation of measurements over delays, which are separated by less than the inverse bandwidth of the signal.

The channel impulse response (CIR) for a particular SoOP is a function of delay τ and time t. The time dependence is due to the satellite motion and the change in the environment (e.g., wind moves the leaves):R=R(τ,t)  [56]The CIR has the LOS peak at τ=0, which corresponds to the true user position.The correlation of this signal with the replica waveform is measured. The correlation as a function of delay and time is given by convolution of CIR with the autocorrelation function of the waveform:C(τ,t)=R(τ,t)*A(τ)=∫R(θ,t)·A(τ−θ)dθ.  [57]Assuming that statistical properties of R(τ,t) are known, pair-wise covariances for values of C(τ,t) can be computed:

$\begin{array}{cc}Q\left({\tau }_i,t_i,{\tau }_j,t_j\right)=E\left\{C\left({\tau }_i,t_i\right)·C^{*}\left({\tau }_j,t_j\right)\right\}=\int \int R\left({\theta }_1,t_i\right)·A\left({\tau }_i-{\theta }_1\right)R^{*}\left({\theta }_2,t_j\right)·A\left({\tau }_j-{\theta }_2\right)d{\theta }_{1}d{\theta }_2& \left[58\right]\end{array}$ Similarly, a noise covariance matrix can be computedη(τ_i,t_i,τ_j,t_j)=E{n(τ_i,t_i)·n*(τ_j,t_j)}.  [59]This computation can be done once for each type of forest and stored at the receiver in the form of a lookup function for different pairs of angles of arrival and delays. (Even though arguments for this covariance matrix include time, and not the angle of arrival per se, the statistics is mostly due to the change in the angle as the time progresses and satellite SoOP moves. From a pre-stored table, the receiver can extract the covariance matrices for a particular satellite pass using the specific geometry.)

The problem can be formulated as follows. Correlations between the signal and the waveform replica can be measured. These measurements comprise a realization {tilde over (C)} of a random process. The receiver position can be estimated from this data.

An estimate for the overall offset of the delay can be denoted with a hat, {circumflex over (τ)}. At a particular time, {circumflex over (τ)}=τ−Δ r· l(t), where Δ r is the error in the user position estimate, and l(t_i) is the LOS direction at time t_i. When the correlations {tilde over (C)} are measured, the true value of the replica delay τ is unknown and therefore the estimate delay {circumflex over (τ)} is used. Thus, measured correlations are in the form {tilde over (C)}({circumflex over (τ)},t). The main concept is that if the position error is zero, then {circumflex over (τ)}=τ, and the measured realization of the random process is most likely in the statistical sense. Thus, the goal is to find such {circumflex over (τ)}, that equation{tilde over (C)}({circumflex over (τ)},t)=C(τ,t)+n(τ,t)  [60]is the most likely realization of the random process for C(τ,t) and n(τ,t).For discrete values of τ_i,t_i, equation [60] can be viewed as an over-determined system of linear equations for C(τ_i,t_i),n(τ_i,t_i). If one long vector of unknowns is formed x={C(τ_i,t_i)|n(τ_i,t_i)} then equation [60] can be written in the form:Mx={tilde over (C)}({circumflex over (τ)},t)  [61]where matrix M has the following structure:

$\begin{array}{cc}M=\left(\begin{array}{ccc}1& & 0\\ & \ddots & \\ 0& & 1\end{array}|\begin{array}{ccc}1& & 0\\ & \ddots & \\ 0& & 1\end{array}\right)& \left[62\right]\end{array}$ A solution of equation [61] is sought while minimizing the following quadratic form (this is what makes the solution the most likely one):F=xWx*  [63]where the covariance matrix W is constructed from the covariance matrices for the signal and for the noise:

$\begin{array}{cc}W=\left(\begin{array}{cc}Q& 0\\ 0& \eta \end{array}\right)& \left[64\right]\end{array}$

This is called a linear equality-constrained least squares problem, and there are efficient numerical routines to solve it. Since {circumflex over (τ)} depends on the user position error, the solution will too be a function of the position error (via equation [61]). Thus, the value of the quadratic form [63] will be a function of the user position error. By varying Δ r, the absolute minimum of F can be found. The minimum value of F will correspond to the position error, which makes measurements the most likely realization of the random process, i.e. the best estimate in the statistical sense.

FIG. 7 illustrates LEAF processing. From as assumed initial pseudorange 700, pseudorange, Doppler and ADR measurements from several navigation signal sources (such as GPS satellites)are collected 710, where measurements are corrupted by diffuse scattering of the signal, e.g. from foliage. Multiple correlators may be used, which estimate correlations of the signal with the signal replica. The correlators may differ by the delays of the signal replica. The measured correlation values are assumed to be generated by the noise and (at least for some correlators) by the direct and scattered signal. Using the statistical properties of noise and the scattered signal (which are preset in the receiver for typical environments), the present disclosure computes the most likely partition of measured values into noise and signal 720. The Bayesian likelihood associated with this partition is calculated 730. By varying the assumed value of the pseudorange 740, the algorithm maximizes the Bayesian likelihood of the noise/signal partition. The pseudorange value which maximizes the likelihood, is the estimated pseudorange 750.

It may be emphasized that the above-described embodiments, particularly any “preferred” embodiments, are merely possible examples of implementations, merely set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiments of the disclosure without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and the present disclosure and protected by the following claims Embodiments of the subject matter and the functional operations described in this specification can be implemented in digital electronic circuitry, or in computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Embodiments of the subject matter described in this specification can be implemented as one or more computer program products, i.e., one or more modules of computer program instructions encoded on a tangible program carrier for execution by, or to control the operation of, data processing apparatus. The tangible program carrier can be a propagated signal or a computer readable medium. The propagated signal is an artificially generated signal, e.g., a machine-generated electrical, optical, or electromagnetic signal that is generated to encode information for transmission to suitable receiver apparatus for execution by a computer. The computer readable medium can be a machine-readable storage device, a machine-readable storage substrate, a memory device, a composition of matter affecting a machine-readable propagated signal, or a combination of one or more of them.

The term “processors or processing” encompasses all apparatus, devices, and machines for processing data, including by way of example a programmable processor, a computer, or multiple processors or computers. The processor can include, in addition to hardware, code that creates an execution environment for the computer program in question, e.g., code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them.

A computer program (also known as a program, software, software application, script, or code) can be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and it can be deployed in any form, including as a stand alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document), in a single file dedicated to the program in question, or in multiple coordinated files (e.g., files that store one or more modules, sub programs, or portions of code). A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network.

The processes and logic flows described in this specification can be performed by one or more programmable processors executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by, and apparatus can also be implemented as, special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application specific integrated circuit).

Processors suitable for the execution of a computer program include, by way of example, both general and special purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor will receive instructions and data from a read only memory or a random access memory or both. The essential elements of a computer are a processor for performing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to receive data from or transfer data to, or both, one or more mass storage devices for storing data, e.g., magnetic, magneto optical disks, or optical disks. However, a computer need not have such devices. Moreover, a computer can be embedded in another device, e.g., a mobile telephone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a Global Positioning System (GPS) receiver, to name just a few.

Computer readable media suitable for storing computer program instructions and data include all forms of non volatile memory, media and memory devices, including by way of example semiconductor memory devices, e.g., EPROM, EEPROM, and flash memory devices; magnetic disks, e.g., internal hard disks or removable disks; magneto optical disks; and CD ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry.

To provide for interaction with a user, embodiments of the subject matter described in this specification can be implemented on a computer having a display device, e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor, for displaying information to the user and a keyboard and a pointing device, e.g., a mouse or a trackball, by which the user can provide input to the computer. Other kinds of devices can be used to provide for interaction with a user as well; for example, input from the user can be received in any form, including acoustic, speech, or tactile input.

Embodiments of the subject matter described in this specification can be implemented in a computing system that includes a back end component, e.g., as a data server, or that includes a middleware component, e.g., an application server, or that includes a front end component, e.g., a client computer having a graphical user interface or a Web browser through which a user can interact with an implementation of the subject matter described is this specification, or any combination of one or more such back end, middleware, or front end components. The components of the system can be interconnected by any form or medium of digital data communication, e.g., a communication network. Examples of communication networks include a local area network (“LAN”) and a wide area network (“WAN”), e.g., the Internet.

The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other.

While this specification contains many specifics, these should not be construed as limitations on the scope of any invention or of what may be claimed, but rather as descriptions of features that may be specific to particular embodiments of particular inventions. Certain features that are described in this specification in the context of separate embodiments can also be implemented in combination in a single embodiment. Conversely, various features that are described in the context of a single embodiment can also be implemented in multiple embodiments separately or in any suitable subcombination. Moreover, although features may be described above as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a subcombination or variation of a subcombination.

Similarly, while operations are depicted in the drawings in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed, to achieve desirable results. In certain circumstances, multitasking and parallel processing may be advantageous. Moreover, the separation of various system components in the embodiments described above should not be understood as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

It may be emphasized that the above-described embodiments, particularly any “preferred” embodiments, are merely possible examples of implementations, merely set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiments of the disclosure without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and the present disclosure and protected by the following claims.

Claims

1. A method to increase accuracy of geolocation of a user device based on timing measurements of RF signals, the method comprising: (a) receiving, at a user device, one or more RF signals transmitted over time from a transmitter;(b) correlating, at the user device, the received signals with a plurality of replica waveforms;(c) creating, at the user device, a first set of estimates of times of arrival of the received signals at the user device;(d) using the first set of estimates to create a plurality of models at the user device, wherein each model assumes at least one of the first set of estimates is affected by multipath or interference;(e) processing, at the user device, the plurality of models to estimate the likelihood of times of arrival for each model;(f) determining, at the user device, a set of line of sight (LOS) measurements for the user device using a model having the largest likelihood;(g) determining, at the user device, a location of the user device as a function of the determined set of LOS measurements.

2. The method of claim 1 wherein creating the plurality of models and estimating their likelihoods uses an Interacting Multiple Models (IMM) method.

3. The method of claim 1 wherein creating the plurality of models uses a Genetic Algorithms method.

4. The method of claim 1 wherein the processing includes using Bayesian Likelihood.

5. The method of claim 1 wherein the step of creating the plurality of models includes the use of likelihood gradients to estimate the likelihoods of the plurality of models, which may be created at the next time interval.

6. The method of claim 1 wherein the processing includes removing a model if its likelihood is lower than a threshold likelihood.

7. The method of claim 6 wherein the threshold likelihood is dynamically adjusted to keep some reasonable number of models.

8. The method of claim 1 wherein the processing includes calculating a priori likelihood for the plurality of models.

9. The method of claim 1 wherein the processing includes calculating posteriori likelihood for the plurality of models.

10. The method of claim 1 wherein the step of creating a first set of estimates includes estimating amplitude and Doppler of the received signals at the user device, and the processing the plurality of models includes testing the plurality of models to estimate the likelihood of amplitude and Doppler measurements for each model.

11. A method of locating a user device having a signal receiver in an environment having a plurality of transmitters, and a plurality of reference stations of known location for receiving signals transmitted by the transmitters, the method comprising: (a) receiving, at a plurality of reference stations, a first signal transmitted by a first one of a plurality of transmitters;(b) estimating, at a reference station in the plurality of reference stations, the location of the first transmitter using time difference of arrival of the first signal at pairs of the plurality of reference stations;(c) transmitting, from the reference station to the user device, signal timing information and estimated pseudo ranges from the first transmitter to the plurality of reference stations;(d) forming, at the user device, a linear combination of multiple time differences of arrival, each time difference of arrival comprising a time of receipt of the first signal at the user device and at a particular reference station of the plurality of reference stations;(e) processing, at the user device, the linear combination of time differences of arrival to compute an estimate of a location of the user device; and(f) repeating steps (a)-(e) for each of the plurality of transmitters.

12. The method of claim 11 further comprising forming navigation equations that express dependence of the estimated pseudo ranges on timing errors at the plurality of reference stations and on positional errors of the plurality of transmitters.

13. The method of claim 12 wherein the step of forming linear combinations cancels the timing and positional errors.

14. A method of estimating a location of a user device in a propagation environment having a plurality of transmitters (a) receiving an RF signal at a user device where the RF signal is transmitted over time from a transmitter;(b) correlating, at the user device, the received RF signal with a plurality of delayed replica waveforms;(c) testing, at the user device, a plurality of estimated distances from the user device to the transmitter;(d) accessing, at the user device, a database of historical statistics for a plurality of correlation values;(e) comparing, at the user device, calculated statistics for each of the plurality of estimated distances with the accessed historical statistics; and(f) estimating, at the user device, a pseudo range between the user device and the transmitter as a function of the comparison.

15. The method of claim 14 wherein the step of comparing includes characterizing the received RF signal as being comprised of noise and a transmitter-generated signal, the noise having propagation delay.

16. The method of claim 14 wherein the steps of (a)-(f) are repeated for the plurality of transmitters, and the location of the user device is determined as a function of the estimated pseudo ranges between the user device and the plurality of transmitters.

17. The method of claim 14 wherein the step of correlating uses a bank of correlators.

18. The method of claim 14 wherein the each of the plurality of replica waveforms is an expected received signal delayed in time.

19. The method of claim 18 wherein the expected received signal is based on an estimated location of the transmitter.

20. The method of claim 14 wherein the location of the user device is estimated based on an inertial navigation system.

21. The method of claim 15, wherein the step of testing an estimated distance includes determining a Bayesian likelihood of having a portion of the received RF signal between the noise and the transmitter-generated signal, the noise having the delay associated with the estimated distance.