Analytical Solution and Satellite Phasing Rules for Designing Dedicated Geosynchronous Orbit Satellite Constellations

2 MATHEMATICAL MODELING

In this section, we derive a dedicated analytical solution for GSO designs in an Earth-centered, Earth-fixed (ECEF) coordinate system. Lee and Hall presented a comprehensive solution for the relative motion of satellites through parametric constellations (PCs) (Lee & Hall, 2021). As depicted in Figure 1, the relative motions are designed as hypocycloids, epicycloids, or their transformed parametric curves based on the correlation between the deferent circle and epicycle. Hypocycloid curves were drawn with the epicycle rolling inside the deferent circles, and epicycloid curves were drawn with the epicycle circle rolling outside.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Hypocycloid and epicycloid curves by deferent circle and epicycle

We define a notation to describe the equations of motion for the relative motion of the satellites. In this study, “c” indicates the chief satellite, “d” denotes the deputy satellite, and the six orbital elements representing their motion are denoted by the general notations [a, e, i, Ω, ω, M₀]. The true anomaly, which is a function of time, is denoted by v, and the orbital radius is denoted by r. The PC theory does not consider perturbations, and the coordinate system represents the relative positions x, y, z of the deputy satellite with respect to the rotating reference frame of the chief satellite using the following generalized equations of motion:

$x=\frac{r_d}{2}\left(\left(1+\cos i_R\right)\cos \left(\Delta {\theta }^{-}+{\psi }_{xy}^{-}\right)+\left(1-\cos i_R\right)\cos \left(\Delta {\theta }^{+}+{\psi }_{xy}^{+}\right)\right)-r_c$ 1a

$y=\frac{r_d}{2}\left(\left(1+\cos i_R\right)\sin \left(\Delta {\theta }^{-}+{\psi }_{xy}^{-}\right)-\left(1-\cos i_R\right)\sin \left(\Delta {\theta }^{+}+{\psi }_{xy}^{+}\right)\right)$ 1b

$z=r_d\sin i_R\sin \left(v_d+{\psi }_z\right)$ 1c

where the terms of the angular frequency $\Delta {\theta }^{-}$ , $\Delta {\theta }^{+}$ and the phase angles ${\psi }_{xy}^{-}$ , ${\psi }_{xy}^{+}$ , ${\psi }_z$ , as well as ϕ_c and ϕ_d, are given as follows:

$\Delta {\theta }^{-}=v_d-v_c$ 2a

$\Delta {\theta }^{+}=v_d+v_c$ 2b

${\psi }_{xy}^{-}=\left({\omega }_d-{\omega }_c\right)-\left({\varphi }_d-{\varphi }_c\right)$ 2c

${\psi }_{xy}^{+}=\left({\omega }_d+{\omega }_c\right)-\left({\varphi }_d+{\varphi }_c\right)$ 2d

${\psi }_z={\omega }_d-{\varphi }_d$ 2e

${\varphi }_c={\tan }^{-1}\frac{\left(\sin \Delta \Omega \sin i_c\sin i_d\right)}{\left(-\cos i_d+\cos i_c\cos i_R\right)}$ 3a

${\varphi }_d={\tan }^{-1}\frac{\left(\sin \Delta \Omega \sin i_c\sin i_d\right)}{\left(\cos i_c-\cos i_d\cos i_R\right)}$ 3b

Here, i_R is the relative inclination between two satellites:

$\cos i_R=\cos i_c\cos i_d+\sin i_c\sin i_d\cos \Delta \Omega$ 4

The relative ascending node ΔΩ is defined as follows:

$\Delta \Omega ={\Omega }_d-{\Omega }_c$ 5

The purpose of this study was to obtain a closed-form solution for designing a GSO from the equations of relative motion described above. The resulting solution represents the relative position of the GSO satellite with respect to a coordinate system in the International Terrestrial Reference Frame (ITRF). For this purpose, the parameters in Equations (2)–(4) were simplified as shown in Table 1.

The angular frequency terms in Table 1 vary as functions of the true anomaly of the satellite orbit and the rotational velocity of the Earth. The nodal drift $\stackrel{.}{\Omega }$ is as follows (Kuiack & Ulrich, 2018):

$\stackrel{.}{\Omega }=-\frac{3}{2}{J}_{2}n\frac{R_e^2}{p^2}\cos i$ 6

