Non-Cooperative LEO Satellite Orbit Determination Using Single Station for Space-Based Opportunistic Positioning

: Space-based opportunistic positioning is a crucial component of resilient positioning, navigation, and timing (PNT) systems, and it requires the acquisition of orbit information for non-cooperative low Earth orbit (LEO) satellites. Traditional methods for orbit determination (OD) of non-cooperative LEO satellites have difficulty in achieving a balance between reliability, hardware costs, and availability duration. To address these challenges, this study proposes a framework for single-station orbit determination of non-cooperative LEO satellites. By utilizing signals of opportunity (SOPs) captured by a single ground station, the system performs initial orbit determination (IOD), precise orbit determination (POD), and orbit prediction (OP), enabling the long-term determination of satellite positions and velocities. Under the proposed framework, the reliability and real-time performance are dependent on the initial orbit determination and the orbit calculation based on the dynamical model. To achieve initial orbit determination, a three-step algorithm is designed. (1


Introduction
Navigation systems provide positioning, navigation, and timing (PNT) services for various applications, playing an important role and having a wide range of applications in various fields [1].Currently, common PNT systems include systems based on wireless navigation signals and inertial navigation systems.For the former, user terminals perform positioning based on dedicated navigation signals and known spatiotemporal information of the signal radiation sources.This requires the establishment of navigation systems such as global navigation satellite system (GNSS), long range navigation (LORAN), and Locata, which have issues such as high infrastructure investment costs, susceptibility to interference, deception, and even destruction.The latter maintains calibrated position information through extrapolation but has drawbacks such as accumulating errors over time and requiring initial alignment for a long time.Opportunistic positioning is a technology that uses non-navigation/non-cooperative signals that can be received to extract positioning information.Opportunistic positioning offers several advantages, including its independence from additional infrastructure construction, resilience against interference, and the absence of error accumulation over time.As such, it effectively compensates for the shortcomings of current commonly used PNT systems [2,3].Opportunistic positioning is divided into terrestrial-based and space-based, with the latter having strong coverage capabilities, enabling seamless global coverage of desert, ocean, polar, and other regions.Currently, there are thousands of low Earth orbit (LEO) satellites in operation, such as Starlink, OneWeb, Iridium, Orbcomm, and Globalstar.Many companies such as SpaceX, Amazon, Samsung, and Boeing plan to launch mega-constellations consisting of tens of thousands of LEO satellites, providing a huge radiation source for space-based opportunistic positioning.Therefore, space-based opportunistic positioning has become a research hotspot in recent years [1,4,5].For the sake of brevity, the opportunistic positioning in the following text specifically refers to space-based opportunistic positioning.
However, opportunistic positioning faces two challenges [6].On the one hand, the architectures of signals of opportunity (SOPs) are often partially or completely undisclosed, making the extraction of observation difficult.Due to its minimal requirements on signal architecture, Doppler is the most common measurement in the field of opportunistic positioning.On the other hand, unlike navigation satellites, such as GNSS satellites, that modulate ephemeris in their signal messages, opportunistic positioning terminals have difficulty in acquiring real-time accurate information on the sources' position and velocity.According to the existing public literature, the methods for obtaining the orbit of noncooperative LEO satellites as radiation sources for opportunistic positioning mostly rely on the two-line element (TLE) published by the North American Aerospace Defense Command (NORAD) [3,7,8].The TLE provides mean orbital elements at an instant by averaging out periodic perturbations and is generated by the orbit determination (OD) system processing approximately 24 h of observations from a globally distributed space surveillance network (SSN) [9,10].It should be noted that the simplified general perturbations 4 (SGP4) model is compatible with TLE, and using SGP4 to propagate the TLE, the position and velocity of the LEO satellite at moments of interest can be obtained [11,12], with the position and velocity accuracy on the order of kilometers and meters per second, respectively [13].
Some scholars have devised various frameworks for non-cooperative LEO satellite orbit determination by integrating TLE with other information.Kassas et al. improved the opportunistic positioning performance by optimizing the position model [14][15][16], LEO ephemeris [6], and orbit propagation model [17].Li et al. combined TLE with the extracted azimuth, elevation, and range measurements from radar, and determined the orbit of Starlink satellites for two days with accuracy at kilometer grade [18].Sakamoto et al. proposed a method for LEO orbit determination using a radio interferometer of small-diameter antennas, where TLE was used for initialization [19].However, in the aforementioned studies, TLE serves as a prerequisite for achieving the objectives.Deng et al. achieved noncooperative LEO satellite orbit determination independent of TLE using several distributed facilities equipped with the SOPs receivers [20].Nevertheless, there are two issues.On the one hand, the hardware of the orbit determination system is complex, requiring temporal and spatial synchronization among stations and large amounts of information transmission.On the other hand, the available time for the obtained orbit results is short, as only about 10 min of data were utilized to estimate an instantaneous orbit, making it difficult to be used for long-term opportunistic positioning.In summary, the existing non-cooperative LEO satellite orbit determination methods are unable to simultaneously achieve optimal performance in terms of reliability, hardware cost, and availability duration.
Initial orbit determination (IOD) is critical for the reliability of the OD system, which is used for the initialization of precise orbit determination (POD) and has a decisive impact on the convergence of the POD process.IOD based on Doppler observations has been intensively investigated.The article published by Patton R. B. includes several relevant statements advocating a "systematic" exploration of a five-parameter state space [21].Li Wenhua et al. proposed a method for single-station single-pass orbit determination employing a random search approach on the six-dimensional orbital elements [22].Nevertheless, these search methods in high-dimensional state spaces entail long computational time.Christian et al. obtained initial orbits through a one-dimensional search approach by exploiting the constraints of circular satellite orbits [23].This approach implicitly assumes that the satellite's eccentricity is zero and neglects the effect of Earth's rotation, making it challenging to ensure the accuracy of the initial orbit.Therefore, it is highly necessary to acquire the initial orbit in a prompt and precise manner.
The part that determines the real-time performance of the OD system is orbit calculation, which refers to obtaining orbital information at the past/future single orbit/multiple time of interest based on instantaneous orbit by utilizing the satellite dynamics models in this study.In most cases, the computation of multiple past orbits is employed for orbit determination, whereas the computation of a single future orbit is used for orbit prediction (OP).With the launch of the mega-constellations, there will be tens of thousands of non-cooperative LEO satellites in the future.In order to obtain real-time orbits of all satellites without increasing parallel computing costs, the computational efficiency of orbit calculation urgently needs to be improved.However, there is a lack of comprehensive information regarding orbit calculation in the field of non-cooperative LEO satellite orbit determination [19,20,24].
The study makes the following contributions.First, a single ground station-based framework is presented for the non-cooperative LEO satellite orbit determination, enabling long-term reliable orbit estimation using a cost-effective hardware platform.Second, within this framework, a three-step initial orbit determination algorithm is proposed to initialize the precise orbit determination process.Finally, two fast orbit calculation algorithms are developed to calculate orbits with low computational complexity.
The remaining sections of the article are organized as follows.Section 2 presents a framework for single-station orbit determination for the non-cooperative LEO satellite for opportunistic positioning.In Section 3, the initial orbit determination is described.Section 4 demonstrates the fast orbit calculation method.Section 5 verifies the effectiveness of the proposed methods by experiments through practical opportunistic positioning experiments.Finally, Section 6 concludes this article.

Overview
In this subsection, we propose a new framework for opportunistic positioning, which is illustrated in Figure 1.
The proposed framework involves three key elements: non-cooperative LEO satellites, a ground station, and a user.Their roles in opportunistic positioning are as follows: non-cooperative LEO satellites serve as radiation sources; a ground station is utilized to acquire precise orbits of the radiation sources; a user performs self-positioning based on observations extracted from signals transmitted from radiation sources and precise orbits transmitted from a ground station.The distinction between a ground station and a user lies in their objectives and outputs.In terms of objectives, a single ground station is dedicated to long-term monitoring of non-cooperative LEO satellites and obtaining their high-precision orbits.The objective of a user is self-positioning.Regarding outputs, the ground station provides continuously updated precise orbits, whereas a user outputs its own position coordinates.
In this proposed framework, the individual target satellites are denoted as satellite 1, . . ., satellite m, . . ., and satellite M, where 1 < m < M, and M is the number of satellites.
A sole ground station is responsible for monitoring the aforementioned satellites, and the scenario in which the satellites are observable from the ground station is represented by the green area.As a consequence of Earth's rotation and the dynamic behavior of the satellites, the ground station receives multiple passes of SOPs transmitted from the satellites.These multi-pass SOPs from satellite m are represented as 1, . . ., Ψ α m .By leveraging the received SOPs, it becomes possible to acquire measurements from multiple passes, thus enabling the iterative determination of the orbital parameters for each satellite.Let k be the number of times the orbit determination has been executed.After updating the kth precise orbit of satellite m, represented as X P,m k , it is disseminated to the user terminal through broadcasting, where the superscript P indicates POD.
When the satellites enter the user's field of visibility, as depicted by the gray area, the positioning algorithm is employed.On one hand, measurements are extracted from the SOPs.On the other hand, the prediction orbits at the measurement epochs are derived based on the most recent precise orbits.Subsequently, the estimation of position coordinates is achieved through the solution of the Doppler position equation, utilizing the prediction orbit and Doppler measurements as known information.It should be emphasized that in cases where a satellite is simultaneously within the field of view of both the user and the ground station, the signals during that pass are utilized not only for positioning purposes but also for the precise determination of the satellite's orbit.

System Composition and Process of Orbit Determination
In light of the absence of differentiation in the system composition and process of orbit determination among all satellites, the system composition and process for satellite m is employed as an illustrative example.This system operates in a closed-loop feedback manner and can be divided into three distinct components: initial orbit determination, precise orbit determination, and orbit prediction.
The detailed process for orbit determination is elucidated by highlighting crucial concepts within the proposed framework, as depicted in In the first POD-OP iteration, the observations for precise orbit determination ỸP,m 1 are acquired through data association utilizing the initial orbit X I,m 0 (as represented by the purple box).The initial values of the orbital state for initialization, X I,m 1 , are derived by propagating the initial orbit X I,m 0 to the time t P,m 1 .By employing the observational data ỸP,m 1 to refine the initial values X I,m 1 , a precise orbit X P,m 1 at the desired time t P,m 1 is obtained (as depicted by the purple dashed box and curve).Before obtaining the second precise orbit X P,m 2 , the orbits at the desired epoch t P,m 1 ≤ t < t P,m 2 are calculated through orbit prediction, utilizing X P,m 1 as the initial values (illustrated by the red solid curve).Subsequently, upon the availability of a new pass, data association is conducted to incorporate the measurements obtained from this new pass.The observations from the Doppler in the new pass, together with the accumulated Doppler, are utilized as observations for the second iteration of precise orbit determination, denoted as ỸP,m 2 .The observations ỸP,m 2 are employed to rectify the initial values X I,m 2 , which represent the prediction orbit derived from the first precise orbit X P,m 1 .Consequently, the second precise orbit X P,m 2 is obtained, as illustrated in the golden-colored section of the diagram.The spatial and temporal information of the satellite at the desired time t P,m 2 ≤ t < t P,m 3 can be obtained through orbit prediction, and it is the current latest precise orbit X P,m 2 that serves as the initial orbital information X P,m 2 , represented by the red solid curve.The aforementioned process of POD-OP is iteratively repeated, with the nth iteration serving as an exemplification.The orbits at the desired epochs through the interval, denoted by t P,m n ≤ t < t P,m n+1 , can be obtained using the current latest precise orbit X P,m n as the initial value.Here, X P,m n represents the corrected orbit derived from X I,m n using the observation ỸP,m n .and t P,m n+1 is the time corresponding to X P,m n+1 which is the precise orbit to be determined at the next iteration.This process is shown by the brown box and curve.
It is worth noting that, in this paper, the duration of collecting Doppler for precise orbit determination is approximately 24 h.The state vector includes the orbit and frequency offset.The parameter estimation algorithm employed is least squares (LS), which is the same as the method used in the third step of the initial orbit determination.For the sake of simplicity, further details regarding the precise orbit determination will not be elaborated on in the subsequent sections.
The characteristics of the proposed framework are summarized as follows.First, the observation information for the orbit determination system is Doppler extracted by a single ground station from the SOPs.This approach facilitates the attainment of longterm and dependable orbit solutions, while keeping hardware costs low.Furthermore, this eliminates the requirement for temporal and spatial synchronization or information transmission among multiple stations.Second, the initialization information for the precise orbit determination is obtained from autonomously acquired initial orbit determination results, thereby ensuring the reliability of the system.Finally, the orbits at the desired time are determined by continuously predicting the orbit using the latest precise orbit as the initial value, guaranteeing the performance of long-term real-time orbit estimation through the system architecture of closed-loop feedback and the process of iterative POD-OP.

Initial Orbit Determination
In order to address the issue of obtaining a reliably accurate orbit for initialization of the precise orbit determination, the initial orbit determination problem is described, and the three-step initial orbit determination algorithm is presented.
Since the orbit determination system employs the same model and algorithm for each satellite in the constellation, the superscript "m" denoting the specific satellite will be omitted in the subsequent text.

Problem Description
The ground station continuously receives the SOPs and extracts the Doppler measurements.The Doppler measurement at the kth epoch is denoted as Y(t k ).The combination of the satellite's position and velocity in the Earth-centered Earth-fixed (ECEF) coordinate system at time t is denoted as , where (•) T denotes the transpose operation, and the superscript s and the subscript F represent the abbreviation of the satellite and ECEF, respectively.Let the frequency offset of the tracking station and satellite at a given epoch be denoted as δ g (t) and δ s (t), respectively.The relative frequency offset is represented as δ g,s (t), where δ g,s (t) = δ g (t) − δ s (t).Let the state vector composed of X s F and δ g,s be denoted as X OD .Neglecting the small Doppler errors caused by the ionosphere and troposphere [25], the relation h k (•) between the Doppler measurement Y(t k ) at a single epoch t k and the state vector X OD is modeled as where ν k denotes the measurement noise, and R g F represents the position of the ground station in the ECEF frame.The terms c and λ are the speed of light and the wavelength of the signal, whereas ∥•∥ 2 indicates the operation of taking the two-norm of a vector.
The relation between orbits at different time can be expressed as follows: where the position R s ECI (t) and velocity V s ECI (t), both in Earth-centered inertial (ECI) frame, are combined to form the orbital information X s ECI (t) = [(R s ECI (t)) T , (V s ECI (t)) T ] T , and the subscript ECI denotes the abbreviation for the ECI coordinate system.It is worth noting that, unlike X s F (t) in the ECEF coordinate system, X s ECI (t) corresponds to the position and velocity in the ECI coordinate system.The ECI and ECEF coordinate systems utilized in this paper are based on the J2000 and WGS84 coordinate systems, respectively.The methods for coordinate transformation are detailed in reference [26].
Because the ground station is equipped with an atomic clock and most clocks of LEO satellites are stable, the rate of change of relative frequency offset can be neglected.Thus, the dynamic equation of relative frequency offsets is modeled as By combining the position-velocity coordinate transformation model between the ECEF and the ECI frame, Equations ( 1) to (3), a measurement equation can be established, which represents the relation between multiple Doppler measurements Y = [Y(t 1 ), . . . ,Y(t k ), . . ., Y(t K )] T and the state vector of the orbit determination system X OD .Here, K denotes the total number of epochs.
It is worth noting that the state vector X OD depends on the orbit determination stage.During the coarse orbit determination process, where an improved search method is used, the accuracy requirements for the coarse orbit are low.Modeling small values of relative frequency offset as unknowns would increase the dimension of the search space, leading to a significantly low computational efficiency.Therefore, the state variables used for searching the coarse orbit consist only of the instantaneous satellite position and velocity.
When employing the LS estimation for the initial and precise orbit determination stages with higher accuracy requirements, incorporating the relative frequency offset not only improves the model's precision but also minimally affects computational efficiency.Consequently, the state vector consists of the instantaneous satellite position and velocity, and frequency offset.
It is worth noting that in the initial stage of LS, the orbital part is provided by IOD, and the frequency offset part is obtained from experience, which is 0 Hz in this study.

Three-Step Initial Orbit Determination Algorithm
To address the challenge of obtaining a reliable initial orbit for non-cooperative low Earth orbit satellites using Doppler observations, a three-step initial orbit determination algorithm is proposed.The algorithm consists of three steps: (1) coarse orbit estimation using an improved search method based on single-pass Doppler measurements, (2) data association to obtain multi-pass Doppler observations, and (3) initial orbit estimation based on coarse orbit and multi-pass Doppler through LS.

Coarse Orbit Estimation
The feasibility of orbit determination based on single-station single-pass Doppler measurements has been proven [21,22].This involves two techniques: (1) extracting a complete pass of Doppler observations from multiple-pass Doppler observations, and (2) orbit determination based on the extracted single pass observations.As the former is relatively straightforward, this study focuses on the latter.
Given that the publicly available document submitted by Motorola to the Federal Communications Commission (FCC) contains three design orbit parameters for the Iridium constellation [27], these three obtainable design orbit parameters X d (the semi-major axis a d , eccentricity e d , inclination i d ) can be used to reduce the dimensionality of optimization variables.
After dimensionality reduction based on known orbit parameters, only three parameters remain unknown: the right ascension of the ascending node Ω, the argument of perigee ω, and the true anomaly ν.Due to the high nonlinearity of the Doppler observations and the aforementioned unknown parameters, directly estimating these three parameters using single-station single-pass Doppler observations is challenging [22].Both the grid search method and particle swarm optimization (PSO) are capable of solving high dimensional nonlinear equations.As the dimensionality of optimization problems increases, the number of local minima in multimodal functions exhibits exponential growth [28].PSO is prone to falling into local optima when its parameters are not appropriate, and its performance is sensitive to control parameters [29], resulting in unstable performance.Grid search determines the optimal solution by searching all possible solutions, and although its computational burden is significant, its performance is not affected by the number of local optima.By combining the advantages of grid search and PSO, an improved search method is proposed, with the following implementation procedure.
First, the right ascension of the ascending node is divided into multiple grid points, expressed as Ω 1 g = [Ω 1 , . . ., Ω n , . . ., Ω N ] T , where the superscript of Ω 1 g represents the number of the search iteration, and the subscript of Ω n represents the order of elements.Second, the optimal variables, including the argument of perigee and the true anomaly for each grid point, are achieved through PSO.Denoting the vector formed by the difference between the actual observations and the theoretical observations calculated based on the estimated orbit as the residual, we select the modulus of the residual vector as the fitness function.In other words, by applying the PSO at every grid point Ω k (1 ≤ k ≤ N), the minimum fitness function F 1 g = [F 1 , . . ., F n , . . ., F N ] T and corresponding optimal variables X 1 g = [X o,1 , . . ., X o,n , . . ., X o,N ] T are obtained for Ω 1 g .Subsequently, by comparing the elements of F 1 g , the range in which the ascending node right ascension of the coarse orbit with the minimum residual modulus is located is determined.Then, the search grid Ω i g is redefined iteratively based on the narrowed range by repeating the aforementioned steps, where the superscript i denotes the iteration number.

Data Association
The purpose of data association here is to determine the most likely Doppler observations that originated from the same target from the acquired Doppler extracted by SOPs.Given the relatively small errors in extrapolating the coarse orbit over a short period, the calculated elevation angle, Doppler, and Doppler rate by the coarse orbit can be approximated as the true values of the satellite for data association.
The specific procedure can be summarized as follows.First, the conversion from the six orbital elements to the ECI coordinate system is performed to obtain the position and velocity (X s I (t I 0 )) c , where the subscript c denotes the coarse orbit.Then, by incorporating the positions of the station, the relation between the elevation angles at the observation epochs and instantaneous orbit h e (•) can be obtained.Additionally, assuming the relative frequency offset is zero, the relation between the theoretical Doppler and instantaneous orbit h(•) is established using the coordinate transformation between ECEF and ECI, Equations (1) to (3).Moreover, the relation between the Doppler rate and instantaneous orbit ḣ(•) is modeled by calculating the quotient of the Doppler difference and the time difference between adjacent epochs.Then, the theoretical elevation, Doppler, and Doppler rate are obtained by substituting the coarse orbit (X s I (t I 0 )) c into the aforementioned three models.Finally, the associated data are obtained, which are the result of satisfying the following three constraints among all observations within a short time period: (1) the calculated elevation angle at the observation epoch is between 0 degrees and 90 degrees, (2) the Doppler residual is less than a specified threshold, and (3) the Doppler rate residual is less than a preset threshold.

Initial Orbit Estimation
After obtaining the multi-arc observations Y I r−1 , the initial orbit can be estimated by refining the coarse orbit through LS based on Newton iteration.For the rth iteration, the estimated state value ( X s I (t I 0 )) r is calculated using the following equation: In this expression,( X s I (t I 0 )) r−1 represents the estimation obtained from the previous iteration, and ( X s I (t I 0 )) 0 refers to the vector composed of the position and velocity coordinates transformed from the coarse orbit (X s I (t I 0 )) c , with the relative frequency offset assigned a value of zero.Here, H I r−1 represent the observation matrix, and y I r−1 represents the residual, which is the difference between the observed Doppler Y I r−1 and the theoretical Doppler h(( X s I (t I 0 )) r−1 ).The least squares algorithm is considered to have converged when the root mean square (RMS) of the positional correction quantities in three dimensions is smaller than a user-specified threshold ε R .

Fast Orbit Calculation
Orbit calculation is employed to obtain the orbital information at other desired times based on the satellite's dynamic conditions using the known instantaneous orbit.The computational time required for orbit calculation significantly impacts the real-time performance of orbit determination.This section presents a description of the orbit calculation problem as well as the algorithms for fast single-orbit calculation (FSOC) and fast multiorbit calculation (FMOC).

Problem Description
The objective of orbit calculation is to acquire the orbital information of a satellite at specific past or future single or multiple time instances, utilizing the position and velocity data at a given moment.Let X s ECI (t 0 ) represent the initial orbit at time t 0 , and let the time sequence of interest corresponding to the computed orbit be arranged in ascending order and denoted as t I = [t 1 , . . ., t k , . . ., t K ] T , where K is the number of elements in the time sequence of interest.The computed orbit sequence is denoted as (X s ECI ) C = [(X s ECI (t 1 )) T , . . ., (X s ECI (t k )) T , . . ., (X s ECI (t K )) T ] T , where the subscript C refers to the computed orbit.The mathematical model for orbit calculation is as follows: where X s ECI (t 0 ) and X s ECI (t k ) represent the position and velocity of the satellite in the ECI coordinate system at time t 0 and t k , respectively.The derivative of the satellite's position and velocity with respect to time is denoted as Ẋs ECI .The dynamic differential equation describing the relation between Ẋs ECI and X s ECI is as follows: In the given context, the relation between the derivative of the orbital state Ẋs ECI and the orbit X s ECI is represented by g RV (•).The satellite's acceleration in the ECI coordinate system is denoted by a s ECI , and the relation between the acceleration a s ECI and the state X s ECI is indicated by g a (•).
The primary perturbations affecting LEO satellites comprise Earth's non-spherical gravity and atmospheric drag.Regarding non-spherical gravity perturbation, with Earth's oblateness being the main factor causing perturbations, this perturbing force only considers the gravitational potential function with fourth-order zonal harmonics [30].As for atmospheric drag perturbation, it is omitted considering the following factor.For LEO satellites positioned at altitudes of approximately 800 km, the maximum error in a 24-h orbit prediction due to neglecting atmospheric drag is approximately 320 m [31].Notably, the sole available real-time orbit data for Iridium satellites are the TLE, characterized by an accuracy on the order of kilometers.Given the objective of achieving kilometer-level orbital accuracy, errors within the range of hundreds of meters in orbit prediction are acceptable.Consequently, in the perturbation force modeling presented in this paper, atmospheric drag is disregarded, and only the fourth-order gravitational potential with zonal harmonic perturbations is considered.

Fast Single-Orbit Calculation Algorithm
Due to the numerical stability exhibited by the Runge-Kutta method, which enables it to handle larger time steps without abrupt divergence, a Runge-Kutta integrator with a relatively large step is chosen as the numerical integration method.This method is employed iteratively to integrate Equation ( 6) and facilitate the computation of the position and velocity of the satellite at the desired time, denoted as X s ECI (t 1 ), based on the known instantaneous orbit X s ECI (t 0 ).The schematic diagram of the fast single-orbit calculation algorithm is depicted in Figure 3.The squares represent the starting time t 0 and ending time t 1 , and the circles denote the intermediate time involved in the process.The integration process is executed to forecast future orbits when t 1 > t 0 , as visually represented by the green color in the diagram.Conversely, when t 1 < t 0 , the integration is carried out to compute past orbits, indicated by the blue color.The algorithm can be summarized into three steps as follows.(1) The integration interval is divided based on the initial time t 0 , final time t 1 , and integration step size δ.The resulting division is denoted as t FSOC = [t 0 , t 0 + T FSOC , . . ., t 0 + iT FSOC , . . ., t 1 ] T , where  6) using the Runge-Kutta method.The integration process begins from the initial orbit X s ECI (t 0 ) and employs the integration limits t 0 + (i − 1)T FSOC and t 0 + iT FSOC for the ith iteration.(3) The integrated orbit obtained in the last iteration X s ECI (t 1 ) is provided, yielding the satellite's position and velocity at the specific future or past time of interest, denoted as t 1 .

