# Patent

## Title
Positioning determinations of receivers

## Abstract

The disclosed method and system is used to determine the position of a user device. The user device can receive data signals and/or carrier signals from orbiting space vehicles. These data signals can be used for positioning calculation and/or track maintenance of the user device. The disclosed method and system can account for time and frequency biases of the user device. For the track maintenance, a Kalman filter state estimator can be extended to include a velocity of the user device.

## Background
RELATED APPLICATION
This application claims the benefit of and priority to U.S. Provisional Application Ser. No. 61/171,009, filed Apr. 20, 2009, and U.S. Provisional Application Ser. No. 61/177,579, filed May 12, 2009, the contents of both of which are incorporated by reference herein in their entirety.
U.S. Pat. No. 7,489,926 issued to Whelan, et al. is hereby incorporated by reference in its entirety in order to more fully describe the state of the art as known to those skilled therein as of the date of the invention described herein.


FIELD
The present disclosure relates to satellite navigation and estimating the position of fixed or moving user devices.
BACKGROUND
Global Positioning Systems (GPS) are widely used to determine the position of a receiver on the earth. The receiver can use position and time data of GPS satellites to calculate its position. Sometimes, however, a GPS receiver or signal is not available, and an alternative method of satellite geolocation may be desirable. The unavailability of the GPS signal may be due several factors. These factors may include, but are not limited to, diminished signal strength in an occluded environment (e.g., buildings, heavy foliage, etc.), deliberate or incidental radio frequency interference, or hardware or software malfunctions resulting in degraded receiver performance.
SUMMARY
The present disclosure relates to a system, apparatus, and method of determining the position of a user device, such as an Iridium receiver, which contains a time bias and a frequency bias. The present disclosure also relates to a system, apparatus, and method of estimating the position and velocity of a user device in addition to estimating the time and frequency bias of that same device.
Some embodiments of the present disclosure provide a method for determining the position of a user device. The user device can receive signals from space vehicles. Based on correlation time and frequency deviation measurements derived from the received signals, the user device can calculate the current position. The signals may not include position information of the space vehicles in order for the user device to determine the current position. A short correlation time of the signals can be sufficient for accurate results.
One or more embodiments of the disclosure provide a method for track maintenance. Track maintenance is the process of estimating a position of a user device if data for positioning calculations becomes unavailable. The method can account for a time bias and a frequency bias. The method can employ an extended Kalman filter state estimator, which can include a velocity of the user device. In some embodiments, the velocity of the user device, time and frequency biases, and statistical errors can be accounted for using only information of latitude and longitude.
Some embodiments of the disclosure provide a user device capable of receiving signals from at least one space vehicle. The user device can sample the received signals in order to compute a correlation time and frequency deviation measurement. Based on these measurements, the user device can determine its position. If the user device does not receive signals from at least one space vehicle, the user device can estimate a position using an extended Kalman filter state estimator, which can include a velocity of the user device.
In one or more embodiments, the method for track maintenance of a user device includes receiving a signal from at least one space vehicle, and obtaining a position of the user device based on the received signal. In some embodiments, the method further includes applying an extended Kalman filter state estimator including the position of the user device and a velocity of the user device. In at least one embodiment, the method also includes obtaining an estimate of a state of the user device for the track maintenance.
In some embodiments, the obtaining of the position of the user device includes a snapshot calculation. In one or more embodiments, the performing of the snapshot calculation and/or the extended Kalman filter state estimator includes a pseudo range measurement model and a pseudo range rate measurement model. (The pseudo range rate is a function of the sum of the Doppler and the frequency deviations of the transmitters and/or the receivers.) In at least one embodiment, the computing of the pseudo range measurement model and the pseudo range rate measurement model includes a nonlinear least squares method.
In one or more embodiments, the performing of the snapshot calculation and/or the extended Kalman filter state estimator includes a time bias and a frequency bias. In some embodiments, the computing of the time bias and the frequency bias includes a linear least squares method. In at least one embodiment, the computing of the time bias includes approximating the time bias with a polynomial. In some embodiments, the computing of the frequency bias includes taking the derivative of the time bias with respect to time. In one or more embodiments, the time bias and the frequency bias are expressed as explicit functions of a position of the user device. In at least one embodiment, expressions of the time bias and the frequency bias modify a nonlinear least squares cost function by removing the biases as independent states, thereby reducing the number of independent variables to those of a position of the user device alone.
In some embodiments, the signal is received from a low earth orbit (LEO) satellite. In one or more embodiments, the LEO satellite is part of the Iridium system. In at least one embodiment, the state of the user device includes a position of the user device, a velocity of the user device, a time bias, and a frequency bias. In one or more embodiments, obtaining the estimate of the state of the user device includes evaluating a latitude and a longitude of the position of the user device.
In one or more embodiments, the method for positioning calculations in an occluded environment includes receiving signals from at least one space vehicle, and computing a pseudo range and a pseudo range rate including a time bias and a frequency bias of the signals. In some embodiments, the method further comprises minimizing a cost function including the pseudo range and the pseudo range rate to obtain a position. In at least one embodiment, the cost function includes a product of the pseudo range and the pseudo range rate. In some embodiments, the computing of the pseudo range and the pseudo range rate includes a nonlinear least squares method.
In some embodiments, the computing of the time bias and the frequency bias includes a linear least squares method. In one or more embodiments, the time bias and the frequency bias are expressed as explicit functions of a position of the user device. In at least one embodiment, expressions of the time bias and the frequency bias modify a cost function by removing the biases as independent states, thereby reducing the number of independent variables to those of a position of the user device alone.
In one or more embodiments, the signal is received from a LEO satellite. In some embodiments, the LEO satellite is part of the Iridium system. In at least one embodiment, the occluded environment includes an indoor environment, a jamming environment, and/or other environment where the signal is degraded by incidental or deliberate radio frequency interference. In some embodiments, the obtaining of the position of the user device is based on latitude and longitude information.
In one or more embodiments, the user device includes an antenna capable of receiving signals from at least one space vehicle. In some embodiments, the user device further includes a computer system capable of processing the received signals to calculate a time bias and a frequency bias, the computer system using the time bias and the frequency bias to determine a position of the user device. In at least one embodiment, the computer system can execute an extended Kalman filter state estimator including a velocity of the user device in order to estimate a position of the user device if the antenna does not receive the signals of the at least one space vehicle. In some embodiments, the computer system further comprises an output device.



DRAWINGS
These and other features, aspects, and advantages of the present disclosure will become better understood with regard to the following description, appended claims, and accompanying drawings where:
FIG. 1 is a prior art block diagram of an indoor positioning system using LEO satellites.
FIG. 2 is a prior art diagram of a differential positioning system using LEO satellites.
FIG. 3a is a prior art graphic representation of a system covariance from LEO and MEMS sources after a first pass.
FIG. 3b is a prior art graphic representation of a system covariance from LEO and MEMS sources after subsequent passes.
FIG. 4 is a prior art block diagram of a tightly coupled LEO inertial integrator.
FIG. 5 is a prior art flowchart to describe a process for deriving a position based upon a LEO signal and an inertial position fix.
FIG. 6 is a schematic illustration of a plurality of space vehicles (SVs) orbiting the earth, in accordance with at least one embodiment of the present disclosure.
FIG. 7 is a schematic illustration of how a position A can be calculated from signals of the plurality of SVs of FIG. 6, in accordance with at least one embodiment of the present disclosure.
FIG. 8 is a schematic illustration of how time inaccuracies can influence the accuracy of the calculated position of FIG. 7, in accordance with at least one embodiment of the present disclosure.
FIG. 9 is a schematic illustration of the signals of FIG. 7 being reflected by buildings, in accordance with at least one embodiment of the present disclosure.
FIG. 10 is a schematic illustration of atmospheric influences on the signals of FIG. 7, in accordance with at least one embodiment of the present disclosure.
FIG. 11 is a flow chart of a process to determine a position of the user device according to one embodiment of the present disclosure.
FIG. 12 is a schematic illustration of a signal transfer between a user device and a space vehicle (SV), in accordance with at least one embodiment of the present disclosure.
FIG. 13 is a contour plot of a cost function obtained from signals of three SVs, in accordance with at least one embodiment of the present disclosure.
FIG. 14 is a contour plot of a cost function obtained from signals of four SVs, in accordance with at least one embodiment of the present disclosure.



