Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Ionospheric Oscillation with Periods of 6–30 Days at Middle Latitudes: A Response to Solar Radiative, Geomagnetic, and Lower Atmospheric Forcing
Next Article in Special Issue
Learning Lightweight and Superior Detectors with Feature Distillation for Onboard Remote Sensing Object Detection
Previous Article in Journal
Recent Advances and Challenges in the Seismo-Electromagnetic Study: A Brief Review
Previous Article in Special Issue
A Robust Star Identification Algorithm Based on a Masked Distance Map
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Self-Organizing Control of Mega Constellations for Continuous Earth Observation

1
College of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China
2
Huzhou Institute of Zhejiang University, Huzhou 313000, China
3
School of Aerospace Engineering, Tsinghua University, Beijing 100084, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(22), 5896; https://doi.org/10.3390/rs14225896
Submission received: 22 October 2022 / Revised: 16 November 2022 / Accepted: 16 November 2022 / Published: 21 November 2022
(This article belongs to the Special Issue CubeSats Applications and Technology)

Abstract

:
This work presents a novel self-organizing control method for mega constellations to meet the continuous Earth observation requirements. In order to decrease the TT&C pressure caused by numerous satellites, constellation satellites are not controlled according to the designed configurations but are controlled with respect to intersatellite constraints. By analyzing the street-of-coverage (SOC) of coplanar constellation satellites, the continuous coverage constraint of the mega constellation is transformed into constraints of the right ascension of ascending node (RAAN) and relative motion bound between every two adjacent coplanar satellites. The proposed continuous coverage constraint can be satisfied by most ongoing or planned mega constellations. Artificial potential functions (APFs) are used to realize self-organizing control. The scale-independent relative orbital elements (SIROEs) are innovatively presented as the self-organizing control variables. Using the Gaussian equations and Lyapunov’s theory, the stability of the APF control in quadratic form is proven, from which it can be concluded that the APF control variables of the controlled satellite should have the same time derivative as the target satellite states under two-body Keplerian motion condition, and SIROEs are ideal choices. The proposed controllers and self-organizing rules are verified in the sub-constellation of the GW-2 mega constellation by simulation. The results demonstrate the goodness in control effect and ground coverage performance.

Graphical Abstract

1. Introduction

In recent years, mega constellations consisting of thousands of satellites have received widespread attention due to their potential to provide high-speed and full-coverage telecommunication services. Commercial companies, such as SpaceX and OneWeb, have made their plans public to deploy mega constellations in low Earth orbits [1]. SpaceX have launched nearly 2000 Starlink satellites and started to provide global communication and Internet services. China plans to construct two Guowang low-Earth-orbit constellations, totaling 12,992 satellites for telecommunication services [2].
Such giant constellations are wonderful choices for performing Earth observation missions. Since they evenly distribute at low Earth orbits all over the globe, they can be used to collect huge amounts of data on gravity [3] and atmosphere [4] for analysis, and they can obtain high–resolution Earth images in real time [5]. Among these remote sensing missions, the priority is to provide continuous global coverage. Limited by the computation load of analytic and grid methods, traditional coverage analysis methods are suitable for constellations with less than 100 satellites [6,7,8]. For the requirement of quick analysis of mega constellations, Dai et al. studied the geometric configuration of constellation projection points on Earth’s surface [9]. The multi-satellite coverage problem is decomposed into a one-satellite coverage problem to solve both continuous and discontinuous coverage problems. Gong et al. proposed a continuous coverage analysis method based on 2D map theory [10]. The proposed method avoids redundant calculation of intersatellite positions and positions between the satellites and ground targets so that the time complexity does not increase significantly with the increase in satellites. Huang et al. presented a systematic method for the design of continuous global coverage Walker and street-of-coverage constellations, by taking seven critical constellation properties (coverage, robustness, self-induced collision avoidance, launch, build-up, station-keeping, and end-of-life disposal) as design criteria [11]. However, the above studies rely on the overall motion information of the constellation, which will result in central ground control.
To that end, this paper innovatively introduces the self-organizing concept of mega constellation control, and the emphasis is on maintaining continuous coverage. Two traditional constellation maintenance methods are the absolute station-keeping method and the relative station-keeping method [12,13,14]. For mega constellations, the former will result in frequent on-orbit maneuvers, while the latter relies on the overall configuration information, which is also inefficient because of the huge number of satellites. In this paper, self-organizing control indicates a decentralized and distributed control mechanism that requires the relative motion information of neighbor satellites only. By that means, constellation control will not be affected by the failure or development of constellation satellites. Artificial potential function (APF) is a commonly used self-organizing control method for its clear mathematical descriptions. Trajectories based on APF control are not planned explicitly but rather are shaped by the potential field in reaction to the changing environment. APF control has shown its potential in many close proximity scenarios, including satellite cluster reconfiguration [15,16], collision avoidance [17], satellite rendezvous and docking [18,19], etc. Spencer improved the APF from Cartesian coordinate form to relative orbital elements (ROEs) form so that the relative orbital dynamics can be reflected in APF control [20]. Our previous work realized reconfiguration, bounded flight, and collision avoidance for satellite clusters by means of ROE-based APF control [21,22]. However, the relative motion distance between constellation satellites is much larger than that of close clusters. The above control methods cannot be applied to mega constellation scenarios.
For the requirement of continuous Earth observation of mega constellations, this paper analyzes the relationship between orbital planes, coplanar adjacent satellites, and coverage performance and proposes a continuous novel coverage constraint that relates to intersatellite motion only. The new constraint can be applied to the self-organizing control of mega constellations. The APF method is adopted to solve the self-organizing control problem in this work. Considering the large range of relative motion between constellation satellites, the scale-independent relative orbital elements (SIROEs) are proposed as the APF control variables. SIROEs have been proposed in our preliminary work [23,24], which have clear geometric interpretations, concise expressions, and high precision for relative motion with near-circular reference orbit and arbitrary intersatellite separation, which are particularly suitable for mega constellation satellites. The stability of the proposed self-organizing control is proven based on the Gaussian equations and Lyapunov’s theory. An ideal APF control variable of the controlled satellite should have the same time derivative as the target satellite states under a two-body Keplerian motion circumstance. Finally, the GW-2 constellation is simulated to verify the proposed self-organizing control method. The results indicate its goodness of control effect and ground-coverage performance.
The continuous coverage constraint is analyzed in Section 2.1. The configuration requirements of mega constellations are given here. The artificial potential functions for a single satellite and the self-organizing control rules for coplanar satellites are studied in Section 2.2. The detailed stability proofs of the proposed artificial potential functions are given in Section 2.3. The performance of self-organizing control for continuous coverage is evaluated based on a GW-2 constellation in Section 3. Conclusions are provided in Section 4.

2. Materials and Methods

2.1. Continuous Coverage Constraint of Mega Constellation

Considering the numerous orbital planes and satellites in one mega constellation, the continuous coverage constraint is redefined in this section. The configuration requirements for mega constellations that satisfy the coverage constraint are proposed. The intersatellite motion constraints during the lifetime are also analyzed.

2.1.1. Definition of Continuous Coverage

