Location Determination Using Doppler and Pseudorange Measurements From Fewer Than Four Satellites

This invention relates to use of signals received from fewer than four satellites to determine location of a signal receiver.

Location determination (LD) systems, such as the Global Positioning System (GPS), the Global Orbiting Navigational Satellite System (GLONASS) and the Low Earth Orbit (LEO) System, receive and analyze timed signals from a plurality of satellites (usually three or more) in order to estimate the location and time offset of an LD signal receiver. The LD signals received may be used to measure pseudoranges from the receiver to each satellite, carrier phases, signal Doppler shifts and the like, but it is usually necessary to receive and analyze these signals from many satellites in order to estimate the location coordinates and time offset or related quantities.

In some situations, acceptable signals from an adequate number of satellites is not available, because of signal interference from an adjacent urban canyon or other structure(s), because of an unfortunate geometric configuration for the contributing satellites, because of partial or full masking of signals transmitted near the horizon, or because of non-availability of signal processing in one or more receiver channels. What is needed is a signal receiver system that can receive and analyze LD signals from fewer than four satellites to estimate location coordinates and receiver time offset, or other relevant variables. Preferably, the system should work with LD signals received from as few as one or two LD signal sources (usually, but not necessarily, satellites) at substantially the same signal receipt time. Preferably, the system should work with an arbitrary satellite configuration, not merely with satellites in equatorial orbits.

These needs are met by the invention, which provides method and apparatus that receives LD signals from fewer than four LD signal sources, measures pseudorange, signal Doppler shift and selected time derivatives of Doppler shift, and estimates the location coordinates and time offset ("location fix coordinates") for the LD signal receiver. In one approach, pseudorange and Doppler shift (no time derivatives) for each of two LD signal sources are measured and analyzed to estimate the location fix coordinates. In another approach, pseudorange and Doppler shift (no time derivatives) for a single moving LD signal source are measured and analyzed at two distinct times to estimate the location fix coordinates. In another approach, any two of the pseudorange and the zeroth, first and second time derivatives of Doppler shift at two distinct times are measured and analyzed to estimate the location fix coordinates. In another approach, pseudorange and at least one of the zeroth, first and second time derivatives of Doppler shift for a single moving LD signal source are measured and analyzed at two distinct times to estimate the location fix coordinates. In another approach, pseudorange, Doppler shift and the first and second time derivatives of Doppler shift are measured and analyzed at a selected time to estimate the location fix coordinates. In another approach, any combination of four of the pseudorange and the zeroth, first and second time derivatives of Doppler shift at two or more distinct times are measured and analyzed to estimate the location fix coordinates. If one or more of the location fix coordinates (x.sub.r,y.sub.r,z.sub.r,b.sub.r) is known, the number of measurements is reduced accordingly.

FIG. 1 illustrates a situation in which the invention can be used.

FIGS. 2-4 are flow charts illustrating methods to practice the invention.

FIG. 5 illustrates analysis of an elliptical satellite orbit.

FIG. 6 illustrates apparatus to practice the invention.

In FIG. 1, an LD signal antenna 12 and receiver 11 (referred to collectively as an "LD receiver") receives LD signals from one or more moving LD signal sources 13A, 13B, 13C, whose location coordinates are known as a function of time t and which move relative to the receiver. The LD receiver 11 rotates with the Earth's surface, has initial location coordinates (x.sub.r0,y.sub.r0,z.sub.r0), and has fixed receiver time offset .DELTA.t.sub.r with b.sub.r =c.multidot..DELTA.t.sub.r, where c is the velocity of light in the ambient medium. The location fix coordinates (x.sub.r,y.sub.r,z.sub.r,b.sub.r) for the LD receiver 11 are to be estimated, using LD signals received from one, two or three LD signal sources.

The LD signal receiver 11 occupies a fixed location relative to the Earth's surface as the Earth rotates on its axis. Expressed as a function of time relative to an absolute coordinate system in space, the location coordinates for the LD receiver 11 can be expressed in a Cartesian system, with origin at the Earth's center of rotation and z-axis coinciding with the Earth's rotation axis, as

where x.sub.r0, y.sub.r0 and z.sub.r0 are the (unknown) location coordinate values of the LD receiver at a selected time t'=t-t0 and .omega..sub.e is the known angular velocity of the Earth.

