Uncertainty Analysis of the Short-Arc Initial Orbit Determination

The solution of the short-arc angles-only orbit determination problem has large uncertainty because the topocentric range is not observable. For a certain angular observation tracklet with measurement noise, there exist numerous potential orbits, all of which are compatible with the observations. However, the solution of a deterministic initial orbit determination algorithm is usually far different from the true orbit, especially for the semimajor axis and the eccentricity. A new sampling method is proposed to describe the probability distribution of the orbit determination solutions. Firstly, a series of orbits are sampled in the semimajor axis - eccentricity plane. A chi-square test method is proposed to select candidate orbits from the sample orbits. The weights of the candidate orbits are calculated to measure their probability being the true orbit. Finally, the kernel density estimation algorithm is used to estimate the probability density function of the true orbits. With some a priori assumptions, the candidate orbits can be further screened, and their weight can be modified. The a priori knowledge can significantly improve the accuracy of the orbit determination solution.


I. INTRODUCTION
The initial orbit determination (IOD) is an old but developing problem, which plays an important role in Space Domain Awareness (SDA). For an unknown resident space object (RSO) firstly detected by an optical or radar sensor, the IOD process can give a preliminary guess of the true orbit, providing valuable guidance information for the follow-up observations. For a known RSO, the IOD process is useful for object identification and tracklet association. Generally, the IOD process is the pre-stage of orbit determination (OD) using a batch or sequential estimation process. Traditional IOD algorithms only use three pairs of measurements. However, modern sensors can gather many pairs of measurements within a few seconds. In this paper, the IOD problem means that the arc used for orbit determination is short, even if there are more than three pairs of measurements are generated.
By the end of 2019. 10.04, there are 19779 RSOs larger than 10cm, including 5181 spacecrafts and 14598 rocket bodies and debris [1]. Far more debris smaller than 10cm exists in the space. Recently, many mega-constellation projects The associate editor coordinating the review of this manuscript and approving it for publication was Nagarajan Raghavan . have been proposed by some companies and institutions. For instance, SpaceX's Starlink constellation will consist of a total of 11943 satellites in the low earth orbit (LEO) [2]. The rapid increase of the RSOs brings unprecedented pressure to the SSA systems. The optical sensors will play a more critical role. To ensure the detection of all RSOs, sensors will produce millions of short-arc observations. ''Short arc'' usually means that the observation duration is less than 1% of the orbital period. The lack of range information leads to great uncertainty of the IOD solution. For a certain tracklet, there are numerous potential orbits compatible with the observations. However, the IOD algorithms usually give a bad solution far from the true orbit. Therefore, instead of a point solution, the corresponding estimation of uncertainty must be given. From the probability point of view, a good IOD process should give the probability density function (PDF) of the true orbit. The PDF of the true orbit will help to perform data association, determine collision probabilities, and initialize a filter for the orbit improvement.
The short-arc orbit determination is also a basic problem in the astronomical field. To determine the orbit of poorly observed asteroids, astronomers developed a kind of orbit determination method based on the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ ''Admissible Region'' (AR) [3]- [7]. AR is defined in the topocentric range-range rate plane, in which each point represents an elliptical orbit. The true orbit of the RSO is searched in the Admissible Region. To estimate the uncertainty quantification for the IOD process, the so-called ranging method is developed. In the ranging method, a series of candidate orbits are sampled by assuming the topocentric ranges (or range rate) at some epochs. Then the probabilities of these candidate orbits being the true orbit are calculated. The probability distribution of the true orbit can be estimated by the candidate orbits. There are two kinds of ranging methods: statistical ranging method [8]- [10] and systematic ranging method [11]- [13]. In the statistical ranging method, the topocentric ranges at two epochs of the observations are sampled. Two pairs of observations with assumed ranges determine a sample orbit. A weight is given to each sample orbit by comparing the errors between virtual observations of the sample orbit and the true observations. The systematic ranging method adopts the technique of the Admissible Region. Candidate orbits are sampled in the range-range rate phase instead of assuming a pair of ranges.
Ansalone et al. [14] and Hinagawa et al. [15] applied the ranging method to the determination of the satellite orbit. Their approaches are similar to the statistical ranging methods. Two assumed ranges determined an orbit, whose virtual observation errors are also determined. They find the optimal IOD orbit by minimizing the observation errors using the genetic algorithm. However, their method only finds a point solution. The uncertainty characteristic of the solution is not analyzed.
Another kind of method to analyze the IOD uncertainty uses the PDF mapping technology. Given the statistical characteristics of the measurement errors, all possible observations constitute an observation domain. Correspondingly, all possible IOD solutions constitute an orbit domain. The IOD process can be regarded as a nonlinear transformation from the observations to the orbit solution. Once the PDF of the observations is given, the PDF of IOD solutions can also be obtained.
Weisman et al. [16], [17] adopted the transformation of variables (TOV) technique to analyze the uncertainty quantification of the IOD solutions obtained by the Herrick-Gibbs approach. The TOV technique allows for the exact mapping of the PDF from the observation domain into the orbit domain. Armellin et al. [18], [19] realized the observationsto-solution (O2S) mapping using the Differential Algebra (DA) tools. The unscented transformation (UT) technique is also used in the uncertainty analysis for the IOD problem [20]. Given the measurement error characteristics, the UT technique is used to estimate the mean and covariance of the solutions. Then the PDF of the solution is described by the Gaussian distribution. Binz et al. [21]- [23] used similar methods and obtained useful results from the ''binary sensors'', demonstrating that low-quality observations generated by low-cost sensors are also useful. However, the solution domain is not necessarily Gaussian after the nonlinear IOD processes, even if the observation domain is Gaussian. This will limit the accuracy of the UT methods.
The rest of this paper is organized as follows. Firstly, we proposed a Laplace-least square (Laplace-LS) IOD method, whose estimation variance is very close to the Cramer-Rao Lower Bound (CRLB). Section 3 describes a new IOD uncertainty analysis approach. A series of orbits are sampled in the semimajor axis-eccentricity plane instead of the range-range rate plane. In section 4, we test our method on three RSOs: a highly elliptical orbit (HEO) satellite, a low earth orbit (LEO) satellite, and a Geosynchronous (GEO) satellite. The a priori knowledge is used to modify the solutions later. Some remarks conclude the paper in section 5.

II. THE LAPLACE-LEAST SQUARE INITIAL ORBIT DETERMINATION
Many methods can solve the angles-only initial orbit determination problem. The Laplace method, the Gauss method, the Double-r method, and the Gooding method are most commonly used [24]- [26]. However, these methods choose only three pairs of observations, with many observations wasted. The improved Laplace method [27], [28] utilizes all observations, which is essential to reduce the influence of measurement error and improve the reliability of the IOD solutions.
The improved Laplace method solves the conditional equations to find an orbit most compatible with the true observations. However, the original statistical characteristics of the measurement errors are distorted by a series of nonlinear transformations. From the perspective of optimal estimation, the improved Laplace method cannot give the optimal IOD solution. Moreover, when the observation arc is very short, the improved Laplace method will give a large-error solution, or even fail to converge.
In this paper, we use the least square method to solve IOD problems. By minimizing the measurement errors directly, this method can give a better solution than traditional IOD methods.

A. ALGORITHM DESCRIPTION
The Laplace-LS IOD method is simple and direct: to minimize the errors between the virtual observations of solution orbit and true observations. Assuming that a sensor has obtained m groups of noisy measurements of right ascensions and declination: The actual observation sequence with measurement noise is defined asỸ Correspondingly, the theoretical observation sequence of the true orbit (without noise) is defined as Supposing the measurement errors follow the normal distribution with variance σ 2 , the observed ascensions and declinations have following statistical characteristics: For any orbit whose state is x 0 at the initial epoch t 0 , there is also a theoretical observation sequence calculated by the dynamic model and the observation model: The weighted square sum of the errors between the theoretical observations and the actual observations indicates the difference between the given orbit x 0 and the actual orbit: The root mean square error (RMSE) of the errors is defined as the target function J (x 0 ), which indicates the average angle error. Then the IOD problem can be transformed into an optimization problem:

B. INITIAL VALUE AND OPTIMIZATION
A reasonable initial guess is very important for an optimization problem. We adopt the same initial value used in the improved Laplace method for our optimization problem. In the improved Laplace method, the conditional equations are comprised of 3m equations: where x 0 , y 0 , z 0 ,ẋ 0 ,ẏ 0 andż 0 are Cartesian components (position and velocity) of the orbit state x 0 . R xi , R yi , and R zi are position components of the sensor at epochs t i . F i and G i are the f and g functions of the RSO at epochs t i . They are decided by x 0 and t i .
Though the f and g functions are changing with x 0 and t i , they are easy to be approximated from their series expansion [27]: Using the approximate f and g functions, (8) comes to linear equations, whose solution is a good initial guess of the optimal solution.
Given an initial guess, there are many off-the-shelf algorithms to solve the optimization problem in (7). We adopt the Levenberg-Marquardt algorithm, which is an efficient and robust algorithm to solve non-linear least-square problems [29].

C. THE CRAMER-RAO LOWER BOUND OF INITIAL ORBIT DETERMINATION PROBLEM
Broadly speaking, the IOD problem is a kind of parameter estimation problem. An IOD algorithm is an estimator of the orbit state x 0 . The Cramer-Rao Lower Bound (CRLB) gives a lower bound for the variance of any unbiased estimator. If the variance of an estimator is close to the CRLB, it is considered a good estimator.
The CRLB is expressed as an inequality: wherex andx are the true value and the estimation of the parameter x respectively. P is the variance of an estimator. F is called the Fisher information, which is expressed as: whereỹ is the actual value of the observation. Specifically, the Fisher information of the IOD problem is (see Appendix)

D. ACCURACY TEST OF THE LAPLACE-LEAST SQUARE IOD METHOD
This section will test the accuracy of the Laplace-LS IOD method by simulation. As a contrast, the accuracy of the improved Laplace method is also analyzed. The simulation parameters are shown in Table 1. The RSO's orbit is described by the Keplerian elements, where a is the semimajor axis, e is the eccentricity, i is the orbit inclination, is the right ascension of ascending node (RAAN), ω is the argument of perigee and f is the true anomaly. The Monte Carlo simulation is necessary to measure the performance of the IOD algorithms because the measurement errors are stochastic. The Laplace-LS IOD algorithm and the improved Laplace method were implemented 1000 times respectively. The covariance of the IOD solutions is calculated and compared with the CRLB respectively. For convenience, only the standard deviations of the orbit elements (square roots of the diagonal elements of the covariance matrix) are listed and compared, as Table 2 shows. The statistical results show that the variance of the Laplace-LS IOD method is close to the CRLB, while the improved Laplace method performance poorly. Figure 1 shows the errors of  the semimajor axis and eccentricity in the 1000 times' simulation. Moreover, the CRLB shows that the large errors of semimajor axis and eccentricity are inevitable in short-arc circumstances, no matter how accurate an estimator is.

III. ORBIT SAMPLING IN THE SEMIMAJOR AXIS -ECCENTRICITY PLANE
The difficulty of short-arc IOD is that there are numerous orbits compatible with the observations. Due to the measurement noises, the least square (LS) orbit (solution of the Laplace-LS IOD method) is not necessarily the true orbit, though it has the smallest angle RMSE. On the contrary, the least square orbit is usually very different from the true orbit, or even an impractical orbit intersecting with the earth. Fortunately, though a and e are hard to determine, the inclinations and the RAANs are easy to determine because they determine the orbital plane's direction in the inertial space.
For the short arc IOD problem, a point solution given by a deterministic IOD algorithm is not credible. More candidate orbits which might be the true orbit must be found to describe the distribution of the true orbit comprehensively. In other words, a good IOD algorithm should give the PDF of the true orbit. In this paper, we propose an orbit sampling method in the semimajor axis -eccentricity plane. The idea of this method is as follows.
1) Sample the semimajor axis and eccentricity uniformly in the (a, e) plane. For each sample point, find a least square orbit with a and e constrained. The least square orbit given a and e is defined as the ''sample orbit''; 2) Set a threshold for the angle RMSE. If the angle RMSE is lower than the threshold, the sample orbit will be accepted as a candidate orbit. Otherwise, it will be rejected; 3) Give a weight for each candidate orbit to measure its probability being the true orbit; 4) Estimate the probability density function using the candidate orbits.

A. THE ORBIT DETERMINATION WITH GIVEN SEMIMAJOR AXIS AND ECCENTRICITY
When the semimajor axis and eccentricity are constrained, it is also probable to find an orbit most consistent with the observations. The optimization variables are reduced to The given semimajor axis, eccentricity, and the solution of this optimization problem constitute the sample orbit.

B. ACCEPTANCE THRESHOLD
For a sample orbit, a threshold is needed to determine whether to accept it as a candidate orbit. By analyzing the error distribution of the theoretical observations, we calculate the threshold based on a χ 2 value. Using this threshold, the probability of accepting the true orbit is greater than 1-α (α is the significance level).
In practical, only one specific observation sequencẽ Y = [α 1 ,δ 1 ,α 2 ,δ 2 , . . . ,α m ,δ m ] can be obtained. The errors between the theoretical true observationsȲ and the actual observationsỸ can be measured by J (x 0 ). Though the true value is unknown, the statistical characteristics of these errors are clear according to (4): Suppose the measurement errors are independent of each other. The sum of the squares of these errors follows the χ 2 -distribution with 2m degrees of freedom: Given a significance level α, the α quantile of χ 2 can be denoted as χ 2 2m (α), and P(χ 2 < χ 2 2m (α)) = 1 − α Then an inequality about J (x 0 ) can be obtained: Therefore, the acceptance threshold can be set to (χ 2 2m (α)/2m) 1/2 σ . When the angle RMSE of a sample orbit is smaller than this threshold, it is accepted as a candidate orbit. This judgment strategy makes sure the true orbit being accepted with a probability of 1-α.

C. WEIGHTS OF THE CANDIDATE ORBITS
For a certain candidate orbit x 0 , the posterior probability density of this orbit can be calculated by the Bayes formula: When no a priori knowledge is available, the distribution of x 0 is uniform. According to the expression of p(Ỹ |x 0 ) (see Appendix), Define the weight of the candidate orbit as w(x 0 ) is proportional to the probability density. Obviously, the least square orbit with the smallest J has the largest weight, i.e., the largest probability.
When the a priori knowledge is available, the weight can be modified. For instance, if an RSO is identified as a GEO object [30], there will be a constraint on the semimajor axis. In this circumstance, we can assume that the semimajor axis of the true orbit follows a Gaussian distribution. The mean is the standard semimajor axis of the GEO a GEO and the standard deviation is σ a . The a priori weight w prior can be defined as Then the sample orbit's weight should be modified to:

D. PROBABILITY DENSITY FUNCTION COMPUTATION
Some candidate orbits and the corresponding weights are obtained after the sampling procedure. They can be used to estimate the probability density function of the IOD solutions. There are many methods to estimate the PDF of the random variables from samples. They are usually classified as parameter estimation methods, semi-parametric estimation methods and non-parametric estimation methods. Due to the nonlinear characteristics of the IOD algorithms, the PDF of the IOD results should not be assumed to follow a certain parametric distribution.
Non-parametric estimation methods can give the PDF directly from the samples. Kernel density estimation (KDE) is a widely used non-parametric density estimation algorithm [31], [32]. It was proposed by Rosenblatt (1955) and Emanuel Parzen (1962), also known as the Parzen Window.
KDE for the one-dimensional random variable x is as follows. Assuming that x 1 , x 2 , . . . , x n are n independent samples of x, then the PDF f (x) can be estimated as: where K (·) is the kernel function, including Gaussian function, Epanechnikov function, Biweight function, etc., and h is the window width. When the samples have different weights, 1/n should be adjusted to each sample's weight w i : For the multi-dimensional random variable X = (x 1 , x 2 , . . . , x d ) T with n independent samples X i = (x i1 , x i2 , . . . , x id ) T (i = 1, 2, . . . , n), the KDE of the PDF is: where h = (h 1 , h 2 , . . . , h d ) is the window width vector. In this paper, the Gaussian function is adopted as the kernel function:

IV. RESULTS AND CONCLUSIONS
Our sample method is tested on three different types of RSOs. The first one is the Radiation Belt Strom Probes-A (RBSP-A), an HEO satellite. The second one is FengYun 3A, an LEO weather satellite. The third one is CHINASAT 9, a GEO communications satellite. Detail orbit parameters are shown in Table 3. The true satellites' orbits are obtained using the Two Line Elements (TLE). The optical observations are simulated assuming the sensor is located in China National Astronomical Observatories (XingLong Observatory).
For the LEO satellite, the observation's duration is 30 seconds with the measurement frequency of 1Hz. The theoretical measurements calculated by the true orbit and sensor's location constitute the true observations. Using the true observations, the solution of IOD algorithms must be the true orbit.
The Gaussian noises with a standard deviation of 3 are added to the true observations. The noisy observations are defined as the actual observations, which are the realization of the measurements with noise. The IOD algorithm and sampling method are implemented on the actual observations.
For the HEO satellite, the observation's duration is set to 7.5 minutes with the frequency of 1/15 Hz. The measurement resolution is set to 0.5 . For the GEO satellite, the observation's length is set to 15 minutes with the frequency of 1/30 Hz. The measurement resolution is also 0.5 .

A. CONFIDENCE BOUNDARY, CONFIDENCE REGION, AND SAMPLE ORBITS
Implement the sampling method in section 3 on the HEO satellite firstly. Figure 2 (a) shows sample orbits whose angle RMSE is smaller than 0.0004 • (1.5 ). They are all very compatible with the actual observations. The red circle represents the least square orbit that has the smallest error. The green diamond represents the true orbit. In Figure 2(a), there is a region where the gradient value is almost zero. All orbits in this region are candidate orbits. Taking J (x 0 ) = (χ 2 2m (α)/2m) 1/2 σ as the boundary of this region, the region will contain the true orbit with a probability of 1-α.
In this paper, the boundary J (x 0 ) = (χ 2 2m (α)/2m) 1/2 σ is defined as the (1-α)− confidence boundary. The corresponding region is defined as the (1-α)−confidence boundary. Figure 2(b) shows the contour plot of Figure 2(a) and some candidate orbits. The significance level α is set to 0.005. The red line represents the corresponding 99.5%-confidence boundary, whose value is 0.0001714 • (0.617 ). The color of each candidate orbit is based on the weight, which has been scaled in the range [0,1].
For the LEO satellite and the GEO satellite, the 99.5%-confidence boundary, 99.5%-confidence region, and the candidate orbits are shown in Figure 3 and Figure 4. For orbits with small eccentricity, the joint distribution of a and e of the candidate orbits is non-Gaussian because the eccentricity is limited in the range [0,1]. For an HEO satellite, the distribution of the candidate orbits may be Gaussian, as Figure 2(b) shows.

B. VALIDATION: THE MONTE CARLO METHOD
This paper uses the weighted candidate orbits to describe the probability distribution of the true orbit. The feasibility of this method can be validated by the Monte Carlo method.
The actual observations are the realization of measurements with noise. In reality, the theoretical observations of the true orbit are unknown, but their probability distribution is clear from (15): In the Monte Carlo method, the unknown theoretical observations of the true orbit are treated as random variables. Their means are the actual observations with noise, and the variance is the same as the measurement errors. Therefore, if we add the Gaussian noises to the actual observations and implement the Laplace-LS IOD algorithm, the solution (the Monte Carlo orbit) may be the true orbit. Repeat the process for some times, all IOD solutions can describe the probability distribution of the true orbit. If our sampling method in the  a-e plane is reasonable, the PDF estimated by the weighted candidate orbits and the Monte Carlo orbits should be similar to each other.
In the validation section, the orbit is described with two kinds of elements, namely the Keplerian elements K = [a e i ω f ] and the Cartesian elements X = [x y z v x v y v z ]. x, y, and z are three position components of the RSO in the International celestial reference frame (ICRF). v x , v y , and v z are the corresponding velocity components. The Cartesian elements have been normalized in this paper. The unit length is 6378.1366km, and the unit vector is 7.905km/s. Figure 5 shows the 2000 Monte Carlo orbits and the weighted candidate orbits projected on different planes. The Keplerian elements are projected on the a−e, i− , and ω−f plane respectively. The Cartesian elements are projected on the x − y, z-v x , and v y -v z plane respectively. The colors of the candidate orbits represent their weights. The weights of the Monte Carlo orbits are the same, so they are all represented by the gray point.
An important conclusion can be summarized in Figure 5. The Keplerian elements of the IOD potential true orbit are not necessary Gaussian because of the elements' value are restrained. The eccentricity is restrained in the range [0,1], VOLUME 8, 2020      we only have to compare the mean and variance calculated by the Monte Carlo orbits with those calculated by the weighted candidate orbits. Table 4 -Table 6 show the validation results for three RSOs. The numbers of candidate orbits for the HEO, LEO, and GEO satellites are 432, 515 and 788 respectively. The numbers of VOLUME 8, 2020 Monte Carlo orbits are all 2000 for three satellites. Comparison results of the means and standard deviations show that the sampling method in this paper is feasible. Moreover, fewer weighted candidate orbit is needed than the unweighted Monte Carlo orbits to estimate the PDF of the true orbit.

C. MODIFICATION OF THE CANDIDATE ORBITS -ADDING THE A PRIORI KNOWLEDGE
The candidate orbits selected and weighted by the angle RMSE J (x 0 ) do not contain any a priori knowledge. In this section, we will discuss how the a priori knowledge can help to modify the candidate orbits, including further screening and the modification of the weights. Figure 5(c) shows that the semimajor axis of some candidate orbits are smaller than the earth radius R e , which is not practical. We should exclude the orbits if they intersect with the earth. The geocentric range of the perigee of an orbit is r p = a(1-e). The condition of rejecting a candidate orbit can be set to:

1) CANDIDATE ORBIT SCREENING
Apply this rejection condition to three test cases. Figure 6 and Figure 7 show the remaining candidate orbits after the screening process in the HEO and LEO test cases. For the 515 HEO candidate orbits, 356 orbits are rejected and 159 orbits are retained. For the 432 LEO candidate orbits, 356 orbits are rejected and 159 orbits are retained. All 788 GEO candidate orbits are retained. Besides, the least square orbit for the HEO and GEO cases are all rejected.

2) WEIGHTS MODIFICATION
Take the GEO satellite for example. Usually, the GEO RSOs are easy to distinguish due to the special characteristics of their orbit. Then the weight can be modified using (23). Setting a GEO = 42166km and σ a = 100km, the modified candidate orbits are shown in Figure 8, whose weights were changed significantly.

3) THE MAXIMUM A POSTERIORI ORBIT
Both the candidate orbit screening and the modification of the weights will lead to the change of the posterior probability density. In this circumstance, the best IOD point solution is no longer the least square orbit, but the maximum a posteriori (MAP) orbit. Equivalently, we can find the MAP orbit by solving the following optimization problem.
For the a priori constraint of the GEO orbit, the optimization problem turns into: For the orbit screening case, w prior should be defined as: w prior = 1, if the condition is met 0, if the condition is not met (32) However, this definition will cause the discontinuity of the objective function. For the convenience of optimization, the screening condition should be transformed into a continuous function. Taking the constrain of r p as example, w prior can be expressed as: When r p is slightly larger than R e +100, w prior approximately equal to 1. When r p is slightly smaller than R e +100, w prior approximately equal to 0. Actually, many similar functions can be defined such as w prior =(1+tanh(r p -R e -100))/2. However, the definition of (33) has an additional benefit. This function is in the exponential form, which helps the simplification of the optimization problem:  the a priori knowledge can significantly improve the accuracy of the IOD point solution.

D. PROBABILITY DENSITY FUNCTION: ESTIMATION AND VALIDATION
This section will estimate the PDF of the IOD result using KDE and the candidate orbits. For the orbit screening process, the orbit weights are unchanged. The PDF estimation result can still be validated by the Monte Carlo orbits if they are screened by the same rejection condition. Figure 9 and Figure10 show the projection of the PDF on the two-dimensional space. The diagonal sub-graphs show the marginal distribution of the orbit state elements x, y, z, v x , v y , and v z , respectively. Other sub-graphs show the joint distribution of two different orbit elements. In the marginal distribution sub-graphs, the solid blue lines show the PDF of a certain variable estimated by the screened candidate orbits. The black dotted lines show the PDF estimated by the screened Monte Carlo orbits. The green dash-dot lines show the value of the true orbit parameters. The magenta dashed line shows the value of the MAP orbit parameters. In the two-dimensional joint distribution sub-graphs, the green diamonds show the position of true orbit, and the magenta stars show the position of the MAP orbit.
The results show that PDF estimated by the candidate orbits and the Monte Carlo orbits are almost identical, indicating the KDE method is feasible to estimate the PDF of the IOD results.
For the GEO satellite, we estimate the PDF of its orbit without and with the a priori knowledge respectively, as Figure 11 shows. The red circles and the bold dashed red line represent the LS orbits. The comparison of Figure 11(a) and (b) show that the a priori knowledge can reduce the IOD uncertainty significantly. The MAP solution is far more accurate than the LS solution.

V. CONCLUSION
In this paper, the problem of short-arc orbit determination is researched. A more accurate IOD method, namely the Laplace-LS orbit determination method, is presented, whose estimation variance is close to the CRLB. For the short-art initial orbit determination problem, a new sampling method in the semimajor axis-eccentricity plane is proposed to describe the probability distribution function of the true orbit. Some candidate orbits are obtained and weighted according to the angle errors. They are efficient to estimate the PDF of the true orbit. When the a priori knowledge is available, the candidate orbits can easily be screened and modified. Then the accuracy of the IOD solution can be significantly improved.
When the detected trajectory is very short or the sensor has limited accuracy, our method can be applied to determine the RSO's orbit comprehensively.
Our method can also be applied to the long-arc IOD problem. However, the uncertainty of the long-arc IOD result will degrade dramatically. A point solution given by a traditional algorithm is usually credible. In this circumstance, methods introduced in section I, such as the uncertainty analysis using the unscented transformation, maybe more straightforward.
Finally, the Fisher information of the initial orbit determination problem is: Her research interests include spacecraft orbit dynamics, the design of constellation, and space mission analysis.
PING JIANG was born in 1996. He received the B.Eng. degree from the University of Science and Technology of China (USTC), Hefei, China, in 2018. He is currently pursuing the master's degree with Space Engineering University, Beijing, China. His research interests include correlation and the orbit determination of space objects. VOLUME 8, 2020