A Simple Three-Dimensional Matrix Method for Global Constellation Intrasatellite Link Topological Design

Abstract

A three-dimensional matrix method is proposed in this paper for Global Navigation Satellite System (GNSS) constellation Intrasatellite Link (ISL) topological design. The rows and columns of proposed matrix contain the information of both constellation orbit planes and satellites in each orbit plane. The third dimension of proposed matrix represents the time sequences during constellation movement. The proposed method has virtues of better advantage of conceptual clarity and computational efficiency, meanwhile, some properties of ISLs in the constellation can be proved easily. At the second part of this paper, a link assignment and optimal routing problem is proposed using three-dimensional topology matrixes, which aimed to minimize the relay hops during the process of data uploading to whole constellation network. Moreover, some practical constrains as antenna beam coverage and relative velocity are considered and analyzed in detail. Finally, some numerical simulations are provided, and the results demonstrated the promising performance of proposed topological method in reduction of computation burden, clear ISL conception, and so forth; the efficiency of provided optimal ISL routing problem is also proved.

1. Introduction

With the rapid progress in GNSS technology, it is now feasible and necessary to build ISL network to increase system functions as communications and range measurements. One of the most important aspects in satellite constellation ISL analysis is the topological design.

In previous researches, Werner [1–3] provided geometry of Virtual Path Connection (VPC) topological representation method for the analysis of the ISL between the satellites. In a series of literatures of Werner, “LEONET,” a Medium Earth Orbit (MEO) satellite constellation has been taken as an example, which is on the basis of proposed topology method. Moreover, the actual needs for global data relay and transmission are considered. By using the idea of Asynchronous Transfer Mode (ATM), Werner also provided a Modified Dijkstra Shortest Path Algorithm (M-DSPA) to achieve the global transmission of data with end-to-end performance. Donner et al. [4] carried out in-depth study of routing algorithm according to the Werner’s topology method, and Lee and Kang [5] also provided an application case of global Search and Rescue (SAR) using Werner’s results.

Chang et al. [6] proposed a graphical method for ISL topology design for a Low Earth Orbit (LEO) communication satellite constellation, and a two-step optimized search routing algorithm was given for global data relay and transmission.

Suzuki et al. [7] provided a topological method that has been successfully applied to the Next-Generation LEO System (NeLS) project which will consist of 120 satellites in constellation with inter satellite link connections. The NeLS can provide global multimedia mobile satellite communication service for handy terminals, which is conducted by the authorities of Japan. In the literature [7], Suzuki et al. considered the use of two-dimensional matrix for ISL topological, while the matrix rows are used for different orbit planes and the matrix columns are used to represent each satellite in corresponding orbit plane. This method has significantly simplified the ISL analyzed and was adopted by some other researchers [8].

As we can see from the analysis above, previous constellation ISL designs are basically based on intuitive graphical analysis for a task-specific mission. Some of the ISL topology analyses are using matrix form, however, those of which are mostly graphical based and without any matrix associated properties analysis. In this paper, the authors proposed a three-dimensional matrix topology analysis method, wherein the rows and columns of the matrix represent the orbital plane and satellite information, and the third dimension represents the time information. By theoretical analysis and numerical simulations, the results show that the proposed method has advantages of computational efficiency and clarity of ISL conceptual.

At the second stage of this paper, an optimal link assignment and routing problem is provided using proposed matrix topology method. The readers will find more details in Section 4.

2. Matrix Model

2.1. The Elements of ISL Topology Matrix

Assuming that a typical global navigation constellation consists of 𝑁 orbital planes, each orbital plane have 𝑀 satellites (to simplify the notation, this section consider a constellation of three orbital planes, with four satellites at each plane; that is 𝑁=3,??𝑀=4). Then the ISL topology matrix model is shown as in Figure 1.

Figure 1 

Demonstration of 3-orbit plane ISL topological.

The rows of the matrix have 𝑁×𝑀 elements, which are divided into three subblocks in accordance with the number of orbit planes, and each subblock is to be allocated according to the number of satellites in each orbital plane. The matrix columns hold the same meaning as rows. The connection between topological matrix and global constellation ISL are as follows: suppose the 𝑗₁th satellite in orbital plane 𝑖₁ established a ISL with the 𝑗₂th satellite in orbital plane 𝑖₂; then the element 𝑖₁−1×𝑀+𝑗₁,𝑖₂−1×𝑀+𝑗₂ (1) of the matrix is 1; the others are 0. According to (1), ISL has another representation of 𝑖₂−1×𝑀+𝑗₂,𝑖₁−1×𝑀+𝑗₁; (2) that is (1) and (2) have the same meaning and mutually matrix transpose. In this paper, for the consideration of conceptual clarity, we use the upper triangular matrix representation, defined as Topological Method one (TM1). That is, the element of min𝑖₁−1×𝑀+𝑗₁,𝑖₂−1×𝑀+𝑗₂max𝑖₁−1×𝑀+𝑗₁,𝑖₂−1×𝑀+𝑗₂ (3) in the matrix is 1; the others are 0.

Based on the definition of the matrix elements above, the following conclusions can be drawn:(1)the diagonal elements of matrix denote each satellite itself, as the blue region in Figure 2;(2)according to the characteristics of the global constellation ISL, the topological matrix can be blocked to some regions, as shown in Figure 2. Area “1” indicates the links between the first orbit plane of global constellation; area “2” indicates the links between the second orbit plane of global constellation; also, area “3” represents the third orbit plane as before; areas “4, 5, and 6,” respectively, denote the links between the orbit plane 1 and 2, 2 and 3, and 3 and 1 of global constellation.

Figure 2 

Demonstration of topological block matrix.

Thus, when the 𝑗₁th satellite in orbit plane 𝑖₁ established an ISL link with the 𝑗₂th satellite in orbit plane 𝑖₂ in the constellation, according to the definition of (3), it can immediately get a topology matrix, as Figure 3 shows. Moreover, the third dimension of the ISL topology matrix denotes the time axis, as in Figure 4.

Figure 3 

Demonstration of 3-orbit plane ISL topological.

Figure 4 

Demonstration of 3-dimensional ISL topological matrix.

In the subsequent content of this paper, for the sake of matrix operations efficiency, some minor changes could be adopted for ISL topology matrix. Supposing the ISL matrix is 𝐴 at one moment, and we call it as Topological Method one (TM1), it could also be given as:(i)Topological Method two (TM2): 𝐵=𝐴+𝐴^𝑇;(ii)Topological Method three (TM3): 𝐶=𝐴+𝐴^𝑇+𝐼.

