Open AccessArticle

Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm

Zichen Xu

^†

Zongnan Li

^†,

Xiaohui Liu

^*,

Zhimin Ji

Qianqian Wu

Hao Liu

and

Chao Wen

Key Laboratory of Satellite Navigation Technology, College of Electronic Science and Technology, National University of Defense Technology, Changsha 410005, China

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Remote Sens. 2024, 16(16), 2958; https://doi.org/10.3390/rs16162958

Submission received: 28 June 2024 / Revised: 9 August 2024 / Accepted: 10 August 2024 / Published: 12 August 2024

(This article belongs to the Special Issue GNSS Positioning and Navigation in Remote Sensing Applications)

Download

Browse Figures

Figure 1
Schematic diagram of an ill−posed problem. "> Figure 2
The sub-figures (a–d) show the linearization errors for initial values set on a spherical surface with the true value as the center and a radius of 200 km at orbit altitudes of 200 km, 400 km, 600 km, and 800 km. The horizontal and vertical coordinates represent the elevation and azimuth angles of the initial values relative to the true position in the ECEF coordinate system. "> Figure 3
The schematic diagram of the Doppler equation in two−dimensional space, where the x and y axes represent the positions of the receiver in meters, and the z axis represents the Doppler frequency shift at different positions. Sub-figure (a) displays a three-dimensional view of the objective function set and the affine set, while sub-figure (b) provides a top-down view. "> Figure 4
The schematic diagram of the objective function for Doppler positioning in a two−dimensional space. Sub-figure (a) displays a three-dimensional view of the objective function set and the affine set, while sub-figure (b) provides a top-down view. "> Figure 5
Relationship between C value and the number of visible satellites. "> Figure 6
Algorithm flowchart. "> Figure 7
The skyplot of visible satellites for the non-convergence case study. "> Figure 8
The convergence status of the algorithm using the LS estimation iteration is shown in the figure above. The left subplot displays the condition number of the iteration matrix, while the right subplot illustrates the loss function during the iteration process. "> Figure 9
The convergence status of the algorithm using the algorithm proposed in this paper is shown in the figure above. The left subplot displays the condition number of the iteration matrix, while the right subplot illustrates the loss function during the iteration process. "> Figure 10
The skyplot of a case study involving local non-uniqueness of solutions. "> Figure 11
The relationship between the logarithmic loss and the number of iterations is depicted in the figure above. The left subplot corresponds to the algorithm proposed in this paper, while the right subplot corresponds to the iteration LS estimation. "> Figure 12
Simulation trajectory. "> Figure 13
The graph of the variation of the number of visible satellites over time. "> Figure 14
The positioning error calculated using the algorithm presented in this paper. "> Figure 15
The velocity error computed using the algorithm proposed in this paper in m/s. "> Figure 16
The clock drift error computed using the algorithm proposed in this paper in s/s. "> Figure 17
The positioning error computed using the iterative LS estimation in s/s. "> Figure 18
The velocity error computed using the iterative LS estimation in s/s. ">

Versions Notes

Abstract

Doppler positioning, as an early form of positioning, has regained significant research interest in the context of Low Earth Orbit (LEO) satellites.Given the LEO mega-constellation scenario, the objective function of Doppler positioning manifests significant nonlinearity, leading to ill-conditioning challenges for prevalent algorithms like iterative least squares (LS) estimation, especially in cases where inappropriate initial values are selected. In this study, we investigate the causes of ill-posed problems from two perspectives. Firstly, we analyze the linearization errors of the Doppler observation equations in relation to satellite orbital altitude and initial value errors, revealing instances where traditional algorithms may fail to converge. Secondly, from an optimization theory perspective, we demonstrate the occurrence of convergence to locally non-unique solutions for Doppler positioning. Subsequently, to address these ill-conditioning issues, we introduce Tikhonov regularization terms in the objective function to constrain algorithm divergence, with a fitted model for the regularization coefficient. Finally, we conduct comprehensive simulation experiments in both dynamic and static scenarios to validate the performance of the proposed algorithm. On the one hand, when the initial values are set to 0, our algorithm achieves high-precision positioning, whereas the iterative LS fails to converge. On the other hand, in certain simulation scenarios, the iterative LS converges to locally non-unique solutions, resulting in positioning errors exceeding 50 km in the north and east directions, several hundred kilometers in the vertical direction, and velocity errors surpassing 120 m/s. In contrast, our algorithm demonstrates typical errors of a position error of 6.8462 m, velocity error of 0.0137 m/s, and clock drift error of 8.3746 ×

10^{- 6}

s/s. This work provides an effective solution to the sensitivity issue of initial points in Doppler positioning and can serve as a reference for the algorithm design of Doppler positioning receivers with LEO mega-constellations.

Keywords:

low earth orbit; Doppler position; ill-posed problems; Tikhonov regularization

1. Introduction

In recent years, the construction of Low Earth Orbit (LEO) satellites has been in a phase of rapid deployment. Currently, more than 6000 LEO satellites provide services, with a planned deployment total exceeding 40,000 satellites [1]. In response to the positioning and navigation needs under Global Navigation Satellite System (GNSS) denial situations caused by complex electromagnetic environments and weak signal power, countries led by the United States have begun researching backup navigation technologies independent of the GNSS [2]. On the one hand, the terrestrial power level of LEO satellite signals surpasses that of GNSS signals by 30 dB [3]. On the other hand, LEO satellite signals exhibit a pronounced Doppler frequency shift, with the associated information being readily accessible [4]. Consequently, positioning technologies based on the Doppler frequency shift of LEO satellites have emerged as a focal point in backup navigation system research, particularly in the context of extensive satellite-based communication services [5].

Research related to satellite Doppler positioning technology can be traced back to 1958. At that time, the Johns Hopkins University Applied Physics Laboratory, commissioned by the U.S. Navy, was developing the Transit satellite navigation system based on Doppler positioning technology [6,7]. From the 1970s to the 1990s, satellite Doppler positioning technology was widely applied in the field of geodesy [8]. During this period, satellite Doppler positioning technology became an important means of long-distance positioning, especially in areas such as oceans and polar regions where other positioning methods were difficult to implement. After 1995, with the widespread application of the Global Positioning System (GPS), Doppler positioning technology was replaced in some fields, but it still maintains unique advantages in certain specific applications. For example, the ARGOS satellite for animal tracking [9] and the Doris system for orbit determination of space probes [10] continue to develop. In recent years, Doppler positioning with LEO mega-constellation has been proposed [11], attracting increasing attention from scholars.

The Doppler positioning of LEO satellites can be divided into three categories. The first category is single-satellite integrated Doppler positioning, exemplified by the TRANSIT navigation system that provided positioning and navigation services to the public as early as 1968 [12]. This system achieved positioning using approximately 2 min of Doppler frequency shift observations, with a single-point positioning accuracy of about 100–200 m [13]. For stationary receivers, sub-meter-level positioning accuracy can be achieved using precise ephemeris determined by the U.S. Defense Mapping Agency and continuous observations over several days. This positioning mode has low accuracy, long single-point positioning time, and long intervals between positioning times, and currently, only the Argos system still uses this positioning method [14,15]. The second category involves positioning using non-cooperative Signals-of-Opportunity (SOPs) [16,17], requiring the processing of a large number of heterogeneous radio signals and the blind estimation of Doppler frequency shift information. In this positioning mode, the satellite’s position, velocity, clock error, and clock drift information are all very ambiguous, and the satellite clock error usually needs to be estimated as unknown parameters [18,19,20]. Neinavaie et al. [11] achieved a horizontal positioning accuracy of 10 m using observations from a complete six arc of Starlink satellites; Khalife and Kassas [21] assumed that the clock drift of satellites and receivers did not change with time and included them in the state estimation; Wang et al. modeled the orbit errors caused by the clock error of LEO satellites [22], and verified the validity of the model using opportunistic signal data from Starlink satellites. The third category involves instantaneous Doppler positioning using cooperative LEO satellite signals [23,24]. If the satellite broadcast ephemerides and clock error products are retrievable, the satellite’s position, velocity, and clock error information can be considered to be known. Shi et al. [25] modeled errors in an LEO satellite ephemeris and atmospheric errors and used iterative LS estimation to simultaneously estimate the user receiver’s position, velocity, and clock drift, achieving simulation results comparable to GNSS positioning. Guo et al. used Extended Kalman Filtering (EKF) and iterative LS estimation for Doppler velocity determination and positioning, and they conducted a comprehensive performance analysis, with simulation results comparable to GNSS positioning and velocity determination [26]. Mark L. Psiaki et al. [27] proposed a new global satellite navigation scheme based solely on carrier Doppler frequency shift positioning and conducted evaluations, demonstrating its practicality in application. Currently, research on determining receiver states using Doppler frequency shift observations from LEO satellites in dynamic trajectory scenarios is still relatively lacking.

