Navigation using carrier Doppler shift from a LEO constellation: TRANSIT on steroids

2.1 Doppler shift from the time derivative of accumulated delta range

The measured carrier Doppler shift equals the negative time derivative of the accumulated delta range divided by the carrier wavelength:

1

where 𝐷^𝑗 is the measured carrier Doppler shift from the 𝑗^th GNSS satellite in Hz, 𝜆 is the nominal wavelength in meters, is the accumulated delta range in meters, and 𝑡_𝑅 is the erroneous receiver clock time. It is related to the true signal reception time 𝑇_𝑅 and the receiver clock offset 𝛿_𝑅 as follows:

2

The accumulated delta range can be modeled as follows:

3

where is the unknown user receiver position vector in Earth-Centered/Earth-Fixed (ECEF) Cartesian coordinates–typically WGS-84 coordinates, is the known orbital position time history of the 𝑗^th GNSS satellite in ECEF coordinates, 𝜔_𝐸 is the Earth’s rotation rate, is the propagation time of the signal from the satellite to the receiver, 𝛿^𝑗 is the satellite clock offset, is the signal delay due to the neutral atmosphere (mostly due to the troposphere), is the time-equivalent carrier phase advance that is caused by dispersive refraction in the ionosphere, and 𝛽 is the bias on the measured beat carrier phase.

Although not explicitly noted in Equation (3), the signal propagation delay implicitly depends on , 𝛿_𝑅, and 𝑡_𝑅 because it must be computed by solving the following delay equation:

4

Note that the ionosphere term acts as a delay in this equation.

The troposphere delay term and the ionosphere advance/delay term also depend on , 𝛿_𝑅, and 𝑡_𝑅. This dependence comes partly through the Line-of-Sight (LOS) unit direction vector from the satellite to the receiver:

5

This unit vector affects the elevation angle through the troposphere and the ionosphere pierce-point location and zenith angle. The receiver location also affects the ionosphere pierce-point location directly.

The 3×3 direction cosines matrix:

6

appears in Equations (3) through (5). It compensates for the rotation of the ECEF coordinate frame while the signal travels from the satellite to the receiver.

2.2 Analytic carrier Doppler shift model

Equations (1) through (6) can be used to derive an analytic carrier Doppler shift model. If one ignores the time derivatives of the troposphere and ionosphere terms, then this model takes the form:

7

after multiplication of Equation (1) on both sides by −𝜆, which transforms to the range-rate equivalent of the carrier Doppler shift. A derivation of this equation is given in the appendix.

The new quantities introduced in this model are the unknown receiver velocity , the unknown receiver clock offset rate 𝑑𝛿_𝑅/𝑑𝑇_𝑅, the known ECEF satellite velocity at the time of transmission (but transformed to the ECEF frame of the time of reception) , which is the time of signal transmission, the satellite clock offset rate , and , which is a term that arises due to the time rate of change of the propagation delay .

The formula for the satellite velocity is:

8

The term takes the form:

9

where is the Earth’s rotation rate vector given in ECEF coordinates.

The carrier Doppler shift model in Equation (7) has some terms in common with the model given in (Aumayer & Petovello, 2015). The present model, however, is more complete and explicit in its characterizations of the effects of clock offset rates, satellite velocity, and the time rate of change of .

The Doppler shift model in Equation (7) has been used extensively by the author to solve for velocity and range-rate-equivalent receiver clock offset rate 𝑐(𝑑𝛿_𝑅/𝑑𝑇_𝑅) by using the carrier Doppler shift measurements from four or more GPS satellites. This velocity solution assumes that the receiver position and the receiver clock offset 𝛿_𝑅 have already been determined from a pseudorange solution. The values of and 𝛿_𝑅 are needed in order to compute the terms , and 𝑑𝛿^𝑗/𝑑𝑇^𝑗 that appear in Equation (7). Peraxis accuracies on the order of 0.01 m/sec RMS have been demonstrated. Velocity biases are negligible.

The accuracy is good enough to enable mapping of a vehicle’s relative route via simple quadrature integration of the velocity solution time history that results from using this model. Relative position accuracies on the order of 1 meter or better have been observed after 5 to 10 minutes of velocity integration. These levels of accuracy were obtained by the author in preparation of a university course and by students who completed the data-processing lab portion of that course using recorded carrier Doppler shift and pseudorange data from a survey-grade Magellan receiver mounted on an automobile. The foregoing statements represent the first public release of these accuracy results.