Walker constellation has been applied in many mega constellation missions because of its simple design and regular configuration. The spatial configuration of the Walker constellation can be represented by the total number of satellites N, the number of orbital planes P, and the phase parameters F, where 2 π F / N represents the phase difference between two satellites in adjacent orbits. All satellites are placed in circular orbits at the same altitude and at the same inclination.
The latitude range that can be covered by the constellation is determined by the designed inclination and orbital altitude, and the distribution density of constellation satellites at higher latitudes is much larger than that at lower latitudes. Hence, it is convincible that continuous coverage for the designed latitude range can be achieved as long as low latitude regions can be covered continuously.
The coverage geometry of a single satellite is shown in Figure 1, where the angular radius of the coverage range is α max , the orbital altitude is H, the Earth’s radius is R e , and the minimum elevation angle is E min . The relationship between α max , H, and E min is given by Equation (1).
α max = arccos R e R e + H cos E min E min
The coverage range of adjacent coplanar satellites overlaps with each other, forming a street-of-coverage (SOC), shown in Figure 2. The width of the SOC σ depends on the coverage range α max and the geocentric angle between adjacent coplanar satellites θ , written in Equation (2).
σ = arccos cos 2 α max cos θ
Since mega Walker constellation satellites are located at the same inclination, the SOC of adjacent orbital planes are approximately parallel to each other at low latitudes, shown in Figure 3. If the SOC of adjacent orbital planes could overlap with each other, continuous coverage can be achieved at low latitudes. It is noted that adjacent orbital planes could be co-rotating planes or counter-rotating planes, where satellites in the co-rotating planes move in the same direction ( Δ Ω is close to zero), while satellites in the counter-rotating planes move in the opposite direction ( Δ Ω is close to π ). Therefore, the relative motion relationships between co-rotating satellites are relatively fixed compared to counter-rotating satellites. From the perspective of self-organization, the adjacent orbital planes in the coverage constraint are specified as co-rotating planes.
To sum up, the continuous coverage constraint of the mega constellation is defined as follows: The continuous coverage can be guaranteed for a designed latitude range if the SOC of every two adjacent co-rotating orbital planes overlaps with each other. Theoretically, the existence of adjacent counter-rotating orbital planes can provide more coverage multiplicity.

2.1.2. Configuration Requirements for Mega Constellations

The mega constellation, containing a total of N satellites, is composed of S satellites evenly spaced on each of P orbital plane. The geocentric angle between adjacent coplanar satellites θ is 2 π / S . According to Equation (2), the width of SOC of one orbital plane σ is
σ = arccos cos 2 α max cos 2 π 2 π S S
where
S π π α max α max
To ensure the SOC of every two adjacent co-rotating orbital planes overlaps with each other there exists
σ 2 π 2 π P P
From Equations (3) and (5), it can be implied that
cos 2 π P cos 2 π S cos 2 α max
Therefore, the configuration requirements for mega constellations that satisfy the proposed continuous coverage constraint are given in Equations (4) and (6).
Several ongoing and planned mega constellation missions are listed in Table 1. Let the minimum elevation angle E min takes 20 . The coverage range α max can be calculated by Equation (1) together with the orbital altitude. Most of these mega constellations are proven to satisfy the configuration requirements for continuous coverage, i.e., Equations (4) and (6).

2.1.3. Intersatellite Motion Constraints

Due to perturbations and initial orbit errors, the coplanar satellites will drift from each other, and thus the RAAN differences between coplanar satellites will get larger over time without consideration of constellation control. The drift range of the RAAN difference between adjacent co-rotating orbital planes δ Ω , and the drift range of the geocentric angle between adjacent coplanar satellites θ under the continuous coverage constraint, are worth investigating.
The SOC geometry of the perturbed coplanar satellites is shown in Figure 4, where
AC = arccos cos 2 α max cos θ
DE = θ
EF = δ Ω
It can be calculated that the coverage range of two adjacent coplanar satellites along the track direction is
DF = arccos cos θ cos δ Ω
The coverage range along the normal direction of the orbital plane is
BC = sin θ cos 2 α max sin 2 θ cos 2 θ cos 2 θ cos 2 2 α max sin 2 δ Ω
To ensure the SOC of every two adjacent co-rotating orbital planes overlap with each other, there exists
BC 2 π P
To ensure the coverage range of every two adjacent coplanar satellites along the track direction overlap with each other, there exists
D F max 2 π S
According to Equations (11) and (12), the drift range of θ can be obtained, which is related to the drift range of δ Ω .
θ min θ θ max
where
A 1 = cos 2 2 π P B 1 = cos 2 2 π P 1 + sin 2 δ Ω + cos 2 2 α max C 1 = cos 2 2 π P sin 2 2 α max sin 2 δ Ω θ min = arcsin B 1 B 1 2 4 A 1 C 1 2 A 1 θ max = arcsin B 1 + B 1 2 4 A 1 C 1 2 A 1
According to Equations (10), (13) and (14), the drift range of δ Ω of every two adjacent satellites can be obtained as
0 δ Ω arcsin B 2 A 2
where
A 2 = 2 sin 2 2 π S cos 2 2 π S cos 2 2 π P cos 2 2 π S cos 2 2 α max + cos 2 2 π P cos 2 2 α max B 2 = sin 2 2 π S cos 2 2 π S cos 2 2 π P cos 2 2 α max
The detailed derivation steps and the proofs of real roots are given in Appendix A.
We denote Ω t as the target RAAN under J 2 perturbation of the designed orbital plane. Then the drift range of Ω i of all coplanar satellites away from Ω t can be written as
δ Ω i = Ω i Ω t 1 2 arcsin B 2 A 2 = δ Ω min δ Ω i = Ω i Ω t 1 2 arcsin B 2 A 2 = δ Ω max
We denote the position vectors of two adjacent coplanar satellites in the geocentric inertial frame as r i and r i + 1 . The constraint of θ can be transformed to the relative distance constraint d i , i + 1 by means of cosine theorem.
d i , i + 1 r i 2 + r i + 1 2 2 r i r i + 1 cos θ min d i , i + 1 r i 2 + r i + 1 2 2 r i r i + 1 cos θ max
Considering that mega constellations are nearly circular orbits, the position vectors in the geocentric inertial frame are substituted by semi-major axis a, given in Equation (20).
d i , i + 1 2 a 2 1 cos θ min = d mint d i , i + 1 2 a 2 1 cos θ max = d maxt
To sum up, the intersatellite motion constraints that every two adjacent coplanar satellites should obey during their lifetime are presented in Equations (18) and (20).

2.2. Self-Organizing Control Using Artificial Potential Functions

In this section, the scale-independent relative orbital elements (SIROEs), as well as the corresponding bound formulations, are introduced first. The relative motion bound denotes the maximum or minimum value of the relative distance between two satellites. Based on SIROEs, the artificial potential functions for single satellite control are designed under the continuous coverage constraint. The self-organizing control rules for the coplanar satellites are presented.

2.2.1. Scale-Independent Relative Orbital Elements and Bound Formulations