The Doppler positioning LEO satellites fundamentally belong to the Doppler frequency shift (DFS) positioning system based on their frequency information. However, due to the high non-linearity between the DFS and target state, compared to the well-established Time of Arrival (TOA) positioning system, the DFS positioning system exhibits complex geometric curves [28,29,30,31]. Moreover, in dynamic scenarios where the simultaneous estimation of receiver position and velocity is required, the equations exhibit even higher-dimensional complex surfaces [32]. The high level of non-linearity often leads to ill posedness in the equation-solving procedure [33,34]. Shi et al. discovered in their research that LEO Doppler positioning is sensitive to the initial iteration values, as in their simulation scenario, where the epoch solutions failed when the initial value error exceeded 200 km [25]. Guo et al. found that larger errors in the initial iteration values resulted in greater velocity errors when using the EKF method [26]. Therefore, investigating the ill-posedness issue in LEO Doppler positioning and obtaining accurate initial values without relying on additional prior information holds significant importance. Currently, research on this issue remains relatively scarce.

Within the framework of instantaneous Doppler positioning using cooperative LEO satellite signals, this study conducts the first theoretical analysis on the origin of ill posedness in LEO satellite Doppler positioning and proposes a Doppler positioning algorithm that is insensitive to initial values. The algorithm’s performance was confirmed through various simulation tests. The structure of this study is as follows: Section 2 introduces the Doppler positioning model, provides an analytical expression for the linearization error of the positioning equation, proves the existence of local minima or saddle points in the equations, and presents specific cases. This section further details the method for determining positioning velocity in dynamic conditions. In Section 3, the algorithm is validated in static and high dynamic scenarios. Finally, Section 4 offers relevant conclusions.

2. Theory and Algorithm

This section begins with the theoretical foundation for using the Doppler frequency shift of LEO satellites to estimate the receiver’s position, velocity, and clock drift simultaneously. It then examines the inherent challenges of Doppler positioning from two different viewpoints. First, we analyze what affects the linearization error in both LS and extended EKF methods to better understand the challenges of linearization. Next, we approach the problem of estimating the receiver’s state as an optimization challenge, focusing on finding the lowest point of the error function. We also explore whether there are multiple local solutions that are not the best. Finally, we suggest a method based on Tikhonov regularization to quickly find good starting values for the process.

2.1. Doppler Positioning Model

The Doppler effect is caused by the relative velocity in the line-of-sight (LOS) direction between the transmitting and receiving ends. The Doppler observation equation can be expressed as follows:

\begin{matrix} f_{d} & = \pm \frac{v_{los}}{c} \cdot f_{c} \\ = - \frac{f_{c}}{c} [(v^{s} - v_{r}) \cdot \frac{p^{s} - p_{r}}{∥p^{s} - p_{r}∥} + c \cdot (δ {\dot{t}}_{r} - δ {\dot{t}}^{s}) \\ + c \cdot d {\dot{R}}_{r}^{s} + {\dot{T}}_{r}^{s} + {\dot{I}}_{r_{r} f}^{s} + d {\dot{E}}_{r}^{s}] + ε_{f} \end{matrix}

(1)

In the equation,

f_{d}

represents the Doppler frequency shift between the satellite and the receiver,

v_{los}

represents relative velocity in the LOS direction, c represents the speed of light, and

f_{c}

represents carrier frequency.

v^{s}

and

p^{s}

are the satellite velocity and position vectors, while

v_{r}

and

p_{r}

are the velocity and position vectors of the receiver;

δ t_{r}

and

δ {\dot{t}}^{s}

correspond to the clock drifts of the receiver and the satellite, respectively;

d {\dot{R}}_{r}^{s}

is the clock drift correction due to relativistic effects;

{\dot{T}}_{r}^{s}

and

i_{r, f}^{s}

represent the tropospheric and ionospheric delay rate errors, and

d {\dot{E}}_{r}^{s}

accounts for the error induced by the Sagnac effect caused by Earth’s rotation.

ε_{f}

denotes the unmodeled error. Throughout this paper,

∥.∥

denotes the L2 norm of a vector. Within the framework of cooperative signals, it can be assumed that the satellite clock bias, clock drift, atmospheric errors, relativistic effects, and Sagnac effect can be eliminated by broadcasting parameters or error models. In this case, the parameters to be estimated are the position and velocity vectors of the receiver, as well as the clock drift of the receiver. In order to facilitate mathematical analysis, the modified Doppler observation equation is given as follows:

{\hat{f}}_{d i} = f_{i} (s) + ε_{f_{i}}

(2)

f_{i} (s) = - \frac{f_{c}}{c} \frac{{(p_{i}^{s} - K^{T} s)}^{T}}{∥p_{i}^{s} - K^{T} s∥} (v_{i}^{s} - Q^{T} s) + f_{c} J^{T} s

(3)

where

{\hat{f}}_{d i}

represents the Doppler frequency shift observation of the ith satellite after error correction. Accordingly,

f_{i} (s)

represents the theoretical Doppler value.

s = {[p_{x}, p_{y}, p_{z}, v_{x}, v_{y}, v_{z}, δ {\dot{t}}_{r}]}^{T}

is the vector to be estimated consisting of the receiver’s position, velocity, and clock drift.

K^{T}, Q^{T}, J^{T}

are selection matrices defined as follows:

K^{T} = {[\begin{matrix} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]}_{7 \times 7}, Q^{T} = {[\begin{matrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]}_{7 \times 7}

J^{T} = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & - 1 \end{matrix}]

Automatically in Equation (3), the ith satellite position vector

p_{i}^{s}

and velocity vector

v_{i}^{s}

are expanded with zero elements to form seven-dimensional vectors. Actually, in the Hilbert space [35],

f_{i} (s)

can be viewed as a function mapping from seven dimensions to one dimension. Without constraints, it can be denoted as follows:

f_{i} (s) : R^{7} \to R^{1}, s \in R^{7}

(4)

To write the equation in vector form at the observation of M satellites, we have the following:

F_{d} = F_{s} + ε_{f}

(5)

where

\{\begin{matrix} {\hat{F}}_{d} = {[{\hat{f}}_{d 1}, {\hat{f}}_{d 2}, \dots, {\hat{f}}_{d M}]}^{T} \\ F_{s} = {[f_{1} (s), f_{2} (s), \dots, f_{M} (s)]}^{T} \\ ε_{f} = {[ε_{f_{1}}, ε_{f_{2}}, \dots, ε_{f_{M}}]}^{T} \end{matrix}\}

(6)

In accordance with the LS criterion, the optimal state estimation can be equivalently achieved through the minimization of the subsequent objective function:

min {({\hat{F}}_{d} - F_{s})}^{T} ({\hat{F}}_{d} - F_{s})

(7)

The current challenge involves an unconstrained optimization problem that seeks to minimize the norm, which is a common occurrence in positioning issues. For unconstrained optimization problems, it is generally necessary to set an initial value. Then, based on the characteristics of the objective function and the scale of the problem, an appropriate optimization method should be chosen. The optimal solution is obtained through iterative computation [36]. In such contexts, the iterative LS is frequently employed, notably in GNSS pseudorange single-point positioning scenarios. However, this method performs optimally only under mild non-linearity of the equations and is highly sensitive to inaccuracies in the initial values [37]. For the application of LEO satellite Doppler positioning, it is more sensitive to initial point errors compared to GNSS single-point positioning. Therefore, designing an algorithm that can address this challenge is particularly important.

2.2. Analysis of Ill-Posed Problems

A problem is considered well-posed when it adheres to the following three fundamental conditions [38,39]:

Existence: A solution exists for the given dataset.
Uniqueness: The solution is unique for the specified dataset.
Stability: Slight variations in the input data lead to only minor changes in the resultant solution.

Conversely, if any of these criteria are not satisfied, the problem is classified as ill posed. The corresponding examples are shown in Figure 1.

Figure 1 shows examples of objective functions for three ill−posed situations: the absence of an optimal solution (Red line), the non-uniqueness of the solutions (Blue line), and instability (Green line). Referencing Figure 1, we employ mathematical language to describe the three factors of the ill-posed problem: In an optimization problem, let the objective function be

f : D \to R

, where D is the domain:

$Not existent \Rightarrow \exists M \in R s . t . \forall x \in D, f (x) > M .$
$Not unique \Rightarrow \exists x_{1}, x_{2} \in D s . t . f (x_{1}) = f (x_{2}) = m$ .
$Not stable \Rightarrow \forall ϵ > 0, \exists δ > 0 such that \forall x, x^{'} \in D, ∥ x - x^{'} ∥ < δ \Rightarrow ∥ f (x) - f (x^{'}) ∥ \geq ϵ .$

In the realm of Doppler positioning systems, ensuring both the existence and uniqueness of a solution depends on having sufficient observations. The inherent instability of the problem makes it ill posed. The significant non-linearity of the Doppler positioning equation causes substantial deviations in the solution with small perturbations or different algorithm choices, leading to instability. This instability in the Doppler positioning computation process manifests in two distinct phenomena:

Algorithmic Challenges: Widely employed algorithms such as the iterative LS and gradient descent may result in non-invertible iteration matrices and unsuccessful solutions attributed to inadequate initial value selection, consequently inducing multicollinearity. Moreover, filtering techniques like the EKF may exhibit varying solution errors in response to modifications in initial values.
Local Non-Uniqueness: The problem may demonstrate locally non-unique solutions, wherein despite algorithmic convergence, the obtained outcomes may still harbor significant errors.

In tackling these challenges, we first analyze the linearization errors of the Doppler observation equations, which are related to the satellite’s orbital altitude and initial value errors, revealing instances where traditional algorithms may fail to converge in Section 2.2.1. Secondly, from the perspective of optimization theory, we demonstrate the existence of the potential for the solution to converge to locally non-unique solutions in the Doppler positioning problem in Section 2.2.2.

2.2.1. Linearization Error

Linearization stands as a prevalent method employed to tackle non-linear equations effectively. In algorithms such as iterative LS and EKF, linearizations of the equations are compulsory. This procedure entails the Taylor expansion of the non-linear function around a specific point, where solely the first-order term is retained, thereby yielding a linearized equation. This resultant equation delineates the hyperplane equation that tangentially intersects the original equation at the designated point. The disparity between the linearized equation and the original equation, denoted as the linearization error, acts as a gauge of the degree of non-linearity inherent in the equation. Taking the derivative of Equation (3), we obtain the following:

\frac{\partial f_{i} (s)}{\partial s} = [\frac{\partial f_{i} (s)}{\partial p_{r}}, \frac{\partial f_{i} (s)}{\partial v_{r}}, \frac{\partial f_{i} (s)}{\partial {\dot{t}}_{r}}]

(8)

where

\frac{\partial f_{i} (s)}{\partial p_{r}} = - \frac{f_{c}}{c} [e_{i}^{p} 〈e_{i}^{p}, e_{i}^{v}〉 - e_{i}^{v}]

(9)

\frac{\partial f_{i} (s)}{\partial v_{r}} = \frac{f_{c}}{c} e_{i}^{p}

(10)

\frac{\partial f_{i} (s)}{\partial {\dot{t}}_{r}} = - f_{c}

(11)

In the Equations (8) and (9), we define

e_{i}^{p} = \frac{p_{i}^{s} - p_{r}}{∥p_{i}^{s} - p_{r}∥}

as the direction cosine (DC) and

e_{i}^{v} = \frac{(v_{i}^{s} - v_{r})}{∥p_{i}^{s} - p_{r}∥}

as the pseudo direction cosine of velocity (PDCV). Therefore, the linearized equation at the initial value

s_{0}

can be represented as follows:

B_{s_{0}}^{i} (s) = f_{i} (s_{0}) + {\frac{\partial_{i} f (s)}{\partial s}|}_{s = s_{0}} (s - s_{0})

(12)

To facilitate comparison, the expression for the relative linearization error is defined as follows:

f_{e r r o r}^{i} (s_{0}) = |\frac{f_{i} (s_{t}) - B_{s_{0}}^{i} (s_{t})}{f_{i} (s_{t})}|

(13)

where

s_{t}

represents the true state of the receiver. Substituting Equations (3) and (8)–(12) into Equation (13) and rearranging, we can obtain the expression for the relative linearization error:

f_{e r r o r}^{i} (s_{0}) = |\frac{Δ f_{i} - \frac{f_{c}}{c} {[e_{i}^{p} 〈e_{i}^{p}, e_{i}^{v}〉 - e_{i}^{v}, e_{i}^{p}, c]}^{T} Δ s}{f_{i} (s_{t})}|

(14)

In Equation (14),

Δ f_{i} = f_{i} (s_{t}) - f_{i} (s_{0})

represents the difference in the Doppler frequency shift observed by the receiver between the true state and the initial state;

Δ s = s - s_{0}

represents the difference in two state vectors. Upon analysis of Equation (14), it is evident that the values of

Δ f_{i}

Δ s

, and

e_{i}^{v}

are solely influenced by the relative positions or velocities of the satellite and the receiver. The distance from the satellite to the receiver can be approximated as the orbital height, primarily affecting the value of the PDCV

e_{i}^{v}

. Therefore, defining

γ_{i} = \frac{1}{∥p_{i}^{s} - p_{r}∥}

as the distance impact factor, the PDCV can be expressed in the following form:

e_{i}^{v} = γ_{i} (v^{s} - v_{r})

(15)

From this, the following conclusions can be drawn. For LEO satellites, the distance impact factor is significant, amplifying slight changes in the relative velocity and consequently causing drastic variations in the linearization errors. Conversely, for medium- to high-orbit satellites, the distance impact factor is only about one-thousandth of that for LEO satellites, resulting in less drastic changes in linearization errors. This observation partially explains that lower orbits are more prone to ill-conditioning issues. For the convenience of drawing and comparing the results, numerical verification of Equation (14) was performed under the condition of zero receiver speed. Taking the actual receiver position

p_{r}

as the center of a sphere, the distance between the initial value

p_{0}

and the actual receiver position is taken as the radius r, and the hypothetical spherical surface S was constructed:

S = \{x \in R^{3} : {∥x - p_{r}∥}_{2} \leq r\}

(16)

Taking r as 200 km, the point

p_{0}

traverses the spherical surface. Subsequently, multiple sets of experiments were conducted with the satellite orbit altitude as the variable, and the relative linear errors obtained are illustrated in Figure 1.

In Figure 2, the range and scale of the vertical axis of the four sub-figures are consistent to facilitate the comparison of differences between different orbital altitudes; the differences in linearization errors at the same orbital altitude in different directions are represented by a heatmap. It can be observed that as the satellite altitude decreased, the relative linearization error increased significantly and changed dramatically. Moreover, under a fixed initial value error distance, the relative linearization error exhibited strong non-linear characteristics, with pronounced variations in relation to the azimuth of the chosen initial values. As shown in Table 1, the specific details are presented.

It should be noted that the convergence interval of the initial values was fundamentally determined by both the optimization algorithm and the objective function. The linearization error can only serve as a reference indicator for the difficulty of convergence; the specific situation is also related to the visible satellite configuration and the orientation of the initial value selection. Reference [25] provides preliminary simulation verification of this issue, demonstrating that in a mega-constellation scenario, if the initial error exceeds 300 km, there will be epochs that cannot converge. However, the reference did not explore its relationship with the constellation orbit altitude, which will be verified and analyzed in detail later in this paper. Additionally, literature reference [40] employed a grid search method to address the sensitivity of initial points. According to the conclusions of this study, by combining satellite orbit altitude, approximate initial value errors and dynamically designing the grid spacing, it is possible to better adapt to actual conditions and reduce computational complexity.

In dynamic scenarios, the conclusions are similar to static situations, with the additional consideration of initial velocity and clock drift. In general, the lower the satellite orbit altitude and the farther the initial state chosen is from the actual state of the receiver, the more likely the first kind of ill-posed problem will occur, leading to algorithm convergence difficulties.

2.2.2. Existence of Local Non-Unique Solutions

The problem of receiver state estimation, as previously described, can be viewed as an optimization problem. If the objective function is non-convex, it may lead to the existence of local non-unique solutions, causing another type of ill-posed problem in position estimation: the algorithm can converge but with significant errors. In this section, we provide a mathematical proof from a mathematical perspective that the Doppler frequency shift observation equation and positioning objective function is non-convex. We rewrite the relevant definitions as follows [36]:

For the function $f : R^{n} \to R^{1}$ , if $\forall x, y \in dom (f), \forall θ \in [0, 1] \to f (θ x + (1 - θ) y) \leq$ $θ f (x) + (1 - θ) f (y)$ , then the function f is convex.

According to the above definition, we define the discriminant function

g (s_{1}, s_{2})

as follows:

g (s_{1}, s_{2}) = f (θ s_{1} + (1 - θ) s_{2}) - [θ f (s_{1}) + (1 - θ) f (s_{2})]

(17)

According to the definition of a convex function, for the observation equation

f_{i} (s)

in this paper, if

g (s_{1}, s_{2})

is always less than 0 or always greater than 0, then

f_{i} (s)

is a convex function. If the range of

g (s_{1}, s_{2})

includes 0, then

f_{i} (s)

is considered non-convex. For a non-convex objective function, the solution obtained after the algorithm converges is not necessarily the global optimal solution, which means that locally non-unique solutions exist. Our proof approach mainly involves scaling

g (s_{1}, s_{2})

using the Cauchy–Schwarz inequality and the Minkowski inequality [41], as well as proving that its maximum value is greater than 0, and the minimum value is less than 0, thereby proving that

f_{i} (s)

is non-convex. Since the observations from each satellite are independent, it can also be demonstrated that the objective function is non-convex.

Specifically, substituting the aforementioned expression into Equation (4) and applying the Cauchy–Schwarz inequality to constrain the first term on the right-hand side of Equation (18), we derive the following:

- \frac{c}{f_{c}} g (s_{1}, s_{2}) ⩽ ∥v^{s} - Q^{T} (θ s_{1} + (1 - θ) s_{2})∥ - [θ e_{1}^{T} (Δ v_{1}^{T} s_{1}) + (1 - θ) e_{2}^{T} (Δ v_{2}^{T} s_{2})]

(18)

where

e_{1}, e_{2}

are direction cosines, with each component’s absolute value being less than or equal to 1; as well,

Δ v_{i} = v^{s} - Q^{T} s_{i}, i \in N

. Then, by utilizing the Minkowski inequality to bound the first term on the right-hand side of Equation (18) and after rearranging, we obtain the following:

g (s_{1}, s_{2}) \geq - \frac{f_{c}}{c} \{θ [∥Δ v_{1}^{T} s_{1}∥ - e_{1}^{T} (Δ v_{1}^{T} s_{1})] + (1 - θ) [∥Δ v_{2}^{T} s_{2}∥ - e_{2}^{T} (v^{s} - Q^{T} s_{2})]\}

(19)

Let the term on the right-hand side of Equation (19) be denoted as

g {(s_{1}, s_{2})}_{min}

. It is important to highlight that

∥Δ v_{i}∥ ⩾ e^{T} (Δ v_{i})

, thus leading to

g {(s_{1}, s_{2})}_{min} \leq 0

and indicating that

inf [g (s_{1}, s_{2})] \leq 0

. Correspondingly, by utilizing the Cauchy–Schwarz inequality to amplify the second and third terms on the right-hand side of Equation (18), we can obtain the upper bound for the discriminant:

g (s_{1}, s_{2}) \leq - \frac{f_{c}}{c} \{θ [e_{1}^{T} (Δ v_{1}^{T} s_{1}) - ∥v^{s} - Q^{T} s_{1}∥] + (1 - θ) [e_{2}^{T} (Δ v_{2}^{T} s_{2}) - ∥Δ v_{2}^{T} s_{2}∥]\}

(20)

This upper bound is strictly greater than or equal to 0, which implies that

sub [g (s_{1}, s_{2})] \geq 0

. Let the right-hand side term of Equation (18) be denoted as

g {(s_{1}, s_{2})}_{max}

. If we do not consider the singular case where the receiver coincides with the satellite positions, then

g (s_{1}, s_{2})

is continuously differentiable. Therefore, we have

\{g (s_{1}, s_{2}) \in [inf (g (s_{1}, s_{2})), sub (g (s_{1}, s_{2}))] / inf (g (s_{1}, s_{2})) \leq 0, sub (g (s_{1}, s_{2})) \geq 0\}

(21)

which means the range of

g (s_{1}, s_{2})

includes a neighborhood of 0. By definition, when

s \in R^{7}

, the function

f (s)

is non-convex, implying the existence of

\exists s \in N

such that

\nabla^{2} f (s) < 0

and the existence of

\exists s \in P

such that

\nabla^{2} f (s) ⩾

0. Here, N represents the set of state variables where the Hessian Matrix

\nabla^{2} f (s)

f (s)

is negative definite, while P represents the set of state variables where

\nabla^{2} f (s)

is positive semi-definite.

Next, we consider the conditions for equality. Based on the Cauchy–Schwarz and Minkowski equality conditions, let

s_{3} = (θ s_{3} + (1 - θ) s_{3})

. It can be deduced that if there exists a set

L

such that

s_{1}, s_{2}, s_{3} \in L

and it satisfies

∥Δ v_{j}^{T} s_{j}∥ = e_{j}^{T} (Δ v_{j}^{T} s_{j}), j = 1, 2, 3

(22)

then

g {(s_{1}, s_{2})}_{\max} = g {(s_{1}, s_{2})}_{min} = 0

(23)

which implies that

g (s_{1}, s_{2}) = 0

This indicates that when

s \in L f (s)

is affine. From this equality case, it can be observed that within the domain, the relationship between the line-of-sight direction vector from satellite to receiver (LOSDSR) and the relative velocity direction vector from satellite to receiver (RVDSR) has a significant impact on the non-linearity of

f (s)

Therefore, let

cos β_{i} = \frac{{(p^{s} - K^{T} s_{i})}^{T} (v^{s} - Q^{T} s_{i})}{∥p^{s} - K^{T} s_{i}∥ ∥v^{s} - Q^{T} s_{i}∥}

(24)

where the physical significance of

β_{i}

is the angle between the LOSDSR and the VDSR. Substituting Equation (23) into Equation (18) and rearranging, we obtain the following:

g (s_{1}, s_{2}) = - \frac{f_{c}}{c} \{∥Δ v_{1}^{θ} + Δ v_{2}^{1 - θ}∥ cos β_{3} - ∥Δ v_{1}^{θ}∥ cos β_{1} - ∥Δ v_{2}^{1 - θ}∥ cos β_{2}\}

(25)

where

Δ v_{1}^{θ} = θ (v^{s} - Q^{T} s_{1}), Δ v_{2}^{1 - θ} = (1 - θ) (v^{s} - Q^{T} s_{2})

Furthermore,

β_{3}

is bounded between

β_{1}

and

β_{2}, cos β_{3} \leq cos β_{1} + cos β_{2}

. By Minkowski’s inequality, it is known that

∥Δ v_{1}^{θ} + Δ v_{2}^{1 - θ}∥ \leq ∥Δ v_{1}^{θ}∥ + ∥Δ v_{2}^{1 - θ}∥

(26)

always holds true. Therefore, the sign of

g (s_{1}, s_{2})

is entirely determined by

β

. Let the discriminant function be defined as

f^{*} = ∥Δ v_{1}^{θ} + Δ v_{2}^{1 - θ}∥ cos β_{3} - ∥Δ v_{1}^{θ}∥ cos β_{1} - ∥Δ v_{2}^{1 - θ}∥ cos β_{2}

(27)

Consequently, the following inference is given:

\exists β \in N^{*} \Rightarrow f^{*} \geq 0, s \in N, \nabla^{2} f (s) < 0; \exists β \in P^{*} \Rightarrow f^{*} < 0, s \in P, \nabla^{2} f (s) ⩾ 0

For a better intuitive understanding, we assume that the relative velocity direction between the satellite and the receiver is

[1, 1]

in a two-dimensional space. After conducting 200,000 Monte Carlo simulations, the following figure was obtained in the subsequent figure.

In Figure 3, the red point set represents an affine set, serving as a reference, while the cyan point set represents the collection of the objective function. Figure 3a displays a three-dimensional view of the objective function set and the affine set, while Figure 3b provides a top-down view. It is clearly observable from the figures that the affine set does not lie within the interior of the objective function set, and their intersection is not a convex set. According to convex optimization theory, if the intersection of a set and an affine set is not convex, then that set cannot be convex either [36]. Therefore, it can be clearly concluded from the figures that the objective function is non-convex, which further corroborates the previous proof. Returning to the optimization problem in this article, taking the second derivative of the objective function (Equation (8)) with respect to the receiver state s yields the following:

\frac{\partial^{2} {({\hat{F}}_{d} - F_{s})}^{T} ({\hat{F}}_{d} - F_{s})}{\partial s^{2}} = 2 \sum_{i = 1}^{M} \nabla^{2} f_{i} (s)

(28)

The addition of matrices does not affect the positive definiteness of a matrix. Therefore, if there are no constraints on the feasible domain of the receiver state vector, i.e.,

s \in R^{7}

, the objective function of this optimization problem is also non-convex. Similarly, in a two-dimensional space, by constructing five satellites with relative velocity vectors

[- 20; 100], [1, 1], [2, 102], [- 2, - 300]

, and

[200, - 30]

, generating 100,000 random points within the feasible domain for Monte Carlo simulation yields the results shown in Figure 4.

Figure 4a displays a three-dimensional view of the objective function set and the affine set, while Figure 4a provides a top-down view. It can be observed that the objective function exhibits significant non-linearity and contains saddle points. Near the global optimal solution, the solution undergoes large fluctuations with small changes in position, demonstrating that the ill posedness of the Doppler positioning is primarily caused by the instability of the problem.

In conclusion, in the unconstrained scenario, Doppler positioning at LEO satellites is a non-convex problem, and the positive definiteness of the Hessian matrix

\nabla^{2} f (s)

of its equations varies with the angle between the relative velocity direction and the line-of-sight direction. According to convex optimization theory, in cases of non-convex functions, the second-order Hessian matrix may become negative definite, indicating the possibility of plateaus, saddle points, or local minima on the function curve. In plateau regions, where both the first and second derivatives are zero, the function lacks a descent direction, posing a significant challenge for optimization algorithms to converge. Similarly, at local minima or saddle points, where the first derivative is zero, the gradient around that point remains zero, preventing the optimization algorithm from reaching the global optimal solution.

2.3. Method

Based on the previous analysis, the Doppler positioning problem for LEO satellites exhibits ill-conditioned characteristics. When the initial value error is significant, traditional methods may encounter issues with an ill-conditioned iteration matrix. Additionally, in cases of sub-optimal satellite configurations, the algorithm is prone to converging towards locally non-unique solutions. Therefore, the most cost-effective approach to addressing this ill-conditioned problem is to obtain a relatively accurate initial value at the beginning of the positioning process. This section will utilize Tikhonov regularization [42] to address this issue.

2.3.1. Tikhonov Regularization

Tikhonov regularization is a widely used technique in linear regression that incorporates a squared penalty term into the loss function for model parameters. This addition helps control the model’s complexity and reduces the risk of overfitting. Numerous studies have extended the application of Tikhonov regularization to solve ill-posed problems. In the context of LEO satellite Doppler positioning, Tikhonov regularization is beneficial from both a problem-solving and computational perspective.

Firstly, the inherently ill-posed nature of Doppler positioning in LEO satellites can lead to multiple solutions or a broad solution space, often worsened by data noise or incomplete information. Tikhonov regularization, by adding a penalty term, balances data fitting with model complexity, thus narrowing the solution space and enhancing stability and reliability. Secondly, Tikhonov regularization generally has a lower computational burden compared to heuristic optimization methods. As a mathematical optimization tool with relatively simple computational procedures, it is both efficient and easy to implement. In scenarios requiring real-time performance, which is crucial for positioning applications, Tikhonov regularization is a suitable choice, aligning well with the demands of real-time processing.

Considering the specific problem discussed in this article, the process is outlined as follows: Assuming there are N Doppler observation values at a certain epoch, substituting the observation value

{\hat{F}}_{d}

and the initial state

s_{0}

into the linearized Equation (12) yields the following linearized system of equations:

{\hat{F}}_{d} = F_{s} + \nabla F^{T}_{s} Δ s + ε

(29)

where

Δ s = (s - s_{0})

, and

\nabla F_{s} = [\frac{\partial f_{1} (s)}{\partial s_{0}}, \dots, \frac{\partial f_{N} (s)}{\partial s_{0}}]

(30)

Upon adding a Tikhonov regularization term to this optimization objective, we obtain a new optimization objective:

min {∥{\hat{F}}_{d} - F_{s} - \nabla F^{T}_{s} Δ s∥}^{2} + λ {∥ Δ s ∥}^{2}

(31)

where

λ

is the regularization coefficient and is typically referred to as the damping factor, and · and its value has a significant impact on the solution of non-linear problems [43]. The role of the regularization term

{∥ Δ s ∥}^{2}

can be understood as restraining the divergence of

Δ s

during the iterative process. From Equation (30), we obtain the iterative formula in the solution process [44]:

s_{i + 1} = s_{i} + {(\nabla F_{s}^{T} \nabla F_{s} + λ I)}^{- 1} \nabla F_{s}^{T} ({\hat{F}}_{d} - F_{s})

(32)

It is evident that, compared to the commonly used iterative LS estimation in Doppler positioning, this method introduces an additional damping factor term

λ I

, and the damping factor

λ

typically needs to be selected based on actual circumstances [45].

2.3.2. Selection of Damping Factor

Due to the presence of local optimal solutions and the strong non-linearity in the Doppler positioning system equations, this is the primary factor leading to the ill conditioning of the iterative matrix in practical calculations [46]. The condition number is a key criterion for determining whether a matrix is ill conditioned [33]. Considering that the damping factor is essentially used to alleviate the ill conditioning of the iterative matrix, the selection of the damping factor should be related to the condition number of the iterative matrix. For a positive semi-definite real matrix

\nabla {F_{s}}^{T} \nabla F_{s} + λ I

, its condition number is given by the following [47]:

cond (\nabla F_{s}^{T} \nabla F_{s} + λ I) = \frac{λ_{max} + λ}{λ_{min} + λ}

(33)

where

λ_{max}

and

λ_{min}

represent the maximum and minimum eigenvalues of the iterative matrix

\nabla {F_{s}}^{T} \nabla F_{s}

, respectively. Let

C = cond (\nabla F_{s}^{T} \nabla F_{s} + λ I)

. Typically,

C

falls within a certain range where the algorithm converges smoothly [48]. Referring to Equation (32), the formula for computing the damping factor is provided.

λ = |\frac{λ_{max} - C λ_{min}}{C - 1}|

(34)

The value of C is intricately linked to the number of visible satellites, as well as the elevation angle, azimuth angle, and velocity of each visible satellite. Among these factors, the number of visible satellites stands out as a parameter that is relatively straightforward to ascertain and exerts the most substantial influence on the value of C. Consequently, this study conducted an extensive array of simulation experiments, utilizing the number of visible satellites as a variable to unravel the underlying patterns governing C. In these experiments, the initial value errors were uniformly set to zero to ensure the invertibility of the iterative matrices throughout the experiments. The findings revealed that the value of C exhibited fluctuations in each epoch’s computation, with the magnitude of these fluctuations varying between the orders of

\pm 10^{6}

\pm 10^{7}

. By computing the average value of C at each epoch and subsequently fitting this average value with a mathematical function, the results are depicted in Figure 5.

The x axis in Figure 5 represents the number of visible satellites, while the y axis denotes the condition number of the iterative matrix. The graph concurrently displays both the average condition number and the fitting curve. In this study, a model incorporating two exponential terms was selected to fit the data. The expression for the fitting function is as follows:

C = 10^{6} (36.116 e^{(- 0.3063 n)} + 963.69 e^{(- 0.0033 n)})

(35)

The goodness of fit for this model is excellent, with an R-squared value reaching 0.9457. Consequently, given the knowledge of the eigenvalues of the iterative matrix and the visibility of satellites, the damping factor can be determined based on Equations (32) and (33).

2.3.3. Algorithm Flow

Based on the above analysis, the Doppler positioning algorithm using Tikhonov regularization can be divided into several parts, as illustrated in Figure 6.

The specific algorithm pseudocode is as follows:

Algorithm 1: Receiver State

s^{*}

Estimation Algorithm

In the context of Doppler positioning in dynamic scenarios, Algorithm 1 is suitable for applications requiring predefined initial values. For instance, employing the algorithm to compute the initial value

s^{*}

enables its utilization as a fresh starting point in iterative methods such as iteration LS estmation or EKF procedures for resolving subsequent epochs.

3. Simulation Experiment and Result Analysis

3.1. Ill-Posed Simulation Case

To elucidate the two initially introduced ill-posed problems more effectively, this section includes two simulation examples showcasing algorithm non-convergence and the presence of locally non-unique solutions. Subsequently, the algorithm introduced in this study is applied in these simulation scenarios to evaluate its efficacy in addressing ill-posed problems. Since this paper primarily investigates the ill-posed problems in real-time Doppler positioning systems, the minor errors in Doppler frequency shift caused by the ionosphere and troposphere are not the main factors, and therefore, they were not fully considered in the simulation experiments.

3.1.1. Algorithm Non-Convergence

Simulation Parameter Table for Non-Convergence Case Study, as shown in Table 2. And the skyplot of visible satellites for this epoch is shown in Figure 7.

The experimental results are shown in Figure 8.

The x axis in the Figure 8 represents the number of iterations, while the y axis represents the logarithm of the condition number and the logarithm of the loss. The loss function (Loss) denotes the correction term

∥s_{m} - s_{m - 1}∥

during the algorithm iteration. Analysis of the graph reveals that with an increase in the number of iterations, both the condition number and the loss exhibited exponential growth. The condition number also displayed significant fluctuations, rendering it challenging to visualize. This phenomenon is attributed to the high non-linearity observed in the Doppler equations for LEO satellites, as discussed in Section 2. In scenarios with substantial initial value errors, the linearization step in the algorithm introduces notable linearization errors, potentially leading to multicollinearity in the iteration matrix and resulting in a markedly high condition number. This, in turn, complicates inversion, leading to algorithm divergence.

By using the algorithm proposed in this paper, the problem of non-convergence in the previously mentioned case can be effectively reduced. Adding a penalty term to the objective function helps lower the condition number of the iteration matrix and limits the divergence of the loss. The experimental results are provided below.

In Figure 9, it can be observed that the condition number of the iteration matrix has been reduced by approximately

10^{7}

times compared to before. Moreover, after 54 iterations, the Loss reached the convergence stopping criteria. The solution took a total of 0.03 s, with position error of 3.949 m, velocity error of 0.071 m/s, and clock bias error of 3.0632 ×

10^{- 5}

. The initial error level already met the requirements for subsequent solution algorithms. This case demonstrates the superior performance of the algorithm proposed in this paper for ill-posed problems of the first kind. In fact, without employing Tikhonov regularization, the performance of the algorithm will heavily depend on the initial values and satellite orbit heights. The table below shows the percentage of convergent epochs corresponding to different orbit heights and initial error distances in the Walker 220/20/10 constellation configuration. In this scenario, the receiver velocity was set to 0, the constellation orbit inclination was 55 degrees, the simulation duration was 300 s, and the epoch interval was 1 s; the initial position’s azimuth and elevation angles were both 45 degrees. This setup ensured that there were at least seven visible satellites in each epoch.

Table 3 illustrates a significant finding indicating the emergence of non-convergent epochs when the satellite orbit altitude dropped below 1000 km. With a continued decrease in orbit altitude and an escalation in initial value errors, the occurrence of convergent epochs diminished. This observed trend aligns with the linearization errors outlined in Figure 2, providing substantial evidence that linearization errors adversely affect the algorithm’s convergence capability. Hence, it is crucial to be vigilant of the detrimental repercussions stemming from the non-linearity of equations in constellations with lower orbits. Nonetheless, the methodology proposed in this study effectively mitigated this scenario, as evidenced by the computational results and average iterations presented in the accompanying table.

As seen in Table 4, the algorithm proposed in this study exhibited a strong adaptability to variations in satellite orbit height. With the total number of satellites in the constellation remaining constant, it is observed that at an orbit height of 400 km, approximately nine visible satellites were present, requiring 77 iterations for convergence. In contrast, when the orbit height increased to 21,500 km, the number of visible satellites could exceed 80, but this came at the cost of a significantly higher iteration count of 605, and the positioning error was reduced. Thus, it can be concluded that both the number of visible satellites and the satellite orbit height play crucial roles in influencing the positioning error of the algorithm.

In summary, the experiments in this section first analyzed the performance of the method proposed in this paper and the iterative LS method in a certain positioning epoch using the condition number of the iteration matrix and the LOSS as indicators. The results show that when the initial error was large, the condition number and LOSS of the method in this paper gradually decreased and eventually stabilized. However, the iterative LS method struggled to converge to the optimal solution. Subsequently, to better explore under what circumstances the traditional iterative LS algorithm cannot converge, comparative experiments with different orbital heights and initial errors were set up in this section, and the results are consistent with those in Reference [25]. The experiments further illustrate that for real-time Doppler positioning applications of LEO satellites, the sensitivity of traditional algorithms to the initial point is mainly related to the constellation orbit altitude. Finally, we analyzed the impact of orbital height on the algorithm proposed in this paper under the same constellation configuration. The experimental results indicate that as the orbital height increased, the performance of the positioning algorithm in this paper decreased. This also demonstrates that the real-time Doppler positioning system needs to rely on LEO satellites to be implemented.

3.1.2. Local Non-Uniqueness of Solutions

This case study exemplifies a scenario in which the local non-uniqueness of solutions arises during the application of the LS estimation for solving a particular epoch. All simulation parameters remained consistent. The skyplot depicting visible satellites for this epoch is depicted in Figure 10.

The experimental results are shown in Figure 9.

In Figure 11, the left plot displays the Loss curve for the iteration LS estimation, while the right plot corresponds to the algorithm proposed in this study. A comparison between the two reveals that the iteration LS estimation method converged after 16 iterations, whereas the algorithm presented in this study achieved convergence after 45 iterations. Despite the steeper decline in the Loss curve of the iteration LS estimation, it converged to a locally non-unique solution. In contrast, the algorithm introduced in this study successfully converged in the vicinity of the optimal solution. Detailed results are provided in the table below.

The results presented in Table 4 demonstrate that, in the context of LEO satellite Doppler positioning, the algorithm proposed in this study exhibited the capability to avoid falling into locally non-unique solutions. It can also be observed that even with 37 visible satellites, the traditional LS algorithm still converged to non-unique solutions. This indirectly confirms that the geometric constraints of the Doppler positioning system are relatively weak, and more constraints are needed in the solution process. As for the impact of geometric configuration on the calculation of Doppler positioning solutions, further in-depth exploration is required in subsequent studies.

3.2. Dynamic Scene Simulation Experiment

In Section 2.2.2, this study rigorously demonstrated that the angle

β

between the line-of-sight direction and the relative velocity direction in Doppler positioning scenarios directly impacts the non-linearity of the objective function. When the receiver is stationary or moving at a significantly lower speed compared to the satellites, the angle

β

can be considered to be influenced solely by the positions and velocities of the satellites; if the receiver is located on the Earth’s surface, the range of

β

values remains relatively stable. However, in high-dynamic scenarios, it is necessary to consider the influence of the receiver’s velocity on the angle

β

. An extreme hypothetical example is when the velocity vector of the receiver perfectly aligns with the velocity vector of the satellite, resulting in a relative velocity of 0 and no Doppler frequency shift. To validate the robustness of the proposed research algorithm under high-dynamic conditions, the simulation scenario was constructed using the TLE file released by the North American Aerospace Defense Command (NORAD) on 25 January 2024 at 00:05:35. The scenario comprises a mixed constellation of OneWeb, Iridium, and Globalstar satellites, totaling 801 LEO satellites. The range of the dynamic experiment path is as follows: Lat: [28.2213, 28.4579]; Lon: [112.9019, 112.9916]; H(m): [0, 3160]. The relevant parameter settings are presented in Table 5.

In the high-dynamic scenario, a 300 s trajectory was simulated with an average receiver velocity of 123.2953 m/s.The simulated trajectory included smooth and natural turns to more closely approximate real-world scenarios, as shown in the figure below.

From Figure 12, it can be seen that the simulated trajectory included ascent, descent, and turning phases, with the receiver’s velocity being comparable to that of the aerial vehicle. In this simulation scenario, the receiver had an average of 38.27 visible satellites, and the relationship between the number of visible satellites and time is shown in the Figure 13.

The errors in the initial state calculation for position, velocity, and clock drift are presented in Figure 14, Figure 15 and Figure 16.

The RMSE statistics of the experimental results are as follows: The positioning error in the north direction is 2.4688 m, with a velocity determination error of 0.00491 m/s; in the east direction, the positioning error is 5.1715 m, with a velocity determination error of 0.01020 m/s; in the down direction, the positioning error is 3.7457 m, with a velocity determination error of 0.007635 m/s. The clock drift error is 8.2776 ×

10^{- 6}

s/s. The level of error is comparable to that of GNSS single-point positioning. Moreover, in the independent solution calculations with over 300 epochs with the carrier moving rapidly, there were no instances of divergence or convergence to local non-unique solutions. To compare with the algorithm presented in this paper, the initial values were also set to the zero vector and using the iterative LS estimation algorithm. The results are statistically shown in Figure 17 and Figure 18.

From the figures, it can be observed that when the initial error was too large, although each epoch converged, the iterative least squares (LS) estimation still failed to converge to the optimal solution, always wandering around local non-unique solutions. The positioning error was in the order of hundreds of kilometers, and the velocity determination error was in the order of thousands of meters per second. It should be noted that the algorithm in this paper is essentially a type of least squares estimation as well, and the iterative LS algorithm is consistent in estimation performance. Therefore, when the initial values are set close to the true values, the accuracy of the algorithm in this paper is almost identical to that of the iterative LS algorithm, with no significant difference.

4. Conclusions

This work focuses on the cutting-edge issue of utilizing Doppler shift information from LEO mega-constellations for real-time positioning. Approaching the problem from an optimization perspective, it proposes an innovative and straightforward Doppler positioning algorithm that effectively addresses the shortcomings of traditional iterative LS estimation algorithms. This work can serve as a reference for the algorithm design of Doppler positioning receivers with LEO mega-constellations.

Specifically, we conducted a thorough theoretical analysis of the ill-posed problems in Doppler positioning of LEO satellites, as well as the sensitivity of Doppler positioning to initial values as identified in References [25,26]. Furthermore, we proposed an algorithm to tackle this issue. The key contributions of this study are as follows:

We analyzed the mechanism behind the sensitivity of initial values from the perspective of linearization errors in the observation equations, concluding that as satellite orbits lower, the linearization errors exhibit more significant fluctuations, rendering traditional iteration LS estimation algorithms more sensitive to initial values.
We have mathematically demonstrated that, under unconstrained conditions, Doppler positioning may suffer from convergence to locally non-unique solutions, and we have provided corresponding simulation scenarios for illustration.
We have introduced the Tikhonov regularization method to address this ill-posed problem, and we have quantitatively modeled the damping factor (regularization coefficient). This approach presents substantial advantages over iterative least squares estimation.

To achieve the goal of real-time positioning using only the Doppler frequency shift of LEO satellites, future research can be conducted in the following areas:

Due to the significant differences between the Doppler positioning system and the pseudorange positioning system, it is necessary to investigate the impact of the ionosphere, troposphere, and multipath effects on the Doppler frequency shift, as well as to find correction methods that can meet real-time requirements. Data from Doris, mentioned in the introduction, can be utilized for related research.
When the satellite orbit is very low, or even in Doppler positioning applications based on the aircraft ADS-B, the method in this paper will still have the issue of sensitivity to initial points. It is necessary to delve deeper into the theoretical knowledge of convex optimization to establish a more comprehensive Doppler positioning method and theory. Additionally, methods such as graph optimization, which are suitable for non-linear systems, can be attempted to better utilize the information connections between previous and subsequent epochs.

Author Contributions

Z.X. conceived the idea and designed the experiments with Z.L., X.L., Z.J., Q.W., H.L. and C.W. Z.X. wrote the main manuscript. Z.L. and X.L. reviewed the paper. All components of this research were carried out under the supervision of X.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China under Grant U20A0193.

Data Availability Statement

The data related to this article are available upon request to the corresponding authors.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Zhang, J.; Cai, Y.; Xue, C.; Xue, Z.; Cai, H. LEO Mega Constellations: Review of Development, Impact, Surveillance, and Governance. Space Sci. Technol. 2022, 2022, 9865174. [Google Scholar] [CrossRef]
Blanchard, W. Would a GNSS need a backup? J. Glob. Position. Syst. 2004, 3, 308–321. [Google Scholar] [CrossRef]
Kassas, Z.Z.M. Navigation from low-earth orbit. In Position, Navigation, and Timing Technologies in the 21st Century: Integrated Satellite Navigation, Sensor Systems, and Civil Applications; IEEE: Piscataway, NJ, USA, 2020. [Google Scholar]
Selvan, K.; Siemuri, A.; Prol, F.S.; Välisuo, P.; Bhuiyan, M.Z.H.; Kuusniemi, H. Precise orbit determination of LEO satellites: A systematic review. GPS Solut. 2023, 27, 178. [Google Scholar] [CrossRef]
Deng, L. Research on Passive Localization Based on Doppler Frequency Shift of Multiple Moving Receivers. Ph.D. Thesis, University of Electronic Science and Technology of China, Chengdu, China, 2019. [Google Scholar]
Houghton, P.A. III—The Future Development of Doppler Navigation. J. Navig. 1958, 11, 130–137. [Google Scholar] [CrossRef]
Moorhen, C.N. II—The Navigational Applications of Doppler Equipments. J. Navig. 1958, 11, 125–130. [Google Scholar] [CrossRef]
Song, C.; Wang, H.; Xie, S. Satellite Doppler Positioning Measurements; Surveying and Mapping Press: Beijing, China, 1987. [Google Scholar]
Lee, J.N.; Jang, Y.-G.; Lee, B.H.; Mun, B.Y. Animal Tracking System Using the Doppler Effect for Single LEO Satellite. J. Korean Soc. Aeronaut. Space Sci. 2006, 34, 61–69. [Google Scholar]
Saunier, J. The DORIS network: Advances achieved in the last fifteen years. Adv. Space Res. 2022, 72, 3–22. [Google Scholar] [CrossRef]
Neinavaie, M.; Khalife, J.; Kassas, Z.M. Acquisition, Doppler Tracking, and Positioning with Starlink LEO Satellites: First Results. IEEE Trans. Aerosp. Electron. Syst. 2022, 58, 2606–2610. [Google Scholar] [CrossRef]
Kershner, R.B.; Newton, R.R. The Transit System. J. Navig. 1962, 15, 129–144. [Google Scholar] [CrossRef]
Black, H.D. Early development of Transit, the Navy navigation satellite system. J. Guid. Control Dyn. 1990, 13, 577–585. [Google Scholar] [CrossRef]
Rozylowicz, L.; Bodescu, F.P.; Ciocanea, C.M.; Gavrilidis, A.A.; Manolache, S.; Matache, M.L.; Miu, I.V.; Moale, I.C.; Nita, A.; Popescu, V.D. Empirical analysis and modeling of Argos Doppler location errors in Romania. PeerJ 2019, 7, e6362. [Google Scholar] [CrossRef] [PubMed]
Levanon, N.; Ben-Zaken, M. Random Error in ARGOS and SARSAT Satellite Positioning Systems. IEEE Trans. Aerosp. Electron. Syst. 1985, AES-21, 783–790. [Google Scholar] [CrossRef]
Garrison, J.L.; Nold, B.; Masters, D.; Brown, C.; Bridgeman, J.; Mansell, J.; Vega, M.; Bindlish, R.; Piepmeier, J.R.; Babu, S.R. A Spaceborne Demonstration of P-Band Signals-of-Opportunity (SoOp) Reflectometry. IEEE Geosci. Remote Sens. Lett. 2023, 20, 3507205. [Google Scholar] [CrossRef]
Florio, A.; Bnilam, N.; Talarico, C.; Crosta, P.; Avitabile, G.; Coviello, G. LEO-Based Coarse Positioning Through Angle-of-Arrival Estimation of Signals of Opportunity. IEEE Access 2024, 12, 17446–17459. [Google Scholar] [CrossRef]
Tan, Z.; Qin, H.; Cong, L.; Zhao, C. Positioning Using IRIDIUM Satellite Signals of Opportunity in Weak Signal Environment. Electronics 2019, 9, 37. [Google Scholar] [CrossRef]
Jiang, M.; Qin, H.; Su, Y.; Li, F.; Mao, J. A Design of Differential-Low Earth Orbit Opportunistically Enhanced GNSS (D-LoeGNSS) Navigation Framework. Remote Sens. 2023, 15, 2136. [Google Scholar] [CrossRef]
Kassas, Z.M.; Khairallah, N.; Kozhaya, S. Ad Astra: Simultaneous Tracking and Navigation With Megaconstellation LEO Satellites. IEEE Aerosp. Electron. Syst. Mag. 2024, 1–19. [Google Scholar] [CrossRef]
Kassas, J.K.N.M. Blind Doppler Tracking from OFDM Signals Transmitted by Broadband LEO Satellites. In Proceedings of the 2021 IEEE 93rd Vehicular Technology Conference (VTC2021-Spring), Helsinki, Finland, 25–28 April 2021. [Google Scholar]
Wang, D.; Qin, H.; Huang, Z. Doppler Positioning of LEO Satellites Based on Orbit Error Compensation and Weighting. IEEE Trans. Instrum. Meas. 2023, 72, 5502911. [Google Scholar] [CrossRef]
Wang, W.; Lu, Z.; Tian, Y.; Bian, L.; Wang, G.; Zhang, L. Doppler-Aided Positioning for Fused LEO Navigation Systems. Aerospace 2023, 10, 864. [Google Scholar] [CrossRef]
Xv, H.; Sun, Y.; Zhao, Y.; Peng, M.; Zhang, S. Joint Beam Scheduling and Beamforming Design for Cooperative Positioning in Multi-beam LEO Satellite Networks. IEEE Trans. Veh. Technol. 2024, 73, 5276–5287. [Google Scholar] [CrossRef]
Shi, C.; Zhang, Y.; Li, Z. Revisiting Doppler positioning performance with LEO satellites. GPS Solut. 2023, 27, 126. [Google Scholar] [CrossRef]
Guo, F.; Yang, Y.; Ma, F.; Zhu, Y.; Liu, H.; Zhang, X. Instantaneous velocity determination and positioning using Doppler shift from a LEO constellation. Satell. Navig. 2023, 4, 9. [Google Scholar] [CrossRef]
Psiaki, M.L. Navigation using carrier Doppler shift from a LEO constellation: TRANSIT on steroids. Navigation 2021, 68, 621–641. [Google Scholar] [CrossRef]
Zhang, B.; Shi, Z. Simulation of FDOA Locating Systems. J. Electron. Opt. Control. 2009, 16, 13–16. [Google Scholar] [CrossRef]
Weinstein, E. Measurement of the differential Doppler shift. IEEE Trans. Acoust. Speech Signal Process. 1982, 30, 112–117. [Google Scholar] [CrossRef]
Tirer, T.; Weiss, A.J. High resolution localization of narrowband radio emitters based on doppler frequency shifts. Signal Process. 2017, 141, 288–298. [Google Scholar] [CrossRef]
Chestnut, P.C. Emitter Location Accuracy Using TDOA and Differential Doppler. IEEE Trans. Aerosp. Electron. Syst. 1982, AES-18, 214–218. [Google Scholar] [CrossRef]
Benzerrouk, H.; Nguyen, Q.; Xiaoxing, F.; Amrhar, A.; Rasaee, A.; Landry, R., Jr. LEO Satellites Based Doppler Positioning Using Distributed Nonlinear Estimation. IFAC PapersOnLine 2019, 52, 496–501. [Google Scholar] [CrossRef]
Tang, L. Research on the Ill-Posed and Solving Methods of Nonlinear Least Squares Problem. Ph.D. Thesis, Central South University of China, Changsha, China, 2011. [Google Scholar]
Deya, A. On ill-posedness of nonlinear stochastic wave equations driven by rough noise. Stoch. Process. Their Appl. 2022, 150, 215–249. [Google Scholar] [CrossRef]
Mathé, P.; Hofmann, B. Tractability of linear ill-posed problems in Hilbert space. J. Complex. 2024, 84, 101867. [Google Scholar] [CrossRef]
Moklyachuk, M. Convex Optimization; Wiley: Hoboken, NJ, USA, 2020. [Google Scholar]
Bakushinsky, A.; Smirnova, A. A study of frozen iteratively regularized Gauss–Newton algorithm for nonlinear ill-posed problems under generalized normal solvability condition. J. Inverse Ill-Posed Probl. 2020, 28, 275–286. [Google Scholar] [CrossRef]
Nicholson, O.J.M. What can be estimated? Identifiability, estimability, causal inference and ill-posed inverse problems. arXiv 2019, arXiv:1904.02826. [Google Scholar]
Latz, J. Bayesian Inverse Problems Are Usually Well-Posed. SIAM Rev. 2023, 65, 831–865. [Google Scholar] [CrossRef]
Bianco, G.M.; Giuliano, R.; Mazzenga, F.; Marrocco, G. Multi-Slope Path Loss and Position Estimation With Grid Search and Experimental Results. IEEE Trans. Signal Inf. Process. Over Netw. 2021, 7, 551–561. [Google Scholar] [CrossRef]
Willem, M. Functional Analysis; Springer Nature: Berlin/Heidelberg, Germany, 2022. [Google Scholar]
Tikhonov, A.N. On the solution of ill-posed problems and the method of regularization. Dokl. Akad. Nauk. 1963, 151, 501–504. [Google Scholar]
Fischer, A.; Cellmer, S.; Nowel, K. Assessment of the double-parameter iterative Tikhonov regularization for single-epoch measurement model-based precise GNSS positioning. Measurement 2023, 218, 113251. [Google Scholar] [CrossRef]
Rastogi, A. Nonlinear Tikhonov regularization in Hilbert scales for inverse learning. J. Complex. 2024, 82, 101824. [Google Scholar] [CrossRef]
Umar, A.O.; Sulaiman, I.M.; Mamat, M.; Waziri, M.Y.; Zamri, N. On damping parameters of Levenberg-Marquardt algorithm for nonlinear least square problems. J. Phys. Conf. Ser. 2021, 1734, 012018. [Google Scholar] [CrossRef]
Eriksson, J.; Wedin, P.A.; Gulliksson, M.E.; Söderkvist, I. Regularization Methods for Uniformly Rank-Deficient Nonlinear Least-Squares Problems. J. Optim. Theory Appl. 2005, 127, 1–26. [Google Scholar] [CrossRef]
Chacha, C.S.; Naqvi, S.M.R.S. Condition Numbers of the Nonlinear Matrix Equation X p − A ∗ e X A = I. J. Funct. Spaces 2018, 2018, 3291867. [Google Scholar] [CrossRef]
Vanderschel, D.J. A Theory of Approximate Inverses for the Solution of Matrix Equations by Iteration. Ph.D. Thesis, Rice University, Houston, TX, USA, 1970. [Google Scholar]

Figure 1. Schematic diagram of an ill−posed problem.

Figure 2. The sub-figures (a–d) show the linearization errors for initial values set on a spherical surface with the true value as the center and a radius of 200 km at orbit altitudes of 200 km, 400 km, 600 km, and 800 km. The horizontal and vertical coordinates represent the elevation and azimuth angles of the initial values relative to the true position in the ECEF coordinate system.

Figure 3. The schematic diagram of the Doppler equation in two−dimensional space, where the x and y axes represent the positions of the receiver in meters, and the z axis represents the Doppler frequency shift at different positions. Sub-figure (a) displays a three-dimensional view of the objective function set and the affine set, while sub-figure (b) provides a top-down view.

Figure 4. The schematic diagram of the objective function for Doppler positioning in a two−dimensional space. Sub-figure (a) displays a three-dimensional view of the objective function set and the affine set, while sub-figure (b) provides a top-down view.

Figure 5. Relationship between C value and the number of visible satellites.

Figure 6. Algorithm flowchart.

Figure 7. The skyplot of visible satellites for the non-convergence case study.

Figure 8. The convergence status of the algorithm using the LS estimation iteration is shown in the figure above. The left subplot displays the condition number of the iteration matrix, while the right subplot illustrates the loss function during the iteration process.

Figure 9. The convergence status of the algorithm using the algorithm proposed in this paper is shown in the figure above. The left subplot displays the condition number of the iteration matrix, while the right subplot illustrates the loss function during the iteration process.

Figure 10. The skyplot of a case study involving local non-uniqueness of solutions.

Figure 11. The relationship between the logarithmic loss and the number of iterations is depicted in the figure above. The left subplot corresponds to the algorithm proposed in this paper, while the right subplot corresponds to the iteration LS estimation.

Figure 12. Simulation trajectory.

Figure 13. The graph of the variation of the number of visible satellites over time.

Figure 14. The positioning error calculated using the algorithm presented in this paper.

Figure 15. The velocity error computed using the algorithm proposed in this paper in m/s.

Figure 16. The clock drift error computed using the algorithm proposed in this paper in s/s.

Figure 17. The positioning error computed using the iterative LS estimation in s/s.

Figure 18. The velocity error computed using the iterative LS estimation in s/s.

Table 1. The convergence rate of epochs in the Walker 220/20/10 scenario.

CR(%)	3000	2000	1000	500	260	100
OA	3000	2000	1000	500	260	100
400	0	0	0	0	0.00	1
600	0.117	0.503	0.893	1	1	1
800	0.583	1	1	1	1	1
1000	1	1	1	1	1	1
21,500	1	1	1	1	1	1

OA means orbit altitude in km; CR means convergence rate; IE means initial error in km.

Table 2. Simulation parameter table for non-convergence case study.

Parameter	Value
Angle mask	$5^{\circ}$
Frequency measurement error	0.1 Hz
Satellite position variance	5 m
Satellite velocity variance	0.1 m/s
Receiver clock drift	1 × $10^{- 6}$
Initial position	0, 0, 0
Initial velocity and clock drift	0
Satellite orbit altitude	400 km
Number of visible satellites	14
Satellite carrier frequency	11.35 × $10^{9}$ Hz
Receiver coordinates (LLA)	28.2213, 112.9913, 20 m
Convergence threshold	1 × $10^{- 3}$

Table 3. The computation results of the algorithm proposed in this paper in the Walker 220/20/10 scenario.

Orbit Altitude (km)	Position Error (km)	Velocity Error (m/s)	Clock Drift Error (s/s)	Iteration Count
400	3.8903	0.0113	7.47 × $10^{- 6}$	77
600	3.5119	0.0087	7.99 × $10^{- 6}$	20
800	4.1474	0.0085	8.47 × $10^{- 6}$	32
1000	3.8119	0.0082	9.01 × $10^{- 6}$	40
21,500	7.5207	0.0073	7.84 × $10^{- 6}$	605

Table 4. Comparison table of solution results.

RMSE Error Terms	Clock Drift s/s	North Position (m)	East Position (m)	Down Position (m)	North Velocity (m)	East Velocity (m)	Down Velocity (m)
Iteration LS	1.34 × $10^{- 5}$	45,316.7	71,039.1	3,668,618.9	694.3	312.8	1813.1
PA	6.60 × $10^{- 6}$	3.01	7.65	4.90	0.030	0.062	0.045

PA means proposed algorithm.

Table 5. Parameters for the dynamic simulation section.

Parameter	Value
Angle mask	$5^{\circ}$
Frequency measurement error	0.1 Hz
Satellite position variance	5 m
Satellite velocity variance	0.1 m/s
Receiver clock drift	Modeled as white noise with variance 1 × $10^{- 6}$
Average receiver velocity	123.2953 m/s

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Xu, Z.; Li, Z.; Liu, X.; Ji, Z.; Wu, Q.; Liu, H.; Wen, C. Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm. Remote Sens. 2024, 16, 2958. https://doi.org/10.3390/rs16162958

AMA Style

Xu Z, Li Z, Liu X, Ji Z, Wu Q, Liu H, Wen C. Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm. Remote Sensing. 2024; 16(16):2958. https://doi.org/10.3390/rs16162958

Chicago/Turabian Style

Xu, Zichen, Zongnan Li, Xiaohui Liu, Zhimin Ji, Qianqian Wu, Hao Liu, and Chao Wen. 2024. "Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm" Remote Sensing 16, no. 16: 2958. https://doi.org/10.3390/rs16162958

APA Style

Xu, Z., Li, Z., Liu, X., Ji, Z., Wu, Q., Liu, H., & Wen, C. (2024). Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm. Remote Sensing, 16(16), 2958. https://doi.org/10.3390/rs16162958

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Doppler Positioning with LEO Mega-Constellation: Equation Properties and Improved Algorithm

Abstract

1. Introduction

2. Theory and Algorithm

2.1. Doppler Positioning Model

2.2. Analysis of Ill-Posed Problems

2.2.1. Linearization Error

2.2.2. Existence of Local Non-Unique Solutions

2.3. Method

2.3.1. Tikhonov Regularization

2.3.2. Selection of Damping Factor

2.3.3. Algorithm Flow

3. Simulation Experiment and Result Analysis

3.1. Ill-Posed Simulation Case

3.1.1. Algorithm Non-Convergence

3.1.2. Local Non-Uniqueness of Solutions

3.2. Dynamic Scene Simulation Experiment

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI