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.
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 , obtain its minimum latitude and longitude encirclement. Then, a group of points in the minimum latitude and longitude encirclement of 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 . 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 , determine whether it is in . If a point is not in , it needs to be removed from . Only the points inside region 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 . For each point , determine whether it is in . Count the number of points in that is also in , denoted by . This step is shown in the fourth subfigure in Figure 1. The green points stand for the points in . The coverage rate of the constellation to the region is calculated, which is denoted by ; then where represents the number of elements in .
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 -fold, the computing time will increase -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 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 . The symbol represents the area of a spherical region . The radius of the Earth is unitized. The area of the spherical disk is as follows:
Denote the longitude and latitude of a point on the Earth’s surface as and as . The spherical distance between the two points, denoted by , can be expressed by the following formula:
The relationship between and its corresponding spherical central angle can be expressed as follows:
The spherical sector with central angle in Figure 2 is written as . The area of is as follows:
The spherical triangle with vertices is denoted by . The area of is
The spherical arch region in Figure 2 is written as . The area of this spherical arch region is
Therefore, we can obtain the area of the spherical arch only if the radius and the two points and 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 . The relationship between a cell and a ground region , which is called the cell-region relation and is written as , can be divided into three categories as follows: (i) is inside of , which is written as (ii) is outside of , which is written as , where represents the complementary set of with respect to the Earth’s surface(iii) is partly inside of , which is written as , where the symbol “” represents a set that is partly contained by another set
Suppose the longitude range and latitude range of are and , respectively. Normalize the radius of the Earth as a unit length; then, the area of the cell 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 is written as .
The intersection area between a cell and a spherical disk is denoted by . The calculation of is a core algorithm in this paper. According to the cell-region relation, it can be classified into four categories as follows: (i)If , then . 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 , then , which can be computed by Equation (8)(iii)If , then (iv)If , the computation process of is very complex, which will be discussed in detail in this section
For , except for very extreme cases, and have two intersection points in general. We analyze a typical case for calculating when .
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. and are the two intersection points. The longitude and latitude coordinates of and 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 , the yellow region , and the green region .
The blue region is a spherical longitude and latitude box. Supposing the value of the longitude interval of and is , then
The yellow region is a typical spherical arch region. If the given coordinates of and , 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, is not a spherical triangle. In obtaining the area of , the spherical triangle with the same vertexes of should be firstly constructed, and then, the extra area is subtracted. We make a great circle by the two points and . The region surrounded by the latitude line and the great circle by the two points and is denoted by . The spherical triangle with the same vertexes of corresponds to the combination of and . The area of the generated spherical triangle can be obtained by the given coordinates , , . The radius of this generated spherical arch region is , and the center angle is . Thus, the area of can be calculated.
Then, 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.
(a)
(b)
(c)
3.4. The Intersection Area of a Ground Region and a Cell
Let be a ground region with an arbitrary shape. The edge of , which can be written as , 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 is written as ; then, .
Additionally, by the conclusion in the former section, there are three types of relationships between and : , , and . If , then 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 crosses the cell(ii)Type-II: more than one circular arc in crosses the cell
According to the spherical geometry, a random spherical circular arc specifying the inner side corresponds to a spherical disk, denoted by . For , if only one circular arc in crosses cell , the intersection area of and is equivalent to that of and .
We can use Equation (11) and its variant formulas to compute .
However, if more than one circular arc in crosses cell , then the intersection area of and cannot be simply computed. To compute the precision result of , cell should be partitioned into a group of subcells.
If is divided into subcells, the th subcell is written as , that is,
Then,
is smaller than ; in general, part of the subcells in belong to type-I, and the others belong to type-II. The intersection area between 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 . The instantaneous coverage of satellite to the ground is denoted by , and the instantaneous coverage of constellation is denoted by .
Then,
For a at a certain time, it corresponds to a ground region. By Subsection 3.4, we can compute the result of . If consists of more than one spherical arc, may belong to type-II with respect to . We can then compute the result of .
For the computation of , we also need to consider the case where 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, (ii)For red cells, (iii)For blue cells, is only partly covered by one satellite ; then, (iv)For green cells, belongs to type-II with respect to . We need to divide into subcells to obtain a precision result
3.6. Computation of Constellation to Ground-Region Coverage Problems
In constellation to ground-region coverage problems, a cell 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 and , respectively. What we need to compute is .
In some combinations of and , the value of can be computed by the conclusion given in former sections. However, in other combinations of and , cannot be computed directly, and should be divided into a group of subcells. The result of under the combinations of and is given in Table 1.
In Table 1, the symbol indicates that needs to be further subdivided in the case of this combination. For other cases, can be directly computed. When can be directly computed, in some cases, , while in other cases, can be calculated by the equation given in former sections.
Then, for the set of , according to Table 1, it can be classified into two subsets and denoted by and , 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:
Generally, is a known quantity; therefore, the core problem for the coverage rate is computing .
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 , which is denoted by .
Step 2. Select a division precision and divide into a group of cell sets .
Step 3. For every cell , compute the results of and . Additionally, we judge whether belongs to type-I or type-II with respect to and , respectively.
Step 4. According to Equation (14) and results in Step 3, classify into and . Then, compute the value of and :
Step 5. If the calculation accuracy does not satisfy the requirement, dividing every cell in into four subcells. For every subcell, use the same process as Step 3 and Step 4 to determine whether the subcell belongs to or . Update the value of and . 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 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 and , respectively. Then,
4. Experimental Results
In this section, we analyze a Walker-Delta constellation with configuration parameters 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 and a latitude interval of . 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 region.
The computational results for the cell analytical method and the traditional net point method are depicted in Figures 6–9. The axis indicates the partition precision, given as the number of cells or net points in 1 km on Earth’s surface. The 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 point division method with an increase of the partition precision.
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 and . 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.
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.
Data Availability
The data used to support the findings of this study are available from the first author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
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).