## Abstract

Utilizing broadband low Earth orbit satellite signals in an opportunistic manner for navigation is becoming increasingly popular. This paper deals with a particularly useful approach for navigation based on satellite signals of opportunity, which uses carrier Doppler-shift observables. We provide analytically derived and simplified formulas for the Jacobian involved in the numerical computation of the navigation solution and derive a global navigation satellite system-like dilution-of-precision metric that can be used to assess accuracy. A numerical study provides preliminary computational results.

## 1 INTRODUCTION

Using signals from low Earth orbit (LEO) satellites as signals of opportunity (SoOps) has recently been proposed as a resilient alternative to global navigation satellite systems (GNSSs), motivated by the fact that signals from LEO satellites have desirable characteristics for positioning, navigation, and timing (PNT) applications. LEO satellites are 25–40 times closer to Earth than their counterparts in GNSS constellations, which reside in medium Earth orbits (MEOs). Reduced orbit radii give rise to two important advantages from a user viewpoint: first, signals received on or close to the Earth’s surface are potentially significantly stronger; second, the orbital periods of LEO satellites are much shorter than those of MEO satellites. Stronger signals imply stronger immunity to jamming, particularly when combined with high spectral distribution. Smaller orbital radii result in higher angular velocities with respect to the receiver, thereby providing rich Doppler measurements that can be used for navigation.

Therefore, it is not surprising that interest in LEO satellite navigation has been increasing, fueled by companies such as OneWeb, SpaceX, Boeing, Samsung, Kepler, Telesat, and LeoSat, planning to launch thousands of LEO satellites in order to provide broadband internet and communication services globally. As a consequence, signals from LEO satellites will become abundant and diverse in frequency and direction. Indeed, SoOps that might be useful for navigation will come in different frequency bands, including the L-band, S-band, Ku-band, and more. The position and velocity of LEO satellites can be approximately calculated for a given epoch by using the two-line-element (TLE) files publicly available from the North American Aerospace Defense Command (Spiridonov et al., 2020) although precise applications may require epehemerides provided by the constellation operator or dedicated ground-control segment. Consequently, exploiting broadband LEO satellite signals in an opportunistic manner for navigation is becoming an increasingly popular approach that will become mature in the near future.

The current work revisits a particularly interesting approach for LEO-satellite SoOp navigation that uses carrier Doppler-shift observables, originally proposed by Psiaki (2021a). The main results include analytically derived and simplified formulas for the Jacobian involved in the numerical computation of the navigation solution, the specification of a GNSS-like dilution-of-precision (DOP) metric that can be used to assess accuracy, and a numerical study showing preliminary results. These developments expand upon the work of Psiaki (2021a), with the DOP formulations and numerical simulations developed via a newly formulated analytical Jacobian matrix. Additionally, the current work uses numerical simulations to examine navigation performance for constellations not previously studied.

This paper is organized as follows. Section 2 contains a review of existing literature. Although LEO navigation is relatively new, several alternative approaches have already been considered. Section 3 concentrates on the algorithm for using carrier Doppler-shift measurements and shows how computations can be improved by using explicit formulas for the Jacobian instead of numerical approximations. Section 4 shows how the DOP can be formulated for the LEO navigation case. Section 5 provides results of numerical simulations of the DOP formulations, the navigation algorithm performance, and a Monte-Carlo comparison of differing options for Jacobian matrices. Section 6 presents conclusions.

## 2 PREVIOUS WORK

The literature on LEO navigation is already quite extensive. Some reports, such as those by Reid et al. (2018), Reid et al. (2016), and Reid et al. (2020), as well as studies on the commercial Xona Ravens constellation, have proposed using special-purpose LEO satellites by tailoring the broadband protocol and infrastructure of the constellation to support precise PNT. Others have favored the opportunistic use of LEO satellite signals designed for communication to compute a navigation solution (Ardito et al., 2019; Khalife & Kassas, 2019; Morales et al., 2018; Pesyna et al., 2012). The latter approach is followed in this paper.

