State Estimation of a Moving Frequency Source From Observations at Multiple Receivers

  • NAVIGATION: Journal of the Institute of Navigation
  • March 2026,
  • 73
  • navi.745;
  • DOI: https://doi.org/10.33012/navi.745

Abstract

The task of position and velocity estimation of a body based on the Doppler effect, whether the body is transmitting or receiving a signal of known or unknown frequency, is a common problem arising in many different application domains. While prior works use nonlinear least-squares methods to solve the Doppler equations for positioning, such an approach requires an initial estimate to perform the initial linearization step. However, in certain scenarios (e.g., initial orbit determination), an initial estimate may not be readily available, and without a good initial estimate, the nonlinear least-squares solution may not converge. This work reformulates the Doppler-based positioning problem as a system of polynomial equations, allowing for a direct solution without any a priori state information. We then leverage homotopy continuation to obtain the global solution to the polynomial system. For an unknown transmitter (or receiver), we show that the data from six or seven receivers (or transmitters) is sufficient in the case of known or unknown frequency, respectively, to recover the unknown state up to finitely many possibilities. This technique provides a structured method for initializing nonlinear least-squares or sequential filters for further refinement of the estimated state, improving their convergence. After a brief development of the mathematics, two simple examples are provided: (1) initial orbit determination of a satellite emitting an electromagnetic signal and (2) position and velocity estimation of a vocalizing dolphin emitting an acoustic signal.

Keywords

1 INTRODUCTION

First described by C. A. Doppler in 1842 (Doppler & Studnica, 1842; Nolte, 2020), the Doppler effect has become a foundational principle in many scientific and engineering applications, including medicine (Kaunitz, 2016), astronomy (Lambourne, 1997), satellite orbit determination (Guier & Weiffenbach, 1958), and submarine localization (Kershner & Newton, 1962). In navigation, Doppler shift measurements enable estimation of a vehicle’s state by observing the frequency change in a signal exchanged between two bodies in relative motion.

One of the earliest examples of Doppler-based navigation was the use of Doppler radar, which was developed in the mid-1940s to provide accurate ground velocity measurements for commercial aircraft, enabling navigation via dead reckoning (Fried, 1993). By the 1960s, the Transit satellite system (Kershner & Newton, 1962; Stansell, 1978) became the first operational satellite-based navigation system using Doppler shift observables from radio frequency signals to enable global positioning. Transit operated until the mid-1990s, broadcasting at frequencies of 150 and 400 MHz, with the position fix computed iteratively using nonlinear least-squares methods (Stansell, 1978). Reliable convergence of the iterative position solution required the initial position estimate to be within a few hundred kilometers (Stansell, 1978). With at least three Doppler measurements obtained at different times, a lost submarine or terrestrial receiver could localize itself on the Earth’s surface with truly worldwide coverage for the first time.

More recently, Doppler-only positioning has been explored using signals of opportunity (SOPs) from large communication satellite constellations in low Earth orbit (LEO). Prior work has investigated batch point solutions for the user position and velocity, along with additional terms for the receiver clock bias (Psiaki, 2021), using SOPs from one or multiple LEO constellations (Kassas et al., 2023; Kozhaya et al., 2024). Given the uncertainty in the position and velocity of uncooperative LEO satellites, prior work has also examined performing simultaneous navigation and tracking using these Doppler observables (Kassas et al., 2024). Additionally, integration of LEO satellite SOP-based observables with inertial units in a tightly coupled filter has been explored for moving receivers (Kassas et al., 2023; Kassas et al., 2024; McLemore & Psiaki, 2022).

For LEO navigation, while state estimation convergence thresholds vary depending on factors such as satellite–receiver geometry, implementations in the literature suggest that initial position errors must generally be within a few hundred kilometers for successful state estimation (Stock et al., 2024). Some studies have demonstrated convergence with errors as large as 3600 km (Kassas et al., 2023), but common initialization strategies, such as those that assume an initial position at Earth’s center (as done for medium Earth orbit-based global navigation satellite systems [GNSSs] (Morton et al., 2021)), have proven insufficient for LEO Doppler positioning (Stock et al., 2024; Van Uytsel et al., 2022).

Beyond satellite-based navigation, Doppler-based localization has been applied in a variety of other fields as well. In marine biology, for example, Doppler-based localization can be used to determine the position and velocity of a vocalizing underwater animal (Gillespie et al., 2020). This method also enables localization of nodes within an underwater sensor network, which is widely used to collect data for underwater exploration, environmental monitoring, and natural disaster prevention (Felemban et al., 2015). More broadly, Doppler-based state estimation has been applied to tracking and navigation for moving vehicles on land, in the air, and in space, with numerous prior studies addressing these problems across various domains. These applications include underwater device detection (Amrita & Mou, 2021; Gong et al., 2020), drone localization and tracking (Chan & Jardine, 1990; Famili & Park, 2020), and initial orbit determination (IOD) (Christian et al., 2022; Ertl, 2023; Patton, 1960).

However, to the best of the authors’ knowledge, most of the existing prior works process Doppler measurements using nonlinear least-squares or sequential filtering (e.g., a Kalman filter or one of its variants), which inherently require a sufficiently accurate initial guess of the unknown state—yet, it is often unclear how such an a priori estimate might be obtained. Indeed, in certain scenarios, such as IOD, a poor initial guess can lead to divergence of the Doppler-based solution when using a traditional, nonlinear least-squares approach (as also demonstrated in Section 4). While a recent study attempted to address this initialization problem for spacecraft IOD (Ertl, 2023), the proposed approach suffered from prohibitive runtimes for the general case.