Scale-independent relative orbital elements are high-precision configuration parameters for relative motion with near-circular reference orbit and arbitrary intersatellite separation. SIROEs have been proposed in our preliminary work [23]. Six independent SIROEs are b 1 , b 2 , b 3 , c , θ 1 , θ 2 , shown in Figure 5 and Figure 6. The relationship between SIROEs and relative position is given in Equation (21). We denote x ˜ , y ˜ , z ˜ as the basic relative motion coordinates. The relative trajectory can be obtained by rotating the basic trajectory by θ 1 , counterclockwise. The basic trajectory projected on the x y plane is the combination of an ellipse (with a semi-major axis of 2 b 1 and a semi-minor axis of b 1 ), a circle (with a radius of b 2 ), and an additional offset (with x-offset of b 3 a cos θ 1 and y-offset of a sin θ 1 ). The relative motion on the z-axis is harmonic.
x y z = cos θ 1 sin θ 1 0 sin θ 1 cos θ 1 0 0 0 1 b 1 cos θ 2 + ψ + b 2 cos 2 ψ + b 3 a cos θ 1 2 b 1 sin θ 2 + ψ b 2 sin 2 ψ + a sin θ 1 c sin ψ = cos θ 1 sin θ 1 0 sin θ 1 cos θ 1 0 0 0 1 x ˜ y ˜ z ˜
Without regard to the minor differences of semi-major axes between constellation satellites, the relative velocities represented by SIROEs are written as
x ˙ x ˙ n n y ˙ y ˙ n n z ˙ z ˙ n n = cos θ 1 sin θ 1 0 sin θ 1 cos θ 1 0 0 0 1 b 1 sin θ 2 + ψ 2 b 2 sin 2 ψ 2 b 1 cos θ 2 + ψ 2 b 2 cos 2 ψ c cos ψ
where n is the mean angular velocity of the reference orbit.
According to Equations (21) and (22), SIROEs can be represented by translational states, given in Equation (23).
c = z 2 + z ˙ z ˙ n n 2 sin ψ = z z c c , cos ψ = z ˙ z ˙ n c n c b 2 = 1 2 a a 2 c 2 b 3 = 1 2 a + a 2 c 2 cos θ 1 = 2 b 3 2 y ˙ + n x + n a + 3 b 2 sin 2 ψ 2 n y x ˙ 2 x ˙ n y 2 n y x ˙ 2 n x + 2 n a + y ˙ 2 y ˙ + n x + n a sin θ 1 = 2 n b 3 2 n x + 2 n a + y ˙ cos θ 1 2 n y x ˙ b 1 = x + a cos θ 1 y sin θ 1 + b 2 cos 2 ψ + b 3 2 + x ˙ n cos θ 1 + y ˙ n sin θ 1 + 2 b 2 sin 2 ψ 2 sin β 1 = x ˙ cos θ 1 + y ˙ sin θ 1 + 2 n b 2 sin 2 ψ n b 1 cos β 1 = x + a cos θ 1 y sin θ 1 + b 2 cos 2 ψ + b 3 b 1 sin θ 2 = sin β 1 cos ψ cos β 1 sin ψ cos θ 2 = cos β 1 cos ψ + sin β 1 sin ψ
According to Equation (21), the squared relative distance between satellites can be written as
w = x 2 + y 2 + z 2 = A 3 s 6 + B 3 s 5 + C 3 s 4 + D 3 s 3 + E 3 s 2 + F 3 s + G 3 1 + s 2 3
where
s = tan ψ 2 , ψ 0 , 2 π
Parameters A 3 , B 3 , C 3 , D 3 , E 3 , F 3 , G 3 are polynomials about SIROEs b 1 , b 2 , b 3 , c , θ 1 , θ 2 . Therefore, A 3 , B 3 , C 3 , D 3 , E 3 , F 3 , G 3 are time-invariant, while s is the only time variable. Differentiating w with respect to s results in
w s = B 3 s 6 + 6 A 3 2 C 3 s 5 + 5 B 3 3 D 3 s 4 + 4 C 3 4 E 3 s 3 + 3 D 3 5 F 3 s 2 + 2 E 3 6 G 3 s 3 + F 3 1 + s 2 4
Let w w s s = 0 , there exists
B 3 s 6 + 6 A 3 2 C 3 s 5 + 5 B 3 3 D 3 s 4 + 4 C 3 4 E 3 s 3 + 3 D 3 5 F 3 s 2 + 2 E 3 6 G 3 s 3 + F 3 = 0
or
s = ±
The real roots of Equation (27) together with s = ± constitute a set of s values that perhaps make the squared relative distance take the maximum or minimum value, which are relative motion bounds. Note that when s = ± , ψ = ± π , the squared relative distance takes w = A 3 .
Assuming that s max is the s value that makes the squared relative distance w take the maximum value and s min is the one that makes w take the minimum value, the maximum bound d i j and the minimum bound d i j for the current time can be written as
d i j = A 3 s max 6 + B 3 s max 5 + C 3 s max 4 + D 3 s max 3 + E 3 s max 2 + F 3 s max + G 3 1 + s max 2 3 , if s max ± d i j = A 3 , if s max = ±
d i j = A 3 s min 6 + B 3 s min 5 + C 3 s min 4 + D 3 s min 3 + E 3 s min 2 + F 3 s min + G 3 1 + s min 2 3 , if s min ± d i j = A 3 , if s min = ±
Relative motion bounds are complicated functions of SIROEs, and SIROEs are functions about translational states, given in Equation (23). Therefore, d i j and d i j are functions about translational states x , y , z , x ˙ , y ˙ , z ˙ .

2.2.2. Design of Artificial Potential Functions for Single Satellite Control

Referring to the concept of relative motion bound, Equation (20) can be rewritten as d i j d maxt , d i j d mint . Together with Equation (18), the proposed continuous coverage constraint relates to RAAN differences and relative motion bounds. Therefore, the control effects will be performed in cases if δ Ω i > δ Ω max or δ Ω i < δ Ω min or d i j > d maxt or d i j < d mint .
The following artificial potential functions are all designed in quadratic form, given in Equation (31).
ϕ i = X i X t T K X i X t
where X i is the control variable. X t is the expected value of X i . K is the control coefficient matrix, which is a positive definite symmetric matrix. In the geocentric inertial frame, the control effect generated from ϕ i is given as
Δ V i = X i T V i ϕ i X i = X i T X ˙ ϕ i X i X i T Y ˙ ϕ i X i X i T Z ˙ ϕ i X i
where V i = X ˙ , Y ˙ , Z ˙ T represents the geocentric velocity vector of satellite i. Hence, the velocity increments of ϕ i are the negative gradient of ϕ i with respect to satellite i’s geocentric velocity vector V i . The stability proof will be given in Section 4.
(1)
RAAN Control
To limit the RAAN difference, the artificial potential function is designed in Equation (33). The control variable takes X i = sin i i sin Ω i , sin i i cos Ω i T . The control coefficient matrix takes K = diag k i Ω , k i Ω . δ Ω i = Ω i Ω t . i t and Ω t represent the target inclination and RAAN values under J 2 perturbation of the designed orbital plane, respectively. Therefore, ϕ i Ω limits the RAAN difference within the expected range and corrects the inclination meanwhile.
ϕ i , i Ω = k i Ω sin i i sin Ω i sin i t sin Ω t 2 + k i Ω sin i i cos Ω i sin i t cos Ω t 2 , if δ Ω i > δ Ω max δ Ω i < δ Ω min 0 , if δ Ω min δ Ω i δ Ω max
According to Equation (32), the velocity increments generated by ϕ i , i Ω in the geocentric inertial frame are
Δ V i , i Ω = X i T V i ϕ i , i Ω X i = ϕ i , i Ω sin i i sin Ω i sin i i sin Ω i X ˙ + ϕ i , i Ω sin i i cos Ω i sin i i cos Ω i X ˙ ϕ i , i Ω sin i i sin Ω i sin i i sin Ω i Y ˙ + ϕ i , i Ω sin i i cos Ω i sin i i cos Ω i Y ˙ ϕ i , i Ω sin i i sin Ω i sin i i sin Ω i Z ˙ + ϕ i , i Ω sin i i cos Ω i sin i i cos Ω i Z ˙
The relationship between X i and V i is given as
sin i i sin Ω i = Y Z ˙ Z Y ˙ Y Z ˙ Z Y ˙ 2 + Z X ˙ X Z ˙ 2 + X Y ˙ Y X ˙ 2 sin i i cos Ω i = Z X ˙ X Z ˙ Y Z ˙ Z Y ˙ 2 + Z X ˙ X Z ˙ 2 + X Y ˙ Y X ˙ 2
In fact, the out-of-plane control can be avoided by setting an appropriate offset of the semi-major axis to ensure the time derivative of Ω is approximately zero under perturbations. Take the GW-2 constellation, for instance. The drift range of δ Ω is ± 12.14 when E min = 20 according to Equation (18). The drift range of δ Ω is relatively large, where the RAAN control will basically not be implemented in two years.
(2)
Relative Motion-bound Control
To realize the restriction of relative motion bound, the artificial potential function is designed in Equation (36). The control variable takes X i = d i j or X i = d i j , and the control coefficient matrix takes the K = k b o u n d .
ϕ i , b o u n d = k b o u n d d i j d mint 2 , if d i j < d mint k b o u n d d i j d maxt 2 , if d i j > d maxt 0 , if d i j d mint d i j d maxt
Satellite j is designated as the reference satellite. d mint and d maxt are determined by Equation (20).
In the geocentric inertial frame, the velocity increments generated by ϕ i , b o u n d can be obtained by calculating the negative gradient of ϕ i , b o u n d with respect to the relative velocity vector v i = x ˙ , y ˙ , z ˙ T , and then transforming them from satellite j’s orbital frame to the geocentric inertial frame.
Δ V i , b o u n d = Q j 1 X i T v i ϕ i , b o u n d X i = Q j 1 ϕ i , b o u n d d i j d i j x ˙ + ϕ i , b o u n d d i j d i j x ˙ ϕ i , b o u n d d i j d i j y ˙ + ϕ i , b o u n d d i j d i j y ˙ ϕ i , b o u n d d i j d i j z ˙ + ϕ i , b o u n d d i j d i j z ˙
d i j , d i j d i j , d i j x ˙ x ˙ , d i j , d i j d i j , d i j y ˙ y ˙ , d i j , d i j d i j , d i j z ˙ z ˙ can be formulated based on Equations (23), (24), (29) and (30). Q j represents the transformation matrix from the geocentric inertial frame to satellite j’s orbital frame.
Q j = sin Ω j cos i j sinu j + cos Ω j cosu j cos Ω j cos i j sinu j + sin Ω j cosu j sin i j sinu j sin Ω j cos i j cosu j cos Ω j sinu j cos Ω j cos i j cosu j sin Ω j sinu j sin i j cosu j sin Ω j sin i j cos Ω j sin i j cos i j
(3)
Semi-Major Axis Control
The differences in semi-major axes between constellation satellites will lead to secular drift, and hence the relative motion-bound constraint will never be satisfied. Moreover, the semi-major axis will keep declining because of the influence of atmospheric drag. From this perspective, semi-major axis control is essential to mega constellation maintenance.
The artificial potential function is designed to be a quadratic function about a i , given in Equation (39). X i = a i . K = k a . a t is the designed value of the orbital semi-major axis. Semi-major axis control will be implemented when the continuous coverage constraint has been satisfied.
ϕ i , a = k a a i a t 2 , if d i j d mint d i j d maxt δ Ω δ Ω min δ Ω δ Ω max 0 , if d i j < d mint d i j > d maxt δ Ω < δ Ω min δ Ω > δ Ω max
The velocity increment of ϕ i , a in the geocentric inertial frame is
Δ V i , a = X i T V i ϕ i , a X i = ϕ i , a a i a i X ˙ ϕ i , a a i a i Y ˙ ϕ i , a a i a i Z ˙
The relationship between a i and V i is given as
a i = μ X ˙ 2 + Y ˙ 2 + Z ˙ 2 2 μ X 2 + Y 2 + Z 2 1
To sum up, the control potential functions for satellite i are given as follows.
ϕ i = k i Ω sin i i sin Ω i sin i t sin Ω t 2 + k i Ω sin i i cos Ω i sin i t cos Ω t 2 , if δ Ω > δ Ω max δ Ω < δ Ω min k b o u n d d i j d mint 2 , if d i j < d mint k b o u n d d i j d maxt 2 , if d i j > d maxt k a a i a t 2 , if δ Ω δ Ω min δ Ω δ Ω max d i j d mint d i j d maxt

2.2.3. Self-Organizing Control Rules for Coplanar Satellites

For three adjacent coplanar satellites, i 1 , i and i + 1 , the self-organizing control rules are defined as follows.
ϕ i = k i Ω sin i i sin Ω i sin i t sin Ω t 2 + k i Ω sin i i cos Ω i sin i t cos Ω t 2 , if δ Ω > δ Ω max δ Ω < δ Ω min k b o u n d d i , i 1 d mint 2 , if d i , i 1 < d mint d i , i + 1 d mint k b o u n d d i , i + 1 d mint 2 , if d i , i 1 d mint d i , i + 1 < d mint k b o u n d d i , i 1 d maxt 2 , if d i , i 1 > d maxt d i , i + 1 d maxt k b o u n d d i , i + 1 d maxt 2 , if d i , i 1 d maxt d i , i + 1 > d maxt k a a i a t 2 , else
To put it simply, the RAAN control will be performed once the RAAN difference between satellite i and the designed value does not fit the δ Ω range. The relative motion-bound control will work only if one of two adjacent coplanar satellites of satellite i (denoted as i 1 and i + 1 ) does not meet the boundary constraint while the other does. For other cases, the semi-major axis control will be implemented, even if d i , i 1 , d i , i 1 , d i , i + 1 , and d i , i + 1 do not all meet the boundary constraint.

2.3. Stability Proof of Satellite Control

In this section, the stability of the proposed self-organizing control is proved based on the Gaussian equations and Lyapunov’s theory. The selection criteria of potential function variables for satellite control are discussed.

2.3.1. Stability Proof of RAAN Control and Semi-Major Axis Control

If there exists one Lyapunov function ϕ that satisfies Equation (44), the control system can be proven to be stable.
ϕ 0 = 0 , ϕ X > 0 , X 0 ϕ ˙ 0 = 0 , ϕ ˙ X 0 , ϕ ˙ X 0 , X 0
In this paper, ϕ i is designated as the Lyapunov function. Referring to Equation (31), it is obvious that ϕ i is positive definite. If the time derivative of ϕ i is negative definite, one can prove that the proposed APF control is stable.
For RAAN control and semi-major axis control, X = X i X t = [ sin i i sin Ω i sin i t sin Ω t , sin i i cos Ω i sin i t cos Ω t , a i a t ] T . Hence, X relates to controlled satellite i only.
According to Equation (31), the time derivative of ϕ i can be written as
d ϕ i d t = ϕ i X i T d X i T d t + ϕ i X t T d X t T d t
The time derivative of X i is
d X i d t = X i R i T d R i d t + X i V i T d V i d t = X i R i T d R i d t + X i V i T μ R i R i 3 + u = d X i d t t w o b o d y + X i V i T u
where R i and V i are the position vector and velocity vector in the geocentric inertial frame, respectively. u represents the control acceleration vector in the geocentric inertial frame. d X i d X i d t d t t w o b o d y represents the time derivative of X i under the two-body Keplerian motion condition.
Since X t does not change under control, the time derivative of X t is
d X t d t = d X t d t t w o b o d y
According to Equation (32), the control effect generated by ϕ i in geocentric inertial frame is the negative gradient of ϕ i with respect to V i = X ˙ , Y ˙ , Z ˙ T .
u = X i T V i ϕ i X i = X 1 X ˙ X 2 X ˙ X n X ˙ X 1 Y ˙ X 2 Y ˙ X n Y ˙ X 1 Z ˙ X 2 Z ˙ X n Z ˙ ϕ i X 1 ϕ i X 2 ϕ i X n
Substituting Equations (46)–(48) into Equation (45), the time derivative of ϕ i can be rewritten as
d ϕ i d t = ϕ i X i T d X i d t t w o b o d y + X i V i T u + ϕ i X t T d X t d t t w o b o d y = ϕ i X i T d X i d t t w o b o d y X i V i T X i T V i ϕ i X i + ϕ i X t T d X t d t t w o b o d y = X i T V i ϕ i X i T X i T V i ϕ i X i + ϕ i X i T d X i d t t w o b o d y + ϕ i X t T d X t d t t w o b o d y = 4 X i T V i K X i X t T X i T V i K X i X t + 2 X i X t T K d X i d t d X t d t t w o b o d y
For RAAN control and semi-major axis control, X i = sin i i sin Ω i , sin i i cos Ω i , a i T stays constant under the two-body Keplerian motion condition without consideration of perturbations. Hence, d X i d X i d t d t t w o b o d y is zero. Similarly, d X t d X t d t d t t w o b o d y is also zero. The time derivative of ϕ i is
d ϕ i d t = 4 X i T V i K X i X t T X i T V i K X i X t
The control coefficient matrix K is positive definite. Therefore, Equation (50) is a negative definite quadratic function when X i X t 0 . ϕ i is positive definite, and the time derivative of ϕ i is negative definite. According to Lyapunov’s second theory, the RAAN and semi-major axis APFs are stable.

2.3.2. Stability Proof of Relative Motion-Bound Control

