A Semi-Analytical and Monte Carlo-Based Phase Dynamic Evolution Approach for LEO Mega-Constellations

Abstract

In recent years, with the reduction of the cost of microsatellites, the development of commercial rockets and the multi-satellite launching technology, the construction of large-scale constellations in low-Earth orbit (Mega-Constellations) has become a development trend. Since the motion of LEO satellites is affected by perturbations such as non-spherical gravitational fields and atmospheric drag, as well as the uncertainty of actuators, measurement systems, and dynamic models, it is easy to cause divergence of constellation configurations. The station-keeping control of the satellites is crucial for the stable operation of the mega-constellation. Aiming at this problem, this paper proposes an uncertainty propagation approach based on semi-analytical and Monte Carlo for LEO Mega-Constellations. Under the assumption that initial uncertainty on the osculating trajectory is Gaussian distribution, through hypothesis testing analysis, the uncertainty propagation simulations of a single satellite suggest that the satellite argument of latitude and the relative phase of co-plane satellites can be both considered as Gaussian distributions with zero means. Multi-group Monte Carlo simulations with product-based least-squares surface fitting establish an approximate mapping between initial and terminal errors. The mapping provides an efficient method for deviation prediction and can be used to design the station-keeping control strategy.

1. Introduction

The development of microsatellites leads to a general trend in Networking with LEO mega-constellations. Many countries and companies have proposed building this type of constellation, such as OneWeb and Starlink in the United States [1]. Of these, the Starlink communications constellation of 12,000 satellites is expected to be constructed around 2025, and accumulated to 42,000 satellites afterward [2].

Compared with the traditional middle and high Earth orbit navigation constellations, the LEO mega-constellations are characterized by low orbit height, dense satellite distribution, large perturbation, and high requirements for station position maintenance. Although the relative phase between satellites will not drift in constellations with the same value of mean elements, due to errors in the measurement data, dynamic model, solution algorithm, actuator, and other factors, the constellation configuration will be out of order in the uncontrolled state, leading to a high risk of collisions within the constellations.

To investigate how much the relative phase between satellites in the same orbital plane in the LEO mega-constellations will be drifted by the perturbations, it is critical to first point out the nominal value of the constellation design. From the semi-analytic theory [3], the satellite orbits can be described by their initial mean elements. If the initial mean elements of the constellation are identical, the relative phases within the constellation will not theoretically diverge. So, the constellation nominal values can be considered a set of mean elements. Accordingly, the problem of phase evolution affected by perturbations can be viewed as the initial uncertainties propagation of the constellation: how to establish the mapping relationship of the terminal state to the initial state. The initial state is the deviation of the position velocity of the two satellites from the nominal design value, and the terminal state is the deviation of the relative phase of the two satellites from the design nominal value after a while.

As widely recognized, orbit propagation is a fundamental problem in astrodynamics. There are three main approaches for modeling orbit propagation [4]: the numerical approach, analytical approach, and semi-analytical approach. The numerical approach [5,6] has high precision and is suitable for scenarios with complex mechanical environments, but consumes more computational resources. The analytical approach [7] is developed on Brouwer’s perturbation solutions [8,9], using an approximate analytical solution instead of the original ordinary differential equations of motion, preserving the main features of the perturbation while neglecting higher-order small terms. The analytical approach makes the calculation much more efficient while losing some accuracy. The semi-analytical approach [10,11] mixes the analytical and numerical approaches’ advantages, classification of perturbation solutions into long-term, long-period, and short-period terms according to perturbation characteristics. Bezděk [12] presents the computational efficiency strength of the semi-analytical approach in long-term dynamics. Golikov [13] presents a dynamic approach of satellite formation flights based on the THEONA semi-analytical theory with arbitrary values of the eccentricity.

For general stochastic dynamical systems, the evolution equation of the probability density over time satisfies the Fokker–Planck equation [14] (FPE), while for orbital dynamical systems with high dimensionality (6 dimensions) and nonlinear perturbations, the FPE solution becomes very tough. Thus, the probability density function (PDF) evolution often requires the assumption of local linearization and only the lower order moments (e.g., mean and covariance matrices) are analyzed. Currently, nonlinear uncertainty propagation approaches include polynomial chaos expansions [15], state transition tensor approaches [16,17], differential algebra approaches [18,19], and Gaussian mixture model approaches [20]. Most of the theories of uncertainty propagation are, however, focused on explaining the behavior of medium and high Earth orbit navigation constellations, seldom involving LEO mega-constellations. A well-known and robust approach used to study the uncertainty propagation problem is the Monte Carlo simulation [21], and the distribution of terminal states is obtained through the steps of randomly taking a large number of samples, high-precision simulation, and statistical analysis. The numerical approach-based Monte Carlo simulation has a long computation time and high hardware requirements. A possible solution to the problem at hand is applying the semi-analytical dynamics model, which does not need to propagate the orbital elements of the satellite or the position and velocity vectors of the satellite to obtain the position information of the satellite at a specific time and therefore has a smaller computational memory and computational load [22].