Multiple SoOp methodologies exist, including time of arrival (TOA), time difference of arrival (TDOA), angle of arrival (AOA), Doppler, and differential Doppler methods (Jafarnia-Jahromi et al., 2012; Neinavaie et al., 2022; Raquet, 2013). One-way TOA methods generally require precise transmitter clocks. LEO satellite clocks, unless disciplined by GNSS, do not possess the needed clock accuracy. Moreover, satellites in LEO constellations need not be tightly time-synchronized with each other (El-Mowafy et al., 2022), at least not to the level required for positioning trilateration. We note that range estimation is easier and more accurate when additional knowledge of the transmitted signal is available, whereas in opportunistic LEO measurements, the code of the transmitted signals may not be public knowledge. “Uncooperative” signals may be successfully decoded (see, e.g., work by Humphreys et al. (2022)) and code-phase ranging made possible, under the assumption that constellation operators will not change parts of the code in the foreseeable future.

Phase ranging and TOA require integer ambiguity resolution via either a more advanced least-squares method or the use of a base station (Subirana et al., 2013; Teunissen, 2017; Zhang et al., 2020). TDOA and differential Doppler methods rely on the use of a base station, whereas AOA methods require receivers with multiple-element antennas (Al Aziz, 2015). In contrast, Doppler methods do not require integer ambiguity resolution (in the absence of a cycle slip), the use of a base station, or multiple-element antennas. Therefore, Doppler-based positioning from LEO satellites appears to be an excellent alternative to GNSS.

Humphreys et al. (2022) explored pseudorange LEO navigation by decoding the broadcasted code from STARLINK and performing code-phase ranging. While not specifically intended for LEO, the work by Morales et al. (2016) explored the use of a pseudorange SoOp to aid an inertial navigation system (INS).

Carrier phase ranging from SoOps of LEO satellites was introduced in the work by Pesyna et al. (2012), in which a standalone carrier phase ambiguity resolution of a time division multiple access signal was attempted. Alternatively, Kassas (2020) explored the use of a base station to resolve the carrier phase ambiguity from a LEO SoOp.

The combination of SoOp pseudorange and Doppler-shift measurements has also been addressed via various methodologies. For instance, Levanon (1999) explored the use of range and range rate for positioning by intersecting a sphere of constant range, a cone of constant range rate, and the Earth’s surface. This method was successfully implemented by Satelles, Inc., in their navigation platform, which claimed 20-m accuracy for a static receiver.^{1}

Kassas (2020), Kassas et al. (2019), Morales et al. (2019), and Ardito et al. (2019) developed simultaneous tracking and navigation (STAN) with a LEO satellite signal framework, which simultaneously navigates and estimates the position of LEO satellites via the use of Doppler measurements and pseudoranges. In the STAN framework, the vehicle utilizes GNSS signals and uncertain information about LEO satellite states for navigation. The vehicle combines GNSS pseudorange mea-surements with the onboard INS while GNSS signals are available and tracks LEO satellites to refine estimates of their states. In the absence of GNSS signals, the vehicle switches to a mode wherein it simultaneously tracks LEO satellites and navigates by integrating SoOp pseudorange and Doppler measurements with its onboard INS.

The DOP associated with LEO SoOPs was initiated by McLemore & Psiaki (2022) under both Doppler shifts and pseudorange measurements. This issue was explored in order to improve the clock bias DOP (CBDOP) obtained when only Doppler shifts are used. Neinavaie et al. (2022) developed a method for using frequency difference-of-arrival SoOP signals from LEO satellites.

Although especially attractive for LEO navigation given the relatively large angular velocity between satellites and a receiver, Doppler shifts have also been proposed to enhance GNSS positioning. For instance, van Diggelen (2009) studied the use of Doppler shifts in an assisted Global Positioning System (A-GPS) context in order to jump-start an initial guess for GNSS positioning and hence reduce the time-to-first-fix under low signal-to-noise scenarios. Lehtinen (2002) and Progri et al. (2020) examined GNSS positioning with Doppler shifts for conditions in which code-pseudoranging was unavailable.

The use of Doppler navigation information as a SoOP from LEO satellites has also been explored. For instance, Khalife & Kassas (2019) examined the receiver architecture needed for Doppler information extraction and developed an extended Kalman filter for SoOp navigation from LEO Doppler shifts, whereas Benzerrouk et al. (2019) examined different estimation techniques for standalone Doppler navigation from the IRIDIUM-next constellation.

In this paper, we follow the alternative approach proposed by Psiaki (2021a), which uses a high-fidelity model of carrier Doppler shifts, to formulate a batch filter for simultaneously computing the eight unknowns of the underlying problem: three-dimensional position and velocity, clock bias, and clock drift.

## 3 LEO DOPPLER-SHIFT BATCH FILTER

