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

Authors: Zhiming Song; Haidong Liu; Guangming Dai; Maocai Wang; Xiaoyu Chen
Journal: International Journal of Aerospace Engineering
 (2021)

---

## 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.


---

## Body

## 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 Figure1. The details of this process are as follows.Step 1.
For a certain target regionR, 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 DR. 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 pointP∈DR, determine whether it is in R. If a point P is not in R, it needs to be removed from DR. 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∈DR, determine whether it is in ΛC. Count the number of points in DR 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
(1)η=NcardDR,where cardDR represents the number of elements in DR.Figure 1
Basic steps of the net point method.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 increasesN-fold, the computing time will increase N2-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 Figure2, 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 2
Schematic diagram 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 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:
(2)mΘ=2π⋅1−cosr.Denote the longitude and latitude of a pointP 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:
(3)d=cos−1cosϕPcosϕQcosθP−θQ+sinϕPsinϕQ.The relationship betweend and its corresponding spherical central angle α can be expressed as follows:
(4)α=2⋅sin−1sind/2sinr.The spherical sector with central angleα in Figure 2 is written as ∇OPQ. The area of ∇OPQ is as follows:
(5)m∇OPQ=1−cosr⋅α.The spherical triangle with verticesO,P,Q is denoted by ΔOPQ. The area of ΔOPQ is
(6)mΔOPQ=2tan−11cosr⋅tanα/2+α−π.The spherical arch region in Figure2 is written as APQΘ. The area of this spherical arch region is
(7)mAPQΘ=m∇OPQ−mΔOPQ.Therefore, we can obtain the area of the spherical arch only if the radiusr 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 byGk. The relationship between a cell Gk and a ground region U, which is called the cell-region relation and is written as ΠGk,U, can be divided into three categories as follows:
(i)
Gk is inside of U, which is written as Gk⊂U(ii)
Gk is outside of U, which is written as Gk⊂UC, where UC represents the complementary set of U with respect to the Earth’s surface(iii)
Gk is partly inside of U, which is written as Gk∣=U, where the symbol “∣=” represents a set that is partly contained by another setSuppose the longitude range and latitude range ofGk are θkL,θkU and ϕkL,ϕkU, respectively. Normalize the radius of the Earth as a unit length; then, the area of the cell Gk is expressed as follows:
(8)mGk=sinϕkU−sinϕkLθkU−θkL.


### 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 indexj is written as Θj.The intersection area between a cellGk and a spherical disk Θ is denoted by mGk∩Θ. The calculation of mGk∩Θ is a core algorithm in this paper. According to the cell-region relation, it can be classified into four categories as follows:
(i)
IfΘ⊂Gk, then mGk∩Θ=mΘ. 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)
IfGk⊂Θ, then mGk∩Θ=mGk, which can be computed by Equation (8)(iii)
IfGk⊂ΘC, then mGk∩Θ=0(iv)
IfGk∣=Θ, the computation process of mGk∩Θ is very complex, which will be discussed in detail in this sectionForGk∣=Θ, except for very extreme cases, Gk and Θ have two intersection points in general. We analyze a typical case for calculating mGk∩Θ when Gk∣=Θ.In Figure3, 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 3
A typical case of a spherical disk intersecting a cell.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 Figure3, the intersection area is divided into three regions, namely, the blue region S1, the yellow region S2, and the green region S3.The blue regionS1 is a spherical longitude and latitude box. Supposing the value of the longitude interval of Q and M is Δθ, then
(9)mS1=sinϕkU−sinϕkL⋅Δθ.The yellow regionS2 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,S3 is not a spherical triangle. In obtaining the area of S3, the spherical triangle with the same vertexes of S3 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 S4. The spherical triangle with the same vertexes of S3 corresponds to the combination of S3 and S4. 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 S4 is π/2−ϕkU, and the center angle is θN−θP. Thus, the area of S4 can be calculated.Then,mGk∩Θ can be expressed as follows:
(10)mGk∩Θ=mS1+mS2+mΔPQN−mS4.Obviously, there are many other cases in which a spherical disk intersects a cell. There exists a total of 16 cases by analysis. Figure4 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.Figure 4
Two typical cases of a spherical disk intersecting a cell.
(a)(b)(c)


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

LetU 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 called spherical circular arcs. A spherical circular arc with index j is written as Aj; then, Aj∈∂U.Additionally, by the conclusion in the former section, there are three types of relationships betweenGk and U: Gk⊂U, Gk⊂UC, and Gk∣=U. If Gk∣=U, then Gk 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 cellAccording to the spherical geometry, a random spherical circular arcAj specifying the inner side corresponds to a spherical disk, denoted by Θj. For U, if only one circular arc in ∂U crosses cell Gk, the intersection area of U and Gk is equivalent to that of Θj and Gk.
(11)mU∩Gk=mΘj∩Gk.We can use Equation (11) and its variant formulas to compute mU∩Gk.However, if more than one circular arc in∂U crosses cell Gk, then the intersection area of U and Gk cannot be simply computed. To compute the precision result of mU∩Gk, cell Gk should be partitioned into a group of subcells.IfGk is divided into N subcells, the jth subcell is written as Gk,j, that is,
(12)Gk=⋃j=1NGk,j.Then,(13)mU∩Gk=∑j=1NmU∩Gk,j.Gk,j is smaller than Gk; in general, part of the subcells in Gk,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 byC. 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,(14)ΛC=⋃S∈CΛS.For aΛS at a certain time, it corresponds to a ground region. By Subsection 3.4, we can compute the result of mGk∩ΛS. If ∂ΛS consists of more than one spherical arc, Gk may belong to type-II with respect to ΛS. We can then compute the result of mGk∩ΛS.For the computation ofmGk∩ΛC, we also need to consider the case where Gk is covered by more than one satellite simultaneously.Consider the situation shown in Figure5. 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,mGk∩ΛC=0(ii)
For red cells,mGk∩ΛC=mGk(iii)
For blue cells,Gk is only partly covered by one satellite S; then, mGk∩ΛC=mGk∩ΛS(iv)
For green cells,Gk belongs to type-II with respect to ΛC. We need to divide Gk into subcells to obtain a precision resultFigure 5
The relationship between the cells and the constellation coverage region.


### 3.6. Computation of Constellation to Ground-Region Coverage Problems

In constellation to ground-region coverage problems, a cellGk 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 ΠGk,ΛC and ΠGk,R, respectively. What we need to compute is mGk∩ΛC∩R.In some combinations ofΠGk,ΛC and ΠGk,R, the value of mGk∩ΛC∩R can be computed by the conclusion given in former sections. However, in other combinations of ΠGk,ΛC and ΠGk,R, mGk∩ΛC∩R cannot be computed directly, and Gk should be divided into a group of subcells. The result of mGk∩ΛC∩R under the combinations of ΠGk,ΛC and ΠGk,R is given in Table 1.Table 1
The result ofmGk∩ΛC∩R under the combinations of ΠGk,ΛC and ΠGk,R.
Gk⊂ΛCCGk∣=ΛC type-IGk∣=ΛC type-IIGk⊂ΛCGk⊂RC0000Gk∣=R type-I0××mGk∩RGk∣=R type-II0×××Gk⊂R0mGk∩ΛC×mGkIn Table1, the symbol × indicates that Gk needs to be further subdivided in the case of this combination. For other cases, mGk∩ΛC∩R can be directly computed. When mGk∩ΛC∩R can be directly computed, in some cases, mGk∩ΛC∩R=0, while in other cases, mGk∩ΛC∩R can be calculated by the equation given in former sections.Then, for the set ofGk, according to Table 1, it can be classified into two subsets and denoted by V1 and V2, respectively, where
(15)V1=GkmGk∩ΛC∩Rcanbedirectlycomputed,V2=GkV1.


### 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:(16)η=mR∩ΛCmR.Generally,mR is a known quantity; therefore, the core problem for the coverage rate is computing mR∩Λ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 regionR, which is denoted by UR.Step 2.
Select a division precisionε and divide UR into a group of cell sets Gk.Step 3.
For every cellGk, compute the results of ΠGk,R and ΠGk,ΛC. Additionally, we judge whether Gk 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, classifyGk into V1 and V2. Then, compute the value of M1 and M2:
(17)M1=∑Gk⊂V1mGk∩ΛC∩R,M2=∑Gk⊂V2mGk.Step 5.
If the calculation accuracy does not satisfy the requirement, dividing every cellGk in V2 into four subcells. For every subcell, use the same process as Step 3 and Step 4 to determine whether the subcell belongs to V1 or V2. Update the value of M1 and M2. 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 thatmGk∩Λ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,
(18)ηinf=M1mR,ηsup=M1+M2mR.


## 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 Figure2, 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 2
Schematic diagram 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 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:
(2)mΘ=2π⋅1−cosr.Denote the longitude and latitude of a pointP 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:
(3)d=cos−1cosϕPcosϕQcosθP−θQ+sinϕPsinϕQ.The relationship betweend and its corresponding spherical central angle α can be expressed as follows:
(4)α=2⋅sin−1sind/2sinr.The spherical sector with central angleα in Figure 2 is written as ∇OPQ. The area of ∇OPQ is as follows:
(5)m∇OPQ=1−cosr⋅α.The spherical triangle with verticesO,P,Q is denoted by ΔOPQ. The area of ΔOPQ is
(6)mΔOPQ=2tan−11cosr⋅tanα/2+α−π.The spherical arch region in Figure2 is written as APQΘ. The area of this spherical arch region is
(7)mAPQΘ=m∇OPQ−mΔOPQ.Therefore, we can obtain the area of the spherical arch only if the radiusr 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 byGk. The relationship between a cell Gk and a ground region U, which is called the cell-region relation and is written as ΠGk,U, can be divided into three categories as follows:
(i)
Gk is inside of U, which is written as Gk⊂U(ii)
Gk is outside of U, which is written as Gk⊂UC, where UC represents the complementary set of U with respect to the Earth’s surface(iii)
Gk is partly inside of U, which is written as Gk∣=U, where the symbol “∣=” represents a set that is partly contained by another setSuppose the longitude range and latitude range ofGk are θkL,θkU and ϕkL,ϕkU, respectively. Normalize the radius of the Earth as a unit length; then, the area of the cell Gk is expressed as follows:
(8)mGk=sinϕkU−sinϕkLθkU−θkL.


## 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 indexj is written as Θj.The intersection area between a cellGk and a spherical disk Θ is denoted by mGk∩Θ. The calculation of mGk∩Θ is a core algorithm in this paper. According to the cell-region relation, it can be classified into four categories as follows:
(i)
IfΘ⊂Gk, then mGk∩Θ=mΘ. 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)
IfGk⊂Θ, then mGk∩Θ=mGk, which can be computed by Equation (8)(iii)
IfGk⊂ΘC, then mGk∩Θ=0(iv)
IfGk∣=Θ, the computation process of mGk∩Θ is very complex, which will be discussed in detail in this sectionForGk∣=Θ, except for very extreme cases, Gk and Θ have two intersection points in general. We analyze a typical case for calculating mGk∩Θ when Gk∣=Θ.In Figure3, 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 3
A typical case of a spherical disk intersecting a cell.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 Figure3, the intersection area is divided into three regions, namely, the blue region S1, the yellow region S2, and the green region S3.The blue regionS1 is a spherical longitude and latitude box. Supposing the value of the longitude interval of Q and M is Δθ, then
(9)mS1=sinϕkU−sinϕkL⋅Δθ.The yellow regionS2 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,S3 is not a spherical triangle. In obtaining the area of S3, the spherical triangle with the same vertexes of S3 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 S4. The spherical triangle with the same vertexes of S3 corresponds to the combination of S3 and S4. 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 S4 is π/2−ϕkU, and the center angle is θN−θP. Thus, the area of S4 can be calculated.Then,mGk∩Θ can be expressed as follows:
(10)mGk∩Θ=mS1+mS2+mΔPQN−mS4.Obviously, there are many other cases in which a spherical disk intersects a cell. There exists a total of 16 cases by analysis. Figure4 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.Figure 4
Two typical cases of a spherical disk intersecting a cell.
(a)(b)(c)


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

LetU 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 called spherical circular arcs. A spherical circular arc with index j is written as Aj; then, Aj∈∂U.Additionally, by the conclusion in the former section, there are three types of relationships betweenGk and U: Gk⊂U, Gk⊂UC, and Gk∣=U. If Gk∣=U, then Gk 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 cellAccording to the spherical geometry, a random spherical circular arcAj specifying the inner side corresponds to a spherical disk, denoted by Θj. For U, if only one circular arc in ∂U crosses cell Gk, the intersection area of U and Gk is equivalent to that of Θj and Gk.
(11)mU∩Gk=mΘj∩Gk.We can use Equation (11) and its variant formulas to compute mU∩Gk.However, if more than one circular arc in∂U crosses cell Gk, then the intersection area of U and Gk cannot be simply computed. To compute the precision result of mU∩Gk, cell Gk should be partitioned into a group of subcells.IfGk is divided into N subcells, the jth subcell is written as Gk,j, that is,
(12)Gk=⋃j=1NGk,j.Then,(13)mU∩Gk=∑j=1NmU∩Gk,j.Gk,j is smaller than Gk; in general, part of the subcells in Gk,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 byC. 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,(14)ΛC=⋃S∈CΛS.For aΛS at a certain time, it corresponds to a ground region. By Subsection 3.4, we can compute the result of mGk∩ΛS. If ∂ΛS consists of more than one spherical arc, Gk may belong to type-II with respect to ΛS. We can then compute the result of mGk∩ΛS.For the computation ofmGk∩ΛC, we also need to consider the case where Gk is covered by more than one satellite simultaneously.Consider the situation shown in Figure5. 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,mGk∩ΛC=0(ii)
For red cells,mGk∩ΛC=mGk(iii)
For blue cells,Gk is only partly covered by one satellite S; then, mGk∩ΛC=mGk∩ΛS(iv)
For green cells,Gk belongs to type-II with respect to ΛC. We need to divide Gk into subcells to obtain a precision resultFigure 5
The relationship between the cells and the constellation coverage region.


## 3.6. Computation of Constellation to Ground-Region Coverage Problems

In constellation to ground-region coverage problems, a cellGk 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 ΠGk,ΛC and ΠGk,R, respectively. What we need to compute is mGk∩ΛC∩R.In some combinations ofΠGk,ΛC and ΠGk,R, the value of mGk∩ΛC∩R can be computed by the conclusion given in former sections. However, in other combinations of ΠGk,ΛC and ΠGk,R, mGk∩ΛC∩R cannot be computed directly, and Gk should be divided into a group of subcells. The result of mGk∩ΛC∩R under the combinations of ΠGk,ΛC and ΠGk,R is given in Table 1.Table 1
The result ofmGk∩ΛC∩R under the combinations of ΠGk,ΛC and ΠGk,R.
Gk⊂ΛCCGk∣=ΛC type-IGk∣=ΛC type-IIGk⊂ΛCGk⊂RC0000Gk∣=R type-I0××mGk∩RGk∣=R type-II0×××Gk⊂R0mGk∩ΛC×mGkIn Table1, the symbol × indicates that Gk needs to be further subdivided in the case of this combination. For other cases, mGk∩ΛC∩R can be directly computed. When mGk∩ΛC∩R can be directly computed, in some cases, mGk∩ΛC∩R=0, while in other cases, mGk∩ΛC∩R can be calculated by the equation given in former sections.Then, for the set ofGk, according to Table 1, it can be classified into two subsets and denoted by V1 and V2, respectively, where
(15)V1=GkmGk∩ΛC∩Rcanbedirectlycomputed,V2=GkV1.


## 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:(16)η=mR∩ΛCmR.Generally,mR is a known quantity; therefore, the core problem for the coverage rate is computing mR∩Λ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 regionR, which is denoted by UR.Step 2.
Select a division precisionε and divide UR into a group of cell sets Gk.Step 3.
For every cellGk, compute the results of ΠGk,R and ΠGk,ΛC. Additionally, we judge whether Gk 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, classifyGk into V1 and V2. Then, compute the value of M1 and M2:
(17)M1=∑Gk⊂V1mGk∩ΛC∩R,M2=∑Gk⊂V2mGk.Step 5.
If the calculation accuracy does not satisfy the requirement, dividing every cellGk in V2 into four subcells. For every subcell, use the same process as Step 3 and Step 4 to determine whether the subcell belongs to V1 or V2. Update the value of M1 and M2. 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 thatmGk∩Λ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,
(18)ηinf=M1mR,ηsup=M1+M2mR.


## 4. Experimental Results

In this section, we analyze a Walker-Delta constellation with configuration parametersT/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 the1km×1km region.The computational results for the cell analytical method and the traditional net point method are depicted in Figures6–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.Figure 6
Results of coverage rate with an increase of the cell precision for target-1.Figure 7
The differences between the upper and lower bound of coverage rate for target-1.Figure 8
Results of coverage rate with an increase of the cell precision for target-2.Figure 9
Comparison of the cell counts from two algorithms with different precisions for target-2.The computation results indicate that the boundary of coverage rate converges to the same value. Figures7 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.Figures10 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 10
Comparison of the computing time from two algorithms with different precision for target-1.Figure 11
Comparison of the computing time from two algorithms with different precision for target-2.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 Figure10, 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.Figures12 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 Gk∈R and Gk∈Λ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 12
The coverage region of simulation for target-1.Figure 13
The coverage region of simulation 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 Figures12 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.



---

*Source: 6679107-2021-06-21.xml*