#### Abstract

With the development of commercial space technology, low-orbit large-scale satellite constellation has shown great development potential in communication services and military applications due to its advantages of low delay, strong signal, high coverage rate, fast communication rate, and low cost of mass production. It has become an urgent problem in the field of satellite constellations to study the perturbation evolution law of low-orbit large-scale constellation initialization and orbit maintenance, and how to maintain and control the deployment and initialization of low-orbit satellite constellations efficiently, economically, and stably. In this paper, the satellite motion law, constellation configuration evolution characteristics, and constellation initialization control strategy are analyzed considering the perturbation of earth’s nonspherical gravity and atmospheric drag under the orbit deviation of satellite orbit elements. Firstly, MonteCarlo simulation was used to simulate the deviation of satellite initial orbit elements, and the evolution law of ascending ascension point and phase angle of different constellations was analyzed and, secondly, establish the initialization phase constellation deployment and the ascending node right ascension deployment of mathematical equations, using Starlink constellation (V1.0-L3, 8) satellite TLE data; initialization of Leo constellation orbit plane is analyzed, by considering the orbit bias satellite simulation data; the same orbital plane satellite initialization phase deployment method is analyzed; finally, this paper provides some suggestions for the future deployment and maintenance control strategy of low-orbit constellation.

#### 1. The Introduction

Leo large constellations in the initial deployment process, due to the existing TT&C systems, satellite launch site, ground control, satellite initial deployment in from the nonspherical earth gravity, the moon and sun, sun light pressure gravity, atmospheric drag, etc., all kinds of perturbation forces affect at the same time, has also been orbit satellite, orbit elements deviation, control the influence of the deviation, and lead to spread constellation configuration. The initial deployment of large satellite constellations has become a key technology, restricting the development of satellite constellations. At present, theoretical analysis of large-scale constellation initialization, orbital maintenance perturbation evolution law, and constellation control strategy results are mainly reflected in the following aspects: in literature 1, it was proposed that the overall constellation configuration of Walker constellation in medium-high orbit would drift under the action of orbit entry deviation, control execution deviation, perturbation force, solar and lunar gravity, etc., and the relative control method was adopted to ensure the stability of the overall constellation geometry configuration and the single star orbit to meet specific requirements. Dalin proposed two constellations keeping methods, absolute and relative, and studied the phase keeping in relative station keeping of Walker constellation in high orbit. In [1], by analyzing the track dynamics characteristics of time-sharing lifting orbit, a multisatellite placement method with phase adjustment realized by time-sharing sequence lifting orbit was proposed. However, in the phase initialization phase, the influence of the orbit deviation on the semimajor axis cannot be ignored under the premise of considering the influence of orbit lifting on the phase angle. On this basis, this paper considers the constellation initialization deployment strategy of the orbit deviation.

The problems of initial deployment and maintenance of large-scale constellation in low orbit are mainly in ascending intersection right ascension and relative phase control of orbital plane stability, which are affected by the existing launch deployment and satellite TT&C system capability conditions. There are certain deviations between satellite orbit acquisition and launching into orbit. The difficulty lies in considering the influence of satellite orbit deviation, control execution deviation, perturbation force, solar and lunar gravity force, and complex control model of constellation configuration, and whether the control method is economical and applicable. Based on TLE data of the third and eighth batch (V1.0-L3 and 8) of Starlink system, the initialization and configuration evolution of this configuration satellite under the perturbation of earth’s nonspherical gravity and atmospheric resistance were analyzed, and the evolution characteristics of parameters related to the initial deployment control strategy of the constellation were verified, providing reference suggestions for future low-orbit constellation deployment.

#### 2. Configuration Perturbation and Evolution Analysis of Low-Orbit Constellation

Low orbit satellites are mainly affected by the perturbation force of earth oblateness and atmospheric resistance, and other perturbations are small and can be ignored.

##### 2.1. Perturbation of the Earth’s Oblateness

If only the gravity of the earth center and the perturbation of the oblateness term are considered, the long-term change rate of the mean root number of the satellite orbit is shown as follows [2]:
where is the right ascension _{d}rift rate of the rising point; is the nonspherical perturbation term of the earth; is the phase angle drift rate; is the radius of the earth; is the half path; is the average angular velocity of the satellite movement; is the orbit inclination, and is the eccentricity. For a group of satellites in Walker constellation with exactly the same orbital height, inclination and eccentricity parameters, the long-term change rates of right ascension and phase are the same, but due to the existing launch deployment and satellite tracking and control system capability conditions, there are certain deviations in satellite orbit acquisition and launch.

For the satellite with an initial deviation of , the relative drift rate of right ascension and phase at the ascending intersection is shown as follows [3]:

Under the influence of perturbation, the relative drift of right ascension and phase angle of ascending intersection satisfies the following relation, as shown in the following:

For low orbit near circular orbit satellites, the initial deviation of the semimajor axis and the inclination angle of the orbit can cause the relative drift of the right ascension and phase angle of the rising point, but the influence of each deviation value on the relative drift is different. The orbit elements of constellation 1 (orbit height 0 ~ 30000 km, orbit inclination 53°, semimajor axis deviation 0 ~ 100 m) are brought into Formula (2), and the relative drift curve of right ascension and phase angle of constellation 1 are obtained; the orbit elements of constellation 2 (orbit height 550 km, orbit inclination 53°, semimajor axis deviation 0 ~ 100 m, inclination deviation 0 ~ 0.1°) are brought into Formula (3), and the relative drift curve of right ascension and phase angle of constellation 2 are obtained; among them, Figure 1 shows the change of the right ascension of constellation 1 satellite with the relative drift of the initial semimajor axis deviation. Figure 2 shows the change of the phase angle of constellation 1 satellite with the relative drift of the initial semimajor axis deviation. Figure 3 shows the change of the right ascension of constellation 2 satellite with the relative drift of the initial semimajor axis deviation and the initial inclination deviation. Figure 4 shows the change of phase angle of constellation 2 satellite with the relative drift of initial semi major axis deviation and initial inclination deviation.

It can be concluded that the main influencing factor of the relative drift of the rising point right ascension is the initial deviation of the orbital inclination. The initial deviation of the inclination of 0.1°will cause the drift of the rising point right ascension of nearly 4°in one year; the main influence factor of phase angle relative drift is the initial deviation of half major axis. The horizontal half major axis deviation of 100 m will cause phase angle drift of nearly 40°in one year.

Due to initial orbit satellite orbit elements exist deviation, due to the initial orbit satellite orbit elements exist deviation, coupled with the constellation orbit relative perturbation difference exists long-term cumulative effects, different constellation orbit relative orbit of the ascending node and within the same orbital plane satellite relative phase angle change, causing the constellation’s geometry spread continuously, thus affecting the coverage, stability and service of the constellation.

##### 2.2. Atmospheric Resistance Perturbation

The variation of semimajor axis and eccentricity caused by atmospheric drag perturbation depends on the atmospheric density and the satellite surface to mass ratio. The atmospheric density decreases with the increase of altitude, and its law decreases exponentially with the increase of altitude. At the same time, it is also affected by season, solar activity, time, and geomagnetic activity. The semimajor axis perturbation of the satellite caused by atmospheric drag mainly depends on the satellite’s surface mass ratio and atmospheric density, and its change rate in the orbital period is shown in the following [2]: where is the atmospheric drag coefficient. is the area-mass ratio of the satellite. is the upwind area of the satellite. is the mass of the satellite. is the atmospheric density.

The optical pressure coefficient is 1.0; the atmospheric drag coefficient is 2.3, and and in high solar activity years. Figure 5 shows the half major axis attenuation of satellites with an orbital height of 1100 km in the high solar activity year, with an average half major axis attenuation of 7 m per day.

#### 3. Monte Carlo Simulation Based Evolution Analysis of Orbit Entry Error Distribution

It is assumed that the intake error and acquisition error of satellite orbit elements follow normal distribution, and the distribution functions of semimajor axis, orbital inclination, eccentricity, ascending point right ascension, and initial phase angle are shown in Table 1. After Monte Carlo method shooting experiment, the deviation distribution of initial orbit elements (, ) in low orbit constellation is obtained in Figures 6 and 7.

Figure 8 shows the distribution of orbit plane phase angle drift evolution, Figure 9 shows the orbital plane distribution of the ascending node right ascension drift evolution, the altitude of 550 km, orbital mass Leo Walker constellation of 53°angle under the condition of considering the deviation of initial orbit, to maintain the relative phase angle within ±1.5°; an average of 90 days shall maintain control at a time; maintenance control is required for an average of 30 days to stay within ±0.1°of right ascension.

#### 4. Research on Constellation Initialization Control Strategy Analysis Method

Constellation initialization, that is, constellation configuration is established. Constellation initialization strategy is to work from the initial orbit to orbit satellite orbit (or anchor) transfer strategy; initialization strategy can be classified into passive mode (orbit drift and perturbation factors), active way (thrust system), and the combination of initiative and passive strategy, mainly completed to work (or berth orbit) phase adjustment and orbital maneuver. The main objective of initialization is to distribute the satellites gathered together during the separation of satellites and arrows on different orbital planes, and evenly distribute multiple satellites in each orbital plane.

##### 4.1. Initialization Analysis of Right Ascension at Ascending Intersection of Orbital Plane

The placement task of multiple stars in different orbital planes can be regarded as the adjustment of right ascension at the relative ascending point of the orbital plane. Its essence is to realize the adjustment of right ascension at the ascending point of the orbital plane by using the precession of the earth nonspherical gravity. The specific operation method is as follows: by lifting the orbital heights of satellites of different orbital planes in the constellation, and by means of time accumulation, the objective of controlling the orbital plane difference between satellites can be realized. To sum up, the maneuver of multiple satellites with one arrow can be attributed to the delay of lifting orbit of satellites with different orbital planes.

Assume that the time delay is to be solved, and the relationship between the right ascension difference of ascending point and the right ascension drift of ascending point of satellites on different orbital planes with time is as follows (6):

Referring to the satellite initialization parameters of V1.0-L3 of Starlink constellation, the third batch of 20 satellites were lifted from 280 km altitude after orbit to 350 km and finally to the target altitude of 550 km, climbing at a speed of 5.9 km/s every day. For the orbital altitude of 550 km and nominal thrust, according to Formula (6), the right ascension bias at the ascending intersection of 20° requires a time delay of about 41.3 days.

Therefore, the deployment time of the 60-star orbital plane of the First Arrow is about 130 days, and the constellation orbital plane deployment is completed 130 days after the satellite is put into orbit. Figures 10 and 11, respectively, show the initial orbital altitude climb and right ascension deployment of v1.0-L3 satellites in the Starlink constellation.

##### 4.2. Analysis of In-Plane Phase Angle Initialization

Multisatellite placement in the same orbital plane can be regarded as the relative phase adjustment problem in the orbital plane. By analyzing the orbital dynamics characteristics of time-sharing lifting rails, a multisatellite placement method with phase adjustment realized by time-sharing sequence lifting rails was proposed in reference [4]. Time-sharing orbit lifting maneuver needs to meet the following conditions: the constellation layout and structure are uniform, and the phase deviation between any adjacent satellites is consistent, that is, the time-sharing orbit lifting maneuver of any two adjacent satellites has the same delay [5].

In the phase initialization stage, it is assumed that the lifting delay to be solved is , and the specific solution process is as follows:

the change of the semimajor axis of the satellite will cause the change of the orbital angular velocity and then cause the phase change. Therefore, the relative change rate of the phase angle can be regarded as a function of the change rate of the semimajor axis and the semimajor axis. As shown in the following formula: where takes the satellite to be maneuvered as the reference; can be regarded as a constant.

For a near-circular orbit, during the control phase of orbit uplift, the change of Satellite speed; the relation expression is as follows:

Simultaneous Formula (7), Formula (8):

Therefore, the waiting time for orbit uplift of adjacent satellites is as follows:

Therefore, the time required for phase deployment of satellites on the same orbital plane is as follows:

is the number of satellites in the same orbital plane.

Referring to the V1.0-L3 satellite initialization parameters of Starlink constellation, the first batch of 20 satellites will directly lift from the altitude of 280 km after entering orbit to the target altitude of 550 km, and climb at a speed of 5.9 km every day. In this case, the phase angle offset of 18°requires a time delay of about 6201 s.

In the phase initialization stage, the semimajor axis modification variables play a major role in the relative drift of the phase angle. The deviation of the semimajor axis in the order of 100 meters will cause a phase drift of more than 40°in one year. Therefore, on the premise of considering the influence of track lifting on the phase angle, the influence of the semi major axis deviation and acquisition error cannot be ignored. This chapter focuses on the semimajor axis as the main influence factor, and analyzes the feasibility and effectiveness of the semimajor axis offset compensation phase angle method by establishing the analytical formula of the semimajor axis offset compensation phase angle considering the semimajor axis deviation.

The relationship between the half major axis offset and the phase angle drift rate is as follows:

is the deviation of the semimajor axis.

Therefore, on the basis of compensating for semimajor axis deviation and aiming at improving on-board computer autonomy, this paper proposes the flow and method of on-board autonomous phase deployment control for satellites in the same orbital plane. The calculation steps are as follows: (1)The half major axis difference and phase difference of satellite I are , and , calculated using satellite 1 as a reference.(2)Input the phase angle offset value and electric thrust acceleration of the satellite.(3)According to the phase angle offset value and the electric thrust acceleration , the orbit lifting of two adjacent satellites needs waiting time .(4)Taking satellite 1 as the benchmark, orbit control is performed on satellite i, and the interval time is, respectively, waited to raise the orbit height of the satellite and realize the relative phase angle adjustment

##### 4.3. Simulation and Analysis

HPOP high precision orbit predictor was used to simulate and analyze two low-orbit Walker constellations of different scales. The perturbation forces of the two constellations were nonspherical perturbation and atmospheric resistance perturbation. The atmospheric density model was Jacchia 70, the atmospheric resistance coefficient was 2.2, and the satellite mass was 260 kg. The propulsion system is equipped with continuous electric thrust, nominal thrust , corresponding acceleration m, and satellite surface to mass ratio 0.003. The initial establishment process of orbital surface phase of the simulation constellation is simulated, and the initial time is 00 : 00 : 00 on January 1, 2022, Beijing time.

Constellation configuration 1 is a large-scale low-orbit Walker constellation with an orbital altitude of 550 km and an orbital inclination of 53°. It is assumed that the initial semimajor axis of each star after the constellation is put into orbit is a random variable of ±3.5 km and the initial phase deviation is ±3°.

Figures 12 and 13 show the orbital height uplift of initial phase deployment of different constellation configurations, respectively.

For the large-scale low orbit Walker constellation with a constellation configuration of 60/3/1, an orbital height of 550 km and an orbital inclination of 53°, assuming that the initial semi major axis of each satellite is a random variable of ±3.5 km and the initial phase deviation is ±3°; it will take about 22 days to realize the initial deployment of 20 satellites on the first orbital plane with a phase difference of 18°between adjacent satellites. For the large-scale low orbit Walker Constellation with a constellation configuration of 100/4/1, an orbital height of 1100 km and an orbital inclination of 88°, assuming that the initial semi major axis of each satellite is a random variable of ±5 km and the initial phase deviation is ±4°, it will take about 9 days to realize the initial deployment of 25 satellites on the first orbital plane with a phase difference of 9°.

#### 5. Conclusion

In this paper, combined with the measured data of Starlink constellation satellite and the simulation data considering the initial orbit deviation, the constellation deployment initialization process is quantitatively analyzed, and the feasibility of ascending right ascension and phase angle deployment strategy is analyzed in a relatively complete way. The simulation data show that for the large-scale low-orbit Walker constellation with an orbital altitude of 550 km and an orbital inclination of 53°, the relative phase angle should be maintained within ±1.5°, and the control should be maintained once in 90 days on average, taking into account the initial orbital deviation. Maintenance control is required for an average of 30 days to stay within ±0.1°of right ascension. For the large-scale low-orbit Walker constellation with constellation configuration of 60/3/1, orbital altitude of 550 km and orbital inclination of 53°, it is assumed that the initial semimajor axis of each star is a random variable of ±3.5 km after the constellation is put into orbit, and the initial phase deviation is ±3°. It takes about 22 days to realize the initial deployment of adjacent satellites with phase difference of 18°on the first orbital plane.

For the large-scale low-orbit Walker constellation with constellation configuration 2 of 200/5/1, orbital altitude of 1100 km and orbital inclination of 86°, it is assumed that the initial semimajor axis of each star is a random variable of ±5 km, and the initial phase deviation is ±1°after the constellation is put into orbit. It takes about 4 days to realize the initial deployment of adjacent satellites with phase difference of 18° on the first orbital plane. This paper can provide reference for the initial deployment and maintenance of large-scale constellation in low orbit, and will follow up the launch and deployment of Starlink system in the future. By analyzing the measured data of in-orbit satellites, relevant conclusions of absolute maintenance control strategy and relative maintenance control strategy are presented.

#### Data Availability

Previously reported data were used to support this study and are available at References [2, 6].

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was funded by the National Natural Science Foundation of China.