Fast Multi-Orbit Calculation Algorithm
The multi-orbit calculation is commonly used for orbit determination in a batch form, computing the satellite's position and velocity at multiple observation epochs.In most cases, the observation epochs are closely spaced, leading to a notable correlation among orbits at different epochs.Directly applying the fast single-orbit calculation algorithm to each observation epoch would disregard the strong correlation among them, leading to a significant amount of redundant computations.Due to the low eccentricity and relatively smooth orbital variations of LEO satellites, a limited set of satellite positions and velocities at specific moments can effectively capture comprehensive orbital information.Therefore, in order to reduce significant redundant computations, a strategy has been introduced.The strategy consists of two steps: (1) obtaining a small set of satellite positions and velocities as intermediate variables X N I through FSOC, and (2) utilizing these intermediate variables to fit the orbits (X s ECI ) C at the desired time of interest t I .The diagram of FMOC is illustrated in Figure 4. Due to the uncertain relationship between the time of interest t and the initial epoch t 0 , three scenarios are considered: (1) all elements of t I are smaller than t 0 , as depicted in the red region, (2) t 0 falls within the interval defined by the maximum and minimum values of t I , as illustrated in the blue region, and (3) all elements of t I are larger than t 0 , as shown in the gray region.Within the elements comprising tI, the one exhibiting the smallest absolute difference in relation to t 0 is identified as t Int1 .This particular element's corresponding index is denoted as Idx 1 and is visually highlighted by the gold dashed line.that exhibits the smallest minimum difference in relation to t 0 .(3) The FSOC is executed with X s ECI (t 0 ) as the initial value to compute the orbit X Idx1 at the specific time t Int1 .(4) The FSOC is performed with X Idx1 as the initial value to obtain the integrated orbit X N I as an intermediate variable.(5) The intermediate variable X N I is employed to fit the sequence of orbital parameters at time t I and, consequently, derive the calculated orbit (X s ECI ) C .