where R_e is the Earth’s radius, J₂ is the oblateness perturbation, n is the mean motion of the satellite orbit, and $p=a(1-e^2)$ The phase angle terms ${\psi }_{xy}^{-}$ ,${\psi }_{xy}^{+}$ are also simplified as functions of ω and Ω of the satellite orbit. Consequently, by substituting the simplified parameters listed in Table 1 in Equation (1), we obtain a closed-form solution for computing the coordinates in the ITRF, a non-inertial frame of reference with its origin at the Earth’s center of mass:

$x=\frac{r}{2}\left[\left(1+\cos i\right)\cos \left(v-\left({\omega }_{\oplus }-\dot{\Omega }\right)t+\omega +\Omega \right)+\left(1-\cos i\right)\cos \left(v+\left({\omega }_{\oplus }-\dot{\Omega }\right)t+\omega -\Omega \right)\right]$ 7a

$y=\frac{r}{2}\left[\left(1+\cos i\right)\sin \left(v-\left({\omega }_{\oplus }-\dot{\Omega }\right)t+\omega +\Omega \right)-\left(1-\cos i\right)\sin \left(v+\left({\omega }_{\oplus }-\dot{\Omega }\right)t+\omega -\Omega \right)\right]$ 7b

$z=r\sin i\sin \left(v+\omega \right)$ 7c

As shown in Equation (7), the phase angle terms of the orbit are determined solely by ω and Ω, which will play important roles in the design of the GSO constellation in the next section. Moreover, Equation (7) is generally used to design the orbits of satellites rotating around Earth. To design the trajectory of the target GSO satellite, only the semi-major axis of the GSO (42,164.17 km) must be considered (Sjoberg et al., 2017).

Figure 2 depicts the orbit of a satellite near the GSO. Figure 2(a) shows the trajectory for a semi-major axis of 41,000 km and an inclination angle of 30°, where the values of the other orbital elements are all zero. The trajectory of the satellite follows a hypocycloid motion, and a loop is formed as it traces outward. In Figure 2(b), with a semi-major axis of 43,000 km, a loop forming inward can be observed, along with an epicycloid motion. In general, the loops are caused by the dynamics of the relative motion between Earth and the satellite, and the shape and number of loops are determined by the rotation rate between the two objects. Additionally, the trajectory featured in Figure 2 has a shape that differs from a figure-eight, owing to the difference in rotation rates between Earth and the satellite. This result indicates that if the rotation rates of the two objects are the same, the trajectory has a figure-eight shape and has characteristics identical to those of the GSO trajectory.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Hypocycloid and epicycloid motions near the GSO (a) Hypocycloid motion near GSO (b) Epicycloid motion near GSO

3 SATELLITE PHASING RULE AND CLOSED-FORM SOLUTIONS

In this section, we aim to derive a satellite phasing rule for designing a GSO constellation with the purpose of having identical and non-identical constellation patterns. First, the identical constellation pattern design, in which all satellites in the constellation have the same ground track, is considered only for Ω and M₀ among the six orbital elements. For GSO satellites, the rotation rates of the Earth and the satellite are the same. This characteristic establishes the following relationship between Ω, the initial position of the Earth’s rotation, and M₀, the initial position of the satellite’s rotation, for maintaining the same ground track:

$M_0=-\Omega$ 8

Therefore, the phasing rule for the same ground track of k GSO satellites can be arranged as follows:

$\begin{array}{cc}{\Omega }_k={\Omega }_1+{\theta }_{{\Omega }_k}\left(k-1\right),& k=1,2,3\dots N\end{array}$ 9

with:

$\begin{array}{cc}{M}_{k0}=M_{10}-{\Omega }_k,& k=1,2,3\dots N\end{array}$ 10

where ${\theta }_{{\Omega }_k}$ is an interval of ascending node spacing.

For a GSO satellite constellation, which consists of the elliptical orbit from Equation (10), the argument of perigee, ω, for satellites with the same eccentricity is added to M_k0. The resulting phasing rule is as follows:

$\begin{array}{cc}{\Omega }_k={\Omega }_1+{\theta }_{{\Omega }_k}\left(k-1\right),& k=1,2,3\dots N\end{array}$ 11a

$\begin{array}{cc}{M}_{k0}=M_{10}-\left({\Omega }_k+\omega \right),& k=1,2,3\dots N\end{array}$ 11b

Next, a non-identical constellation pattern for GSO satellites is made possible by adding Δє from the desired ground track interval and the initial mean anomaly. Thus, the phasing rule is given by the following:

$\begin{array}{cc}{M}_{k0}=M_{10}-\left({\Omega }_k+\omega \right)+\Delta \epsilon ,& k=1,2,3\dots N\end{array}$ 12

It is well known that the trajectory of a GSO satellite with an inclination has a figure-eight shape. Therefore, Equation (12) can place figure-eight shapes in non-identical patterns on the longitude in intervals of Δє on the Earth’s surface.

As a result, by substituting the phasing rule from Equation (11) into Equation (7), we obtain the following closed-form solution for which N_s satellites have the same GSO pattern:

$x_k=\frac{r}{2}\left(\left(1+\cos i\right)\cos \left(v_k-{\omega }_{\oplus }t+\omega +{\Omega }_k\right)+\left(1-\cos i\right)\cos \left(v_k+{\omega }_{\oplus }t+\omega -{\Omega }_k\right)\right)$ 13a

$y_k=\frac{r}{2}\left(\left(1+\cos i\right)\sin \left(v_k-{\omega }_{\oplus }t+\omega +{\Omega }_k\right)-\left(1-\cos i\right)\sin \left(v_k+{\omega }_{\oplus }t+\omega -{\Omega }_k\right)\right)$ 13b

$z_k=r\sin i\sin \left(v_k+\omega \right)$ 13c

where k = 1, 2, 3, …, N_s and the effect of the short-term perturbation by J₂ on the GSO is assumed to be small and is thus ignored. The closed-form solution in Equation (13) can directly describe the ground tracks of N_s satellites on the Earth’s surface; in particular, it is observed that the three orbital elements Ω, M₀, and ω play an important role in defining the GSO constellation configuration.

By evaluating the effect on the GSO ground track for each orbital element in Equation (13), it can be seen that the change in the shape of the GSO ground track is determined only by ω, the phase value of the z axis. Figure 3(a) shows the change in ground tracks as a function of ω values from 0º to 30º for seven GSO satellites with the same orbital elements. This change affects the relative dynamics between the ground station on the Earth’s surface and the satellite, which is closely related to the performance (GDOP, navigation errors, etc.) of the satellite constellation. The next critical orbital element is Ω, which controls the satellite spacing of the GSO satellite constellation. In the case of GDOP, which generally measures GNSS performance, good results are obtained with equally spaced satellite placement. However, the performance of the satellite constellation varies depending on the location of the ground station and the spacing between satellites. Figure 3(b) shows equally spaced and non-equally spaced deployments for four satellites with i = 30º and ω = 45º, and the GDOP performance of the two satellite constellations is shown in Table 2. In some cases, the GDOP performance of the four satellites is better in non-equally spaced placements than in equally spaced placements for a specific ground station, as shown in Table 2. This result emphasizes that there are various GSO satellite constellation design cases, that evenly spaced constellations do not always have better performance than non-equally spaced constellations, and that better satellite constellation geometries can be found through simulations with the proposed solution. Thus far, no dedicated solution related to GSO satellite constellation design has been presented in the literature; accordingly, Equations (11)–(13) have academic novelty, providing a means to easily analyze various constellation designs and verify overall GDOP performance by simply substituting orbital elements into the closed-form solution.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Change in GSO constellation according to ω and Ω values (a) Change in GSO ground track shape (b) Equally spaced (circles) and non-equally spaced (squares) placement

4 GDOP PERFORMANCE

This section examines the theoretical background for evaluating the performance of various GSO satellite constellations via the GDOP and the algorithm for applying the closed-form solution. The purpose of satellite navigation is to provide users with the time and actual position. However, the expected location and time differ from the user’s actual position and time owing to various factors, including atmospheric errors and the geometry between the receiver and satellite. Navigation system errors can be broadly classified into two types. The first type is the statistical user equivalent range error (UERE) that occurs during the distance measurement process. The second error type is the GDOP error caused by the geometric arrangement of the navigation satellites. The GDOP reduces the accuracy of the time and position calculated at the receiver end. In contrast, the UERE includes errors due to satellites, receiver clock biases, atmospheric propagation, and multipaths. Here, errors caused by the geometric characteristics of the satellite receiver can be minimized by maintaining a good geometric arrangement between the satellite and the user (Sharp et al., 2009; Langley, 1999; El-Hakim, 2023).

The GDOP can be used to determine how to deploy satellites to provide accurate positioning and time services to users. Figure 4 shows a flowchart for obtaining the GDOP performance for various constellation configurations by applying the closed-form solution proposed in this study. First, the orbital elements of the first satellite are set and the total number of satellites (N_s) is determined. The GDOP performance for a specific ground station is calculated based on considerations such as the location of the ground station and the spacing $\left({\theta }_{{\Omega }_k}\right)$ between satellites. Under the same orbital element conditions, the closed-form solution simply calculates the position of each satellite distributed by the phasing rule of Ω_k and M_k0 while minimizing the computational load. Here, we examine the steps for calculating the GDOP and obtain the unit vector for the position vector between the satellite and desired receiver using Equation (13):

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Implementation of the GDOP performance algorithm for assessing a satellite constellation via the application of a closed-form solution

$\left(\frac{x_k-x}{{\rho }_k},\frac{y_k-y}{{\rho }_k},\frac{z_k-z}{{\rho }_k}\right)$ 14

where:

${\rho }_k=\sqrt{{(x_k-x)}^2+{(y_k-y)}^2+{(z_k-z)}^2}$ 15

In Equation (15), x, y, and z are the ECEF coordinates of the user receiver, and x_k, y_k, and z_k are the ECEF coordinates of the k-th satellite. Matrix A for the pseudorange measurement residual equations is formulated as follows:

$\left(\begin{array}{cccc}\frac{x_1-x}{{\rho }_1}& \frac{y_1-y}{{\rho }_1}& \frac{z_1-z}{{\rho }_1}& -1\\ \frac{x_2-x}{{\rho }_2}& \frac{y_2-y}{{\rho }_2}& \frac{z_2-z}{{\rho }_2}& -1\\ \frac{x_3-x}{{\rho }_3}& \frac{y_3-y}{{\rho }_3}& \frac{z_3-z}{{\rho }_3}& -1\\ ⋮& ⋮& ⋮& ⋮\end{array}\right)$ 16

If the element in the fourth column of matrix A is –1, the time dilution of precision (TDOP) is calculated appropriately. Based on matrix A, the covariance matrix Q can be obtained, which is used to determine the GDOP, position dilution of precision (PDOP), and TDOP as follows:

$Q={\left(A^{T}A\right)}^{-1}$ 17

and:

$\left(\begin{array}{cccc}{\sigma }_x^2& {\sigma }_{xy}& {\sigma }_{xz}& {\sigma }_{xt}\\ {\sigma }_{xy}& {\sigma }_y^2& {\sigma }_{yz}& {\sigma }_{yt}\\ {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_z^2& {\sigma }_{zt}\\ {\sigma }_{xt}& {\sigma }_{yt}& {\sigma }_{zt}& {\sigma }_t^2\end{array}\right)$ 18

PDOP, TDOP, and GDOP are given as follows:

$PDOP=\sqrt{{\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2}$ 19

$TDOP=\sqrt{{\sigma }_t^2}$ 20

$GDOP=\sqrt{tr Q}$ 21

The GDOP represents the comprehensive uncertainty of time and position due to the satellite constellation geometry and is used to determine the performance index of the space segment design of the satellite system. If the GDOP value is small, the uncertainty of the user’s position is also small. The GDOP value changes as the satellite position changes with respect to the ground location. Therefore, the overall mean value of the GDOP can be used as a measure of constellation performance.

5 SIMULATION RESULTS

In the previous section, a closed-form solution and satellite phasing rules for designing a ground track and satellite constellation were proposed based on the geometric principles of Ptolemaic deferent–epicycle and spring–mass systems. In this section, the point and regional coverage obtained for different GDOP values are evaluated to analyze the performance of GSO constellations.

5.1 Point Coverage

We analyzed the performance of two satellite constellation patterns using the analytical solution and phasing rules presented above. A performance analysis was performed by comparing the GDOP values. The GSO constellations were designed such that seven satellites (k = 7) have identical ground tracks, according to Equation (13). Table 3 lists the specifications of the two types of GSO constellations, which have different eccentricities/arguments of perigee. The spacing of the ascending nodes is 51.4º.

View this table:

• View inline
• View popup

TABLE 3

Specifications for Two Types of GSO Constellations

As shown in Figure 5, the GSO I constellation has seven satellites evenly distributed in a circular orbit and a symmetrical structure with a ground station on the equator. The GDOP values present sinusoidal curve characteristics. The GSO II constellation shown in Figure 6 has an eccentricity of 0.1 and an argument of perigee of 90º. Therefore, the figure-eight-shaped crossing point descends toward the south. The distribution of GDOP values exhibits an irregular curve for the equatorial ground station. However, the mean GDOP values of the two GSO constellations were 2–4 for a ground station located at the equator, indicating excellent performance.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Constellation and GDOP performance for a circular GSO (a) Circular GSO I pattern (b) Time vs. GDOP values

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Constellation and GDOP performance for a noncircular GSO (a) Noncircular GSO II pattern (b) Time vs. GDOP values

5.2 Regional Coverage

This section compares the performance of regional GDOP values as the ground station location is varied, aiming to analyze the local coverage for various GSO constellation patterns. The GSO constellation patterns are divided into identical patterns, in which all satellites have the same ground track, and non-identical patterns, in which the satellites have different ground tracks. In this study, for convenience, the number of satellites in an identical pattern is denoted as σ_i, and the number of satellites in a non-identical pattern is denoted as σ_n. For the regional coverage analysis, we determined the GDOP value while constantly changing the ground station latitude based on the longitude of the GSO crossing point.

Two types of GSO constellations with seven satellites were compared and analyzed for the case of σ_i = 7/σ_n = 0, which consists of only identical patterns, and σ_i = 3 × 2/σ_n = 1, which consists of two patterns combined. Table 4 presents the specifications of the satellite constellations used for the comparative analysis.

View this table:

• View inline
• View popup

TABLE 4

Specifications for Two Types of GSO Constellations

Figure 7 presents the results for case I, as described in Table 4, with an identical circular GSO constellation. The circularly identical GSO pattern shown in Figure 7(a) implies that the GDOP values of the symmetrically distributed ground stations are symmetric about the equator. The best GDOP value is observed for the location at the equator; as the distance from the equator increases, the GDOP performance becomes worse. Figure 7(b) illustrates three GSO patterns, with the figure-eight in the middle being the trajectory of one satellite and the figure-eights on both sides consisting of three satellites. Compared with the results shown in Figure 7(a), the GDOP performance is improved at the Equator, but the overall deviation from the ground station is notably large. This result indicates that, in the case of a circular GSO, maintaining an identical constellation pattern for all satellites helps improve the GDOP performance.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

Regional GDOP distribution for a circular pattern (a) Only circular orbit pattern (b) Combined circular orbit pattern

Figure 8 presents the results for case II, as described in Table 4, with a non-identical GSO constellation. For the non-identical pattern in Figure 8(a), the crossing point of the figure-eight moves to the southern hemisphere owing to the eccentricity and the argument of perigee. The variation in GDOP values is markedly large, especially in the northern hemisphere, where there are no crossing points, and the GDOP performance is very poor. Figure 8(b) displays a constellation with combined patterns of a GSO in an elliptic orbit. Although the GDOP performance is marginally better than that shown in Figure 8(a), it is still worse than that of the circular pattern. The QZSS and KPS are representative examples that leverage a noncircular orbit pattern. For the constellation patterns of these navigation satellite systems, engineers can improve the performance of the GDOP by moving the crossing point of the elliptical orbit toward the region of each respective country in the northern hemisphere.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Constellation and GDOP performance for a noncircular GSO (a) Only elliptical orbit pattern (b) Combined elliptical orbit pattern

Figure 9 presents the results from Figures 7 and 8 according to the latitude of the ground station. Overall, for latitudes from –30º to 30º, the GDOP value remains at 4 or less, showing excellent performance. For all constellation patterns, the ground station at the equator exhibits the best GDOP performance. In summary, the identical circular pattern is stable and exhibits the best performance in regional areas compared with other patterns. Even for circular orbits, the GDOP performance gap between the northern and southern hemispheres widens as it shifts toward a non-identical pattern. In addition, the elliptical orbit pattern moves the crossing point of the figure-eight to the southern or northern hemisphere, depending on the value of the argument of perigee, which causes the GDOP performance to be biased toward one region.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

Regional GDOP values as a function of the ground station latitude

6 CONCLUSIONS

The dedicated tool proposed in this study worked effectively and obtained meaningful results for analyzing and evaluating the GDOP performance of point and regional coverage of various GSO satellite constellations. Specifically, the point coverage was used to compare the GDOP performance of circular and noncircular GSO patterns for seven evenly distributed satellites. Here, the GDOP performance exhibited average GDOP values of 2–4 for ground stations located at the equator, irrespective of the GSO constellation. To evaluate the regional coverage performance, the regional GDOP distribution was analyzed for various GSO patterns. The ground station located on the equator displayed the best GDOP performance. Moreover, for a satellite group with a given number of satellites, the ideal constellation pattern showed better performance than that with multiple ground track patterns. Additionally, the crossing point location on the figure-eight trajectory contributes to the improved GDOP performance. The ideal crossing point location should correspond to a ground station location at any point on the equator.

DATA AVAILABILITY

The authors declare that all data relevant to this study are available within the paper.

DECLARATION OF COMPETING INTERESTS

None.