The dependence of the terms , and on position and receiver clock offset is a liability when using carrier Doppler shift measurements for velocity estimation. A pseudorange-based navigation solution is needed to overcome this liability.

In the present context, however, this dependence is an asset. It makes position and receiver clock offset simultaneously observable with velocity and clock offset rate. The minimum number of needed signals grows from four to eight in this situation because of the increased number of unknowns that must be estimated. As will be shown in the section on generalized GDOP analysis, this method’s accuracy depends on the rapidity with which the GNSS satellites pass over the receiver. That is why this technique will not work well with MEO GNSS satellites even when carrier Doppler shift data are available from eight or more of them.

2.3 Carrier Doppler shift model based on finite differencing of accumulated delta range

Although the Doppler shift model in Equation (7) works well for velocity and clock offset rate estimation when position and clock offset are known, it is not clear that it will work well for the present application. It’s performance may be lacking due to its neglect of the time derivatives of the troposphere and ionosphere terms.

Therefore, an alternate carrier Doppler shift model has been developed for use in the joint estimation of position, receiver clock offset, velocity shift, and receiver clock offset rate. One could develop an analytic model that included the effects of the non-zero time derivatives of the troposphere and ionosphere terms. The complexity of such a derivation, however, would be considerable. Instead, a suitable model has been developed based on a 5-point finite difference calculation that numerically computes the time derivative of the accumulated delta range.

This model starts with the assumption that there already exists a function which computes the model of the accumulated delta range given in Equation (3). The accumulated delta range model’s calculations are assumed to include solution of propagation delay Equation (4) for and model-based calculation of the troposphere and ionosphere terms and . Thus, it captures all of the ways in which depends on , 𝛿_𝑅, and 𝑡_𝑅.

The 5-point finite difference formula for the accumulated delta range rate uses the accumulated delta range model to perform the following calculation:

10

where Δ𝑡_𝑅 is the nominal finite difference interval of the erroneous receiver clock time 𝑡_𝑅 and:

11

is shorthand for the unknown receiver clock offset rate. The true reception time’s corresponding finite difference interval is , which accounts for the presence of the factor in several parts of this formula. Note that this 5-point finite-difference formula only contains four weighted terms because it uses a weighting of zero for the middle term, (i.e., the term with zero time offset).

The accumulated delta range rate model in Equation (10) is substituted into Equation (1) to yield the carrier Doppler shift model that forms the basis of point solutions for position, clock offset, velocity, and clock offset rate. This finite-difference model has been employed successfully by the author in a navigation system that uses carrier Doppler shift data from the Iridium constellation in addition to data from other sensors. A finite difference interval in the range 0.1 sec ≤ Δ𝑡_𝑅 ≤ 0.25 sec works well.

2.4 Simplified analytic carrier Doppler shift model for use in GDOP analysis

The preceding Doppler shift models are too complicated for use in a Geometric Dilution of Precision (GDOP) analysis. A GDOP analysis only needs to consider the dominant ways in which the unknowns affect the measurements. Its model must characterize how errors in the measurements are related to errors in the estimates of the unknown quantities. Second-order effects, while important in the calculation of an actual navigation solution, typically have a negligible impact on a GDOP calculation.

Several small effects have been neglected in order to simplify the carrier Doppler shift model for use in the GDOP analysis. First, the effects of the troposphere and the ionosphere have been omitted. They do not have a large impact on the sensitivity of the measurements to small perturbations of the unknowns. Omission of the atmospheric terms allows Equation (7) to be the starting point for developing a simplified model.

Three additional simplifications substitute the value 1 in place of three terms in Equation (7) that nearly equal 1. The three substitutions are:

12

The first two of these approximations are typically accurate to seven or more significant digits due to the stability of the transmitter and receiver clocks. The last of these approximations is accurate to more than four significant digits due to the low speed of the satellite relative to the speed of light.

Given the foregoing approximations, the simplified carrier Doppler shift model takes the form:

13

