Random Error Estimation and Propagation Analysis for Satellites ’ Initial Positions

. 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 a ﬀ ects space safety, network accuracy, and operating e ﬃ ciency. 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 e ﬀ ective, and good at accuracy.


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 (Bei-Dou 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][16][17][18][19][20], data error propagation model [21][22][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.

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

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Õ, 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 Ω. The field angle between the perigee and the ascending node with respect to the geocenter is called the perigee argument and is marked as ω. 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.

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. OXYZ is the geocentric inertial coordinate system, as shown in Figure 1.
In the geocentric orbital coordinate system, the coordinate origin point is located at geocentric O, while the OU and OV axes are in the orbit plane of the satellite. Simultaneously, the OU axis points to the perigee, and the OW axis coincides with the normal vector of the satellite orbital 2 International Journal of Aerospace Engineering plane. The axes OV, OU, and 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Õ and the radius a, as shown in Figure 2. Draw a vertical line of the OU axis through S and intersect with OU at point H. Then, extend HS intersection auxiliary circle toS, connect pointÕ and pointS, and the included angle between lineÕS and lineÕU is the eccentric anomaly E. The angle between the lineÕS andÕU is the true anomaly, marked as f . By Equation (1), point C can be obtained. The included angle between the lineÕU andÕ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: 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: where μ is the Kepler constant, whose value is 3:986 × 10 5 km 3 /s 2 .

Equation of Satellite
Motion. The satellite's orbital coordinate is usually expressed by six orbital elements, written as the vector: 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]: The three-dimensional coordinate expression of a satellite position in the geocentric orbital coordinate system is the following: where U is the component of a satellite orbital position in the direction of the OU axis, V is the component of a satellite's position in the direction of the OV axis, and W is the component in the direction of the OW axis. The geocentric inertial coordinate system can be written as follows: x y z 2 6 6 4 where R 1 , R 2 , and R 3 are rotation matrices expressed as follows: The elements on the main diagonal in D XX are variances of random variables a, e, i, Ω, ω, 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 D XX can be written as follows: : ð11Þ Since they are subjected to the random error of the orbital position of the satellite, the vector μ X and the variance matrix D XX are the digital characteristics of the 6dimensional 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.

Random Error Expression of the Initial Orbital Position
in the Geocentric Inertial Coordinate System. Suppose that the expected satellite's initial orbital position isÃðx,ỹ,zÞ in the geocentric inertial coordinate system. However, the actual position is Aðx, y, zÞ, as shown in Figure 3.
The deviation value Δ of the actual initial position can be expressed as follows: where Δ x , Δ y , and Δ z are the components of the spatial position error Δ on the X, Y, and Z axes, respectively.
sin E cos ω sin i According to the definition of mean and variance, According to the linear principle of the mean value and by combining Equations (13) and (14), where EðΔ 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 σ 2 A , it can be expressed as follows: where the magnitude of σ 2 A is independent of the coordinate system and σ A is the satellite's orbital position error [27].

Construction of Satellite's Orbital Position Error Model
The satellite's initial position variance σ 2 A 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.

The Function for the Error Ellipsoid
Model. The satellite's initial position in orbit is presented as a normal threedimensional distribution in the geocentric inertial coordinate system. Moreover, its density function of the joint distribution can be expressed as follows [28]: wherex,ỹ, andz are the mathematical expectations of the initial position on the X, Y, and Z axes. Respectively, D rr is the uncertainty matrix of the satellite's initial position in orbit under the geocentric inertial coordinate system, and jD rr j is the determinant of the uncertainty matrix. The matrix D rr can be expressed as follows: where σ 2 i ði = x, y, zÞ is the position variance of the coordinate axis direction and σ ij ði, j = x, y, zÞ is the covariance between the coordinate axes.

Calculation of Uncertainty
Matrix D rr . The covariance matrix D XX of six orbital elements is expressed as Equation (11). Therefore, the uncertainty matrix D rr can be calculated according to the covariance propagation rate [29]: A (x, y, z) Δ Figure 3: Error in the satellite's initial position.  where each element of the matrix 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 μ X , as shown in Equation (9). Matrix A can be expressed as follows: Elements in the matrix A are partial derivatives of the function with respect to each variable. These values are obtained using the approximate valueã,ẽ,ĩ,Ω,ω,Ẽ 0 in place of a, e, i, Ω, ω, E 0 , whereã,ẽ,ĩ,Ω,ω,Ẽ 0 are constants. The accuracy is relatively high when the values a,ẽ,ĩ,Ω,ω,Ẽ 0 are very close to their respective valuesã,ẽ,ĩ, Ω,ω,Ẽ 0 .

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: 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.
Then, Equation (21) can be written as follows: 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 OUVW, as shown in Figure 7.
In Figure 7, OU, OV, and OW are the principal axes of the ellipsoid family. According to Equations (19) and (20), D rr is a real symmetric matrix, and there must be an orthogonal matrix M that satisfies the following [30]: Equation (23) can be transformed, and an inverse is taken on both sides to obtain the following: Therefore, The following expression is obtained by combining Equations (21) and (24).
Additional calculation and expansion yield the following: Equation (29) represents the error ellipsoid equation in the principal axis coordinate system. The lengths of the half-axis are k ffiffiffiffi ffi , where k is the above-defined amplification factor and λ 1 , λ 2 , and λ 3 are the eigenvalues of the covariance matrix D rr . The length of the corresponding ellipsoid axes can be determined by employing the covariance matrix and the value of k.

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 OXYZ to the ellipsoidal principal axis coordinate system OUVW. If the Euler angle rotates under the Z − Y − X compliance, the rotation matrix is as follows: where θ z , θ y , and θ 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 ðθ z Þ, R Y 1 ðθ y Þ and R X 2 ðθ x Þ are corresponding rotation matrices [31], as shown in Figure 8. According to Equation (23), the corresponding orthonormal eigenvectors of D rr can be easily calculated and denoted by M. Based on linear algebra [32], we can choose M such that M ij is equal to R ij . So that the Euler angles θ z , θ y , and θ x can be calculated.

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: and substitute Equation (32) into Equation (31). Thus, the following expression can be obtained: where C k is a sphere of radius k/√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]: where k is the ellipsoidal amplification factor mentioned above. If the covariance matrix is determined, the ratio between the three axes ffiffiffiffi ffi λ 1 p , ffiffiffiffi ffi λ 2 p , and ffiffiffiffi ffi λ 3 p 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.
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.

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,  International Journal of Aerospace Engineering eccentricity, and initial eccentric anomaly. Combining Equations (1), (2), and (4) yields the following: 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: Parameters ∂E 1 /∂a, ∂E 1 /∂e, ∂E 1 /∂i, ∂E 1 /∂Ω, ∂E 1 /∂ω, and ∂E 1 /∂E 0 are the partial derivatives of Equation (8) with respect to each variable. These values can be obtained by using approximate valuesã,ẽ,ĩ,Ω,ω,Ẽ 0 instead of a, e, i, Ω, ω, E 0 . Moreover, the accuracy is relatively high when are very close to their respective values of a, e, i, Ω, ω, E 0 . The run time t is discretized to t 1 , t 2 , t 3 , ⋯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: 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.

Example Analysis and Simulation Verification
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      Error ellipsoids are shown in Figure 10. 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.

Random Initial Orbital Position Based on the Monte
Carlo Simulation. According to the following equation [34], the number of samples can be determined.
where N is the number of samples, σ is the standard deviation of the random variable, ε is the simulation error of the Monte Carlo method, and Z α is the α quantile on the standard normal distribution. Let σ = 0:1, Z α = 3, and ε = 0:002, then N ≥ 22500. In this experiment, 25000

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. 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.
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.

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: θ x ∈ ½0 ∘ , 360 ∘ Þ, θ y ∈ ½−180 ∘ , 180 ∘ Þ, and θ z ∈ ½0 ∘ , 360 ∘ Þ. 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.    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.
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.  14 International Journal of Aerospace Engineering 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

15
International Journal of Aerospace Engineering error range response of an ellipsoidal volume is shown indirectly in Figure 17.
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.

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.

Data Availability
Data and materials for replication of the case studies will be available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.

16
International Journal of Aerospace Engineering