Doppler Positioning Using Multi-Constellation LEO Satellite Broadband Signals as Signals of Opportunity

2 DOPPLER-SHIFT POSITIONING MODELS

In the following representations, the unknown position and velocity vectors of the receiver (r) at the true time of reception (t = t_r – δt_r) are denoted as r_r(t) = [x_r, y_r, z_r]^T and vr(t)=[x˙r,y˙r,z˙r]T , where t_r is the erroneous receiver clock time and δt_r is the receiver clock offset. The state vector of the satellite at the time of burst transmission (t_r – δt_r – τ) includes the satellite position vector r^s(t_r – δt_r – τ) = [x^s, y^s, z^s]^T and velocity vector νs(tr−δtr−τ)=[x˙s,y˙s,z˙s]T , where τ is the satellite-to-receiver time of flight for each burst. Both the satellite and receiver state vectors are given in the Earth-centered Earth-fixed (ECEF) frame.

The general Doppler-shift observation (DrS) model can be derived from the derivatives of the accumulated delta range Δρ^s of satellite (s) with respect to the erroneous receiver clock time (t_r) as follows:

−λfDrs=dΔρsdtr=d(‖rr(t)−Rrs(tr−δtr−τ)‖+c(δtr−δts)+Trs−Ir,fs+εDrs)dtr 1

where λ_f is the transmitted carrier wavelength for frequency f and c is the speed of light. The satellite clock offset is denoted by δt^s. Tropospheric and ionospheric errors are denoted by T and I, respectively. εDrs represents the unmodeled errors and noise of the observations.

Denoting the angular velocity of the Earth with ω_e and assuming that the Earth’s rotational axis is aligned with the Z-axis, the rotation matrix R is defined as follows:

R=[cos(ωeτ)sin(ωeτ)0−sin(ωeτ)cos(ωeτ)0001] 2

This matrix applies the ECEF frame rotation due to the Earth’s rotation during the signal travel time and transforms the satellite state vector to the ECEF frame at the time of reception.

Defining the range between satellite and receiver as ρ = ||r_r(t) – Rr^s(t_r – δt_r – τ)|| and considering that the change in the satellite clock error is not solely a function of the satellite time but also adjusted by the signal travel time, the partial derivatives in Equation (1) are written as follows:

−λfDrs=(dρdt+c(dδtrdt−dδtsdts(1−dτdt))+dTrsdt−dIr,fsdt)dtdtr 3

The time derivative of the range in Equation (3) is calculated from the following equation:

dρdt=rr(t)−rs(t)‖rr(t)−rs(t)‖·(drr(t)dt−drs(t)dt)=e^rs(t)·(vr(t)−drs(t)dt) 4

where e^rS(t)=rr(t)−rs(t)‖rr(t)−rs(t)‖ is the line-of-sight vector from the satellite to the receiver.

The satellite position after the Earth rotation correction is r^s(t) = Rr^s(t_r – δt_r – τ), where its derivative is derived as follows:

drs(t)dt=(dRdτd(τ≅ρ/c)dt)·rs(t−τ)+Rdrs(t−τ)dt=(ωe[−sin(ωeτ)cos(ωeτ)0−cos(ωeτ)−sin(ωeτ)0001]1cdρdt)·rs(t−τ)+Rvs(t−τ)(1−1cdρdt) 5

After defining Ω_e = [0,0, ω_e]^T and performing some rearrangements, Equation (5) can be written in the following compact form:

drs(t)dt=(Ωe×Rrs(t−τ))1cdρdt+Rνs(t−τ)(1−1cdρdt) 6

Substituting Equation (6) into Equation (4) yields the following:

dρdt=e^rs(t)·(vr(t)−[(Ωe×Rrs(t−τ))1cdρdt+Rvs(t−τ)(1−1cdρdt)])=e^rs(t)·(vr(t)−Rvs(t−τ))−e^rs(t)·(Ωe×Rrs(t−τ)−Rvs(t−τ))(1cdρdt) 7

This range derivative is then calculated by factoring and applying rearrangement steps as follows:

dρdt=e^rs(t)·(vr(t)−Rvs(t−τ))1−1ce^rs(t)·(Rvs(t−τ)−(Ωe×Rrs(t−τ))) 8

where the numerator is the projection of the relative velocity between the receiver and the satellite onto the line-of-sight vector and the denominator accounts for the signal travel time correction, incorporating the satellite’s velocity and the Earth’s rotation effects.