It can be proved that those three methods above hold the equivalent meaning (which will be proved in Section 3.2).

2.2. Some Characteristics of ISL Matrix

Definition 1. For a typical Walker 𝑁_𝑤/𝑀_𝑤/𝑃_𝑤 constellation, by using the agile scanning beam ISL design, one traversal cycle of constellation is defined as each satellite establish a link with other satellites once during a cycle period. The upper triangular elements of ISL topology matrix will be fully filled with ones by timing accumulate.

Proof. According to the definition of ISL topology matrix, under any epoch 𝑘, the ISL topology matrix has elements 𝑎_𝑖𝑗=1, which denotes the establishment of ISL between some specific pair of satellites. By timing accumulation of one traversal cycle, if the ISL topology matrix in a constellation appears, 𝑎_𝑖𝑗≠1𝑖,𝑗=2,3,…,𝑁_𝑤×𝑀_𝑤,𝑖≠𝑗. (4)
Then it indicates that a particular pair of satellites ISL disconnected or redundant, by the timing accumulation of the ISL topology matrix, the upper triangular elements would be fully filled with ones.

Lemma 2. ISLs between different orbit planes are established in pairs.

Proof. For a typical Walker 𝑁_𝑤/𝑀_𝑤/𝑃_𝑤 constellation, suppose the Right Ascension of Ascending Node (RAAN) and the mean anomaly of the 𝑗th satellite in orbital plane 𝑖 are Ω_𝑖𝑗=2𝜋𝑀_𝑤𝑀_𝑖−1,𝑀_𝑖=1,2,…,𝑀_𝑤,𝑢_𝑖𝑗=mod2𝜋𝑁_𝑤/𝑀_𝑤𝑗−1+2𝜋𝑀_𝑤𝑃_𝑤𝑀_𝑖−1,2𝜋, (5) where mod denotes the remainder operation.
Note that all the satellites in a Walker constellation are geometrical symmetry, and any satellite in the constellation can geometrically be representative of all others. Let us take a look at the first satellite in first orbit plane (𝑖=1,𝑗=1), subsequently referred to as the reference satellite. The phase difference of the 𝑗₁th satellite in the second orbit plane relative to the reference satellite is Δ𝑢_𝑗1=2𝜋𝑁_𝑤/𝑀_𝑤𝑗₁−1+2𝜋𝑀_𝑤𝑃_𝑤. (6)
Also, the phase difference of the 𝑗₂th satellite in the orbit plane 𝑀_𝑤 relative to the reference satellite is Δ𝑢_𝑗2=2𝜋𝑁_𝑤/𝑀_𝑤𝑗₂−1+2𝜋𝑀_𝑤𝑃_𝑤𝑀_𝑤−1. (7)
Obviously, when Δ𝑢_𝑗1+Δ𝑢_𝑗2=𝑗2𝜋¹+𝑗₂−2𝑁_𝑤/𝑀_𝑤+𝑃_𝑤 (8) holds integer reflects the 𝑗₁th and the 𝑗₂th satellite in different orbit plane which establish the geometrically equal ISLs to reference satellite. That is, ISLs between different orbit planes established in pairs.

Conclusion 1. For a typical Walker 𝑁_𝑤/𝑀_𝑤/𝑃_𝑤 constellation, by using one agile scanning beam at each satellite, the overall hop number of one traversal cycle is 𝑀_𝑤𝑀_𝑤−12𝑁_𝑤𝑀_𝑤+𝛼, (9) where 𝑁𝛼=^𝑤𝑀_𝑤,𝑁_𝑤𝑀_𝑤𝑁=Oddnumber_𝑤𝑀_𝑤𝑁−1,^𝑤𝑀_𝑤=Evennumber. (10)

Proof. The hop number for agile scanning beam analysis is divided into two different types: interplane (ISL’s between satellites in adjacent/corotating planes) and intraplane (ISL’s connecting satellites in the same orbit plane) [1].(1)For the orbit plane number 𝑀_𝑤, due to the ISLs between different orbit plane in pairs, the total combination of interplane ISLs gives 𝐶²_𝑀𝑤=𝑀_𝑤!2!𝑀_𝑤!=𝑀−2^𝑤𝑀_𝑤−12. (11) The ISLs between each pair of interplane have situations of (𝑁_𝑤/𝑀_𝑤) types; then the interplane ISLs have number of 𝑀_𝑤𝑀_𝑤−12𝑁_𝑤𝑀_𝑤 (12) hops.(2)For the intraplane ISLs, each orbital plane has satellites of (𝑁_𝑤/𝑀_𝑤), due to the ISLs between satellites in coorbit plane in pairs. Then at any epoch in a traversal cycle, one has the following:(a)coorbit plane with odd number satellites: each orbit plane has one satellite bye, traversing the hop count 𝑁_𝑤/𝑀_𝑤−12𝑁×2+2=^𝑤𝑀_𝑤 (13)(b)coorbit plane with even number satellites: each orbit plane get the full pairing, traversing the hop count 𝑁_𝑤𝑀_𝑤−1. (14) Then the overall hop number of one traversal cycle is 𝑀_𝑤𝑀_𝑤−12𝑁_𝑤𝑀_𝑤+𝛼. (15)

Corollary 3. According to Conclusion 1 and Figure 2, for a typical Walker 𝑁_𝑤/𝑀_𝑤/𝑃_𝑤 constellation, by using one agile scanning beam ISL at each satellite, the upper triangular elements of topology matrix will accumulating filled with ones during one ISL traversal cycle. The interplane matrix area of Figure 2 will be completely filled and the number of elements is 𝑀_𝑤𝑀_𝑤−12𝑁_𝑤𝑀_𝑤². (16)
Also, the intraplane matrix area of Figure 2 will be completely filled with number of elements 𝑁_𝑤2𝑁_𝑤𝑀_𝑤−1. (17)

