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

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


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][2][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 2 International Journal of Navigation and Observation 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.

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.
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  1 th satellite in orbital plane  1 established a ISL with the  2 th satellite in orbital plane  2 ; then the element of the matrix is 1; the others are 0. According to (1), ISL has another representation of 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 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.
Thus, when the  1 th satellite in orbit plane  1 established an ISL link with the  2 th satellite in orbit plane  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.
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).

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, 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 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  1 th satellite in the second orbit plane relative to the reference satellite is Also, the phase difference of the  2 th satellite in the orbit plane   relative to the reference satellite is Obviously, when holds integer reflects the  1 th and the  2 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 where 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 The ISLs between each pair of interplane have situations of (  /  ) types; then the interplane ISLs have number of 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 (b) coorbit plane with even number satellites: each orbit plane get the full pairing, traversing the hop count Then the overall hop number of one traversal cycle is 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 Also, the intraplane matrix area of Figure 2 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 combinations; considering the orbit plane number   , the coorbit ISLs have the total number of hops as 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 Equation (23) = (24); coorbit ISLs is proved.
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 ): 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: where  =  represents the number of satellites at each orbit plane.By using matrix multiplication, the ()  () at each epoch results in Since each of its elements are all ones, so 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: International Journal of Navigation and Observation 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 has the form of diagonal; the diagonal elements of the matrix and the eigenvalues are 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].

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

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.First of all, here we give a brief introduction of some properties of the topology matrix and linear space: Let R be the real number region, with elements in lowercase letters , ,  . ... On R; a collection of all  ×  matrixes constitute a linear space on the basis of matrix addition and matrix multiplication, which is denoted as R × .
Obviously, the Topological Method one is a spanned upper triangular matrix linear space on the base vectors E  , (,  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 R × on real number region R; this   International Journal of Navigation and Observation mapping reflecting the transform from ISLs is intact to ISLs failure recombination.Based on the analysis in this section, taking a base E  from linear space R × , let  be the ISL traversal vector in R × ; Then we have The question now is how to calculate the coordinates ( 1 ,  2 , . . .,   )  of the base E  under ()?Here we write linear transformation as follows: where ( 1 ,  2 , . . .,   )  are the coordinates of (  ) on the base E  . Let This matrix is the linear transformation on the base E  , wherein the column  elements are the coordinates of (  ) on the base E  .
The analysis of ISLs redundancy situation and restructure is similar as above, which is not provided here.
Clearly, mapping  12 ,  13 is the linear transformation in linear space R × .Meanwhile, for a set of specific selected base E  , a unique matrix can be obtained using linear transformation  12 ,  13 in R × .Therefore, we create a oneto-one relationship between the matrix and the linear transformation; that is, TM 1/2/3 are equal.[9].Commonly used computation efficient evaluation algorithm on modern computers is flop (floating octal point); one flop is a floatingpoint 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.

Analysis of Computational Efficiency
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 ,  ∈ R × , the product  =  is also upper triangular, that is, for any  <  or  < , there are     = 0, so Then, calculations of   ( ≤ ) require 2( −  + 1) flops.The amount of the product calculation of two matrices is 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.

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][11][12][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 endto-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 constraint models in the optimal ISL routing problem are introduced in detail hereafter, and the overall optimization process will finally be provided.

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).
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 Generally, the elevation angle is defined in interval of [−90 deg, 90 deg]; formula (44) can be solved as Since the circular orbit approximation of GNSS can be adopted, then we have || = ||.According to the spatial International Journal of Navigation and Observation geometry relationship as in Figure 6, one can obtain the following: Then the ISL antenna beam coverage condition of two satellites is 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 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.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 The coordinates of any point in orbit plane 2 are The Thus, when the mean anomaly () ∈ [0, ], the visible arcs of orbit plane 2 from satellite  are Similarly, the visible arcs of orbit plane 2 from satellite S 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 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.

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 where Therefore, the second condition for ISL establish between two satellites, is where V 0 is the relative velocity range which depends on the performance of satellite communication payload.

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.

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.8km, 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 = 15 deg,  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.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.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.

Simulation without ISL Antenna Configuration Constraint.
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.

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.

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 realtime constellation status and dynamic calculations.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.
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 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.

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.

Figure 5 :
Figure 5: Demonstration of one ISL failure.
"diff-orbit" denotes the different orbit plane.* * "ISL types" denotes the ISL connection types for different orbit during one traversal cycle.

Figure 7 :
Figure 7: The ISL antenna beam coverage analysis.

Figure 10 :
Figure 10: Data uploading to visible satellite from ground.

Figure 11 :
Figure 11: Satellites that have received the uploading data.

Figure 12 :
Figure 12: Path searching for possible ISL establishment.

Figure 13 :Figure 14 :Figure 15 :
Figure 13: Satellites that have received data by relay process.

Table 1 :
Analysis of ISL hops and link number.