Considering the satellite and receiver clock drift as δt˙s and δt˙r , defining the time derivatives of the ionosphere and troposphere as I˙r,fs and T˙rs , and substituting Equation (8) into Equation (3) with some rearrangements results in the following Doppler-shift observation model:

−λfDrs=[e^rs(t)·(vr(t)−Rvs(t−τ))(1+cδt˙s1−1ce^rs(t)·(Rvs(t−τ)−(Ωe×Rrs(t−τ))))+c(δt˙r−δt˙s)+T˙rs−I˙r,fs](1+δt˙r)−1 9

This compact form integrates all of the components necessary for computing the Doppler shift observed at the receiver, accounting for the motion and time variations of the receiver and satellite, environmental effects, and adjustments for clock drifts. Ignoring the atmospheric term results in a model similar to that developed by Psiaki (2021), who compensated for the ignored terms using a five-point stencil approach. The analytical implementation of his approach in a batch filter was presented by Baron et al. (2024). Instead of using the finite differencing approach, we develop multi-frequency, multi-constellation models with detailed elements provided in the following for both absolute and differential positioning, designed for real-time applications.

Each broadband constellation broadcasts various frequency ranges for its downlink transmission. For instance, Starlink satellites use a combination of Ku-band (10.7–12.7 GHz) and Ka-band (17.7–19.7 GHz) frequencies. These satellites employ advanced frequency reuse and beamforming techniques to optimize spectrum utilization and mitigate interference. Therefore, the received burst from each satellite could fall within specified frequency ranges. To accommodate these considerations in the model, the following vectorial representation is provided.

First, we define Dr=[Dr,1s,…,Dr,fs]T , where Dr,f=[−λfDr,f1,…,−λfDr,fm] for a constellation of 1 to m satellites transmitting bursts across various frequency ranges. This handles the transmission of different frequency ranges for each satellite in the model. We define e_f = [1,…, 1]_1×f and the following parameters:

ur=efT⊗[e^r1·(vr−Rv1)(1+cδt˙11−1ce^r1·(Rv1−(Ωe×Rr1))),…,e^rm·(vr−Rvm)(1+cδt˙m1−1ce^rm·(Rvm−(Ωe×Rrm)))]T(1+δt˙r)−1δt˙=efT⊗[c(δt˙r−δt˙1),…,c(δt˙r−δt˙m)]T(1+δt˙r)−1T˙r=efT⊗[T˙r1,…,T˙rm]T(1+δt˙r)−1I˙r=[I˙r,1T,…,I˙r,fT]TwhereI˙r,f=[I˙r,f1,…,I˙r,fm](1+δt˙r)−1εDr=efT⊗[εDr1,…,εDrm]T 10

The vectorial representation of the Doppler-shift observation model for absolute positioning, considering the entire constellations transmitting various frequency ranges, is then defined as follows:

E(Dr)=ur+δt˙+T˙r−I˙rQDr=diag(εDr) 11

Here, E represents the expected value, and QDr is the covariance matrix of the observations, which distinguish the deterministic and stochastic components of the model, respectively. ⊗ denotes the Kronecker product (Schott, 2016). Equation (11) is applicable for the multi-constellation case in which each satellite transmits bursts at various frequency ranges. The size of each term in this equation is mf × 1, and the size of QDr is mf × mf. The unknown parameters in absolute positioning via the developed model include the state vector, which consists of the user position, velocity, and clock drift (x(t)=[rr,vr,δt˙r]) . The estimation process requires linearization of the observation model through partial derivatives to construct the design matrix A. Based on the Doppler-shift observation model in Equation (9), a sequential estimation approach (either least-squares estimation or Kalman filtering) is employed, where the design matrix A represents the sensitivity of the observations with respect to the state parameters. The design matrix A is defined as the Jacobian matrix of the observation equation A=∂(−λfDrs)/∂x . The partial derivatives incorporate complex geometric relationships between receiver and satellite state vectors, accounting for their respective positions and velocities, as well as the Earth’s rotational effects. These relationships are available in the terms N_{i=0,…, 7} defined in the appendix and demonstrate how changes in each variable affect the observed Doppler shift. The derivation of these partial derivatives is presented in the following, establishing the mathematical foundation for the estimation process.

3 BROADBAND BURST SIMULATION

Satellites in broadband constellations deliver data to users in rapid, short bursts rather than in a continuous stream. This approach is primarily applied to enhance bandwidth utilization, reduce power consumption, and increase system flexibility in accommodating more users. This intermittent availability of data must be considered in positioning applications. Additionally, as demonstrated in the packet internet groper analyses of Starlink by Huston (2024), the Starlink signals fluctuate every 15 s because of satellite tracking and handover processes; thus, the user equipment tracks each satellite in 15-s intervals before possibly switching to another satellite in a complex handover process.