This paper aims to develop an overarching framework to reveal the uncertainty propagation law of the phase of LEO mega-constellations. The perturbation solution is modeled in Section 2. Section 3 describes the Monte Carlo uncertainty propagation approach. Section 4 carries out a numerical simulation. Conclusions are drawn in Section 5.

2. Modeling the Orbital Dynamics

In the geocentric equatorial inertial coordinate systems, the orbit motion of a satellite can be reduced to a perturbed two-body problem, and the motion equation is a set of ordinary differential equations [23]

$$\left\{\begin{array}{c}\ddot{r}=F_0\left(r\right)+F_{\epsilon }(r,\dot{r},t;\epsilon ),\hfill \\ r\left(t_0\right)=r_0,\dot{r}\left(t_0\right)={\dot{r}}_0,\hfill \end{array}\right.$$

(1)

where r and $\dot{r}$ are the position vectors and velocity vectors of the satellite, t is the propagation time, $F_0$ is the two-body motion acceleration of the central body, $F_{\epsilon }$ is the perturbation acceleration, and $\epsilon$ is a small parameter. Utilizing asymptotic expansion and other approaches, the power series form solution of small parameter $\epsilon$ can be constructed [24], such as

$$\kappa \left(t\right)={\kappa }^{\left(0\right)}\left(t\right)+\epsilon \beta \left(t\right)+{\epsilon }^2{\beta }^2\left(t\right)+\cdots +{\epsilon }^l{\beta }^l\left(t\right),$$

(2)

where $\kappa ={(a,e,i,\Omega ,\omega ,M)}^{\mathrm{T}}$ are the classical orbital elements: semimajor axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and mean anomaly at epoch, respectively. By substituting the formal solution into the original perturbed motion equation and comparing the coefficients of the same power ${\epsilon }^l$ at both ends, the perturbation solutions of each order can be obtained. By the average approach, the perturbation solutions can be divided into secular terms, long-period terms, and short-period terms according to forms of presentation.

$$\left\{\begin{array}{cc}\hfill \kappa \left(t\right)& =f({\overline{\kappa }}_0,t)=\overline{\kappa }\left(t\right)+{\kappa }_{\mathrm{lp}}\left(t\right)+{\kappa }_{\mathrm{sp}}\left(t\right),\hfill \\ \hfill \overline{\kappa }\left(t\right)& ={\overline{\kappa }}_0+(\delta n+{\kappa }_{\sec }^{\left(1\right)}+{\kappa }_{\sec }^{\left(2\right)}+\cdots )(t-t_0).\hfill \end{array}\right.$$

(3)

Among them, $t_0$ is the epoch time, $\kappa \left(t\right)$ are the osculating elements at time $\kappa \left(t\right)$, ${\overline{\kappa }}_0$ is the initial mean elements, $\overline{\kappa }\left(t\right)$ is the mean elements at time t, $\delta ={[0,0,0,0,0,1]}^T$, n is the mean motion of satellite in two body, ${\kappa }_{\sec }(t-t_0)$ makes up the secular terms, ${\kappa }_{\mathrm{lp}}\left(t\right)$ consists of the long period terms, and ${\kappa }_{\mathrm{sp}}\left(t\right)$ consists of the short-period terms. Besides this, Equation (3) describes the identity mapping between mean elements and the osculating elements, which can be solved by numerical algorithm such as Newton iteration [25].

LEO constellations are mainly affected by atmospheric drag and non-spherical Earth gravity. The first-order and second-order secular term of $J_2$, the first-order long-period term of $J_2$, the first-order short-period term of $J_2$, and the secular term of atmospheric drag are involved in the semi-analytical dynamic model. It should be noted that in order to avoid the abnormal matrix caused by the large numerical difference among the Keplerian elements, distance and time variables need to be normalized by $\left[R_e\right]$ and $\left[{(R{e}^3/\mu )}^{1/2}\right]$, where $R_e$ is the Earth’s equatorial radius and $\mu$ is the Earth’s gravitational constant. Thus, the mean motion n is expressed by $a^{-3/2}$. Detailed semi-analytical solutions are recorded in [26].

2.1. Earth Non-Spherical Gravity Perturbation Solution

The bulge at equator exerts a pull on the satellite toward the equatorial plane. Whether the satellite is above or below the equatorial plane, the orbital plane will move toward the equator and make a precession motion [27]. The main perturbation solutions of $J_2$ are listed below.

The first-order secular terms of $J_2$ are

$$\left\{\begin{array}{c}{a}_{\sec }^{\left(1\right)}=0,e_{\sec }^{\left(1\right)}=0,i_{\sec }^{\left(1\right)}=0,\hfill \\ {\Omega }_{\sec }^{\left(1\right)}=-{\displaystyle \frac{3J_2}{2p^2}}cosin,\hfill \\ {\omega }_{\sec }^{\left(1\right)}={\displaystyle \frac{3J_2}{2p^2}}\left(2-{\displaystyle \frac{5}{2}}{sin}^{2}i\right)n,\hfill \\ M_{\sec }^{\left(1\right)}={\displaystyle \frac{3J_2}{2p^2}}\left(1-{\displaystyle \frac{3}{2}}{sin}^{2}i\right)\sqrt{1-e^2}n.\hfill \end{array}\right.$$

(4)

The second-order secular terms of $J_2$ are

$$\left\{\begin{array}{cc}\hfill a_{\sec }^{\left(2\right)}& =0,e_{\sec }^{\left(2\right)}=0,i_{\sec }^{\left(2\right)}=0,\hfill \\ \hfill {\Omega }_{\sec }^{\left(2\right)}& =-{\left(\frac{3J_2}{2p^2}\right)}^{2}cosi\left[\left(\frac{3}{2}+\frac{1}{6}{e}^2+\sqrt{1-e^2}\right)-{sin}^{2}i\left(\frac{5}{3}-\frac{5}{24}{e}^2+\frac{3}{2}\sqrt{1-e^2}\right)\right]n,\hfill \\ \hfill {\omega }_{\sec }^{\left(2\right)}& ={\left(\frac{3J_2}{2p^2}\right)}^2\left[\left(4+\frac{7}{12}{e}^2+2\sqrt{1-e^2}\right)-{sin}^{2}i\left(\frac{103}{12}+\frac{3}{8}{e}^2+\frac{11}{2}\sqrt{1-e^2}\right)+{sin}^{4}i\left(\frac{215}{48}-\frac{15}{32}{e}^2+\frac{15}{4}\sqrt{1-e^2}\right)\right]n,\hfill \\ \hfill M_{\sec }^{\left(2\right)}& ={\left(\frac{3J_2}{2p^2}\right)}^2\sqrt{1-e^2}\left[\frac{1}{2}{\left(1-\frac{3}{2}{sin}^{2}i\right)}^2\sqrt{1-e^2}+\left(\frac{5}{2}+\frac{10}{3}{e}^2\right)\right.-{sin}^{2}i\left(\frac{19}{3}+\frac{26}{3}{e}^2\right)\hfill \\ \hfill & +{sin}^{4}i\left(\frac{233}{48}+\frac{103}{12}{e}^2\right)+\frac{e^4}{1-e^2}\left(\frac{35}{12}-\frac{35}{4}{sin}^{2}i+\frac{315}{32}{sin}^{4}i\right)]n.\hfill \end{array}\right.$$

(5)

The first-order short-period terms of $J_2$ are

$$\left\{\begin{array}{cc}\hfill a& {}_{\mathrm{sp}}^{\left(1\right)}=\frac{3J_2}{2a}\left\{\frac{2}{3}(1-\frac{3}{2}{sin}^{2}i)\left[\frac{a^3}{r^3}-{(1-e^2)}^{-3/2}\right]+{sin}^{2}i\frac{a^3}{r^3}cos\left(2f+2\omega \right)\right\},\hfill \\ \hfill e& {}_{\mathrm{sp}}^{\left(1\right)}=\left(\frac{3J_2}{2p^2}\right)\{\frac{1}{3}(1-\frac{3}{2}{sin}^{2}i)\left[e\left(\frac{1}{1+\sqrt{1-e^2}}+\sqrt{1-e^2}\right)+cosf\left(3(1+ecosf)+{(cosf)}^2\right)\right]\hfill \\ \hfill & +\frac{1}{2}{sin}^{2}i\left[\left(e+cosf\left(3(1+ecosf)+{(ecosf)}^2\right)\right)cos2(f+\omega )-\left(1-e^2\right)\left(cos(f+2\omega )+\frac{1}{3}cos(3f+2\omega )\right)\right]\},\hfill \\ \hfill i& {}_{\mathrm{sp}}^{\left(1\right)}\left(t\right)=\left(\frac{3J_2}{2p^2}\right)sin2i\left\{\frac{e}{4}cos(f+2\omega )+\frac{1}{4}cos(2f+2\omega )+\frac{e}{12}cos(3f+2\omega )\right\},\hfill \\ \hfill \Omega & {}_{\mathrm{sp}}^{\left(1\right)}\left(t\right)=-\left(\frac{3J_2}{2p^2}\right)cosi\{(f-M+esinf)-\frac{1}{2}\left[esin(f+2\omega )+sin(2f+2\omega )+\frac{e}{3}sin(3f+2\omega )\right]\},\hfill \\ \hfill \omega & {}_{\mathrm{sp}}^{\left(1\right)}\left(t\right)=-cosi{\Omega }_{\mathrm{sp}}^{\left(1\right)}\left(t\right)+{\left[{\omega }_{\mathrm{sp}}^{\left(1\right)}\left(t\right)\right]}_1,\hfill \\ \hfill [\omega & {}_{\mathrm{sp}}^{\left(1\right)}{\left(t\right)]}_1=\frac{1}{e}\left(\frac{3J_2}{2p^2}\right)\left\{\left(1-\frac{3}{2}{sin}^{2}i\right)\left[(f-M+esinf)e+\left(1-\frac{e^2}{4}\right)sinf+\frac{e}{2}sin2f+\frac{e^2}{12}sin3f\right]\right.\hfill \\ \hfill & +{sin}^{2}i[-\left(\frac{1}{4}-\frac{7}{16}{e}^2\right)sin(f+2\omega )+\frac{3}{4}esin2(f+\omega )+\left(\frac{7}{12}+\frac{11}{48}{e}^2\right)sin(3f+2\omega )\hfill \\ \hfill & +\frac{3}{8}esin(4f+2\omega )+\frac{e^2}{16}(sin(5f+2\omega )+sin(f-2\omega ))\left]\right\},\hfill \\ \hfill M& {}_{\mathrm{sp}}^{\left(1\right)}\left(t\right)=-\sqrt{1-e^2}{\left[{\omega }_{\mathrm{sp}}^{\left(1\right)}\left(t\right)\right]}_1+\left(\frac{3J_2}{2p^2}\right)\sqrt{1-e^2}\hfill \\ \hfill & \left\{\left(1-\frac{3}{2}{sin}^{2}i\right)(f-M+esinf)+{sin}^{2}i\left[\frac{3}{4}esin(f+2\omega )+\frac{3}{4}sin(2f+2\omega )+\frac{1}{4}esin(3f+2\omega )\right]\right\}.\hfill \end{array}\right.$$

(6)

The first-order long-period terms of $J_2$ are

$$\left\{\begin{array}{cc}\hfill a_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =0,\hfill \\ \hfill e_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =\left(\frac{3J_2}{2p^2}\right)\frac{2{sin}^{2}i}{4-5{sin}^{2}i}\left(\frac{7}{24}-\frac{5}{16}{sin}^{2}i\right)(1-e^2)ecos2\omega ,\hfill \\ \hfill i_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =-\left(\frac{3J_2}{2p^2}\right)\frac{sin2i}{4-5{sin}^{2}i}\left(\frac{7}{24}-\frac{5}{16}{sin}^{2}i\right)e^{2}cos2\omega ,\hfill \\ \hfill {\Omega }_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =-\left(\frac{3J_2}{2p^2}\right)\frac{cosi}{{\left(4-5{sin}^{2}i\right)}^2}\left(\frac{7}{3}-5{sin}^{2}i+\frac{25}{8}{sin}^{4}i\right)e^{2}sin2\omega ,\hfill \\ \hfill {\omega }_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =-\left(\frac{3J_2}{2p^2}\right)\frac{1}{{\left(4-5{sin}^{2}i\right)}^2}[{sin}^{2}i\left(\frac{25}{3}-\frac{245}{12}{sin}^{2}i+\frac{25}{2}{sin}^{4}i\right)\hfill \\ \hfill & -e^2\left(\frac{7}{3}-\frac{17}{2}{sin}^{2}i+\frac{65}{6}{sin}^{4}i-\frac{75}{16}{sin}^{6}i\right)]sin2\omega ,\hfill \\ \hfill M_{\mathrm{lp}}^{\left(1\right)}\left(t\right)& =\left(\frac{3J_2}{2p^2}\right)\frac{\sqrt{1-e^2}}{4-5{sin}^{2}i}{sin}^{2}i\left[\left(\frac{25}{12}-\frac{5}{2}{sin}^{2}i\right)-e^2\left(\frac{7}{12}-\frac{5}{8}{sin}^{2}i\right)\right]sin2\omega .\hfill \end{array}\right.$$

(7)

By substituting short-period terms (Equation (6)) and long-period terms (Equation (7)) into the first equation in Equation (3), the mean elements and osculating elements conversion can be realized.

2.2. Perturbation Solution of Atmospheric Drag

The secular terms of atmospheric drag are

$$\left\{\begin{array}{cc}\hfill a_{\sec ,\mathrm{drag}}& =-B_1{a}^{2}n\left[(I_0+2e{I}_1)+Ccos2\omega I_2+\nu z_0^2\left(\frac{3}{4}{I}_0-I_1+\frac{1}{4}{I}_2\right)\right],\hfill \\ \hfill e_{\sec ,\mathrm{drag}}& =-B_{1}an[\left(\frac{e}{2}{I}_0+I_1+\frac{e}{2}{I}_2\right)+\frac{C}{2}cos2\omega (I_1+I_3)\hfill \\ \hfill & +\nu z_0^2\left(-\frac{1}{2}{I}_0+\frac{7}{8}{I}_1-\frac{1}{2}{I}_2+\frac{1}{8}{I}_3\right)],\hfill \\ \hfill i_{\sec ,\mathrm{drag}}& =-\frac{1}{4}{B}_{2}asini(I_0+cos2\omega I_2),\hfill \\ \hfill {\Omega }_{\sec ,\mathrm{drag}}& =-\frac{1}{4}{B}_{2}asin2\omega I_2,\hfill \\ \hfill {\omega }_{\sec ,\mathrm{drag}}& =-cosi{\Omega }_{\sec ,\mathrm{drag}}-B_{1}an\frac{1}{e}Csin2\omega \left[\frac{e}{4}{I}_0-\frac{1}{2}{I}_1-e{I}_2+\frac{1}{2}{I}_3+\frac{3}{4}e{I}_4\right],\hfill \\ \hfill M_{\sec ,\mathrm{drag}}& =-({\omega }_{\sec ,\mathrm{drag}}+cosi{\Omega }_{\sec ,\mathrm{drag}})-\frac{3n}{4a}{a}_{\sec ,\mathrm{drag}}(t-t_0),\hfill \end{array}\right.$$

(8)

where $\nu$ is the atmospheric density scale height rate, $I_m\left(z\right)$ is the Bessel function of the first type of virtual variables, and variable z abides by

$$z=\frac{ae}{H_{p_0}}.$$

(9)

Other intermediate quantities are given as follows:

$$I_m\left(z\right)=\sum _{k⩾0}\frac{1}{k!(m+k)!}{\left(\frac{z}{2}\right)}^{m+2k},$$

(10)

$$B_1=\frac{C_{D}S}{m}{\rho }_{p_0}{F}^{2}exp\left(-\frac{1}{H_{p_0}}(a-a_0+a_0{e}_0)-Ccos2{\omega }_0\right),$$

(11)

$$B_2=\frac{C_{D}S}{m}{\rho }_{p_0}F{n}_{e}exp\left(-\frac{1}{H_{p_0}}(a-a_0+a_0{e}_0)-Ccos2{\omega }_0\right),$$

(12)

$$F=1-\frac{r_{p_0}{n}_e}{v_{p_0}}cos{i}_0,C=\frac{1}{2}\left(\frac{{\epsilon }_e}{H_{p_0}}{r}_{p_0}\right){sin}^2{i}_0,$$

(13)

where S is the cross-sectional reference area, m is the satellite mass, $C_D$ is the drag coefficient, $n_e$ is the Earth spin rate, the elements $(a_0,e_0,i_0,{\omega }_0)$ are the initial Keplerian elements, ${\rho }_{p_0}$ and $H_{p_0}$ are atmospheric density and the density scale height of the initial perigee of the satellite orbit, and $r_{p_0}$ and $v_{p_0}$ are the position and velocity at the same point.

3. Uncertainty Propagation Approach of Constellation Phase

The uncertainty propagation problem satisfying the dynamic constraints Equation (1) can be described as: Given the initial state ${\xi }_0$ and the PDF $p(t_0,{\xi }_0)$ of the initial state deviation $\delta {\xi }_0$, find the PDF $p(t,\xi )$ of the state $\xi \left(t\right)$ of the dynamical system at any moment or its statistical moments of each order. Generally only the first two order statistical moments of the deviation $\delta \xi$, i.e., the mean and covariance matrix, are of interest.

Because controllers, thrusters, and observers have a direct influence on the instantaneous state of the satellite, and the terminal state quantity we are interested in is the mean elements of the satellite, the samples are generated on the satellite’s osculating trajectory, and the terminal samples appear around the mean trajectory. Therefore, the initial state errors are modeled as a probability distribution of three-dimensional position/velocity components $x_0={[r,v]}^{\mathrm{T}}$ of the two co-plane satellites at the moment $t_0$. Since the short period effect is removed by averaging in the dynamics model, only long and long-period effects are retained, so the different initial positions of satellites in one orbital plane do not influence the evolution process. It means that the relative phase evolution of any two co-plane satellites can represent the divergence of this whole orbital plane. The relative phase $\Delta \varphi ={\varphi }_1\left(t\right)-{\varphi }_2\left(t\right)$ of the two satellites after time t is chosen as the terminal state quantity. The schematic diagram of two co-plane satellites’ uncertainty propagation is shown in Figure 1. Since the orbits of LEO mega-constellations tend to be near-circular, the argument of latitude $\varphi =\omega +M$ is used to characterize the phase angle of a satellite. The evolution of the initial state to the terminal state involves, in turn, the transformation of the position-velocity components to the osculating elements, the transform between osculating elements and mean elements, and the semi-analytical propagation. The whole physical process shows strong nonlinearities, so the Monte Carlo is an appropriate analysis approach.

3.1. Semi-Analytical-Based Monte Carlo Simulation Design

Without loss of generality, we assume that the initial state quantities of each co-planar satellite in the satellite orbit coordinate system follow a three-dimensional Gaussian distribution with zero mean $X_0\sim N(0,\Sigma )$, and the components are independent of each other. Then, the joint Gaussian PDF is

$$p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{3/2}{|\Sigma |}^{1/2}}exp\left[-\frac{1}{2}{x}^{\mathrm{T}}{\Sigma }^{-1}x\right],$$

(14)

where the covariance matrix of the initial position velocity errors ${\Sigma }_r$ and ${\Sigma }_v$ are

$$\left\{\begin{array}{cc}\hfill {\Sigma }_r& =diag({\sigma }_x^2,{\sigma }_y^2,{\sigma }_z^2),\hfill \\ \hfill {\Sigma }_v& =diag({\sigma }_{v_x}^2,{\sigma }_{v_y}^2,{\sigma }_{v_z}^2),\hfill \end{array}\right.$$

(15)

where diag represents a diagonal matrix with the entries on the diagonal.

To further simplify the problem, Monte Carlo simulations are performed in two hypothetical cases: Case I, where both two satellites have only initial position errors, and the position standard deviation is ${\sigma }_x={\sigma }_y={\sigma }_z={\sigma }_r$; Case II, where both two satellites have only velocity errors, and the velocity error is noted as ${\sigma }_{v_x}={\sigma }_{v_y}={\sigma }_{v_z}={\sigma }_v$.

A number of N samples are generated, and the samples are then propagated over time individually, and the a posteriori uncertainty distribution is reconstructed from the integrated samples. Multiple sets of position (velocity) standard deviations ${\sigma }_r$ (${\sigma }_v$) with orbital propagation time t are taken, and multiple Monte Carlo simulations are performed to count the first two order moments of the relative phase distribution, which are the mean and standard deviation. Combining hypothesis testing and surface fitting, an approximate mapping between the initial covariance and the terminal covariance can be obtained.

It should be noted that since the relative phase is a directional statistic, and $0^{\circ }$ and ${360}^{\circ }$ are identical angles, so that for example ${180}^{\circ }$ is not a sensible mean of $1^{\circ }$ and ${359}^{\circ }$. To avoid this error, the directional mean $E\left(\theta \right)$ and the directional standard deviation $\sigma \left(\theta \right)$ of the directional quantity $\theta$ should be calculated by Equation (16) [28]

$$\left\{\begin{array}{cc}\hfill E\left(\theta \right)& ={tan}^{-1}(\overline{S}/\overline{C}),\hfill \\ \hfill \sigma \left(\theta \right)& =\sqrt{-2ln\overline{R}},\hfill \end{array}\right.$$

(16)

where the center of mass $(\overline{C},\overline{S})$ of the directional statistic is

$$\overline{C}=\frac{1}{N}\sum _{j=1}^{N}cos{\theta }_i,\overline{S}=\frac{1}{N}\sum _{j=1}^{N}sin{\theta }_i,$$

(17)

and the mean resultant length is

$$\overline{R}={({\overline{C}}^2+{\overline{S}}^2)}^{1/2}.$$

(18)

3.2. Hypothesis Testing of Monte Carlo Simulations Results

The overall terminal phase quantities obtained from the Monte Carlo simulation are denoted as $X_t$. The Gaussian distribution remains normal after propagation by the linear or linearized dynamic equations [29], while the two-body motion plays a dominant role in the orbital dynamics, and the $J_2$ regimens and the atmospheric drag regimens are small quantities. Thus, the terminal quantities can be considered to obey the same Gaussian distribution $X_t\sim N({\mu }_t,{\sigma }_t)$. In order to test whether the mean ${\mu }_t$ of the relative phase distribution is zero, a statistical hypothesis test is required [30]. Two mutually opposing hypotheses are made for $X_t$

$$H_0:{\mu }_t=0,H_1:{\mu }_t\ne 0.$$

(19)

The test statistic distinguishing the first hypothesis from the alternative hypothesis is z value and is given by

$$z=\frac{E(\Delta \varphi )-0}{\sigma (\Delta \varphi )/\sqrt{N}}.$$

(20)

The significance level $\alpha$ is given and the standard Gaussian distribution table is queried to determine the rejection domain. If $\left|z\right|<z_{\alpha /2}$ is satisfied, the hypothesis $H_0$ is accepted (the opposing hypothesis $H_1$ is rejected), i.e., it is considered that, ideally, two co-plane satellites will not diverge or cluster under perturbations; if $\left|z\right|\ge z_{\alpha /2}$, the hypothesis $H_1$ is accepted (the hypothesis $H_0$ is rejected), i.e., it is considered that ideally the non-uncertainty case, the relative phase will also have drifted.

3.3. Product-Based Least-Squares Surface FITTING Approach

Multiple sets of Monte Carlo simulations with different orbital propagation times t and position (velocity) standard deviations ${\sigma }_r$ (${\sigma }_v$) are selected. The time series is denoted as $\mathbf{t}$ (days), the position or velocity standard deviation series is denoted as $\sigma$ ($\mathrm{m}$ or $\mathrm{m}/\mathrm{s}$), the standard deviation of the phase difference obtained by the Monte Carlo approach is denoted as $\mathit{\varphi }$ ($\mathrm{rad}$), and $\mathit{t},\mathit{\sigma },\mathit{\varphi }$ are all row vectors, which approximately form a three-dimensional surface. The set of functions $\left\{{\phi }_p\left(t\right){\psi }_q\left(\sigma \right)\right\}$ is used as the basis function for the product-type least-squares fit [31], where

$$\begin{array}{cc}\hfill & {\phi }_p\left(t\right)=t^p,p=0,1,\dots ,\mathrm{N},\hfill \\ \hfill & {\psi }_q\left(\sigma \right)={\sigma }^q,q=0,1,\dots ,\mathrm{M}.\hfill \end{array}$$

(21)

Construct a surface with $c_{pq}$ as parameter, which can be expressed as

$$\Gamma (t,\sigma )=\sum _{q=0}^{\mathrm{N}}\sum _{p=0}^{\mathrm{M}}{c}_{pq}{\phi }_p\left(t\right){\psi }_q\left(\sigma \right).$$

(22)

where, the coefficient $c_{pq}$ is a matrix form given by

$$\mathit{C}={\left({\mathit{B}}^{\mathrm{T}}\mathit{B}\right)}^{-1}{\mathit{B}}^{\mathrm{T}}\mathit{\varphi }\mathit{G}{\left({\mathit{G}}^{\mathrm{T}}\mathit{G}\right)}^{-1},$$

(23)

with

$$\begin{array}{cc}\hfill \mathit{B}& =[{\mathit{\phi }}_0\left({\mathit{t}}^{\mathrm{T}}\right),\dots ,{\mathit{\phi }}_{\mathrm{N}}\left({\mathit{t}}^{\mathrm{T}}\right)],\hfill \\ \hfill \mathit{G}& =[{\mathit{\psi }}_0\left({\sigma }^{\mathrm{T}}\right),\dots ,{\mathit{\psi }}_{\mathrm{M}}\left({\mathit{\sigma }}^{\mathrm{T}}\right)].\hfill \end{array}$$

(24)

In this way, an approximate function of the standard deviation of the initial position velocity and the standard deviation of the phase difference can be obtained.

4. Numerical Simulations and Results

We will now demonstrate the method for a specific problem. Because the Starlink constellation is a typical LEO mega-constellation with great research significance, we choose the orbit and mechanical environmental parameters of Starlink for simulation. Take into consideration the satellite with mean Keplerian elements as shown in

Table 1. Initial mean elements.

Elements	Values
a	6921 $\mathrm{km}$
e	0.0001
i	53 deg
$\Omega$	10 deg
$\omega$	10 deg
M	60 deg

and mechanical environmental parameters as shown in

Table 2. Analytical model coefficients [32].

Parameters	Values	Explanation
$J_2$	$1082.62668355\times {10}^{-6}$	zonal harmonic
$n_e$	$7.2921158553\times {10}^{-5}$ $\mathrm{rad}/\mathrm{s}$	Earth spin rate
$R_e$	6378.137 $\mathrm{km}$	Earth’s equatorial radius
$\mu$	$3.98600436\times {10}^{14}$ ${\mathrm{m}}^3/{\mathrm{s}}^2$	Earth’s gravitational constant
$\nu$	0.1	variability of scale height
${\epsilon }_e$	1/298.257	Earth flattening
H	68.7 $\mathrm{km}$	scale height
${\rho }_0$	$2.34\times {10}^{-13}$ $\mathrm{kg}/{\mathrm{m}}^3$	atmospheric density
$C_D$	2.2	drag coefficient
m	227 $\mathrm{kg}$	mass of satellite
S	1.2 ${\mathrm{m}}^2$	cross-sectional reference area of satellite

to conduct Monte Carlo simulations. Orbital and structural parameters of Starlink satellites can be found in SpaceX’s FCC documents (such as 36 FCC Rcd 7995 (11)(2021)). In our simulations, the semi-analytical approach shows a strong advantage in computing efficiency. On a PC equipped with i7-1065G7 CPU 1.30 GHz, the numerical and the semi-analytical approach respectively take 4.258 $\mathrm{s}$ and 0.103 $\mathrm{s}$ under $J_2$ and atmospheric drag for the 7-day orbit prediction.

Without loss of generality, the initial position and velocity uncertainty of the satellite is assumed to be Gaussian distributed, and are defined in the geocentric inertial coordinate system. The mean is zero and the covariance matrix is

$$\left\{\begin{array}{cc}\hfill {\Sigma }_r& =\mathrm{diag}(100,100,100){\left(\mathrm{m}\right)}^2,\hfill \\ \hfill {\Sigma }_v& =\mathrm{diag}(0.01,0.01,0.01){(\mathrm{m}/\mathrm{s})}^2.\hfill \end{array}\right.$$

(25)

Four thousand points are distributed to be consistent with the initial uncertainty, generated and propagated using the averaged $J_2$ and atmosphere drag dynamics, referring to the second equation of Equation (3). The secular terms of $J_2$ and atmosphere drag are the same as Equations (4), (5), and (8). The conversion between osculating elements and mean elements abides by the first equation of Equation (3), and short-period terms and long-period terms are described in Equations (6) and (7).

Figure 2 shows initial points distribution before and after conversion from osculating elements to mean elements. The mean track propagated in averaged $J_2$ and atmosphere drag dynamics is marked as a green line, while the osculating track propagated in full dynamics—which includes short-period and long-period terms besides secular terms—is marked as an orange line. To better display the size of the uncertainty ellipsoid, the absolute position of the points is subtracted from the nominal position, and denoted by $\delta$. It can be seen that it is the conversion between osculating elements and mean elements that makes the distribution of sample points show strong nonlinearity at initial moments. Figure 3 shows points distribution after different revolutions up to 100 orbits, which is about 6.2 days. The samples are gradually correlated with the mean trajectory as the propagation time progresses. After 1 orbit period, the shape of points takes on a shuttle shape. Hence, it is reasonable to assume that the position nonlinearity becomes weaker after several revolutions. The argument of latitude of a satellite is Gaussian distribution after going through the same revolutions.

Table 3. Relative phase statistics and hypothesis testing.

Propagation Time	Mean (rad)	Standard Deviation	z Value	Accept ${\mathit{H}}_0$
0.3 orbit	$-2.442\times {10}^{-6}$	$7.656\times {10}^{-5}$	−2.017	No
0.5 orbit	$-3.693\times {10}^{-6}$	$1.223\times {10}^{-4}$	−1.909	Yes
1 orbit	$-6.819\times {10}^{-6}$	$2.403\times {10}^{-4}$	−1.795	Yes
5 orbits	$-3.184\times {10}^{-5}$	$1.196\times {10}^{-3}$	−1.683	Yes
10 orbits	$-6.317\times {10}^{-5}$	$2.393\times {10}^{-3}$	−1.669	Yes

shows the first two order moments of the relative phase of two co-plane satellites in the constellation, z value of hypothesis test Equation (20) and whether $H_0$ should be accepted are also presented in this table.

There is only a ${15}^{\circ }$ difference in relative position between the two satellites, and other initial mean elements of them are the same. When the significance level of hypothesis test $\alpha$ is given as 0.05, the z value of the two-sided test is 1.96. Thus, the hypothesis $H_0$ is rejected in the beginning and is accepted after 1 orbit, because the standard deviation is growing faster than the mean over time. As a consequence, the phase of two co-plane satellites can also be considered as a Gaussian distribution with zero mean as propagation time goes by.

In the multi-group MC simulation, the total number of samples per group is set to be 4000, the propagation time is chosen as $\mathit{t}=\{0.5,1,\dots ,5\}$ days, the position standard deviation is ${\sigma }_r=\{40,80,\dots ,1000\}$ $\mathrm{m}$, and the velocity standard deviation is ${\sigma }_v=\{0.04,0.08,\dots ,1\}$ $\mathrm{m}/\mathrm{s}$. There are a total of 625 groups in each case and are parallel computed.

The standard deviation surfaces of the relative phase in two cases are shown in Figure 4, and the least-squares surface fitting parameters are shown in

Table 4. Least-squares fitting coefficient matrix.

Coefficient	${\mathit{C}}_{\mathbf{r}}$			${\mathit{C}}_{\mathbf{v}}$
	1	$\mathit{y}$	${\mathit{y}}^{\mathbf{2}}$	1	$\mathit{y}$	${\mathit{y}}^{\mathbf{2}}$
1	$-3.163\times {10}^{-7}$	$-1.249\times {10}^{-4}$	$-8.969\times {10}^{-7}$	$1.368\times {10}^{-6}$	$6.991\times {10}^{-6}$	$-1.846\times {10}^{-7}$
x	$3.477\times {10}^{-9}$	$5.896\times {10}^{-5}$	$5.945\times {10}^{-9}$	$-1.092\times {10}^{-5}$	$5.278\times {10}^{-2}$	$4.280\times {10}^{-7}$
$x^2$	$-3.491\times {10}^{-12}$	$-7.843\times {10}^{-10}$	$-8.068\times {10}^{-12}$	$1.794\times {10}^{-5}$	$-2.711\times {10}^{-5}$	$-8.750\times {10}^{-7}$

. Taking $\alpha =0.01$, the z value of the two-sided test is 2.58, and the acceptance rate of hypothesis $H_0$ reaches 96% in both two scenarios, indicating that our assumption is logical.

Since the fitting relationship of the initial standard deviation to the terminal standard deviation has been obtained, we can easily make deviation forecasts for controllers with known error magnitudes, or infer the control and observation error of a constellation backward from its deviation evolution. For example, if two co-plane satellites have an initial uncertainty of 100 $\mathrm{m}$ standard deviation in three-dimensional position, after 20 orbits, the relative phase appears an uncertainty of $0.{438}^{\circ }$ standard deviation, which could be also be caused by 0.109 $\mathrm{m}/\mathrm{s}$ three-dimensional velocity errors.

5. Conclusions

In this paper, an uncertainty propagation analysis approach based on the semi-analytical dynamic approach, Monte Carlo simulation, and least-squares fitting is proposed. The overall approach shows significant benefits in terms of improving the computational efficiency of uncertainty propagation and is suitable for LEO mega-constellations scenarios. On this basis, numerical tests confirm that the probability density will become highly nonlinear after the conversion between osculating elements and mean elements, and the positions distribution of the satellite fits its mean trajectory after the initial nonlinear moment about one orbit period, leading to the phase of two co-plane satellites can be considered as a Gaussian distribution with zero mean as propagation time goes by. Ulteriorly, multi-group Monte Carlo simulations under the effect of initial normal position or velocity error and hypothesis testing verify the assumption. Finally, the initial standard deviation, propagation time, and terminal standard deviation can be described by a binary equation, which can estimate the phase error or infer initial uncertainties.

The major limitation of the present study is that the initial position and velocity errors are considered separately in two cases, but the two coexist in practice, and the other is that our approach is only suitable for low Earth orbit scenarios.