GNSS and SoOp Doppler navigation are based on using variations of the Doppler-shift model (see, e.g., works by van Diggelen (2009), Lehtinen (2002), Khalife & Kassas (2019), Benzerrouk et al. (2019), and Morales-Ferre et al. (2020)). This model can be expressed as follows:

1

Here, *D ^{j}* denotes the Doppler shift from satellite

*j*,

*λ*is the wavelength, is the user-satellite unit line-of-sight (LOS) vector,

**v**

^{j}is the satellite velocity vector,

*c*is the speed of light, is the user clock drift calculated by differentiating the user clock bias

*δ*with respect to the true received signal time

_{R}*T*, is the satellite clock drift calculated by differentiating the satellite clock bias

_{R}*δ*with respect to signal transmission time

^{j}*T*, and

^{j}**v**is the user velocity. In this paper, unless otherwise mentioned, all position, velocity, and acceleration vectors are resolved in Earth-centered Earth-fixed (ECEF) coordinates. The standard use of this model (Angrisano et al., 2022; Petovello, 2015; Psiaki, 2021b) assumes that the position of the user is either known or has been estimated from pseudorange measurements, and the velocity is then computed. In contrast, in this paper, we will follow the approach of Psiaki (2021a), with the following eight equations in this section emanating from that reference. The current work utilizes the simultaneous dependence of the Doppler-shift observable on user position, velocity, clock bias, and clock bias rate. This approach is practical because of the large number of satellites that can be simultaneously observed in LEO constellations, together with the previously mentioned relatively large angular velocity with respect to the observer. Thus, we consider the eight-dimensional state vector:

2

where **r** and **v** are the respective user position and velocity vectors in ECEF coordinates and *δ _{R}* and are the user clock bias and drift, respectively. Now, let us consider the connection between the range-rate equivalent Doppler shift (RREDS) and the derivative of the accumulated delta range (ADR):

3

The ADR can be modeled by the following equation:

4

where **r**^{j} is the satellite position vector computed from ephemeris or TLE data, *ω _{E}* is the Earth’s rotation rate, is the propagation time from the satellite to the receiver, is the signal delay due to the troposphere, is the signal advance due to the ionosphere, and

*β*is the bias on the measured beat carrier phase. The time propagation is given as follows:

5

The matrix , which corrects the satellite’s position at the time of reception for Earth’s rotation during signal propagation, is given by the following:

6

By substituting Equation (4) into Equation (3) and neglecting the ionospheric and tropospheric terms, Psiaki (2021a) proposed a more accurate model:

7

where:

is the velocity of the satellite at the time of signal reception and:

is a term arising from the time rate of change of the propagation delay. **ω**_{E} is Earth’s rotational angular velocity vector in ECEF coordinates.

Because Equation (7) includes the effect of the satellite clock bias rate and propagation time on the relative LOS velocity vector, the model in Equation (7) is more accurate than the model in Equation (1), although Psiaki (2021a) noted that this model neglects the time derivatives of the tropospheric and ionospheric terms. These ionospheric and tropospheric error terms are functions of the receiver time-dependent user positions and clock biases, and therefore, they must be differentiated as such. As an alternative, Psiaki (2021a) proposed using a finite-difference approximation:

8

and then implementing a batch filter minimizing the performance index:

9

This approach is referred to herein as the Psiaki batch filter (PBF) approach. Because the cost function is nonlinear in the variables, the problem can be solved numerically, e.g., by using a recursive least-squares approach.

### 3.1 PBF Implementation

The PBF algorithm can be used to minimize Equation (9) by applying Algorithm 1. In Algorithm 1, a descent direction may be computed by using the Jacobian of the vector **f**(**x**_{k}). The most straightforward way to do this is by using perturbations to compute a numerical approximation, as shown in Equation (10). Although relatively easy to implement, a numerical approximation has two main disadvantages. First, a numerical approximation may become time-consuming, especially when a relatively large number of measurements are available. Second, a careful tuning of the perturbation parameters is required to compute a reasonable approximation:

10

A more sophisticated approach would be to compute an analytic formula for the Jacobian by direct differentiation. Although this procedure requires lengthy calculations, once this step has been completed, this approach becomes efficient from a computational viewpoint, and, moreover, is highly accurate. Our subsequent numerical studies also show this approach to be superior in terms of convergence to a solution. The interested reader is referred to Appendix A for a statement and derivation of the analytical Jacobian matrix.

In practical applications, one may trade-off accuracy for computational complexity, with a preference for computing a modified descent direction by using a simplified Doppler-shift model. For instance, the following approximation to the Jacobian can be used:

11

This Jacobian was derived by Psiaki (2021a) as a Jacobian matrix for the simplified Doppler-shift model given in Equation (1). Implementing this approximation is relatively straightforward when compared with the fully analytical formula or the numerical approximation. The derivative of the unit vector of geometric range between the satellite and receiver can be written as follows:

12

We note that in Equation (12), the time derivative of the propagation delay is not considered because, during the derivation of the analytical Jacobian matrix in Appendix A.1.1.1, it was found that the magnitude of is on the order of 10^{−9} s/m. Therefore, the time propagation is considered constant for the purposes of this derivation. The acceleration of the *j*-th satellite in ECEF coordinates can be obtained as follows:

13

where is the skew-symmetric matrix equivalent of the Earth angular velocity vector transformed from Earth-centered inertial to ECEF coordinates (Groves, 2013). Alternatively, more detailed analytical formulas for can be obtained from the TLE using formulas by Zhang et al. (2006) or via simple numerical differentiation. All three alternatives for , along with Equation (12) as an expression for , are possibilities for implementation of the simplified Jacobian matrix PBF.

## 4 DOP ANALYSIS

Given uncertain measurements, the DOP quantifies the expected accuracy of each computed variable as a function of the measurement variance. Assuming measurements with zero-mean Gaussian distributions with differing variances, the uncertainty on each PBF measurement can be modeled by *λσ _{dopp}*, i.e., the uncertainty in the RREDS measurements. A physically meaningful DOP can then be computed by scaling the output variables so that the matrix relating measurements to variables becomes non-dimensional (Psiaki, 2021a), thereby allowing an “apples and oranges” comparison of the different errors (e.g., position vs. velocity). To that end, the analytical expression for the Jacobian shown in Appendix A will be used. Given the state in Equation (2), the analytical Jacobian from Equation (A26) can be restated as follows:

14

Once the optimization algorithm has converged, the RREDS measurement noise and state errors are related via the following:

15

Now, let *S* be a diagonal matrix of scaling factors. We can then write the following:

16

The rules for a geometric DOP (GDOP) analysis are stated below.

A GDOP analysis requires that all elements of the Jacobian matrix be non-dimensional. Examination of Equation (14) leads to the following observations: the first three columns have units of ; the fourth column has units of ; the fifth through seventh columns are dimensionless; the eighth column has units of m/s.

A GDOP analysis requires that all estimation errors have the same units and that those units be identical to the units of the measurement errors. Examination of the error state in Equation (15) shows that Δ

**r**has units of meters, Δ*δ*has units of seconds, Δ_{R}**v**has units of , and has units of .

The above rules indicate that all error states must be scaled to have units of . Additionally, we must simultaneously divide the corresponding column in the Jacobian matrix in Equation (14) so that it is dimensionless. The following manipulations are identical to those performed by Psiaki (2021a).

### 4.1 Scaling of

The simplest example is the case in which , as it has units. Therefore, the last column of *H* in Equation (14) must be divided by *c*, and the error state must be multiplied by *c*.

### 4.2 Scaling of Δr

One approach for rescaling Δ**r** is to multiply Δ**r** by the upper bound on the corresponding column in the *H* matrix in Equation (14). If we examine the first three columns in *H*, we see that they are the time rate of change of the unit vector from the satellite to the receiver. The maximum of this value is called the maximum LOS sweep rate, which is given by the following (Psiaki, 2021a):

17

where *R _{E}* is the radius of the Earth,

*a*is the semimajor axis of the satellite’s orbit, and

_{orb}*μ*is Earth’s gravitational constant. Thus, we can divide the first three columns of

*H*by

*γ*and multiply the error Δ

**r**by

*γ*.

### 4.3 Scaling of *δ*_{R}

_{R}

The fourth column of *H* in Equation (14) is the negative acceleration in the direction of the unit vector between the satellite and the receiver. This value is referred to as the range acceleration, and its maximum is given as follows (Psiaki, 2021a):

18

### 4.4 GDOP Formulas

Consequently, we can synthesize the following relation:

19

where:

20

We can calculate the GDOP as follows:

21

When *N* satellites are used, a value must be chosen for *a _{orb}*; the average or minimum value can be chosen.

### 4.5 Additional DOP Definitions

In this subsection, we provide additional DOP definitions for individual components of the state. For the DOP associated with position and velocity, the values should be related to the expected accuracy in the north–east–down (NED) reference frame. Therefore, we slightly modify *A _{GDOP}* from Equation (20):

22

where is the rotation matrix from the ECEF reference frame to the NED reference frame. All DOP quantities depend on the following quantity:

23

The different DOP quantities are listed below.

Position DOP (PDOP):

24

Horizontal Position DOP (HDOP):

25

Clock Bias DOP (CBDOP):

26

Velocity DOP (VDOP):

27

Clock bias rate DOP (CBRDOP):

28

To calculate the expected output accuracy of a state element, the DOP value is utilized according to Table 1. Based on the relation in Table 1, there are several issues to consider with regard to the ability of obtaining a good expected accuracy of a state element (a smaller DOP is desired for better accuracy). Additionally, having lower noise for *λσ _{Dopp}* improves accuracy. Lastly, a large

*η*and a large

*γ*also improve accuracy. LEO satellites have much larger values for

*η*and

*γ*, and therefore, their use with the PBF yields a better expected accuracy relative to higher-orbit satellites.

## 5 NUMERICAL SIMULATIONS

In this section, we provide results from three simulations. The first is a DOP simulation. This simulation calculates the expected DOP and, thereby, the expected accuracy of the PBF for any constellation for any time and receiver location. The second simulation is a navigation simulation, which is also capable of propagating satellite constellations to a desired time and location. However, the navigation simulation uses the receiver time and location as well as satellite information to simulate RREDS measurements subject to noise. The simulation then uses the measurements along with the satellite orbits to perform PNT processing and obtain a navigation solution via the PBF. These first two simulations complement each other; the first calculates the expected accuracy using DOP formulations, whereas the second examines the actual PBF accuracy when RREDS measurements are subject to noise. Lastly, a Monte-Carlo performance comparison for different formulations of the PBF Jacobian matrices is performed.

### 5.1 DOP Simulation

In the DOP simulation, we simulated the expected DOP and accuracy for a stationary receiver located in Ramat Aviv, Israel, on June 14, 2022, at 14:59:41 for a duration of 10 min. Our code used a TLE file from the same day to propagate a multitude of constellations via the MATLAB satellite communications toolbox. Our simulation subsequently used the satellite positions and velocities and the receiver location as inputs to the *A _{GDOP}* matrix from Equation (22). Subsequently, DOP values were calculated via Equations (24)–(28). We obtained the expected accuracy using equations from Table 1. We simulated a RREDS measurement error of 0.1 m/s. Each constellation had its own simulation, and we obtained DOP values for every time point during the simulation for each state member. In Table 2, we show the results of the expected position accuracy for all analyzed constellations.

From Table 2, it is evident that, with a measurement noise level of 0.1 m/s, STARLINK has the best expected horizontal accuracy. GPS and GNSS constellations will have poorer accuracy because of their small *η* and *γ* values. LEO satellites have much larger values for *η* and *γ*, and therefore, their expected accuracy is better.

### 5.2 Navigation Simulation

For our navigation simulation, we used TLEs and the MATLAB satellite communications toolbox to propagate satellite constellations to a desired time, and from those satellites that were in view of a ground-truth location, we generated RREDS measurements subject to Gaussian noise using Equation (8). Our satellite data, from a multitude of constellation and receiver ground-truth data, were used as inputs into the ADR (Equation (4)). We also provided ground-truth clock biases *δ _{R}* and clock bias rates . It should be noted that atmospheric effects were not modeled in the generated measurements. We then finite-differenced the ADR values according to Equation (8) and added 0.1-m/s random Gaussian noise to the simulated measurements. When at least eight satellites were in view of the receiver, we had the receiver use Algorithm 1 to estimate its own location, velocity, and clock data. We then compared the estimated state to the ground truth used to generate measurements. We again simulated a receiver located in Ramat Aviv, Israel, on June 14, 2022, at 14:59:41 for a duration of 10 min. This simulation helps determine how well a LEO SoOp receiver would perform.

The results of the simulation are shown in Table 3, which shows that, again, STARLINK has the best performance. These results agree with the order of magnitude of the DOP results for the expected accuracy in position shown in Table 2. We define accuracy as the root mean square (RMS) of the magnitude of the error vector. Additionally, in the absence of real signals from LEO constellations, this simulation gives a reasonable approximation of the real expected accuracy of the PBF for any constellation at any time and receiver location.

### 5.3 Comparison of Jacobian Matrices

The numerical Jacobian in Equation (10), the analytical Jacobian in Equation (A26), and the simplified Jacobian in Equation (11) were tested using Monte-Carlo conditions as input into Algorithm 1 for the PBF for 100 trials according to Sections 5.3.1, 5.3.2, and 5.3.2. The purpose of this analysis was to investigate the feasibility of using alternative Jacobian matrices in the PBF.

#### 5.3.1 Truth Parameter Sampling

We selected truth locations according to the same distributions used by Psiaki (2021a). Truth locations were randomly selected across the Earth’s surface, with an equal probability distribution over the entire globe. Truth altitudes were chosen randomly from a flat distribution within the range of 0–9,144 m above sea level (30,000 ft). Truth receiver clock offsets were sampled from a flat distribution between −0.25 s and +0.25 s. The individual ECEF components of truth velocities were randomly sampled from a Gaussian distribution with a mean of 0 and a standard deviation of 137 m/s, which results in a 10% probability that the truth velocity magnitude exceeds the speed of sound. Truth receiver clock offset rates were sampled from a Gaussian distribution with a mean value of 0 and a standard deviation of 3.336 × 10^{−9} s/s, which corresponds to a standard deviation of 1 m/s for the truth range-rate-equivalent receiver clock offset rates (where represents the offset rate).

#### 5.3.2 Initial Solution Sampling

We selected initial solution guesses for the PBF according to the following distributions. ECEF initial position errors were sampled from a uniform distribution between 143 and 151 km. The filter was initialized with zero velocity for each ECEF component. The filter was initialized with a receiver clock offset of 0 s. The receiver was initialized with a receiver clock offset rate of 0 s/s for all 100 cases.

#### 5.3.3 Miscellaneous Conditions

This simulation was conducted with the following miscellaneous conditions. The elevation mask used in the simulation was set to 7.5°. To model the satellite orbits, the simulation utilized TLE data obtained from June 14, 2022. The specific time point covered by the simulation was set to 14:59:41. Lastly, because the purpose of this analysis is to compare the performance of the Jacobian matrices, a noise signal with 0.1-meter-per-second (mps) 1-*σ* was added to the generated RREDS measurements used in the simulation.

#### 5.3.4 Results

From Table 4, it is clear that the analytical Jacobian matrix given in Equation (A26) has the best performance in some categories. All of the Jacobian matrices had an equivalent number of unconverged trials, and the analytical and simplified Jacobian matrices required similar numbers of iterations to reach convergence. This analysis also shows that the simplified Jacobian matrix in Equation (11) can feasibly be used with the PBF in Algorithm 1, with only a small degradation in performance, in some categories, relative to the analytical Jacobian. In some cases, the simplified Jacobian matrix even outperformed the analytical Jacobian matrix. The results from Table 4 are significant because they show that there are no downsides to using the simplified Jacobian. There are no outliers in any of the simulation cases with a much greater number of required iterations, outright failure to converge, or convergence to a much less accurate solution compared with the other two Jacobians. We note that in both the navigation and DOP simulations detailed above, a measurement error of 0.1 m/s was chosen. It remains unclear whether RREDS measurements can be obtained with this level of accuracy from the main lobe or side lobe of LEO satellite transmissions at the specified 7.5° elevation mask.

### 6 CONCLUSIONS

This paper provides a comprehensive analysis of Doppler navigation with LEO satellites, presenting a useful approach for navigation using satellite SoOps and demonstrating a method for assessing accuracy using a GNSS-like DOP metric. Results on numerical simulations of the DOP formulations, the navigation algorithm performance, and a Monte-Carlo comparison of differing options for Jacobian matrices are presented, which provide valuable insights into the implementation and accuracy of this approach. This paper also provides a review of existing literature and discusses the algorithm for using carrier Doppler-shift measurements. Finally, the paper concludes that the simplified Jacobian is a viable option for implementing the PBF. This work provides closer insight for researchers and practitioners into Doppler navigation with LEO satellites.

## HOW TO CITE THIS ARTICLE

Baron, A., Gurfil, P., & Rotstein, H. (2024). Implementation and accuracy of Doppler navigation with LEO satellites. *NAVIGATION, 71*(2). https://doi.org/10.33012/navi.649

## APPENDIX A DERIVATION OF AN ANALYTICAL JACOBIAN

### A.1 Derivation

The Jacobian matrix is defined as the derivative of **f(x)** in Algorithm 1 with respect to each element of the state defined in Equation (2). In practice, because RREDS measurements do not depend on elements of the state, the Jacobian can be defined as the derivative of the measurement model defined in Equation (8) with respect to each state element. We denote the measurement model as follows:

A1

A2

A3

A4

The ADR is given in Equation (4).

#### A.1.1 Term 1

Substituting the finite differences from Term 1 of Equation (A1) into Equation (4), we have the following:

A5

Here, we define the following:

A6

and denote:

A7

##### A.1.1.1 Derivative of Term 1 with Respect to r

We differentiate Equation (A5) with respect to **r**:

A8

The term requires a numerical approximation. However, is on the order of 10^{−9} s/m whereas is on the order of 10^{5} m/s; thus, their product is on the order of 10^{−4}. Consequently, the term can be considered negligible. Moreover, . We have the following for the derivative:

T the *z*-coordinate of will be zero. Additionally, experimental tests show that in the *x* and *y* coordinates, the product is on the order 10^{2} m/s. When is multiplied by , which is on the order of 10^{−9} s/m, the product is on the order of 10^{−7}. This term is deemed negligible. With these simplifications, Equation (A8) reduces as follows:

A9

##### A.1.1.2 Derivative of Term 1 with Respect to δ_{R}

We start with the following:

A10

The following term:

must be evaluated numerically, but again, is on the order of 10^{−2} s/s whereas is on the order of 10^{−17} s/s; thus, their product is on the order of 10^{−19} s/s. Therefore, the term is negligible when compared with the −1 term. is on the order of 10^{−17} s/s and, when multiplied by , which is on the order of 10^{2} m/s, the product is considered negligible. Noting that is the velocity of the satellite at time *a*_{1} and using the simplifications above, we can rewrite Equation (A10) as follows:

A11

We realize that the dot product of the LOS unit vector and the satellite velocity is the component of RREDS arising from the satellite velocity. We denote this as follows:

A12

With this notation, Equation (A10) can be simplified as follows:

A13

##### A.1.1.3 Derivative of Term 1 with Respect to v

Here, we have the following:

A14

Following Psiaki (2021a), we assume here that .

##### A.1.1.4 Derivative of Term 1 with Respect to

Again, following Psiaki (2021a), we assume that and approximate as (**r** – 2**v**Δ*t _{R}*) in the first term in the unit vector below. Therefore, we can write the following:

A15

We can see experimentally that is on the order of 10^{5} m/s and . Their product leaves the left-hand term in Equation (A15) on the order of 10^{5} m. Thus, the right-hand term in Equation (A15) is clearly dominant. Therefore, we can approximate .

#### A.2 Jacobian Statement

The above process is repeated for Terms 2, 3, and 4, and a five-point finite difference of results is performed:

A16

Consequently, the analytical Jacobian matrix for the PBF is stated in MATLAB notation:

A17

A18

A19

A20

where:

A21

A22

A23

A24

and:

A25

In Equations (A16)–(A24), as mentioned by Psiaki (2021a), a finite-difference interval of 0.1 s ≤ Δ*t _{R}* ≤ 0.25 s can be applied. A compact version of the analytical Jacobian can be written as follows:

A26

## Abbreviations

- ADR
- accumulated delta range
- AOA
- angle of arrival
- CBDOP
- clock bias dilution of precision
- CBRDOP
- clock bias rate dilution of precision
- DOP
- dilution of precision
- ECEF
- Earth-centered Earth-fixed
- GDOP
- geometric dilution of precision
- GNSS
- global navigation satellite system
- GPS
- Global Positioning System
- HDOP
- horizontal position dilution of precision
- INS
- inertial navigation system
- LEO
- low Earth orbit
- LOS
- line of sight
- MEO
- medium Earth orbit
- mps
- meters-per-second
- NED
- north–east–down
- PDOP
- position dilution of precision
- PBF
- Psiaki batch filter
- PNT
- positioning, navigation, and timing
- RREDS
- range-rate equivalent Doppler shift
- SoOp
- signal of opportunity
- STAN
- simultaneous tracking and navigation
- TDOA
- time difference of arrival
- TLE
- two-line element
- TOA
- time of arrival
- VDOP
- velocity dilution of precision

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.