Random Error Estimation and Propagation Analysis for Satellites’ Initial Positions

Abstract

An increasing number of satellites are currently being launched into orbit to work in the form of clusters or constellations. However, the initial orbit position is accompanied by random errors, which will propagate during their running. Therefore, the orbit precision of the satellites directly affects space safety, network accuracy, and operating efficiency. Hence, accurate and fast random error estimation is essential to improve satellite control. The traditional method will take much time and cost, and it is associated with complex calculations or low accuracy, especially for large-scale constellations. In this paper, a random error evaluation model based on the ellipsoid is proposed. It can be used to estimate initial positions and error propagation for any orbit satellites. By comparing with the experiment results using the Monte Carlo method, it is proved that the proposed model is relatively simple, highly effective, and good at accuracy.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Positioning the satellite into orbit is accompanied by the uncertainty that spreads during the satellite running and affects its efficiency. Random and systematic errors mainly cause the uncertainty of orbit positioning during satellite production, design, manufacturing, and launch. Systematic errors can be reduced, but random errors are caused by accidental factors that are difficult to predict and control. In practice, satellites with high requirements, high accuracy, and high value are usually corrected or adjusted when they deviate from the normal orbital position. However, since their correction ability is limited, the cost of orbit correction may be much higher than that of satellites for small or micro-nanosatellites. Moreover, it is also difficult to conduct extensive orbit correction for large-scale satellite clusters. However, accurately and quickly estimating the random errors of the orbital position can help calculate the position error distribution of the satellite in orbit, evaluate the possibility of collision, and analyze the deviation of coverage or communication area of a satellite, constellation, or cluster. Consequently, correct control strategies can be specified.

Space target monitoring and orbit theory analysis are the main methods for determining the orbit of the satellite and the orbit error. Space target monitoring discovers, tracks, and measures the motion parameters of space targets in real-time through various observation methods, determining their orbital characteristics. The US Space Surveillance Network (SSN) continuously updates orbit data using the SGP4 (simplified general perturbations)/SDP4 (simplified deep space perturbations) orbit theoretical model by cataloging and tracking the satellites in orbit [1]. NORAD (North American Aerospace Defense Command) also regularly publishes TLE (two-line elements) of most space targets by monitoring them. These orbit determination methods require a perfect observation system and an accurate and efficient prediction model. However, they have the disadvantages of poor reliability, high hardware cost, and long orbit determination time. Nevertheless, TLE may become private and make these methods unreliable.

Some experts have proposed convenient methods for different types of satellites. Knogl et al. [2] used LEO (low earth orbit) satellite communication channels to locate GEO (geosynchronous earth orbit) satellites accurately. Vighnesam et al. [3] investigated the systematic method of determining the operating orbit of IRS (Indian remote sensing) satellites through the position difference. Yi et al. [4] and Liu et al. [5] used onboard GPS (global positioning system) and BDS2 (BeiDou satellite) to accurately determine the relative orbit of China’s TH-2 (Tian Hui) satellite formation. Li et al. [6] analyzed the accuracy of the BDS3 orbits using SLR (satellite laser ranging), while Gu et al. [7] analyzed the principle of GPS precise orbit determination evaluation by the same method. Existing investigations determine the orbital position of the satellite by receiving satellite signals from the measuring station [8, 9]. Larki et al. [10] determined the satellite orbit based on the gradient method. Based on historical orbit data, Bai [11] investigated the orbit prediction error and the collision probability of space objects. According to the literature review, most authors investigated satellite orbital positioning using other equipment. BDS or other equipment used in orbit determination may malfunction and are characterized by complex systems, high costs, and poor reliability.

Uncertainty can also be analyzed from the perspective of mathematics. Methods such as spatial estimation [12, 13], real positioning error and error model expression [14], description of spatial data error [15–20], data error propagation model [21–23], and data error detection and correction [24, 25] were proposed by mathematical statistics. These methods have high efficiency and low cost. However, few studies have been conducted on the description of errors in data and the applicability analysis in aerospace.

In recent years, many satellites have been launched into orbit in constellations, clusters, or groups. It is meaningful to study the orbital position laws of satellites subjected to the initial random error and find the certainty hidden in chance, such as determining the positioning accuracy for navigation satellites, calculating the coverage for observation satellites, analyzing the collision risk for space targets, and conducting mission planning for constellations.

In this paper, an ellipsoid model was proposed to estimate the random error and propagation analysis for the initial positions of satellites. The error model is provided by the analytical method, and the results are numerically verified. The remainder of this paper is organized as follows. In the second part, the background material of satellite dynamics is introduced. In part three, the position uncertainty of the satellite launch into the orbital position is calculated according to the uncertainty matrix of the six orbital elements of the satellite. Then, the three coordinate axis directions under the geocentric inertial coordinate system are calculated by the propagation rate of the covariance matrix. The theoretical model of the equal probability density surface is constructed, and the probability calculation formula is provided for the satellite launch into the orbit in a certain error range. Based on these values, the transfer of the satellite’s initial orbital position error is further derived. In part four, the orbital positions of LEO, MEO (medium earth orbit), and IGSO (inclined geosynchronous orbit) satellites are simulated under random errors using STK (satellite tool kit) software and the Monte Carlo method. The last part provides conclusions and research prospects.

2. Background Material

In this part, relevant symbols and concepts will be introduced. Furthermore, satellite space coordinate systems and satellite dynamics knowledge will be explained.

2.1. Concepts and Symbols

Mark the satellite as $S$, the center of the earth (geocenter) as $O$, and the distance from the satellite to the geocenter as $r$. It is considered that the satellite follows an elliptic orbit. The center of the ellipse is marked as $\tilde{O}$, and the point closest to the geocenter is called perigee, which is marked as $P_0$. The point farthest from the geocenter is called apogee, which is denoted as $B$. Half of the major axis of an elliptical orbit is called the satellite orbit semimajor axis, which is marked as $a$. The ratio of the distance between the two focuses of the ellipse and the major axis is called satellite orbit eccentricity, marked as $e$. The angle between the orbital plane and the earth’s equatorial plane is called orbital inclination, marked as $i$. The satellite orbit (except for that with zero inclination degree) has two focal points between the satellite orbit and the earth’s equatorial plane. The satellite running arc segment from the southern hemisphere through the equatorial plane to the northern hemisphere is called the ascending segment. The point crossing the equatorial plane is called ascending node and is marked as $N$. The angle between the vernal equinox, the ascending node, and the geocenter is called the right ascension of ascending node, marked as $\Omega$. The field angle between the perigee and the ascending node with respect to the geocenter is called the perigee argument and is marked as $\omega$. The included angle between the orbital perigee and the satellite and the geocenter is called the true anomaly, marked as $f$. The average angular velocity of the satellite is marked as $n$, and the time when the satellite passes the perigee is recorded as $t_0$. In addition, to express the motion of the satellite, the eccentric anomaly $E$ and the mean anomaly $M$ are introduced, which will be explained in Section 2.2.

2.2. Coordinate System with Respect to the Satellite

The coordinate system is crucial for accurately expressing the position and describing the motion of the satellite. In this research, the geocentric inertial coordinate system and the geocentric orbital coordinate system are used.

The geocentric inertial coordinate system selects the geocenter ($O$) as the coordinate origin. The $X$-axis points to the vernal equinox in the equatorial plane, the $Z$-axis coincides with the earth’s rotation axis, and the $Y$-axis forms the right-hand Cartesian coordinate system with the $X$-axis and the $Z$-axis in the equatorial plane. $\text{OXYZ}$ is the geocentric inertial coordinate system, as shown in Figure 1.