The first, second and third time derivatives of the location coordinates for the LD signal receiver 11 are found to be

The LD signal source(s) is/are assumed to be moving, with known location coordinates as a function of time. For definiteness, each LD signal source is assumed to be carried on a satellite that revolves around the Earth in a circular orbit of radius a with an orbital plane that includes the center of the Earth. Each satellite rotates with a known characteristic angular velocity .omega..sub.s, in an orbital plane that is described by a known azimuthal angle .phi. and by a known angle .theta. that this plane makes with the equatorial line on the Earth. A satellite thus has a time-varying location given by the Cartesian coordinates

The time variable t" in Eqs. (13)-(15) is a translated time variable t"=t-t(j), to account for the difference in location of each satellite (j). The parameters a, .omega..sub.s, .phi. and .theta. also vary with the satellite (j). The first, second and third time derivatives of the location coordinates of the satellite are found to be

A computer 14 is programmed to handle the computations.

The square of the measured pseudorange is first formed. If the vectors r.sub.r (t) and r.sub.s (t) represent the respective vectors from the coordinate origin to the LD signal receiver (i) and to the satellite (j), the (corrected) pseudorange equation for a selected receiver (i) and a selected satellite (j) becomes ##EQU1## The unknowns in Eq. (25) are the location coordinates x.sub.r0, y.sub.r0, z.sub.r0 and the receiver time offset b.sub.r (optionally determined). After simplification, the time derivative of the separation distance (squared) becomes ##EQU2##

More generally, an nth power of the quantity (.DELTA.r.multidot..DELTA.v) becomes

The first time derivative of Doppler shift involves the quantity ##EQU3##

The second derivative of Doppler shift involves the quantity ##EQU4##

For a first body whose motion can be described by the usual special relativistic equations, the frequency Doppler shift of a signal emitted by that body and received by a second body is given by

(C. Moller, The Theory of Relativity, Oxford Univ. Press, 1952, pp. 58-62), where .DELTA.r and .DELTA.v are the first body-second body location vector difference and velocity vector difference, respectively, and .nu..sub.0 is the signal frequency as transmitted by the first body.

Using the cosine law, the quantity

in the denominator in Eq. (45) can be replaced by an approximation

where a.sub.e is the Earth's radius and .psi. is the estimated angle between a vector directed toward the LD signal receiver and a vector directed toward the satellite whose signal is being received (.pi./2.ltoreq..psi..ltoreq..pi.). The numerical value .vertline..DELTA.r.vertline. may be written as .chi.(a-a.sub.e), where a=26,560 km, a.sub.e =7,040 km and 1.ltoreq..chi..ltoreq..chi.(max)={(a+a.sub.e)/(a-a.sub.e)}.sup.1/2 =1.312 for GPS orbits, and a=25,510 km, a.sub.e =7,040 km and 1.ltoreq..chi..ltoreq..chi.(max)=1.328 for GLONASS orbits, for example. Preferably, the quantity .vertline..DELTA.r.vertline. is estimated using a measured pseudorange value.

The first and second time derivatives of the frequency Doppler shift between the first body and the second body are given by

In a first embodiment of the invention, a pseudorange (optional) and Doppler shifts .DELTA..nu.(Doppler;t;i;j) are measured for three different (non-collinear) satellites, j=1, 2, 3, at a selected time t=t1. From Eqs. (26), the three Doppler shift relations become ##EQU5## where the quantities in the 3.times.3 matrix M and in the 3.times.1 column vector G are calculated or measured. Equation (40) is then inverted to provide the receiver (initial) location coordinates ##EQU6## The matrix M(t1;i;1,2,3) is invertible, in part because the three satellites are non-collinear. Equations (1), (2) and (3) can be applied to estimate the receiver location coordinates on the rotating Earth at any other time t. If the receiver clock offset b.sub.r is also needed (optional), the pseudorange PR(t1;i;j) is also measured, and Eq. (25) is applied to determine the offset value br.

In a second embodiment, three or four of pseudorange and the zeroth, first and second time derivatives of the Doppler shift for a single satellite (j) are measured and/or calculated for a selected time t=t1. The three Eqs. (36), (39) and (40), supplemented by Eqs. (26), (32) and (34), can be re-expressed in the forms ##EQU7## where M'(t1;i;j)=(M'.sub.k1) is a 3.times.3 matrix with entries