For relative motion-bound control, X = X i X t = d i j d maxt , d i j d mint T . Different from RAAN and semi-major axis control, the relative motion-bound control relates to the controlled satellite i and its adjacent coplanar satellite j. Therefore, the time derivative of X i will be different. Here we give the stability proof of relative motion-bound control under self-organizing control rules.
Assuming that ϕ i = k bound d i , i + 1 d maxt 2 . It can be inferred from Equation (43) that d i , i 1 d maxt , and d i , i + 1 > d maxt . In this case, X i = d i , i + 1 . ϕ i is positive definite. The time derivative of X i can be written as
d X i d t = X i R i T d R i d t + X i R i + 1 T d R i + 1 d t + X i V i T d V i d t + X i V i + 1 T d V i + 1 d t = d X i d t t w o b o d y + X i V i T u i + X i V i + 1 T u i + 1 = X i V i T u i + X i V i + 1 T u i + 1
The time derivative of X t is
d X t d t = d X t d t t w o b o d y = 0
The control effect working on satellite i is
u i = X i T V i ϕ i X i
Since d i , i 1 d maxt and d i , i + 1 > d maxt , the control APF working on satellite i + 1 is
ϕ i + 1 = k bound d i , i + 1 d maxt 2 = ϕ i , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 d maxt 0 , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 > d maxt
The control effect working on satellite i + 1 is
u i + 1 = X i T V i + 1 ϕ i X i , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 d maxt 0 , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 > d maxt
Hence, the time derivative of ϕ i is
d ϕ i d t = ϕ i X i T d X i d t + ϕ i X t T d X t d t = ϕ i X i T X i V i T u i + X i V i + 1 T u i + 1 = ϕ i X i T X i V i T X i T V i ϕ i X i ϕ i X i T X i V i + 1 T X i T V i + 1 ϕ i X i , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 d maxt ϕ i X i T X i V i T X i T V i ϕ i X i , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 > d maxt = 4 k bound 2 d i , i + 1 d maxt 2 d i , i + 1 V i T d i , i + 1 V i 4 k bound 2 d i , i + 1 d maxt 2 d i , i + 1 V i + 1 T d i , i + 1 V i + 1 , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 d maxt 4 k bound 2 d i , i + 1 d maxt 2 d i , i + 1 V i T d i , i + 1 V i , i f d i , i + 1 > d maxt a n d d i + 1 , i + 2 > d maxt
Since k bound > 0 , the time derivative of ϕ i is negative definite. To that end, the stability of the relative motion-bound control under self-organizing rules is proven.

2.3.3. Selection Criteria of Quadratic Potential Function Variables

In past studies, the control variables of quadratic APFs can be position states X , Y , Z [16,18,25], velocity states X ˙ , Y ˙ , Z ˙ [19], relative orbital elements [20,21,22], classical orbital elements [26], etc. The stabilities of quadratic APFs in terms of different control variables are discussed as follows.
(1)
Using R and V as Control Variables
In this case, X i = R i ; V i , X t = R t ; V t . The control effects aim to make the geocentric position and velocity of the controlled satellite converge to that of the target one. The difference between the time derivative of X i and that of X t under the two-body condition is
d X i d t d X t d t t w o b o d y = d R i d t d R t d t d V i d t d V t d t = V i V t μ R i R i 3 + μ R t R t 3
When the relative distance between the controlled satellite and the target one is much closer than the orbital altitude, the difference in velocity and two-body gravity in Equation (57) is minor and is ignored. Referring to Equation (49), d X i d X i d t d t d X t d X t d t d t t w o b o d y = 0 , and the time derivative of ϕ i is negative definite. The control effect is stable. However, the orbital dynamics differ a lot with the increase in the intersatellite distance. d X i d X i d t d t d X t d X t d t d t t w o b o d y can be positive or negative, and the stability cannot be guaranteed anymore.
Therefore, quadratic APF control in terms of R and V is stable only if the relative distance is small, or if the gravity is weak, for example, on the ground or in deep space.
(2)
Using Orbital Elements as Control Variables
In this case, X i = a i , e i , i i , Ω i , ω i , M i T , X t = a t , e t , i t , Ω t , ω t , M t T . Under two-body Keplerian motion condition, it can be inferred that
d X i d t d X t d t t w o b o d y = d a i d a i d t d t d a t d a t d t d t d e i d e i d t d t d e t d e t d t d t d i i d i i d t d t d i t d i t d t d t d Ω i d Ω i d t d t d Ω t d Ω t d t d t d ω i d ω i d t d t d ω t d ω t d t d t d M i d M i d t d t d M t d M t d t d t = 0 0 0 0 0 n i n t
According to Equation (49), if the slow variables of orbital elements, i.e., a , e , i , Ω , ω , are chosen as APF control variables, the control effect is stable.
(3)
Using Relative Orbital Elements as Control Variables
In this case, X i = b 1 , b 2 , b 3 , c , θ 1 , β 1 , ψ T , X t = b 1 t , b 2 t , b 3 t , c t , θ 1 t , β 1 t , ψ t T . There exists
d X i d t d X t d t t w o b o d y = d b 1 d b 1 d t d t d b 1 t d b 1 t d t d t d b 2 d b 2 d t d t d b 2 t d b 2 t d t d t d b 3 d b 3 d t d t d b 3 t d b 3 t d t d t d c d c d t d t d c t d c t d t d t d θ 1 d θ 1 d t d t d θ 1 t d θ 1 t d t d t d β 1 d β 1 d t d t d β 1 t d β 1 t d t d t d ψ d ψ d t d t d ψ t d ψ t d t d t = 0 0 0 0 0 n i n t n i n t
Therefore, relative orbital elements b 1 , b 2 , b 3 , c , θ 1 are ideal control variables that ensure control stability. When a i = a t , β 1 , ψ can be chosen as control variables.
To sum up, the selection criteria of quadratic APF control variables are given as follows. If the time derivative of X i is the same as that of X t under the two-body condition, the corresponding APF control is stable.

3. Results and Discussion

In this section, the self-organizing control under continuous coverage constraints is simulated. The GW-2 constellation is designated as the mega constellation example. The control consumption is assessed, and the coverage performance is analyzed.

3.1. Self-Organizing Control of Two Satellites

The two-satellite case is simulated, where one satellite, named the chief, performs semi-major axis control only, and the other, named the deputy, performs self-organizing control to keep the relative motion bound within the expected range.
The simulation parameters are listed in Table 2.
The relative motion-bound changing with time is shown in Figure 7, where the blue, black, red, and green curves represent the maximum bound under control, the minimum bound under control, the expected maximum bound, and the expected minimum bound, respectively. It can be found that both the maximum bound and the minimum bound are constrained within the expected range.
Due to the existence of atmospheric drag, the mean semi-major axis drops by 10 m every day without semi-major control, shown in Figure 8. It can be seen that the relative motion-bound control makes the semi-major axis of the deputy converge to that of the chief as well. That is because the expected bound values are fixed, and the fixed bounds exist only if the semi-major axis difference takes zero. Figure 9 shows the semi-major axes of the chief satellite and the deputy under semi-major axis control. The results indicate that the semi-major axes of the two satellites have been maintained around the expected value.
The control effect of the deputy is shown in Figure 10. The accumulated velocity increment is 1.98 m/s, and the single-step velocity increment is less than 0.02 m/s, which is minor enough for mega constellation control.

3.2. Self-Organizing Control of Coplanar Satellites