[figure(s) omitted; refer to PDF]

In the geocentric orbital coordinate system, the coordinate origin point is located at geocentric $O$, while the $\text{OU}$ and $\text{OV}$ axes are in the orbit plane of the satellite. Simultaneously, the $\text{OU}$ axis points to the perigee, and the $\text{OW}$ axis coincides with the normal vector of the satellite orbital plane. The axes $\text{OV}$, $\text{OU}$, and $\text{OW}$ form the right-hand Cartesian coordinate system in the satellite orbital plane, as shown in Figure 1.

For simplicity, an auxiliary circle is produced with the center $\tilde{O}$ and the radius $a$, as shown in Figure 2. Draw a vertical line of the $\text{OU}$ axis through $S$ and intersect with $\text{OU}$ at point $H$. Then, extend $\text{HS}$ intersection auxiliary circle to $\tilde{S}$, connect point $\tilde{O}$ and point $\tilde{S}$, and the included angle between line $\tilde{\text{O}}\tilde{\text{S}}$ and line $\tilde{\text{O}}\text{U}$ is the eccentric anomaly $E$. The angle between the line $\tilde{\text{O}}\tilde{\text{S}}$ and $\tilde{\text{O}}\text{U}$ is the true anomaly, marked as $f$. By Equation (1), point $C$ can be obtained. The included angle between the line $\tilde{\text{O}}\text{U}$ and $\tilde{\text{O}}\text{C}$ is the mean anomaly, which is different from the true anomaly and eccentric anomaly and does not have geometric interpretation. However, it can be directly used in the Kepler equation, as shown in the following: $$\begin{array}{}\text{(1)}& M=nt-t_0+M_0,\end{array}$$where $t_0$ is the start time, $t$ is the elapsed time, $M_0$ is the mean anomaly at $t_0$, and $n$ is the average angular velocity of the satellite expressed as follows: $$\begin{array}{}\text{(2)}& n=\sqrt{\frac{\mu }{a^3}},\end{array}$$where $\mu$ is the Kepler constant, whose value is $3.986\times 10^5\text{k}{\text{m}}^3/{\text{s}}^2$.

[figure(s) omitted; refer to PDF]

2.3. Equation of Satellite Motion

The satellite’s orbital coordinate is usually expressed by six orbital elements, written as the vector: $$\begin{array}{}\text{(3)}& \mathbf{X}={a,e,i,\Omega ,\omega ,f}^T.\end{array}$$

Regarding the satellite’s six orbital elements, only true anomaly changes with time, and its change equation can be expressed indirectly through the eccentric and mean anomalies. According to orbital dynamics, the relationship between the true anomaly, eccentric anomaly, and mean anomaly is as follows [26]: $$\begin{array}{}\text{(4)}& \begin{array}{c}\tan \frac{f}{2}=\frac{\sqrt{1+e}}{\sqrt{1-e}}\tan \frac{E}{2}\hfill \\ M=E-e\sin E\hfill \end{array}.\end{array}$$

The three-dimensional coordinate expression of a satellite position in the geocentric orbital coordinate system is the following: $$\begin{array}{}\text{(5)}& \begin{array}{c}U=a\cos E-ae,\hfill \\ V=a\sqrt{1-e^2}\sin E\hfill \\ W=0,\hfill \end{array},\end{array}$$where $U$ is the component of a satellite orbital position in the direction of the $\text{OU}$ axis, $V$ is the component of a satellite’s position in the direction of the $\text{OV}$ axis, and $W$ is the component in the direction of the $\text{OW}$ axis. The geocentric inertial coordinate system can be written as follows: $$\begin{array}{}\text{(6)}& \begin{array}{c}x\\ y\\ z\end{array}={\mathbf{R}}_3{\mathbf{R}}_2{\mathbf{R}}_1\begin{array}{c}U\\ V\\ W\end{array},\end{array}$$

where ${\mathbf{R}}_1$, ${\mathbf{R}}_2$, and ${\mathbf{R}}_3$ are rotation matrices expressed as follows: $$\begin{array}{c}\text{(7)}& {\mathbf{R}}_1=\begin{array}{ccc}\cos \omega & -\sin \omega & 0\\ \sin \omega & \cos \omega & 0\\ 0& 0& 1\end{array},{\mathbf{R}}_2=\begin{array}{ccc}1& 0& 0\\ 0& \cos i& -\sin i\\ 0& \sin i& \cos i\end{array},\\ {\mathbf{R}}_3=\begin{array}{ccc}\cos \Omega & -\sin \Omega & 0\\ \sin \Omega & \cos \Omega & 0\\ 0& 0& 1\end{array}.\end{array}$$

Hence, $$\begin{array}{}\text{(8)}& \begin{array}{c}x\\ y\\ z\end{array}=\begin{array}{ccc}\cos \Omega \cos \omega -\sin \Omega \sin \omega \cos i& -\cos \Omega \sin \omega -\sin \Omega \cos \omega \cos i& \sin \Omega \sin i\\ \sin \Omega \cos \omega +\cos \Omega \sin \omega \cos i& -\sin \Omega \sin \omega +\cos \Omega \cos \omega \cos i& \sin \Omega \sin i\\ \sin i\sin \omega & \cos \omega \sin i& \cos i\end{array}\begin{array}{c}U\\ V\\ W\end{array}=\begin{array}{c}a\cos E-e\cos \Omega \cos \omega -\sin \Omega \sin \omega \cos i-a\sqrt{1-e^2}\sin E\cos \Omega \sin \omega +\sin \Omega \cos \omega \cos i\\ a\cos E-e\sin \Omega \cos \omega +\cos \Omega \sin \omega \cos i-a\sqrt{1-e^2}\sin E\sin \Omega \sin \omega -\cos \Omega \cos \omega \cos i\\ a\cos E-e\sin i\sin \omega +a\sqrt{1-e^2}\sin E\cos \omega \sin i\end{array}.\end{array}$$

2.4. Expression for the Initial Orbital Position of a Satellite

Satellite orbital positions usually adopt different methods on various occasions and application backgrounds. If the impact of random factors during launching is considered, the satellite’s initial orbital position is a random variable. The mathematical expectation of its orbital elements vector can be written as follows: $$\begin{array}{}\text{(9)}& {\mu }_X=\begin{array}{c}Ea\\ Ee\\ Ei\\ E\Omega \\ E\omega \\ Ef\end{array}=\begin{array}{c}\tilde{a}\\ \tilde{e}\\ \tilde{i}\\ \tilde{\Omega }\\ \tilde{\omega }\\ \tilde{f}\end{array},\end{array}$$where $Ea$, $Ee$, $Ei$, $E\Omega$, $E\omega$, and $Ef$ are the expectation of the semimajor axis, eccentricity, inclination, right ascension of ascending, perigee argument, and true anomaly of the satellite orbit, respectively. Let the values of $Ea$, $Ee$, $Ei$, $E\Omega$, $E\omega$, and $Ef$ are $\tilde{a}$, $\tilde{e}$, $\tilde{i}$, $\tilde{\Omega }$, $\tilde{\omega }$, and $\tilde{f}$, respectively.