Proof. (1) For interplane ISLs, each orbit plane has (𝑁_𝑤/𝑀_𝑤) satellites. It will be (𝑁_𝑤/𝑀_𝑤) ISLs between different orbit planes at one epoch. According to the hop number of interplane ISLs during one ISL traversal cycle from (12), the interplane ISLs has the number of 𝑀_𝑤𝑀_𝑤−12𝑁_𝑤𝑀_𝑤². (18)
Considering the geometric relationships of interplane area of Figure 2, the number of matrix elements of the interplane area is multiplying elements number of each square by number of squares: 𝑁_𝑤𝑀_𝑤²×𝑀_𝑤𝑀_𝑤−12. (19)
Equation (18) = (19); interplane ISLs proved.
(2) For intraplane ISLs, each orbital plane has (𝑁_𝑤/𝑀_𝑤) satellites; due to the paired establish of ISLs between satellites in coorbit, each epoch will have 𝑁floor^𝑤2𝑀_𝑤 (20) combinations, wherein floor is a numeric integer arithmetic.
According to the orbit planes 𝑀_𝑤 and hop number of coorbit 𝛼, the coorbit ISLs have the hop number of 𝑁floor^𝑤2𝑀_𝑤𝑀_𝑤𝛼. (21)
An alternate way to prove it is as follows: each orbital plane has (𝑁_𝑤/𝑀_𝑤) satellites, due to the paired establish of ISLs between satellites in coorbit, the coorbit ISLs number at each epoch have 𝐶²_{𝑁𝑤/𝑀𝑤}=𝑁_𝑤/𝑀_𝑤!2!𝑁_𝑤/𝑀_𝑤!=𝑁−2^𝑤/𝑀_𝑤𝑁_𝑤/𝑀_𝑤−12=𝑁_𝑤2𝑀_𝑤𝑁_𝑤𝑀_𝑤−1 (22) combinations; considering the orbit plane number 𝑀_𝑤, the coorbit ISLs have the total number of hops as 𝑁_𝑤2𝑀_𝑤𝑁_𝑤𝑀_𝑤−1⋅𝑀_𝑤=𝑁_𝑤2𝑁_𝑤𝑀_𝑤−1. (23)
This also proved that the formula (21) equals formula (23).

According to the geometric relationship of coorbit areas in Figure 2, the number of matrix elements in coorbit areas equals the upper triangular matrix elements of each block multiplying the number of matrix 12×𝑁_𝑤𝑀_𝑤×𝑁−11+^𝑤𝑀_𝑤−1×𝑀_𝑤=𝑁_𝑤2𝑁_𝑤𝑀_𝑤−1. (24)

Equation (23) = (24); coorbit ISLs is proved.

Conclusion 2. For a given ISL topology matrix 𝐴(𝑘) at epoch 𝑘, it always can be (1) transformed into subunit matrixes with same column sequence using 𝐴(𝑘)=𝐴(𝑘)^𝑇𝐴(𝑘) or??(2) transformed into subunit matrixes with same row sequence of using 𝐴(𝑘)=𝐴(𝑘)𝐴(𝑘)^𝑇.

Proof. The ISL topology matrix 𝐴(𝑘) at any epoch 𝑘 is upper triangular, according to the difference of orbit plane, the following general form of matrix block can be obtained (omitting epoch symbols 𝑘): 𝐴𝐴=¹¹𝐴₁₂𝐴₁₃⋯𝐴_1𝑁𝐴₂₂𝐴₂₃⋯𝐴_2𝑁𝐴₃₃⋯𝐴_3𝑁𝐴⋱⋮_𝑁𝑁, (25) where 𝑁 is the orbit plane number of constellation. It can be seen from above that 𝐴_𝑖𝑖,𝑖=1,2,3…𝑁 denote the coorbit ISLs topology matrix, and the others are different orbit ISLs topology matrix. According to the definition of topology matrix, the different orbit ISLs submatrix 𝐴_𝑖𝑗 can be written as the following general form: 𝐴_𝑖𝑗=𝑎_𝑝+1𝑎_𝑝+2⋱𝑎_𝑞𝑎₁𝑎₂⋱𝑎_𝑝, (26) where 𝑞=𝑀 represents the number of satellites at each orbit plane. By using matrix multiplication, the 𝐴(𝑘)^𝑇𝐴(𝑘) at each epoch results in 𝐴^𝑇_𝑖𝑗𝐴_𝑖𝑗=𝑎²_𝑝+1𝑎²_𝑝+2⋱𝑎²_𝑞𝑎²₁𝑎²₂⋱𝑎²_𝑝. (27)
Since each of its elements are all ones, so 𝑎²_𝑖=1,𝑖=1,2,…,𝑝,𝑝+1,…,𝑞. (28)
Then the subunit matrixes with same column sequence can be obtained. For the coorbit ISL submatrix 𝐴_𝑖𝑖, according to the ISL topology matrix definition, it can be written as the following general form: 𝐴_𝑖𝑗=0𝑎_𝑝+10𝑎₁0⋱⋱𝑎_𝑝00𝑎_𝑝+20⋱⋱𝑎_𝑞0. (29)
The result of the multiplication between any two submatrices 𝐴^𝑇_𝑖𝑖𝐴_𝑖𝑖 at each epoch is similar as (27). Then the ISL topology matrix 𝐴(𝑘) can be transformed into subunit matrixes with same column sequence using 𝐴(𝑘)^𝑇𝐴(𝑘).

Proposition 2 can be proved similar to proposition 1, which is not provided here.

Conclusion 3. During a traversal cycle of constellation using one ISL agile scanning antenna beam, the cumulative matrix 𝐴_batch=^𝑛_𝑘=1𝐴(𝑘) (30) has the form of diagonal; the diagonal elements of the matrix and the eigenvalues are diag𝐴_batch=[],𝜆𝐴0,1,2,3,…,𝑛−1_batch=[0,1,2,3,…,𝑛−1], (31) where 𝑛 is the hop number during one traversal cycle using one agile scanning beam.

Proof. The matrix 𝐴(𝑘) is the subunit matrix of 𝐴(𝑘) with the same column sequence at arbitrary epoch 𝑘, and the accumulation of the matrix 𝐴(𝑘) will fully fill the square with dimension of 𝑀𝑁×𝑀𝑁 in one traversal cycle of constellation. So, the accumulative matrix 𝐴_batch is substantially the accumulate of columns, with the form of diagonal matrix and diagonal elements diag(𝐴_batch)=[0,1,2,3,…,𝑛−1]. For the matrix in form of 𝐴_batch, the eigenvalues can be easily obtained as 𝜆(𝐴_batch)=[0,1,2,3,…,𝑛−1].

3. Some Applications

3.1. Analysis of ISL Hops

This section presents the analysis of optimal coorbit/different orbit ISL hops and the total ISL numbers under different satellite constellation conditions. It is assumed that each satellite has the capability of establishing only one agile scanning beam ISL during simulation. The selected constellations include practical engineering constellations and the theoretical ones as in Table 1.

Table 1 

Analysis of ISL hops and link number.

Here we take the Walker77/7/1 constellation as an example to demonstrate the computation method of total ISL hops and the meaning of each item in Table 1, as shown in the following: (32)

3.2. Cases of ISL Failure and Redundancy

It will be inevitable that in the event of ISLs failure during GNSS operation, supposing that one satellite failure occurred, the row and column 𝑖,𝑗 of the topology matrix will be vacant as shown in Figure 5. This section provides the analysis of ISL reconfiguration and restructure in this case, as well as ISL redundancy situation.

Figure 5 

Demonstration of one ISL failure.

First of all, here we give a brief introduction of some properties of the topology matrix and linear space: Let 𝐑 be the real number region, with elements in lowercase letters 𝑎,𝑏,𝑐…. On 𝐑; a collection of all 𝑛×𝑛 matrixes constitute a linear space on the basis of matrix addition and matrix multiplication, which is denoted as 𝐑^𝑛×𝑛.

It can be proved that one group of bases in space 𝐑^𝑛×𝑛 is 𝐄_𝑖𝑗=0⋮00⋮0(𝑖),𝑗0⋯010⋯0 (33) where 𝑖=1,2,…,𝑛;𝑗=1,2,…,𝑛.

Obviously, the Topological Method one is a spanned upper triangular matrix linear space on the base vectors 𝐄_𝑖𝑗,𝑖,𝑗satisfy(3). (34)

In case of ISL failure, as shown in Figure 5, the row and column 𝑖,𝑗 of the matrix vanished, some special treatment is needed in ISLs reorganization.

Here we assumed that the mapping 𝑓 is a linear transformation of linear space 𝐑^𝑛×𝑛 on real number region 𝐑; this mapping reflecting the transform from ISLs is intact to ISLs failure recombination. Based on the analysis in this section, taking a base 𝐄_𝑖𝑗 from linear space 𝐑^𝑛×𝑛, let 𝜶 be the ISL traversal vector in 𝐑^𝑛×𝑛; 𝜶=𝑥₁𝐄₁₁+𝑥₂𝐄₁₂+⋯+𝑥_𝑛𝐄_𝑖𝑗. (35)

Then we have 𝑓(𝜶)=𝑦₁𝐄₁₁+𝑦₂𝐄₁₂+⋯+𝑦_𝑛𝐄_𝑖𝑗. (36)

The question now is how to calculate the coordinates (𝑦₁,𝑦₂,…,𝑦_𝑛)^𝑇 of the base 𝐄_𝑖𝑗 under 𝑓(𝜶)? Here we write linear transformation as follows: 𝑓𝜶₁=𝑎₁₁𝜶₁+𝑎₂₁𝜶₂+⋯+𝑎_𝑛1𝜶_𝑛𝑓𝜶₂=𝑎₁₂𝜶₁+𝑎₂₂𝜶₂+⋯+𝑎_𝑛2𝜶_𝑛⋮𝑓𝜶_𝑛=𝑎_1𝑛𝜶₁+𝑎_2𝑛𝜶₂+⋯+𝑎_𝑛𝑛𝜶_𝑛, (37) where (𝑎_1𝑗,𝑎_2𝑗,…,𝑎_𝑛𝑗)^𝑇 are the coordinates of 𝑓(𝜶_𝑗) on the base 𝐄_𝑖𝑗.

Let 𝐴=𝑎₁₁𝑎₂₁⋯𝑎_𝑛1𝑎₁₂𝑎₂₂⋯𝑎_𝑛2𝑎⋮⋮⋱⋮_1𝑛𝑎_2𝑛⋯𝑎_𝑛𝑛. (38)

This matrix is the linear transformation on the base 𝐄_𝑖𝑗, wherein the column 𝑗 elements are the coordinates of 𝑓(𝜶_𝑗) on the base 𝐄_𝑖𝑗.

The analysis of ISLs redundancy situation and restructure is similar as above, which is not provided here.

Corollary 4. Topological Methods one/two/three (TM 1/2/3) are equal.

Proof. In linear space 𝐑^𝑛×𝑛, assume that the vector 𝜶∈𝐑^𝑛×𝑛 denotes the coordinates of the matrix 𝐴 in base 𝐄_𝑖𝑗 by using TM 1. Here we define the transform mapping from TM 1 to TM 2, and TM 1 to TM 3 are 𝑓₁₂=𝜶+𝜶^𝑇,𝑓₁₃=𝜶+𝜶^𝑇+𝐼. (39)
For any vector 𝜶,𝜷∈𝐑^𝑛×𝑛, we have 𝑓₁₂(𝜶+𝜷)=(𝜶+𝜷)+(𝜶+𝜷)^𝑇=𝜶+𝜶^𝑇+𝜷+𝜷^𝑇=𝑓₁₂(𝜶)+𝑓₁₂(𝜷),𝑓₁₃(𝜶+𝜷)=(𝜶+𝜷)+(𝜶+𝜷)^𝑇+2𝐼=𝜶+𝜶^𝑇+𝐼+𝜷+𝜷^𝑇+𝐼=𝑓₁₃(𝜶)+𝑓₁₃(𝜷). (40)
For any 𝑎∈𝐑, we have 𝑓₁₂(𝑎𝜶)=𝑎𝜶+𝑎𝜶^𝑇=𝑎𝜶+𝜶^𝑇=𝑎𝑓₁₂(𝜶),𝑓₁₃(𝑎𝜶)=𝑎𝜶+𝑎𝜶^𝑇+𝑎𝐼=𝑎𝜶+𝜶^𝑇+𝐼=𝑎𝑓₁₃(𝜶). (41)
Clearly, mapping 𝑓₁₂,𝑓₁₃ is the linear transformation in linear space 𝐑^𝑛×𝑛. Meanwhile, for a set of specific selected base 𝐄_𝑖𝑗, a unique matrix can be obtained using linear transformation 𝑓₁₂,𝑓₁₃ in 𝐑^𝑛×𝑛. Therefore, we create a one-to-one relationship between the matrix and the linear transformation; that is, TM 1/2/3 are equal.

3.3. Analysis of Computational Efficiency [9]

Commonly used computation efficient evaluation algorithm on modern computers is flop (floating octal point); one flop is a floating-point operation. Matrix flops calculation can be realized by adding the basic operation loops in one algorithm. In the cases of matrix multiplication and addition, the basic algorithm is 𝑎_𝑖𝑘𝑏_𝑘𝑗+𝑐_𝑖𝑗, which has two flops. This calculation requires 𝑚𝑛𝑝 times which were obtained by a simple loop count, so, generally the calculation of the matrix multiplication and addition needs 2𝑚𝑛𝑝 flops.

In this paper, the matrix topology representation is in the form of upper triangular matrix; here we give the calculation efficiency analysis.

For any upper triangular matrix 𝐴,𝐵∈𝐑^𝑛×𝑛, the product 𝐶=𝐴𝐵 is also upper triangular, that is, for any 𝑘<𝑖 or 𝑗<𝑘, there are 𝑎_𝑖𝑘𝑏_𝑘𝑗=0, so 𝑐_𝑖𝑗=^𝑗_𝑘=𝑖𝑎_𝑖𝑘𝑏_𝑘𝑗. (42)

Then, calculations of 𝑐_𝑖𝑗(𝑖≤𝑗) require 2(𝑗−𝑖+1) flops. The amount of the product calculation of two matrices is ^𝑛𝑛_𝑖=1𝑗=12=𝑗−𝑖+1_{𝑛𝑖=1𝑛−𝑖+1𝑗=1}≈2𝑗_𝑛𝑖=12𝑗−𝑖+1²2=^𝑛_𝑖=1𝑖²≈𝑛³3. (43)

That is, the amount of calculation for upper triangular matrix multiplication is only 1/6 of the general dense matrix, which clearly has higher computational efficiency.

4. 3D Topology Matrix Based ISL Routing Problem

Currently, ISL intersatellite communication has been widely used in new generation of global navigation satellite system, and operational data can be exchanged through ISL multihops [10–13]. It has obvious advantages in improving communication performance, reducing communication costs, and so forth. At the same time, due to the significant reduction in transmission delay and the increases in transmission data rate, it is particularly suitable for data exchange using ISL carrying a large amount of information including voice and multimedia transmission.

However, unlike the traditional network system, the ISL nodes in the network are moving constantly, and the ISLs on-off switch frequently with the change of relative position between satellites. The satellite in the network needs to establish a new link when the original one disconnects to accomplish the data transmission mission. Therefore, finding an optimal ISL routing is a crucial task in the GNSS ISL network design, and it directly affects the ISL network performance.

Previous research results in ISL routing problem include the following. Harathi et al. [14] provided an autonomous linking assignment algorithm based on the Line-Of-Sight (LOS) visibility of the satellites in constellation. Literatures [15, 16] also proposed some link assignment and routing algorithms with object of network robustness. Moreover, by using the theory of Finite State Automation (FSA), the link assignment problem in LEO satellite networks can be treated as a set of link assignment problems in fixed topology networks [17]. However, the optimization problems above are used for performance improvement for satellite communication networks. The optimal ISL routing algorithm designed for GNSS ISL network which is used for both intersatellite measurement and communication applications has little results published.

Here we provide a tailored optimal link assignment and routing problem for a GNSS ISL system below.

An optimal problem includes two parts: optimal object and constraints.

Previous research results about optimal ISL routing problems have been focused on the object of short path for end-to-end communications, including the critical up/downlink between space segment (satellites) and earth segment (mobile users and gateway stations), as well as the routing protocol model [1, 2]. However, those literatures are particularly used for analysis of ISL network for communication satellites, and it is not suitable for GNSS ISL network. A functional GNSS ISL network is mainly severing for the whole system stable operation and critical data dissemination; that is, at any time, the data from control segment on ground can achieve constellation traversal by visible satellites and minimum number of ISL relay between satellites. Therefore, this paper considered the minimum ISL relay hops between satellites as the optimal routing object which is used for the data upload from ground station to whole constellation.

Some constraints should be considered in an optimal ISL routing problem such as linking on-off switch conditions, linking propagation distance, and linking time. Those of constraints directly affect the ISL network performance. Generally, ISLs establishment between satellites in constellation relies on two factors: (1) the couple satellites are under each other’s ISL antenna beam coverage; (2) the projection of relative velocity on the direction of LOS between couple satellites is within the given range due to the payload specification. Based on the analysis above, here we propose the optimal ISL routing problem expressed as follows:?Finding?min?(relay hops)?such as?(1) ISL antenna beam coverage; (2) relative velocity constraint.

The constraint models in the optimal ISL routing problem are introduced in detail hereafter, and the overall optimization process will finally be provided.

4.1. Analysis of Antenna Beams

This section provides the analysis of the antenna beam coverage constraints. First, definition of elevation angle and related physical quantities are given; then, ISL visibility criterion between couple satellites is introduced.

Figure 6 shows the geometry quantity associated with the ISL elevation angle. 𝛼 is vertical plane perpendicular to the line 𝑂𝐴 which is from satellite 𝐴 to the center of the earth 𝑂. The elevation angle 𝐸 that is from satellite 𝐴 to satellite 𝐵 is defined as the angle between line 𝐴𝐵 and plane 𝛼. The center of the sphere angle 𝜃_𝐴𝐵 is called ISL angle (geocentric angle).

Figure 6 

Geometry of ISL elevation angle.

According to the geometrical relationship in Figure 6 and the definition of ISL elevation angle, we can see that if the LOS of two satellites is not obscured by the earth, the elevation angle should satisfy 𝜋sin2≥𝑅−𝐸^𝑒𝑅_𝑒.+ℎ (44)

Generally, the elevation angle is defined in interval of [−90deg,90deg]; formula (44) can be solved as 𝐸≤𝐸_max=𝜋2𝑅−arcsin^𝑒𝑅_𝑒.+ℎ (45)

Since the circular orbit approximation of GNSS can be adopted, then we have |𝑂𝐴|=|𝑂𝐵|. According to the spatial geometry relationship as in Figure 6, one can obtain the following: 𝐸(𝑡)=𝜃_𝐴𝐵(𝑡)2. (46)

Then the ISL antenna beam coverage condition of two satellites is 2𝐸_min≤𝜃_𝐴𝐵(𝑡)≤2𝐸_max. (47)