Experimental Results
Because Iridium is a widely studied non-cooperative LEO constellation, it is used as an example to demonstrate the effectiveness of the proposed method.In this section, the efficacy of the proposed method was assessed through practical opportunistic localization experiments, accompanied by a comparative analysis of performance across various orbit calculation methods.
It should be noted that the high-precision orbit of Iridium is not public, and, currently, the only available real-time orbit is the TLE with an accuracy of kilometers.TLE can be used to roughly evaluate the orbit determination performance, but cannot accurately evaluate the actual orbit determination error.In order to verify the effectiveness of opportunistic positioning based on the proposed method, a stationary receiver is used as the positioning terminal, and positioning is carried out based on the orbit determination results and observations extracted by the next pass of the satellites whose orbits are estimated.

Experimental Environment and Basic Settings
The experiment on non-cooperative orbit determination of LEO satellites was conducted using Iridium satellites as a representative case due to frequent utilization as opportunistic signal emitters.
As for the hardware of the orbit determination system, it was situated at the New Main Building of Beihang University, with approximate geographic coordinates of 116.3 degrees longitude and 40.0 degrees latitude, and an elevation of 96.6 m, as depicted in Figure 5.The satellite antenna is installed on the rooftop to capture SOPs, and the signal processing center is located within the building and with dedicated computers developed by the Communication Navigation and Testing Laboratory (CNT Lab) for long-term continuous orbit determination of Iridium NEXT satellites with NORAD identifications 42807 and 43927.When it comes to the software, the orbit determination system initially extracts the Doppler observations, followed by orbit determination using these observations.Subsequently, the effectiveness of the obtained orbit is validated through the process of two position experiments.The epochs for observation, orbit determination, and positioning are presented in Table 1.
Regarding the basic settings, during the initial orbit determination stage, the particle swarm optimization algorithm [32] employs a population size of 30, and the maximum number of iterations is set to 200.The learning factors, denoted as c 1 and c 2 , are both set to 1.5, and the inertia weight is fixed at 0.5.The thresholds for Doppler residuals and Doppler rate of change residuals related to data association are empirically set to 1500 Hz and 300 Hz/s, respectively.The three-dimensional root mean square threshold ε R used for position correction is set to 1 m based on experience.
In terms of orbit calculation, the proposed algorithm is adopted.The numerical integration step size δ is 30 s, and extrapolation factors a 1 and a 2 are both set to 0.5; the fitting method employed is cubic spline interpolation.Obtaining the signal structure feature is a prerequisite for extracting Doppler measurements.As for the Iridium constellation, its subscriber link frequency is divided into duplex and simplex operation channels, whose frequencies are 1616.0-1626.0MHz and 1626.0-1626.5 MHz, respectively.By controlling the transmitted spot beams, Iridium NEXT realizes a channel for users.In addition, the signal structure of Iridium NEXT downlink transmission in the seventh channel mainly consists of pilot, BPSK, and QPSK modulation, which can be received every 4.32 s.