The variance of six orbital elements of the vector $\mathbf{X}$ can be expressed by ${\mathbf{D}}_{\text{XX}}$. $$\begin{array}{}\text{(10)}& {\mathbf{D}}_{\text{XX}}=\begin{array}{cccccc}{\sigma }_a^2& {\sigma }_{ae}& {\sigma }_{ai}& {\sigma }_{a\Omega }& {\sigma }_{a\omega }& {\sigma }_{af}\\ {\sigma }_{ea}& {\sigma }_e^2& {\sigma }_{ei}& {\sigma }_{e\Omega }& {\sigma }_{e\omega }& {\sigma }_{ef}\\ {\sigma }_{ia}& {\sigma }_{ie}& {\sigma }_i^2& {\sigma }_{i\Omega }& {\sigma }_{i\omega }& {\sigma }_{\mathit{\text{if}}}\\ {\sigma }_{\Omega a}& {\sigma }_{\Omega e}& {\sigma }_{\Omega i}& {\sigma }_{\Omega }^2& {\sigma }_{\Omega \omega }& {\sigma }_{\Omega f}\\ {\sigma }_{\omega a}& {\sigma }_{\omega e}& {\sigma }_{\omega i}& {\sigma }_{\omega \Omega }& {\sigma }_{\omega }^2& {\sigma }_{\omega f}\\ {\sigma }_{fa}& {\sigma }_{fe}& {\sigma }_{fi}& {\sigma }_{f\Omega }& {\sigma }_{f\omega }& {\sigma }_f^2\end{array}.\end{array}$$

The elements on the main diagonal in ${\mathbf{D}}_{\text{XX}}$ are variances of random variables $a,e,i,\Omega ,\omega ,$ and $f$, while the other elements express the covariance that describes the correlation between any two orbital elements. At the initial time, the six orbital elements are independently designed, i.e., the covariance is zero. Therefore, the variance matrix of the vector ${\mathbf{D}}_{\text{XX}}$ can be written as follows: $$\begin{array}{}\text{(11)}& {\mathbf{D}}_{\text{XX}}=\begin{array}{cccccc}{\sigma }_a^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_e^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_i^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{\Omega }^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{\omega }^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_f^2\end{array}.\end{array}$$

Since they are subjected to the random error of the orbital position of the satellite, the vector ${\mu }_X$ and the variance matrix ${\mathbf{D}}_{\text{XX}}$ are the digital characteristics of the 6-dimensional normal random vector. Once the satellite is launched into orbit, the true anomaly is affected by eccentricity, semimajor axis, and other perturbation factors. Therefore, the covariance is not zero.

2.5. Random Error Expression of the Initial Orbital Position in the Geocentric Inertial Coordinate System

Suppose that the expected satellite’s initial orbital position is $\tilde{A}\tilde{x},\tilde{y},\tilde{z}$ in the geocentric inertial coordinate system. However, the actual position is $Ax,y,z$, as shown in Figure 3.

[figure(s) omitted; refer to PDF]

The deviation value $\Delta$ of the actual initial position can be expressed as follows: $$\begin{array}{}\text{(12)}& \begin{array}{c}{\Delta }_x=x-\widetilde{x,}\hfill \\ {\Delta }_y=y-\tilde{y}\hfill \\ {\Delta }_z=z-\tilde{z},\hfill \end{array},\end{array}$$where ${\Delta }_x$, ${\Delta }_y$, and ${\Delta }_z$ are the components of the spatial position error $\Delta$ on the $X$, $Y$, and $Z$ axes, respectively. $$\begin{array}{}\text{(13)}& {\Delta }^2={\Delta }_x^2+{\Delta }_y^2+{\Delta }_z^2.\end{array}$$

According to the definition of mean and variance, $$\begin{array}{}\text{(14)}& \begin{array}{c}E{\Delta }_x^2=E{x-\tilde{x}}^2=E{x-E\tilde{x}}^2={\sigma }_x^2,\hfill \\ E{\Delta }_y^2=E{y-\tilde{y}}^2=E{y-E\tilde{y}}^2={\sigma }_y^2\hfill \\ E{\Delta }_z^2=E{z-\tilde{z}}^2=E{z-E\tilde{z}}^2={\sigma }_z^2.\hfill \end{array},\end{array}$$