DESCRIPTION
A system, apparatus, and method are disclosed for determining the position of a user device, such as an Iridium receiver, which contains a time bias and a frequency bias. Additionally, a system, apparatus, and method are disclosed that allow for estimating the position and velocity of a user device in addition to estimating the time and frequency bias of that same device.
Discussion of Prior Art
One important example of a prior art method for estimating a precise position of a user device from signals from a low-earth orbit (LEO) satellite is discussed in detail as follows. This method is disclosed in U.S. Pat. No. 7,489,926, issued to Whelan, et al. FIGS. 1 through 5 provide a thorough understanding of the various embodiments of this prior art method taught in Whelan, et al.
By way of overview, this prior art method for estimating a precise position of a user device from signals from a low earth orbit (LEO) satellite includes receiving at least one carrier signal at a user device, where each carrier signal is being transmitted from a distinct LEO satellite. The user device processes the carrier signals to obtain a first carrier phase information. Then, the user device recalls an inertial position fix derived at an inertial reference unit. Finally, the user device derives a position of the user device based on the inertial position fix and the first carrier phase information.
FIG. 1 illustrates a system in which Iridium (or other LEO) satellites 12, 14 are used to provide ranging systems to a user in conjunction with one or more reference stations 16, 18. One of the advantages of using Iridium is that it is able to produce a signal that is much stronger than that produced by GPS satellites. Depending on circumstances, the Iridium satellite can be configured to provide users with approximately 20 dB to 40 dB or more received power than GPS.
Positioning using a single ranging source in a three-dimensional dynamic environment with Iridium differs significantly from previous positioning systems in that single ranging sources have been limited to two-dimensioned (2D) resolution on an idealized surface. With the Navy Navigation Satellite System known as TRANSIT, for example, the user was only able to make quasi-static, two-dimensional measurements that were limited in accuracy. Normally, a minimum of four operational TRANSIT satellites were needed to provide the required frequency of precise navigation fixes.
GPS now provides at least four ranging sources simultaneously in order to enable instantaneous, three-dimensional (3D) positioning. However, GPS has a low-power signal that limits operations indoors or in conditions of heavy jamming. A fundamental advantage of the system is that it simultaneously addresses the limitations of its predecessors, providing dynamic, three-dimensional, accurate position fixes, even indoors or in the presence of jamming. Augmented positioning using Iridium should be able to achieve suitable performance that is limited principally by the effects of ambient multipath.
A ground support infrastructure is used to provide differential reference measurements. In the present embodiment, a reference station 16 receives signals from satellites 12 and 14 using reference equipment. Such reference equipment can be functionally identical to user equipment 20 differing only in that the local position of a receiving antenna is precisely known by a survey or other conventional means including GPS positioning.
Differential reference measurement involves the cooperation of at least two receivers, the reference station 16, and the user equipment 20. The cooperation of at least two receivers relies upon a signal 24 received at both the reference station 16 and the user equipment 20, which are both degraded by virtually the same errors. The cooperation is possible on earth when the signals pass through virtually the same slice of atmosphere containing the same obstructions to signals 26. To occur on the surface of the earth, the user equipment 20 and the reference station 16, generally, can be separated by fewer than approximately a thousand kilometers. Where such geometry is present, the signal 24 that reaches both of the user equipment 20 and the reference station 16 will have traveled through the same obstacles 26 or will be augmented by the same pattern of jamming.
The reference station 16 provides real-time measurements of the Iridium clock. A data message 22, which, in the present embodiment, is transmitted over Iridium from the reference station 16 to the user receiver 20, provides a real-time range correction to each measurement to account for both Iridium clock errors and atmospheric effects including obstacles 26 or jamming. Since the reference station 16 has no way of knowing which of the many available satellites the user receiver 20 might be using to calculate its position, the reference receiver 16 quickly runs through all the visible satellites, such as satellite 14, and then computes the error attendant to its signal 28. The corrections necessary to bring the calculated result into line with the known local position of the reference station are then transmitted on any suitable band with adequate confidence in the jamming environment to the user equipment in association with time references to establish near real time correction. In general, navigation performance degrades as the separation between user and reference station gets greater due to attendant differences in obstacles 26 or jamming the signal 24 experiences.
Where a second reference station 18 is suitably close, the second reference station 18 can perform the same calculations on the signal 28 as the first reference station 16 yielding a second correction factor from, for instance, the satellite 14, thereby allowing the user equipment to achieve greater precision by averaging or other suitable means of harmonizing the error calculation.
FIG. 2 depicts a block diagram for the present system architecture for a positioning system 30 that uses Iridium or other LEO satellites. Each component of the positioning system 30 is driven from the same master clock, which is a precise time standard 40. A synthesizer 38 creates each of the requisite coherent sine wave and clock signals for each component based upon a clock signal that is fed to the synthesizer 38 from the precise time standard 40 through a data bus 42.
An antenna 32 is configured to receive transmissions from the Iridium or other LEO satellites as the present embodiment is configured and is optimized for L band reception. An Iridium receiver 34 receives a raw signal that is received at the antenna 32 and compares it with the signal generated by the synthesizer 38 and presented to the receiver 34 at a data bus 48. By comparing the signal at the data bus 48 with the transmission received at the antenna 32, the Iridium receiver 34 presents data sufficient to compute a position solution.
An augmented position solution is calculated using an inertial measurement unit 36 that receives a clock signal from the synthesizer 38. Measuring acceleration with the inertial measurement unit 36 in the present embodiment is accomplished by accelerometers oriented in three orthogonal axes and measuring angular rate about each such axis to compute attitude accurately relative to a vertical axis accomplish accurate attitude sensing. Attitude and other parameters or orientation and motion of the user are derived from the data produced by the accelerometers and rate sensors within the common assembly. In the presently preferred embodiment, the accelerometers are MEMS inertial sensors.
Measuring acceleration with the inertial measurement unit 36 in the present embodiment augments the system to provide a system that anticipates the next position of the user. Optionally, the position solutions derived by use of the inertial measurement unit 36 may be harmonized with earlier solutions to gain a self-testing ability and to reduce a radius of error in the calculation of the position with the inertial measurement unit 36.
Three-dimensional positioning and filtering using Iridium operates over time scales of about 10 minutes, which is much less than the 84-minute Schuler period. The Schuler period is the period for a simple undamped pendulum with a length equivalent to the radius of the earth, and has been used to correct traditional inertial navigation equipment for the curved movement of a spot on the surface of the earth. Therefore, the inertial unit needs to be capable of providing relative position measurements whose accuracy is significantly better than the filtered range measurement accuracy of the Iridium signal.
With MEMS inertial sensors of sufficient performance, degradation due to the ambient multipath of an indoor environment will dominate the overall system-level accuracy. The total system accuracy will start out in the 4-meter range representing one sigma in standard deviation. Advanced signal processing techniques applied to the Iridium signal significantly reduce indoor multipath error. In outdoor applications with an unobstructed view of the sky, the accuracy will be considerably better, which is limited mostly by the performance of the inertial reference unit.
In at least one embodiment, the method and apparatus create a Secure Iridium Broadcast Signal. Although the Iridium signal is technically a TDMA signal, the superposition of several sub-bands together formulate a high-powered signal to appear more like a secure CDMA signal. With such a formulation, a navigation user knows the code in advance to be able to make use of it. If the pulse patterns that make up the secure Iridium Broadcast Signal are programmed correctly, the high-power signal would appear like the secure Y-code signal of GPS or its equivalent for processing.
The systems architecture for the indoor case driven by multipath imposes an implicit requirement on the total position bias of about 1 meter after 10 minutes of coasting. The limiting inertial parameter is likely to come from the gyro-rate bias stability or angle random-walk error. A higher performance inertial system is required if the system is to be used outdoors for high-accuracy and integrity navigation.
The computer 54 serves to tie all the Iridium ranging measurements together, especially when there is only a single ranging source in view at any given time. “High accuracy” means position errors at the centimeter level. “High integrity” is a safety related term that means that there is enough redundant information present in the form of excess satellite ranging measurements to determine if there is an error in the positioning system. Such capability can be used to alert an operator of the system when that system should not be used for navigation. High performance navigation employs the carrier phase of the LEO satellite to attain raw range measurements precise to the centimeter level.
Because the system will often be measuring only one ranging source at a time, it is desirable that a precise frequency standard be used. Two types of frequency standards are available for this purpose: an ovenized quartz crystal oscillator and an atomic rubidium frequency standard. An ovenized quartz crystal as long as the Allan variance at 600 seconds (10 minutes) does not exceed 10−11. This corresponds to about 2 meters of position error over the Iridium pass, which is significantly less than the multipath error on the Iridium signal. If additional accuracy is needed, a compact, ruggedized rubidium standard should be used. The corresponding Allan variance is 10−13, corresponding to a position error of about 2 cm over the 10-minute interval.
Raw position solutions from the Iridium receiver 34 through a data bus 50 and acceleration measurements from the inertial measurement unit 36 through a data bus 52 are fed into a computer 54 which executes a Kalman filter to process the measurements into final solutions. The Kalman filter is a set of mathematical equations that provides an efficient computational (recursive) solution of the least-squares method. The filter is very powerful in several aspects: it supports estimations of past, present, and even future states, and it can do so even when the precise nature of the modeled system is unknown.
The Kalman filter estimates a process by using a form of feedback control: the filter estimates the process state at some time and then obtains feedback in the form of (noisy) measurements. As such, the equations for the Kalman filter fall into two groups: time update equations and measurement update equations. The time update equations are responsible for projecting forward (in time) the current state and error covariance estimates, which are used to obtain the a priori estimates for the next time step. The measurement update equations are responsible for the feedback, i.e. for incorporating a new measurement into the a priori estimate in order to obtain an improved a posteriori estimate.
Since raw position solutions from the Iridium receiver 34 through a data bus 50 and acceleration measurements from the inertial measurement unit 36 through a data bus 52 that are fed into the computer 54 are measurements of the same phenomenon (i.e. movement in space), the measurements are related in the system modeled by the Kalman filter 157 (FIG. 4).
Depending on the circumstances, not all states (such as yaw attitude) will necessarily be observable at all times. However, because of the orbit geometry of Iridium, the system design ensures that the position component of output will effectively always be observable to within the accuracy of the Iridium ranging measurements.
There are two fundamental modes of operation of this system. The first is based on code phase measurements. Inside of a building, there are many sources of multipath; so using the carrier is not especially feasible. However, LEO satellites provide an abundance of geometry, as shown in FIG. 3, along with significantly higher broadcast power that is useful for penetrating physical barriers. The code ranging measurements can be combined using this geometry to solve for reasonably accurate position, using the inertial navigation unit to bridge measurements made at different times.
The second mode of operation is based on carrier phase measurements. If carrier phase measurements are made outdoors, it is possible to obtain a clean line of sight to the LEO satellites, and therefore, achieve centimeter-level positioning accuracy. The same abundance of geometry, as shown in FIG. 3, enables these precision measurements to be combined into high accuracy and high integrity position solutions, thereby again using the inertial navigation unit to bridge measurements made at different times.
FIG. 3 shows a typical geometry pass from the standpoint of the user. The Iridium satellites fly in an arc over an interval of several minutes. Multipath will generally be the largest error source. The Iridium carrier phase can be used to drive the ranging error to be arbitrarily small—potentially to centimeter level—when the user has a clear view of the sky. Unfortunately, raw ranging errors will tend to increase to roughly 20-30 m working indoors. Because the Iridium satellite spans a large-angle arc in the sky, it should be possible to take advantage of the spatial diversity in order to average down much of this indoor multipath error. By analogy with experimental GPS performance, it is possible to predict what Iridium performance is likely to be by scaling the parameters. The correlation of time between Iridium measurements is estimated to be about 10 seconds; meaning that over a 10-minute pass, the receiver can gather roughly 60 “independent” measurements. Therefore, the ranging accuracy may perhaps be improved to roughly 4 meters (dividing the raw ranging accuracy by the square root of 60).
As shown in FIG. 3, a cold start initialization 60 uses a trajectory of the first Iridium satellite pass to define a local section of the Iridium orbit sphere 64 having a zenith 62 relative to the position of the user 68. Inertial navigation by the inertial measurement unit 36 yields a positional covariance after the first pass 66 as shown relative to the position of the user 68. The rapidly changing angle of the orbit of the LEO satellite in the LEO satellite orbit sphere 64 allows for a rapid convergence of the position estimate allowable by means of the LEO satellite in its orbit sphere 64.
The system structure resembles a tightly coupled GPS-Inertial unit. However, as shown in FIG. 4, the system 100 is intended to process as few as a single range measurement at a time using a Kalman filter 150. For dynamic applications, a MEMS Inertial Reference Unit (IRU) 102 is coupled to the system and subjected to error preprocessing in the error-preprocessing unit 105. In more demanding applications, an inertial-grade IRU may be desirable.
A general model for a suitable IRU 102 includes a strap down inertial navigation system. Strap down inertial navigation systems are rigidly fixed to the moving body. Therefore, strap down inertial reference units move with the body, thereby their gyros experience and measure the same changes in angular rate as the body in motion. The strap down inertial reference unit contains accelerometers to measure changes in linear rate in terms of the body's fixed axes. The body's fixed axes serve as a moving frame of reference as opposed to the constant inertial frame of reference. The navigation computer uses the gyros' angular information and the accelerometers' linear information to calculate the body's 3D motion with respect to an inertial frame of reference.
The IRU 102 senses inertial acceleration, which it outputs as rotational acceleration. The rotational acceleration vector information is fed into an error preprocessor 105. The inertial error preprocessor 105 corrects pre-calibrated parameters, including scale factor and alignment errors. Next, the corrected measurements pass through the time update blocks 108 and 111, including the addition of the accelerometer and gyro bias states and the integration of the strap down IRU 102 measurements into position, velocity, and attitude vectors.
At an in-phase coordinate processor 114 and a quatrature coordinate processor 117, a vector translation, xl, and attitude motion, represented by the 3×3 attitude rotation matrix A, of the user platform. With prior knowledge of the antenna mounting lever arm b 120, with respect to the body frame of the user platform, it is possible to use the inertial signal output from the in-phase coordinate processor 114 and the quatrature coordinate processor 117 to project the antenna motion into the line of sight of the satellite, ŝ at a processor 126. The output of the processor 126 is a complex, real-time phase correction. The phase correction is to subtract out short-term user motion and enable long integration times on a LEO signal, when such a LEO signal is available.
On the LEO receiver side of a LEO (in the case of the present embodiment, an IRIDIUM) receiver 132 receives a carrier signal from the LEO satellite. In a present embodiment, a second carrier signal received at a reference ground-station in proximity to the user device is also received in association with the precise position of the ground-station position at an optional data link 135. The second carrier signal insures a rapid integration of the carrier signal from the LEO satellite and further enables operation of an LEO error preprocessor 138.
As with the inertial side, the LEO error preprocessor 138 corrects pre-calibrated parameters, including scale factor and alignment errors. Additionally, the LEO error preprocessor 138 corrects propagation-induced errors based upon the information received at the optional data link 135. The error processor 138 applies corrections such as for atmospheric/ionospheric effects, time tag alignments, and blending code and carrier.
Bias state time update blocks 141, 144, 147, and 151 apply the scalar receiver clock and clock bias estimates to the raw measurements. A further bias block 154, uses the output of the processor 126 to subtract out short-term user motion and enable long integration times on a LEO signal, when such a LEO signal is available. The corrected LEO position is ready for feeding into the Kalman filter 157. In the present embodiment, the computer 54 executes a 17-state Kalman filter estimator to solve for: position (3 axes), velocity (3 axes), accelerometer bias (3 axes), attitude (3 axes), gyro bias (3 axes), clock bias, and clock drift.
A covariance time updater 160 propagates a state covariance estimate. The estimated inertial position, which is projected into the line of sight of each given LEO satellite by the processor 126, is compared with the measured range to the LEO satellite to form the measurement update to the Kalman filter 157.
Referring to FIG. 5, a method 200 is provided for estimating a precise position of a user device in a satellite-based navigation system. At a block 201, a user device receives transmitted carrier signals from a set of LEO satellites. At a block 204, the user device processes the carrier signals to obtain user carrier phase information including geometrically diverse user carrier phase information from the set of LEO satellites. At a block 207, the user device recalls an inertial position fix. At a block 210, the precise position of the user device is determined based on the inertial position fix and the user carrier phase information. At a block 213, the user device derives user carrier information from the set of LEO satellites based upon the inertial position to resolve integer cycle ambiguities in the user carrier phase information.
In a preferred embodiment, the method 200 includes tracking the carrier signals at a reference station in order to obtain reference carrier phase information. The reference carrier phase information includes geometrically diverse reference carrier phase information from the set of LEO satellites. At a block 216, the user device refines the accuracy of the position calculation based upon the reference carrier phase information. In a preferred embodiment, the method further comprises estimating an approximate user position and clock offset using code phase signals received from a set of navigational satellites.
Preferably, differential code phase techniques are used to improve the accuracy of the initial estimate. The preferred embodiment of the method also includes additional advantageous techniques such as: compensating for frequency dependent phase delay differences between carrier signals in user and reference receiver circuits, reading navigation carrier information and LEO carrier information within a predetermined time interval selected in dependence upon an expected motion of the user receiver and the LEO signal sources, calibrating LEO oscillator instabilities using navigation satellite information, compensating for phase disturbances resulting from a bent pipe LEO communication architecture, compensating for oscillator instabilities in the user and reference receivers, predicting present reference carrier phase information based on past reference carrier phase information, and monitoring the integrity of the position calculation.
Depending on the circumstances, not all states (such as yaw attitude) will necessarily be observable at all times. Because of the orbit geometry of Iridium—specifically the rapid large-angle overhead motion—the system ensures that, upon convergence, the position component of output will effectively always be observable to within the accuracy of the Iridium ranging measurements. If high-performance carrier ranging is to be carried out, an optional float bias state is added, one for each LEO satellite, as shown in FIG. 4, to account for the integer cycle ambiguity.
Discussion of Present Disclosure
In the following description, numerous details are set forth in order to provide a more thorough description of the system. It will be apparent, however, to one skilled in the art, that the disclosed system may be practiced without these specific details. In the other instances, well known features have not been described in detail so as not to unnecessarily obscure the system.
FIG. 6 illustrates a plurality of space vehicles (SVs) 600 orbiting the earth 610. The plurality of SVs 600 may include, but is not limited to, satellites, high-flying unmanned air vehicles (UAV), or any other suitable air vehicles. A user device on the earth may include an antenna and a computer system to receive signals from the plurality of SVs 600 and process the signals in order to determine a position of the user device. In some embodiments, the user device can include an output device; for example, a display to visualize the position and/or a connection bus to transfer position data to an external computer system.
In some embodiments, the signals can be modulated or pure carrier signals, while, in other embodiments, the signals can include position data of the SVs 600. In some embodiments, satellites of the global positioning system (GPS) can continuously transmit data signals including their position in space and time signals when the data signals have been sent. Time differences between the GPS satellites and the user device (i.e., the time it takes for the signal to travel from the GPS satellites to the user device) can be used to calculate the position of the user device with respect to each GPS satellite. Consequently, the user device's position on the earth can be determined according to a reference frame. In some embodiments, the reference frame can be the earth centered earth fixed (ECEF) coordinate system.
In some embodiments, the accuracy of the computed position of the user device on the earth can depend on the accuracy of the time signal. As shown in FIG. 7, if the time signals of the user device and the GPS satellites are perfectly synchronized, the user device is capable of determining the travel time of the signals and an exact location A. If the time signals of the user device and the GPS satellites are only slightly offset, as shown in FIG. 8, the position of the user device can be estimated within a certain range, which is denoted by the triangle defined by points B. Unsynchronized time signals can also result from reflections of the signals 900, as shown in FIG. 9, and/or atmospheric/weather conditions 1000, as shown in FIG. 10.
In some embodiments, the user device can process carrier signals from LEO satellites to determine the position of the user device. In some embodiments, the LEO satellites are communication satellites, such as those used in the Iridium system. In some embodiments, the LEO satellites can transmit at a higher intensity than the GPS satellites. As a result, the carrier signals can be received by the user device, even if signals from the GPS satellites become unavailable.
In some embodiments, the position of the user device can be determined from carrier signals of the SVs using a snapshot calculation. The user device can process the carrier signals to obtain carrier phase information. The carrier phase information can include a phase signal, an average slope of the phase signal, a time bias, and a frequency bias. In some embodiments, the carrier phase information can be obtained from an iterative process. The carrier signals can be sampled periodically and/or continuously. In some embodiments, a correlation time of the carrier signals of about ten seconds can suffice to obtain reliable carrier phase information. The snapshot calculation can use the carrier phase information to compute the current position of the user device.
If the carrier signals from the SVs become unavailable, an earlier position can be used to estimate the position of the user device, also called track maintenance. In some embodiments, the track maintenance can be used to extrapolate an estimate for a future position from a current position. The track maintenance can include an extended Kalman filter (EKF) state estimator. The EKF state estimator can model a dynamic evolution of a random process (i.e., arbitrary movements of the user device). In some embodiments, the EKF state estimator can include a velocity of the user device. In some embodiments, the velocity of the user device can be a ground velocity, while in other embodiments the velocity can include velocities faster than the speed of sound. In some embodiments, the EKF state estimator can include an inertial system if the velocity of the user device is greater than the speed of sound. The EKF state estimator can reliably predict a position of the user device if a positioning calculation cannot be performed.
The track maintenance of some embodiments can take a state of the user device into account. In one embodiment, the state of the user device can include a position of the user device, a velocity of the user device, a time bias of the user device, and/or a frequency bias of the user device.
In some embodiments, a hierarchy of different methods to determine the position of the user device can be employed. The hierarchy can include a method using signals from GPS satellites, a method using signals from SVs other than GPS satellites, and/or a method for track maintenance. FIG. 11 illustrates a flow chart 1100 of the process to determine the position of the user device according to one embodiment of the invention. In some embodiments, the hierarchy of the positioning methods can be specific to an application or a mission of the user device. In FIG. 11, it is first determined whether a GPS signal is received 1110. If a GPS signal is received, then the position of the user device is determined 1120. However, if a GPS signal is not received, then it is determined if an Iridium signal is received 1130.
If an Iridium signal is not received, then the position of the user device is estimated based on Kalman filter computation, which includes a velocity of the user device 1140. However, if an Iridium signal is received, then the time bias and frequency bias are computed 1150. After the time bias and frequency bias are computed, the cost function, which is a measure of the expected measurement error as a function of the estimated position, is then computed 1160. And, after the cost function is computed, the position of the user device is determined by minimizing the cost function 1170.
FIG. 12 illustrates a signal transfer between the user device and a distinct SV. FIG. 12 can be used to explain the snapshot calculations. A vector Xu, can point from an origin O to the user device {dot over (X)}, which can be also called a receiver. During a known transmission time t, the SV can be located at a position, which can be described by the vector Xt,SV. The SV can send a signal to the user device {dot over (X)} so that the user device can uniquely identify the SV. The distance between the user device and the SV during the known transmission time t can be called range rt,SV. The range rt,SV can be regarded as the length of the difference between the vector Xu and the vector Xt,SV. During the known transmission time, the SV can travel at a speed {dot over (X)}t,SV so that a range rate {dot over (r)}t,SV can be defined as




  
    
      r
      .
    
    
      t
      ,
      SV
    
  
  ⁢
  
    =
    ^
  
  ⁢
  
    
      -
      
        
          X
          .
        
        
          t
          ,
          SV
        
        T
      
    
    ⁢
    
      
        
          
            X
            u
          
          -
          
            X
            
              t
              ,
              SV
            
          
        
        
          r
          
            t
            ,
            SV
          
        
      
      .
    
  




The accuracy of the range rt,SV and the range rate {dot over (r)}t,SV can be compromised by a number of influences. In some embodiments, the method to compute the position of the user device can include a time bias bt,u, a frequency bias {dot over (b)}t,u, and random measurement errors εt,ηt. The time bias bt,u can be regarded as a time drift and the frequency bias {dot over (b)}t,u can be regarded as a frequency offset. A pseudo range ρt,SV and a pseudo range rate {dot over (ρ)}t,SV can be defined for each SV as follows:

ρt,SV=rt,SV(Xu)+bt,u+εt,

{dot over (ρ)}t,SV={dot over (r)}t,SV(Xu)+{dot over (b)}t,u+ηt 

wherein the random measurement errors εt,ηt can be specific to the user device and the SV.

The time bias bt,u can be approximated by a polynomial in time defined by




  
    b
    
      t
      ,
      u
    
  
  ⁢
  
    =
    ^
  
  ⁢
  
    
      ∑
      
        n
        =
        0
      
      N
    
    ⁢
    
      
        q
        n
      
      ⁢
      
        
          t
          n
        
        /
        
          
            n
            !
          
          .
        
      
    
  




The order of the polynomial N can be any integer number. In some embodiments, the order to the polynomial can be governed by the available compute power of the user device. The coefficients qn of the polynomial can be independent of time. In some embodiments, the coefficients qn of the polynomial can be constant. The frequency bias {dot over (b)}t,u can be a derivative of the time bias bt,u with respect to time resulting





  
    in
    ⁢
    
      
    
    ⁢
    
      
        b
        .
      
      
        t
        ,
        u
      
    
  
  ⁢
  
    =
    ^
  
  ⁢
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        q
        n
      
      ⁢
      
        
          t
          
            n
            -
            1
          
        
        /
        
          
            
              (
              
                n
                -
                1
              
              )
            
            !
          
          .
        
      
    
  