Location of ground station
Since there are unmodulated pilot data in the seventh channel with a carrier frequency of 1626.270833MHz, the received signal frequency can be extracted by the fast Fourier transform (FFT) algorithm on the downlink pilot data.By calculating the difference between the received signal frequency and the transmitted signal frequency, the Doppler frequency offset can be obtained.

Initial Orbit Determination
The orbit derived from TLE serves as a reference orbit for assessing performance.The initial orbit determination process for satellites 42807 and 43927 was conducted employing the three-step method.Figure 6 illustrates the outcomes for each step of satellite 42807.Satellite 43927 yielded identical outcomes in the initial step, followed by the acquisition of the initial orbit through the second and the third steps.
From Figure 6a, it is apparent that both satellites demonstrate two distinct points of extremity across the entire search range of 0 degrees to 360 degrees.Among these points, the extremity between 120 degrees and 160 degrees possesses a smaller value, indicating that this particular range corresponds to the minimum.Consequently, a coarse orbit determination is carried out by conducting a localized search within the range of 120 degrees to 160 degrees.According to Figure 6b, the optimal ascending node right ascension is 136 degrees, and the argument of perigee and the true anomaly are obtained using the PSO algorithm based on the coarse orbit estimation in the three-step method under the ascending node right ascension parameter of 136 degrees.
By incorporating these specific right ascension values with the optimal true argument of periapsis and true anomaly, coarse orbits can be derived.Utilizing these coarse orbits, Doppler observations are obtained for an additional arc within a 1.5-h interval, resulting in the generation of multi-arc Doppler observations for the same satellite.These observations are illustrated in Figure 6c.The LS method is subsequently applied iteratively using the multi-arc Doppler observations and the initial values of the coarse orbit.The root mean square error (RMSE) of position and velocity for each iteration is presented in Figure 6d, demonstrating that, after three iterations, a reasonably precise initial orbit can be estimated.

Orbit Calculation
In order to assess the numerical stability of various single-orbit calculation methods, a comparative analysis of computed orbit RMSE is conducted using different integration steps between the SFOC method and the Adams-Bashforth method.The results of this comparison are illustrated in Figure 7, where the Kepler equation analytical model is employed as the reference standard.As depicted in Figure 7a, it can be observed that with a one-day extrapolation, the position and velocity RMSE of the SFOC method increase as the step size increases; however, the rate of divergence is relatively gradual.For instance, when using a 30-s integration step, the RMSE for position is on the order of 10 m, while the RMSE for velocity is on the order of 10 m and 0.01 m/s.In contrast, the RMSE for position and velocity in the Adams-Bashforth method demonstrates rapid divergence as the step size increases, as illustrated in Figure 7b.The errors exceed or attain magnitudes of 10 9 m and 10 5 m/s, respectively, after a 20-s step.Therefore, it can be concluded that the FSOC method has more numerical stability compared to the Adams-Bashforth method.
To compare the disparities in results and efficiency between multi-orbit calculation methods, the orbit calculation involved in the first iteration of POD for Iridium 42807 is taken as an illustrative example.The deviations in the computed orbit positions and velocities are compared for the specific step sizes using the FMOC and the Adams-Bashforth methods.These comparisons are presented in Figure 8.In this context, the step size for the FMOC is set to 30 s, which aligns with the step size utilized in the orbit determination system.On the other hand, the step size for the Adams-Bashforth method is determined empirically as the maximum value that guarantees the convergence of the integrator, which equals 16.7 s.Based on the results from Figure 8, it can be deduced that, for an extrapolation period of approximately one day, employing the specified step sizes, the maximum discrepancy in computed orbit positions and velocities between the two methods is approximately 45 m and 0.045 m/s, respectively.Moreover, according to the recorded data, the computational time for the FMOC method amounts to 51 s, whereas for the Adams-Bashforth method, it is 23,279 s.Consequently, in comparison to the Adams-Bashforth method, the FMOC method demonstrates comparable accuracy while reducing the computational time by 99.8%.It is worth noting that this experiment was conducted using MATLAB.If implemented in C language, it will further reduce the runtime.

Precise Orbit Determination and Orbit Prediction
To obtain accurate orbit estimations, multiple-pass Doppler observations are acquired through data association.The Doppler observations utilized for the precise determination of the orbits of Iridium satellites 42807 and 43927 are depicted in Figure 9. Based on statistical analysis, the time intervals between consecutive visible arcs are distributed within two intervals: 1.4-1.6 h and 9.4-11.4h.
Figure 10 illustrates the RMSE in the precise and prediction orbits of satellites 42807 and 43927, which are obtained through precise orbit determination based on Doppler observations and orbit prediction utilizing the latest precise orbit information.It can be inferred from the figure that the proposed framework enables the continuous and long-term determination of satellite positions and velocities.The determined orbits exhibit position and velocity errors within 1 km and 1 m/s, respectively, when compared to the orbital positions and velocities derived from TLE.As TLE provides orbital accuracy at the kilomter and meters per second levels, the proposed method achieves the same level of accuracy.

Opportunistic Positioning
Two opportunistic localization experiments were conducted using TLE and precise orbit by the proposed method, with the corresponding duration listed in Table 1.The sky map for each experiment is shown in Figure 11.The positioning RMSE using TLE and the precise orbit are presented in Table 2.It can be observed from Table 2 that positioning with the precision orbit yields higher accuracy in positioning compared to the utilization of TLE-derived orbit.

Conclusions
This paper introduces a novel method for orbit determination for non-cooperative LEO satellites using a single ground station, with application to opportunistic positioning.The proposed method offers several advantages, including high reliability, low hardware cost, and the elimination of the need for a unified space-time reference among ground stations and information transmission.To enable real-time opportunistic positioning based on reliable satellite orbits, a new orbit determination framework is proposed.In this framework, ground stations utilize Doppler measurements from SOPs as the observation data.The initial orbit or the previously updated precise orbit is iteratively refined by incorporating newly acquired measurements and the accumulated Doppler observations.The most recent precise orbit serves as the initial value for orbit prediction, enabling the determination of positions and velocities at the desired time of interest.To obtain a reliable initial orbit, a three-step initial orbit determination algorithm is proposed.(1) An improved search method is employed to estimate a coarse orbit based on single-pass Doppler measurements.
(2) Data association is performed based on the coarse orbit to obtain multi-pass observations.(3) The initial orbit is estimated by refining the coarse orbit with the associated data.Furthermore, two fast orbit calculation algorithms are proposed.First, the numerical stability of the Runge-Kutta method is leveraged to reduce the number of integration times.Second, the strong temporal correlation in LEO orbits is exploited to minimize redundant computations.As a result, the proposed algorithms achieve high computational efficiency in orbit calculation.The effectiveness was verified through positioning experusing real observations.The results demonstrated comparable precision between the autonomously determined orbits and those obtained using TLE.This method provides long-term reliable orbit determination using a single station, achieving accuracy on par with TLE.These findings have significant implications for enhancing the reliability of opportunistic positioning.

Figure 1 .
Figure 1.Diagram of the non-cooperative LEO satellite orbit determination using a single station for opportunistic positioning framework.

Figure 2 .Figure 2 .
Figure 2.An illustration of the important concepts.

Figure 3 .
Figure 3. Diagram of the fast single-orbit calculation algorithm.
and ⌈•⌉ represents the operation of ceiling rounding.(2) The algorithm performs the calculation of the orbit for each epoch t FSOC by iteratively integrating Equation (

Figure 4 .
Figure 4. Diagram of the fast multi-orbit calculation algorithm.FMOC consists of five steps: (1) The time t N I of the numerical integrated orbit X N I is determined based on the desired time of interest t I , step size δ, and extrapolation factors a 1 and a 2 .The time interval t FMOC is defined as t FMOC = [t I (1) − a 1 δ, t I (1) − a 1 δ + T FMOC , . . . ,t I (K) + a 2 δ + T FMOC ] T , where T FMOC = T N FMOC .The value of T is determined as T = tI(K) − tI(1) + (a 1 + a 2 )δ, where t I (K) and t I (1) represent the last and first elements of t I , respectively.The extrapolation factors a 1 and a 2 are positive numbers.The number of steps for the FMOC is given by N FMOC = ⌈T⌉ δ .(2) We identify the element within t N I

Figure 5 .
Figure 5. Orbit determination system for experiment.

Figure 6 .
Figure 6.Results of the initial orbit determination for each step of Iridium 42807.(a) Shows the curve of the optimal fitness function varying with the RAAN across the entire search range.(b) Shows the curve of the optimal fitness function varying with the RAAN within the narrowed search range.(c) Shows the multiple-arc Doppler observations obtained through data association.(d) Shows the RMSE of position and velocity for each iteration based on LS.

Figure 7 .Figure 8 .
Figure 7.The RMSE of the calculated orbit's position and velocity using different methods as a function of extrapolation time.(a) Shows the RMSE of the calculated orbit's position and velocity by the FSOC as a function of extrapolation time.(b) Shows the RMSE of the calculated orbit's position and velocity by Adams-Bashforth as a function of extrapolation time.

Figure 9 .
Figure 9. Doppler observations used for the precise determination of the orbits of Iridium satellites 42807 and 43927.(a) Shows Doppler observations used for the precise orbit determination of Iridium 42807.(b) Shows Doppler observations used for the precise orbit determination of Iridium 43927.

Figure 10 .
Figure 10.The RMSE in the precise and prediction orbits of satellites 42807 and 43927.(a) Shows the RMSE in the precise and prediction orbits of satellites 42807.(b) Shows the RMSE in the precise and prediction orbits of satellites 43927.

Figure 11 .
Figure 11.Sky map for each experiment.(a) Shows the sky map for the first experiment.(b) Shows the sky map for the second experiment.

Table 1 .
The epochs for observation extraction and positioning experiments.

Table 2 .
The positioning RMSE using TLE and the precise orbit.