In this work, we explore how simultaneous Doppler shift measurements may be used to solve the state initialization problem via homotopy continuation, in the case of both known and unknown frequency at the source. A polynomial formulation of the problem allows the analyst to find the global solution to the minimal problem in the absence of a priori state information, thus providing a means of obtaining the initial guess required by earlier works. If the frequency at the source is known, six Doppler measurements suffice to determine the unknown three-dimensional position and velocity up to finitely many possibilities. If the transmitted frequency is unknown, seven measurements are needed, and the unknown frequency can be estimated as well. This method is especially useful in the case of over-determined systems, enabling existing filter-based solutions to solve these systems utilizing the initial guess provided by the approach outlined in this work.

The use of simultaneous observations relies on the assumption of synchronized observers. Loose synchronization may lead to biases in the Doppler frequency measurement owing to unmodeled clock drift, as well as incorrect temporal matching of the measurements collected by different receivers due to unmodeled clock biases. As a result, the applicability of the proposed technique is generally limited to scenarios for which zero or negligible clock drift is a reasonable assumption for the observers (e.g., ground stations measuring the signal emitted by an orbiting satellite). The experimental results in Section 4 provide further details on the tested Doppler measurement errors for which algorithm convergence was observed. In the case that the clock drift of a transmitter being localized is significant or unknown, we introduce a problem formulation that simultaneously estimates the transmitted frequency, along with the transmitter’s position and velocity, to account for this scenario.

The remainder of this paper is organized as follows. In Section 2, we provide the mathematical problem formulation as a polynomial expression, while in Section 3, we describe our solution method using homotopy continuation. In Section 4, we illustrate this method using several examples with applications ranging from underwater localization to Earth-based IOD. Finally, after concluding this work in Section 5, we additionally provide Macaulay2 source code for the proposed approach with homotopy continuation in the Appendix.

2 POLYNOMIAL FORMULATION

The approach presented in this work applies to both transmitter state estimation (tracking) and receiver state estimation (navigation). Because the Doppler equations remain the same in both cases — differing only in the interpretation of variables — we present our methodology as a tracking problem for clarity. However, the same method can be directly applied to navigation by reinterpreting the variables, without any modifications to the approach.

In this scenario, let us consider a transmitter at an unknown location r3 moving at an unknown velocity v3 and emitting a signal with a frequency f>0. We assume that the signal is detected by a receiver, placed at position ri moving at velocity vi. We assume that all receivers are synchronized with each other (e.g., a network of ground monitors that use the Global Positioning System to maintain time synchronization). Owing to the Doppler effect, the frequency observed by the receiver fi will be related to the frequency at the source through the following relation:

fi=(1ρ˙ic)f1

Here, c is the signal propagation speed in the given medium, and ρ˙i is the range rate, or the derivative of the range ρi between the transmitter and the receiver, represented as follows:

ρi=(rir)T(rir)2

By taking the derivative of Equation (2), we derive the following expression for the range rate ρ˙i in terms of the receiver and transmitter position and velocity:

ρ˙i=1ρi(rir)T(viv)3

A few simple manipulations lead to the polynomial expression we seek. Rearrangement of Equations (1)(3) allows us to write the product ρiρ˙i as follows:

ρiρ˙i=(rir)T(rir)(c(ffi)f)=(rir)T(viv)4

By squaring the middle and right-hand term of Equation (4), we obtain a polynomial expression relating the known parameters (ri,vi,fi) to the unknown parameters (r,v,f):

c2(ffi)2(rir)T(rir)f2[(rir)T(viv)]2=05

For numerical reasons, it is sometimes convenient to scale the frequencies by some common constant factor. To facilitate this and improve numerical performance, we can divide the previous expression by the frequency f2 and introduce the parameter k=1/f:

c2(1kfi)2(rir)T(rir)[(rir)T(viv)]2=06

The repetition of Equation (6) (or, equivalently, of Equation (5)) for seven receivers produces a polynomial system of seven equations with seven unknowns. As previously discussed in Section 1, this formulation holds under the assumption of zero clock drift for the receivers. Clock drift would introduce a bias term for each receiver, producing an under-determined system of equations. We note that if the transmit frequency is known, but the clock drift at the source is not, Equation (6) allows us to easily incorporate this indeterminate term in the problem by treating f as an unknown.

Note that the c2 coefficient in the left-most term of Equation (6) will likely cause performance degradation for applications with high propagation speeds (e.g., speed of light) on a finite-precision computer. In such cases, scaling of the parameters may provide considerable performance improvements (Higham, 2002). In particular, note that any scaling applied to the position vectors and to the frequency f does not change the number or location of the solutions.

If the transmit frequency f is known, the signal from six transmitters may be used to solve for the remaining six unknowns. In this case, the issues arising from c2 may be avoided by recognizing that the polynomial equation may be directly written in terms of the measured range rate ρ˙i instead of fi, yielding the following:

(rir)T(rir)ρ˙i2[(rir)T(viv)]2=07

Note that in the form of a navigation problem, Equation (6) and consequently Equation (7) remain unchanged. In this case, the symbols ri,vi, would be used to represent the known position and velocity of the transmitters, while r,v would represent the unknown position and velocity of the receiver.

3 NUMERICAL SOLUTION WITH HOMOTOPY CONTINUATION

The polynomial system formed by a set of frequency observations—each represented by its own instantiation of Equation (6)—can be solved in many different ways. Techniques such as Newton’s method may only find a single solution to the system or may fail to converge altogether, especially if the initial guess is poor or the system is highly nonlinear. Another family of techniques is homotopy continuation (Morgan, 1987; Sommese & Wampler, 2005), which tracks all solutions to the system during its solution process. This approach has much better convergence performance with a thoughtfully chosen initial setup, as we will discuss in more depth in this section. In this work, we use an implementation of parameter homotopy provided by the software package NAG4M2 (Leykin, 2011) in the computer algebra system Macaulay2 (Grayson & Stillman, 2025).

Homotopy continuation is a numerical method used to solve a square system of N equations with N unknowns. This method starts with a system of equations, known as the start system, for which the solutions are already know. The method then gradually deforms this system, tracking how the solutions change, until it transforms into the target system – the system we actually want to solve. Once this system deformation and tracking procedure is complete, because the final system coincides with the target system, the final system solutions also correspond to that of the desired problem. When the target consists of polynomial equations, a commonly used homotopy approach is parameter homotopy. In this type of homotopy, the deformation process corresponds to changing the parameters or coefficients of the (polynomial) system of equations, and a careful choice of the start system improves the numerical stability and convergence of this technique. Limitations in this approach arise in the presence of limited observability of the unknown state, e.g., a receiver and transmitter with negligible relative motions.

3.1 Solving a Polynomial System of Doppler Equations With Parameter Homotopy

Let x be the vector of unknowns, which, for our problem context, is a six- or seven-dimensional vector with an unknown position r and velocity v of the transmitter (or receiver) and a potentially unknown transmitted frequency f (or k). The system of N = 6 or 7 equations for our problem context corresponds to a set of N copies of Equation (6) for N receivers (or transmitters). For example, for the case of unknown frequency f, we can more explicitly define a family of systems T as follows:

T(r1:N,v1:N,f1:N):={c2(1kf1)2(r1r)T(r1r)[(r1r)T(v1v)]2=0c2(1kfi)2(rir)T(rir)[(rir)T(viv)]2=0c2(1kfN)2(rNr)T(rNr)[(rNr)T(vNv)]2=08

We notice that this system depends on some parameters, in particular, the locations and velocities of the receivers (or transmitters) r1:N,v1:N, and the Doppler frequency measurements f1:N. Once the parameters are set, they determine the coefficients of the monomials. When the values of these parameters match the actual known values from the observations (i.e., r1:N,v1:N, and f1:N correspond to the actual observers’ state and frequency measurements), the resulting system from the family of equations becomes the target system F that we aim to solve.

For parameter homotopy continuation, we select the start system S to belong to the same family as the target system S that we desire to solve, i.e., selecting ST. Whereas the target system F is characterized by the parameters of the problem we are interested in solving and therefore changes every time we wish to solve a different state estimation problem, the start system S is chosen and solved once and for all. The parameters and solutions of S can be stored and used in every other Doppler-based state estimation problem the user wishes to solve. The parameters of the start system S need not be physically meaningful; thus, they are not required to represent a realistic geometry of observation. Once the parameters have been selected, we can compute the solutions using a technique based on monodromy (see the work by Duff et al. (2019)). Section 3.3 provides further details on how we obtain a generalizable start system.

Next, beginning from the known solutions to the start system S, homotopy continuation follows a path in the space of the family of systems, transitioning from the start system to the target system while deforming the corresponding solutions. This transition is typically accomplished by using an adaptive predictor–corrector method, as depicted in Figure 1. Given that there are multiple solutions to any system in Equation (8), each solution defines a continuation path that must be tracked. Specifically, we track the solutions of the following system:

H(t)=(1t)S+tF9

as the continuation parameter t is varied from 0 to 1 using a shooting method. Note that t = 0 corresponds to the start system, whereas t = 1 corresponds to the target system we are interested in solving. To estimate the solutions for the current intermediate system H(t), we differentiate Equation (9) with respect to t:

H(t)xx(t)+H(t)(t)=010

Given the roots of H at a certain t0, we estimate the solutions at t0+δt by numerically integrating x(t) from Equation (10), thereby obtaining the guess x˜t0+δt. This guess is then refined using a corrector method (such as Newton’s method) to solve the intermediate system represented by Equation (9) at t0+δt, with the solution of the numerical integration serving as the initial guess. This procedure is schematically illustrated in Figure 1 for the case of a one-dimensional system. A short example code that shows the setup of the polynomial system and the necessary calls to perform the tracking in Macaulay2 is provided in the Appendix.

FIGURE 1

Depiction of the adaptive predictor–corrector method along the continuation path from the start system S to the target system F

Given the solution x0 of the system H at t = 0, we can obtain an initial guess of the solution x˜Δt at the following step t = Δt by numerical integration of Equation (10) (red arrow). This solution is then refined with a corrector method (blue arrow) to obtain the new estimate xΔt. Repetition of this procedure will ultimately lead to the solution of the system H at t = 1, which coincides with the target system F.

3.2 Benefits of Parameter Homotopy Over Total-Degree Homotopy

The construction of the family of polynomial systems is a key factor in this approach. This aspect is critical because the number of solutions of a generic start system corresponds to the number of continuation paths that must be tracked. For instance, many solvers implement the idea of total-degree homotopy. Total-degree homotopy, an approach that works for any square system, is attractive because it results in a simple start system whose solutions have a simple expression. However, the number of continuation paths equals the total degree in Bézout’s theorem (see chapter 18 of the work by Harris (1992)), which is often much larger than the number of target solutions we seek. The consequence is that total-degree homotopy often results in more paths than are necessary. Tracking these excess paths (which are divergent) significantly increases the computational cost. Reliance on total-degree homotopy likely led to exceptionally long runtimes in a similar approach attempted recently (Ertl, 2023). Depending on the application, using parameter homotopy (i.e., a start system within the same family as the target system) can provide a considerable speedup in comparison to total-degree homotopy.

3.3 Selecting a Start System and Exploiting Scenario-Specific Symmetries

A given system in the family T has multiple solutions, but the number of solutions remains constant for different choices of sufficiently generic parameters. The number of solutions is explored further in Sections 3.3.1 and 3.3.2. Note that the path associated with a single solution of the start system can lead to just one solution of the target system. As a result, our objective is to choose the parameters of the start system S so as to avoid special configurations (e.g., reflecting some special geometry of observation). Using parameters that reflect such configurations would result in fewer solutions and produce a start system that would be usable only in part when that special configuration is not satisfied by the observations. Because our solution space is complex, we selected the parameters via random sampling from complex numbers. We solved the start system using an approach based on monodromy (Duff et al., 2019).

We now apply these ideas to different scenarios with the objective of showing how application-specific symmetries may be exploited to reduce the computational cost required to track the number of continuation paths. Here, we consider two types of symmetry. Symmetries in the geometry of the observation affect the structure of the polynomial system T itself and can reduce the dimensionality of the solution set in the complex space. The second type of symmetries arises in the solution set itself. By establishing relationships between groups of solutions, we can identify a smaller subset of independent solutions. Only these independent solutions need to be determined by path tracking. The dependent solutions can be evaluated by leveraging the determined relationships. Both of these types of symmetries can significantly reduce the computational cost of solving the state estimation problem.

After introducing the setup for the general case of receivers and transmitters in generic relative motion, we analyze the simplified scenario in which an unknown transmitter is moving, but the known receivers are stationary in some common reference frame. This is the situation encountered, for example, in the localization of an underwater device emitting a signal received by multiple anchored sensors, but also in the IOD framework, when the receivers are bounded to move with the Earth’s surface. The analysis will also consider the dimension of the solution set for each problem formulation.

3.3.1 General Case: Navigation of a Receiver With Multiple Moving Transmitters

Let us consider the polynomial system given by Equation (6) specialized for seven different transmitters with an unknown transmit frequency f. According to Bézout’s theorem, we can determine an upper limit of the number of solutions for a polynomial system by multiplying the total degree of each polynomial in the system. In this case, we can expect up to 47 = 16, 384 solutions in the form of complex-valued seven-tuples (r, v, f ). However, the actual number of solutions is often much smaller than that specified by Bézout’s theorem. As a result, in practice, researchers often use techniques based on the Gröbner basis to determine the actual number of solutions (Bose & Bose, 1995; Cox et al., 2015). With this approach, we found that the actual number of solutions is 672. Using total-degree homotopy, which requires tracking the number of solutions identified by Bézout’s theorem, almost all paths diverge, and many hours (on contemporary laptop computers) are required to converge to the 672 solutions. Instead, parameter homotopy allows us to track only 672 solutions and solve the desired polynomial system more quickly. We can repeat a similar reasoning for the case of known frequency f, where the total-degree homotopy requires the tracking of 46 = 4, 096 solutions as six-tuples (r, v), whereas the parameter homotopy allows us to track the actual number of 128 solutions.

We have experimentally demonstrated that most of the 672 (for unknown f ) or 128 (for known f ) complex solutions can be eliminated without using observations from any additional transmitters. In fact, even if all of the solutions satisfy the polynomial system obtained from Equation (6), not all satisfy the original (i.e., non-squared) relation of Equation (4). However, it is not guaranteed that this approach will eliminate all but one candidate. If an additional observation is available, we can discard the remaining “wrong” solutions by assessing the residuals of Equation (4) written for the additional transmitter. The screening procedure described here, combined with a simple consideration of the geometry of the observations, was sufficient to uniquely identify the desired solution in all of the experiments we performed.

3.3.2 Special Case: Localization of a Transmitter With Multiple Stationary Receivers

Although the general case encompasses all classes of observer motion (both moving and stationary), the special case of stationary receivers (vi = 0) simplifies the polynomial systems and permits a more efficient solution. In particular, setting vi = 0 causes some of the terms in the original polynomial to vanish. Furthermore, we observe that if (r, v, f ) is a solution to the polynomial system, then (r, −v, f ) must also be a solution. Counting the number of solutions for this special case, we found that the number of solutions of the polynomial system with stationary observers is 296 for an unknown transmit frequency f. However, owing to the symmetry of the problem, we can discard half of these solutions and track only 148 paths. Furthermore, if the frequency at the source is also known, then we only need to track 24 paths for the reduced system. This reduction in the number of paths increases the speed of the solver.

Localization of a Transmitter With Multiple Terrestrial Receivers

The problem of receivers located on the Earth’s surface and rotating with it produces an observing scenario analogous to that of stationary receivers. In fact, as also derived by Liu et al. (2019), the contribution of the Coriolis acceleration to the Doppler shift measurement is zero.

To understand this, consider the Earth-centered Earth-fixed (ECEF) frame. Let us assume that we have multiple observers located on the Earth’s surface (potentially at different altitudes). In addition, assume that we have an orbiting transmitter, whose state we wish to estimate. We will now show that the same form of the Doppler equation holds in the ECEF and in the Earth-centered inertial (ECI) frame.

If we let v and vF be the ECI and ECEF velocities of the transmitter, respectively, we have the following:

v=vF+ω×r11

where ω is the angular velocity of the Earth. Because the observer lies on the Earth’s surface, its velocity is as follows:

vi=ω×ri12

Substituting these expressions for v and vi in Equation (4), we obtain the following:

ρiρ˙i=(rir)T(ω×rivFω×r)=013

corresponding to the following:

ρiρ˙i=(rir)T[ω×(rir)vF]14

Because the following holds:

(rir)T[ω×(rir)]=015

then Equation (14) is as follows:

ρiρ˙i=(rir)T(vF)16

which is the same as Equation (4) when vi = 0.

4 EXPERIMENTS

The efficacy of the proposed solution is demonstrated through two illustrative examples. The first example examines IOD for an Earth-based satellite transmitter using Doppler frequency observations from terrestrial receivers. We consider both scenarios in which the satellite’s transmitted frequency is known or unknown. The second example focuses on underwater localization of a dolphin emitting a kilohertz-level sonar whistle underwater (Oswald et al., 2003).

4.1 Example 1: IOD With Terrestrial Receivers

This first experiment consists of the IOD of a satellite emitting an electromagnetic signal. The receivers are located at different altitudes above the Earth’s surface and move along circles of latitude. We assume that the terrestrial receivers are synchronized, which can be readily achieved through the GNSS to maintain alignment with Coordinated Universal Time. The orbital elements (Bate et al., 2020) for the transmitting spacecraft are assumed as follows:

a=12,000 kme=0.1i=20Ω=200ω=20v=0

where a is the semi-major axis, e is the eccentricity, Ω is the right ascension of the ascending node, i is the inclination, ω is the argument of the periapsis, and v is the true anomaly at epoch. The satellite is assumed to transmit a signal with frequency f = 9 GHz.

These orbital elements may be converted to the three-dimensional inertial position and velocity using standard equations found in most modern astrodynamics textbooks (Bate et al., 2020; Curtis, 2020; Vallado, 2007):

r=[8349.56732.81263.4]kmv=[3.97214.54172.0478]km/s17

The geometry of the observation is shown in Figure 2.

FIGURE 2

Observers (red) and the transmitting satellite (black) with the corresponding direction of motion

We now assume that the ground stations measure the frequencies fi using a frequency-locked loop (FLL) (Hegarty & Kaplan, 2005). The tracking loop measurement error, mainly driven by thermal noise σt,f and the dynamic stress error fe, may be modeled as follows (Hegarty & Kaplan, 2005):

σf,i=σt,f+13fe,i18

Letting T be the prediction integration time, C/N0 the received signal’s carrier-to-noise-density, and Bn the carrier loop noise bandwidth, the thermal noise in the FLL can be modeled as follows (Hegarty & Kaplan, 2005):

σt,f=12πT4FBnC/N0[1+1TC/N0]19

where F =1 for C/N0 values above the threshold associated with the 1σ tracking threshold of 1/(12T), with F =2 otherwise. Assuming a second-order FLL, the dynamic stress error may be calculated as follows (Hegarty & Kaplan, 2005):

fe,i=1360ω02d3ρidt320

where we remind the reader that ρi is the range between the satellite and the ground station. The parameter ω0 can be calculated using the relation Bn = 0.53ω0 (Hegarty & Kaplan, 2005). Assuming a carrier loop bandwidth of 2 Hz and an integration time of 5 ms, we calculated the dynamic stress error by numerically evaluating the third derivative of the range (i.e., the second derivative of the range rate) at the time of the observation for each ground station and used the resulting σf,i value. With this setup, we ran multiple Monte Carlo simulations of 1, 000 runs each for values of C/N0 ranging between 15 and 36 dB-Hz.

4.1.1 Baseline Method: Solving for the State With Standard Nonlinear Least Squares

The objective of the homotopy-based solution is to provide an initial guess that other filter-based approaches can then utilize when incorporating additional information (e.g., further Doppler and non-Doppler measurements, dynamics). Therefore, we start by examining the performance of a standard framework for nonlinear least squares.

In this section, we will show that with the conventional method (using a standard nonlinear least-squares approach), a reasonable initial guess on the unknown state is often not sufficient to guarantee convergence. This result highlights the need for a more structured initialization solution (e.g., the homotopy-based approach proposed in this work). While this example is considered for the case of a known frequency at the source, the procedure is analogous in the case of unknown f. Following standard approaches, we will work with the range rate measurement ρ˙i, instead of the frequency measurement fi.

We begin by considering the expression of ρ˙i in Equation (3), from which we may construct the following vector function:

F(r,v)=F(x)=[1ρ1(r1r)T(v1v)1ρN(rNr)T(vNv)]21

where N is again the number of observations. The Jacobian JF can be analytically determined by taking the partial derivatives of F with respect to the state. By letting x0 be an initial guess of the unknown state, the optimal solution is estimated by minimizing the difference between the calculated F vector and the measurements. With Newton’s method, at each step k, the estimate x0 is updated as follows:

x0,k=x0,k1JF(x0,k1)1(F(x0,k1)[ρ˙1ρ˙N])22

To solve the problem with the baseline approach using nonlinear least squares, an initial guess is needed. Without taking advantage of the initialization procedure shown in this paper, we attempt to make a reasonable initial guess. If multiple ground stations, located at different latitude and longitude values, are able to receive a signal from the satellite, an equatorial orbit represents a sensible starting point. However, two choices need to be made: the altitude of the orbit and the direction of the velocity. While some fortuitous choices of orbit radius (e.g., r = 10, 000 km) and velocity direction will enable convergence of the least-squares solution in all of the Monte Carlo simulations, choosing the wrong direction of motion or a less fortunate radius (e.g., r = 12, 000 km) is sufficient to prevent the solver from converging to the solution in most cases. In fact, if we assume the physically plausible initial guess of rT = [9192.6, 7713.5, 0] and vT = [−3.7046, 4.415, 0], only one of the 8,000 total runs of the Monte Carlo analysis converges to the true solution when solving the nonlinear least squares in the ECEF frame.

Thus, we observe that with the conventional nonlinear least-squares approach, reasonable initial guesses can lead to divergence. This finding highlights the importance of a good initial guess, such as the one provided by the homotopy continuation method described in this paper.

4.1.2 Solving for the State With Homotopy Continuation

Using the approach described in this paper, we can solve for the unknown state without any initial guesses, avoiding the disadvantages pointed out in Section 4.1.1.

In the ECEF frame, we use the setup of Section 3.3.2 to determine the 24 independent complex-valued candidate solutions. Homotopy tracking was performed on a machine with an Intel Core i9-13900HX (32 CPUs), 32 GB of RAM. The run-time values are reported in Table 1.

View this table:
TABLE 1 Mean and Standard Deviation of the Tracking Time Utilized to Recover the 24 Independent Solutions of the IOD Problem with Terrestrial Receivers for a Monte Carlo Simulation of 1,000 Runs

Because we are interested in real-valued solutions, we can begin filtering out undesired solutions by discarding the complex ones. Once this step is completed, recalling the sign symmetry described in Section 3.3.2, we double the solution set by recalling that, for each pair (r, v), the symmetric solution (r, −v) should also be considered.

As described in Section 3.3.1, not all of these remaining solutions satisfy Equation (7). Finally, the evaluation of Equation (7) for an additional receiver is usually sufficient to identify the unique desired solution.

It is important to note that this procedure might not be sufficient to confidently discard all of the undesired solutions. When the residuals of two or more solutions were similar, we utilized a visibility constraint to rule out the undesired solution. In particular, we only considered candidate position solutions for which the line-of-sight vector forms a positive cosine with the receiver antenna’s pointing direction—that is, the line-of-sight vector lies within the hemisphere defined by the antenna’s boresight.

The output of the homotopy-based approach is reported in Tables 2 and 3, as well as Figures 3 and 4. When the least-squares solution of Section 4.1.1 was initialized with the output of the homotopy-based approach, the percentage of runs converging to the true solution rose to 100%, potentially enabling, at the same time, an improvement in the solution estimate by simply including additional measurements. These results demonstrate the level of Doppler error that the homotopy-based approach can tolerate. In particular, even at the lowest tested C/N0 of 15 dB-Hz, which corresponds to a 1σ Doppler error of 61 Hz, the method consistently converged across all 1,000 Monte Carlo trials, suggesting robustness in the tested scenario to relative clock drift errors up to a comparable magnitude.

View this table:
TABLE 2 Error in the Position Estimate for Monte Carlo Simulations of 1,000 Runs at Various Levels of C/N0 Under the Assumption of Known Transmitted Frequency Obtained Through Homotopy Results are reported in the radial (R), along-track (T), and cross-track (C) directions. std: standard deviation
View this table:
TABLE 3 Error in the Velocity Estimate for Monte Carlo Simulations of 1,000 Runs at Various Levels of C/N0 Under the Assumption of Known Transmitted Frequency Obtained Through Homotopy Results are reported in the radial (R), along-track (T), and cross-track (C) directions.
FIGURE 3

Visualization of the position error for the Monte Carlo simulation represented in Table 2

FIGURE 4

Visualization of the velocity error for the Monte Carlo simulation represented in Table 3

Finally, we performed one last experiment assuming that the frequency transmitted by the orbiting satellite was unknown. Using a value of C/N0 = 30 dB-Hz, we ran a Monte Carlo simulation and obtained the results shown in Figure 5. For comparison with the case of known frequency, we also report the mean and standard deviation of the estimated state obtained through homotopy in the radial, along-track, and cross-track frame RTC:

FIGURE 5

Results of a Monte Carlo simulation of 1,000 runs with C/N0 = 30 dB-Hz for the IOD of a satellite emitting a signal with unknown frequency f

Over all Monte Carlo runs, we observe consistent convergence to the correct position solution, with 3σ errors of < 10 km in each dimension.

ΔrRTC=[1.82×102±1.681.10×102±2.895.69×102±1.02]kmΔvRTC=[4.06×104±6.42×1023.50×104±4.39×1021.71×105±1.97×103]km/sΔf=(4.04±1.73×102) Hz

We observe that the position and velocity errors are marginally higher compared with the known-frequency case. Nonetheless, in both scenarios, the positioning errors of the IOD solution remain within kilometer-level accuracy.

4.2 Example 2: Underwater Localization With Stationary Receivers

The second illustrative experiment is the simultaneous position and velocity estimation of a vocalizing Tursiops truncatus (the common bottlenose dolphin) swimming within a volume equipped with a network of acoustic sensors (e.g., hydrophones). We assume that the receivers are stationary, with well-known positions from prior surveys. Recognizing that these dolphins have typical swim speeds of 1–4 m/s (Wursig & Wursig, 1979) and their whistles have a frequency of approximately 7.4–17.2 kHz (Oswald et al., 2003), we simulated a dolphin with the following state that is unknown a priori to the sensor network:

r=[5.235.2815.00]mv=[1.381.530.22]m/sf=15 kHz

When provided perfect measurements (i.e., no noise), we recover the dynamical state (position and velocity) of the source to machine precision in less than 0.5 s of runtime for unknown f (N = 7 observations) and in less than 0.1 s for known f (N = 6 observations), in this case of stationary receivers. Owing to the simplicity of the problem’s mathematical structure, an approach similar to that of Hruby et al. (2022) may lead to an exceptionally fast solver with an execution time of less than a millisecond per problem instance. Implementation of such a speed-up was not pursued further here, as the current sub-second runtimes were sufficiently short for our proof-of-concept purposes.

We also conducted a Monte Carlo analysis to study the effect of imperfect receiver state knowledge and noisy frequency measurements. Figure 6 shows the results from 1,000 numerical experiments for both a known (N = 6) and unknown (N = 7) whistle frequency. In each experiment, we applied additive, zero-mean Gaussian noise to each of the observers’ positions (σr = 1.5 cm) in each dimension and to each of the observed frequencies (σf = 0.1 Hz).

FIGURE 6

State errors in position (left column) and velocity (right column) for the example of a vocalizing dolphin (underwater application)

Frequency estimate errors are shown at the bottom for the case for N = 7 receivers. Results were obtained from a 1,000-trial Monte Carlo analysis. Over all Monte Carlo runs, we observe consistent convergence to the correct position solution, with 3σ errors of < 2 km in each dimension.

5 CONCLUSIONS

The position and velocity of an object emitting an acoustic or electromagnetic signal may be inferred from only the apparent frequency as measured by a multitude of known receivers. Similarly, the position and velocity of a receiver may be inferred from the measurement of the frequency of a signal emitted by a multitude of transmitters. This state estimation problem may be written as a system of polynomials, which permits an efficient solution via homotopy. The frequency at the source may be either known or unknown. We highlight the importance of using parameter homotopy, instead of total-degree homotopy, and provide a strategy to filter out multiple solutions.

The approach solves the state estimation problem without using any initial guess on the unknown state. This feature makes this approach a good candidate for initialization of other state estimation techniques that are tailored to incorporate additional measurements and more complex dynamical modeling.

Two example applications show the feasibility of this procedure, both in perfect (noise-free) and imperfect (noisy) scenarios. The proposed method provided a noisy estimate of the transmitter state consistent with the amount of noise in the measurements, both in the case of known and unknown frequency.

HOW TO CITE THIS ARTICLE:

Mancini, M., Mina, T., Leykin, A., & Christian, J.A. (2026). State estimation of a moving frequency source from observations at multiple receivers. NAVIGATION, 73. https://doi.org/10.33012/navi.745

ACKNOWLEDGMENTS

This work was partially supported by the Air Force Office of Scientific Research award FA95502310512 and the National Science Foundation Division of Mathematical Sciences award 2001267. The authors thank Rickey Huang for valuable discussions related to the number of polynomial solutions and the implementation of Macaulay2 code. We also thank Dr. Timothy Duff for assistance with developing the monodromy code used to acquire the solutions to our start system.

APPENDIX

This appendix provides example code using Macaulay2, which is a popular software system for solving polynomial systems and other problems in algebraic geometry. For more information on Macaulay2, see the work by Grayson and Stillman (2025).

The following code example shows how to build a polynomial system based on multiple instantiations of Equation (6). Note that it is assumed that the velocity provided is already scaled by the signal’s speed c. The variables are named on the basis of the navigation framework.

Given the polynomial system G produced by the above code, we may write another script to track the solutions using parameter homotopy continuation. It is assumed that the parameters describing the start and target system have already been loaded in the 1 × (7 nObs + 1) vectors startParameters and targetParameters. The same holds for the solutions of the start system, contained in the matrix startSolutions, of dimension (number of unknowns × number of solutions). In practice, we obtain startSolutions using MonodromySolver, a package of Macaulay2 based on the work by Duff et al. (2019).

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.

REFERENCES

  1. Amrita, D., & Mou, D. (2021). On accurate localization of sensor nodes in underwater sensor networks: A Doppler shift and modified genetic algorithm based localization technique. Evolutionary Intelligence, 14(11), 119131. https://doi.org/10.1007/s12065-019-00343-1
  2. Bate, R., Mueller, D., White, J., & Saylor, W. (2020). Fundamentals of astrodynamics (2nd ed.). Dover.
  3. Bose, N., & Bose, N. (1995). Gröbner bases: An algorithmic method in polynomial ideal theory. Springer.
  4. Chan, Y.-T., & Jardine, F. (1990). Target localization and tracking from Doppler-shift measurements. IEEE Journal of Oceanic Engineering, 15(3), 251257. https://doi.org/10.1109/48.107154
  5. Christian, J., Ertl, C., Horneman, K., & Lovell, A. (2022, January). Doppler-only initial orbit determination for an orbiting transmitter. In AIAA SciTech Forum. https://doi.org/10.2514/6.2022-1000
  6. Cox, D. A., Little, J., & O’Shea, D. (2015). Ideals, varieties, and algorithms: An introduction to computational algebraic geometry and commutative algebra (4th ed.). Springer. https://doi.org/10.1007/978-3-319-16721-3
  7. Curtis, H. (2020). Orbital mechanics for engineering students (4th ed.). Elsevier.
  8. Doppler, C. A., & Studnica, F. J. (1842). Ueber das farbige licht der doppelsterne und einiger anderer gestirne des himmels. Prag K. Böhm.
  9. Duff, T., Hill, C., Jensen, A., Lee, K., Leykin, A., & Sommars, J. (2019). Solving polynomial systems via homotopy continuation and monodromy. IMA Journal of Numerical Analysis, 39, 14211446. https://doi.org/10.1093/imanum/dry017
  10. Ertl, C. A. (2023). Solutions and analysis of multivariate polynomial systems for geolocation and initial orbit determination [PhD dissertation, Rensselaer Polytechnic Institute].
  11. Famili, A., & Park, J.-M. J. (2020, May). ROLATIN: Robust localization and tracking for indoor navigation of drones. In 2020 IEEE Wireless Communications and Networking Conference (WCNC), 16. https://doi.org/10.1109/WCNC45663.2020.9120619
  12. Felemban, E., Shaikh, F. K., Qureshi, U. M., Sheikh, A. A., & Qaisar, S. B. (2015). Underwater sensor network applications: A comprehensive survey. International Journal of Distributed Sensor Networks, 11(11), 896832. https://doi.org/10.1155/2015/896832
  13. Fried, W. R. (1993). History of Doppler radar navigation. NAVIGATION, 40(2), 121136. https://doi.org/10.1002/j.2161-4296.1993.tb02299.x
  14. Gillespie, D., Palmer, L., Macaulay, J., Sparling, C., & Hastie, G. (2020). Passive acoustic methods for tracking the 3D movements of small cetaceans around marine structures. PLoS ONE, 15(5), e0229058. https://doi.org/10.1371/journal.pone.0229058
  15. Gong, Z., Li, C., Jiang, F., & Zheng, J. (2020). AUV-aided localization of underwater acoustic devices based on Doppler shift measurements. IEEE Transactions on Wireless Communications, 19(4), 22262239. https://doi.org/10.1109/TWC.2019.2963296
  16. Grayson, D. R., & Stillman, M. E. (2025). Macaulay2, a software system for research in algebraic geometry. http://www2.macaulay2.com
  17. Guier, W. H., & Weiffenbach, G. C. (1958). Theoretical analysis of Doppler radio signals from Earth satellites. Nature, 181, 15251526. https://doi.org/10.1038/1811525a0
  18. Harris, J. (1992). Algebraic geometry. A first course (Vol. 133). Springer. https://doi.org/10.1007/978-1-4757-2189-8
  19. Hegarty, C., & Kaplan, E. (2005). Understanding GPS principles and applications (2nd ed.).
  20. Higham, N. J. (2002). Accuracy and stability of numerical algorithms. SIAM.
  21. Hruby, P., Duff, T., Leykin, A., & Pajdla, T. (2022, June). Learning to solve hard minimal problems. In Proc. of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 55325542. https://doi.org/10.1109/CVPR52688.2022.00545
  22. Kassas, Z., Kozhaya, S., Saroufim, J., Kanj, H., & Hayek, S. (2023). A look at the stars: Navigation with multi-constellation LEO satellite signals of opportunity. Inside GNSS. https://insidegnss.com/a-look-at-the-stars-navigation-with-multi-constellation-leo-satellite-signals-of-opportunity/
  23. Kassas, Z. M., Khairallah, N., & Kozhaya, S. (2024). Ad astra: Simultaneous tracking and navigation with megaconstellation LEO satellites. IEEE Aerospace and Electronic Systems Magazine, 39(9), 4671. https://doi.org/10.1109/MAES.2023.3267440
  24. Kaunitz, J. D. (2016). The Doppler effect: A century from red shift to red spot. Digestive Diseases and Sciences, 61, 340341. https://doi.org/10.1007/s10620-015-3998-9
  25. Kershner, R. B., & Newton, R. R. (1962). The Transit system. Journal of Navigation, 15(2), 129144. https://doi.org/10.1017/S0373463300035943
  26. Kozhaya, S., Saroufim, J., & Kassas, Z. M. (2024, December). Opportunistic positioning with Starlink and OneWeb LEO Ku-band signals. In 2024 IEEE International Conference on Wireless for Space and Extreme Environments (WiSEE) (pp. 112117). https://doi.org/10.1109/WiSEE61249.2024.10850237
  27. Lambourne, R. (1997). The Doppler effect in astronomy. Physics Education, 32(1), 34. https://doi.org/10.1088/0031-9120/32/1/017
  28. Leykin, A. (2011). Numerical algebraic geometry. Journal of Software for Algebra and Geometry, 3, 510. https://doi.org/10.2140/jsag.2011.3.5
  29. Liu, Q., Drake, S. P., & Anderson, B. (2019). Mapping target location from Doppler data. ArXiv. https://doi.org/10.48550/arXiv.1903.04979
  30. McLemore, B., & Psiaki, M. L. (2022). Navigation using Doppler shift from LEO constellations and INS data. IEEE Transactions on Aerospace and Electronic Systems, 58(5), 42954314. https://ieeexplore.ieee.org/document/9743612
  31. Morgan, A. (1987). Solving polynomial systems using continuation for engineering and scientific problems. Prentice Hall Inc.
  32. Morton, Y. J., van Diggelen, F., Spilker Jr, J. J., Parkinson, B. W., Lo, S., & Gao, G. (2021). Position, navigation, and timing technologies in the 21st century: Integrated satellite navigation, sensor systems, and civil applications (Vol. 1). John Wiley & Sons.
  33. Nolte, D. D. (2020). The fall and rise of the Doppler effect. Physics Today, 73(3), 3035. https://doi.org/10.1063/PT.3.4429
  34. Oswald, J. N., Barlow, J., & Norris, T. F. (2003). Acoustic identification of nine delphinid species in the eastern tropical pacific ocean. Marine Mammal Science, 19(1), 20037. https://doi.org/10.1111/j.1748-7692.2003.tb01090.x
  35. Patton, R. B. (1960). Orbit determination from single pass Doppler observations. IRE Transactions on Military Electronics, MIL-4(2/3), 336344. https://doi.org/10.1109/IRET-MIL.1960.5008246
  36. Psiaki, M. L. (2021). Navigation using carrier Doppler shift from a LEO constellation: TRANSIT on steroids. NAVIGATION, 68(3), 621641. https://doi.org/10.1002/navi.438
  37. Sommese, A. J., & Wampler, C. W., II. (2005). The numerical solution of systems of polynomials. World Scientific Publishing Co.
  38. Stansell, T. A. (1978). The TRANSIT navigation satellite system. Magnavox.
  39. Stock, W., Schwarz, R. T., Hofmann, C. A., & Knopp, A. (2024). Survey on opportunistic PNT with signals from LEO communication satellites. IEEE Communications Surveys & Tutorials, 27(1), 77107. https://doi.org/10.1109/COMST.2024.3406990
  40. Vallado, D. (2007). Fundamentals of astrodynamics and applications (3rd ed.). Microcosm Press.
  41. Van Uytsel, W., Janssen, T., Halili, R., & Weyn, M. (2022). Exploring positioning through pseudoranges using low earth orbit satellites. International Conference on P2P, Parallel, Grid, Cloud and Internet Computing, 278287. https://doi.org/10.1007/978-3-031-19945-5_28
  42. Wursig, B., & Wursig, M. (1979). Behavior and ecology of the bottlenose dolphin (tursiops truncatus) in the south atlantic. Fishery Bulletin, 77(2), 399412.
Loading
Loading
Loading
Loading