A nonlinear least squares method can be employed to estimate the pseudo range ρt,SV and the pseudo range rate {dot over (ρ)}t,SV by minimizing a cost function Jcost, which can be defined as:



  
    
      
        J
        cost
      
      ⁢
      
        =
        ^
      
      ⁢
      
        
          1
          2
        
        ⁢
        
          ∑
          
            
              
                [
                
                  
                    Z
                    
                      t
                      ,
                      SV
                    
                  
                  -
                  
                    W
                    
                      t
                      ,
                      SV
                    
                  
                  -
                  
                    
                      Ψ
                      t
                    
                    ⁢
                    
                      q
                      u
                    
                  
                
                ]
              
              T
            
            ⁢
            
              
                R
                
                  t
                  ,
                  SV
                
                
                  -
                  1
                
              
              ⁡
              
                [
                
                  
                    Z
                    
                      t
                      ,
                      SV
                    
                  
                  -
                  
                    W
                    
                      t
                      ,
                      SV
                    
                  
                  -
                  
                    
                      Ψ
                      t
                    
                    ⁢
                    
                      q
                      u
                    
                  
                
                ]
              
            
          
        
      
    
    ,
    
      

    
    ⁢
    with
  




  
    
      
        Z
        
          t
          ,
          SV
        
      
      =
      
        [
        
          
            
              
                ρ
                
                  t
                  ,
                  SV
                
              
            
          
          
            
              
                
                  ρ
                  .
                
                
                  t
                  ,
                  SV
                
              
            
          
        
        ]
      
    
    ,
    
      

    
    ⁢
    
      
        W
        
          t
          ,
          SV
        
      
      =
      
        [
        
          
            
              
                r
                
                  t
                  ,
                  SV
                
              
            
          
          
            
              
                
                  r
                  .
                
                
                  t
                  ,
                  SV
                
              
            
          
        
        ]
      
    
    ,
    
      

    
    ⁢
    
      
        Ψ
        t
      
      =
      
        [
        
          
            
              
                1
                ⁢
                
                  
                
                ⁢
                t
                ⁢
                
                  
                
                ⁢
                ...
                ⁢
                
                  
                
                ⁢
                
                  
                    t
                    N
                  
                  /
                  
                    N
                    !
                  
                
              
            
          
          
            
              
                0
                ⁢
                
                  
                
                ⁢
                1
                ⁢
                
                  
                
                ⁢
                ...
                ⁢
                
                  
                
                ⁢
                
                  
                    t
                    
                      N
                      -
                      1
                    
                  
                  /
                  
                    
                      (
                      
                        N
                        -
                        1
                      
                      )
                    
                    !
                  
                
              
            
          
        
        ]
      
    
    ,
    
      

    
    ⁢
    
      
        q
        u
      
      =
      
        
          [
          
            
              q
              o
            
            ⁢
            
              
            
            ⁢
            
              q
              1
            
            ⁢
            
              
            
            ⁢
            ...
            ⁢
            
              
            
            ⁢
            
              q
              N
            
          
          ]
        
        T
      
    
    ,
  



and Rt,SV being a measurement error covariance matrix or other weighting matrix.

In some embodiments, a linear least squares method can be employed to estimate the time bias bt,u and the frequency bias {dot over (b)}t,u. The estimate for the time bias can be labeled bt,u and the estimate for the frequency bias can be labeled {dot over (b)}t,u. To calculate the estimates, the cost function Jcost can be reduced to two dimensions, latitude and longitude, because the coefficients qn are linearly dependent on the pseudo range ρt,SV and the pseudo range rate {dot over (ρ)}t,SV.
For reducing the order of the cost function Jcost to two dimensions, a linear least squares sub-problem can be solved. To this end, the vector qu containing the coefficients qn can be re-written as:



  
    
      
        q
        u
      
      ⁢
      
        =
        ^
      
      ⁢
      
        
          q
          z
        
        -
        
          q
          w
        
      
    
    ,
    
      

    
    ⁢
    wherein
  




  
    
      q
      z
    
    ⁢
    
      =
      ^
    
    ⁢
    
      
        
          (
          
            
              ∑
              
                τ
                ,
                SV
              
            
            ⁢
            
              
                Ψ
                τ
                T
              
              ⁢
              
                R
                
                  τ
                  ,
                  SV
                
                
                  -
                  1
                
              
              ⁢
              
                Ψ
                τ
              
            
          
          )
        
        
          -
          1
        
      
      ⁢
      
        (
        
          
            ∑
            
              τ
              ,
              SV
            
          
          ⁢
          
            
              Ψ
              τ
              T
            
            ⁢
            
              R
              
                τ
                ,
                SV
              
              
                -
                1
              
            
            ⁢
            
              Z
              
                τ
                ,
                SV
              
            
          
        
        )
      
    
  




  and




  
    
      q
      w
    
    ⁢
    
      =
      ^
    
    ⁢
    
      
        
          (
          
            
              ∑
              
                τ
                ,
                sv
              
            
            ⁢
            
              
                Ψ
                τ
                T
              
              ⁢
              
                R
                
                  τ
                  ,
                  sv
                
                
                  -
                  1
                
              
              ⁢
              
                Ψ
                τ
              
            
          
          )
        
        
          -
          1
        
      
      ⁢
      
        
          (
          
            
              ∑
              
                τ
                ,
                sv
              
            
            ⁢
            
              
                Ψ
                τ
                T
              
              ⁢
              
                R
                
                  τ
                  ,
                  sv
                
                
                  -
                  1
                
              
              ⁢
              
                W
                
                  τ
                  ,
                  sv
                
              
            
          
          )
        
        .
      
    
  



Replacing the vector q, in the cost function Jcost yields in a reduced order cost function {tilde over (J)}cost in term of measurements {tilde over (Z)}, model residuals {tilde over (W)}, and the error covariance matrix {tilde over (R)} correlating the variances of {tilde over (Z)}. The bias and bias rate estimates are then recovered as




  
    
      [
      
        
          
            
              
                b
                _
              
              
                t
                ,
                u
              
            
          
          
            
              (
              
                X
                u
              
              )
            
          
        
        
          
            
              
                b
                
                  _
                  .
                
              
              
                t
                ,
                u
              
            
          
          
            
              (
              
                X
                u
              
              )
            
          
        
      
      ]
    
    =
    
      
        
          Z
          _
        
        t
      
      -
      
        
          W
          _
        
        t
      
    
  
  ,




with Zt{circumflex over (=)}Ψtqz and Wt{circumflex over (=)}Ψtqw. In some embodiments, the reduction of order of the cost function Jcost can save a substantial amount of computation time.

Neglecting altitude, the position of the user device Xu can be a function of latitude y1 and longitude y2 alone:



  
    
      
        X
        u
      
      =
      
        
          N
          ⁡
          
            (
            
              y
              1
            
            )
          
        
        ⁡
        
          [
          
            
              
                
                  
                    cos
                    ⁡
                    
                      (
                      
                        y
                        1
                      
                      )
                    
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      
                        y
                        2
                      
                      )
                    
                  
                
              
            
            
              
                
                  
                    cos
                    ⁡
                    
                      (
                      
                        y
                        1
                      
                      )
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      
                        y
                        2
                      
                      )
                    
                  
                
              
            
            
              
                
                  
                    
                      (
                      
                        b
                        /
                        a
                      
                      )
                    
                    2
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      
                        y
                        1
                      
                      )
                    
                  
                
              
            
          
          ]
        
      
    
    ,
    
      
    
    ⁢
    
      

    
    ⁢
    with
  




  
    
      
        N
        ⁡
        
          (
          
            y
            1
          
          )
        
      
      ⁢
      
        =
        ^
      
      ⁢
      
        a
        
          
            1
            -
            
              
                (
                
                  1
                  -
                  
                    
                      (
                      
                        b
                        /
                        a
                      
                      )
                    
                    2
                  
                
                )
              
              ⁢
              
                
                  sin
                  2
                
                ⁡
                
                  (
                  
                    y
                    1
                  
                  )
                
              
            
          
        
      
    
    ,
  



a being the major axis of the earth ellipsoid and b being the minor axis of the earth ellipsoid according to the WGS-84 earth surface model (see FIG. 12).

In some embodiments, the latitude y1 and the longitude y2 can be evaluated with an iterative method, like a Newton Raphson, by solving yk+1=yk−Dk(yk)∇{tilde over (J)}cost(yk), with y being a vector having y1 and y2 as components, k being an iteration step, Dk being a gradient deflection matrix and ∇{tilde over (J)}cost being the gradient of the reduced order cost function {tilde over (J)}cost.
The gradient deflection matrix Dk can be the inverse Hessian of the reduced order cost function




  
    
      J
      ~
    
    cost
  
  ,
  
    
  
  ⁢
  
    
      
        D
        k
      
      ⁡
      
        (
        
          y
          k
        
        )
      
    
    =
    
      
        [
        
          
            
              ∂
              2
            
            
              
                ∂
                2
              
              ⁢
              y
            
          
          ⁢
          
            
              
                J
                ~
              
              cost
            
            ⁡
            
              (
              y
              )
            
          
        
        ]
      
      
        y
        =
        yk
      
      
        -
        1
      
    
  
  ,




which can be required to be evaluated at each iterative step k. As a result, the evaluation of Dk can be computationally intensive.

FIGS. 13 and 14 illustrates surface plots of the reduced order cost function {tilde over (J)}cost computed with signals from three different SVs and four different SVs, respectively. Because the surfaces are well behaved, the gradient deflection matrix of some embodiments can be approximated by




  
    
      D
      k
    
    ⁡
    
      (
      
        y
        k
      
      )
    
  
  ⁢
  
    =
    ^
  
  ⁢
  
    
      [
      
        
          ∑
          
            t
            ,
            sv
          
        
        ⁢
        
          
            
              ∂
              
                
                  
                    W
                    ~
                  
                  
                    t
                    ,
                    sv
                  
                  T
                
                ⁡
                
                  (
                  
                    
                      X
                      u
                    
                    ⁡
                    
                      (
                      y
                      )
                    
                  
                  )
                
              
            
            
              ∂
              y
            
          
          ⁢
          
            
              R
              ~
            
            
              t
              ,
              sv
            
            
              -
              1
            
          
          ⁢
          
            
              ∂
              
                
                  
                    W
                    ~
                  
                  
                    t
                    ,
                    sv
                  
                
                ⁡
                
                  (
                  
                    
                      X
                      u
                    
                    ⁡
                    
                      (
                      y
                      )
                    
                  
                  )
                
              
            
            
              ∂
              
                y
                T
              
            
          
        
      
      ]
    
    
      y
      =
      yk
    
    
      -
      1
    
  




involving only first-order derivatives with respect to y. As a result, the computation time for the snapshot calculation can be significantly reduced.

The iterative step can be initialized with a first guess of the position of the user device. In some embodiments, a weighted average of one or more SV beam positions projected onto the earth can be used as the initial guess, while in other embodiments, a weighted average of one or more SV positions projected onto the earth can serve as the initial guess. Once the iterative step converges with a desired accuracy, the snapshot calculation can be completed resulting in an (instantaneous) position X, of the user device along with estimates of the time and frequency biases.
To predict and/or estimate a position of the user device, the EKF state estimator can track the (changing) position of the user device. A state vector xt for the user device can be defined as xt{circumflex over (=)}[Xt,u, {dot over (X)}t,u, bt,u, {dot over (b)}t,u]T, with Xt,u being the position Xu at time instance t, {dot over (X)}t,u being the velocity of the user device at time instance t, and bt,u and {dot over (b)}t,u defined as above for the snapshot calculation. In some embodiments, the EKF state estimator can track a plurality of user devices simultaneously.
In some embodiments, the EKF state estimator can include four distinct sub-steps, a state extrapolation step, a covariance extrapolation step, a state update step and covariance update step. The state vector itself xt, can be modeled as a random process according to the equation xt=Φt,τxτ+ξr, wherein Φt,τ is a state transition matrix ξτand is the process noise with covariance Qt,τ.
The covariance matrix Qt,τ can be used to tune the EKF state estimator to achieve good agreement. The covariance matrix Qt,τ can be regarded as a process noise or diffusion matrix. Because the state vector xt accounts for a position and a velocity of the user device, a neglected acceleration of the user device can be compensated in the covariance matrix Qt,τ. In some embodiments, the covariance matrix Qt,τ can be defined as:



  
    
      
        Q
        
          t
          ,
          τ
        
      
      =
      
        [
        
          
            
              
                
                  
                    Γ
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  3
                
                /
                3
              
            
            
              
                
                  
                    Γ
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  2
                
                /
                2
              
            
            
              
                0
                
                  3
                  ×
                  1
                
              
            
            
              
                0
                
                  3
                  ×
                  1
                
              
            
          
          
            
              
                
                  
                    Γ
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  2
                
                /
                2
              
            
            
              
                Γ
                ⁡
                
                  (
                  
                    t
                    -
                    τ
                  
                  )
                
              
            
            
              
                0
                
                  3
                  ×
                  1
                
              
            
            
              
                0
                
                  3
                  ×
                  1
                
              
            
          
          
            
              
                0
                
                  1
                  ×
                  3
                
              
            
            
              
                0
                
                  1
                  ×
                  3
                
              
            
            
              
                
                  
                    
                      σ
                      d
                      2
                    
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  3
                
                /
                3
              
            
            
              
                
                  
                    
                      σ
                      d
                      2
                    
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  2
                
                /
                2
              
            
          
          
            
              
                0
                
                  1
                  ×
                  3
                
              
            
            
              
                0
                
                  1
                  ×
                  3
                
              
            
            
              
                
                  
                    
                      σ
                      d
                      2
                    
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  2
                
                /
                2
              
            
            
              
                
                  σ
                  d
                  2
                
                ⁡
                
                  (
                  
                    t
                    -
                    τ
                  
                  )
                
              
            
          
        
        ]
      
    
    ,
    
      

    
    ⁢
    with
  




  
    
      Γ
      =
      
        
          [
          
            
              u
              1
            
            ⁢
            
              u
              2
            
            ⁢
            
              u
              3
            
          
          ]
        
        ⁢
        
          
            Λ
            ⁡
            
              [
              
                
                  u
                  1
                
                ⁢
                
                  u
                  2
                
                ⁢
                
                  u
                  3
                
              
              ]
            
          
          T
        
      
    
    ,
    
      
    
    ⁢
    
      Λ
      =
      
        2
        ⁢
        
          
            τ
            m
          
          ⁡
          
            [
            
              
                
                  
                    σ
                    xy
                    2
                  
                
                
                  0
                
                
                  0
                
              
              
                
                  0
                
                
                  
                    σ
                    xy
                    2
                  
                
                
                  0
                
              
              
                
                  0
                
                
                  0
                
                
                  
                    σ
                    z
                    2
                  
                
              
            
            ]
          
        
      
    
    ,
  



ui being unit vectors, τm being a maneuver time constant, σd being an uncertainty in the bias acceleration, σxy being a planar maneuver uncertainty, and σz being a vertical maneuver uncertainty.

The unit vectors [u1 u2 u3] can be oriented to a local tangent plane centered at the estimate of the position of the user device Xt,u. In some embodiments, the unit vectors ui can be evaluated at the initial position estimate computed by the snapshot calculation and remain constant throughout the track maintenance. In some embodiments, a component of the covariance matrix Qt,τ and/or the uncertainties σi can be empirically determined.
Continuing with the state extrapolation step, the state transition matrix Φt,τ can be defined in some embodiments as




  
    
      Φ
      
        t
        ,
        τ
      
    
    ⁢
    
      =
      ^
    
    ⁢
    
      [
      
        
          
            
              I
              
                3
                ×
                3
              
            
          
          
            
              
                (
                
                  t
                  -
                  τ
                
                )
              
              ⁢
              
                I
                
                  3
                  ×
                  3
                
              
            
          
          
            
              0
              
                3
                ×
                3
              
            
          
          
            
              0
              
                3
                ×
                3
              
            
          
        
        
          
            
              0
              
                3
                ×
                3
              
            
          
          
            
              I
              
                3
                ×
                3
              
            
          
          
            
              0
              
                3
                ×
                3
              
            
          
          
            
              0
              
                3
                ×
                3
              
            
          
        
        
          
            0
          
          
            0
          
          
            1
          
          
            
              (
              
                t
                -
                τ
              
              )
            
          
        
        
          
            0
          
          
            0
          
          
            0
          
          
            1
          
        
      
      ]
    
  
  ,




wherein I is the identity matrix and (t−τ) is an extrapolation time interval. In some embodiments, the state transition matrix Φt,τ can be modeled as a linear prediction for the state vector xt.

In some embodiments, the measurement model can be substantially identical to the model described with respect to the snapshot calculation above, i.e., Zt,SV. A measurement vector zt can be defined as



  
    
      
        z
        t
      
      =
      
        
          [
          
            
              
                
                  ρ
                  
                    t
                    ,
                    u
                    ,
                    SV
                  
                
              
            
            
              
                
                  
                    ρ
                    .
                  
                  
                    t
                    ,
                    u
                    ,
                    SV
                  
                
              
            
          
          ]
        
        =
        
          
            h
            ⁡
            
              (
              
                x
                1
              
              )
            
          
          +
          
            η
            t
          
        
      
    
    ,
    
      

    
    ⁢
    wherein
  




  
    
      h
      ⁡
      
        (
        
          x
          t
        
        )
      
    
    ⁢
    
      =
      ^
    
    ⁢
    
      [
      
        
          
            
              
                r
                
                  t
                  ,
                  u
                  ,
                  sv
                
              
              +
              
                b
                
                  t
                  ,
                  u
                
              
            
          
        
        
          
            
              
                
                  r
                  
                    .
                    
                      
                    
                  
                
                
                  t
                  ,
                  u
                  ,
                  sv
                
              
              +
              
                
                  b
                  .
                
                
                  t
                  ,
                  u
                
              
            
          
        
      
      ]
    
  



is a function which estimates the predicted measurement based only on the present state estimate. In the previous expression, rt,u,sv, {dot over (r)}t,u,sv are the estimated range and range rate between the user and the SV and ηt is the measurement noise with covariance Rt.

For the covariance extrapolation step, a prediction {circumflex over (x)}t− of the state vector xt can be computed with {circumflex over (x)}t−=Φt,τ{circumflex over (x)}τ−. A prediction for the covariance matrix Pt for the state vector xt can be evaluated by Pt−=Φt,τPτ−Φt,τT+Qt,τ. In the above equations, the superscript ‘−’ denotes that the estimates of xt and Pt have not been updated by the most recent measurement value. The covariance update step (denoted by the superscript ‘+’) can use this data to compute a new estimate {circumflex over (x)}t+ for the state vector xt. The new estimate {circumflex over (x)}t+ can be computed with {circumflex over (x)}t+={circumflex over (x)}t−+Pt−HT({circumflex over (x)}t−)(H({circumflex over (x)}t−)Pt−HT({circumflex over (x)}t−)+Rt)−1 (zt−h({circumflex over (x)}t−). The covariance matrix Pt can be updated with Pt+=Pt−−Pt−HT({circumflex over (x)}t−)(H({circumflex over (x)}t−)Pt−HT({circumflex over (x)}t−)+Rt)−1H({circumflex over (x)}t−)Pt−, with H(xt) being the Jacobian matrix of h(xt). The Jacobian matrix H(xt) can be calculated with




  
    
      H
      ⁡
      
        (
        
          x
          t
        
        )
      
    
    ⁢
    
      =
      ^
    
    ⁢
    
      [
      
        
          
            
              S
              
                t
                ,
                u
                ,
                sv
              
              T
            
          
          
            
              0
              
                1
                ×
                3
              
            
          
          
            1
          
          
            0
          
        
        
          
            
              
                (
                
                  
                    d
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        X
                        .
                      
                      
                        t
                        ,
                        u
                        ,
                        sv
                      
                    
                  
                  -
                  
                    
                      
                        r
                        .
                      
                      
                        t
                        ,
                        u
                        ,
                        sv
                      
                    
                    ⁢
                    
                      S
                      
                        t
                        ,
                        u
                        ,
                        sv
                      
                    
                  
                
                )
              
              T
            
          
          
            
              S
              
                t
                ,
                u
                ,
                sv
              
              T
            
          
          
            0
          
          
            1
          
        
      
      ]
    
  
  ,




wherein st,u,sv{circumflex over (=)}dXt,u,sv/rt,u,sv is a unit vector pointing from the SV to the user device, d{dot over (X)}t,u,sv{circumflex over (=)}({dot over (X)}t,u−{dot over (X)}t,sv) is a relative velocity between the user device and the SV, and {dot over (r)}t,u,sv{circumflex over (=)}d{dot over (X)}t,u,svTst,u,sv is the range rate from the SV to the user device.

In some embodiments, the covariance extrapolation step and the covariance update step can be regarded as a predictor-corrector step method for the state vector xt and the covariance matrix Pt. To compute the state vector xt of the user device, the EKF state estimator of some embodiments can use different extrapolation time intervals. In some embodiments, the correlation time of a few seconds can suffice to obtain accurate results for the state vector xt.
In some embodiments, the track maintenance can give reliable results in occluded environments. Types of occluded environments can include environments in which signals from other positioning systems are jammed or otherwise unavailable. In some embodiments, the occluded environments can include indoor environments, like buildings, basements, tunnels, caves, and undersea locations. In some embodiments, the accuracies of the snapshot calculations obtained for an Iridium system can be about 1,200 meters for three SVs and an occluded environment, 400 meters for four SVs and an occluded environment, and 32 meters for four SVs and a clean environment.
Different objects and advantages of the disclosure are described. It is understood that not necessarily all such objects or advantages may be achieved in accordance with any particular example of the disclosure. Those skilled in the art will recognize that the disclosure may be embodied or carried out in a manner that achieves or optimizes one advantage or group of advantages as taught herein without necessarily achieving other objects or advantages as may be taught or suggested.
While the method and system have been described in terms of what are presently considered to be the most practical and preferred examples, it is to be understood that the disclosure need not be limited to the disclosed examples. It is intended to cover various modifications and similar arrangements included within the spirit and scope of the claims, the scope of which should be accorded the broadest interpretation so as to encompass all such modifications and similar structures. The present disclosure includes any and all examples of the following claims.

## Claims

1. A method for track maintenance of a user device comprising:
receiving a signal from at least one space vehicle;
obtaining a position of the user device based on the received signal;
applying an extended Kalman filter state estimator including the position of the user device and a velocity of the user device; and
obtaining an estimate of a state of the user device for the track maintenance,
wherein obtaining the position of the user device includes a snapshot calculation, and
wherein performing at least one of the snapshot calculation and the extended Kalman filter state estimator includes a pseudo range measurement model and a pseudo range rate measurement model.

2. The method of claim 1, wherein computing the pseudo range measurement model and the pseudo range rate measurement model includes a nonlinear least squares method.

3. The method of claim 1, wherein performing at least one of the snapshot calculation and the extended Kalman filter state estimator includes a time bias and a frequency bias.

4. The method of claim 3, wherein computing the time bias and the frequency bias includes a linear least squares method.

5. The method of claim 3, wherein computing the time bias includes approximating the time bias with a polynomial.

6. The method of claim 3, wherein computing the frequency bias includes taking a derivative of the time bias with respect to time.

7. The method of claim 3, wherein the time bias and the frequency bias are expressed as explicit functions of a position of the user device.

8. The method of claim 3, wherein expressions of the time bias and the frequency bias modify a nonlinear least squares cost function by removing the biases as independent states, thereby reducing the number of independent variables to those of a position of the user device alone.

9. The method of claim 1, wherein the signal is received from a low earth orbit (LEO) satellite.

10. The method of claim 9, wherein the LEO satellite is part of the Iridium system.

11. The method of claim 1, wherein the state of the user device includes a position of the user device, a velocity of the user device, a time bias, and a frequency bias.

12. The method of claim 1, wherein obtaining the estimate of the state of the user device includes evaluating a latitude and a longitude of the position of the user device.