Another important factor is the restriction of the minimum elevation angle of the received signal for LEO broadband satellites (Pachler et al., 2021). For instance, optimal performance for Starlink is achieved for signals at elevation angles higher than 40°, but it has been suggested that this angle be decreased to 25° (Cakaj, 2021). Signals at higher elevation angles experience lower absolute Doppler shifts but higher Doppler rates. While this can simplify initial signal acquisition, the higher Doppler rates create additional challenges for continuous tracking. However, higher elevation angles still provide benefits in terms of signal quality (i.e., higher C/N₀) and reduced atmospheric signal delays.

Moreover, broadband constellations typically do not broadcast signals simultaneously to the same user from multiple satellites. Instead, they employ phased array antennas to dynamically steer their beams. As a satellite moves out of view of a user terminal, the connection is handed off to another satellite coming into view. These advanced techniques ensure continuous coverage, maintain a stable connection, and mitigate interference as satellites move across the sky. Each user terminal communicates with a specific satellite at any given time using assigned frequencies and time slots.

Therefore, to simulate realistic scenarios, the following assumptions are considered for Doppler-shift frequency simulation from the tested constellations to the user:

• 1) The sample interval for each burst is set to 15 s with a 5-s disruption, referred to in this paper as the “short-burst limitation.”

• 2) The elevation mask angle is set to 30°, referred to as the “elevation angle limitation.”

• 3) The receiver is capable of receiving signals from various constellations simultaneously.

• 4) The limitation applied to broadband LEO signals to avoid interference with broadband geostationary orbit services (Iannucci & Humphreys, 2020) is ignored. Therefore, no interruption to the signals due to this limitation is considered in the simulation.

The orbits of 6,336 satellites from the selected constellations (Starlink, OneWeb, Iridium) are propagated using the initialized simplified general perturbations (SGP4) model and approximated orbits estimated from two/three-line element (TLE or 3LE) files (Vallado & Crawford, 2008). Although satellite clock errors can be a significant source of error (El-Mowafy et al., 2023), they are not considered in the simulation because of a lack of detailed information about the oscillators of these commercial satellites. In real-world scenarios, satellite clock errors should be either estimated or modeled if sufficient information is available. Additionally, receiver clocks are not simulated because their errors are estimated epoch-wise (i.e., epoch-by-epoch) without the application of a specific receiver clock model, making the estimation independent of the simulated clocks. Orbital errors using TLE values can reach several kilometers, significantly impacting positioning accuracy. Onboard precise orbit determination (POD) could mitigate this issue, contingent upon precise corrections via space links (Allahvirdi-Zadeh et al., 2021; Allahvirdi-Zadeh, El-Mowafy, & Wang, 2024), power-efficient approaches (Wang et al., 2020; Allahvirdi-Zadeh & El-Mowafy, 2022, 2024; Allahvirdi-Zadeh, Wang et al., 2022), and developed software packages (Allahvirdi-Zadeh, El-Mowafy, McClusky et al., 2024). Figure 1 illustrates the selected constellations simulated for February 1, 2024. Each constellation transmits across a range of frequencies. In the simulation, random frequencies are selected for each satellite within the frequency range of each constellation. Refer to Table 1 for the specifications of these constellations.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1

Simulated constellations: Starlink (cyan), OneWeb (pink), Iridium (green)