The analysis of ISL antenna coverage above is only of the schematic formulas, which cannot be used directly in practice. Here we provide the analytical analysis of specific visible satellites’ mean anomaly under the conditions of ISL antenna beam coverage.

Due to the relative motion between satellites in different orbital plane in constellation, the specific time of ISL visibility depends on the constellation configuration and time 𝑡; here we provide an analytical method to calculate visibility according to time 𝑡.

Figure 7 shows the geometric relationship between satellites in different orbital plane. Orbit plane 1 and plane 2 intersect at point 𝑃 in the northern hemisphere. Considering a polar coordinate with axis 𝑂𝑃 that is from center of the earth 𝑂 to polar point 𝑃, the positive direction of the polar coordinate is complied with the satellite motion direction (as shown by the arrows in Figure 7). Obviously, the angle of the satellite in this coordinate is called mean anomaly. Note that the ascending node angle of the orbit plane in the polar coordinate is 𝜂, which is a determined parameter by the constellation configuration. According to the GNSS ephemeris parameters, we can see that the mean anomaly 𝜑(𝑡) of satellite at a time 𝑡 is 𝜑(𝑡)=𝑡−𝑡_ref𝜇𝑎³+𝜔+𝜂, (48) where 𝜇 is the product of the gravitational constant and the mass of Earth, 𝑡_ref is reference time of ephemeris, a is orbital semimajor axis, 𝜔 is perigee at time 𝑡_ref, and 𝜂 is the polar coordinate angle of ascending node in the corresponding orbit plane.

Figure 7 

The ISL antenna beam coverage analysis.

The shaded spatial region of Figure 7 illustrates the ISL antenna scan area (visible area) of satellite 𝑆 in orbit plane 1. Arc 𝐴𝐵 and 𝐶𝐷 of orbital plane 2 are under the visible area of the satellite 𝑆, called visual arcs hereafter, which can be represented by the corresponding geocentric angle. Bellowing provide the analysis of four angles of point 𝐴,𝐵,𝐶,𝐷 in polar coordinates in orbit plane 2.

Establish Cartesian Coordinates which are centered at the Earth, as in Figure 7. The fundamental plane is orbit plane 1; the vector 𝑂𝑥 is directed from the Earth’s center to satellite 𝑃. 𝑂𝑧 is normal to the fundamental plane, positive in the direction of the reader, and 𝑂𝑦 completes the setup. Then the coordinates of point 𝑆 are 𝑆=𝑎cos𝜑_𝑆sin𝜑_𝑆0. (49)

The coordinates of any point in orbit plane 2 are 𝑇=𝑎cos𝜑_𝑇sin𝜑_𝑇cos𝛾sin𝛾. (50)

The geocentric angle ∠𝑆𝑂𝑇 satisfies cos𝜃_𝑇𝑆=cos𝜑_𝑆cos𝜑_𝑇+sin𝜑_𝑆sin𝜑_𝑇cos𝛾 (51) which can be solved as cos𝜃_𝑇𝑆=cos²𝜑_𝑆+sin²𝜑_𝑆cos²𝛾cos𝜑_𝑇−tansin𝜑_𝑆cos𝛾cos𝜑_𝑆. (52)

The mean anomaly of point 𝑇 is 𝜑_𝑇=tansin𝜑_𝑆cos𝛾cos𝜑_𝑆±arccoscos𝜃_𝑇𝑆cos²𝜑_𝑆+sin²𝜑_𝑆cos²𝛾. (53)

Accordingly, the boundary of ISL antenna beam coverage has four points 𝜑_𝐴=tansin𝜑_𝑆cos𝛾cos𝜑_𝑆+arccoscos2𝐸_mincos²𝜑_𝑆+sin²𝜑_𝑆cos²𝛾,𝜑_𝐵=tansin𝜑_𝑆cos𝛾cos𝜑_𝑆+arccoscos2𝐸_maxcos²𝜑_S+sin²𝜑_𝑆cos²𝛾,𝜑_𝐶=tansin𝜑_𝑆cos𝛾cos𝜑_𝑆−arccoscos2𝐸_mincos²𝜑_𝑆+sin²𝜑_𝑆cos²𝛾,𝜑_𝐷=tansin𝜑_𝑆cos𝛾cos𝜑_𝑆−arccoscos2𝐸_maxcos²𝜑_𝑆+sin²𝜑_𝑆cos²𝛾. (54)

Thus, when the mean anomaly 𝜑(𝑡)∈[0,𝜋], the visible arcs of orbit plane 2 from satellite 𝑆 are [0,2𝜋]𝜑(𝑡)∈2𝐸_min,𝜋−2𝐸_min𝐼𝜋−2𝐸_max,2𝐸_max,𝐸_𝐴,𝐸_𝐵∪2𝜋−𝐸_𝐶,2𝜋−𝐸_𝐷𝜑(𝑡)∉2𝐸_min,𝜋−2𝐸_min,𝜑(𝑡)∉𝜋−2𝐸_max,2𝐸_max,𝐸_𝐴,2𝜋−𝐸_𝐷𝜑(𝑡)∉2𝐸_min,𝜋−2𝐸_min,𝜑(𝑡)∈𝜋−2𝐸_max,2𝐸_max,0,𝐸_𝐵∪2𝜋−𝐸_𝐶,2𝜋𝜑(𝑡)∈2𝐸_min,𝜋−2𝐸_min,𝜑(𝑡)∉𝜋−2𝐸_max,2𝐸_max. (55)

Similarly, the visible arcs of orbit plane 2 from satellite 𝑆 can be obtained when the mean anomaly 𝜑(𝑡)∈[𝜋,2𝜋], which are not provided here for short.