Equation (46) is then inverted to provide the receiver (initial) location coordinates ##EQU8## The matrix M'(t1;i;1,2,3) should be verified to be invertible by computation of its determinant. Equations (1), (2) and (3) can be applied to estimate the present receiver location coordinates on the rotating Earth at any time t. Equation (25) is applied to estimate the time offset b.sub.r for the receiver, if this quantity is needed. Equations (41), (46) and (71)-(80) are exact if separation distance .vertline..DELTA.r(t1;i;j).vertline. can be measured or otherwise determined. If receiver clock offset b.sub.r is also needed, the pseudorange PR(t1;i;j) is also measured or calculated, and Eq. (25) is applied to determine the offset value b.sub.r.

In a third embodiment of the invention, pseudorange PR(t;i;j) is measured, preferably for at least two satellites, j=j1 and j=j2, and Eq. (25) is rewritten for each of these two satellites in the forms

Equation (58) is subtracted from Eq. (59), and the result is rearranged to yield ##EQU9##

Equation (60) is rewritten in terms of the unknown quantities x.sub.r0, y.sub.r0 and z.sub.r0 as ##EQU10##

At least two of the following six quantities are also measured: the zeroth, first and second time derivatives of Doppler shift for satellite j=j1;

and the zeroth, first and second time derivatives of the Doppler shift for j=j2.

From Eqs. (26) and (46), the zeroth derivative of Doppler shift for a signal received from satellite j may be rearranged to

where the quantities m.sub.x (t), m.sub.y (t), m.sub.z (t) are set forth in Eqs. (27)-(29). From Eq. (39), a first selected combination of zeroth and first time derivatives of the Doppler shift may be expressed as

From Eqs. (40), a second selected combination of zeroth, first and second time derivatives of the Doppler shift may be expressed as

Each of Eqs. (68), (69) and (70) can be expressed in the general form

where the index k=0, 1 and 2 corresponds, respectively, to the zeroth time derivative of Doppler shift, the first selected combination of zeroth and first time derivatives of Doppler shift, and the second selected combination of zeroth, first and second time derivatives of Doppler shift, respectively, for a signal received at the receiver (i) from the satellite (j). Because Doppler shift measurements from two satellites (j=j1 and j=j2) are available here, Eq. (71) represents any of six Doppler shift measurement combinations: three for the satellite j=j1 and three for the satellite j=j2.

Two of the six relations (71) are chosen for index values k=k1 and k=k2 and for the satellite indices j=j1 and j=j2, respectively, where j1 may be the same as, or may be different from, j2. If j1=j2, the indices k1 and k2 must differ; if j1.noteq.j2, the index k1 may be the same as, or may be different from, the index k2. These two relations are rewritten in the form ##EQU11## and Eq. (72) is inverted to express x.sub.r0 and y.sub.r0 in terms of the unknown z.sub.r0 :

For simplicity in notation, Eqs. (73) and (74) are rewritten as

with the coefficients C.sub.x, D.sub.x, C.sub.y and D.sub.y thus defined.

Equations (76) and (77) are now used in Eq. (66) to express the unknown z.sub.r0 in terms of the remaining unknown, b.sub.r, ##EQU12## For simplicity of notation, Eq. (77) is rewritten as

Equations (76), (77) and (79) are then substituted in Eq. (25) to produce a quadratic equation ##EQU13## for the remaining unknown, the clock offset b.sub.r, where t'=t-t0 and t"=t-t(j). The quadratic relation (80) is then solved for the clock offset b.sub.r ; only the more realistic value of b.sub.r, close to zero, is accepted here. This value of b.sub.r is used in Eq. (79) to determine the value of z.sub.r0, and the value for z.sub.r0 is used in Eqs. (76) and (77) to obtain the values for x.sub.r0 and y.sub.r0. Note that the value for clock offset b.sub.r is, or can be, determined in this approach for the third embodiment, whether or not this value is of any intrinsic interest. In other approaches for this third embodiment, a quadratic equation can be developed in terms of one of the other unknowns, x.sub.r0, y.sub.r0 or z.sub.r0, using a similar development, and the value of the clock offset b.sub.r may not need to be computed explicitly.

