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

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 e ﬃ ciency signi ﬁ cantly 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 di ﬀ erent types of coverage problems are conducted, and the results reveal the correctness and high e ﬃ ciency of the proposed analytical method.


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 constellationto-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 thirdorder 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][19][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][22][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.

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 ∈ 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 Λ C . For each point P ∈ D R , determine whether it is in Λ C . Count the number of points in D R that is also in Λ C , denoted by N. This step is shown in the fourth subfigure in Figure 1. The green points stand for the points in Λ C . The coverage rate of the constellation to the region is calculated, which is denoted by η; then where card ðD R Þ represents the number of elements in D R . 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.

Cell Area Analytical Method
In this section, we propose a new method for efficiently and exactly calculating the coverage area for the constellationto-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 Θ and denote the spherical circle by ∂Θ.
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.
For the spherical circle ∂Θ, 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 m½U represents the area of a spherical region U. The radius of the Earth is unitized. The area of the spherical disk is as follows: Denote the longitude and latitude of a point P on the Earth's surface as ðθ P , ϕ P Þ and Q as ðθ Q , ϕ Q Þ. The spherical distance between the two points, denoted by d, can be expressed by the following formula: The relationship between d and its corresponding spherical central angle α can be expressed as follows: The spherical sector with central angle α in Figure 2 is written as ∇ OPQ . The area of ∇ OPQ is as follows: The spherical triangle with vertices O, P, Q is denoted by The spherical arch region in Figure 2 is written as A PQ ðΘÞ. The area of this spherical arch region is Therefore, we can obtain the area of the spherical arch only if the radius r and the two points P and Q are known.

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 3 International Journal of Aerospace Engineering relation and is written as ΠðG k , UÞ, can be divided into three categories as follows: where U C represents the complementary set of U with respect to the Earth's surface where the symbol "| =" represents a set that is partly contained by another set Suppose the longitude range and latitude range of G k are ½θ L k , θ U k and ½ϕ L k , ϕ U k , respectively. Normalize the radius of the Earth as a unit length; then, the area of the cell G k is expressed as follows: 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 Θ j . The intersection area between a cell G k and a spherical disk Θ is denoted by m½G k ∩ Θ. The calculation of m½G k ∩ Θ is a core algorithm in this paper. According to the cell-region relation, it can be classified into four categories as follows: However, in general, the size of a cell is much smaller than that of a spherical disk. Therefore, this category need not be considered the computation process of m½G k ∩ Θ is very complex, which will be discussed in detail in this section For G k | = Θ, except for very extreme cases, G k and Θ have two intersection points in general. We analyze a typical case for calculating m½G k ∩ Θ when G k | = Θ.
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].
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 Δθ, then 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 π/2 − ϕ U k , and the center angle is θ N − θ P . Thus, the area of S 4 can be calculated.
Then, m½G k ∩ Θ can be expressed as follows: 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.

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 ∂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 Additionally, by the conclusion in the former section, there are three types of relationships between G k and U: G k ⊂ U, G k ⊂ U C , and G k | = U. If G k | = 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 ∂U crosses the cell (ii) Type-II: more than one circular arc in ∂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 Θ j . For U, if only one circular arc in ∂U crosses cell G k , the intersection area of U and G k is equivalent to that of Θ j and G k .
We can use Equation (11) and its variant formulas to compute m½U ∩ G k .
However, if more than one circular arc in ∂U crosses cell G k , then the intersection area of U and G k cannot be simply computed. To compute the precision result of m½U ∩ G k , cell G k should be partitioned into a group of subcells.
If G k is divided into N subcells, the jth subcell is written as G k,j , that is, Then, 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.

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 Λ S , and the instantaneous coverage of constellation C is denoted by Λ C . Then, For a Λ S at a certain time, it corresponds to a ground region. By Subsection 3.4, we can compute the result of m½G k ∩ Λ S . If ∂Λ S consists of more than one spherical arc, G k may belong to type-II with respect to Λ S . We can then compute the result of m½G k ∩ Λ S .
For the computation of m½G k ∩ Λ 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 ∩ Λ C = 0 (ii) For red cells, m½G k ∩ Λ C = m½G k (iii) For blue cells, G k is only partly covered by one satellite S; then, m½G k ∩ Λ C = m½G k ∩ Λ S (iv) For green cells, G k belongs to type-II with respect to Λ C . We need to divide G k into subcells to obtain a precision result Figure 4: Two typical cases of a spherical disk intersecting a cell. 5 International Journal of Aerospace Engineering 3.6. Computation of Constellation to Ground-Region Coverage Problems. In constellation to ground-region coverage problems, a cell G k has two important properties, that is, the relationship between the cell and the ground region and the relationship between the cell and the constellation coverage regions, which are written as ΠðG k , Λ C Þ and ΠðG k , RÞ, respectively. What we need to compute is m½G k ∩ Λ C ∩ R.
In some combinations of ΠðG k , Λ C Þ and ΠðG k , RÞ, the value of m½G k ∩ Λ C ∩ R can be computed by the conclusion given in former sections. However, in other combinations of ΠðG k , Λ C Þ and ΠðG k , RÞ, m½G k ∩ Λ C ∩ R cannot be computed directly, and G k should be divided into a group of subcells. The result of m½G k ∩ Λ C ∩ R under the combinations of ΠðG k , Λ C Þ and ΠðG k , RÞ is given in Table 1.
In Table 1, the symbol × indicates that G k needs to be further subdivided in the case of this combination. For other cases, m½G k ∩ Λ C ∩ R can be directly computed. When m½G k ∩ Λ C ∩ R can be directly computed, in some cases, m½G k ∩ Λ C ∩ R = 0, while in other cases, m½G k ∩ Λ C ∩ R can be calculated by the equation given in former sections.
Then, for the set of fG k g, according to Table 1, it can be classified into two subsets and denoted by V 1 and V 2 , respectively, where 3.7. Basic Steps of the Cell Area Analytical Method. Computing the coverage rate is a common indicator in coverage problems. The coverage rate is defined as follows: International Journal of Aerospace Engineering Generally, m½R is a known quantity; therefore, the core problem for the coverage rate is computing m½R ∩ Λ C .
According to the preview discussion, we can propose the cell area analytical method. The basic steps of this method are as follows: Step 1. Obtain the minimum latitude and longitude box of the region R, which is denoted by U R .
Step 2. Select a division precision ε and divide U R into a group of cell sets fG k g.
Step 3. For every cell G k , compute the results of ΠðG k , RÞ and ΠðG k , Λ C Þ. Additionally, we judge whether G k belongs to type-I or type-II with respect to R and Λ C , respectively.
Step 4. According to Equation (14) and results in Step 3, classify fG k g into V 1 and V 2 . Then, compute the value of M 1 and M 2 : Step 5. If the calculation accuracy does not satisfy the requirement, dividing every cell G k in V 2 into four subcells. For every subcell, use the same process as Step 3 and Step 4 to determine whether the subcell belongs to V 1 or V 2 . Update the value of M 1 and M 2 . Repeat the process until the calculation accuracy satisfies the requirement.
Step 6. Obtain the result of coverage rate. Because for any precision, there are some cells that m½G k ∩ Λ C ∩ R cannot directly compute, so the exact coverage rate cannot be obtained. However, we can then compute the upper and lower results of the coverage rate. The upper and lower results of the coverage rate are written as η sup and η inf , respectively. Then,

Experimental Results
In this section, we analyze a Walker-Delta constellation with configuration parameters ðT/P/FÞ = ð36/4/1Þ in a geopotential model of Earth. It includes the point mass effect as well as the dominant effect of the asymmetry in the gravitational field, which is particularly useful for modeling "ideal" accurately maintained orbit (note that the proposed algorithm can be applied to any orbit propagator model) [30]. All the satellites in the constellation are in a circular orbit and have the same altitude of 1300 km and the same inclination of 45°. The right ascension of the ascending node (RAAN) and the true anomaly of the first satellite in the constellation are both 0°at 2020-01-01 00:00:00. Every satellite carries a simple conic sensor with a beam angle of 45°. Two different ground targets are considered as follows: Target-1: a latitude and longitude box with a longitude interval of ½0°, 80° and a latitude interval of ½10°, 50°. It is a regular spherical region.
Target-2: the contiguous United States is a concave polygon of irregular shape, and its boundary consists of 3369 latitude and longitude points. It is an arbitrarily shaped spherical polygon.
For the cell area analytical method, the partition accuracy is the most important factor affecting the result precision. We define the cell precision as follows: Definition 2 Cell precision. Cell precision is defined as the square root of the cell number in the 1 km × 1 km region.
The computational results for the cell analytical method and the traditional net point method are depicted in Figures 6-9. The x axis indicates the partition precision, given as the number of cells or net points in 1 km on Earth's surface. The y axis shows the instantaneous coverage rate and gaps of the coverage boundary for targets, separately in Figures 6 and 8 and in Figures 7 and 9.
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 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 highprecision 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 ∈ R and G k ∈ Λ 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 Coverage rate difference Figure 9: Comparison of the cell counts from two algorithms with different precisions for target-2.    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.

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.

Data Availability
The data used to support the findings of this study are available from the first author upon request.