This is the model that will be used to develop a generalized GDOP analysis. It is suitable only for that purpose. It should not be employed to compute point solutions using this paper’s new navigation method. It should not even be used for determining velocity and receiver clock offset rate when position and clock offset are known from a pseudorange-based navigation solution.

5.1 Truth-model simulation parameters and cases

The nonlinear least-squares batch filter has been tested using truth-model simulation data from two sets of cases. The cases all assume the 2825-satellite final Starlink constellation that is listed in the final line of Table 1. The satellites are assumed to yield a valid carrier Doppler shift measurement if they lie above a 7.5^o elevation mask at the truth receiver location. This mask angle yields between 89 and 189 visible satellites for the set of cases that have been considered. The batch filter’s models for the troposphere and ionosphere terms have been assumed to match the truth troposphere and ionosphere terms.

The range-rate-equivalent measurement error standard deviation is assumed to be m/sec for all visible satellites. This level of range-rate error is commonly achieved for GPS signals. Whether it will be achievable for any given LEO constellation depends on the received signal’s carrier-to-noise ratio, form of the frequency discriminator, length of the interval used to compute accumulations for input to the carrier frequency discriminator, and PLL or FLL bandwidth (if a PLL or FLL is applied to a sequence of discriminator outputs).

Truth locations have been selected randomly over the surface of the Earth with equal probability density over the entire globe. The truth altitudes have been selected randomly from a flat distribution in the range 0 to 9,144 meters above sea level (30,000 ft). The truth receiver clock offsets have been sampled from a flat distribution between −0.25 sec and +0.25 sec. The truth velocities’ individual ECEF component have been sampled randomly from a Gaussian distribution with a mean of 0 and a standard deviation of 137 m/sec. This yields a 10% probability that the truth velocity magnitude is above the speed of sound. The truth receiver clock offset rates have been sampled from a Gaussian distribution with a mean value of 0 and a standard deviation of 3.336 × 10⁻⁹ seconds/second, which yields a standard deviation of 1 m/sec for the truth range-rate-equivalent receiver clock offset rates .

Two sets of 100 cases have been considered. For one set of cases, the satellite ephemerides and transmitter clock frequency offsets 𝑑𝛿^𝑗/𝑑𝑇^𝑗 are assumed to be perfectly known by the batch filter. The other set of cases assumes that there are knowledge errors in these quantities on the part of the batch filter.

The filter’s satellite position vectors for 𝑗 = 1,…,𝑁 are assumed to have per-axis knowledge errors that are sampled from a Gaussian distribution with a mean of 0 and a standard deviation of 2 meters. The filter’s satellite velocity vectors for 𝑗 = 1,…,𝑁 are assumed to have per-axis knowledge errors that are sampled from a zero-mean Gaussian with standard deviation equal to 0.002 meters/sec. The filter’s satellite clock frequency offsets 𝑑𝛿^𝑗/𝑑𝑇^𝑗 for 𝑗 = 1,…,𝑁 are assumed to have knowledge errors that are sampled from a zero-mean Gaussian distribution with a standard of 3.3 × 10⁻¹¹ seconds/second.

These levels of orbit and clock accuracy can be achieved if the satellites use an onboard GNSS receiver for real-time navigation (Montenbruck et al., 2009). If the system must operate independently of other GNSS, then achieving these levels of orbit and clock frequency accuracy could be a challenge.

5.2 Batch filter initialization

The 200 considered cases have been initialized with errors that are sized to provide significant tests of the filter’s ability to converge during normal operation. The initial position errors range from 143 to 151 km in the horizontal direction. The filter has been initialized with an altitude of 0 meters, a receiver clock offset of 0 seconds, and a receiver clock offset rate of 0 seconds/second in all 200 cases.

Thus, the initial altitude errors form a flat distribution in the range −9,144 m to 0 m. The initial clock offset errors form a flat distribution in the range ±0.25 sec. The initial per-axis velocity errors form a zero-mean Gaussian distribution with a standard deviation of 137 m/sec. The initial clock offset rate errors form a zero-mean Gaussian distribution with a standard deviation of 3.336 × 10⁻⁹ seconds/second.

In a true cold start, the initial horizontal position error could be much larger than 143 to 151 km. It might be as large as 20,000 km. Future work should examine the batch filter’s ability to converge from such large errors. If it has trouble, then there likely exist strategies that could be used to augment the filter in order to ensure convergence.

5.3 Batch filter performance results

The performance of the batch filter on the two sets of simulated cases is summarized in Table 2. The upper line of results applies to the 100 cases which assume that the filter has perfect knowledge of the satellites’ ephemerides and clock offset rates. The second line of results applies to the 100 cases that consider the effects of realistic knowledge errors between the truth and filter knowledge of the ephemerides and the transmitter clock offset rates. Note that the quantities , , and in some of the table’s column headers indicate the batch filter’s optimal estimates of , , and 𝛿_𝑅. The quantities , , and 𝛿_{𝑅𝑡𝑟𝑢𝑡ℎ} indicate their corresponding truth values from the simulations.

The first, second, and third results columns of Table 2 document the batch filter’s convergence performance. The filter successfully converged within 17 Gauss-Newton iterations for all of the perfect-ephemerides cases and within 14 iterations for all of the erroneous-ephemerides cases, as indicated by the first and third columns.

The average number of iterations to converge was between 8 and 9, as indicated by the second column. Actual convergence usually occurred in 4-5 iterations. The inflated iteration count statistics in Table 2 are the result of having used overly stringent termination criteria for purposes of deciding when to stop the Gauss-Newton iterations and declare that they have converged to a solution.

The fourth and fifth results columns of Table 2 indicate good position estimation performance. The RMS and peak position estimation error magnitudes are 1.35 and 4.16 m respectively in the perfect-ephemerides case. These error statistics increase modestly to 2.27 m RMS and 5.43 m peak when using the erroneous ephemerides. This level of absolute positioning accuracy is comparable with pseudorange-based MEO GNSS.

The sixth and seventh columns of Table 2 indicate good velocity estimation performance. The RMS accuracy is on the order of 0.01 m/sec, the peak error is only 0.0435 m/sec for the 100 cases with erroneous ephemerides, and the 100 cases with perfect ephemerides have better accuracy. This performance is similar to that of pseudorange-based MEO GNSS when using carrier Doppler shift only for purposes of determining velocity and receiver clock offset rate.

The eighth and ninth columns of Table 2 show that the receiver clock offset is the weak part of this system. The RMS and peak errors are on the order of a fraction of a msec, (i.e., in the range 0.0002 to 0.0009 sec). These error statistics are orders of magnitude worse than for pseudorange-based GNSS. As will be shown in the analysis of GDOP, the reason for this poorer timing accuracy is that the system’s ability to estimate time is tied directly to the maximum value of the acceleration of the distance between the receiver and the satellite. This maximum is not nearly large enough to support pseudorange-like timing accuracy.

Table 2 does not report the accuracy of the estimated receiver clock offset rate. The RMS and peak errors in the range-rate-equivalent receiver clock offset rate , are both smaller than 0.01 m/sec for both sets of 100 cases. This accuracy is better than that of the velocity estimates. This level of performance is comparable to that of pseudorange-based MEO GNSS when using carrier Doppler shift to estimate only this quantity and the velocity.

6.1 Linearized relationship between measurement errors and point-solution errors

The generalized GDOP analysis for this system uses a similar approach to the standard GDOP analysis of the pseudorange navigation solution. It starts by developing a linearized relationship between the errors in the carrier Doppler shift measurements and the errors in the point solution’s estimated quantities. This relationship is developed using the approximate carrier Doppler shift model in Equation (13). It is:

22

where Δ𝐷^𝑗 for 𝑗 = 1,…,𝑁 are the carrier Doppler shift measurement errors, is the error in the estimated ECEF position vector, Δ𝛿_𝑅 is the error in the estimated receiver clock offset, is the error in the estimated ECEF velocity, and is the error in the estimated receiver clock offset rate. The ECEF acceleration of the 𝑗^th satellite that is used in this formula is:

23

Two related challenges must be addressed before the model in Equation (22) can be used in a generalized GDOP analysis. The first challenge concerns the 𝑁×8 coefficient matrix on the right-hand side of this equation. Its first three columns have units of 1/sec, its fourth column has units of m/sec², its fifth through seventh columns are non-dimensional, and its eighth column has units of m/sec. A GDOP analysis needs to work with a coefficient matrix that is non-dimensional.

The second challenge is that the estimation errors in , , , and all have different units. The units of are meters, the units of are seconds, the units of are m/sec, and is non-dimensional. A GDOP analysis needs all of the estimation errors to have the same units. It also needs their common units and the units of the measurement errors to be identical.

Suppose one defines the range-rate-equivalent measurement error for the 𝑗^th satellite to be −𝜆Δ𝐷^𝑗. Then all 𝑁 of the range-rate-equivalent measurement errors have units of m/sec. The estimation error also has units of m/sec. It can be left in its current form. The estimation errors in , , and , on the other hand, need to be rescaled in order to have m/sec units.

The needed rescaling is straightforward in the case of . This receiver clock offset rate error can be rescaled through multiplication by 𝑐 in order to define the range-rate-equivalent receiver clock offset rate error . The selection of appropriate rescaling factors for and Δ𝛿_𝑅, on the other hand, requires careful thought.

6.2 Maximum Line-of-Sight (LOS) sweep rate and maximum range acceleration

A reasonable approach to rescaling the errors and Δ𝛿_𝑅 is to multiply them by upper bounds on the magnitudes of their corresponding columns in the 𝑁×8 coefficient matrix of Equation (22). The first three columns, the columns, contain the time rates of change of the unit direction vectors that point from the satellites to the receiver. The fourth column, the Δ𝛿_𝑅 column, contains the negatives of the accelerations of the scalar ranges from the satellites to the receiver that are caused by satellite motion. One can compute simple upper bounds on the magnitudes of these quantities. These bounds provide the needed scaling factors.

Consider the relative geometry between a satellite, its circular orbit, and a receiver, as depicted in Figure 1. The quantity 𝑎_𝑜𝑟𝑏 is the semi-major axis of the satellite’s orbit. The quantity 𝑅_𝐸 is the radius of the Earth. It is also the nominal distance from the center of the Earth to the receiver. The figure depicts the unit direction vector that points from the satellite to the receiver , and the scalar range from the satellite to the receiver 𝜌. The needed scaling parameters are an upper bound on and an upper bound on under the assumption of a stationary receiver.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1

Geometry of a satellite pass over a receiver that is relevant to computing maximum LOS sweep rate and maximum range acceleration

Upper bounds can be computed under the simplifying assumptions of a non-rotating Earth, a circular Keplerian satellite orbit, and a stationary receiver. In this case, the maxima of the two quantities occur if the orbit passes directly overhead of the receiver, and they both occur exactly at the time when the satellite is directly overhead. A straightforward analysis of this case can be used to show that:

24

25

The scaling factor 𝛾 has units of 1/sec (or, equivalently, radians/sec). The initial nondimensional factor in its formula is a function of the distance ratio 𝑅_𝐸/𝑎_𝑜𝑟𝑏. This factor increases as 𝑎_𝑜𝑟𝑏 approaches 𝑅_𝐸 because its denominator gets to be significantly smaller than 1.

The second term in the 𝛾 formula is the mean motion of the satellite orbit. It also increases as 𝑎_𝑜𝑟𝑏 decreases. The proposed system’s position accuracy tends to improve as 𝛾 increases, which will be discussed in the next subsection. Therefore, if all else is equal, then better performance can be achieved by making 𝑎_𝑜𝑟𝑏 smaller, (i.e., by using a LEO constellation for such a system).

The scaling factor 𝜂 has units of m/sec². The initial nondimensional factor in its formula is also a function of the distance ratio 𝑅_𝐸/𝑎_𝑜𝑟𝑏, and it also increases as 𝑎_𝑜𝑟𝑏 approaches 𝑅_𝐸, both because its denominator becomes much smaller than 1 and because its numerator increases. The second term in the 𝜂 formula is the 1/𝑟² gravitational acceleration at the satellite altitude. This factor also increases as 𝑎_𝑜𝑟𝑏 decreases.

The next subsection will discuss how the proposed system’s receiver clock offset accuracy tends to improve as 𝜂 increases. Again, if all else is equal, then better performance can be achieved by this system by making 𝑎_𝑜𝑟𝑏 smaller.

6.3 GDOP analysis using dimensional rescaling parameters

The rescaling parameters 𝛾 and 𝜂 can be used to redefine the position and receiver clock offset estimation errors so that they have units of m/sec. The position error in Equation (22) is replaced by . The units of 𝛾 are 1/sec. Therefore, the units of are m/sec. The receiver clock offset error Δ𝛿_𝑅 in Equation (22) is replaced by 𝜂Δ𝛿_𝑅. The units of 𝜂 are m/sec². Therefore, the units of 𝜂Δ𝛿_𝑅 are also m/sec.

These rescalings allow the linearized error relationship in Equation (22) to be rewritten in the following equivalent form:

26

which uses the following 𝑁×8 nondimensional coefficient matrix:

27

All of the elements of this matrix will have magnitudes on the order of 1. Its last four columns equal the entire corresponding matrix in the GDOP analysis of a pseudorange solution for position and clock offset. In the present analysis, they apply the velocity/clock-offset-rate part of the solution.

If the 𝑁 satellites have a multiplicity of semi-major axes, for 𝑗 = 1,…,𝑁, then a single value must be used to compute a single 𝛾 and a single 𝜂 for use in this GDOP analysis. One could use the average value or the minimum value in Equations (24) and (25).

This analysis concludes by computing its generalized nondimensional scalar GDOP value in much the same way as the standard GDOP analysis of a pseudorange navigation solution:

28

The scaling factors employed in this analysis can be used to develop rules for inferring the precision of the estimated quantities , , , and . These rules are listed in Table 3.

An examination of the scaling rules in the right-hand column shows that each precision has the correct units: meters for the precision of , seconds for the precision of 𝛿_𝑅, meters/sec for the precision of , and seconds/second for the precision of . Each quantity’s precision scales linearly with GDOP. Thus, it is always better to have a smaller GDOP, all else being equal. Each quantity’s precision also scales linearly with 𝜆𝜎_{𝐷𝑜𝑝𝑝}. Thus, it is always better to have a smaller range-rate-equivalent measurement error standard deviation. The precision of scales inversely with β, and the precision of 𝛿_𝑅 scales inversely with 𝜂. Therefore, it is better to have a large 𝛾 in order to get good position precision and a large 𝜂 in order to get good timing precision, all else being equal.

This analysis indicates that LEO orbits are better than higher orbits because they yield larger values of 𝛾 and 𝜂. Note, however, that the use of LEO orbits becomes practical only if a sufficient number of satellites are visible to each receiver so that GDOP will also be low.

Using MEO orbits would make 𝛾 and 𝜂 very small, and such a system would have very poor position and timing precision. The values of GDOP might be reasonably low for a MEO constellation, but the scaling parameters would be on the order of 𝛾 = 2×10⁻⁴ rad/sec and 𝜂 = 0.2 m/sec². A GDOP equal to 1 and a measurement error standard deviation of 𝜆𝜎_{𝐷𝑜𝑝𝑝} = 0.01 m/sec would yield a position precision of 1×0.01×(1/0.0002) = 50 m and a clock offset precision of 1×0.01×(1/0.2) = 0.05 sec.

Another significant result of this analysis is the importance of the geometric diversity and magnitudes of the normalized Line-of-Sight (LOS) rate vectors . In addition to geometric diversity for the LOS direction vectors themselves, , their time rates of change also must point in a variety of directions, and they must have sufficient magnitudes relative to the maximum possible sweep rates in order to achieve a low GDOP in this generalized analysis.

6.4 Validation of GDOP analysis using batch filter covariance

The foregoing GDOP analysis is valid if the two nondimensional matrices in the following formula are nearly equal:

29

where the 8×8 diagonal scaling matrix used in this formula is:

30

Recall that P_𝑥𝑥 is the batch filter’s approximate estimation error covariance matrix from Equation (21). In order for the GDOP derivation to be valid, near equality between the two matrices in Equation (29) must hold true if , a constant that does not depend on 𝑗, i.e., if all 𝑁 satellites have identical carrier Doppler shift measurement error standard deviations.

These two matrices have been compared for a representative example case. Their elements differ by, at most, a value equal to the largest element multiplied by 10⁻⁴. Furthermore, all eight of their eigenvalues agree to four significant digits or better as do all eight of their eigenvector directions. Therefore, the simplified carrier Doppler shift model of Equation (13) has yielded a GDOP analysis with good fidelity.

8.1 Limited downlink signal footprint of each satellite

A significant concern about the constellations as presently configured has to do with downlink antenna gain patterns. These constellations are communication systems. It does not make sense for neighboring satellites to have large overlaps in their downlink coverage areas. Large overlaps would waste transmission power and would cause unneeded challenges to the problem of servicing multiple independent data transmission requests within the system’s limited spectrum.

The specific beam patterns and overlap possibilities vary for different constellations. One constant among all the constellations is the following: It is unlikely that any one user will see strong main-lobe downlink signals simultaneously from eight or more satellites of a given constellation.

Therefore, the system envisioned will be impossible to implement without some additional strategies to see enough satellites. One strategy might be to receive signals from downlink antenna side lobes. Given the high data rates envisioned by these systems, the main lobes will have a lot of power. Therefore, the side lobes may have sufficient power to enable signal acquisition and determination of carrier Doppler shift.

In this case, however, it will likely be necessary to know a given signal’s modulated bit pattern, at least for short bursts of data. The modulated pattern would be used as a pseudo PRN code and would allow a long enough coherent accumulation interval to enable signal acquisition and determination of carrier Doppler shift.

A second strategy could be to accept the limitation of seeing few satellites simultaneously and try to see enough different satellites over a short time window to get a good navigation fix. The use of an INS might help to navigate using a time series of carrier Doppler shift data. This concept has been studied in (McLemore & Psiaki, 2020) and shows promise.

A third strategy could be to modify a constellation’s satellites to carry special navigation beacon antennas with wide gain patterns. It would be allowable to put the beacon signals in a different frequency band than the given satellite’s communication signals. Each satellite could broadcast a narrowband beacon signal because there would be no need for a wideband ranging code. It might be allowable to broadcast a beacon signal intermittently with a short duty cycle, perhaps one 0.1 sec burst every second.

8.2 High frequencies of downlink signals

Current plans for these large constellations call for downlink signal frequencies in the range 11 GHz–40 GHz. Signals at these frequencies can experience significant attenuation by foliage and significant scattering by rain. It would be better to have signals in a lower frequency range. If a constellation owner decided to add a beacon signal with a wide-gain-pattern antenna, it would be best if that signal were at a significantly lower frequency than the currently planned downlink communication signals.

8.3 Additional error sources

The present study has not considered the degradation of accuracy that might be caused by multipath errors. A thorough analysis of the impact of multipath on carrier-Doppler-shift-based navigation should be carried out for this concept. Given the good performance that the author has experienced with velocity estimation based on GPS carrier Doppler shift measurements, there is reason to believe that the impacts of multipath will not be too deleterious, but the issue needs study.

A great amount of research has been carried out on the subject of multipath-induced errors, but almost all of it considers the effects on pseudorange and beat carrier phase. Xu and Rife (2019) is the only work known to this author that specifically considers multipath effects on carrier Doppler shift. Its methods and results are not directly applicable to the question of how multipath would affect the present system. Therefore, a new type of multipath study will be needed to address this issue.

Another potential source of accuracy degradation is the mismatch between the true troposphere and ionosphere effects on carrier Doppler shift and the models for these effects that will be used by the batch filter. There is reason to hope that the impact of atmospheric mis-modeling will be small. Again, this hope comes from the good results for velocity estimation from carrier Doppler shift that have been achieved when completely ignoring atmospheric effects. Nevertheless, this issue must be carefully studied.

Most studies of troposphere and ionosphere effects consider only pseudorange and beat carrier phase. Graziani et al. (2009) examines a method for estimating and correcting troposphere effects of carrier Doppler shift measurements from deep space probes. (Klobuchar, 1996) briefly discusses ionosphere effects on carrier Doppler shift. However, neither of these works analyze the effects of troposphere or ionosphere modeling errors in a way that is directly relevant to the present system. To this author’s knowledge, no other published works address these questions in ways that would be relevant to this study. Therefore, research on these topics must develop new methods to address these concerns.

8.4 Signal design and signal processing for multiple access operation

The necessity to distinguish signals from 2,000 or more satellites, perhaps with 100 or more of them arriving simultaneously, raises issues of multiple access. If the satellites are to carry special beacons, then a way must be found to design signals so that they can all fit into a relatively narrow bandwidth and yet all be distinguishable. Therefore, effective multiple access strategies will be needed. The possibilities of code-division, frequency division, and time-division multiple access could be considered. Signal processing strategies for such multiple access signals would have to be developed in tandem with their signal structures.

If the existing signals are to be used, including side-lobe signals, then ways will have to be devised to distinguish among the many side-lobe signals that will be present in weak form. It is not clear whether this will be feasible. This problem is worth considering for the following reason: If it can be solved, then this navigation concept might be able to use a constellation such as OneWeb without requiring hardware modifications.

A solution to the side-lobe problem might involve sending special short data messages on the constellation’s downlinks that serve as PRN codes to enable multiple access and long coherent integration of weak side-lobe signals. If the constellation did not cooperate by broadcasting the needed data messages, then it might be possible to listen to main-lobe signals from a network of ground stations and provide after-the-fact information about modulated data. Such a system, however, would require a means to disseminate information about modulated data, which would involve high-bandwidth data broadcasts to user receivers. It would also introduce latency to their navigation solutions.

8.5 Batch filter convergence from large initial errors

The point-solution batch filter has been tested for its convergence capability from initial horizontal errors on the order of 150 km. A cold start capability would require assured convergence to the optimal solution from a first guess that is much farther away from the true position than 150 km. There likely exist good ways to ensure a large region of convergence of the batch filter, but additional work may need to be done to devise and implement a suitable strategy for ensured convergence.

8.6 Use of signals from multiple constellations

This study has considered signals that all come from one or the other of the several constellations that are being planned or fielded. It might make sense to use signals from multiple constellations. Such an approach could alleviate some of the GDOP problems that have been noted in the previous section.

McLemore and Psiaki (2020) have already considered this issue to a limited extent. Their proceeding paper examines possible use of the Iridium constellation in conjunction with either OneWeb, Starlink, or Kuiper. In a multi-constellation approach, it is not necessary that any one of the constituent constellations be able to provide full 4D navigation. Thus, the use of the existing Iridium constellation is reasonable within a multi-constellation system.

The usefulness of multiple constellations hinges on various factors such as signal ground footprints, the resulting impact on GDOP, and the practicality of designing and building receivers to cover the various frequency bands used by the constellations. These issues should all be studied in hopes of finding a multi-constellation solution that is practical.

8.7 Determining accurate satellite ephemerides and clock frequency offsets independent of existing GNSS

The present study has considered the impact of satellite ephemeris and clock frequency errors on navigation solution accuracy. The levels of errors considered, 2 m RMS per-axis position error, 0.002 m/sec RMS per-axis velocity error, and 3.3× 10⁻¹¹ seconds/sec clock frequency error, may be challenging to achieve for a stand-alone system that uses neither onboard GNSS receivers nor onboard atomic clocks.

Some sort of ground-control segment would have to be set up. It would operate a global network of receivers at accurately surveyed locations that have very stable clocks. They could estimate orbits and satellite transmitter clock frequency offsets based on the same carrier Doppler shift signals that support user equipment navigation. Such a system may be able to achieve the stated levels of accuracy if designed properly. The nearest existing equivalent based on on-way carrier Doppler shift is the DORIS orbit determination system (Jayles & Costes, 2004).

This study provides hope that the accuracies stated above might be achievable within a period when data are available. Whether they could be achieved over a needed forward prediction interval is an open question. The largest challenge might be clock frequency prediction.

It might be possible to address this challenge using constellation cross-links, if available, to distribute a frequency standard to the entire constellation that originates from radio-frequency links between a few satellites and a small number of ground stations. An alternative approach might be to link each satellite to a ground station often. Such an approach would reduce the needed prediction interval for clock frequency, but it would require a large network of ground stations.

It may be possible to achieve good user navigation accuracy with less accuracy of the ephemerides and transmitter clock frequencies than has been assumed in this paper. Perhaps the user equipment filter could consider states to allow for ephemeris and transmitter clock frequency errors.

The achievable accuracies of the ephemerides and the transmitter clock frequency offsets for an independent system and their impact on navigation error are worthy subjects for future research. These will be very important topics if independence from MEO GNSS is a requirement of such a system.