The coplanar satellites are selected as satellites of one orbit in GW-2 sub-constellation with an orbital altitude of 1145 km and inclination of 60 . The designed configuration of the GW-2 constellation has been presented in Table 1. There are 48 orbital planes and 36 coplanar satellites in the GW-2 sub-constellation total. The simulation parameters are listed in Table 3. The initial errors of δ Ω , δ ω , δ f are relatively large, which simulate the drift of Ω , ω , f under perturbations without control.
The angular radius of the coverage range is calculated as α max = 17.19 . According to Equation (18), the RAAN drift should be within δ Ω min = 12.14 and δ Ω max = 12.14 . The drift range of the relative motion bound can be calculated using Equations (15) and (20) based on the extreme value of δ Ω .
The initial distribution of the coplanar satellites in the geocentric inertial system is shown in Figure 11, while the final distribution is shown in Figure 12. The red dots represent satellites that do not satisfy the continuous coverage constraint, named unqualified satellites, whereas the green dots represent satellites that satisfy the continuous coverage constraint. The number of unqualified satellites changing with time is shown in Figure 13. The unqualified satellites are reduced from 11 to 0 due to the self-organizing control, which verifies the feasibility of the proposed control method.
The differences in the maximum relative motion bounds away from their corresponding constraints are shown in Figure 14. The differences should be negative theoretically and they actually are. Hence the maximum bound control has not been performed during the simulation time. The differences in the minimum relative motion bounds away from their corresponding constraints are shown in Figure 15, from which it can be found that all minimum bounds are controlled to be larger than the constraints. The occurrence of the unqualified satellite at 40 days in Figure 13 can be understood by referring to Figure 15. The small drift in the minimum bound away from the constraint due to perturbations is identified as unqualified.
The semi-major axis differences away from target value a t are shown in Figure 16. The maximum drift of the semi-major axis at the final moment is 54.27 m. The RAAN differences away from the target RAAN under J 2 perturbation of the designed orbital plane are shown in Figure 17. The maximum drift of RAAN at the final moment is 9.73 . The control effects of coplanar satellites in the geocentric inertial system are shown in Figure 18. The exact velocity increments have no value of reference because the initial orbit conditions are dramatically bad. However, Figure 16 reveals the convergence of the proposed self-organizing control.

3.3. Self-Organizing Control of Mega Constellation

This paper focuses on self-organizing control to keep continuous coverage within the designed latitude zones, whereas collision avoidance control has not been considered so far. Since the proposed continuous coverage constraint relates to coplanar satellites only, the self-organizing control of different orbital planes is independent of each other.
The initial simulation parameters are listed in Table 4. In order to demonstrate the effect of relative motion-bound control, the initial errors of δ Ω , δ ω , δ f are set as δ Ω 5 , 5 deg , δ ω 1 , 1 deg , δ f 1 , 1 deg . The initial and final distributions of the mega constellation are shown in Figure 19 and Figure 20, respectively. Satellites that satisfy the continuous coverage constraint are colored green, while others are colored red. It can be observed that unqualified satellites turn into qualified satellites under self-organizing control. The number of unqualified satellites is reduced from 248 to 0, shown in Figure 21. The control effects in the geocentric inertial system changing with time are shown in Figure 22, which reveals the convergence of the proposed self-organizing control.
The time histories of the maximum bound difference, minimum bound difference, semi-major axis difference, and RAAN difference are given in Figure 23, Figure 24, Figure 25 and Figure 26. All the minimum relative motion bounds converge to expected constraints automatically under control. The semi-major axis differences are restricted within 98.02 m, and the maximum drift of RAAN at the final moment is 5.46 , which is far away from the δ Ω drift constraint.
Figure 27 shows the coverage performance of the GW-2 sub-constellation at the final moment, where yellow dots represent virtual ground stations that can be single-covered, and red dots represent ground stations that cannot be single-covered. The results indicate that the GW-2 sub-constellation can provide continuous coverage within ± 78 along latitudes. Figure 28 shows the multiplicity of coverage of the GW-2 sub-constellation, where the blue curve denotes the initial condition and the green curve denotes the final condition. It can be found that the GW-2 constellation always provides continuous coverage within ± 78 latitude zones during the simulation time.

4. Conclusions

In view of the numerous satellites in mega constellations, this paper innovatively introduces the self-organizing concept of mega constellation control. The emphasis is on realizing continuous Earth observation within the designed latitude zones so that the combination of multiple sub-constellations can provide global remote sensing services. A novel continuous coverage constraint has been presented that relates to intersatellite motion only, including the RAAN drift range and relative motion-bound drift range between two adjacent coplanar satellites. The configuration requirements for mega constellations have been studied, and most ongoing or planned mega constellations satisfy the proposed coverage constraint. Considering RAAN, the relative motion-bound and semi-major axis constraints, the corresponding artificial potential functions as well as the self-organizing rules have been designed. As a highlight, the scale-independent relative orbital elements are used as control variables in mega constellation control because these elements have high precision for near-circular orbits and arbitrary intersatellite ranges. The stability of satellite control using quadratic artificial potential functions has been proven. One key finding is that if the time derivative of the control variable is the same as that of the expected value under two-body Keplerian motion condition, the corresponding quadratic artificial potential function is stable. The self-organizing control simulations of coplanar satellites and the GW-2 sub-constellation have been carried out. The simulation results indicate the goodness of control effect for the small single-step velocity increment and low accumulated control consumption, and the mega constellation satisfies continuous coverage under control.

Author Contributions

Conceptualization: Y.X., Y.Z. and Z.W.; funding acquisition: Z.W. and L.F.; Methodology: Y.X.; validation verification: Y.X. and Y.H.; writing—original draft: Y.X.; writing—review and editing: Y.X. and Y.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China of funder grant number No. U20B2056 and the Huzhou Distinguished Scholar Program of funder grant number ZJIHI-KY0016. The APC was funded by Zhejiang University.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

According to Equations (11) and (12),
cos BC cos 2 π P
cos 2 2 α max sin 2 θ cos 2 2 π P sin 2 θ cos 2 θ sin 2 δ Ω cos 2 θ + cos 2 2 α max sin 2 δ Ω
Let τ = sin 2 θ , 1 τ = cos 2 θ ,
cos 2 π P τ 2 + cos 2 2 π P 1 + sin 2 δ Ω cos 2 2 α max τ cos 2 2 π P sin 2 2 α max sin 2 δ Ω 0
Let
A 1 = cos 2 2 π P B 1 = cos 2 2 π P 1 + sin 2 δ Ω + cos 2 2 α max C 1 = cos 2 2 π P sin 2 2 α max sin 2 δ Ω
According to Equation (6), it can be proven that
cos 2 π P cos 2 α max
Therefore,
A 1 0 B 1 cos 2 2 α max sin 2 δ Ω 0 C 1 0
The range of τ can be obtained by Equation (A3)
B 1 B 1 2 4 A 1 C 1 2 A 1 τ B 1 + B 1 2 4 A 1 C 1 2 A 1
According to Equation (A6), it is obvious that
B 1 + B 1 2 4 A 1 C 1 2 A 1 B 1 B 1 2 4 A 1 C 1 2 A 1 0
Therefore, τ = sin 2 θ 0 can be proven, yet B 1 2 4 A 1 C 1 0 depends on the δ Ω value.
The range of sin θ is
B 1 B 1 2 4 A 1 C 1 2 A 1 sin θ B 1 + B 1 2 4 A 1 C 1 2 A 1
The range of cos θ is
2 A 1 + B 1 + B 1 2 4 A 1 C 1 2 A 1 cos θ 2 A 1 + B 1 B 1 2 4 A 1 C 1 2 A 1
According to Equations (10) and (13),
cos D F min cos 2 π S
cos θ min cos 2 π S cos δ Ω
B 1 2 4 A 1 C 1 2 A 1 sin 2 2 π S 2 A 1 cos 2 2 π S sin 2 δ Ω B 1
B 1 2 4 A 1 C 1 2 A 1 sin 2 2 π S + 2 A 1 cos 2 2 π S sin 2 δ Ω + B 1 2
Equation (A14) ensures that θ is a real value, given in Equations (A9) and (A10). Equation (A14) can be simplified as
A 2 sin 2 δ Ω B 2 0
where
A 2 = 2 sin 2 2 π S cos 2 2 π S cos 2 2 π P cos 2 2 π S cos 2 2 α max + cos 2 2 π P cos 2 2 α max B 2 = sin 2 2 π S cos 2 2 π S cos 2 2 π P cos 2 2 α max
According to Equation (6),
A 2 2 sin 2 2 π S cos 2 2 α max cos 2 2 π S cos 2 2 α max + cos 2 2 π P cos 2 2 α max = cos 2 2 π P + cos 2 2 π S 2 cos 2 2 α max = sin 2 2 π P + sin 2 2 π S cos 2 2 α max 0 B 2 = sin 2 2 π S cos 2 2 π S cos 2 2 π P cos 2 2 α max 0
δ Ω is a positive variable. Therefore, the range of δ Ω is
0 δ Ω arcsin B 2 A 2

