Research on Initial Orbit Determination of Near-Circular Orbit Space Objects Based on the Admissible Region Method

The orbital energy ε of an Earth-orbiting object must satisfy the two-body motion condition. From the geocentric inertial position vector r and the velocity vector ṙ, and using the topocentric spherical coordinate representation, the energy equation is derived:

After simplification, an explicit relationship involving ρ and $\dot{\rho }$ is obtained (DeMars et al. 2012):

Where w₁ is the projection coefficient of the observer velocity onto the line-of-sight direction, and F(ρ) is a function of ρ (Zhao 2020). Given the energy ε, one can solve for:

Using the elliptical orbit constraint, solution curve boundaries and their interior region D₁: ε < 0 can be generated.

Using the relation between semi-major axis a and energy, ε = –μ / (2a), combined with the energy equation:

Substituting into the energy constraint equation yields the solution boundary curve corresponding to the semi-major axis upper limit, with its interior region satisfying D₂: a ≤ a_max.

From the angular momentum h = r × ṙ and the eccentricity formula:

A fourth-order polynomial equation in ρ and $\dot{\rho }$ is derived:

The coefficients a_i are functions of ρ. Given an eccentricity upper limit e_max, solving this equation yields the boundary, with its interior satisfying D₃: e ≤ e_max.

The intersection of the above constraint regions, D₁ ∩ D₂ ∩ D₃ (i.e., the energy, semi-major axis, and eccentricity constraint regions), forms a compact two-dimensional solution domain. As shown in Fig. 1, the blue curve represents D₁; the green curve represents D₂; the red curve represents D₃; and the asterisk marks the objects true solution. This domain contains the objects true orbital parameters (ρ, $\dot{\rho }$), providing a prior solution space for subsequent probabilistic modeling (e.g., uniform distribution assumption) or filtering estimation. It can be seen from Fig. 1 that the feasible region corresponding to the eccentricity constraint is divided into two areas: one contains the true solution, while the other constitutes an erroneous solution that needs to be eliminated in automated processing.

Fig. 1. Example of an admissible region construction.

Download Original Figure

For space-based observation of near-circular orbit objects, the AR algorithm can produce erroneous solutions, which severely impacts the subsequent determination of orbital parameters. The cause for the appearance of these erroneous solutions is briefly explained below. For a circular orbit (eccentricity zero), the eccentricity constraint simplifies significantly. Eq. (9) can be reduced to:

Combining with the semi-major axis and energy relationship, it further simplifies to:

Furthermore, due to the circular orbit, h = a²ω, where ω is the mean angular velocity, yielding:

Eq. (13) establishes the relationship between semi-major axis and mean motion under the circular orbit assumption. This forms the basis for the single-parameter initial orbit determination method discussed in the following section.

First, assume an initial semi-major axis â. The theoretical mean motion n can be obtained via Keplers law:

Next, the computed value of mean motion n₂ is given by â and r_i:

From Eqs. (4) and (5), the objective function can be constructed as:

when f(â) = 0, the function converges, and the resulting â is the solution. In the above equation, r_i is unknown and needs to be obtained using ρ_i. Utilizing the geometric relationship between the telescope and the object, ρ_i can be expressed in terms of â and other known quantities:

where z is the zenith distance, whose geometric relationship is illustrated in Fig. 2.

Fig. 2. Geometry of the single-parameter method.

Download Original Figure

Although Eq. (13) suggests a one-to-one correspondence between a and ω under ideal conditions, the practical solution of Eq. (16) can yield multiple roots. This occurs because $n_1\left(\widehat{a}\right)=\sqrt{\mu /a^3}$ is strictly decreasing in â, while n₂(â) is generally non-monotonic. The non-monotonicity of n₂(â) arises from the dependence of r₁ and r_end on â via the slant range ρ (Eqs. (17)–(19)), and the fact that the dot product r₁ · r_end involves the cosine of the angle between two geocentric position vectors. As â varies, this angle can change non-linearly, causing n₂(â) to increase or decrease in a complex manner. Consequently, the decreasing curve of n₁(â) may intersect the non-monotonic curve of n₂(â) at two distinct points, yielding two different â values that satisfy f(â) = 0. One corresponds to the true orbit; the other is a mathematically valid but physically incorrect spurious solution—referred to as an erroneous solution.

In essence, the apparent angular rate derived from angle-only measurements depends not only on â but also on the instantaneous line-of-sight geometry. Different combinations of â and geometry can produce the same observed angular rate, leading to multiple mathematical solutions that satisfy all physical constraints. This is the core reason for the multi-solution problem.

In the space-based scenario, the platform itself moves along a Keplerian orbit, introducing additional geometric variations that make the dependence of n₂(â) on â more intricate. This increases the chance of multiple intersections, as illustrated by the examples in Section 2.3.3.

Existing space-based long-range surveillance platforms primarily operate in near-circular orbits, and small-eccentricity orbits offer excellent performance in terms of the number of observed objects and observation arc length for low-Earth orbit multi-object surveillance (Zhao 2020). The space-based observation discussed in this paper is also based on long-range near-circular orbit observation platforms. The key difference between space-based and ground-based observation lies in the different motion laws followed by the observation platform itself. The ground-based platform is subject to Earth’s rotation, whereas the space-based platform, like the observed object, follows Keplers laws. This distinction makes the single-parameter algorithm more prone to generating erroneous solutions due to geometric configurations when applied to space-based observation data.

For the objective function Eq. (17), when the semi-major axis value equals the observation platforms own a_obs, the function also equals zero. This solution is considered a trivial solution and can be filtered out by applying a slant range threshold. However, beyond this, consider the following geometric configuration: the object orbit and the observation platform orbit are coplanar, with velocity directions aligned. In this case, within the reasonable range of the independent variable, both $\widehat{\text{n}}$ and n are monotonically decreasing functions, and two solutions exist that satisfy the objective function equaling zero (the true solution and the trivial solution). If the orbital elements of the two satellites remain unchanged, but their velocity directions are opposite, within the same value range, $\widehat{\text{n}}$ remains unchanged, but n exhibits a minimum point, and the extremum is zero. The geometric configuration of the observation line of sight changes continuously, and the corresponding objective function also changes. To intuitively understand this characteristic, Fig. 3 shows a hypothetical space-based observation scenario, and Fig. 4 shows the corresponding function variation.

Fig. 3. Hypothetical space-based observation scenario.

Download Original Figure

Fig. 4. Variation of the objective function under different orbital plane angles.

Download Original Figure

As can be seen from Fig. 4, during the process of the observation platforms orbital plane changing from 0° to 180°, geometric configurations exist that cause the objective function to produce erroneous solutions. These solutions conform to physical laws and are difficult to distinguish from the true solution. Fig. 5 demonstrates the corresponding effect of combining the AR and the single-parameter algorithm. The x-axis in the figure is a shared coordinate axis. It can be observed that both algorithms yield two solutions, one true and one erroneous. However, the solutions from the two algorithms are not entirely consistent. This discrepancy arises because the meanings of the solutions differ: one represents the slant range in the AR algorithm, and the other represents the semi-major axis in the single-parameter algorithm. The figure shows that erroneous solutions corresponding to the single-parameter algorithm exist within the AR. Such erroneous solutions are not limited to the example geometric configuration. Therefore, a method to eliminate erroneous solutions is needed to enable the AR algorithm to be better applied to space-based angle-only initial orbit determination.

Fig. 5. Correspondence between admissible region (AR) and single-parameter algorithm solutions.

Download Original Figure

In the initial orbit determination process, to effectively identify and eliminate erroneous solutions, this paper constructs multiple sets of observation data with varying geometric configurations by changing the observation arc length, start/end times, or a combination of both. Based on the key characteristic that the true solution is consistent across different observation geometries, whereas erroneous solutions are sensitive to such geometric changes, to effectively distinguish between them, geometric configurations with significant differences need to be constructed to amplify the distinction in solution set responses. To this end, this paper systematically selected six observation combinations aimed at maximizing geometric diversity. The selection primarily considered two key dimensions: arc length (shortest, medium, longest) and the relative position of the combinations central time within the observation time window (earliest, middle, latest). Consequently, the following six pairs were chosen:

Based on the method in Section 2.3.2, the orbital semi-major axis can be estimated using any two observation points. Therefore, six different data combinations can be extracted from one observation arc segment, solved separately to form a set of semi-major axis solutions {a_1,1, a_2,1, ..., a_6,1}. If multiple solution groups exist, all solution sets can be represented as{a_1,1, a_2,1, ..., a_6,1}, ..., {a_1,m, a_2,m, ..., a_6,m}. It should be noted that when the number of data points is insufficient to select point pairs according to the above scheme, all possible combination pairs can be used as the selection scheme.

The true solution typically does not fluctuate significantly with changes in geometric configuration, whereas erroneous solutions are sensitive to changes in observation conditions. To characterize this difference, the variance within each solution set can be calculated to evaluate its dispersion (relative variance, as shown in Eq. (20), is used to avoid the influence of the expected value, where s̅ represents the mean of the solution set). Solution sets with smaller variance have higher stability and are more likely to correspond to the true orbit; solution sets with larger variance are judged as erroneous and excluded.

Furthermore, the existence of trivial solutions cannot be ignored. Especially when the zenith angle is less than 90°, their numerical values are less affected by geometric configuration, making it difficult to effectively eliminate them using the aforementioned variance analysis. Therefore, this paper incorporates a slant range ρ judgment threshold: if the solved slant range is below a positive number close to zero, it is considered a trivial solution, and that solution set is directly excluded.

To validate the rationality of the proposed six-group observation point selection scheme, the following comparative experiment was designed. The experiment aimed to systematically compare this scheme with two control schemes—the “scheme using only two points near the start and end times” (2-group scheme) and the “scheme using all possible point pairs” (full combination scheme)—in terms of performance when generating single-parameter solution sets.

The experiment utilized simulation data from three representative orbit types (LEO, MEO, GEO), with example orbital parameters shown in Table 1. For each orbit type, 100 repeated experiments were conducted. In each repeated experiment, an independent set of Gaussian random noise with a mean of 0 and a standard deviation of 3 arcseconds was generated and added to the simulated angle observation sequence to simulate the randomness of actual observation errors. Subsequently, for the same noise-added observation arc segment, the above three point-pair selection strategies were applied in sequence to solve the orbital semi-major axis based on the single-parameter method, forming the corresponding solution set S. The relative variance RV (measuring stability), mean, and average single computation time for each strategys solution set were recorded.

The experimental results are shown in Table 2. In terms of solution stability, for all orbit types, the relative variance of the solution sets from the 6-group scheme was significantly lower than that of the 2-group scheme, indicating smaller dispersion and stronger stability. Simultaneously, its stability was comparable to that of the full combination scheme. In terms of computational efficiency, the time consumption of the 6-group scheme, although slightly higher than the 2-group scheme, was far lower than the full combination scheme. For example, in the LEO scenario, the six-group scheme took 0.29 seconds, while the full combination scheme took 56.40 seconds, representing an improvement in computational efficiency of over two orders of magnitude.

Overall, the six-group scheme achieves an optimal balance between stability and efficiency, significantly improving solution set stability (compared to the 2-group scheme) while its computational cost is far lower than the full combination scheme. Therefore, adopting this scheme as the core strategy for constructing multi-geometry-configuration solution sets is reasonable and efficient.

Random noise in optical observation data is a key factor affecting the accuracy and solution set stability of initial orbit determination algorithms. To suppress noise, eliminate gross errors, and provide a smooth and reliable observation sequence for subsequent calculations, this paper employs a cubic polynomial fitting method to smooth the right ascension and declination observation sequences separately. Given the short observation arc, the objects celestial coordinates change smoothly. Using the time offset relative to the first observation epoch as the independent variable, the following fitting model is constructed:

The polynomial coefficients are obtained through least squares estimation, and then the smoothed values ${\widehat{\text{α}}}_{\text{i}}$, ${\widehat{\text{δ}}}_{\text{i}}$ at each epoch are calculated. To eliminate the influence of gross errors, an iterative processing strategy is adopted: after calculating the residuals at each observation point, outliers are identified and eliminated based on the Pauta criterion (3σ criterion). Polynomial fitting is then repeated using the remaining points, iterating until no gross errors remain, ultimately obtaining a clean, smooth observation sequence.