Based on the analysis above, for a typical Walker 𝑁_𝑤/𝑀_𝑤/𝑃_𝑤 constellation, let 𝑆_𝑖𝑗 denote the 𝑗th satellite in orbital plane 𝑖, if the mean anomaly of 𝑆_𝑖𝑗 is known, we can determine any satellite 𝑆_𝑘𝑙 in other orbital plane according to the constellation geometry. Suppose we have mean anomaly 𝜑_𝑖𝑗(𝑡) of satellite 𝑆_𝑖𝑗; then the angles 𝜑_𝑖1(𝑡), 𝜑_𝑘1(𝑡), 𝜑_𝑘𝑙(𝑡) of satellites 𝑆_𝑖1(𝑡), 𝑆_𝑘1(𝑡), 𝑆_𝑘𝑙(𝑡) are 𝜑_𝑖1(𝑡)=mod𝜑_𝑖𝑗(𝑡)−2𝜋𝑁_𝑤/𝑀_𝑤,𝜑𝑗−1,2𝜋_𝑘1(𝑡)=mod𝜑_𝑖1(𝑡)+2𝜋𝑀_𝑤𝑃_𝑤(),𝜑𝑘−𝑖,2𝜋_𝑘𝑙(𝑡)=mod𝜑_𝑘1(𝑡)+2𝜋𝑁_𝑤/𝑀_𝑤(𝑙−1),2𝜋, (56) where 𝑗,𝑙=2,3,…,𝑁_𝑤/𝑀_𝑤.

In summary, the procedures of calculating visible satellites in orbit plane 𝑘 from satellite 𝑆_𝑖𝑗 at time 𝑡 are as follows:(1)calculate mean anomaly of 𝑆_𝑖𝑗 by formula (48);(2)calculate the visible arcs of orbit plane 𝑘 from satellite 𝑆_𝑖𝑗 by the formula (55);(3)calculate the mean anomaly 𝜑_𝑘𝑙(𝑡) of each satellite in orbit plane 𝑘 at time 𝑡 by formula (56);(4)collect all the visible satellites according to the analysis of visible arcs and mean anomaly of each satellite.

4.2. Analysis of Relative Velocity

Assume that satellite 𝑚 is in constellation with inclination 𝑖, mean anomaly 𝑢_𝑚, and RAAN Ω_𝑚. The differential orbital element parameters with any other satellite 𝑛 in constellation are Δ𝑢_𝑚𝑛,ΔΩ_𝑚𝑛. Then we can obtain the projection of the relative velocity on the direction of LOS as 𝑣_𝑟=−𝑣sin2𝑢^′_𝑚+Δ𝑢^′_𝑚𝑛1−cosΔΩ^′_𝑚𝑛, (57) where 𝑢^′_𝑚=𝑢_𝑚−arctantanΔΩ_𝑚𝑛2cos𝑖,Δ𝑢^′_𝑚𝑛=Δ𝑢_𝑚𝑛+2arctantanΔΩ_𝑚𝑛2cos𝑖,cosΔΩ^′_𝑚𝑛=cos²𝑖+sin²𝑖cosΔΩ_𝑚𝑛. (58)

Therefore, the second condition for ISL establish between two satellites, is 𝑣_𝑟<𝑣₀, (59) where 𝑣₀ is the relative velocity range which depends on the performance of satellite communication payload.

4.3. Summary of Optimal Routing Problem

Based on the analysis of object and constrains for optimal routing problem. The overall optimization process can now be expressed as in Figure 8.

Figure 8 

Flow diagram of optimal ISL routing problem.

5. Simulation of Optimal Routing Problem

According to the introduction of optimal ISL routing problem in Section 4, here we provide some numerical simulations in this section to verify the effectiveness of proposed optimization procedure in this paper.

Suppose we use Walker24/3/2 constellation in simulation, with orbital inclination of 55?deg, orbital altitude of 21528.8?km, intersection period of 12?h?53?min, and regression cycle of 13 laps/7 days. This constellation achieved availability of 65.18% and 94.52% with conditions of PDOP = 2 and positioning error less than 10?m, respectively, after statistics analysis. Moreover, the distribution difference of constellation’s largest PDOP around the world is small, and the average positioning accuracy is about 6?m, which is more practical for engineering use. Control segment is selected as Beijing Station with coordinates of latitude 39°54'27'' north, longitude 116°23'17'' east, and 5?deg elevation angle constraint when observing zenith.

The ISL antenna beam range is assumed to be 𝐸_min=15deg, 𝐸_max=70deg. The maximum relative velocity is set to 𝑣₀=±3km/s considering ISL’s communication payload.

The initial epoch of simulation is set to April 8, 2014, 141500 (GMT +08:00). Simulation duration is set to one regression period of 7 days.

The actual configuration of ISL antennas in each satellite is a practical engineering problem during the process of optimal ISL link assessment. Different optimal results will be obtained with different types of antenna configurations. Therefore, this section provides two types of simulation scenarios: one without ISL antenna configuration constraints and one with.

5.1. Simulation without ISL Antenna Configuration Constraint

Figures 9–13 show the simulation results without ISL antenna configuration constraint. Those figures are arranged to two columns: the actual space location of constellation at each moment (left (a)) combined with the topological matrix (right (b)). Figure 9(a) shows the position of each satellite in constellation at initial epoch (satellite numbers 1–24 are the batch serial sequence of each 8 satellites in 3 orbit plane); Figure 9(b) denotes the corresponding topological matrix. According to the optimization process, first of all, the ground station implement data uploading to visible satellite zenith, as shown in Figure 10(a). Then, those visible satellites as in part of the constellation have received the uploading data from ground, as the red dots shown in Figure 11(a). Also, the topological matrix has changed in some elements, as the red dots in Figure 11(b), and the red stars in Figure 11(b) denote the satellites in constellation that has not received data from ground yet.

(a)

(b)

Figure 9 

Constellation states at initial epoch.

(a)

(b)

Figure 10 

Data uploading to visible satellite from ground.

(a)

(b)

Figure 11 

Satellites that have received the uploading data.

(a)

(b)

Figure 12 

Path searching for possible ISL establishment.

(a)

(b)

Figure 13 

Satellites that have received data by relay process.

Thereafter, entering the phase of data relay, as shown in Figure 12(a), the uploading data from ground begin to implement the second transmit to the unreceived satellites. The dashed blue lines denote the every possible ISL routing to the visible satellites relying on relay satellites (red ones), and Figure 12(b) provided the corresponding topology matrix which also performs the ISL routing search. Eventually, the whole constellation network achieved uploading data traversal. The blue dots in Figure 13(a) denote the data relay receiving satellites which are invisible from ground station at first time, and Figure 13(b) illustrates the final topology matrix.

A variety of assess criterion could be adopted for the optimal ISL routing problem. Here we use the percentage of special time compared with total simulation time which achieved the uploading data from ground transmitted to whole constellation using just ONE relay. By using numerical statistics after simulation, the ISL routing can achieve 100% of the time that realize uploading data traversal using just ONE relay, under the condition of no ISL antenna configuration constraint; clearly, it is a pretty good result impersonality.

5.2. Simulation with ISL Antenna Configuration Constraint

The simulation scenario in Section 5.1 is based on the assumption of no ISL antenna configuration constraint, which is not practical for engineering applications. Due to the limitation of satellite platform equipment, the actual number of antenna beams is generally less than the possible ones. Then, proper selection should be implemented during the navigation constellation operation, and the assignment of ISL resources must be considered.

This section assumes that each satellite is configured with two agile scanning ISL antenna beams. In each epoch, those two beams select the geometrical optimal paths and establish ISL links. Other simulation scenarios and parameters are identical to the above section. The illustrations of topological matrixes are not provided here for conciseness.

Figures 14, 15, 16, 17, 18, and 19 show the simulation results of using just 2 agile scanning ISL beams in each satellite. The dashed blue lines in Figure 17 denote the path searching process relying on the visible satellites from ground station. The solid blue lines in Figure 18 represent the actual using paths (note only two ISL paths at most could be established relying on relay satellites, with the standards of optimal ISL geometry), and ultimately realized the uploading data traversal of the whole network. With numerical statistics, the ISL routing can achieve 89% of the time that realize uploading data traversal using just one relay, under the condition of 2 ISL antenna configuration constraint.

Figure 14 

Initial epoch.

Figure 15 

Data upload.

Figure 16 

Data received (1st).

Figure 17 

Path searching.

Figure 18 

Data relaying (1st).

Figure 19 

Data received (2nd).

5.3. Some Suggestions of Using Consistent ISL Beams

Generally, two fundamental demands should be satisfied when using ISL routing algorithm: ISL ranging measurements and communications between satellites. The former emphasize on the number of established ISLs. The more the ISLs measurements acquired, the more the accuracy results that will be got. The latter emphasize on the quality of established ISLs, such as propagation delay or the carrier-to-noise ratio. Therefore, how to efficiently use of the limited on-board antenna beams that take into account both demands is the difficulty of the optimal routing algorithm.

The link assignment and routing algorithms used in previous two simulation scenarios are based on the real-time constellation status and dynamic calculations. The ISL antennas have to capture and track the calculated couple satellite with speed and accuracy at each epoch, which is almost impossible to achieve in engineering. According to the cyclical nature of the Walker constellation, routing assignment calculated in advance and stored in each satellite cannot cope with the situations as satellite failures and system abnormalities. Therefore, seeking for the combination of consistent antenna beams with agility ones is a more feasible solution in practical.

The consistent ISL antenna beams are the links between two satellites, in the same or different orbit planes, that can be seen from each other with no block during the entire constellation regression cycle. Those links can be analytically calculated in advance. Those ISL beams are commonly realized by using reflector antenna in engineering. With proper configuration of ISL antennas on satellite platform and the usage of high reliability shaft, uninterrupted tracking can be achieved with established linking satellite and intersatellite ranging and data transfer can be performed.

The consistent ISL beams have merits of high precision ranging and high rate of data transfer speed, and so forth, which are the preferred solutions for practical use. On the other hand, usage of ISL agile scanning beams has to consider the time sequence design in advance, which provides poor support to real-time intersatellite data relay.

Figure 20 provided the simulation results of data relaying index under different conditions of consistent/agile ISL beams configurations in spacecraft. Suppose we have one satellite platform which allows up to four ISL antennas onboard. The other simulation scenarios and parameters are the same as in Section 5.1. The 𝑥-axis of Figure 20 denotes the number of agile scanning beams onboard, and the 𝑦-axis denotes the number of consistent ISL beams onboard. The 𝑧-axis denotes the percentage of total simulation time that realized uploading data traversal using just ONE relay from ground.

Figure 20 

Data relaying index under different conditions.

As we can see through the simulation data in Figure 20, under the same conditions, the more the ISL antenna beams equipped onboard of satellite platform, the better the relaying index. Such as under the same conditions of consistent ISL link number (as we can see those bars of same color in Figure 20), increasing the number of agile scanning beams can significantly improve relaying index, and vice versa. According to the illustrations in Figure 20, it is clear that four ISL antenna beams should be equipped on satellite platform, if one wishes to achieve the best relaying index. However, there are five options for four ISL beams configurations (as the sequential numeral of 1–5 in Figure 20):?option (1): 4 consistent linking beams;?option (2): 3 consistent linking beams +1 agile scanning beam;?option (3): 2 consistent linking beams +2 agile scanning beams;?option (4): 1 consistent linking beam +3 agile scanning beams;?option (5): 4 agile scanning beams.

According to the numerical statistics after simulation, the relaying indexes using those five options above are 93.1%; 98.1%; 81.2%; 73.1%; 82.3%. As we can see from the statistical data, the best relaying index can be obtained using option (2). The explanations can be given as follows: some agile scanning beams equipped onboard have merits of flexibility and effectively; however, proper time sequences have to be designed and assigned during the actual constellation network operation, which restricted the ISL network real-time relaying performance. Using the combination of consistent linking beams with agile scanning ones as option (2) can inspire the best performance of both antennas. Moreover, the whole ISL network exhibit robustness under conditions of ISL linking failure and redundancy, which can implement network reconfiguration nicely and quickly. Clearly, option (2) is a more feasible ISL network for actual use.

6. Conclusions

In this paper, a three-dimensional matrix topology representation method is proposed. This method compared with the previous ones, despite the increasing in matrix dimensions, possesses the merits of physical concept clarity, upper triangular matrix with sparse elements, smaller amount of computation, and the easier proof of ISL network properties. Moreover, an optimal link assignment and routing problem is provided, which considers the minimization of data relay times during the whole constellation operation. Some engineering problems as ISL antenna beam coverage and relative velocity constrains are also analyzed. A variety of scenarios are simulated for proposed three-dimensional topology matrix and routing problem. The results demonstrated that the proposed method has advantage of conceptual clarity and computational efficiency; moreover, the effectiveness of provided optimal ISL routing problem is also proved.

The development of ISL technology in future will be focused on the following aspects: supporting to the autonomous navigation and constellation system operation; supporting to the system differential and integrity monitoring; and critical data transmitting to the whole network under critical conditions, as well as the rapid restructuring of ISL network during satellite failure or redundant. The application of three-dimensional matrix, presented in this paper, to those areas will be provided in authors’ consequent researches.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.