References

  1. Del Portillo, I.; Cameron, B.G.; Crawley, E.F. A technical comparison of three low earth orbit satellite constellation systems to provide global broadband. Acta Astronaut. 2019, 159, 123–135. [Google Scholar] [CrossRef]
  2. Williams, A.; Hainaut, O.; Otarola, A.; Tan, G.H.; Rotola, G. Analysing the impact of satellite constellations and ESO’s role in supporting the astronomy community. arXiv 2021, arXiv:2108.04005. [Google Scholar]
  3. Boggio, M.; Colangelo, L.; Virdis, M.; Pagone, M.; Novara, C. Earth Gravity In-Orbit Sensing: MPC Formation Control Based on a Novel Constellation Model. Remote Sens. 2022, 14, 2815. [Google Scholar] [CrossRef]
  4. Wang, Z.; Zhang, Y.; Wen, G.; Bai, S.; Cai, Y.; Huang, P.; Han, D.; He, Y. Atmospheric Density Model Optimization and Spacecraft Orbit Prediction Improvements Based on Q-Sat Orbit Data. arXiv 2021, arXiv:2112.03113. [Google Scholar]
  5. Mastro, P.; Masiello, G.; Serio, C.; Pepe, A. Change Detection Techniques with Synthetic Aperture Radar Images: Experiments with Random Forests and Sentinel-1 Observations. Remote Sens. 2022, 14, 3323. [Google Scholar] [CrossRef]
  6. Lang, T. Low Earth orbit satellite constellations for continuous coverage of the mid-latitudes. In Proceedings of the Astrodynamics Conference, San Diego, CA, USA, 29 July–31 July 1996; p. 3638. [Google Scholar]
  7. Santos, M.; Shapiro, B. Relating satellite coverage to orbital geometry. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, HI, USA, 18–21 August 2008; p. 6609. [Google Scholar]
  8. Hackett, T.M.; Bilén, S.G.; Bell, D.J.; Lo, M.W. Geometric approach for analytical approximations of satellite coverage statistics. J. Spacecr. Rocket. 2019, 56, 1286–1299. [Google Scholar] [CrossRef]
  9. Dai, G.; Chen, X.; Wang, M.; Fernández, E.; Nguyen, T.N.; Reinelt, G. Analysis of satellite constellations for the continuous coverage of ground regions. J. Spacecr. Rocket. 2017, 54, 1294–1303. [Google Scholar] [CrossRef] [Green Version]
  10. Gong, Y.; Zhang, S.; Peng, X. Quick coverage analysis of mega Walker Constellation based on 2D map. Acta Astronaut. 2021, 188, 99–109. [Google Scholar] [CrossRef]
  11. Huang, S.; Colombo, C.; Bernelli-Zazzera, F. Multi-criteria design of continuous global coverage Walker and Street-of-Coverage constellations through property assessment. Acta Astronaut. 2021, 188, 151–170. [Google Scholar] [CrossRef]
  12. Arnas, D.; Casanova, D.; Tresaco, E. Relative and absolute station-keeping for two-dimensional–lattice flower constellations. J. Guid. Control. Dyn. 2016, 39, 2602–2604. [Google Scholar] [CrossRef]
  13. Chen, Y.; Zhao, L.; Liu, H.; Li, L.; Liu, J. Analysis of configuration and maintenance strategy of LEO walker constellation. J. Astronaut. 2019, 40, 1296–1303. [Google Scholar]
  14. Li, J.; Hu, M.; Wang, X.; Li, F.; Xu, J. Analysis of configuration offsetting maintenance method for LEO Walker constellation. Chin. Space Sci. Technol. 2021, 41, 38. [Google Scholar]
  15. Izzo, D.; Pettazzi, L. Autonomous and distributed motion planning for satellite swarm. J. Guid. Control. Dyn.s 2007, 30, 449–459. [Google Scholar] [CrossRef]
  16. Nag, S.; Summerer, L. Behaviour based, autonomous and distributed scatter manoeuvres for satellite swarms. Acta Astronaut. 2013, 82, 95–109. [Google Scholar] [CrossRef]
  17. Hu, Q.; Dong, H.; Zhang, Y.; Ma, G. Tracking control of spacecraft formation flying with collision avoidance. Aerosp. Sci. Technol. 2015, 42, 353–364. [Google Scholar] [CrossRef]
  18. McInnes, C.R. Autonomous proximity manoeuvring using artificial potential functions. ESA J. 1993, 17, 159–169. [Google Scholar]
  19. Cao, L.; Qiao, D.; Xu, J. Suboptimal artificial potential function sliding mode control for spacecraft rendezvous with obstacle avoidance. Acta Astronaut. 2018, 143, 133–146. [Google Scholar] [CrossRef]
  20. Spencer, D.A. Automated trajectory control using artificial potential functions to target relative orbits. J. Guid. Control. Dyn. 2016, 39, 2142–2148. [Google Scholar] [CrossRef]
  21. Wang, Z.; Xu, Y.; Jiang, C.; Zhang, Y. Self-organizing control for satellite clusters using artificial potential function in terms of relative orbital elements. Aerosp. Sci. Technol. 2019, 84, 799–811. [Google Scholar] [CrossRef]
  22. Xu, Y.; Wang, Z.; Zhang, Y. Bounded flight and collision avoidance control for satellite clusters using intersatellite flight bounds. Aerosp. Sci. Technol. 2019, 94, 105425. [Google Scholar] [CrossRef]
  23. Xu, Y.; Wang, Z.; Zhang, Y. Modeling and control of scale-independent relative orbital elements for near-circular orbits. Acta Astronaut. 2022, 198, 642–658. [Google Scholar]
  24. Jiang, C.; Wang, Z.; Zhang, Y. Decomposition analysis of spacecraft relative motion with different inter-satellite ranges. Acta Astronaut. 2019, 163, 56–68. [Google Scholar]
  25. An, M.; Wang, Z.; Zhang, Y. Self-organizing control strategy for asteroid intelligent detection swarm based on attraction and repulsion. Acta Astronaut. 2017, 130, 84–96. [Google Scholar] [CrossRef]
  26. Sun, Y.; Shen, H. The control of mega-constellation at low earth orbit based on TLE. Acta Astronaut. 2020, 42, 156. [Google Scholar]
Figure 1. Satellite coverage geometry.
Figure 1. Satellite coverage geometry.
Remotesensing 14 05896 g001
Figure 2. Street-of-coverage geometry of one orbital plane.
Figure 2. Street-of-coverage geometry of one orbital plane.
Remotesensing 14 05896 g002
Figure 3. Streets-of-coverage geometries of adjacent orbital planes.
Figure 3. Streets-of-coverage geometries of adjacent orbital planes.
Remotesensing 14 05896 g003
Figure 4. Street-of-coverage geometry of the perturbed coplanar satellites. D and E are two sub-satellite points. A and C are two intersections of two coverage circle. Point A, B, C form a right triangle. Point D, E, F form a right triangle.
Figure 4. Street-of-coverage geometry of the perturbed coplanar satellites. D and E are two sub-satellite points. A and C are two intersections of two coverage circle. Point A, B, C form a right triangle. Point D, E, F form a right triangle.
Remotesensing 14 05896 g004
Figure 5. Illustration of the scale-independent relative orbital elements ( x y plane).
Figure 5. Illustration of the scale-independent relative orbital elements ( x y plane).
Remotesensing 14 05896 g005
Figure 6. Illustration of the scale-independent relative orbital elements ( x z plane).
Figure 6. Illustration of the scale-independent relative orbital elements ( x z plane).
Remotesensing 14 05896 g006
Figure 7. Relative motion-bound changing with time.
Figure 7. Relative motion-bound changing with time.
Remotesensing 14 05896 g007
Figure 8. Semi-major axis time histories without semi-major axis control.
Figure 8. Semi-major axis time histories without semi-major axis control.
Remotesensing 14 05896 g008
Figure 9. Semi-major axis time histories under semi-major axis control.
Figure 9. Semi-major axis time histories under semi-major axis control.
Remotesensing 14 05896 g009
Figure 10. Control effects changing with time.
Figure 10. Control effects changing with time.
Remotesensing 14 05896 g010
Figure 11. The initial distribution of the coplanar satellites in the geocentric inertial system. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.
Figure 11. The initial distribution of the coplanar satellites in the geocentric inertial system. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.
Remotesensing 14 05896 g011
Figure 12. The final distribution of the coplanar satellites in the geocentric inertial system.
Figure 12. The final distribution of the coplanar satellites in the geocentric inertial system.
Remotesensing 14 05896 g012
Figure 13. Number of unqualified satellites in the coplanar satellites.
Figure 13. Number of unqualified satellites in the coplanar satellites.
Remotesensing 14 05896 g013
Figure 14. The differences in maximum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are maximum bound differences.
Figure 14. The differences in maximum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are maximum bound differences.
Remotesensing 14 05896 g014
Figure 15. The differences in minimum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are minimum bound differences.
Figure 15. The differences in minimum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are minimum bound differences.
Remotesensing 14 05896 g015
Figure 16. Semi-major axis differences away from the target value.
Figure 16. Semi-major axis differences away from the target value.
Remotesensing 14 05896 g016
Figure 17. RAAN differences away from the target RAAN under J 2 perturbation. Red lines represent δ Ω min and δ Ω max respectively. Blue lines represent RAAN differences.
Figure 17. RAAN differences away from the target RAAN under J 2 perturbation. Red lines represent δ Ω min and δ Ω max respectively. Blue lines represent RAAN differences.
Remotesensing 14 05896 g017
Figure 18. Control effects of coplanar satellites in the geocentric inertial system.
Figure 18. Control effects of coplanar satellites in the geocentric inertial system.
Remotesensing 14 05896 g018
Figure 19. Initial distribution of the mega constellation. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.
Figure 19. Initial distribution of the mega constellation. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.
Remotesensing 14 05896 g019
Figure 20. Final distribution of the mega constellation.
Figure 20. Final distribution of the mega constellation.
Remotesensing 14 05896 g020
Figure 21. Number of unqualified satellites in the mega constellation.
Figure 21. Number of unqualified satellites in the mega constellation.
Remotesensing 14 05896 g021
Figure 22. Control effects of constellation satellites.
Figure 22. Control effects of constellation satellites.
Remotesensing 14 05896 g022
Figure 23. Maximum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are maximum bound differences.
Figure 23. Maximum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are maximum bound differences.
Remotesensing 14 05896 g023
Figure 24. Minimum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are minimum bound differences.
Figure 24. Minimum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are minimum bound differences.
Remotesensing 14 05896 g024
Figure 25. Semi-major axis differences of constellation satellites.
Figure 25. Semi-major axis differences of constellation satellites.
Remotesensing 14 05896 g025
Figure 26. RAAN differences of constellation satellites. Red lines represent δ Ω min and δ Ω max respectively. Blue lines represent RAAN differences.
Figure 26. RAAN differences of constellation satellites. Red lines represent δ Ω min and δ Ω max respectively. Blue lines represent RAAN differences.
Remotesensing 14 05896 g026
Figure 27. Coverage performance of the GW-2 sub-constellation.
Figure 27. Coverage performance of the GW-2 sub-constellation.
Remotesensing 14 05896 g027
Figure 28. Multiplicity of coverage of the GW-2 sub-constellation at the initial and final moments.
Figure 28. Multiplicity of coverage of the GW-2 sub-constellation at the initial and final moments.
Remotesensing 14 05896 g028
Table 1. Summary of several mega constellation missions.
Table 1. Summary of several mega constellation missions.
ConstellationAltitudeInclinationPlanesSatellites per PlaneWhether Continuous Coverage Constraint Is Satisfied
Starlink Genaration 1550 km53 7222Yes
540 km53.2 7222Yes
570 km70 3620Yes
560 km97.6 658No
560 km97.6 443No
OneWeb1120 km87.9 1840Yes
1150 km53 3250Yes
1110 km53.8 3250Yes
GW-A59590 km85 1630No
600 km50 4050Yes
508 km55 6060Yes
GW-21145 km30 4836Yes
1145 km40 4836Yes
1145 km50 4836Yes
1145 km60 4836Yes
Table 2. Simulation parameters for the two-satellites case.
Table 2. Simulation parameters for the two-satellites case.
ParametersValues
Simulation Time60 days
Step Length600 s
Initial Conditions a c , e c , i c , Ω c , ω c , M c = R e + 550 km , 0 , 53 , 30 , 40 , 50
a d , e d , i d , Ω d , ω d , M d = a c , e c + 0.01 , i c , Ω c + 1 , ω c + 1 , M c + 10
Perturbations J 2 , atmospheric drag, lunar gravity, solar gravity
Control Coefficients k bound = 2 × 10 8 , k a = 2 × 10 8
Ballistic Coefficient C D A C D A m m = 0.01
Expected Range d maxt = d max 1 km , d mint = d min + 1 km , a t = R e + 500 km
Table 3. Simulation parameters for the coplanar satellites case.
Table 3. Simulation parameters for the coplanar satellites case.
ParametersValues
Simulation Time60 days
Step Length600 s
Designed Configuration a t = R e + 1145 km
e t = 0
i t = 60
Ω t = 0
ω t = 0
Initial Orbit Errors δ a = 100 , 100 m , uniform
δ e = 0 , 0.1 , uniform
δ i = 0.1 , 0.1 deg , uniform
δ Ω = 10 , 10 deg , uniform
δ ω = 1 , 1 deg , uniform
δ f = 1 , 1 deg , uniform
Perturbations J 2 , atmospheric drag, lunar gravity, solar gravity
Control Coefficients k bound = 2 × 10 8 , k a = 2 × 10 8
Ballistic Coefficient C D A C D A m m = 0.01
Minimum Elevation Angle E min = 20
Table 4. Simulation parameters for the mega-constellation case.
Table 4. Simulation parameters for the mega-constellation case.
ParametersValues
Simulation Time60 days
Step Length600 s
Designed Configuration a t = R e + 1145 km
e t = 0
i t = 60
P = 48
S = 36
Initial Orbit Errors δ a = 100 , 100 m , uniform
δ e = 0 , 0.1 , uniform
δ i = 0.1 , 0.1 deg , uniform
δ Ω = 5 , 5 deg , uniform
δ ω = 1 , 1 deg , uniform
δ f = 1 , 1 deg , uniform
Perturbations J 2 , atmospheric drag, lunar gravity, solar gravity
Control Coefficients k bound = 2 × 10 8 , k a = 2 × 10 8
Ballistic Coefficient C D A C D A m m = 0.01
Minimum Elevation Angle E min = 20
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Xu, Y.; Zhang, Y.; Wang, Z.; He, Y.; Fan, L. Self-Organizing Control of Mega Constellations for Continuous Earth Observation. Remote Sens. 2022, 14, 5896. https://doi.org/10.3390/rs14225896

AMA Style

Xu Y, Zhang Y, Wang Z, He Y, Fan L. Self-Organizing Control of Mega Constellations for Continuous Earth Observation. Remote Sensing. 2022; 14(22):5896. https://doi.org/10.3390/rs14225896

Chicago/Turabian Style

Xu, Yun, Yulin Zhang, Zhaokui Wang, Yunhan He, and Li Fan. 2022. "Self-Organizing Control of Mega Constellations for Continuous Earth Observation" Remote Sensing 14, no. 22: 5896. https://doi.org/10.3390/rs14225896

APA Style

Xu, Y., Zhang, Y., Wang, Z., He, Y., & Fan, L. (2022). Self-Organizing Control of Mega Constellations for Continuous Earth Observation. Remote Sensing, 14(22), 5896. https://doi.org/10.3390/rs14225896

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop