Cell Area-Based Method for Analyzing the Coverage Capacity of Satellite Constellations

Abstract

Constellation-to-ground coverage analysis is an important problem in practical satellite applications. The classical net point method is one of the most commonly used algorithms in resolving this problem, indicating that the computation efficiency significantly depends on the high-precision requirement. On this basis, an improved cell area-based method is proposed in this paper, in which a cell is used as the basic analytical unit. By calculating the accuracy area of a cell that is partly contained by the ground region or partly covered by the constellation, the accurate coverage area can be obtained accordingly. Experiments simulating different types of coverage problems are conducted, and the results reveal the correctness and high efficiency of the proposed analytical method.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

With the rapid development of science and technology, satellite technology plays an increasingly important role in mobile communications [1], remote sensing detection [2], military reconnaissance [3], disaster monitoring [4], navigation [5], and other fields. With the deepening of the application field, the complexity of space missions has gradually increased, and an increasing number of space missions can no longer be completed by a single satellite. Constellations composed of many satellites according to a specific space configuration have become an inevitable development trend.

In satellite applications, ground coverage is a very important component. Constellation-to-ground coverage is accomplished by continuously or intermittently serving from various types of sensors loaded in constellation satellites. The coverage capacity of a constellation (i.e., the coverage rate and the coverage quality) over a long time period will be important issues in constellation applications. This kind of problem is called “the constellation-to-ground coverage problem.” The problem model ignores the ground features and considers the whole Earth’s surface to be spherical or ellipsoid without fluctuations. In most satellite applications, the model can support engineering needs.

The constellation-to-ground coverage problem is complicated. To date, many scholars have carried out relevant research on this topic. References [6, 7] used probabilistic statistical models to analyze the problem and obtained the change rule of the spatial and temporal coverage characteristics in the sense of probability. In Reference [8], a judgment theorem for global complete coverage was proposed based on spherical geometry relations. In References [9, 10], a performance analysis method based on the transformation group was proposed that uses crystallographic point group theory to describe it. In Reference [11], the numerical fitting method was used to fit the coverage characteristics to obtain the fitting function of the constellation-to-ground coverage problem. Reference [12] transformed the coverage problem into a differential equation and solved the function by a numerical integration method. In References [13, 14], the ground coverage problem of the Flower constellation was described by using a matrix, and the spatial configuration of the constellation was isomorphic to a second- or third-order square matrix to realize the design and performance analysis of the constellation. In Reference [15], some feature points of the coverage process were obtained, and the coverage performance of the whole region was characterized by the performance of these feature geometries. In Reference [16], an interpolation algorithm was used to calculate the cumulative coverage of agile satellites over a period of time. In Reference [17], the coverage of satellites in a period of time is regarded as a family of curves, the envelope of the family of curves is solved, and the correlation between the envelope and the boundary of the coverage is obtained. All these methods construct a mathematical model of the coverage problem and then analyze the mathematical characteristics of the problem to represent the coverage problem. These algorithms have some common characteristics. First, these methods generally only model the characteristics of a certain aspect of the constellation-to-ground coverage problem. Second, in the process of mathematical modeling, some approximate processing is often needed. Therefore, the obtained result is mostly an approximate result and is often a qualitative result. Third, this kind of method adopts mathematical modeling, so the result has some resolvable characteristics.

The net point method is a traditional method for solving the constellation-to-ground coverage problem [18–20]. Almost all kinds of coverage problems can be solved by the net point method. The streets-of-coverage method is a common method for solving satellite coverage problems [21–23]. It approximates the coverage range of the satellite in a period of time as a coverage street to analyze the performance of the constellation in a period of time. The spherical triangulation method divides the Earth’s surface into a group of spherical Delaunay triangles with satellites as nodes to perform performance analysis [24, 25]. It has been widely used in constellation performance analysis. The spherical polygon division method [26] is based on the dual vino diagram of the spherical triangulation method, and it can be used to analyze arbitrarily shaped regions. The latitude band method [27] divides the target region into a group of latitude bands for analysis, and the algorithm is efficient. The feature area analysis method divides the two-dimensional sphere into a series of regions based on the trajectories of the satellites in the constellation [28, 29]. These types of methods involve the decomposition of the regional target and have several common characteristics. First, they divide the sphere or region into a set of specific graphs by a division strategy. Second, in general, the coverage problem can be solved for a variety of characteristics, and the accuracy of the calculated results is high. Third, these methods are more numerical and can obtain more numerical information but less analytical information than other methods.

In this paper, an improved algorithm is proposed, which is called the cell area analytical method based on the net point method. On the basis of keeping the basic operation of the net point method, the cell area analytical method takes the cell as the core of the whole algorithm and calculates the precise coverage area of each cell by the spherical geometric relation. Then, the characteristics of the satellite’s ground coverage are obtained.

The basic organizational structure of this paper is as follows: the first chapter is the introduction, and the second chapter introduces the net point method and the analysis of its existing problems. In the third chapter, an improved cell method is proposed to analyze the flow of the algorithm. The fourth chapter is the simulation experiment, and the last chapter is the conclusion.

2. Basis of the Net Point Method

The basic steps of the net point method are shown in Figure 1. The details of this process are as follows.

Step 1.

For a certain target region $R$, obtain its minimum latitude and longitude encirclement. Then, a group of points in the minimum latitude and longitude encirclement of $R$ is sampled according to certain rules, such as the system sampling method. The set of the sampling points is called net points and is denoted by $D_R$. This basic step is shown in the second subfigure in Figure 1, in which the black points represent all of the net points.

Step 2.

For each point $P\in D_R$, determine whether it is in $R$. If a point $P$ is not in $R$, it needs to be removed from $D_R$. Only the points inside region $R$ are retained. This step is shown in the third subfigure in Figure 1, in which the red and blue points are the inside and outside net points of the region, respectively.

Step 3.

Let the coverage region of the constellation be ${\Lambda }_C$. For each point $P\in D_R$, determine whether it is in ${\Lambda }_C$. Count the number of points in $D_R$ that is also in ${\Lambda }_C$, denoted by $N$. This step is shown in the fourth subfigure in Figure 1. The green points stand for the points in ${\Lambda }_C$. The coverage rate of the constellation to the region is calculated, which is denoted by $\eta$; then $$\begin{array}{}\text{(1)}& \eta =\frac{N}{\mathrm{card}{D}_R},\end{array}$$where $\mathrm{card}{D}_R$ represents the number of elements in $D_R$.

[figure omitted; refer to PDF]

In the classical net point division method, the coverage state of the net point is regarded as that of the corresponding net grid. That is, if a satellite constellation can cover the net point, it is considered that the constellation can cover the whole corresponding net grid; otherwise, it is considered that the constellation has no coverage for the corresponding net grid. This strategy is called the 0-1 judgment strategy. According to the principle of the algorithm and the statistics of the experimental results, if the precision of the results increases $N$-fold, the computing time will increase $N^2$-fold. Thus, this method is usually unsuitable for very large regions under high accuracy requirements [26].

From a numerical simulation point of view, the net point division method chooses the net point as the basic unit of analysis. However, the generated net points are just a set of single point on the sphere and have no any spatial structure, making the low-precision calculation results cannot provide any guidance for the high-precision calculation results. Hence, it is difficult to calculate iteratively.

This method exhibits a low computational efficiency and little reliability of the results. Thus, an improved method called the cell area analytical method (CAAM) is proposed in this paper. The improvement of the CAAM mainly focuses on two points: one is replacing the net point with a cell as the basic analytical unit, and the other is changing the 0-1 judgment strategy into the real intersection area calculation strategy.

3. Cell Area Analytical Method

In this section, we propose a new method for efficiently and exactly calculating the coverage area for the constellation-to-ground coverage problem, which is called the cell area analytical method.

3.1. Calculation of the Area of the Spherical Arch Region

For a spherical circle, the inner region is called a spherical disk. That is, the spherical circle is the boundary of the spherical disk. Denote the spherical disk by $\Theta$ and denote the spherical circle by $\partial \Theta$.

As shown in Figure 2, the circle in the figure represents a spherical disk, and the red line segment represents a great arc segment on the sphere. The minor region where a spherical disk is separated by a spherical great arc is called a spherical arch region, which corresponds to the yellow area in Figure 2. In this section, we derive the formula for the area of the spherical arch region.

[figure omitted; refer to PDF]

For the spherical circle $\partial \Theta$, the spherical radius, which means the spherical distance from the center of the spherical circle to a random point on the spherical circle, is written as $r$. The symbol $mU$ represents the area of a spherical region $U$. The radius of the Earth is unitized. The area of the spherical disk is as follows: $$\begin{array}{}\text{(2)}& m\Theta =2\pi \cdot 1-\cos r.\end{array}$$

Denote the longitude and latitude of a point $P$ on the Earth’s surface as ${\theta }_P,{\varphi }_P$ and $Q$ as ${\theta }_Q,{\varphi }_Q$. The spherical distance between the two points, denoted by $d$, can be expressed by the following formula: $$\begin{array}{}\text{(3)}& d={\cos }^{-1}\cos {\varphi }_P\cos {\varphi }_Q\cos {\theta }_P-{\theta }_Q+\sin {\varphi }_P\sin {\varphi }_Q.\end{array}$$

The relationship between $d$ and its corresponding spherical central angle $\alpha$ can be expressed as follows: $$\begin{array}{}\text{(4)}& \alpha =2\cdot {\sin }^{-1}\frac{\sin d/2}{\sin r}.\end{array}$$

The spherical sector with central angle $\alpha$ in Figure 2 is written as ${\nabla }_{OPQ}$. The area of ${\nabla }_{OPQ}$ is as follows: $$\begin{array}{}\text{(5)}& m{\nabla }_{OPQ}=1-\cos r\cdot \alpha .\end{array}$$

The spherical triangle with vertices $O,P,Q$ is denoted by ${\Delta }_{OPQ}$. The area of ${\Delta }_{OPQ}$ is $$\begin{array}{}\text{(6)}& m{\Delta }_{OPQ}=2{\tan }^{-1}\frac{1}{\cos r\cdot \tan \alpha /2}+\alpha -\pi .\end{array}$$

The spherical arch region in Figure 2 is written as $A_{PQ}\Theta$. The area of this spherical arch region is $$\begin{array}{}\text{(7)}& m{A}_{PQ}\Theta =m{\nabla }_{OPQ}-m{\Delta }_{OPQ}.\end{array}$$

Therefore, we can obtain the area of the spherical arch only if the radius $r$ and the two points $P$ and $Q$ are known.

3.2. Definition of the Cell and Its Basic Properties

The definition of cell in this paper is as follows:

Definition 1 Cell.

A cell is the spherical region surrounded by two latitude lines and two longitude lines. Additionally, it has the same length of longitude interval and latitude interval.

The length of the longitude interval is called the cell width.

Denote a cell by $G_k$. The relationship between a cell $G_k$ and a ground region $U$, which is called the cell-region relation and is written as $\Pi G_k,U$, can be divided into three categories as follows:

(i) $G_k$ is inside of $U$, which is written as $G_k\subset U$

(ii) $G_k$ is outside of $U$, which is written as $G_k\subset U^C$, where $U^C$ represents the complementary set of $U$ with respect to the Earth’s surface

(iii) $G_k$ is partly inside of $U$, which is written as $G_k\mid =U$, where the symbol “$\mid =$” represents a set that is partly contained by another set

Suppose the longitude range and latitude range of $G_k$ are ${\theta }_k^L,{\theta }_k^U$ and ${\varphi }_k^L,{\varphi }_k^U$, respectively. Normalize the radius of the Earth as a unit length; then, the area of the cell $G_k$ is expressed as follows: $$\begin{array}{}\text{(8)}& m{G}_k=\sin {\varphi }_k^U-\sin {\varphi }_k^L{\theta }_k^U-{\theta }_k^L.\end{array}$$

3.3. Intersection Area of a Spherical Disk and a Cell

A spherical disk is defined as the spherical region surrounded by a spherical circle, and a spherical disk with index $j$ is written as ${\Theta }_j$.

The intersection area between a cell $G_k$ and a spherical disk $\Theta$ is denoted by $m{G}_k\cap \Theta$. The calculation of $m{G}_k\cap \Theta$ is a core algorithm in this paper. According to the cell-region relation, it can be classified into four categories as follows:

(i) If $\Theta \subset G_k$, then $m{G}_k\cap \Theta =m\Theta$. However, in general, the size of a cell is much smaller than that of a spherical disk. Therefore, this category need not be considered

(ii) If $G_k\subset \Theta$, then $m{G}_k\cap \Theta =m{G}_k$, which can be computed by Equation (8)

(iii) If $G_k\subset {\Theta }^C$, then $m{G}_k\cap \Theta =0$

(iv) If $G_k\mid =\Theta$, the computation process of $m{G}_k\cap \Theta$ is very complex, which will be discussed in detail in this section

For $G_k\mid =\Theta$, except for very extreme cases, $G_k$ and $\Theta$ have two intersection points in general. We analyze a typical case for calculating $m{G}_k\cap \Theta$ when $G_k\mid =\Theta$.

In Figure 3, We use a diagram to show one of the relations between a spherical circle and a cell. All of the elements in the diagram represent relations on the sphere. The square represents a cell, and the circle represents a spherical disk. $P$ and $Q$ are the two intersection points. The longitude and latitude coordinates of $P$ and $Q$ can be computed by spherical geometry [27].

[figure omitted; refer to PDF]

According to the characteristics of the intersection region, we compute the intersection area between the spherical disk and the cell by dividing this region into subregions. In Figure 3, the intersection area is divided into three regions, namely, the blue region $S_1$, the yellow region $S_2$, and the green region $S_3$.

The blue region $S_1$ is a spherical longitude and latitude box. Supposing the value of the longitude interval of $Q$ and $M$ is $\Delta \theta$, then $$\begin{array}{}\text{(9)}& m{S}_1=\sin {\varphi }_k^U-\sin {\varphi }_k^L\cdot \Delta \theta .\end{array}$$

The yellow region $S_2$ is a typical spherical arch region. If the given coordinates of $P$ and $Q$, its area can be obtained according to Equations (2)–(7).

Since the upper and the lower boundaries of the cell are latitude lines, which are spherical small arcs in general, $S_3$ is not a spherical triangle. In obtaining the area of $S_3$, the spherical triangle with the same vertexes of $S_3$ should be firstly constructed, and then, the extra area is subtracted. We make a great circle by the two points $P$ and $N$. The region surrounded by the latitude line and the great circle by the two points $P$ and $N$ is denoted by $S_4$. The spherical triangle with the same vertexes of $S_3$ corresponds to the combination of $S_3$ and $S_4$. The area of the generated spherical triangle can be obtained by the given coordinates $P$, $N$, $Q$. The radius of this generated spherical arch region $S_4$ is $\pi /2-{\varphi }_k^U$, and the center angle is ${\theta }_N-{\theta }_P$. Thus, the area of $S_4$ can be calculated.

Then, $m{G}_k\cap \Theta$ can be expressed as follows: $$\begin{array}{}\text{(10)}& m{G}_k\cap \Theta =m{S}_1+m{S}_2+m{\Delta }_{PQN}-m{S}_4.\end{array}$$

Obviously, there are many other cases in which a spherical disk intersects a cell. There exists a total of 16 cases by analysis. Figure 4 shows the division of the intersection region under three other circumstances. In any case, the intersection region can be treated in a similar way by dividing them into one or sometimes two spherical rectangular regions (such as the blue region in Figures 3 and 4), a spherical arched region (such as the yellow region in Figures 3 and 4), a spherical triangle (such as the combination region or the difference region of the green and red region in Figures 3 and 4), and an arched region between the latitude zone and the spherical great circle (such as the red region in Figures 3 and 4). In different cases, the calculation process is slightly different, but the overall method is similar.

[figures omitted; refer to PDF]

3.4. The Intersection Area of a Ground Region and a Cell

Let $U$ be a ground region with an arbitrary shape. The edge of $U$, which can be written as $\partial U$, consists of or can be approximated as a group of great circular arcs and small circular arcs. The great circular arcs and small circular arcs are collectively called spherical circular arcs. A spherical circular arc with index $j$ is written as $A_j$; then, $A_j\in \partial U$.

Additionally, by the conclusion in the former section, there are three types of relationships between $G_k$ and $U$: $G_k\subset U$, $G_k\subset U^C$, and $G_k\mid =U$. If $G_k\mid =U$, then $G_k$ can be classified into two categories with respect to $U$, which are called type-I and type-II.

(i) Type-I: only one circular arc in $\partial U$ crosses the cell

(ii) Type-II: more than one circular arc in $\partial U$ crosses the cell

According to the spherical geometry, a random spherical circular arc $A_j$ specifying the inner side corresponds to a spherical disk, denoted by ${\Theta }_j$. For $U$, if only one circular arc in $\partial U$ crosses cell $G_k$, the intersection area of $U$ and $G_k$ is equivalent to that of ${\Theta }_j$ and $G_k$. $$\begin{array}{}\text{(11)}& mU\cap G_k=m{\Theta }_j\cap G_k.\end{array}$$

We can use Equation (11) and its variant formulas to compute $mU\cap G_k$.

However, if more than one circular arc in $\partial U$ crosses cell $G_k$, then the intersection area of $U$ and $G_k$ cannot be simply computed. To compute the precision result of $mU\cap G_k$, cell $G_k$ should be partitioned into a group of subcells.

If $G_k$ is divided into $N$ subcells, the $j$th subcell is written as $G_{k,j}$, that is, $$\begin{array}{}\text{(12)}& G_k=\bigcup _{j=1}^N{G}_{k,j}.\end{array}$$

Then, $$\begin{array}{}\text{(13)}& mU\cap G_k=\sum _{j=1}^{N}mU\cap G_{k,j}.\end{array}$$

$G_{k,j}$ is smaller than $G_k$; in general, part of the subcells in $G_{k,j}$ belong to type-I, and the others belong to type-II. The intersection area between $U$ and subcells that belong to type-I can be computed, and subcells that belong to type-II should repeat this process until the calculation results meet the precision requirements.

By this process, we can obtain the intersection area between the grid and the region with arbitrary precision.

3.5. Intersection Area of the Constellation and Cell

A constellation is a set of satellites that is denoted by $C$. The instantaneous coverage of satellite $S$ to the ground is denoted by ${\Lambda }_S$, and the instantaneous coverage of constellation $C$ is denoted by ${\Lambda }_C$.

Then, $$\begin{array}{}\text{(14)}& {\Lambda }_C=\bigcup _{S\in C}{\Lambda }_S.\end{array}$$

For a ${\Lambda }_S$ at a certain time, it corresponds to a ground region. By Subsection 3.4, we can compute the result of $m{G}_k\cap {\Lambda }_S$. If $\partial {\Lambda }_S$ consists of more than one spherical arc, $G_k$ may belong to type-II with respect to ${\Lambda }_S$. We can then compute the result of $m{G}_k\cap {\Lambda }_S$.

For the computation of $m{G}_k\cap {\Lambda }_C$, we also need to consider the case where $G_k$ is covered by more than one satellite simultaneously.

Consider the situation shown in Figure 5. The two circles represent the instantaneous coverage of two satellites, which intersect with each other, and the square grids represent the cells. According to the relations between the cell and the instantaneous coverage of the constellation, the cells can be divided into four categories: white cells, red cells, blue cells, and green cells, as shown in Figure 5. A white cell is a cell that is not covered by any satellite. A red cell is a cell that at least one satellite in the constellation can completely cover. A blue cell is a cell partially covered by only one satellite in the constellation, while a green cell is a cell partially covered by more than one satellite, and each satellite partially covers the cell.

(i) For white cells, $m{G}_k\cap {\Lambda }_C=0$

(ii) For red cells, $m{G}_k\cap {\Lambda }_C=m{G}_k$

(iii) For blue cells, $G_k$ is only partly covered by one satellite $S$; then, $m{G}_k\cap {\Lambda }_C=m{G}_k\cap {\Lambda }_S$

(iv) For green cells, $G_k$ belongs to type-II with respect to ${\Lambda }_C$. We need to divide $G_k$ into subcells to obtain a precision result

[figure omitted; refer to PDF]

The computation results indicate that the boundary of coverage rate converges to the same value. Figures 7 and 9 show the rates of convergence of the two targets are completely different. When the cell precision is increased by 10 times, the difference between the upper and lower bounds decreases by 100 times for target-1, while it only decreases by approximately 10 times for target-2. The reason is that for target-2, the large number of boundaries causes many cells to be classified as type-II.

Figures 10 and 11 represent the computing time of the proposed cell area analytical method and the traditional net point division method with an increase of the partition precision.

[figure omitted; refer to PDF]

Further results show that when the cell precision is low, the cell area analytical method needs to take longer computing time than the traditional net point method. That is because the cell area analytical method takes much longer time to compute a single cell than the numerical method to compute a signal net point. However, with the high-precision calculation requirement based on the increase of the cell partition precision, the computing time of the traditional net point method increases sharply and quickly exceeds than that of the cell analytical method. Additionally, in Figure 10, the computing time of the coverage rate with the cell area analytical method increases slightly. It is because only a very few numbers of cells belong to type-II, and the number of cells that need to be calculated increases slightly with an increase of the high accuracy requirement. Figure 11 shows that for an arbitrarily shaped polygon bounded by many boundaries, the computing time of the cell area analytical method increases in linear form. The computing time of the traditional net point method increases in the form of the square with the increases of the cell precision both in target-1 and in target-2.

Figures 12 and 13 represent the simulation of the coverage region for the two targets. The circles in Figures 12 and 13 represent the instantaneous coverage region of the satellites in the constellation. The red cells represent the cells in which $G_k\in R$ and $G_k\in {\Lambda }_C$. The blue cells represent the type-I cells, while the green cells represent the type-II cells. Of course, because all the green cells need to be divided, in the final results, all green cells have been divided into very small cells.

[figure omitted; refer to PDF]

We can see that the cell area analytical method can solve all types of coverage regions, including regular regions and irregular regions. For ground regions with fewer boundaries, the efficiency of the algorithm is extremely high, but for ground regions with a high number of boundaries, the efficiency of the algorithm is slightly lower than that of regions with fewer boundaries but still very high compared with the traditional net point method.

In Figures 12 and 13, we can see that there are many cells classified as type-II because more than one boundary crosses the cell. To obtain the exact area of the cell belonging to the region, the cell must divide into subcells. This is why the efficiency of the algorithm is slightly lower than that of regions with fewer boundaries.

5. Conclusion

In this paper, based on the net point method, an improved method is proposed, which we call the cell area analytical method. Compared with the traditional net point method, the improved cell area analytical method can accurately compute the coverage area of the partial coverage cell. By a simulation experiment, we can see that the method has high computational efficiency, especially for regions with fewer boundaries.

However, the cell area analytical method still has some disadvantages; that is, it can only solve the instantaneous coverage problem, while the net point method can also solve the accumulative coverage problem and continuous coverage problem.

Additionally, this algorithm has room for improvement. That is, for ground regions with many boundaries, the spherical polygon area formula can be used to calculate the exact area of the intersecting regions, but we need to consider many complicated situations.

Acknowledgments

This work is supported by the National Key R & D Program of China under Grant No. 2016YFB0501001, the National Natural Science Foundation of China under Grant No. 62006214, the 13th Five-year Pre-research Project of Civil Aerospace in China, the China Postdoctoral Science Foundation under Grant No. 2019TQ0291, the Aeronautical Science Fund under Grant No. 2018ZCZ2002, the Opening Fund of Key Laboratory of Geological Survey and Evaluation of Ministry of Education under Grant No. GLAB2019ZR04, and the Open Research Project of the Hubei Key Laboratory of Intelligent Geo-Information Processing (KLIGIP-2018B06 and KLIGIP-2019B07).

References

[1] K. Ohata, K. Kobayashi, K. Nakahira, M. Ueba, "Broadband and scalable mobile satellite communication system for future access networks," Acta Astronautica, vol. 57 no. 2-8, pp. 239-249, DOI: 10.1016/j.actaastro.2005.03.028, 2005.

[2] J. E. Patino, J. C. Duque, "A review of regional science applications of satellite remote sensing in urban settings," Computers, Environment and Urban Systems, vol. 37,DOI: 10.1016/j.compenvurbsys.2012.06.003, 2013.

[3] B. Penné, C. Tobehn, M. Kassebom, H. Lübberstedt, "The Earth observation program at OHB-system," Acta Astronautica, vol. 56 no. 1-2, pp. 17-23, DOI: 10.1016/j.actaastro.2004.09.045, 2005.

[4] G. Santilli, C. Vendittozzi, C. Cappelletti, S. Battistini, P. Gessini, "CubeSat constellations for disaster management in remote areas," Acta Astronautica, vol. 145, pp. 11-17, DOI: 10.1016/j.actaastro.2017.12.050, 2018.

[5] B. Parkinson, "The origins, evolution and future of satellite navigation," 34th Aerospace Sciences Meeting and Exhibit, 1996.

[6] R. D. Luders, "Satellite networks for continuous zonal coverage," ARS Journal, vol. 31 no. 2, pp. 179-184, DOI: 10.2514/8.5422, 1961.

[7] Y. Seyedi, S. M. Safavi, "On the analysis of random coverage time in mobile LEO satellite communications," IEEE Communications Letters, vol. 16 no. 5, pp. 612-615, DOI: 10.1109/LCOMM.2012.031912.112323, 2012.

[8] J. E. Draim, "Three- and four-satellite continuous-coverage constellations," Journal of Guidance, Control, and Dynamics, vol. 8 no. 6, pp. 725-730, DOI: 10.2514/3.20047, 1985.

[9] G. V. Mozhaev, "The problem of the continuous earth coverage and the kinematically regular satellite networks," Cosmic Research, vol. 11, 1973.

[10] G. V. Mozhaev, "Capabilities of kinematically regular satellite systems with symmetry groups of the second type in the problem of continuous single coverage of the earth," Cosmic Research, vol. 43 no. 3, pp. 205-212, DOI: 10.1007/s10604-005-0036-2, 2005.

[11] A. Modiri, L. Mohammady, "Mathematical prediction of sun-synchronous polar LEO satellite visions for earth stations," IEEE 10th International Conference on Advanced Communication Technology, vol. 3, pp. 1559-1563, .

[12] P. Sengupta, S. R. Vadali, K. T. Alfriend, "Satellite orbit design and maintenance for terrestrial coverage," Journal of Spacecraft and Rockets, vol. 47 no. 1, pp. 177-187, DOI: 10.2514/1.44120, 2010.

[13] J. J. Davis, M. E. Avendaño, D. Mortari, "The 3-D lattice theory of flower constellations," Celestial Mechanics and Dynamical Astronomy, vol. 116 no. 4, pp. 339-356, DOI: 10.1007/s10569-013-9494-7, 2013.

[14] S. Lee, D. Mortari, "2-D lattice flower constellations for radio occultation missions," Frontiers in Aerospace Engineering, vol. 2 no. 2, pp. 79-90, 2013.

[15] Z. Song, X. Hu, M. Wang, G. Dai, "Judgement theorems and an approach for solving the constellation-to-ground coverage problem," Mathematical Problems in Engineering, vol. 2018 no. 2, 2018.

[16] H. Wang, C. Han, S. Liu, Y. Sun, H. Zhang, "Adaptive algorithm to determine the coverage belt for agile satellite with attitude maneuvers," IEEE 10th International Conference on Mechanical and Aerospace Engineering (ICMAE), .

[17] M. Zuo, G. Dai, L. Peng, M. Wang, "An envelope curve-based theory for the satellite coverage problems," Aerospace Science and Technology, vol. 100 no. 2020, article 105750, 2020.

[18] J. J. Morrison, A System of Sixteen Synchronous Satellites for World Wise Navigation and Surveillance, 1973.

[19] H. Li, D. Li, Y. Li, "A multi-index assessment method for evaluating coverage effectiveness of remote sensing satellite," Chinese Journal of Aeronautics, vol. 31 no. 10, pp. 2023-2033, DOI: 10.1016/j.cja.2018.05.015, 2018.

[20] Z. Song, G. Dai, M. Wang, X. Chen, "A novel grid point approach for efficiently solving the constellation-to-ground regional coverage problem," IEEE Access, vol. 6, pp. 44445-44458, DOI: 10.1109/ACCESS.2018.2864181, 2018.

[21] L. Rider, "Optimized polar orbit constellations for redundant earth coverage," Journal of the Astronautical Sciences, vol. 33, pp. 147-161, 1985.

[22] W. S. Adams, L. Rider, "Circular polar constellations providing continuous single or multiple coverage above a specified latitude," Journal of the Astronautical Sciences, vol. 35 no. 2, pp. 155-192, 1987.

[23] M. Yang, X. Dong, M. Hu, "Design and simulation for hybrid LEO communication and navigation constellation," IEEE Chinese Guidance, Navigation and Control Conference (CGNCC), pp. 1665-1669, .

[24] A. H. Ballard, "Rosette constellations of earth satellites," IEEE Transactions on Aerospace and Electronic Systems, vol. 16 no. 5, pp. 656-673, 1980.

[25] J. G. Walker, Continuous Whole-Earth Coverage by Circular-Orbit Satellite Patterns, 1977.

[26] G. Dai, X. Chen, M. Wang, E. Fernández, T. N. Nguyen, G. Reinelt, "Analysis of satellite constellations for the continuous coverage of ground regions," Journal of Spacecraft and Rockets, vol. 54 no. 6, pp. 1294-1303, DOI: 10.2514/1.A33826, 2017.

[27] M. Wang, X. Luo, G. Dai, X. Chen, "Application of latitude stripe division in satellite constellation coverage to ground," International Journal of Aerospace Engineering, vol. 2016 no. 10, 2016.

[28] R. Lüders, L. Ginsberg, "Continuous zonal coverage-a generalized analysis," Mechanics and Control of Flight Conference, Mechanics and Control of Flight Conference, .

[29] Q. Chen, Y. Bai, L. Chen, Z. Pang, "Design of LEO constellations providing internet services based on soc method," MATEC Web of Conferences, vol. 114, .

[30] X. Chen, G. Dai, G. Reinelt, M. Wang, "A semi-analytical method for periodic earth coverage satellites optimization," IEEE Communications Letters, vol. 22 no. 3, pp. 534-537, DOI: 10.1109/LCOMM.2017.2780107, 2018.