Two users, one static and one kinematic, are included in the simulation. The static receiver, CUT000AUS (or briefly CUT0), is a continuously operating reference station (CORS) of the international GNSS service network at Curtin University, Western Australia. Ground truth for CUT0 is computed using long-baseline GNSS processing techniques through the Geoscience Australia AusPOS service (https://www.ga.gov.au/scientific-topics/positioning-navigation/geodesy/auspos). The kinematic receiver is a Trimble R12i placed on top of a car, which was driven around the Curtin University campus for approximately 10 min. The trajectory of the reference car positions (shown by red lines in Figure 2) is derived from post-processed kinematic solutions, with CUT0 serving as the reference station. The velocity of the car is obtained through numerical differentiation. Figure 2 illustrates these scenarios on a map.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

Kinematic test trajectory shown on Curtin campus map (left), static receiver: CUT0 (top right), and kinematic receiver set up (bottom right)

To gain insight into the atmospheric loss factors and their impact on positioning using SoPs, the simulation considers the following attenuations. These factors are calculated based on models recommended by the International Telecommunication Union for Radio-Wave Propagation (ITU-R P), available at https://www.itu.int/pub/R-REC.

Signal attenuation due to atmospheric layers includes tropospheric scintillation-induced signal fluctuations, modeled according to the ITU-R P618-13 guidelines. Additionally, attenuation through the ionospheric layer is accounted for by using the ITU-R P531-13 model, considering propagation loss influenced by the total electron content. These calculations are supplemented by data from the space weather access point (APF107).

Signal attenuation due to atmospheric absorption is based on the ITU-R P676-9 model, which employs ray tracing along the propagation path. The atmosphere is segmented into concentric shells, with attenuation computed empirically along each segment and aggregated to determine the total attenuation.

Signal attenuation due to rain is evaluated using the ITU-R P618-13 rain model, which assesses signal degradation as it travels through precipitation. This degradation primarily stems from water molecule absorption and varies with both frequency and elevation angle. Rain-induced signal loss escalates at higher frequencies and lower ground elevation angles because of the increased atmospheric path distances through rain-affected layers. Furthermore, rain contributes to higher antenna noise temperatures P_n. The annual rainfall rates are determined using the ITU-R P618-13 rain model, with an annual rain probability distribution of 21%–23% considered for Perth, Western Australia.

Signal attenuations due to clouds and fog are modeled by employing ITU-R P840-7. This model covers a comprehensive data set spanning 0° to 360° longitude and -90° to +90° latitude, with a spatial resolution of 1.125° in both latitude and longitude. The model provides interpolated total columnar content of reduced cloud liquid water applicable to any location on Earth’s surface. Rain loss modeling is included by utilizing annual average rainfall data adjusted for the mean surface temperature using insights from the ITU-R P1510 model. Additionally, the model calculates the depolarization effect on communication links due to rain.

Figure 3 (top) illustrates the attenuation values affecting the transmitted signal from a Starlink satellite with an assigned 11.7-GHz frequency (Starlink_1006_44712) received at CUT0 during the test day. Free space loss exhibits the highest values, reaching approximately 200 dB, primarily influenced by the distance between the satellite and receiver (Ippolito Jr, 2017). Conversely, losses due to clouds and fog are minimal, typically below 1 dB. Similar patterns are observed for other satellites and in the kinematic scenario, albeit with varying degrees of fluctuation depending on the elevation angle, as depicted in the bottom panel of Figure 3. The impact of free space loss is relatively stable, primarily reflecting the satellite–receiver distance.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Signal loss values for Starlink_1006_44712 (top) (UTC: coordinated universal time) and loss values with respect to the signal elevation angle (bottom)

To account for the impact of signal attenuations on the observations, we selected an EIRP value of 30 dBW for the transmitting antenna. This value represents the power radiated by the satellite in the direction of the receiver, taking into account the antenna gain. Although Starlink and OneWeb satellites typically have higher EIRP values and Iridium has a lower EIRP, this value is an average selection for the simulation. To define the quality of the receiving system, we set G_r / T to 20 dB/K for the receiving antenna, which is an average value for a high-quality receiving station. For the purpose of our simulations, we approximate the Doppler-shift variance using the closed-loop tracking model (Won & Pany, 2017):

σDr2=14π2T2BfC/N0(1+1TC/N0) 28

where B_f = 25 Hz is the tracking loop bandwidth, T = 5 ms is the integration time, and C / N₀ is converted from dB-Hz to a linear scale. Although this variance model with the chosen values is specifically applicable to closed-loop tracking scenarios, which may not be feasible in short-burst conditions, it serves as a reasonable approximation for our simulation purposes. Alternative tracking strategies, such as open-loop or batch-processing techniques, might be more suitable for actual implementations with intermittent signals.

Figure 4 depicts the number of available satellites for both static and kinematic receivers. Although the entire testing period for CUT0 spans 2 h to encompass one revolution of the tracked satellites, the time displayed in Figure 4 is aligned with the kinematic test duration (~10 min) for consistency. Three scenarios are evaluated: no limitations (i.e., no mask angle or short-burst limitation), inclusion of the elevation mask angle restriction, and inclusion of both the elevation angle and short-burst limitations. In the hypothetical no-limitation scenario, the receiver can access signals from all tested constellations. On average, CUT0 has 99 satellites available under these conditions, whereas the kinematic receiver has access to an average of 89 satellites. Implementing a 30° elevation mask angle limitation reduces the average number of satellites to 81 for CUT0 and 69 for the kinematic receiver. Further applying short-burst limitations decreases these numbers to 59 and 47, respectively. Despite the reduced number of satellites, these quantities still enable effective SoP positioning and facilitate the application of advanced techniques to enhance positioning accuracy.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Number of available Starlink satellites for CUT0 (top) and for the kinematic receiver (bottom)

To illustrate the impact of the applied limitations, Figure 5 shows observations from 10 satellites of the Starlink constellation. Applying the elevation mask angle limitation results in the removal of approximately 24% of observations across all satellites. This percentage increases to approximately 50% when the short-burst limitation is also applied.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

Doppler-shift frequency of 10 selected Starlink satellites when limitations are applied

Positioning using Doppler-shift observations based on Equations (11) and (21) was performed on the simulated observations to validate the developed models in both absolute and differential modes. The results are discussed in the following section.

4 POSITIONING RESULTS

This section presents the results of positioning using the developed Doppler-shift observation models. Results for the absolute positioning mode are presented first. Figure 6 compares the estimated coordinates derived from the model against the ground truth for each receiver. The positioning results shown in this figure are from the simulation described in Section 3. The testing duration for CUT0 spans 2 h to observe the impact of satellite orbital dynamics (i.e., for orbital periods of ~1.5 h) on positioning accuracy. The initial values for the state position vector are set more than 4,000 km away from the true values, with the same values considered for all simulations. The initial velocity vector and the receiver clock offset are also set to zero. Among the applied limitations, the elevation mask angle limitation has a relatively minor effect on accuracy changes, primarily because of the lesser reduction of observations (~24%) compared with the addition of the short-burst limitation, as explained in the previous section. However, the elevation mask angle limitation is expected to have a more significant impact in real-world scenarios because of the presence of actual tropospheric delays. Despite the small sample number of observations, which is sufficient for solving Equation (11), discontinuities in observations can still affect the accuracy of SoP-based positioning. This degradation is more pronounced in the kinematic case.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Comparing SoP-based positioning with the ground truth for CUT0 (top) and the kinematic receiver (bottom)

The RMSE values from this comparison are given in Table 2. For the static receiver (CUT0) under no limitations, the 3D RMSE is below 3.14 m. Introducing Doppler-shift observation restrictions based on elevation angle increases the 3D RMSE to 3.90 m, with increases observed in all dimensions. When both elevation angle and signal burst limitations are applied, the 3D RMSE increases to 4.32 m, indicating reduced positioning accuracy, primarily due to discontinuities in the observations.

For the kinematic receiver in absolute positioning mode, the no-limitation scenario results in a 3D RMSE of 3.81 m, slightly worse than the static case and with slightly increased errors in all directions. Applying the elevation angle limitation deteriorates the 3D RMSE to 3.98 m. Trends similar to those of the static case are observed for the kinematic mode, where positioning errors increase with elevation angle restrictions and discontinuous observations. This effect could be more severe in high-density urban areas, where sky visibility is limited. For example, limiting the number of satellites to less than 10 increased positioning errors to approximately 10 m in each direction.

To fully leverage broadband constellations, a higher number of available satellites is required. This can be achieved by designing receivers and antennas capable of collecting multi-constellation observations. Additionally, constellation operators might consider deploying omnidirectional beacons on satellites, which would eliminate the restrictions posed by the availability of a limited number of satellites and allow multi-satellite observation by a single user.

The accuracy levels provided above might be sufficient for coarse navigation applications. However, considering the actual orbital errors and clock instabilities of LEO satellites, similar to those discussed by Allahvirdi-Zadeh, Awange et al. (2022), would significantly impact real-time positioning based on Doppler-shift observations. In such cases, applying network-based POD along with simultaneous tracking and navigation frameworks or using differential positioning scenarios should be considered to reduce the impact of orbital and atmospheric errors and to increase positioning accuracy.

Table 3 provides the accuracy of the kinematic receiver positions calculated in the differential mode using the developed model described in the previous section. No matter which limitations are applied, the accuracy of all position components is better than 1 m. The no-limitation case outperformed the other cases, with the accuracy of all directions remaining below 30 cm. The 3D RMSE increases slightly to 53 cm when the number of satellites is limited based on their elevation angles. The same trend is observed for intermittent signal bursts, degrading the 3D accuracy to 62 cm. This level of accuracy is beneficial for various precise positioning applications. However, in real-world scenarios, the number of satellites might be reduced by natural environmental barriers in high-density urban areas. To simulate such conditions, the number of satellites was restricted to less than 10 in a testing scenario. The positioning accuracy in three directions was 0.7, 1.08, and 1.10 m, demonstrating possible positioning accuracy when the differential model in Equation (21) is applied.