According to the linear principle of the mean value and by combining Equations (13) and (14), $$\begin{array}{}\text{(15)}& E{\Delta }^2=E{\Delta }_x^2+E{\Delta }_y^2+E{\Delta }_z^2={\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2,\end{array}$$where $E{\Delta }^2$ is the theoretical average value of the square of the satellite’s initial position error. If the point-seat variance of $A$ is marked as ${\sigma }_A^2$, it can be expressed as follows: $$\begin{array}{}\text{(16)}& {\sigma }_A^2={\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2,\end{array}$$

where the magnitude of ${\sigma }_A^2$ is independent of the coordinate system and ${\sigma }_A$ is the satellite’s orbital position error [27].

3. Construction of Satellite’s Orbital Position Error Model

The satellite’s initial position variance ${\sigma }_A^2$ can be used to evaluate its accuracy, but it cannot show the directional difference. It is necessary to consider the directional deviation of the satellite’s orbital position when discussing the space collision safety, ground coverage, or efficiency of a satellite or a cluster. Therefore, it is necessary to determine the positional vector difference of the satellite in a three-dimensional space. The set of equal variance position points of the satellite’s initial position forms a closed surface representing errors in all directions. However, this surface is untypical and inconvenient to be drawn. However, its shape is similar to an ellipsoid, as shown in Figure 4. For convenience of analysis, the aforementioned shape is approximated by an ellipsoid, denoted as an error ellipsoid with equal variance (Figure 5). For simplicity, this ellipsoid will be referred to as an error ellipsoid. The parameters $U$, $V$, and $W$ are the three principal axes of the ellipsoid.

[figure(s) omitted; refer to PDF]

3.1. The Function for the Error Ellipsoid Model

The satellite’s initial position in orbit is presented as a normal three-dimensional distribution in the geocentric inertial coordinate system. Moreover, its density function of the joint distribution can be expressed as follows [28]: $$\begin{array}{}\text{(17)}& fx,y,z=\frac{1}{{2\pi }^{32}{{\mathbf{D}}_{\text{rr}}}^{12}}\exp -\frac{1}{2}x-\tilde{x},y-\tilde{y},z-\tilde{z}{\mathbf{D}}_{\text{rr}}^{-1}{x-\tilde{x},y-\tilde{y},z-\tilde{z}}^T,\end{array}$$

where $\tilde{x}$, $\tilde{y}$, and $\tilde{z}$ are the mathematical expectations of the initial position on the $X$, $Y$, and $Z$ axes. Respectively, ${\mathbf{D}}_{\text{rr}}$ is the uncertainty matrix of the satellite’s initial position in orbit under the geocentric inertial coordinate system, and ${\mathbf{D}}_{\text{rr}}$ is the determinant of the uncertainty matrix. The matrix ${\mathbf{D}}_{\text{rr}}$ can be expressed as follows: $$\begin{array}{}\text{(18)}& {\mathbf{D}}_{\text{rr}}=\begin{array}{c}\begin{array}{ccc}{\sigma }_x^2& {\sigma }_{x\text{y}}& {\sigma }_{\text{xz}}\end{array}\\ \begin{array}{ccc}{\sigma }_{\text{xy}}& {\sigma }_y^2& {\sigma }_{\text{yz}}\end{array}\\ \begin{array}{ccc}{\sigma }_{\text{zx}}& {\sigma }_{\text{zy}}& {\sigma }_z^2\end{array}\end{array},\end{array}$$where ${\sigma }_i^{2}i=x,y,z$ is the position variance of the coordinate axis direction and ${\sigma }_{ij}i,j=x,y,z$ is the covariance between the coordinate axes.

3.2. Calculation of Uncertainty Matrix $D_{\text{rr}}$

The covariance matrix ${\mathbf{D}}_{\text{XX}}$ of six orbital elements is expressed as Equation (11). Therefore, the uncertainty matrix ${\mathbf{D}}_{\text{rr}}$ can be calculated according to the covariance propagation rate [29]: $$\begin{array}{}\text{(19)}& {\mathbf{D}}_{\text{rr}}=\mathbf{A}{\mathbf{D}}_{\text{XX}}{\mathbf{A}}^T,\end{array}$$

where each element of the matrix $\mathbf{A}$ is the partial derivative approximate of position coordinates $x,y,z$ with respect to each variable of six orbital elements in the geocentric inertial coordinate system at the expected value ${\mu }_X$, as shown in Equation (9). Matrix $\mathbf{A}$ can be expressed as follows: $$\begin{array}{}\text{(20)}& \mathbf{A}={\begin{array}{cccccc}\frac{\partial x}{\partial a}& \frac{\partial x}{\partial e}& \frac{\partial x}{\partial i}& \frac{\partial x}{\partial \Omega }& \frac{\partial x}{\partial \omega }& \frac{\partial x}{\partial f}\\ \frac{\partial y}{\partial a}& \frac{\partial y}{\partial e}& \frac{\partial y}{\partial i}& \frac{\partial y}{\partial \Omega }& \frac{\partial y}{\partial \omega }& \frac{\partial y}{\partial f}\\ \frac{\partial z}{\partial a}& \frac{\partial z}{\partial e}& \frac{\partial z}{\partial i}& \frac{\partial z}{\partial \Omega }& \frac{\partial z}{\partial \omega }& \frac{\partial z}{\partial f}\end{array}}_{\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },\tilde{f}}.\end{array}$$

Elements in the matrix $\mathbf{A}$ are partial derivatives of the function with respect to each variable. These values are obtained using the approximate value $\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0$ in place of $a,e,i,\Omega ,\omega ,E_0$, where $\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0$ are constants. The accuracy is relatively high when the values $\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0$ are very close to their respective values $\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0$.

3.3. Calculating the Length of Error Ellipsoids Three Axes

According to the joint distribution density function described by Eq. (17), the same point of probability density in the three-dimensional normal distribution space can be expressed as follows: $$\begin{array}{}\text{(21)}& x-{\Delta }_x,y-{\Delta }_y,z-{\Delta }_z{\mathbf{D}}_{\text{rr}}^{-1}{x-{\Delta }_x,y-{\Delta }_y,z-{\Delta }_z}^T=k^2,\end{array}$$where $k$ is the amplification factor of the ellipsoid since different probability density points form various ellipsoids. According to linear algebra, Equation (21) is exactly the expression of a similar ellipsoidal family. If a value of $k$ is provided, an ellipsoid can be obtained, as shown in Figure 6.

[figure(s) omitted; refer to PDF]

Let $x-{\Delta }_x=\tilde{x},y-{\Delta }_y=\tilde{y},z-{\Delta }_z=\tilde{z}$. Then, Equation (21) can be written as follows: $$\begin{array}{}\text{(22)}& \tilde{x},\tilde{y},\tilde{z}{\mathbf{D}}_{\text{rr}}^{-1}{\tilde{x},\tilde{y},\tilde{z}}^T=k^2.\end{array}$$

Equation (22) expresses a similar ellipsoid family with $O$ as the ellipsoid center under the geocentric inertial coordinate system. Coordinate conversion is performed for convenience in analysis. Furthermore, the expression is obtained in the principal axis coordinate system $\text{OUVW}$, as shown in Figure 7.

[figure(s) omitted; refer to PDF]

In Figure 7, $\text{OU}$, $\text{OV}$, and $\text{OW}$ are the principal axes of the ellipsoid family. According to Equations (19) and (20), $D_{\text{rr}}$ is a real symmetric matrix, and there must be an orthogonal matrix $\mathbf{M}$ that satisfies the following [30]: $$\begin{array}{}\text{(23)}& {\mathbf{M}}^T{\mathbf{D}}_{\text{rr}}\mathbf{M}=\begin{array}{ccc}{\lambda }_1& & \\ & {\lambda }_2& \\ & & {\lambda }_3\end{array}.\end{array}$$

Equation (23) can be transformed, and an inverse is taken on both sides to obtain the following: $$\begin{array}{}\text{(24)}& {\mathbf{M}}^T{\mathbf{D}}_{\text{rr}}^{-1}\mathbf{M}=\begin{array}{ccc}{\lambda }_1^{-1}& & \\ & {\lambda }_2^{-1}& \\ & & {\lambda }_3^{-1}\end{array}.\end{array}$$

Therefore, $$\begin{array}{}\text{(25)}& {\mathbf{D}}_{\text{rr}}^{-1}=\mathbf{M}\begin{array}{ccc}{\lambda }_1^{-1}& & \\ & {\lambda }_2^{-1}& \\ & & {\lambda }_3^{-1}\end{array}{M}^T.\end{array}$$

The following expression is obtained by combining Equations (21) and (24). $$\begin{array}{}\text{(26)}& {{\mathbf{M}}^T\begin{array}{c}\tilde{x}\\ \tilde{y}\\ \tilde{z}\end{array}}^T\begin{array}{ccc}{\lambda }_1^{-1}& & \\ & {\lambda }_2^{-1}& \\ & & {\lambda }_3^{-1}\end{array}{\mathbf{M}}^T\begin{array}{c}\tilde{x}\\ \tilde{y}\\ \tilde{z}\end{array}=k^2,\end{array}$$

By substituting $$\begin{array}{}\text{(27)}& {\mathbf{M}}^T\begin{array}{c}\tilde{x}\\ \tilde{y}\\ \tilde{z}\end{array}=\begin{array}{c}u\\ v\\ w\end{array},\end{array}$$the following expression can be obtained: $$\begin{array}{}\text{(28)}& {\begin{array}{c}u\\ v\\ w\end{array}}^T\begin{array}{ccc}{\lambda }_1^{-1}& & \\ & {\lambda }_2^{-1}& \\ & & {\lambda }_3^{-1}\end{array}\begin{array}{c}u\\ v\\ w\end{array}=k^2.\end{array}$$

Additional calculation and expansion yield the following: $$\begin{array}{}\text{(29)}& \frac{u^2}{{\lambda }_1}+\frac{v^2}{{\lambda }_2}+\frac{w^2}{{\lambda }_3}=k^2.\end{array}$$

Equation (29) represents the error ellipsoid equation in the principal axis coordinate system. The lengths of the half-axis are $k\sqrt{{\lambda }_1}$, $k\sqrt{{\lambda }_2}$, and $k\sqrt{{\lambda }_3}$, where $k$ is the above-defined amplification factor and ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are the eigenvalues of the covariance matrix $D_{\text{rr}}$. The length of the corresponding ellipsoid axes can be determined by employing the covariance matrix and the value of $k$.

3.4. Determining Axial Directions of the Error Ellipsoid

It is known that the axial directions of the error ellipsoid can be described by the Euler angles which are rotation angles from the geodetic coordinate system $\text{OXYZ}$ to the ellipsoidal principal axis coordinate system $\text{OUVW}$. If the Euler angle rotates under the $Z-Y-X$ compliance, the rotation matrix is as follows: $$\begin{array}{}\text{(30)}& R={R_{\text{ij}}}_{3\times 3}=R_Z{\theta }_z{R}_{Y_1}{\theta }_y{R}_{X_2}{\theta }_x,\end{array}$$where ${\theta }_z$, ${\theta }_y$, and ${\theta }_x$ are the rotation angles around the $Z$-axis, the $Y$-axis after one rotation, and the $X$-axis after two sequentially rotations, and $R_Z{\theta }_z$, $R_{Y_1}{\theta }_y$ and $R_{X_2}{\theta }_x$ are corresponding rotation matrices [31], as shown in Figure 8.

[figure(s) omitted; refer to PDF]

According to Equation (23), the corresponding orthonormal eigenvectors of ${\mathbf{D}}_{\text{rr}}$ can be easily calculated and denoted by $\mathbf{M}$. Based on linear algebra [32], we can choose $\mathbf{M}$ such that ${\mathbf{M}}_{\text{ij}}$ is equal to ${\mathbf{R}}_{\text{ij}}$. So that the Euler angles ${\theta }_z$, ${\theta }_y$, and ${\theta }_x$ can be calculated.

3.5. Probability of Satellite Orbit Position within a Certain Error Range

According to the probability density shown in Equation (17) and the error ellipsoid shown in Equation (29), the probability of the satellite’s initial position in the error ellipsoid can be determined as follows: $$\begin{array}{}\text{(31)}& P={\iiint }_{B_k}\frac{1}{{2\pi }^{3/2}\sqrt{{\lambda }_1{\lambda }_2{\lambda }_3}}\exp -\frac{1}{2}\frac{u^2}{{\lambda }_1}+\frac{v^2}{{\lambda }_2}+\frac{w^2}{{\lambda }_3}\text{dudvdw}.\end{array}$$

Let $$\begin{array}{}\text{(32)}& \begin{array}{c}\tilde{u}=\frac{u}{\sqrt{2{\lambda }_1}},\hfill \\ \tilde{v}=\frac{v}{\sqrt{2{\lambda }_2}}\hfill \\ \tilde{w}=\frac{w}{\sqrt{2{\lambda }_3}},\hfill \end{array},\end{array}$$and substitute Equation (32) into Equation (31). Thus, the following expression can be obtained: $$\begin{array}{}\text{(33)}& P=\frac{1}{{\pi }^{3/2}}{\iiint }_{C_k}\exp -{\tilde{u}}^2-{\tilde{v}}^2-{\tilde{w}}^2\text{d}\tilde{\text{u}}\text{d}\tilde{\text{v}}\text{d}\tilde{\text{w}},\end{array}$$where $C_k$ is a sphere of radius $k/\surd 2$. Then, the probability that the satellite’s initial position is within the ellipsoid $B_k$ is equivalent to the probability of falling into the sphere $C_k$. The approximate probability that the satellite in orbit falls within a certain error range can be obtained as follows [33]: $$\begin{array}{}\text{(34)}& P=\frac{4}{\sqrt{2\pi }}\sum _{n=0}^{+\infty }{-1}^n·\frac{1}{n!}·\frac{1}{2n+3}·\frac{k^{2n+3}}{2^{n+1}},\end{array}$$

where $k$ is the ellipsoidal amplification factor mentioned above. If the covariance matrix is determined, the ratio between the three axes $\sqrt{{\lambda }_1}$, $\sqrt{{\lambda }_2}$, and $\sqrt{{\lambda }_3}$ of the error ellipsoid is determined. Since the amplification ratio of the three-axis lengths is different for various $k$ values, the probability of the satellite’s initial orbital position being within the error ellipsoid is also modified. The relationship between the amplification factor $k$ and probability $P$ of the satellite’s initial orbital position in the error ellipsoid is shown in Figure 9.

[figure(s) omitted; refer to PDF]

As shown in Figure 9(b), the satellite’s initial orbit position is within four times the ellipsoidal axis length under random error. The probability of a satellite’s position within different error limits can also be obtained.

3.6. Error Transfer of the Initial Position of the Satellite

According to the satellite’s orbit dynamics, only the eccentric anomaly of the satellite’s six orbital elements will change from time to time, affected by the semimajor axis, eccentricity, and initial eccentric anomaly. Combining Equations (1), (2), and (4) yields the following: $$\begin{array}{}\text{(35)}& E_1-e\sin E_1=\sqrt{\frac{\mu }{a^3}}·t_1+E_0-e\sin E_0,\end{array}$$where $E_1$ is the eccentric anomaly at the time $t_1$ and $E_0$ is the eccentric anomaly at the initial time $t_0$. According to the covariance propagation law, the covariance of the eccentric anomaly from the initial time $t_0$ to time $t_1$ can be expressed as follows: $$\begin{array}{}\text{(36)}& {\sigma }_{E_1}^2=\tilde{\mathbf{A}}\begin{array}{cccccc}{\sigma }_a^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_e^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_i^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{\Omega }^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{\omega }^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{E_0}^2\end{array}{\tilde{\mathbf{A}}}^T,\end{array}$$where $$\begin{array}{}\text{(37)}& \tilde{\mathbf{A}}={\begin{array}{cccccc}\frac{\partial E_1}{\partial a}& \frac{\partial E_1}{\partial e}& \frac{\partial E_1}{\partial i}& \frac{\partial E_1}{\partial \Omega }& \frac{\partial E_1}{\partial \omega }& \frac{\partial E_1}{\partial E_0}\end{array}}_{\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0}.\end{array}$$

Parameters $\partial E_1/\partial a$, $\partial E_1/\partial e$, $\partial E_1/\partial i$, $\partial E_1/\partial \Omega$, $\partial E_1/\partial \omega$, and $\partial E_1/\partial E_0$ are the partial derivatives of Equation (8) with respect to each variable. These values can be obtained by using approximate values $\tilde{a},\tilde{e},\tilde{i},\tilde{\Omega },\tilde{\omega },{\tilde{E}}_0$ instead of $a,e,i,\Omega ,\omega ,E_0$. Moreover, the accuracy is relatively high when are very close to their respective values of $a,e,i,\Omega ,\omega ,E_0$. The run time $t$ is discretized to $t_1,t_2,t_3,\cdots t_n$, and eccentric anomaly variance at time $t_n$ is recursively calculated from $t_0$ to obtain the error at any time $t$. Finally, the covariance matrix of six satellite orbital elements at time $t_n$ is obtained as follows: $$\begin{array}{}\text{(38)}& D_{\text{X}{\text{X}}_n}=\begin{array}{cccccc}{\sigma }_a^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_e^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_i^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{\Omega }^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{\omega }^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{E_n}^2\end{array}.\end{array}$$

Then, calculate the ellipsoid axes length, the Euler angles, and satellite position distribution probability according to Sections 3.3, 3.4, and 3.5, respectively.

4. Example Analysis and Simulation Verification

Three representative types of satellites are selected as an example: an LEO satellite with a circular orbit ($S-1$), an MEO satellite with an ellipsoid circle ($S-2$), and an IGSO satellite with a near-circular orbit ($S-3$). Random errors in the satellite’s expected initial orbit positions under the impact of random factors are simulated and verified if the error ellipsoid can accurately express random errors by comparing the Monte Carlo simulation results. Then, the transfer of errors is also simulated based on the error ellipsoid theory.

4.1. Related Satellites Expected Six Orbital Elements and Covariance

The initial six expected orbital elements of the three satellites are shown in Table 1.

Table 1

Satellite’s expected initial six orbital elements.

	$a$ (km)	$e$	$i$(°)	$\Omega$ (°)	$\omega$ (°)	$f$ (°)
$$S-1$$	6904.14	0	97.5	0	0	60
$$S-2$$	26553.4	0.740969	63.4	240.377	270.0	0
$$S-3$$	42167.2	0.002100	54.8	211.400	167.1	201.3
Standard deviation	2	0.0002	0.05	0.03	0.03	0.03

According to Equation (11), the covariance ${\mathbf{D}}_{\text{XX}}$ of six orbital elements can be expressed as follows: $$\begin{array}{}\text{(39)}& {\mathbf{D}}_{\text{XX}}=\begin{array}{cccccc}{2}^2& & & & & \\ & 0.0002^2& & & & \\ & & {\frac{0.05\pi }{180}}^2& & & \\ & & & {\frac{0.03\pi }{180}}^2& & \\ & & & & {\frac{0.03\pi }{180}}^2& \\ & & & & & {\frac{0.03\pi }{180}}^2\\ & & & & & \end{array}.\end{array}$$

The expected initial orbital positions in the geocentric inertial coordinate system can be obtained according to Equation (8), as shown in Table 2.

Table 2

The expected initial orbital positions.

	$x$ (km)	$y$ (km)	$z$ (km)
$$S-1$$	3452.1	-780.4	5928.0
$$S-2$$	-2677.2	1522.3	-6150.1
$$S-3$$	-33821.8	-24813.0	5043.4

4.2. The Random Initial Orbital Position Errors Based on Error Ellipsoid

The partial equation derivative with respect to the orbital six elements can be calculated according to Equations (8) and (20). The partial derivative matrix $\mathbf{A}$ is also acquired. Then, according to Equation (19), the uncertainty matrix ${\mathbf{D}}_{\text{rr}}$ under the geocentric inertial coordinate system can be obtained as follows: $$\begin{array}{}\text{(40)}& \begin{array}{c}{\mathbf{D}}_{\text{r}{\text{r}}_{S-1}}=\begin{array}{ccc}20.8885& 1.9628& -9.2987\\ 1.9628& 30.1972& 2.2432\\ -9.2987& 2.2432& 10.1869\end{array},\hfill \\ {\mathbf{D}}_{\text{r}{\text{r}}_{\text{S}-2}}=\begin{array}{ccc}33.0537& -2.5663& -2.6298\\ -2.5663& 29.9996& 1.4953\\ -2.6298& 1.4953& 29.9866\end{array},\hfill \\ {\mathbf{D}}_{\text{r}{\text{r}}_{S-3}}=\begin{array}{ccc}390.5257& -377.3623& 320.5532\\ -377.3623& 515.2274& -317.6344\\ 320.5532& -317.6344& 650.0961\end{array}.\hfill \end{array}\end{array}$$

The error ellipsoid functions of the satellite’s initial orbital position error can be obtained by substituting ${\mathbf{D}}_{\text{rr}}$, $r_x$, $r_y$, and $r_z$ into Equation (17). The orthogonal matrix $\mathbf{M}$ is then obtained as follows: $$\begin{array}{}\text{(41)}& \begin{array}{c}{\mathbf{M}}_{S-1}=\begin{array}{ccc}-0.5000& 0.8493& 0.1691\\ 0.1130& -0.1296& 0.9851\\ -0.8586& -0.5117& 0.0312\end{array},\hfill \\ {\mathbf{M}}_{S-2}=\begin{array}{ccc}0.3892& -0.4943& 0.7773\\ -0.2213& -0.8693& -0.4420\\ 0.8942& 0& -0.4478\end{array}\hfill \\ {\mathbf{M}}_{S-3}=\begin{array}{ccc}0.8005& -0.3021& -0.5177\\ 0.5872& 0.5681& 0.5766\\ -0.1199& 0.7655& -0.6322\end{array}.\hfill \end{array},\end{array}$$

The corresponding diagonal matrix $\lambda$ (km) can be obtained via Equation (23). Then, $\lambda$ can be calculated as follows: $$\begin{array}{}\text{(42)}& \begin{array}{c}{\lambda }_{S-1}=\begin{array}{ccc}4.4767& & \\ & 26.1907& \\ & & 30.6052\end{array},\hfill \\ {\lambda }_{S-2}=\begin{array}{ccc}28.4717& & \\ & 28.5404& \\ & & 36.0279\end{array},\hfill \\ {\lambda }_{S-3}=\begin{array}{ccc}65.6803& & \\ & 287.8779& \\ & & 1202.2910\end{array}.\hfill \end{array}\end{array}$$

Lastly, the error vector of the satellite’s initial orbital position can be calculated. The calculation results for the amplification factor $k$ equal to one are shown in Table 3. According to Equation (34), when $k$ is 2.8, the satellite’s initial orbital positions will be in the ellipsoid range of approximately 95%. Error ellipsoids are shown in Figure 10.

Table 3

Results based on error ellipsoid analysis.

Parameters	$S-1$	$S-2$	$S-3$
The ellipsoid center (km)	(3452.1, -780.4, 5928.0)	(-2677.2, 1522.3, -6150.1)	(-33822, -24813, 5043)
Minor axis length (km)	2.1158	5.3359	8.1043
Middle axis length (km)	5.1177	5.3423	16.9670
Major axis length (km)	5.5322	6.0023	34.6741
$${\theta }_x,{\theta }_y,{\theta }_z{}^{\circ }$$	(342.7417, 277.8334, 253.1135)	(72.0364, 27.4769, 79.0675)	(327.5203, 306.5469, 283.5594)

[figure(s) omitted; refer to PDF]

Table 3 and Figure 10 demonstrate that the error ellipsoids of the three satellites differ in size, shape, and direction because they are mainly determined by six orbital elements and initial covariance matrices. Although these three types of satellites have the same covariance, the values of the six orbit elements are different from each other. So, the error ellipsoids for $S-1$, $S-2$, and $S-3$ are different in size, shape, and direction. Moreover, the changes of the error ellipsoid at every position are very complex, and we can easily get the random error of satellites by the error ellipsoid model.

4.3. Random Initial Orbital Position Based on the Monte Carlo Simulation

According to the following equation [34], the number of samples can be determined. $$\begin{array}{}\text{(43)}& N\ge {\frac{Z_{\alpha }·\sigma }{\epsilon }}^2,\end{array}$$

where $N$ is the number of samples, $\sigma$ is the standard deviation of the random variable, $\epsilon$ is the simulation error of the Monte Carlo method, and $Z_{\alpha }$ is the $\alpha$ quantile on the standard normal distribution. Let $\sigma =0.1$, $Z_{\alpha }=3$, and $\epsilon =0.002$, then $N\ge 22500$. In this experiment, $25000$ samples are chosen. The results of the Monte Carlo simulation using the same six orbital elements and the same covariance matrix are shown in Figure 11.

[figure(s) omitted; refer to PDF]

Figure 11 shows the initial positions of the satellites based on the Monte Carlo method, represented by black points, and the respective ellipsoids derived from Figure 12. The figure illustrates that the initial position distribution areas exhibit ellipsoidal shapes on the three planes of the coordinate system. Additionally, the probability is higher in the central area compared to the edge area. These features are consistent with the error ellipsoid model proposed in this paper.

[figure(s) omitted; refer to PDF]

4.4. Comparison with the Error Ellipsoids and the Simulation Results Based on the Monte Carlo Method

According to Figures 10 and 11, the error ellipsoids are very similar to the simulation results of the Monte Carlo method, including the shape, size, and rotation angles. To further verify the rationality of the error ellipsoid for expressing random errors, the initial positions of the satellite simulation are compared with the range of the error ellipsoid, and the probability of the satellite’s orbital initial position in the ellipsoid is calculated with a certain $k$, as shown in Table 4. The probability of the satellite’s initial orbital position estimated by the error ellipsoid model is compared with the simulation results for the entire range of the amplification factor $k$ from 0 to 4, as shown in Figure 12. At $k=0$, the cumulative probability distribution is consistently zero for both the error ellipsoid model and the Monte Carlo method. For certain values of $k$, the two methods have slight errors. However, as the size of the ellipsoid increases with $k$, the probabilities of the satellite’s initial orbital position within the ellipsoid increase and eventually reach 1 for both methods, resulting in a greater degree of similarity in the cumulative probability distribution.

Table 4

Probability of the satellite’s initial position distribution in a certain $k$.

Value of $k$	Distribution probability based on error ellipsoid theory (%)	The results of the Monte Carlo simulation probability (%)
		$S-1$	$S-2$	$S-3$
1.00	19.87	22.64	22.55	21.89
1.50	47.78	50.95	50.21	50.48
2.00	73.85	76.25	75.57	75.50
2.50	89.99	90.98	90.90	90.99
2.80	95.06	95.86	95.78	95.44

According to Table 4 and Figure 12, the data from the two methods about three types of satellites are all very close at a certain $k$ value. The probability of the Monte Carlo simulation is slightly bigger than those of the error ellipsoid results. Since the Monte Carlo experiments were conducted many times, the results are close to the actual result. Therefore, it can be considered that the error ellipsoid is close to the actual one, and the distribution probability is slightly lower than the actual probability. This indicates that the error ellipsoid theoretical model can reliably estimate the satellite’s orbital position under random error.

Figure 13 shows the calculation time spent by the error ellipsoid theoretical model and the Monte Carlo method.

[figure(s) omitted; refer to PDF]

According to the figure, it can be found that the calculation time of the Monte Carlo method increases linearly with the improvement of simulation accuracy, while the proposed method in this paper remains almost unchanged. When the simulation error control is 0.2 %, the Monte Carlo method calculation takes 1000 times longer than the error ellipsoid model method. Therefore, the method proposed in this article has higher computational efficiency.

4.5. Promotion and Application of the Error Ellipsoid Theory

Moreover, the error ellipsoid theory can analyze a satellite’s initial position error range under any orbital element and covariance. Figure 14 shows the error ellipsoid axial length and direction of the satellite’s initial position at different true anomalies for the LEO, MEO, and IGSO satellites. In fact, the Euler angles and the lengths of the axes are all gradually changing with the true anomaly. Specifically, in order to well display the directions of the error ellipsoid axes, the ranges of the Euler angles are specified as follows: ${\theta }_x\in 0^{\circ },360^{\circ }$, ${\theta }_y\in -180^{\circ },180^{\circ }$, and ${\theta }_z\in 0^{\circ },360^{\circ }$. When the Euler angle exceeds the value range, it will be equivalently converted to the specified value range, and it will cause a dramatic change in the figure.

[figure(s) omitted; refer to PDF]

According to Figures 14 and 15, different change rules accompany various orbits. If one of the six orbital elements changes, the size, shape, and angles of the axes’ inclination of the error ellipsoid will change. These changes are very complex and difficult to regulate. However, they are all symmetric about the major or minor axes of the orbit.

[figure(s) omitted; refer to PDF]

According to the calculation and the conducted analyses, the accuracy of the error propagation with time will drop rapidly using the error ellipsoid theory. However, it can roughly estimate the error range of the satellite’s position in orbit. The length changes of the three axes for the LEO satellite are shown in Figure 16.

[figure(s) omitted; refer to PDF]

According to Figure 16, the minor axis changes periodically. The maximum value is approximately 2.43 km and appears at the true anomaly of 0° and 180°, while the minimum value is roughly 2 km and appears at the true anomaly of 90° and 270°. Similar change regulation is observed for the middle and major axes. During a specific period, the maximum values appear at the true anomaly of 90° and 270°, and the minimum values appear at 0° and 180°. Furthermore, the major axis will become longer over time, while different orbital elements will lead to various change laws, complicating the entire process.

The ellipsoidal volume is introduced to explain the variation of satellite orbit position error range. The size of the error range response of an ellipsoidal volume is shown indirectly in Figure 17.

[figure(s) omitted; refer to PDF]

As shown in Figure 17, the volume of the error ellipsoid increases gradually. In each period, the smallest point is found at 0° and 180° of the true anomaly, and the largest position is found at 90° and 270° of the true anomaly.

According to calculation and analysis, the accuracy will decrease rapidly with time. Therefore, the error ellipsoid theoretical model can only be used to estimate the position error roughly.

5. Conclusions

In this paper, the rationality of the error ellipsoid for describing the positioning error of satellites in orbit was deduced and demonstrated. According to the uncertainty matrix of the satellite’s six orbital elements under random error, the equal probability density surface of the random error of the initial orbital position was approximated as an ellipsoid. The length and the direction of the three axes of the ellipsoid were determined, and the calculation method of the probability of the satellite entering the orbit within a certain error range was provided. According to the actual case of aerospace engineering, the experiments of launching LEO, MEO, and IGSO satellites into orbit under the influence of random factors were simulated using the Monte Carlo method and STK software. Thus, the actual distribution of the experimental satellite’s initial orbit positions was obtained. A comparison with the probability distribution of the orbital position under the error ellipsoid model showed that the analysis results of the ellipsoid model were consistent with the simulation results. Consequently, it can be concluded that the error ellipsoid theory can be used to estimate the random orbital position error. Lastly, the satellite’s initial position error propagation in orbit can be simply simulated using the error ellipsoid model.

Accurately estimating the random error of the satellite in orbit is of great significance, and it can calculate the collision safety during satellites and the coverage effect on the ground. The contribution of this paper consists in the rapid, simple, and reliable model of error ellipsoid to determine the satellites’ orbital positions. This model is very useful for large-scale constellations. This study is based on six independent variances of orbital elements, and an approximate treatment is made during error transmission. However, the variances are complex, and the accuracy of error transmission will decrease significantly with time. In the next stage, the error of the satellite’s orbit and its impact on the efficiency can be further investigated via the ellipsoid theory and the error propagation characteristics.

Acknowledgments

This work is supported by the Natural Science Foundation of China under Grant nos. 42271391 and 62006214, the Joint Funds of Equipment Pre-Research and Ministry of Education of China under Grant no. 8091B022148, the 14th Five-year Pre-research Project of Civil Aerospace in China, and the Hubei excellent young and middle-aged science and technology innovation team plan project under Grant no. T2021031.

References

[1] F. R. Hoots, P. W. Schumacher, R. A. Glover, "History of analytical orbit modeling in the U.S. space surveillance system," Journal of Guidance, Control and Dynamic, vol. 27 no. 2, pp. 174-185, DOI: 10.2514/1.9161, 2004.

[2] J. S. Knogl, P. Henkel, C. Gunter, "Precise positioning of a geostationary data relay using LEO satellites," 53rd International Symposium ELMAR-2011, pp. 325-328, .

[3] N. V. Vighnesam, A. Sonney, P. K. Soni, "CARTOSAT-1 orbit determination system and achieved accuracy during early phase," AIAA/AAS Astrodynamics Specialist Conference, pp. 511-522, .

[4] B. Yi, D. Gu, K. Shao, B. Ju, H. Zhang, X. Qin, X. Duan, Z. Huang, "Precise relative orbit determination for Chinese TH-2 satellite formation using onboard GPS and BDS2 observations," Remote Sensing, vol. 13 no. 21,DOI: 10.3390/rs13214487, 2021.

[5] M. Liu, Y. Yuan, J. Ou, G. Yang, "Precise orbit determination and precision comparison for FY3C and FY3D with receiver antenna GPS and BDS PCV using spaceborne BDS and GPS observation data," Advance in Space Research, vol. 71 no. 1, pp. 375-389, DOI: 10.1016/j.asr.2022.08.070, 2023.

[6] D. Li, X. Cheng, X. Jia, W. Yang, "Precision analysis of BDS-3 satellite orbit by using SLR data," China Satellite Navigation Conference (CSNC) 2019 Proceedings. CSNC 2019, pp. 389-399, DOI: 10.1007/978-981-13-7759-4_35, 2019.

[7] D. Gu, X. Tu, D. Yi, "System error calibration for GPS precise orbit determination with SLR data," Journal of National University of Defense Technology, vol. 30 no. 6, pp. 14-18, 2008.

[8] Z. Li, X. Yang, G. Ai, H. Si, R. Qiao, C. Feng, "A new method for determination of satellite orbits by transfer," Science in China Series G: Physics, Mechanics & Astronomy, vol. 52 no. 3, pp. 384-392, 2009.

[9] R. Deng, H. Qin, H. Li, D. Wang, H. Lyu, "Non-cooperative LEO satellite orbit determination based on single pass Doppler measurements," IEEE Transactions on Aerospace and Electronic Systems, vol. 59,DOI: 10.1109/TAES.2022.3194977, 2022.

[10] M. H. A. Larki, M. Maboodi, H. Bolandi, "Satellite orbit determination based on gradient method," 2012 24th Chinese Control and Decision Conference (CCDC), pp. 2554-2559, DOI: 10.1109/ccdc.2012.6244407, .

[11] X. Bai, Research on Orbital Prediction Error and Collision Probability of Space Objects, [Ph.D. Thesis], 2013.

[12] N. Cressie, J. Kornak, "Spatial statistics in the presence of location error with an application to remote sensing of the environment," Statistical Science, vol. 18 no. 4, pp. 436-456, DOI: 10.1214/ss/1081443228, 2003.

[13] J. Gabrosek, N. Cressie, "The effect on attribute prediction of location uncertainty in spatial data," Geographical Analysis, vol. 34 no. 3, pp. 262-285, DOI: 10.1111/j.1538-4632.2002.tb01088.x, 2002.

[14] G. Arbia, D. A. Griffith, R. P. Haining, "Spatial error propagation when computing linear combinations of spectral bands," Environmental and Ecological Statistics, vol. 10 no. 3, pp. 375-396, DOI: 10.1023/A:1025167225797, 2003.

[15] C. R. Ehlschlaeger, "Representing multiple spatial statistics in generalized elevation uncertainty models£°moving beyond the variogram," International Journal of Geographical Information Science, vol. 16 no. 3, pp. 259-285, DOI: 10.1080/13658810110099116, 2002.

[16] F. J. Aguilar, F. Carvajal, "Effects of terrain morphology£¬sampling density£¬and interpolation methods on grid DEM accuracy," Photogrammetric Engineering and Remote Sensing, vol. 71 no. 7, pp. 805-816, DOI: 10.14358/PERS.71.7.805, 2005.

[17] F. J. Aguilar, M. A. Aguilar, F. Agüera, "Accuracy assessment of digital elevation models using a non-parametric approach," International Journal of Geographical Information Science, vol. 21 no. 6, pp. 667-686, DOI: 10.1080/13658810601079783, 2007.

[18] F. J. Aguilar, J. P. Mills, "Accuracy assessment of LIDAR-DERIVED digital elevation models," The Photogrammetric Record, vol. 23 no. 122, pp. 148-169, DOI: 10.1111/j.1477-9730.2008.00476.x, 2008.

[19] Z. Li, "Mathematical models of the accuracy of digital terrain model surfaces linearly constructed from square gridded data," Photogrammetric Record, vol. 14 no. 82, pp. 661-674, DOI: 10.1111/j.1477-9730.1993.tb00776.x, 1993.

[20] Z. Li, "A comparative study of the accuracy of digital terrain models (DTMs) based on various data models," ISPRS Journal of Photogrammetry and Remote Sensing, vol. 49 no. 1,DOI: 10.1016/0924-2716(94)90051-5, 1994.

[21] G. J. Hunter, M. F. Goodchild, "Modeling the uncertainty of slope and aspect estimates derived from spatial databases," Geographical Analysis, vol. 29 no. 1, pp. 35-49, DOI: 10.1111/j.1538-4632.1997.tb00944.x, 1997.

[22] P. C. Kyriakidis, A. M. Shortridge, M. F. Goodchild, "Geostatistics for conflation and accuracy assessment of digital elevation models," International Journal of Geographical Information Science, vol. 13 no. 7, pp. 677-707, DOI: 10.1080/136588199241067, 1999.

[23] C. Q. Zhu, W. Z. Shi, Q. Li, G. Wang, T. C. K. Cheung, E. Dai, G. Y. K. Shea, "Estimation of average DEM accuracy under linear interpolation considering random error at the nodes of TIN model," International Journal of Remote Sensing, vol. 26 no. 24, pp. 5509-5523, DOI: 10.1080/10245330500169029, 2005.

[24] C. Lopez, "Locating some types of random errors in digital terrain models," International Journal of Geographical Information Science, vol. 11 no. 7, pp. 677-698, DOI: 10.1080/136588197242149, 1997.

[25] M. Albani, B. Klinkenberg, "A spatial filter for the removal of striping artifacts in digital elevation models," Photogrammetric Engineering and Remote Sensing, vol. 69 no. 7, pp. 755-765, DOI: 10.14358/PERS.69.7.755, 2003.

[26] J. Zhao, Orbital Dynamics of Spacecraft, 2011.

[27] S. Wang, Error Theory and Surveying Adjustment, 2015.

[28] L. Werner, Probability Theory, 2016.

[29] G. Dai, X. Chen, M. Zuo, L. Peng, M. Wang, Z. Song, "The influence of orbital element error on satellite coverage calculation," International Journal of Aerospace Engineering, vol. 2018,DOI: 10.1155/2018/7547128, 2018.

[30] Q. Xie, Advanced Algebra, 2022.

[31] H. Li, E. Bai, M. Mi, Y. Yan, "Identification of Euler angles of permanent magnet spherical motor rotor based on hall sensors array," Measurement, vol. 199, article 111500,DOI: 10.1016/j.measurement.2022.111500, 2022.

[32] Z. Chen, Higher Algebra and Analytic Geometry, 2008.

[33] H. Vazquez-Leal, R. Castaneda-Sheissa, U. Filobello-Nino, A. Sarmiento-Reyes, J. Sanchez Orea, "High accurate simple approximation of normal distribution integral," Mathematical Problems in Engineering, vol. 2012,DOI: 10.1155/2012/124029, 2012.

[34] T. Zou, L. Zhao, "A method for estimating sample size of Monte Carlo method in accident reconstruction," China Safety Science Journal, vol. 23 no. 5, pp. 22-26, 2013.