To recapitulate, this third embodiment uses two pseudorange measurements, for satellites j=j1 and j=j2, and uses two of the following six combinations of measurements: (1) the zeroth derivative of Doppler shift for satellite j=j1 and/or for satellite j=j2; (2) a first selected combination of the zeroth and first time derivatives of Doppler shift for satellite j=j1 and/or for satellite j=j2; and (3) a second selected combination of the zeroth, first and second time derivatives of Doppler shift for satellite j=j1 and/or for satellite j=j2, where the index values j1 and j2 are chosen independently. An alternate version uses one pseudorange value to obtain b.sub.r (optional) and three of the six combinations to obtain x.sub.r0, y.sub.r0 and z.sub.r0.

Thus far, the development has focused primarily on measurement and analysis of Doppler shift and its time derivatives. This development is extended to pseudorange and its time derivatives by the following observations. Pseudorange PR, after removal of error terms such as satellite clock offset, satellite noise, receiver noise, ionospheric propagation time delay, tropospheric propagation time delay and multipath errors, is approximately equal to the magnitude of the radius vector difference,

The first time derivative of pseudorange becomes

which is the Doppler shift, apart from a constant of proportionality, -c/.nu..sub.0.

Thus, measurement of the zeroth, first and second time derivatives of Doppler shift is substantially the same as measurement of the first, second and third time derivatives, respectively, of pseudorange. One can thus work with multiple time derivatives of pseudorange and/or with multiple time derivatives of Doppler shift (including the zeroth derivative) within the same formalism. For ease of reference, pseudorange measurement values and Doppler shift measurement values from a selected satellite to an LD signal receiver are collectively referred to as "location signal measurement" (LSM) values.

FIG. 2 is a flow chart illustrating the first embodiment of the invention. In step 21, at least one LSM value (pseudorange or Doppler shift) is determined for three non-collinear satellites. In step 23 (optional), LSM values for the three satellites (e.g., entries of the matrix M(t1;i;1;2;3) in Eq. (42)) are determined. In step 25, the initial location coordinates x.sub.r0, y.sub.r0 and z.sub.r0 are determined, using Eq. (45) or equivalent. In step 27 (optional), the receiver time offset b.sub.r is determined, using the pseudorange measurement in Eq. (25) or some other approach. In step 29 (optional), the present location coordinates x.sub.r (t), y.sub.r (t) and z.sub.r (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations.

FIG. 3 is a flow chart illustrating the second embodiment. In step 31, three or four LSM and/or LSM derivative values (pseudorange, zeroth, first and second time derivatives of Doppler shift) are determined for a single satellite. In step 33, LSM and derivative values for the entries of the matrix M'(t1;i;j) in Eq. (46) are used to determine the initial location coordinates x.sub.r0, y.sub.r0 and z.sub.r0, using Eqs. (46) or (57). In step 35 (optional), the receiver time offset b.sub.r is determined, using Eq. (25) or some other approach. In step 37 (optional), the present location coordinates x.sub.r (t), y.sub.r (t) and z.sub.r (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations.

FIG. 4 is a flow chart illustrating the third embodiment of the invention. In step 41, a first LSM value (preferably pseudorange) for each of two satellites (j=j1 and j=j2) is measured. In step 43, a first linear relation (corresponding to Eq. (66)) is obtained between the four unknowns x.sub.r0, y.sub.r0, z.sub.r0, and b.sub.r. In step 45, third and fourth LSM values (two of the six Doppler shift time derivative combinations) are measured, and second and third linear relations are obtained between three of the unknowns x.sub.r0, y.sub.r0 and z.sub.r0. In step 47, a quadratic relation involving only one of the four unknowns x.sub.r0, y.sub.r0, z.sub.r0, and b.sub.r is obtained, using the first, second and third linear relations, and a (realistic) solution is determined. In step 49, the remaining initial receiver location fix coordinates, x.sub.r0, y.sub.r0, z.sub.r0, and b.sub.r, (or of three location fix coordinates, if the value of b.sub.r or of one of the other location fix coordinates is not of interest) are computed. In step 51 (optional), the present location coordinates x.sub.r (t), y.sub.r (t) and z.sub.r (t) are determined, using Eqs. (1), (2) and (3) or some other suitable equations

For any embodiment, if three or fewer receiver location fix coordinates are to be estimated, for example, because the remaining coordinates are known, the total number of zeroth, first and second time derivatives of Doppler shift is reduced correspondingly. For example, if only two location coordinates are estimated, only two of the three Doppler shift derivatives need be measured from a single satellite, or from one Doppler shift derivative from each of two satellites. Measurement of pseudorange from one satellite plus Doppler shift(s) allows determination of clock offset. If pseudoranges are measured from two or more satellites, the number of Doppler shift derivative measurements is reduced by two, as shown in the preceding development.

The time derivatives for an LSM value can be estimated by at least two approaches. In a first, straightforward approach, LSM values L(t;i,j) are measured from signals received at a receiver (i) from a selected LD signal source (j) at a sequence of times t=.tau.1, .tau.2, .tau.3, .tau.4, . . . that satisfy .tau.1<.tau.2.ltoreq..tau.3<.tau.4< . . . A first time derivative of the LSM value at t=.tau.2 is defined as ##EQU14## where the dependence upon receiver index (i) and upon LD signal source index (j) are suppressed here to simplify the notation. A second time derivative of LSM value at a time t=.tau.2'=(.tau.2+t3)/2 is defined as ##EQU15## Users of this first approach to determining LSM value time derivatives may experience difficulties with reduced signal-to-noise (S/N) ratios because of the definitions adopted for the derivatives in Eqs. (83) and (84).

In a second approach, the S/N ratio is not reduced as severely, but the required computation is increased. In a first variation of the second approach, the LSM value is measured at several distinct values t=.tau.k(k=1, 2, 3, 4,. . . ), with .tau.1<.tau.2<.tau.3<.tau.4<. . . , and LSM values for intermediate values (.tau.1<t<.tau.4) are defined by interpolation as

and the time derivatives at t=.tau.2 become

Polynomials of degree higher than the second, third or fourth degree polynomials used in Eqs. (85) and (92) can also be used to determine the first and second derivatives of the LSM values, by an obvious extension of the development in Eqs. (85)-(90).

In a third approach, M not-necessarily-consecutive values (M.gtoreq.3) of LSM values L(.tau.m) are measured and a polynomial in t of degree M'(.ltoreq.M-2) is chosen to pass as closely as possible, according to a selected criterion, to each of the measured LSM values. For example, the polynomial may be ##EQU17## and the criterion may be least squares minimization of the collective error associated with the M measured values, by suitable choices of the polynomial coefficients p.sub.k (k=0, 1, . . . , M'): ##EQU18## Equations (93) provide a linear system of M'+1 relations in the polynomial coefficients p.sub.j, which can be solved by standard linear algebra techniques. The first and second derivatives of LSM values are then estimated as ##EQU19## which can be evaluated for any selected time t (with arbitrary time origin).

Whether the first approach, the second approach, the third approach, or some other reasonable approach is used here to estimate the first and second derivatives of LSM values from (only) the measured LSM values, the estimates for these derivatives will each be a linear combination of a plurality of measured values L(.tau.k) (k=1, 2, 3, 4, . . . ). At least two such measured values and at least three measured values of LSM values are required to estimate the first derivative and the second derivative, respectively.

Another method for obtaining equations relating LSM values and their time derivatives to (stationary) receiver location coordinates makes use of the circular rotation of the receiver about the Earth's rotational axis and of the circular rotation of the satellite (LD signal source) about the Earth. If r.sub.r (t;i) and r.sub.s (t;j) represent the present radius vectors (relative to the Earth's center or another common origin) of the receiver (i) and the satellite (j), respectively, the time derivatives of these quantities may be expressed as

where .omega..sub.r or and .omega..sub.s are the respective angular velocity vectors (constant, with dimensions sec-.sup.-1) associated with the rotation of the (stationary) receiver about the Earth's axis and with the rotation of the satellite about the Earth. Using the notations introduced in Eqs. (25), (30), (33) and (35), the zeroth, first and second Doppler shift derivatives are expressed as ##EQU20##

Each of the Eqs. (101)-(103) is a scalar equation, linear in the three unknown components (x.sub.r (t;i),y.sub.r (t;i),z.sub.r (t;i)) of the vector r.sub.r, with coefficients that are either constant or are known functions of t. The left hand sides of the Eqs. (101)-(103) are measurable combinations of the zeroth, first and second derivatives of Doppler shift (multiplied in some instances by the quantity .vertline..DELTA.r(t;i;j).vertline.). These three scalar equations are expressed in matrix form as

Here, the solutions are the present values (x.sub.r (t;i),y.sub.r (t;i),z.sub.r (t;i)) of the three location coordinates, not the "initial" values (x.sub.r (t0;i),y.sub.r (t0;i),z.sub.r (t0;i)) of these coordinates. Note that no restrictions are imposed on the angular velocity vectors .omega..sub.r and .omega..sub.s, except that each is referenced to the same origin, such as the Earth's center.

The orbit(s) of the satellite(s), assumed to be circular in Eqs. (13)-(15), may be extended to an elliptical orbit in the first embodiment, illustrated in FIG. 5, by the following development. The ellipse has a semi-major axis (x') of length a, a semi-minor axis (parallel to y') of length b, and an eccentricity .epsilon. given by

where b<a. The radius length .rho. of a line segment from an ellipse focus F0 to its intersection with the ellipse curve is given by

where .phi. is the angle the line segment makes with the y'-axis in FIG. 5. (K. R. Symon, Mechanics, Addison Wesley, Cambridge, Mass., 1953, pp. 111-114). Note that the azimuthal angle .phi. is defined with reference to the y'-axis rather than to the more conventional x'-axis. The analogues of Eqs. (13)-(15) for an elliptical orbit, instead of a circular orbit, are verified to be

where the "anchoring" body is located at an ellipse focus FO (FIG. 5) and .kappa. is a measure of the angular momentum of the orbiting satellite. Using a result found in G. Petit Bois, Tables of Indefinite Integrals, Dover Publs. New York, 1961, p. 118, Eq. (131) can be integrated to yield an implicit relation between .phi. and t for an elliptic orbit:

In some situations, it may not be possible to obtain sufficiently accurate measurements on signals from one or more satellites to provide values for pseudorange, Doppler shift and two time derivatives of Doppler shift (or, equivalently, of pseudorange and three time derivatives of pseudorange). For example, it may be possible only to obtain values of pseudorange and the first and second time derivatives of pseudorange for a satellite at the receiver. In this situation, receiver location can be estimated, using an approximation. Alternatively, this approximation may be used if one of the location fix coordinates, receiver clock offset b.sub.r, is known to be small and is not of much interest. In this approximation, one returns to Eq. (25) and re-expresses this as an approximation, ##EQU23## where presence of the receiver clock offset b.sub.r is ignored. The linear relations (73) are then used to express the quantities x.sub.r0 and y.sub.r0 in terms of z.sub.r0, and Eq. (133) is restated as a quadratic equation in the remaining unknown, here z.sub.r0. ##EQU24## The more realistic solution z.sub.r0 for Eq. (134) is then accepted and used in Eqs. (73)-(74) to determine the remaining initial location coordinates x.sub.r0 and y.sub.r0. Equations (1)-(3) can be used, as before, to determine the values of the x-, y- and z-coordinates of the receiver at some other selected time t.

FIG. 6 is a schematic view of apparatus suitable for practicing the invention and may be a console or be handheld. The apparatus includes an LD signal antenna 12 that receives and feeds LD signals to an LD signal receiver 11 for processing and determination of the LSM values. The LSM values are then passed to a computer 14, which may be part of the LD signal receiver 11 or may operate independently. The computer 14 includes a microprocessor 15 and associated memory module 16 that are programmed to perform the computations on the pseudorange values and Doppler shift time derivative values to determine the initial receiver location fix coordinate values (x.sub.r0,y.sub.r0 z.sub.r0,b.sub.r), receiver coordinate values (x.sub.r (t),y.sub.r (t),z.sub.r (t)) at some other selected time, elliptical orbit parameters, or any other relevant numerical quantity. The computer 14 optionally includes a data entry module 17, such as a keyboard or light pen, a visual and/or audible display 18 for the information computed or measured, and a power module to provide power for the other components shown in FIG. 6.