To systematically evaluate the effectiveness of this preprocessing method, a controlled experiment was designed. The experiment used the same sensor and object orbital parameters (LEO, MEO, GEO) as in Section 2.4.1. Under noise levels of 0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8 arcseconds, the following three processing methods were compared: 1) No noise reduction (N); 2) The polynomial fitting method of this paper polynomial fitting (PF); 3) Iterative Reweighted Least Squares (IRLS). Each experimental group underwent 100 Monte Carlo repetitions, using the relative variance (RV) of the solution set as the evaluation metric, and distinguishing between correct solutions (small RV and orbital parameter error < 200 km) and erroneous solutions.

The experimental results are shown in Figs. 6–8 (corresponding to LEO, MEO, GEO, respectively). The horizontal axis represents noise level, and the vertical axis represents relative variance (exponential coordinates for 0–1 arcsecond, linear coordinates for 1–8 arcseconds). Blue, green, and red lines represent no noise reduction, polynomial fitting, and IRLS methods, respectively. Circular markers indicate erroneous solutions, square markers indicate correct solutions, and star markers indicate conditions under which the relative variance of the correct solution is smaller than that of the erroneous solution, meaning the algorithm can effectively eliminate erroneous solutions.

Fig. 6. Effect of denoising methods on solution set relative variance for an LEO object.

Download Original Figure

Fig. 7. Effect of denoising methods on solution set relative variance for an MEO object.

Download Original Figure

Fig. 8. Effect of denoising methods on solution set relative variance for a GEO object.

Download Original Figure

Analysis of the charts reveals that: across all orbit types and noise levels, both polynomial fitting and IRLS methods significantly reduce the relative variance of both correct and erroneous solutions, with the relative variance of correct solutions generally lower than that of erroneous solutions. Especially when noise ≥ 3″, the relative variance of correct solutions under the polynomial fitting method remains stable, while that of erroneous solutions increases significantly, creating a clear distinction. This indicates that the preprocessing method of this paper not only effectively smooths observation noise and eliminates gross errors but also enhances the discriminability of the solution sets, providing a reliable basis for subsequent erroneous solution elimination based on relative variance. Overall, the polynomial fitting method, while ensuring computational efficiency, achieves noise reduction and solution set separation effects similar to the IRLS method, and is therefore selected as the preprocessing method in this paper.

Addressing the erroneous solution problem in the AR algorithm for space-based angle-only initial orbit determination, this paper proposes an improved calculation method incorporating data preprocessing. The core idea of the algorithm is to distinguish between true and erroneous solutions based on the difference in stability of solution results under different geometric configurations after preprocessing the observation data. A pseudo-code description of the complete calculation flow of this method is shown in Algorithm (Fig. 9).

Fig. 9. Algorithm: RV-SPSS Method

Download Original Figure

This flow enhances the quality of input data through a front-end data preprocessing module. Subsequently, the multi-geometry-configuration solution strategy amplifies the difference in statistical stability between true and erroneous solutions. Finally, through dual-threshold criteria based on relative variance and slant range, effective elimination of erroneous and trivial solutions is achieved. This method systematically integrates observation data purification, multi-solution path exploration, and statistical decision-making, significantly improving the reliability and robustness of space-based near-circular orbit initial orbit determination.

This set of experiments aims to construct a high-precision, highly realistic simulated observation environment to verify the core capability of the RV-SPSS method to identify and eliminate erroneous solutions under complex, near-actual working conditions. The data generation follows the process below: First, 50 LEO objects with continuous ground-based photoelectric observation data between January and June 2023 are selected. Their historical TLEs are used as initial orbits, and corresponding ground-based angle measurement data are used for precise orbit determination and prediction, generating high-precision orbital ephemerides (position accuracy better than 10 m) within the reference period. This study sets the space-based observation platform in LEO orbit, as it can perform wide-area surveillance of medium and high-altitude objects and is more likely to obtain very-short-arc observation data for a large number of space objects. This configuration aligns with the survey mission mode of space-based space surveillance, and existing research indicates that LEO observation platforms have advantages in associating multi-object observation data (Fujimoto & Scheeres 2013). This scenario can effectively validate the erroneous solution suppression capability of the proposed algorithm in surveillance mode. Subsequently, the SWARM satellite, which has global observation data and laser ranging support, is selected as the simulated space-based platform to generate high-precision platform orbits (orbit determination root mean square (RMS) ~0.7 m) for the same period. Finally, the precise ephemerides of the object and platform are differenced to calculate the objects right ascension and declination observation sequence relative to the platform. Gaussian white noise with zero mean and a standard deviation of 2 arcseconds is added to simulate measurement errors. This value is chosen to be representative of the typical accuracy of space-based optical sensors, which is reported to be in the range of 0.8–1.5 arcseconds for modern star trackers (Agostinelli et al. 2024).

The orbit extrapolation and precise orbit determination processes employed high-fidelity dynamic models; specific parameter settings are shown in Table 3. Furthermore, the simulation set the arc length to 80 s. Fig. 10 shows the parameter information for the simulated arcs. It should be noted that although actual space-based optical observation is also subject to various constraints such as Earth shadow, solar stray light, and camera field of view, these factors primarily affect object visibility and data quality and have a minor impact on the observation geometry and erroneous solution generation mechanism central to this study; therefore, they were not included in the simulation.

Fig. 10. Distribution of orbital parameters in the simulated data. INC, inclination; ECC, eccentricity; SMA, semi-major axis.

Download Original Figure

To test the algorithms applicability and erroneous solution suppression capability across different orbital configurations, this set of experiments constructs a large-scale random simulation scenario. The observation platform is set in a typical near-circular Low Earth Orbit. Object orbits are randomly generated for 300 space objects based on the parameter ranges for LEO, MEO, and GEO objects in Table 4, ensuring that objects are not occulted by Earth relative to the platform and are within the preset field of view during the simulation period. Based on two-body dynamics considering J2 perturbation, numerical integration generates precise orbits for objects and platforms, from which noise-free angle observation sequences are calculated. Finally, the measurement noise is modeled as zero-mean Gaussian white noise with a standard deviation of 2 arcseconds. The observation arc length is set to 80 s. This scenario aims to systematically verify the broad effectiveness and stability of the erroneous solution exclusion mechanism of the RV-SPSS method when facing a large number of randomly generated orbital objects and different observation arcs.

Fig. 11 shows a scatter plot of the solution set distribution based on TLE-derived data. The horizontal axis represents the mean of the solution sets SMA, and the vertical axis represents the solution sets RV. The plot distinguishes three cases: “Correct Solution,” “Erroneous Solution,” and “Large-Error Solution.” “Correct Solution” refers to the solution with SMA error < 200 km and the smallest relative variance; “Erroneous Solution” refers to solutions that are not the one with the smallest relative variance; “Large-Error Solution” refers to the solution that, although having the smallest relative variance, is not the one with the smallest error.

Fig. 11. Scatter plot of solution set distribution from historical data simulation.

Download Original Figure

It can be seen from the figure that erroneous solutions and large-error solutions are more dispersed in the relative variance-semi-major axis plane, while correct solutions are relatively concentrated, verifying the insensitivity of the true solution to observation geometric configurations. Among a total of 269 observation arcs, 239 exhibited erroneous solutions, accounting for approximately 88.8%.

Table 5 compares the performance of different IOD methods on TLE-derived data, including the original AR-based single-parameter method, the proposed RV-SPSS method, as well as the Gooding and improved Gooding (I-Gooding) methods (Yin et al. 2023). While the I-Gooding method shows improvement over the original Gooding method, both algorithms demonstrate limited performance in the test scenarios which include a mix of LEO, MEO, and GEO objects. Gooding-based methods are known to be sensitive to measurement errors, and their orbital constraint effectiveness diminishes for MEO and GEO objects over short observation arcs, leading to degraded stability. In our statistical analysis, the probability of large errors (SMA error > 500 km, as a proportion of errors > 200 km) for the Gooding and I-Gooding methods reached 90.9% and 88.0%, respectively. In contrast, the erroneous solution convergence ratio for the original AR method was 88.8%. The proposed RV-SPSS method, through multi-direction solution set variance screening and slant range threshold exclusion, drastically reduced this ratio to below 5.5%. Simultaneously, it increased the percentage of SMA errors less than 100 km from 69.1% to 80.0% and less than 200 km from 76.6% to 90.5%. These results demonstrate that the proposed method, by specifically targeting geometry-induced erroneous solutions, offers superior robustness and accuracy compared to both the original AR method and recent Gooding-based approaches in a space-based very-short-arc IOD context.

To evaluate the robustness of the proposed RV-SPSS method under more challenging observation conditions, we performed additional experiments with varying levels of measurement noise. Three additional noise levels (3″, 4″, and 6″) were tested, while keeping all other simulation parameters unchanged. The results are summarized in Table 5. It can be seen that the erroneous-solution convergence ratio remains below 7% for 3″ noise and below 20% even for 6″ noise, while the percentage of semi-major axis errors below 200 km stays above 87% for 3″ noise and above 69% for 6″ noise. These sensitivity tests confirm that the RV-SPSS algorithm maintains its effectiveness under degraded data quality, demonstrating its practical applicability in realistic space-based surveillance scenarios.

The experimental results clearly indicate that the erroneous solution problem in the AR algorithm for space-based observation is pervasive. As shown in Table 6, without erroneous solution suppression measures, the AR algorithms convergence ratios to erroneous solutions in random simulation experiments for LEO, MEO, and GEO objects were 52.8%, 64.6%, and 17.7%, respectively. This data confirms that even for objects with randomly distributed orbital parameters, due to the complex relative motion geometry between the space-based observation platform and the objects, the generation of mathematically valid but physically unrealistic erroneous solutions by the AR algorithm under angle-only observation conditions is a frequent occurrence, particularly more prominent in MEO and LEO orbit scenarios.

Further analysis of the distribution of erroneous solutions, as illustrated in Fig. 12 using large-scale random orbit simulations (300 objects across LEO, MEO, and GEO), provides direct experimental validation of the core theoretical principle established in Section 2.3—namely, that true solutions are invariant under geometric variations, while erroneous solutions are highly sensitive to such changes. This theoretical distinction is clearly reflected in Fig. 12 from two perspectives. First, examining the horizontal axis (semi-major axis values), erroneous solutions are widely dispersed across a remarkably broad range, from below 10,000 km to over 46,000 km. This dispersion validates the theoretical prediction that erroneous solutions are geometry-dependent: different targets have inherently different observation geometries, and the resulting erroneous solutions therefore take on correspondingly different values. Even for GEO objects, whose orbital altitudes are relatively concentrated, erroneous solutions span an interval from 10,000 km to 46,000 km, demonstrating that these erroneous solutions cannot be eliminated by simple numerical constraints based on prior orbit type assumptions. Second, examining the vertical axis (relative variance of solution sets), correct solutions consistently cluster in the low relative variance region (below 10^–6), demonstrating their stability across different geometric configurations. In contrast, erroneous solutions predominantly concentrate in the high relative variance region (above 10^–5). This distribution directly validates two key theoretical claims: (1) the relative variance criterion serves as an effective discriminator for distinguishing true solutions from erroneous ones, and (2) erroneous solutions are indeed sensitive to changes in observation geometry, manifesting as high statistical dispersion when evaluated across different geometric configurations. Critically, these patterns hold consistently across LEO, MEO, and GEO objects in the large-scale random simulation, confirming that the theoretical principle—and the screening strategy based upon it—is broadly applicable rather than limited to specific orbital regimes.

Fig. 12. Scatter plot of solution set distribution from random orbit simulation.

Download Original Figure

This inherent geometry-dependence of erroneous solutions means that relying solely on prior orbit type assumptions and setting numerical constraints is insufficient to effectively identify and eliminate them. Therefore, the RV-SPSS method proposed in this paper, by constructing multi-geometry-configuration solution sets and discriminating based on their statistical stability, provides a theoretically grounded and more effective approach to solving such erroneous solution screening problems.

The proposed RV-SPSS method demonstrated excellent erroneous solution suppression capability and algorithm robustness in the random simulation experiments. After applying this method, the erroneous solution convergence ratios for LEO, MEO, and GEO objects significantly decreased to 5.5%, 2.0%, and 5.4%, respectively. This improvement indicates that by constructing multi-geometry-configuration solution sets and using their relative variance for stability discrimination, erroneous solutions arising from changes in observation geometry can be effectively identified and eliminated. It is particularly noteworthy that the erroneous solution ratio for MEO objects drastically dropped from 64.6% to 2.0%, demonstrating the methods strong adaptability and robustness to the complex relative motion geometry in the medium Earth orbit region. Similarly, the method shows stable erroneous solution elimination effects for LEO and GEO objects, further verifying its wide applicability across different orbital altitudes and motion characteristics.