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

Authors: Xiaoliang Wang; Yansong Meng; Zhe Su; Gang Wang; Lixin Zhang; Deren Gong
Journal: International Journal of Navigation and Observation
 (2014)

---

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


---

## Body

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


## 2. Matrix Model

### 2.1. The Elements of ISL Topology Matrix

Assuming that a typical global navigation constellation consists ofN orbital planes, each orbital plane have M satellites (to simplify the notation, this section consider a constellation of three orbital planes, with four satellites at each plane; that is N
=
3,  M
=
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 haveN
×
M 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 j
1th satellite in orbital plane i
1 established a ISL with the j
2th satellite in orbital plane i
2; then the element
(1)
[
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
]
of the matrix is 1; the others are 0. According to (1), ISL has another representation of
(2)
[
(
i
2
-
1
)
×
M
+
j
2
,
(
i
1
-
1
)
×
M
+
j
1
]
;
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
(3)
[
min
⁡
{
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
}
max
⁡
{
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
}
]
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 Figure2;
(2)
according to the characteristics of the global constellation ISL, the topological matrix can be blocked to some regions, as shown in Figure2. 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 thej
1th satellite in orbit plane i
1 established an ISL link with the j
2th satellite in orbit plane i
2 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 isA at one moment, and we call it as Topological Method one (TM1), it could also be given as:(i)
Topological Method two (TM2):B
=
A
+
A
T;
(ii)
Topological Method three (TM3):C
=
A
+
A
T
+
I.It can be proved that those three methods above hold the equivalent meaning (which will be proved in Section3.2).


### 2.2. Some Characteristics of ISL Matrix

Definition 1.
For a typical WalkerN
w
/
M
w
/
P
w 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 epochk, the ISL topology matrix has elements a
i
j
=
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,
(4)
a
i
j
≠
1
(
i
,
j
=
2,3
,
…
,
N
w
×
M
w
,
i
≠
j
)
.
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 WalkerN
w
/
M
w
/
P
w constellation, suppose the Right Ascension of Ascending Node (RAAN) and the mean anomaly of the jth satellite in orbital plane i are
(5)
Ω
i
j
=
2
π
M
w
(
M
i
-
1
)
,
M
i
=
1,2
,
…
,
M
w
,
u
i
j
=
mod
⁡
[
2
π
N
w
/
M
w
(
j
-
1
)
+
2
π
M
w
P
w
(
M
i
-
1
)
,
2
π
]
,
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 (i
=
1
,
j
=
1), subsequently referred to as the reference satellite. The phase difference of the j
1th satellite in the second orbit plane relative to the reference satellite is
(6)
Δ
u
j
1
=
2
π
N
w
/
M
w
(
j
1
-
1
)
+
2
π
M
w
P
w
.
Also, the phase difference of thej
2th satellite in the orbit plane M
w relative to the reference satellite is
(7)
Δ
u
j
2
=
2
π
N
w
/
M
w
(
j
2
-
1
)
+
2
π
M
w
P
w
(
M
w
-
1
)
.
Obviously, when(8)
Δ
u
j
1
+
Δ
u
j
2
2
π
=
(
j
1
+
j
2
-
2
)
N
w
/
M
w
+
P
w
holds integer reflects the j
1th and the j
2th 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 N
w
/
M
w
/
P
w constellation, by using one agile scanning beam at each satellite, the overall hop number of one traversal cycle is
(9)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
+
α
,
where
(10)
α
=
{
N
w
M
w
,
N
w
M
w
=
Odd
number
N
w
M
w
-
1
,
N
w
M
w
=
Even
number
.Proof.
The hop number for agile scanning beam analysis is divided into two different types:interplane (ISL’s between satellites in adjacent/corotating planes) andintraplane(ISL’s connecting satellites in the same orbit plane) [1].(1)
For the orbit plane numberM
w, due to the ISLs between different orbit plane in pairs, the total combination of interplane ISLs gives
(11)
C
M
w
2
=
M
w
!
2
!
(
M
w
-
2
)
!
=
M
w
(
M
w
-
1
)
2
.
The ISLs between each pair of interplane have situations of (
N
w
/
M
w
) types; then the interplane ISLs have number of
(12)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
hops.
(2)
For the intraplane ISLs, each orbital plane has satellites of(
N
w
/
M
w
), 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(13)
[
(
(
N
w
/
M
w
)
-
1
)
2
]
×
2
+
2
=
N
w
M
w
(b)
coorbit plane with even number satellites: each orbit plane get the full pairing, traversing the hop count(14)
N
w
M
w
-
1
.
Then the overall hop number of one traversal cycle is
(15)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
+
α
.Corollary 3.
According to Conclusion 1 and Figure2, for a typical Walker N
w
/
M
w
/
P
w 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
(16)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
2
.
Also, the intraplane matrix area of Figure2 will be completely filled with number of elements
(17)
N
w
2
(
N
w
M
w
-
1
)
.Proof.
(
1
) For interplane ISLs, each orbit plane has (
N
w
/
M
w
) satellites. It will be (
N
w
/
M
w
) 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
(18)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
2
.
Considering the geometric relationships of interplane area of Figure2, the number of matrix elements of the interplane area is multiplying elements number of each square by number of squares:
(19)
(
N
w
M
w
)
2
×
M
w
(
M
w
-
1
)
2
.
Equation (18) = (19); interplane ISLs proved.
(
2
) For intraplane ISLs, each orbital plane has (
N
w
/
M
w
) satellites; due to the paired establish of ISLs between satellites in coorbit, each epoch will have
(20)
floor
(
N
w
2
M
w
)
combinations, whereinfloor is a numeric integer arithmetic.
According to the orbit planesM
w and hop number of coorbit α, the coorbit ISLs have the hop number of
(21)
floor
(
N
w
2
M
w
)
M
w
α
.
An alternate way to prove it is as follows: each orbital plane has(
N
w
/
M
w
) satellites, due to the paired establish of ISLs between satellites in coorbit, the coorbit ISLs number at each epoch have
(22)
C
N
w
/
M
w
2
=
(
N
w
/
M
w
)
!
2
!
(
(
N
w
/
M
w
)
-
2
)
!
=
(
N
w
/
M
w
)
(
(
N
w
/
M
w
)
-
1
)
2
=
N
w
2
M
w
(
N
w
M
w
-
1
)
combinations; considering the orbit plane number M
w, the coorbit ISLs have the total number of hops as
(23)
N
w
2
M
w
(
N
w
M
w
-
1
)
·
M
w
=
N
w
2
(
N
w
M
w
-
1
)
.
This also proved that the formula (21) equals formula (23).According to the geometric relationship of coorbit areas in Figure2, the number of matrix elements in coorbit areas equals the upper triangular matrix elements of each block multiplying the number of matrix
(24)
1
2
×
(
N
w
M
w
-
1
)
×
(
1
+
N
w
M
w
-
1
)
×
M
w
=
N
w
2
(
N
w
M
w
-
1
)
.Equation (23) = (24); coorbit ISLs is proved.Conclusion 2. For a given ISL topology matrix A
(
k
) at epoch k, it always can be (
1
) transformed into subunit matrixes with same column sequence using A
-
(
k
)
=
A
(
k
)
T
A
(
k
) or  (
2
) transformed into subunit matrixes with same row sequence of using A
~
(
k
)
=
A
(
k
)
A
(
k
)
T.Proof.
The ISL topology matrixA
(
k
) at any epoch k is upper triangular, according to the difference of orbit plane, the following general form of matrix block can be obtained (omitting epoch symbols k):
(25)
A
=
[
A
11
A
12
A
13
⋯
A
1
N
A
22
A
23
⋯
A
2
N
A
33
⋯
A
3
N
⋱
⋮
A
N
N
]
,
where N is the orbit plane number of constellation. It can be seen from above that A
i
i
,
i
=
1,2
,
3
…
N 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 A
i
j can be written as the following general form:
(26)
A
i
j
=
[
a
p
+
1
a
p
+
2
⋱
a
q
a
1
a
2
⋱
a
p
]
,
where q
=
M represents the number of satellites at each orbit plane. By using matrix multiplication, the A
(
k
)
T
A
(
k
) at each epoch results in
(27)
A
i
j
T
A
i
j
=
[
a
p
+
1
2
a
p
+
2
2
⋱
a
q
2
a
1
2
a
2
2
⋱
a
p
2
]
.
Since each of its elements are all ones, so(28)
a
i
2
=
1
,
i
=
1,2
,
…
,
p
,
p
+
1
,
…
,
q
.
Then the subunit matrixes with same column sequence can be obtained. For the coorbit ISL submatrixA
i
i, according to the ISL topology matrix definition, it can be written as the following general form:
(29)
A
i
j
=
[
0
a
p
+
1
0
a
1
0
⋱
⋱
a
p
0
0
a
p
+
2
0
⋱
⋱
a
q
0
]
.
The result of the multiplication between any two submatricesA
i
i
T
A
i
i at each epoch is similar as (27). Then the ISL topology matrix A
(
k
) can be transformed into subunit matrixes with same column sequence using A
(
k
)
T
A
(
k
).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
(30)
A
batch
=
∑
k
=
1
n
A
-
(
k
)
has the form of diagonal; the diagonal elements of the matrix and the eigenvalues are
(31)
diag
⁡
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]
,
λ
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]
,
where n is the hop number during one traversal cycle using one agile scanning beam.Proof.
The matrixA
-
(
k
) is the subunit matrix of A
(
k
) with the same column sequence at arbitrary epoch k, and the accumulation of the matrix A
(
k
) will fully fill the square with dimension of M
N
×
M
N in one traversal cycle of constellation. So, the accumulative matrix A
batch is substantially the accumulate of columns, with the form of diagonal matrix and diagonal elements diag
⁡
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]. For the matrix in form of A
batch, the eigenvalues can be easily obtained as λ
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
].


## 2.1. The Elements of ISL Topology Matrix

Assuming that a typical global navigation constellation consists ofN orbital planes, each orbital plane have M satellites (to simplify the notation, this section consider a constellation of three orbital planes, with four satellites at each plane; that is N
=
3,  M
=
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 haveN
×
M 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 j
1th satellite in orbital plane i
1 established a ISL with the j
2th satellite in orbital plane i
2; then the element
(1)
[
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
]
of the matrix is 1; the others are 0. According to (1), ISL has another representation of
(2)
[
(
i
2
-
1
)
×
M
+
j
2
,
(
i
1
-
1
)
×
M
+
j
1
]
;
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
(3)
[
min
⁡
{
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
}
max
⁡
{
(
i
1
-
1
)
×
M
+
j
1
,
(
i
2
-
1
)
×
M
+
j
2
}
]
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 Figure2;
(2)
according to the characteristics of the global constellation ISL, the topological matrix can be blocked to some regions, as shown in Figure2. 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 thej
1th satellite in orbit plane i
1 established an ISL link with the j
2th satellite in orbit plane i
2 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 isA at one moment, and we call it as Topological Method one (TM1), it could also be given as:(i)
Topological Method two (TM2):B
=
A
+
A
T;
(ii)
Topological Method three (TM3):C
=
A
+
A
T
+
I.It can be proved that those three methods above hold the equivalent meaning (which will be proved in Section3.2).


## 2.2. Some Characteristics of ISL Matrix

Definition 1.
For a typical WalkerN
w
/
M
w
/
P
w 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 epochk, the ISL topology matrix has elements a
i
j
=
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,
(4)
a
i
j
≠
1
(
i
,
j
=
2,3
,
…
,
N
w
×
M
w
,
i
≠
j
)
.
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 WalkerN
w
/
M
w
/
P
w constellation, suppose the Right Ascension of Ascending Node (RAAN) and the mean anomaly of the jth satellite in orbital plane i are
(5)
Ω
i
j
=
2
π
M
w
(
M
i
-
1
)
,
M
i
=
1,2
,
…
,
M
w
,
u
i
j
=
mod
⁡
[
2
π
N
w
/
M
w
(
j
-
1
)
+
2
π
M
w
P
w
(
M
i
-
1
)
,
2
π
]
,
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 (i
=
1
,
j
=
1), subsequently referred to as the reference satellite. The phase difference of the j
1th satellite in the second orbit plane relative to the reference satellite is
(6)
Δ
u
j
1
=
2
π
N
w
/
M
w
(
j
1
-
1
)
+
2
π
M
w
P
w
.
Also, the phase difference of thej
2th satellite in the orbit plane M
w relative to the reference satellite is
(7)
Δ
u
j
2
=
2
π
N
w
/
M
w
(
j
2
-
1
)
+
2
π
M
w
P
w
(
M
w
-
1
)
.
Obviously, when(8)
Δ
u
j
1
+
Δ
u
j
2
2
π
=
(
j
1
+
j
2
-
2
)
N
w
/
M
w
+
P
w
holds integer reflects the j
1th and the j
2th 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 N
w
/
M
w
/
P
w constellation, by using one agile scanning beam at each satellite, the overall hop number of one traversal cycle is
(9)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
+
α
,
where
(10)
α
=
{
N
w
M
w
,
N
w
M
w
=
Odd
number
N
w
M
w
-
1
,
N
w
M
w
=
Even
number
.Proof.
The hop number for agile scanning beam analysis is divided into two different types:interplane (ISL’s between satellites in adjacent/corotating planes) andintraplane(ISL’s connecting satellites in the same orbit plane) [1].(1)
For the orbit plane numberM
w, due to the ISLs between different orbit plane in pairs, the total combination of interplane ISLs gives
(11)
C
M
w
2
=
M
w
!
2
!
(
M
w
-
2
)
!
=
M
w
(
M
w
-
1
)
2
.
The ISLs between each pair of interplane have situations of (
N
w
/
M
w
) types; then the interplane ISLs have number of
(12)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
hops.
(2)
For the intraplane ISLs, each orbital plane has satellites of(
N
w
/
M
w
), 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(13)
[
(
(
N
w
/
M
w
)
-
1
)
2
]
×
2
+
2
=
N
w
M
w
(b)
coorbit plane with even number satellites: each orbit plane get the full pairing, traversing the hop count(14)
N
w
M
w
-
1
.
Then the overall hop number of one traversal cycle is
(15)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
+
α
.Corollary 3.
According to Conclusion 1 and Figure2, for a typical Walker N
w
/
M
w
/
P
w 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
(16)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
2
.
Also, the intraplane matrix area of Figure2 will be completely filled with number of elements
(17)
N
w
2
(
N
w
M
w
-
1
)
.Proof.
(
1
) For interplane ISLs, each orbit plane has (
N
w
/
M
w
) satellites. It will be (
N
w
/
M
w
) 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
(18)
M
w
(
M
w
-
1
)
2
(
N
w
M
w
)
2
.
Considering the geometric relationships of interplane area of Figure2, the number of matrix elements of the interplane area is multiplying elements number of each square by number of squares:
(19)
(
N
w
M
w
)
2
×
M
w
(
M
w
-
1
)
2
.
Equation (18) = (19); interplane ISLs proved.
(
2
) For intraplane ISLs, each orbital plane has (
N
w
/
M
w
) satellites; due to the paired establish of ISLs between satellites in coorbit, each epoch will have
(20)
floor
(
N
w
2
M
w
)
combinations, whereinfloor is a numeric integer arithmetic.
According to the orbit planesM
w and hop number of coorbit α, the coorbit ISLs have the hop number of
(21)
floor
(
N
w
2
M
w
)
M
w
α
.
An alternate way to prove it is as follows: each orbital plane has(
N
w
/
M
w
) satellites, due to the paired establish of ISLs between satellites in coorbit, the coorbit ISLs number at each epoch have
(22)
C
N
w
/
M
w
2
=
(
N
w
/
M
w
)
!
2
!
(
(
N
w
/
M
w
)
-
2
)
!
=
(
N
w
/
M
w
)
(
(
N
w
/
M
w
)
-
1
)
2
=
N
w
2
M
w
(
N
w
M
w
-
1
)
combinations; considering the orbit plane number M
w, the coorbit ISLs have the total number of hops as
(23)
N
w
2
M
w
(
N
w
M
w
-
1
)
·
M
w
=
N
w
2
(
N
w
M
w
-
1
)
.
This also proved that the formula (21) equals formula (23).According to the geometric relationship of coorbit areas in Figure2, the number of matrix elements in coorbit areas equals the upper triangular matrix elements of each block multiplying the number of matrix
(24)
1
2
×
(
N
w
M
w
-
1
)
×
(
1
+
N
w
M
w
-
1
)
×
M
w
=
N
w
2
(
N
w
M
w
-
1
)
.Equation (23) = (24); coorbit ISLs is proved.Conclusion 2. For a given ISL topology matrix A
(
k
) at epoch k, it always can be (
1
) transformed into subunit matrixes with same column sequence using A
-
(
k
)
=
A
(
k
)
T
A
(
k
) or  (
2
) transformed into subunit matrixes with same row sequence of using A
~
(
k
)
=
A
(
k
)
A
(
k
)
T.Proof.
The ISL topology matrixA
(
k
) at any epoch k is upper triangular, according to the difference of orbit plane, the following general form of matrix block can be obtained (omitting epoch symbols k):
(25)
A
=
[
A
11
A
12
A
13
⋯
A
1
N
A
22
A
23
⋯
A
2
N
A
33
⋯
A
3
N
⋱
⋮
A
N
N
]
,
where N is the orbit plane number of constellation. It can be seen from above that A
i
i
,
i
=
1,2
,
3
…
N 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 A
i
j can be written as the following general form:
(26)
A
i
j
=
[
a
p
+
1
a
p
+
2
⋱
a
q
a
1
a
2
⋱
a
p
]
,
where q
=
M represents the number of satellites at each orbit plane. By using matrix multiplication, the A
(
k
)
T
A
(
k
) at each epoch results in
(27)
A
i
j
T
A
i
j
=
[
a
p
+
1
2
a
p
+
2
2
⋱
a
q
2
a
1
2
a
2
2
⋱
a
p
2
]
.
Since each of its elements are all ones, so(28)
a
i
2
=
1
,
i
=
1,2
,
…
,
p
,
p
+
1
,
…
,
q
.
Then the subunit matrixes with same column sequence can be obtained. For the coorbit ISL submatrixA
i
i, according to the ISL topology matrix definition, it can be written as the following general form:
(29)
A
i
j
=
[
0
a
p
+
1
0
a
1
0
⋱
⋱
a
p
0
0
a
p
+
2
0
⋱
⋱
a
q
0
]
.
The result of the multiplication between any two submatricesA
i
i
T
A
i
i at each epoch is similar as (27). Then the ISL topology matrix A
(
k
) can be transformed into subunit matrixes with same column sequence using A
(
k
)
T
A
(
k
).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
(30)
A
batch
=
∑
k
=
1
n
A
-
(
k
)
has the form of diagonal; the diagonal elements of the matrix and the eigenvalues are
(31)
diag
⁡
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]
,
λ
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]
,
where n is the hop number during one traversal cycle using one agile scanning beam.Proof.
The matrixA
-
(
k
) is the subunit matrix of A
(
k
) with the same column sequence at arbitrary epoch k, and the accumulation of the matrix A
(
k
) will fully fill the square with dimension of M
N
×
M
N in one traversal cycle of constellation. So, the accumulative matrix A
batch is substantially the accumulate of columns, with the form of diagonal matrix and diagonal elements diag
⁡
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
1
]. For the matrix in form of A
batch, the eigenvalues can be easily obtained as λ
(
A
batch
)
=
[
0,1
,
2,3
,
…
,
n
-
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 Table1.Table 1
Analysis of ISL hops and link number.
Constellations WalkerN
w
/
M
w
/
P
w
Different orbit plane
Coorbit plane
Hops/cycle
ISLs/cycle
Remark
ISL numbers
Hops (multiplication of two items)
Hopsα
ISL numbers
Pairs of ISLs/epoch
ISLs/diff-orbit*
Combinations of diff-orbit
ISL types**
ISLs/coorbitfloor
(
N
w
/
2
M
w
)
Orbit plane numberM
w
Walker9/3/2
1
3
3
3
3
1
3
12
36
Walker18/6/2
3
3
3
5
3
1
6
18
153
Walker24/3/2
1
8
8
3
7
4
3
31
276
Compass
Walker24/6/1
3
4
4
5
3
2
6
23
276
GPS
Walker25/5/1
2
5
5
5
5
2
5
30
300
Walker27/9/0
4
3
3
9
3
1
9
30
351
Walker72/6/1
3
12
12
5
11
6
6
71
2556
Walker77/7/1
3
11
11
7
11
5
7
88
2926
Iridium
Walker66/6/1
3
11
11
5
11
5
6
66
2145
Iridium
Walker120/10/1
5
12
12
9
11
6
10
119
7140
NeLS
Walker288/12/1
6
24
24
11
23
12
12
287
459504
Teledesic
*“diff-orbit” denotes the different orbit plane.
**“ISL types” denotes the ISL connection types for different orbit during one traversal cycle.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 Table1, 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 columni
_
,
j
_ 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: LetR be the real number region, with elements in lowercase letters a
,
b
,
c
…. On R; a collection of all n
×
n matrixes constitute a linear space on the basis of matrix addition and matrix multiplication, which is denoted as R
n
×
n.It can be proved that one group of bases in spaceR
n
×
n is
(33)
E
i
j
=
[
0
⋮
0
0
⋯
0
1
0
⋯
0
0
⋮
0
]
(
i
)
,
(
j
)
where i
=
1,2
,
…
,
n
;
j
=
1,2
,
…
,
n.Obviously, the Topological Method one is a spanned upper triangular matrix linear space on the base vectors(34)
E
i
j
,
(
i
,
j
satisfy
(
3
)
)
.In case of ISL failure, as shown in Figure5, the row and column i
_
,
j
_ of the matrix vanished, some special treatment is needed in ISLs reorganization.Here we assumed that the mappingf is a linear transformation of linear space R
n
×
n on real number region R; this mapping reflecting the transform from ISLs is intact to ISLs failure recombination. Based on the analysis in this section, taking a base E
i
j from linear space R
n
×
n, let α be the ISL traversal vector in R
n
×
n;
(35)
α
=
x
1
E
11
+
x
2
E
12
+
⋯
+
x
n
E
i
j
.Then we have(36)
f
(
α
)
=
y
1
E
11
+
y
2
E
12
+
⋯
+
y
n
E
i
j
.The question now is how to calculate the coordinates(
y
1
,
y
2
,
…
,
y
n
)
T of the base E
i
j under f
(
α
)? Here we write linear transformation as follows:
(37)
f
(
α
1
)
=
a
11
α
1
+
a
21
α
2
+
⋯
+
a
n
1
α
n
f
(
α
2
)
=
a
12
α
1
+
a
22
α
2
+
⋯
+
a
n
2
α
n
⋮
f
(
α
n
)
=
a
1
n
α
1
+
a
2
n
α
2
+
⋯
+
a
n
n
α
n
,
where (
a
1
j
,
a
2
j
,
…
,
a
n
j
)
T are the coordinates of f
(
α
j
) on the base E
i
j.Let(38)
A
_
=
[
a
11
a
21
⋯
a
n
1
a
12
a
22
⋯
a
n
2
⋮
⋮
⋱
⋮
a
1
n
a
2
n
⋯
a
n
n
]
.This matrix is the linear transformation on the baseE
i
j, wherein the column j elements are the coordinates of f
(
α
j
) on the base E
i
j.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 spaceR
n
×
n, assume that the vector α
∈
R
n
×
n denotes the coordinates of the matrix A in base E
i
j by using TM 1. Here we define the transform mapping from TM 1 to TM 2, and TM 1 to TM 3 are
(39)
f
12
=
α
+
α
T
,
f
13
=
α
+
α
T
+
I
.
For any vectorα
,
β
∈
R
n
×
n, we have
(40)
f
12
(
α
+
β
)
=
(
α
+
β
)
+
(
α
+
β
)
T
=
α
+
α
T
+
β
+
β
T
=
f
12
(
α
)
+
f
12
(
β
)
,
f
13
(
α
+
β
)
=
(
α
+
β
)
+
(
α
+
β
)
T
+
2
I
=
α
+
α
T
+
I
+
β
+
β
T
+
I
=
f
13
(
α
)
+
f
13
(
β
)
.
For anya
∈
R, we have
(41)
f
12
(
a
α
)
=
a
α
+
a
α
T
=
a
(
α
+
α
T
)
=
a
f
12
(
α
)
,
f
13
(
a
α
)
=
a
α
+
a
α
T
+
a
I
=
a
(
α
+
α
T
+
I
)
=
a
f
13
(
α
)
.
Clearly, mappingf
12
,
f
13 is the linear transformation in linear space R
n
×
n. Meanwhile, for a set of specific selected base E
i
j, a unique matrix can be obtained using linear transformation f
12
,
f
13 in R
n
×
n. 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 isa
i
k
b
k
j
+
c
i
j, which has two flops. This calculation requires m
n
p times which were obtained by a simple loop count, so, generally the calculation of the matrix multiplication and addition needs 2
m
n
p 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 matrixA
,
B
∈
R
n
×
n, the product C
=
A
B is also upper triangular, that is, for any k
<
i or j
<
k, there are a
i
k
b
k
j
=
0, so
(42)
c
i
j
=
∑
k
=
i
j
a
i
k
b
k
j
.Then, calculations ofc
i
j
(
i
≤
j
) require 2
(
j
-
i
+
1
) flops. The amount of the product calculation of two matrices is
(43)
∑
i
=
1
n
‍
∑
j
=
1
n
2
(
j
-
i
+
1
)
=
∑
i
=
1
n
‍
∑
j
=
1
n
-
i
+
1
2
j
≈
∑
i
=
1
n
2
(
j
-
i
+
1
)
2
2
=
∑
i
=
1
n
i
2
≈
n
3
3
.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.


## 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 Table1.Table 1
Analysis of ISL hops and link number.
Constellations WalkerN
w
/
M
w
/
P
w
Different orbit plane
Coorbit plane
Hops/cycle
ISLs/cycle
Remark
ISL numbers
Hops (multiplication of two items)
Hopsα
ISL numbers
Pairs of ISLs/epoch
ISLs/diff-orbit*
Combinations of diff-orbit
ISL types**
ISLs/coorbitfloor
(
N
w
/
2
M
w
)
Orbit plane numberM
w
Walker9/3/2
1
3
3
3
3
1
3
12
36
Walker18/6/2
3
3
3
5
3
1
6
18
153
Walker24/3/2
1
8
8
3
7
4
3
31
276
Compass
Walker24/6/1
3
4
4
5
3
2
6
23
276
GPS
Walker25/5/1
2
5
5
5
5
2
5
30
300
Walker27/9/0
4
3
3
9
3
1
9
30
351
Walker72/6/1
3
12
12
5
11
6
6
71
2556
Walker77/7/1
3
11
11
7
11
5
7
88
2926
Iridium
Walker66/6/1
3
11
11
5
11
5
6
66
2145
Iridium
Walker120/10/1
5
12
12
9
11
6
10
119
7140
NeLS
Walker288/12/1
6
24
24
11
23
12
12
287
459504
Teledesic
*“diff-orbit” denotes the different orbit plane.
**“ISL types” denotes the ISL connection types for different orbit during one traversal cycle.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 Table1, 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 columni
_
,
j
_ 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: LetR be the real number region, with elements in lowercase letters a
,
b
,
c
…. On R; a collection of all n
×
n matrixes constitute a linear space on the basis of matrix addition and matrix multiplication, which is denoted as R
n
×
n.It can be proved that one group of bases in spaceR
n
×
n is
(33)
E
i
j
=
[
0
⋮
0
0
⋯
0
1
0
⋯
0
0
⋮
0
]
(
i
)
,
(
j
)
where i
=
1,2
,
…
,
n
;
j
=
1,2
,
…
,
n.Obviously, the Topological Method one is a spanned upper triangular matrix linear space on the base vectors(34)
E
i
j
,
(
i
,
j
satisfy
(
3
)
)
.In case of ISL failure, as shown in Figure5, the row and column i
_
,
j
_ of the matrix vanished, some special treatment is needed in ISLs reorganization.Here we assumed that the mappingf is a linear transformation of linear space R
n
×
n on real number region R; this mapping reflecting the transform from ISLs is intact to ISLs failure recombination. Based on the analysis in this section, taking a base E
i
j from linear space R
n
×
n, let α be the ISL traversal vector in R
n
×
n;
(35)
α
=
x
1
E
11
+
x
2
E
12
+
⋯
+
x
n
E
i
j
.Then we have(36)
f
(
α
)
=
y
1
E
11
+
y
2
E
12
+
⋯
+
y
n
E
i
j
.The question now is how to calculate the coordinates(
y
1
,
y
2
,
…
,
y
n
)
T of the base E
i
j under f
(
α
)? Here we write linear transformation as follows:
(37)
f
(
α
1
)
=
a
11
α
1
+
a
21
α
2
+
⋯
+
a
n
1
α
n
f
(
α
2
)
=
a
12
α
1
+
a
22
α
2
+
⋯
+
a
n
2
α
n
⋮
f
(
α
n
)
=
a
1
n
α
1
+
a
2
n
α
2
+
⋯
+
a
n
n
α
n
,
where (
a
1
j
,
a
2
j
,
…
,
a
n
j
)
T are the coordinates of f
(
α
j
) on the base E
i
j.Let(38)
A
_
=
[
a
11
a
21
⋯
a
n
1
a
12
a
22
⋯
a
n
2
⋮
⋮
⋱
⋮
a
1
n
a
2
n
⋯
a
n
n
]
.This matrix is the linear transformation on the baseE
i
j, wherein the column j elements are the coordinates of f
(
α
j
) on the base E
i
j.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 spaceR
n
×
n, assume that the vector α
∈
R
n
×
n denotes the coordinates of the matrix A in base E
i
j by using TM 1. Here we define the transform mapping from TM 1 to TM 2, and TM 1 to TM 3 are
(39)
f
12
=
α
+
α
T
,
f
13
=
α
+
α
T
+
I
.
For any vectorα
,
β
∈
R
n
×
n, we have
(40)
f
12
(
α
+
β
)
=
(
α
+
β
)
+
(
α
+
β
)
T
=
α
+
α
T
+
β
+
β
T
=
f
12
(
α
)
+
f
12
(
β
)
,
f
13
(
α
+
β
)
=
(
α
+
β
)
+
(
α
+
β
)
T
+
2
I
=
α
+
α
T
+
I
+
β
+
β
T
+
I
=
f
13
(
α
)
+
f
13
(
β
)
.
For anya
∈
R, we have
(41)
f
12
(
a
α
)
=
a
α
+
a
α
T
=
a
(
α
+
α
T
)
=
a
f
12
(
α
)
,
f
13
(
a
α
)
=
a
α
+
a
α
T
+
a
I
=
a
(
α
+
α
T
+
I
)
=
a
f
13
(
α
)
.
Clearly, mappingf
12
,
f
13 is the linear transformation in linear space R
n
×
n. Meanwhile, for a set of specific selected base E
i
j, a unique matrix can be obtained using linear transformation f
12
,
f
13 in R
n
×
n. 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 isa
i
k
b
k
j
+
c
i
j, which has two flops. This calculation requires m
n
p times which were obtained by a simple loop count, so, generally the calculation of the matrix multiplication and addition needs 2
m
n
p 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 matrixA
,
B
∈
R
n
×
n, the product C
=
A
B is also upper triangular, that is, for any k
<
i or j
<
k, there are a
i
k
b
k
j
=
0, so
(42)
c
i
j
=
∑
k
=
i
j
a
i
k
b
k
j
.Then, calculations ofc
i
j
(
i
≤
j
) require 2
(
j
-
i
+
1
) flops. The amount of the product calculation of two matrices is
(43)
∑
i
=
1
n
‍
∑
j
=
1
n
2
(
j
-
i
+
1
)
=
∑
i
=
1
n
‍
∑
j
=
1
n
-
i
+
1
2
j
≈
∑
i
=
1
n
2
(
j
-
i
+
1
)
2
2
=
∑
i
=
1
n
i
2
≈
n
3
3
.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.Figure6 shows the geometry quantity associated with the ISL elevation angle. α is vertical plane perpendicular to the line O
A which is from satellite A to the center of the earth O. The elevation angle E that is from satellite A to satellite B is defined as the angle between line A
B and plane α. The center of the sphere angle θ
A
B is called ISL angle (geocentric angle).Figure 6
Geometry of ISL elevation angle.According to the geometrical relationship in Figure6 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
(44)
|
sin
(
π
2
-
E
)
|
≥
R
e
R
e
+
h
.Generally, the elevation angle is defined in interval of[
-
90
deg
⁡
,
90
deg
⁡
]; formula (44) can be solved as
(45)
E
≤
E
max
⁡
=
π
2
-
arcsin
R
e
R
e
+
h
.Since the circular orbit approximation of GNSS can be adopted, then we have|
O
A
|
=
|
O
B
|. According to the spatial geometry relationship as in Figure 6, one can obtain the following:
(46)
E
(
t
)
=
θ
A
B
(
t
)
2
.Then the ISL antenna beam coverage condition of two satellites is(47)
2
E
min
⁡
≤
θ
A
B
(
t
)
≤
2
E
max
⁡
.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 timet; here we provide an analytical method to calculate visibility according to time t.Figure7 shows the geometric relationship between satellites in different orbital plane. Orbit plane 1 and plane 2 intersect at point P in the northern hemisphere. Considering a polar coordinate with axis O
P that is from center of the earth O to polar point P, 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 φ
(
t
) of satellite at a time t is
(48)
φ
(
t
)
=
(
t
-
t
ref
)
μ
a
3
+
ω
+
η
,
where μ is the product of the gravitational constant and the mass of Earth, t
ref is reference time of ephemeris, a is orbital semimajor axis, ω is perigee at time t
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 Figure7 illustrates the ISL antenna scan area (visible area) of satellite S in orbit plane 1. Arc A
B
^ and C
D
^ of orbital plane 2 are under the visible area of the satellite S, called visual arcs hereafter, which can be represented by the corresponding geocentric angle. Bellowing provide the analysis of four angles of point A
,
B
,
C
,
D in polar coordinates in orbit plane 2.Establish Cartesian Coordinates which are centered at the Earth, as in Figure7. The fundamental plane is orbit plane 1; the vector O
x is directed from the Earth’s center to satellite P. O
z is normal to the fundamental plane, positive in the direction of the reader, and O
y completes the setup. Then the coordinates of point S are
(49)
S
=
a
[
cos
⁡
φ
S
sin
φ
S
0
]
.The coordinates of any point in orbit plane 2 are(50)
T
=
a
[
cos
⁡
φ
T
sin
φ
T
cos
⁡
⁡
γ
sin
γ
]
.The geocentric angle∠
S
O
T satisfies
(51)
cos
⁡
θ
T
S
=
cos
⁡
φ
S
cos
⁡
φ
T
+
sin
φ
S
sin
φ
T
cos
⁡
γ
which can be solved as
(52)
cos
⁡
θ
T
S
=
cos
⁡
2
φ
S
+
si
n
2
φ
S
cos
⁡
2
γ
cos
⁡
(
φ
T
-
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
)
.The mean anomaly of pointT is
(53)
φ
T
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
±
arccos
cos
⁡
θ
T
S
cos
⁡
2
φ
S
+
si
n
2
φ
S
cos
⁡
2
γ
.Accordingly, the boundary of ISL antenna beam coverage has four points(54)
φ
A
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
+
arccos
cos
⁡
2
E
min
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
B
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
+
arccos
cos
⁡
2
E
max
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
C
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
-
arccos
cos
⁡
2
E
min
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
D
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
-
arccos
cos
⁡
2
E
max
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
.Thus, when the mean anomalyφ
(
t
)
∈
[
0
,
π
], the visible arcs of orbit plane 2 from satellite S are
(55)
[
0,2
π
]
φ
(
t
)
∈
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
I
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
E
A
,
E
B
]
∪
[
2
π
-
E
C
,
2
π
-
E
D
]
φ
(
t
)
∉
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∉
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
E
A
,
2
π
-
E
D
]
φ
(
t
)
∉
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∈
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
0
,
E
B
]
∪
[
2
π
-
E
C
,
2
π
]
φ
(
t
)
∈
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∉
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
.Similarly, the visible arcs of orbit plane 2 from satelliteS can be obtained when the mean anomaly φ
(
t
)
∈
[
π
,
2
π
], which are not provided here for short.Based on the analysis above, for a typical WalkerN
w
/
M
w
/
P
w constellation, let S
i
j denote the jth satellite in orbital plane i, if the mean anomaly of S
i
j is known, we can determine any satellite S
k
l in other orbital plane according to the constellation geometry. Suppose we have mean anomaly φ
i
j
(
t
) of satellite S
i
j; then the angles φ
i
1
(
t
), φ
k
1
(
t
), φ
k
l
(
t
) of satellites S
i
1
(
t
), S
k
1
(
t
), S
k
l
(
t
) are
(56)
φ
i
1
(
t
)
=
mod
⁡
[
φ
i
j
(
t
)
-
2
π
N
w
/
M
w
(
j
-
1
)
,
2
π
]
,
φ
k
1
(
t
)
=
mod
⁡
[
φ
i
1
(
t
)
+
2
π
M
w
P
w
(
k
-
i
)
,
2
π
]
,
φ
k
l
(
t
)
=
mod
⁡
[
φ
k
1
(
t
)
+
2
π
N
w
/
M
w
(
l
-
1
)
,
2
π
]
,
where j
,
l
=
2,3
,
…
,
N
w
/
M
w.In summary, the procedures of calculating visible satellites in orbit planek from satellite S
i
j at time t are as follows:(1)
calculate mean anomaly ofS
i
j by formula (48);
(2)
calculate the visible arcs of orbit planek from satellite S
i
j by the formula (55);
(3)
calculate the mean anomalyφ
k
l
(
t
) of each satellite in orbit plane k at time t 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 satellitem is in constellation with inclination i, mean anomaly u
m, and RAAN Ω
m. The differential orbital element parameters with any other satellite n in constellation are Δ
u
m
n
,
Δ
Ω
m
n. Then we can obtain the projection of the relative velocity on the direction of LOS as
(57)
v
r
=
-
v
sin
⁡
(
2
u
m
′
+
Δ
u
m
n
′
)
(
1
-
cos
⁡
⁡
Δ
Ω
m
n
′
)
,
where
(58)
u
m
′
=
u
m
-
arctan
(
tan
Δ
Ω
m
n
2
cos
⁡
⁡
i
)
,
Δ
u
m
n
′
=
Δ
u
m
n
+
2
arctan
(
tan
Δ
Ω
m
n
2
cos
⁡
⁡
i
)
,
cos
⁡
⁡
Δ
Ω
m
n
′
=
cos
⁡
2
i
+
si
n
2
i
cos
⁡
Δ
Ω
m
n
.Therefore, the second condition for ISL establish between two satellites, is(59)
|
v
r
|
<
v
0
,
where v
0 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 Figure8.Figure 8
Flow diagram of optimal ISL routing problem.


## 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.Figure6 shows the geometry quantity associated with the ISL elevation angle. α is vertical plane perpendicular to the line O
A which is from satellite A to the center of the earth O. The elevation angle E that is from satellite A to satellite B is defined as the angle between line A
B and plane α. The center of the sphere angle θ
A
B is called ISL angle (geocentric angle).Figure 6
Geometry of ISL elevation angle.According to the geometrical relationship in Figure6 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
(44)
|
sin
(
π
2
-
E
)
|
≥
R
e
R
e
+
h
.Generally, the elevation angle is defined in interval of[
-
90
deg
⁡
,
90
deg
⁡
]; formula (44) can be solved as
(45)
E
≤
E
max
⁡
=
π
2
-
arcsin
R
e
R
e
+
h
.Since the circular orbit approximation of GNSS can be adopted, then we have|
O
A
|
=
|
O
B
|. According to the spatial geometry relationship as in Figure 6, one can obtain the following:
(46)
E
(
t
)
=
θ
A
B
(
t
)
2
.Then the ISL antenna beam coverage condition of two satellites is(47)
2
E
min
⁡
≤
θ
A
B
(
t
)
≤
2
E
max
⁡
.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 timet; here we provide an analytical method to calculate visibility according to time t.Figure7 shows the geometric relationship between satellites in different orbital plane. Orbit plane 1 and plane 2 intersect at point P in the northern hemisphere. Considering a polar coordinate with axis O
P that is from center of the earth O to polar point P, 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 φ
(
t
) of satellite at a time t is
(48)
φ
(
t
)
=
(
t
-
t
ref
)
μ
a
3
+
ω
+
η
,
where μ is the product of the gravitational constant and the mass of Earth, t
ref is reference time of ephemeris, a is orbital semimajor axis, ω is perigee at time t
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 Figure7 illustrates the ISL antenna scan area (visible area) of satellite S in orbit plane 1. Arc A
B
^ and C
D
^ of orbital plane 2 are under the visible area of the satellite S, called visual arcs hereafter, which can be represented by the corresponding geocentric angle. Bellowing provide the analysis of four angles of point A
,
B
,
C
,
D in polar coordinates in orbit plane 2.Establish Cartesian Coordinates which are centered at the Earth, as in Figure7. The fundamental plane is orbit plane 1; the vector O
x is directed from the Earth’s center to satellite P. O
z is normal to the fundamental plane, positive in the direction of the reader, and O
y completes the setup. Then the coordinates of point S are
(49)
S
=
a
[
cos
⁡
φ
S
sin
φ
S
0
]
.The coordinates of any point in orbit plane 2 are(50)
T
=
a
[
cos
⁡
φ
T
sin
φ
T
cos
⁡
⁡
γ
sin
γ
]
.The geocentric angle∠
S
O
T satisfies
(51)
cos
⁡
θ
T
S
=
cos
⁡
φ
S
cos
⁡
φ
T
+
sin
φ
S
sin
φ
T
cos
⁡
γ
which can be solved as
(52)
cos
⁡
θ
T
S
=
cos
⁡
2
φ
S
+
si
n
2
φ
S
cos
⁡
2
γ
cos
⁡
(
φ
T
-
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
)
.The mean anomaly of pointT is
(53)
φ
T
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
±
arccos
cos
⁡
θ
T
S
cos
⁡
2
φ
S
+
si
n
2
φ
S
cos
⁡
2
γ
.Accordingly, the boundary of ISL antenna beam coverage has four points(54)
φ
A
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
+
arccos
cos
⁡
2
E
min
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
B
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
+
arccos
cos
⁡
2
E
max
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
C
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
-
arccos
cos
⁡
2
E
min
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
,
φ
D
=
tan
sin
φ
S
cos
⁡
γ
cos
⁡
φ
S
-
arccos
cos
⁡
2
E
max
⁡
cos
⁡
2
φ
S
+
sin
2
φ
S
cos
⁡
2
γ
.Thus, when the mean anomalyφ
(
t
)
∈
[
0
,
π
], the visible arcs of orbit plane 2 from satellite S are
(55)
[
0,2
π
]
φ
(
t
)
∈
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
I
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
E
A
,
E
B
]
∪
[
2
π
-
E
C
,
2
π
-
E
D
]
φ
(
t
)
∉
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∉
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
E
A
,
2
π
-
E
D
]
φ
(
t
)
∉
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∈
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
,
[
0
,
E
B
]
∪
[
2
π
-
E
C
,
2
π
]
φ
(
t
)
∈
[
2
E
min
⁡
,
π
-
2
E
min
⁡
]
,
φ
(
t
)
∉
[
π
-
2
E
max
⁡
,
2
E
max
⁡
]
.Similarly, the visible arcs of orbit plane 2 from satelliteS can be obtained when the mean anomaly φ
(
t
)
∈
[
π
,
2
π
], which are not provided here for short.Based on the analysis above, for a typical WalkerN
w
/
M
w
/
P
w constellation, let S
i
j denote the jth satellite in orbital plane i, if the mean anomaly of S
i
j is known, we can determine any satellite S
k
l in other orbital plane according to the constellation geometry. Suppose we have mean anomaly φ
i
j
(
t
) of satellite S
i
j; then the angles φ
i
1
(
t
), φ
k
1
(
t
), φ
k
l
(
t
) of satellites S
i
1
(
t
), S
k
1
(
t
), S
k
l
(
t
) are
(56)
φ
i
1
(
t
)
=
mod
⁡
[
φ
i
j
(
t
)
-
2
π
N
w
/
M
w
(
j
-
1
)
,
2
π
]
,
φ
k
1
(
t
)
=
mod
⁡
[
φ
i
1
(
t
)
+
2
π
M
w
P
w
(
k
-
i
)
,
2
π
]
,
φ
k
l
(
t
)
=
mod
⁡
[
φ
k
1
(
t
)
+
2
π
N
w
/
M
w
(
l
-
1
)
,
2
π
]
,
where j
,
l
=
2,3
,
…
,
N
w
/
M
w.In summary, the procedures of calculating visible satellites in orbit planek from satellite S
i
j at time t are as follows:(1)
calculate mean anomaly ofS
i
j by formula (48);
(2)
calculate the visible arcs of orbit planek from satellite S
i
j by the formula (55);
(3)
calculate the mean anomalyφ
k
l
(
t
) of each satellite in orbit plane k at time t 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 satellitem is in constellation with inclination i, mean anomaly u
m, and RAAN Ω
m. The differential orbital element parameters with any other satellite n in constellation are Δ
u
m
n
,
Δ
Ω
m
n. Then we can obtain the projection of the relative velocity on the direction of LOS as
(57)
v
r
=
-
v
sin
⁡
(
2
u
m
′
+
Δ
u
m
n
′
)
(
1
-
cos
⁡
⁡
Δ
Ω
m
n
′
)
,
where
(58)
u
m
′
=
u
m
-
arctan
(
tan
Δ
Ω
m
n
2
cos
⁡
⁡
i
)
,
Δ
u
m
n
′
=
Δ
u
m
n
+
2
arctan
(
tan
Δ
Ω
m
n
2
cos
⁡
⁡
i
)
,
cos
⁡
⁡
Δ
Ω
m
n
′
=
cos
⁡
2
i
+
si
n
2
i
cos
⁡
Δ
Ω
m
n
.Therefore, the second condition for ISL establish between two satellites, is(59)
|
v
r
|
<
v
0
,
where v
0 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 Figure8.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 Section4, 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 beE
min
⁡
=
15
deg
⁡, E
max
⁡
=
70
deg
⁡. The maximum relative velocity is set to v
0
=
±
3
km
/
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

Figures9–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.Constellation states at initial epoch.
(a)
(b)Data uploading to visible satellite from ground.
(a)
(b)Satellites that have received the uploading data.
(a)
(b)Path searching for possible ISL establishment.
(a)
(b)Satellites that have received data by relay process.
(a)
(b)Thereafter, entering the phase of data relay, as shown in Figure12(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 Section5.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.Figures14, 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.Figure20 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 x-axis of Figure 20 denotes the number of agile scanning beams onboard, and the y-axis denotes the number of consistent ISL beams onboard. The z-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 Figure20, 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.


## 5.1. Simulation without ISL Antenna Configuration Constraint

Figures9–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.Constellation states at initial epoch.
(a)
(b)Data uploading to visible satellite from ground.
(a)
(b)Satellites that have received the uploading data.
(a)
(b)Path searching for possible ISL establishment.
(a)
(b)Satellites that have received data by relay process.
(a)
(b)Thereafter, entering the phase of data relay, as shown in Figure12(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 Section5.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.Figures14, 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.Figure20 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 x-axis of Figure 20 denotes the number of agile scanning beams onboard, and the y-axis denotes the number of consistent ISL beams onboard. The z-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 Figure20, 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.



---

*Source: 502158-2014-03-30.xml*