A System-Level Engineering Approach for Preliminary Performance Analysis and Design of Global Navigation Satellite System Constellations
A System-level Engineering Approach for Preliminary Performance Analysis and Design of Global Navigation Satellite System Constellations
Marco Nugnes , Camilla Colombo , Massimo Tipaldi Abstract – This paper presents a system-level engineering approach for the preliminary coverage performance analysis and the design of a generic Global Navigation Satellite System (GNSS) constellation. This analysis accounts for both the coverage requirements and the robustness to transient or catastrophic failures of the constellation. The European GNSS, Galileo, is used as reference case to prove the effectiveness of the proposed tool. This software suite, named GNSS Coverage Analysis Tool (G-CAT), requires as input the state vector of each satellite of the constellation and provides the performance of the GNSS constellation in terms of coverage. The tool offers an orbit propagator, an attitude propagator, an algorithm to identify the visibility region on the Earth’s surface from each satellite, and a counter function to compute how many satellites are in view from given locations on the Earth’s surface. Thanks to its low computational burden, the tool can be adopted to compute the optimal number of satellites per each orbital plane by verifying if the coverage and accuracy requirements are fulfilled under the assumption of uniform in-plane angular spacing between coplanar satellites.
Keywords : constellation design, Earth coverage, Galileo satellites, GNSS, space engineering tool Nomenclature α Right ascension [deg] G α Greenwich meridian right ascension [deg] T α Temperature gradient [K/m] β Ballistic coefficient [m /kg] γ Pitch angle [deg] δ Declination [deg] r ∆ Satellite relative distance w.r.t. the Sun [km] t ∆ Time interval [s] ε Void dielectric constant [F/m] φ Satellite longitude [deg] λ Wavelength [m] κ Geometric signal bending [m] µ ⊕ Earth gravitational parameter [km /s ] η Antenna efficiency ϑ Azimuthal angle [deg] ρ Atmospheric density [kg/m ] ψ Roll angle [deg] ν Angular spacing [deg] ω Satellite angular velocity [rad/s] p ω Pericentre anomaly [deg] ω ⊕ Earth angular velocity [rad/s] Ω Right ascension of the ascending node [deg] ζ Yaw angle [deg] a semi-major axis [km] A Rotation matrix b Leverage vector w.r.t. the pole [m] B Magnetic field vector [T]
BER
Bit Error Rate c Speed of light [m/s] a c Absorption coefficient d c Diffuse reflection coefficient s c Specular reflection coefficient D C Drag coefficient d Antenna diameter [m] zdry d Dry tropospheric delay [m] tropo d Tropospheric delay [m] zwet d Wet tropospheric delay [m] D Atmospheric drag per unit mass [N/kg] e Water vapour pressure [mbar] s e Surface water vapour pressure [mbar] b E Energy per bit
EIRP
Equivalent Isotropic Radiated Power [dB] f Signal frequency [Hz] F Solar radiation force [N] e F Solar radiation power per unit surface [W/m ] FCC
Federal Communications Commission
FSL
Free Space Losses [dB]
G-CAT
GNSS Coverage Analysis Tool x R G Receiver antenna gain [dB] x T G Transmitter antenna gain [dB]
M. Nugnes, C. Colombo, M. Tipaldi
GEO
Geostationary Orbit
GNSS
Global Navigation Satellite System
GPS
Global Positioning System H Signal bandwidth [Hz] i Orbital inclination [deg] iono i Ionospheric delay [m]
ITU
International Telecommunication Union J Second zonal coefficient L Distance between the receiver and satellite [m] gas L Gas attenuation [dB] rain L Rain attenuation [dB] x T L Transmitter line losses [dB]
LEO
Low Earth Orbit J Inertia matrix k ′ Refractivity coefficient [K/mbar] k Wet refractivity coefficient [K /mbar] B k Boltzmann constant [J/K] d k Derivative gain [Hz] p k Proportional gain m Satellite mass [kg] M Mean anomaly [deg] s n Satellite mean angular velocity [rad/s] N Signal noise power [dB] p Semi-latus rectum [km] m p Residual magnetic induction P Atmospheric pressure [mbar] s P Surface pressure [mbar] x T P Transmitter power [dB] q Elemental charge [C] r Relative position vector [km] R Inertial satellite position vector [km] R ⊕ Earth equatorial radius [km] d R Dry air specific gas constant [J/KgK] S Satellite surface [m ] S Static moment [kgm] T Disturbing torque [Nm] e T Equivalent noise temperature [K] m T Mean temperature of water vapour [K] s T Surface temperature [K]
TEC
Total Electron Content u Control action [Nm] in v Satellite inertial velocity vector [km/s] rel v Satellite relative velocity vector [km/s] V Magnetic field potential
WGS
World Geodetic System z ′ Zenith angle at the ionospheric point [deg] I. Introduction
Satellite constellations signed the beginning of a new era in terms of space mission design philosophy [1]. While in the second half of the last century monolithic spacecraft were deployed to accomplish the mission requirements, from the beginning of the new millennium the tendency changed to allocate the instrumentations and the tasks of space (scientific) missions over more than one spacecraft. The list of the future space missions involving satellite constellations demonstrates why it is useful to develop engineering approaches for the analysis of their performance. One of the most demanding missions exploiting telecommunication satellite systems has been recently proposed by SpaceX to shorten the communication time among internet users on Earth and space-faring satellites, speeding up surfing speeds. In November 2018, SpaceX received the approval for launching an additional bunch of 7518 satellites reaching the number of 11943 satellites [2] and the first 240 satellites were already injected in orbit. SpaceX is not the only company aiming to provide internet broadband services to individual users. In March 2018, OneWeb asked the approval to authorise 1260 more satellites to be added to the approved 720 ones [3], while Telesat received approval for a 117 satellite network in September 2018 [4]. Another mission involving 78-108 Low Earth Orbit (LEO) satellites is LeoSat, which is addressed to enterprises, employing big data transactions with lower latency. Indeed, this is a service that Geostationary Orbit (GEO) satellites cannot provide and the company plans to enter this new market [5]. The Federal Communications Commission (FCC) which oversees the regulations regarding telecommunication systems is adapting their rules to the incoming constellation deployment. Indeed, it gave the approval to SpaceX under the condition of full constellation deployment in nine years. If SpaceX fails to reach full deployment in that time, its authorised number of satellites will shrink to the number already in orbit [6]. This means that in the next decades several satellite constellations will be deployed. Besides communication and internet services, satellite constellations have been exploited also for Earth observation through remote sensing techniques (for instance, the study of the ionosphere total electron content parameter [7] and the global ocean altimetry [8]). This paper addresses Global Navigation Satellite System (GNSS) constellations, whose space segment consists of a group of navigation satellites in orbit around the Earth placed in an orbit positioning to get a minimum requirement for global coverage. The ground segment involves the ground stations and antennas and it is responsible for managing the constellation of navigation satellites, controlling the core functions of the navigation missions as well as determining the integrity information [9]. The user segment consists of a receiver converting the navigation signal of the satellites to position and velocity at a specific time. Another application involving GNSS is security. Indeed, the international entities for the transportation sector exploit the high precision determination offered by GNSS constellations to air
M. Nugnes, C. Colombo, M. Tipaldi traffic management to prevent the occurrence of accidents. In literature, the GNSS performance analysis has been addressed by examining specific engineering aspects. For instance, Sünderhauf et al. [10] developed a method to mitigate the multipath effect with the aim of improving the GNSS positioning determination of a generic user. Multipath effects consist in the repetition of the same signal shifted in time due to the reflection of the signal itself coming from different paths [11]. This phenomenon is more dangerous in large cities, where the reflection of the navigation signals is considerably higher. Coulot et al. [12] proposed an optimisation method of the GNSS performance by analysing the on-ground reference stations. Their work relies on the employment of algorithms to develop the future ground stations network to improve the communication with the GNSS space segment and, consequently, the performance of the system. Pan et al. [13] investigated the GNSS positioning performance in partially obstructed environments. By using a carrier double-difference model, they analysed the influence of the satellite-pair geometry on the correlation among the carrier phase observation equations, and proposed a method for determining the position of the reference satellites by using the minimum condition number rather than the maximum elevation. Tadic et al. [14] presented a web-based GNSS performance analysis and simulation tool offering a cross comparison module to benchmark devices from an established reference. Their work focuses on the modelling of the signal propagation and uses statistical data collected from existing GNSS. All the mentioned works addressed the GNSS constellation performance by considering the enhancement of only one single aspect of these complex systems, such as the navigation signal multipath phenomenon. This paper presents a system-level engineering methodology and the corresponding tool named GNSS Coverage Analysis Tool (G-CAT) to perform the design of GNSS or generic satellite constellations, by means of a coverage performance analysis starting from the initial orbital parameters of the space segment. The tool can suit different types of pointing scenarios (like geocentric, geodetic or a generic pointing altered by the effect of orbital and attitude perturbations), and gives the possibility to analyse not only the nominal configuration, but also scenarios where one or more satellites fail both for a transient period or permanently [15]. Finally, the robustness to failures associated with the accuracy in the position determination is a driver for the optimisation of the number of satellites per each orbital plane by running the tool few numbers of times. This manuscript is an extended version of [16], and explores more deeply all the aspects of the proposed methodology, together with the models embedded and their assumptions. The paper is organised as follows. In Sec. II the G-CAT functional model is described with its inputs and outputs. Sec. III presents the modelling of each function contained in the G-CAT functional model. Sec. IV describes the G-CAT graphical user interface. Sec. V presents the results obtained in case of Galileo and Global Positioning System (GPS) nominal operational scenarios and in some failure conditions. Sec. VI shows how to optimise the number of satellites per each orbital plane. Finally, Sec. VII concludes the paper and outlines the future directions.
II.
G-CAT Functional Model Overview
This section introduces the G-CAT functional model shown in the block diagram of Fig. 1. G-CAT contains a section dedicated to the space segment, a section related to the signal propagation, and one for the coverage analysis.
Fig. 1. G-CAT functional model block diagram
II.1.
Space Segment
The space segment block consists of an attitude and orbit propagator. The orbit propagator works by taking as input the initial state vector of each satellite of the constellation in terms of classic orbital elements and a propagation time expressed in hours. A standard two-body problem propagator or a perturbed propagator can be chosen. An attitude propagator is also available, which takes as input all the physical properties and the initial conditions for the Euler angles and the angular velocities of each spacecraft of the constellation. The attitude propagator will be presented in Sec. III as it is important to understand its modelling and the type of output provided to the other functions. Attitude perturbations are already embedded in the models, and the propagation time is equal to the time used for the coverage analysis.
II.2.
Signal in Space Path
The second functional block deals with the signal propagation starting from the navigation antenna up to the ground surface. All the models embedded in the tool rely on approximations of the main factors affecting the navigation signal propagation. The signal propagation requires the modelling of the interaction of an electromagnetic wave with the atmosphere, precisely two layers: the ionosphere and the troposphere.
M. Nugnes, C. Colombo, M. Tipaldi
Several phenomena like refraction, scintillation, and multipath need to be considered. For the sake of simplicity, only the refraction is considered by the tool. Both the ionosphere and the troposphere can cause the refraction of the signal for different reasons, with the result of the signal deviation from the straight path and the introduction of a delay in the reception of the signal. The ionosphere is responsible for the deflection of the signal due to the presence of charged particles which interact with the electromagnetic wave. On the other side, the troposphere also introduces a refraction and a curvature in the signal propagation, which is related to the interaction with atmosphere molecules like water vapour.
II.3.
Coverage Analysis
The last functional block determines the coverage region of each satellite given the spacecraft position vector and the direction of the navigation line of sight. This information is derived, respectively, from the orbit propagator, whose output is the spacecraft position vector in Cartesian coordinates, and from the attitude propagator. These two inputs are provided by two different propagators since it is assumed that the orbit propagation is decoupled from the satellite attitude. To determine the coverage region of each satellite, the coverage analysis functional block uses a new refined analytical approach, which models the Earth’s surface as an oblate ellipsoid of rotation [17]. Such block also contains the function in charge of the verification of the constellation coverage. Starting from the knowledge of each coverage area associated to a single spacecraft, a function verifies whether the points located on a grid computed from the Earth’s surface are within or not the coverage area and counts the number of satellites visible from the regions on the Earth’s surface. This way, it is possible to derive the coverage and the degree of accuracy according to the number of visible satellites. The coverage analysis gives as output four coverage indices: • “Red index”, i.e., the number of points on the grid covered by less than 4 satellites. • “Yellow index”, i.e., the number of points on the grid covered by exactly 4 satellites. • “Green index”, i.e., the number of points on the grid covered by more than 4 satellites. • “Global index”, i.e., the number of points on the grid covered by at least 4 satellites and it is equal to the sum of the “yellow index” and “green index”. These indices have been introduced by considering the resolution of the general global positioning problem, also known as trilateration process. Starting from the measurement of the pseudo-distance between a satellite and a receiver, a single measurement provides the position of the user on the surface of a sphere. Two measurements reduce the size of the problem into a planar one since the user’s position this time is on the region described by a circle which is the intersection of the two spheres having radii equal to the pseudo-distances between the satellite and the receivers. Finally, three measurements give as output the intersection of two circles, which allows choosing the right position between two points. Even if this ambiguity can be solved, the number of measurements is not enough due to the presence of the time-bias caused by the different accuracy between the on-board satellite’s clock and the receiver’s clock. The result is that four measurements are needed to derive unequivocally the user’s position on the Earth’s surface [18]. III.
G-CAT Functional Model Description
In this section, the modelling approach of the G-CAT functional models is presented. Each function can be configured via specific panels of the G-CAT graphical user interface, which will be presented in Sec. IV.
III.1.
Orbit Panel Functions
The “orbit panel” belongs to the functional block of the space segment. The default orbit propagator implements the standard two-body Keplerian orbit propagator, but the user has the possibility to select an orbit propagator which considers perturbations. The default propagator models in a simplified way the two most important perturbations for GNSS constellations: the non-sphericity of the Earth and the solar radiation pressure. Atmospheric drag can be neglected since the altitude of a generic GNSS is above 28000 km. For future constellations (e.g., Starlink, OneWeb), for which the altitude is under 1000 km, the effect of the atmosphere will be included. The non-sphericity of the Earth is modelled by considering only the secular effect of the J zonal harmonic on the right ascension of the ascending node, the pericentre anomaly, and the mean anomaly. The formulations for the secular rates of the orbital parameters are [19]: ( ) s n R J ip ⊕ Ω = − (1) ( ) ( ) sec sp n R J ip ω ⊕ = − (2) ( ) ( ) s n R JM ip ⊕ = − − (3) where sec Ω , sec p ω , sec M are the secular rates of the right ascension of the ascending node, the pericentre anomaly, and the mean anomaly respectively, s n represents the mean motion of the satellite, R ⊕ is the Earth mean equatorial radius, J is the second zonal harmonic coefficient, p is the semi-latus rectum and i is M. Nugnes, C. Colombo, M. Tipaldi the orbital inclination. The secular rates are used to determine the linear variation with time of the corresponding orbital parameters: t Ω = Ω + Ω ∆ (4) p p p t ω ω ω= + ∆ (5) M M M t = + ∆ (6) The solar radiation pressure is modelled by considering only the effect on the semi-major axis [19]: aa F r µ ⊕ = − ∆ (7) with a and a representing the semi-major axis and its rate respectively, µ ⊕ is the Earth gravitational constant, r ∆ is the satellite relative distance with respect to the Sun, and F is the solar radiation force. Also in this case, the semi-major axis rate is used to determine the linear variation with time of the semi-major axis. This kind of variation is not secular, but it is a cyclic perturbation with zero net effect in one solar year if the satellite is assumed to be always in sunlight: a a a t = + ∆ (8) The default perturbed orbital propagator is quite simple. Its execution does not require a significant computational load. All the satellites are assumed as concentrated in their centre of mass, therefore the orbit propagation can be decoupled from the attitude. The time scales for the orbit propagation of the satellite are extremely larger than the ones used for the attitude propagator. In Sec. IV, the advantage of the decoupling between the attitude and the orbit propagator is explained for both the short and long period analysis of the constellation. III.2.
Attitude Propagator
The attitude propagator is one of the most important functions of the tool since it allows to generalise the direction of the navigation signal and to simulate more realistic orbital scenarios. For this reason, a brief description of the modelling functions and related assumptions will be provided. In this case, each satellite is modelled as a rectangular parallelepiped rigid body. The attitude representation adopted for the constellation is the one using Euler angles and direct cosine matrices. Euler angles are used for the inputs to the propagator while all the processing is performed with the direct cosine matrices since they present the advantages to be more robust to numerical integration and do not present relevant singularities. The block diagram showed in Fig. 2 summarises the data flow within the attitude propagator. The starting point of the attitude propagator is the resolution of the dynamics through the Euler equations which needs as input the disturbances, given by the attitude perturbations, and the control actions of the actuators together with the initial angular velocity.
Fig. 2. Attitude propagator functional block diagram
The result of the dynamics is the angular velocity at the current time that is combined with the initial attitude configuration to solve the attitude kinematics in terms of direct cosine matrices. The attitude state in terms of direct cosine matrix and angular velocity is known at the current time and the corrective control action to track the desired state is implemented. This control action is added to the satellite attitude perturbations which depend only on the Cartesian position of the spacecraft and the next cycle begins solving again the Euler equations with the updated variables. From the knowledge of the initial conditions in terms of angular velocity , , x y z ω ω ω = ω , the Euler equations [20] are solved to compute the value of the angular velocity at each instant with respect to the Earth-centred inertial frame: ( ) ( ) ( ) xx x yy zz y z xyy y zz xx x z yzz z xx yy x y z J J J uJ J J uJ J J u ω ωωω ωωω ωω = − + = − + = − + (9) with xx J , yy J , zz J representing the diagonal terms of the inertia matrix J expressed in the principal axes reference frame of the satellite and x u , y u , z u are the components of the torques generated from the attitude perturbations and spacecraft control actuators. After having computed the attitude dynamics, the attitude kinematics can be retrieved. This is done by solving the following differential equation [21]: [ ] / / / B L B L B L ∧ = − A ω A (10) where / B L A is the direct cosine matrix expressing the rotation matrix from the satellite orbital frame coincident with the local frame and denoted with the letter “ L ”, and the satellite body frame denoted with the letter “ B” , while the [ ] / B L ∧ ω is the skew-symmetric matrix built from the knowledge of the angular velocity of the spacecraft orbital frame expressed in the body frame: M. Nugnes, C. Colombo, M. Tipaldi [ ] / z yB L z xy x ω ωω ωω ω ∧ − = − ω (11) The direct cosine matrix / B L A is a classic “321” rotation matrix [22]: where ψ , γ and ζ are the roll, pitch and yaw rotation angles of the body reference system with respect to the local reference frame. The angular velocity vector / B L ω is obtained as the result of two contributions: / / / / B L B N B L L N = − ω ω A ω (12) where / B N ω and / L N ω are the angular velocity vectors of the body reference frame and local reference frame with respect to the inertial reference system denoted with the letter “ N ”. One of the disadvantages in using direct cosine matrices for solving the attitude kinematics is the loss of the orthogonality property due to numerical errors after the numerical integration. This problem can be solved by using a Gram-Schmidt algorithm [23] which makes the matrix orthogonal at each instant. Another problem associated with “321” rotation matrices is the singularity occurring for values of the pitch angle equal to or . However, such values for the pitch angle are never achieved under the action of the attitude control system and so they can be neglected in the propagation processing. The last model worth mentioning is the one related to the control action, u . The formulation used is the following [24]: ( ) [ ] // / Td e p e e e e L NB N B N k k ∨ ∧ = − − − −+ ∧ u J ω J A A J ω A ωω Jω (13) with d k and p k representing the derivative and proportional gain respectively while e ω and e A are the angular velocity and attitude matrix error with respect to the desired state. The default desired state is the configuration where the spacecraft is aligned along the navigation antenna line of sight, coincident with the z -axis of the body reference frame. For this reason, the desired angular velocity should allow only rotations around the z -axis and the desired attitude matrix is the one associated with a yaw angle equal to zero. This type of control is called responsive control since it has the property to cancel the first term in the Euler equations (9) and to make the dynamics resolution faster. Of course, this method presents also the disadvantage to make slower the process if the values of the inertia matrix are not accurate in their determination [24]. After this processing, the z -direction of the body reference frame, which is assumed to be aligned with the navigation signal direction, is determined and in this way, it is possible to simulate the moving line of sight. The attitude propagator models four kinds of perturbations: atmospheric drag, gravity gradient, solar radiation pressure and magnetic torque. The atmospheric drag is the reactive force related to the atmosphere nearby the spacecraft. Even if the density is very low with respect to the Earth’s surface level, there is a small contribution. This effect becomes insignificant for altitudes higher than 10000 km. To determine the value of the torque, the expression of the drag resistance per unit mass is written: rel rel v βρ= D v with D SCm β = (15) where m is the mass of the spacecraft, D C is the drag coefficient which is in around 2-3 for common spacecraft, S is the wet area of the spacecraft, rel v is the velocity vector relative to the atmosphere, ρ is the atmospheric density and β is the ballistic coefficient, parameter that characterises the response of the spacecraft to drag. In G-CAT it is required to insert the value of the ballistic coefficient which is known from the design and geometry of the satellite, while the density is computed from the height of the satellite through an exponential model, and the relative velocity is derived from the knowledge of the inertial velocity as follows: rel in ⊕ = − ∧ v v ω R (16) The inertial velocity, in v , and the position vector in Cartesian coordinates, R , come from the orbital propagator while the [ ]
0, 0, 2 π 86400 ⊕ = ω is the angular velocity of the Earth in the geocentric equatorial frame neglecting the nutation and precession of the rotational axis. The torque associated to the atmospheric drag is, finally, given by: ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) / cos cos cos sin( ) sinsin cos cos sin sin cos cos sin sin( ) sin cos sinsin sin cos sin cos cos sin sin sin cos cos cos B L γ ζ γ ζ γζ ψ ζ γ ψ ζ ψ ζ γ ψ γ ψζ ψ ζ γ ψ ζ ψ ζ γ ψ γ ψ− = − + + + − + A (14) M. Nugnes, C. Colombo, M. Tipaldi drag = ∧
T b D (17) where b is the vector associated to the leverage between the application point of the force and the pole. The leverage is assumed to be 10% of the length of the spacecraft. The gravity gradient is the difference in terms of gravity acceleration of different parts of the spacecraft since the spacecraft is not considered as a dot point but as a 3D rigid body. To derive the expression of the gravity gradient it is necessary to express first the torque acting on an infinitesimal element of mass, dm . ( ) dmd µ ⊕ = − ∧ ++ T r r Rr R (18) where r is the relative position vector of the infinitesimal mass with respect to the body reference frame, while R is the position vector of the spacecraft centre of mass. The total torque acting on the spacecraft is the integral over the whole body of the infinitesimal torque: ( ) B dm µ ⊕ = − ∧ ++ ∫ T r r Rr R (19) It is possible to linearise the expression of the denominator since the magnitude of r is much smaller than R . Taking the series expansion of the denominator in the integral close to zero: R 1 3 R − − ⋅ + ≈ −
R rr R (20) Substituting the linearised form in Eq. (19): ( ) B dm µ ⊕ ⋅ = − ∧ − + ∫ R rT r r R (21) Considering the first addendum: ( ) B dm ∧ + = − ∧ ∫ r r R R S (22) where S is the static moment of the spacecraft. The first cross product has result equal to 0 and in the second term the position vector R is independent from the infinitesimal mass and it can be taken outside the integration. If the reference system of the body frame has the origin in the satellite’s centre of mass like the one considered, the static moment is equal to 0. The other term in Eq. (21) gives the final expression of the gravity gradient torque: ( )( ) B dm µ ⊕ = ⋅ ∧ ∫ T r R r R (23) This integral can be further developed. Indeed, considering the body reference frame, those vectors assume these expressions: ˆ ˆˆ x y z = + + r r θ k (24) ˆ ˆˆ c c c = + +
R r θ k (25) where [ ] , , c c c are the direction cosines of the radial direction. Substituting inside Eq. (23): ( ) B yc zcxc yc zc zc xc dmxc yc µ ⊕ − = + + − − ∫ T (26) Evaluating all the terms under the integral sign, and considering that the reference frame consists of principal axes of inertia: ( ) ( ) ( ) zz yyxx zzyy xx J J c cJ J c cJ J c c µ ⊕ − = − − T (27) Therefore, if one of the principal axes is aligned with the radial direction the torque is zero because only one of the direction cosines is non-zero. Solar radiation illuminating the outer surface of a satellite generates a pressure, that in turn generates a force and a torque around the centre of mass of the satellite. The main sources of electromagnetic radiation are the direct solar radiation, the solar radiation reflected by the Earth (or by any other planet) and the radiation directly emitted by the Earth. Solar radiation intensity varies with the inverse square distance from the source, so that direct solar radiation is almost constant for geocentric orbits, independent from the orbit radius since the distance to the Sun does not change significantly. On the contrary, the reflected radiation and the Earth radiation intensity are strongly dependent on the orbit radius. The average pressure due to radiation can be evaluated as: e FP c = (28) where c is the speed of light and e F is the solar radiation power per unit surface. Radiation forces can be modelled assuming that part of the incident radiation is absorbed, part is reflected in a specular way and part reflected with diffusion. Coefficients a c , d c , s c are respectively the coefficient of absorption, diffuse M. Nugnes, C. Colombo, M. Tipaldi reflection and specular reflection, and are constrained by the following equality: a d s c c c + + = (29) Merging the three components due to the three possible effects, the total force due to radiation is evaluated as: ( ) ( ) ( ) ˆ ˆ ˆˆ ˆ1 22 ˆ3 j j j s s jd j PS c cc = − ⋅ − + ⋅+
F O N O O NN (30) with ˆ O denoting the unit direction of the satellite-Sun vector and ˆ N the unit direction of the normal to each spacecraft surface j . The total torque will be the summation of the single torque on each spacecraft surface. ˆ j jj = = ∧ ∑ T r F (31) Again, the leverage of the solar radiation pressure has been considered 10% of the length of the spacecraft. The last source of possible attitude perturbations is offered by the interference of the geomagnetic field with the electronics onboard of the spacecraft. The expression for the magnetic torque is: m = ∧ T p B (32) where m p is the residual magnetic induction due to parasitic currents in the satellite and B is the Earth’s magnetic field. Normally, m p is an undesired effect while B is always somehow present. There are models that, given the satellite position, allow evaluating its components. Magnetic torques on a satellite therefore do not depend on in its inertia properties but rather on its position and attitude. So, the next open point is how the magnetic field is evaluated. The magnetic field B can be modelled in a precise way as the gradient of a scalar potential V , that is V = −∇ B . Normally, V is modelled as series expansion of spherical harmonics [21]: ( ) ( ( ) ) ( )
11 0 ( , , ) cosrsin nk n mnn mm mn n
RV r R g mh m P ϑφ φφ ϑ +⊕⊕ = = = + ∑ ∑ (33) where r B represents the radial component, positive outwards, B ϑ the co-elevation component, positive if directed towards south, and B φ the azimuth component, positive towards east. Coefficients mn g and mn h are evaluated under the assumption that polynomials mn P are normalised according to Schmidt [21]: ( ) ( ) ( ) mmn P d n π δϑ ϑ ϑ − = + ∫ (34) where m δ is the Kronecker delta, ji δ = if i j = or ji δ = if i j ≠ . In alternative, a normalisation due to Gauss [21] is possible, according to which: ,, m n mn n m P S P = (35) ( ) ( ) ( )
10 2, mn m n m nS n m n m δ − − − = + − (36) Coefficients , n m S are independent from r , ϑ and φ can be evaluated from the Gaussian coefficients mn g and mn h , following the rule: , , n m mn m n g S g = , , n m mn m n h S h = S = ( ) ,0 1,0 n n nS S n − −= for n ≥ ( ) ( )( )
11 2, , 1 mn m n m n mS S n m δ − + − + = + for m ≥ In a similar way, it is possible to get recursive formulas for the polynomials , n m P : P = ( ) , 1, 1 sin n n n n P P ϑ − − = ( ) , 1, , 2, cos n m n m n m n m P P K P ϑ − − = − , n m K = for n = ( )( )( )
12 1 2 3 n m n mK n n − −= − − for n > With these models, the components of the magnetic field, B , in the [ ] , , r ϑ φ reference frame are: ( ) ( ) ( ( ) ) ( ) nk n n mr n mn m n m RB n g mrh m P φφ ϑ +⊕= = = + + ∑ ∑ (37)
M. Nugnes, C. Colombo, M. Tipaldi ( ) ( ( ) ) ( ) cossin nk n n mn m n mn m RB g mr Ph m ϑ φϑφ ϑ +⊕= = = − ∂+ ∂ ∑ ∑ (38) ( ) ( ) ( ( ) ) ( ) nk n n mn mn m n m RB m g mrh m P φ φϑ φ ϑ +⊕= = = − − + ∑ ∑ (39) In a geocentric inertial reference frame, the components of the magnetic field assume the form: ( ) ( ) ( ) ( ) ( ) cos sin cos sin x r B B B B ϑ φ δ δ α α= + − (40) ( ) ( ) ( ) ( ) ( ) cos sin sin cos y r
B B B B ϑ φ δ δ α α= + + (41) ( ) ( ) sin cos z r
B B B ϑ δ δ= − (42) where δ is the satellite declination ( ) δ π ϑ= − and α is the right ascension, linked to the longitude of the satellite by the relation: G α φ α= + (43) with G α the Greenwich meridian right ascension. III.3.
Failure Function
This function allows turning “on” or “off” the navigation signal of a satellite for a given period. Indeed, it is possible to shut down the satellite for the whole simulation or to disable it for a transient period, and then turn it on again. This gives the possibility to simulate realistic scenarios, such as the execution of recovery procedures for the navigation payload and its navigation signal. During this time interval, the satellite is not operating properly, and it should be considered as not existing in the constellation. Another possible failure scenario is the shutdown of the navigation payload due to under-voltage or over-voltage problems. In such cases, the navigation payload operational status switches from the nominal mode into a safe mode, where it stays until the problem has been solved. The failure function gives the possibility to consider many cases where the outage of the navigation signal occurs and to check how the performance of the GNSS constellation is affected.
III.4.
Coverage Function
The coverage function is in charge of determining the coverage area where a satellite is visible (see Fig. 3), by starting from the satellite position vector and the navigation signal direction and modelling the Earth’s surface as a sphere or an oblate ellipsoid of rotation according to the World Geodetic System (WGS-84) representation [25]. The reader can refer to [17] for more details on the equations used to compute the coverage area for an oblate ellipsoid of rotation. Such equations simplify to the ones of a perfectly spherical Earth if the values of the equatorial and polar radius of the Earth are the same.
Fig. 3. Geometry for the field of view in the Earth-centred reference frame
III.5.
Signal Propagation Function
The tool offers the possibility to have a first preliminary estimation of the link budget starting from the signal properties and the antenna’s characteristics. Moreover, it gives as output a figure of merit on the delay in the signal propagation due to the interaction with the atmosphere. It is important to highlight that these computations are decoupled from the coverage analysis since they are external and add simply more information. The link budget is one of the most important quantities during the design phase of a space mission giving information about the quality of the telecommunication signal. For navigation satellites where the signal represents the payload of the mission, it becomes of crucial importance. The link budget results as output of a series of operations the signal to noise ratio, starting from a fixed set of variables like carrier frequency, antenna gain, spacecraft height, and internal hardware. Once the signal to noise ratio is evaluated, according to the type of codification and modulation used for the signal, the quality of the communication is verified. The starting point of the link budget computation is the carrier frequency, car f , of the navigation signal which must be approved from the International Telecommunication Union (ITU). Then, according to the type of antenna onboard and the hardware it is possible to evaluate the Equivalent Isotropic Radiated Power (EIRP) summarising the quality of the signal exiting from the M. Nugnes, C. Colombo, M. Tipaldi transmitter. The EIRP is evaluated using the following formula [26]: x x x
T T T
EIRP P G L = + − (44) where x T P is the power of the transmitter, x T G is the antenna gain and x T L represents the line losses associated to the transmitter. It has to be highlighted that in the previous formula all the quantities are expressed in decibel ( dB ). The value of the antenna gain x T G is fixed once the physical properties of the antenna have been established: x T dG π ηλ = (45) with d representing the diameter of the antenna, λ the wavelength of the signal and η the antenna efficiency. This quantity is non-dimensional and must be converted in dB in order to be used:
10 log x xdB
T T
G G = (46) Once the EIRP is evaluated, the next step is to compute the losses existing for the distance between the user receiver on the Earth’s surface and the antenna of the satellite. This kind of losses are called Free Space Losses (FSL) and may be estimated using the following formulation [26]: FSL L λπ = (47) where L represents the distance between the user receiver and the antenna. Also, in this case the conversion in dB is necessary. At this stage the signal arrives at the boundary of the atmosphere and so new losses must be introduced. For this purpose, some statistical plots modelling the value of the loss associated to an atmospheric phenomenon exist. The atmospheric losses may be distinguished in two main components: gas losses and rain losses. The gaseous losses are associated to the resonance of the frequency of the signal with the frequency of vibration of the molecules in the atmosphere, in particular molecular oxygen, O , and water vapour, H O . It is important to select a value for the carrier far away from the absorption peaks as it is shown in Fig. 4.
Frequency (GHz) -5 -4 -3 -2 -1 S pe c i f i c a tt enua t i on ( d B / k m ) OxygenWater vapourTotal
Fig. 4. Specific attenuation losses due to atmospheric gases versus the navigation signal frequency
The second source is rain and for each location on the Earth a certain probability of rain in one year is evaluated. This way, it is possible to have a figure of merit for the expected loss in one year as shown in Fig. 5. Both the models for the specific attenuation due to gases and rain have been taken from [27]. -3 -2 -1 Percentage of time (%) R a i n A tt enu t a t i on ( d B ) Fig. 5. Specific attenuation losses due to rain versus the percentage of rain time for different navigation signal frequencies
The last part in the signal to noise ratio computation is to evaluate the properties of the receiving antenna. Also, in this case there is a gain that can be computed in the same way done for the transmitter. However, an additional source of disturbance is present which is represented by the noise power which depends on the bandwidth and on the physical temperature of the system. Generally, an equivalent noise temperature e T for each component of the receiver is introduced to express noise power. The temperatures are summed to get an equivalent system noise temperature in such way that the noise power is expressed as: B e
N k HT = (48) M. Nugnes, C. Colombo, M. Tipaldi with B k the Boltzmann constant equal to − × and H the bandwidth of the signal. The final quantity that is evaluated is the energy per bit to noise power ratio, b E N , which is strictly related to the signal to noise ratio [26]: x b rain gas R E EIRP FSL L L G NN = − − − + − (49) At this stage according to the type of modulation, the Bit Error Rate (BER) is derived as function of the energy per bit to noise power ratio following [28]. The BER is an index of the quality of the signal as it is shown in Fig. 6. Indeed, in order to have a good communication it is expected that the BER is at least
10 bit s − . -5 0 5 10 15 20 25 E b /N (dB) -6 -5 -4 -3 -2 -1 E s t i m a t ed BE R BPSKQPSK8PSK16PSK32PSK4QAM16QAM64QAM256QAM
Fig. 6. Estimated Bit Error Rate (BER) versus signal over noise ratio for different typed of codification
The modelling of the ionosphere is performed using first-order approximations taken from [29]. The ionosphere is the upper part of the atmosphere where charged particles are mixed with neutral particles. The charged particles are created by photoionisation caused by incoming UV and X-radiation from the Sun: gas molecules are heated, and electrons are liberated then. The rate of this ionisation depends on the density of gas molecules and the intensity of the radiation. In the neutral atmosphere charged particles are practically absent, since the created charged particles are recombined rapidly due to the high density of particles. In the ionosphere, however, only the charged particles can influence the propagation of radio waves. Mainly, the free electrons affect the propagation since the free ions are much heavier than the electrons. The interaction of the ionosphere with the navigation signal can be decomposed into four sources of delay called first-order, second-order, third-order delay, and geometrical signal bending denoted respectively as (1) i , (2) i , (3) i and κ [29]. The second and third-order delays are often referred to as the ionospheric higher-order terms. The total ionospheric delay can be written as: (1) (2) (3) iono i i i i κ= + + + (50) All the formulations for the evaluation of the different sources for the ionospheric delay are taken from [29]. The first-order delay, (1) i , which is also the most significant for the ionospheric delay computation is approximated as follows: (1) 2 Ai TECf = with
22 0 e qA m π ε= (51) where f is the frequency of the navigation signal, TEC is the Total Electron Content representing the integral of the electron density in the ionosphere, q − = × is the elemental charge, e m − = × is the electron mass and ε − = × is the void dielectric constant. The second-order delay, (2) i , is associated to geomagnetic field: ( ) (2) 3 cos2 e qAi TECf m ϑπ≈ B (52) with ( ) cos ϑ B representing the projection of the magnetic field vector along the signal path. The third-order delay, (3) i , can be formulated as follows: e Ai N TECf τ≈ (53) where ,max e N is the maximum electron density in the ionosphere and τ represents a shape factor that is equal to 0.66. Finally, the geometric signal bending is evaluated by using the following approximation: ( ) tan8 e A z N TECf κ η′≈ (54) with z ′ representing the zenith angle at the ionospheric point, that is the intersection between the line of sight connecting the satellite with the user receiver and the ionosphere boundary modelled as a thin-uniform shell located at 350 km from the Earth’s surface. The tropospheric delay represents the time shift introduced by the neutral part of the atmosphere due to the introduction of a geometric curvature of the navigation signal because of the variation in the refraction index in the different layers of the troposphere. The modelling of the tropospheric delay is carried out M. Nugnes, C. Colombo, M. Tipaldi also in this case by decomposing it into two main terms in order to distinguish the effects due to the dry gases (primarily nitrogen and oxygen) and the water vapour. The formulation of the dry tropospheric delay is given by [30]: z dhyd sm k Rd Pg − = (55) where k = is the refractivity constant of the dry gases, d R = is the specific gas constant for dry air, m g is the gravity at the centre of mass of the column of air and s P = mbar is the surface pressure. The second part of the tropospheric delay which is denoted as wet delay can be evaluated with the following equation [30]: ( ) ( )
10 1 m dzwet m k k T Rd eg − ′ += Λ + (56) ( ) T dm s m
RT T g α = − Λ + ( ) s s Pe e P λ+ = with m T the mean temperature of water vapor, ( ) w d k k k M M ′ = − = a modified refractivity coefficient, k = × the wet refractivity coefficient, Λ an empiric coefficient depending on the local conditions, T α the temperature gradient with respect to the surface height, s e = mbar the surface water vapour pressure, s T the surface temperature, e the water vapour pressure and P the tropospheric pressure that can be computed using the exponential model of the standard atmosphere. The values of the coefficients used in G-CAT are T α − = × and Λ =
The total tropospheric delay is get from the summation of the two contributions: z z ztropo hyd wet d d d = + (57) The expressions mentioned before are valid for the zenith directions. If the direction is different from the zenith, the effect in the new direction is obtained by considering mapping functions that require further modelling and they are not implemented in G-CAT.
IV.
G-CAT Graphical User Interface and Configuration
In the following section, the G-CAT user graphical interface (see Fig. 7) is presented. The procedures to configure the G-CAT tool, run a simulation, and analyse the related results are also provided. The G-CAT tool configuration consists in defining all the parameters needed to set up a specific simulation scenario. The first step is to insert the initial state vector in terms of orbital elements for each satellite of the constellation. The tool can handle up to 30 satellites at the current time (in line with the number of satellites of existing GNSS), and a tab is dedicated to each satellite. After the selection of the type of orbit propagator, the propagation time can be inserted. This time can be different from the coverage analysis time. This way, the orbit and attitude propagations are performed as part of the coverage performance analysis. More specifically, the
Fig. 7. G-CAT graphical user interface
M. Nugnes, C. Colombo, M. Tipaldi two propagations are decoupled and carried out for different integration times. Indeed, the orbit perturbations time scales are much larger than the attitude ones. This means that the user can start an orbit propagation for a longer time and the tool gives as output the evolution of the orbit parameters with each time instant. At this stage, the user can change the initial orbital parameters selected with the one obtained from the orbit propagation. This way, it is possible to simulate the attitude behaviour for a small time interval by starting from the orbital parameters of the constellation at a different epoch to check the effects of the orbital perturbations on the performance. Once the orbit propagation is concluded, the attitude panel can be filled in the same way of the orbital one. For each satellite a tab is available where the geometrical characteristics can be inserted, together with the initial attitude conditions: • Principal moment of inertia • Length, width and height, assuming that each satellite is a parallelepiped • Ballistic coefficient • Reflectance and absorptance coefficients • Gains for the control law • Initial Euler angles and angular velocities. There is also the possibility to select the type of attitude perturbations to be included in the propagation. The failure panel is a box where the user can indicate whether a satellite is in nominal or in safe mode, and in the latter case, the duration of such off-nominal condition as explained in Sec. III. The signal in space panel is another secondary box divided in “payload properties”, where the information about the navigation signal should be filled to compute the figures of merit for the signal to noise ratio, and the “atmosphere interaction”, where the tropospheric and ionospheric delays are computed. As already mentioned in Sec. III, the evaluations of the two delays are decoupled from the coverage analysis because the G-CAT tool calculates a first-order estimation for those variables. Finally, in the coverage panel, the user can choose the type of modelling for the Earth’s surface (sphere or oblate ellipsoid of rotation), the input for the navigation signal, the type of pointing, the set of points to be analysed, and the type of output. The navigation signal is modelled as a conical field of view by having as axis of symmetry the navigation signal direction. For this reason, the coverage area is derived by considering either the conical half-aperture angle or the minimum elevation angle from where the satellites are visible on the Earth’s surface [17]. A small panel offers the possibility to perform both the local and the global coverage analysis. In the first case, the user can insert the geographic coordinates of the location whose coverage has to be verified. In the latter case, the tool discretises the Earth’s surface as a grid, and verifies if every point on such a grid is covered. The user has also the possibility to choose three different types of pointing, according to the mission requirements: geocentric, geodetic and perturbed. The geocentric pointing has the satellite always looking towards the centre of the Earth. The geodetic pointing has the satellite pointing along the local normal direction, which is different from the geocentric direction in the oblate ellipsoid case. Finally, the perturbed pointing starts from an initial direction, and then considers the evolution over time of the navigation signal direction caused by the attitude perturbations. Only in the case of perturbed pointing selection the attitude propagator is involved in the coverage analysis together with the orbit propagator. Indeed, the other two pointing modes imply a fixed direction for the navigation antenna and no attitude propagation is required. The user can also select if the output of the analysis should be plotted on a 3D representation of the Earth’s globe or on a planar one. The G-CAT tool flags with a red cross all the regions on the Earth’s surface which are not covered for at least one time instant during the whole simulation. The last parameter to be set before the tool can start the simulation is the time for the coverage performance. Simulations are not time-consuming processes and can last from less than 5 minutes (in case of geocentric or geodetic pointing) up to 10 minutes (in case of perturbed pointing) using an Intel® Core™ i7-4710HQ CPU @2.50 GHz. Fig. 8 shows how the G-CAT functions are executed when the coverage performance analysis is launched. The tool receives as input all the initial data regarding the
Fig. 8. Order of execution of the G-CAT functions for the GNSS coverage performance analysis
M. Nugnes, C. Colombo, M. Tipaldi orbit and attitude configuration of each satellite of the constellation. The user can decide to analyse the full operative constellation or to select some transient or permanent failures. The first function executed by the tool is the orbital propagator which computes the position of all the spacecraft for the desired time. The attitude propagator (if needed) combines this information with the initial attitude configuration to determine the navigation antenna line of sight, that is the direction aligned along the signal transmitted from the payload onboard each satellite. The last function performs the coverage analysis and uses as input the position vector of each satellite and the navigation antenna line of sight, computed by the two propagators. Such function determines first the coverage area associated to the constellation, counts the number of times a point on the Earth’s surface is inside the coverage region of each satellite, and according to the final sum, the four indices defined in Sec. III are computed. V. GNSS Coverage Analysis and Results with the G-CAT Tool
This section presents the results obtained by the G-CAT tool when configured with the nominal operational conditions of the GPS and Galileo constellations. These results are used to validate the tool design and implementation. In the second part of the section some failure scenarios are considered to show all the functionalities of the tool.
V.1.
GPS Nominal Configuration
The GPS is the first GNSS constellation in orbit for navigation and communication purposes. The nominal configuration is a Walker Delta pattern 30/6/1, which involves 30 satellites placed in 6 equally spaced orbital planes each one having an inclination of 55° with respect to the equatorial plane. The minimum number for obtaining GPS information is 24 [31], but the additional satellites provide global coverage with a refined precision positioning. They can also be used as spare satellite. The revolution period of each satellite is half sidereal day. All this information is used to configure the G-CAT tool together with a geocentric pointing. The coverage performance analysis can be run for a single satellite revolution. The inputs for the coverage panel part are the following: • Ellipsoidal representation of the Earth’s surface. • Minimum elevation angle equal to 5°. • Geocentric pointing. The results of the tool for the green and global coverage indices for a full operative capability of 30 and 24 satellites are shown in Fig. 9 and Fig. 10, respectively. In both the configurations, the coverage is global. However, as for the case of 24 satellites, it is not guaranteed that each point on the Earth’s surface is visible by more than 4 satellites. Moreover, the additional satellites can improve the failure robustness of the GPS constellation.
Fig. 9. Green and global coverage indices for GPS considering 30 operative satellites
Fig. 10. Green and global coverage indices for GPS considering 24 operative satellites
V.2.
Galileo Nominal Configuration
The European GNSS Galileo is a 24/3/1 Walker Delta pattern constellation. For this reason, there are 24 satellites placed in 3 equally spaced orbital planes each one having an inclination of 56° [32]. The nominal configuration includes 6 additional satellites to be used as spare satellites in case of catastrophic failures. The revolution period of each satellite is slightly higher (15 hours) and a global coverage is required. It is possible to repeat the same procedure performed for the GPS constellation by configuring similar options for the simulation: • Ellipsoidal representation of the Earth’s surface. • Minimum elevation angle equal to 5°. • Geocentric pointing. As shown in in Fig. 11, the coverage is global and there is also a higher percentage of regions on the Earth’s
M. Nugnes, C. Colombo, M. Tipaldi surface visible from more than 4 satellites for an amount of time longer than the GPS case.
Fig. 11. Galileo coverage indices in nominal configuration
V.3.
Failure Case
It is interesting to analyse what happens when there are off-nominal configurations. The first failure scenario involves a catastrophic failure of the first satellite of the constellation. In order to set a catastrophic failure, the option “payload safe mode” has to be flagged to simulate the transition from the nominal mode into a configuration where there is the outage of the navigation signal (see Sec. III). The Galileo case study with the following G-CAT configuration is considered: • Ellipsoidal representation of the Earth’s surface. • Half-aperture angle equal to 12°. • Geodetic pointing. • Catastrophic failure of the first satellite of the constellation.
Fig. 12. 2D representation of critical regions for failure case
The outcome of the simulation is presented both in 2D version (see Fig. 12) and 3D version (see Fig. 14). The red crosses identify the points on the Earth’s surface uncovered for at least one time instant during the whole simulation. This is the criterion for checking the global coverage for the GNSS constellation.
Fig. 13. Galileo coverage indices for failure case
Even if the coverage indices of Fig. 13 already show that a global coverage is not achieved, the 2D and the 3D representations of the terrestrial globe allow a direct association between the critical points and their geographic locations on the Earth’s surface.
Fig. 14. 3D representation of critical regions for failure case
V.4.
Failure Case
The G-CAT tool allows analysing the GNSS coverage performance with both a fixed direction and a moving line of sight. In this operational scenario, the effect of the attitude perturbations in the analysis for the Galileo constellation has been included, thus the following inputs have been set in the G-CAT user interface: • Ellipsoidal representation of the Earth’s surface.
M. Nugnes, C. Colombo, M. Tipaldi • Half-aperture angle equal to 12°. • Perturbed pointing. • Catastrophic failure of the first satellite of the constellation.
The results are presented in Fig. 15, where the yellow circles highlight the difference with respect to Fig. 12. Indeed, the perturbed pointing introduces differences in the coverage analysis that cannot be clearly understood by looking only at the coverage index plots because of their similarity (see Fig. 13 and Fig. 16). Even though from an average point of view the two types of pointing seem to show similar coverage performance, it is important to use the perturbed pointing and the 2D/3D representations to identify which are the actual affected regions by the attitude perturbations.
Fig. 15. 2D representation of critical regions for failure case
V.5.
Failure Case
The last operational case addresses the occurrence of a transient failure. For example, there is the possibility that a reconfiguration of the navigation signal makes the satellite unavailable for a limited time. The coverage performance changes for this small portion of time and even if nothing varies in terms of average percentage index, the actual coverage is affected. In this case, in the failure panel of the tool, the “payload safe mode” has to be flagged, but this time the “recovery time” shall be less than the time used for the coverage analysis simulation. There is also the possibility to shift the time instant at which the failure begins by filling the “failure time” box. The same set of inputs has been considered with the only difference that the first satellite of the constellation is unavailable for 20 minutes: • Ellipsoidal representation of the Earth’s surface. • Half-aperture angle equal to 12°. • Perturbed pointing. • First satellite of the constellation not operative for 20 minutes.
Fig. 16. Galileo coverage indices for failure case
The outcome of the simulation is shown in Fig. 17 and Fig. 18. As expected, the global coverage is not achieved during the failure recovery time, which is much shorter than the orbital period. This is shown by the red cross of Fig. 18.
Fig. 17. Galileo coverage indices for failure case
The failure recovery procedure is defined by the ground stations, and it is important to execute it properly to minimise its effect on the coverage performance requirements. For instance, uncrowded regions can be addressed by such recovery procedure, instead of densely populated areas.
M. Nugnes, C. Colombo, M. Tipaldi
Fig. 18. 2D representation of critical regions for failure case
VI.
GNSS Design using G-CAT
In this section, the G-CAT tool is used to design the optimal number of satellites per orbital plane given the orbit parameters as well as the coverage and accuracy requirements. Again, Galileo has been used as reference for discussing the results provided by the tool. Basically, the design process consists in the execution of many simulations considering each time a different number of satellites. This pragmatic approach can be applied under the assumption of having the same number of satellites in all the planes with uniform angular spacing between coplanar satellites. Thanks to the low computational time associated to each simulation, it is possible to store the results of all the simulations and plot them to choose the best GNSS constellation configuration. This way, the number of satellites per orbital plane can be optimised.
The first parameter to be analysed is the global coverage index to determine the minimum number of satellites in order to achieve a global coverage. This is shown in Fig. 19, where it is clear that this minimum number is equal to 6 satellites per orbital plane associated to an angular spacing of 60°. At this stage, the obvious question to be asked is why the actual configuration of Galileo constellation involves 8 satellites per orbital plane when there are other cheaper configurations giving the same result in terms of global coverage. The answer has to be investigated in the robustness of the satellite constellation to failures. For this reason, three different failure scenarios can be considered: • Case 1: complete loss of the first satellite of the first orbital plane. • Case 2: complete loss of the first satellite of the first orbital plane and of the second satellite of the second orbital plane. • Case 3: a mix of transient failures affecting the first and second satellites of the first orbital plane, the third satellite of the second orbital pane and the second satellite of the third orbital plane.
40 45 50 55 60 65 70 75 80 85 90020406080100120
Fig. 19. Mean global coverage index for different Galileo constellation configurations
The two parameters under investigation for the optimisation process are the global index and the green index for the reasons explained hereafter.
VI.1.
Global Index Optimisation
The first parameter to be optimized is the global index. By looking at Fig. 20, it is evident that the configuration with 6 satellites per plane cannot fulfil the coverage requirements. The new optimal configuration is given by an angular spacing equal to ν = : it corresponds to 7 satellites per orbital plane. This is a proper configuration since, in all the three failure scenarios, the global index keeps close to 100%, which is the expected result to have a global coverage.
40 45 50 55 60 65 70 75 80 85 90020406080100120
Fig. 20. Mean global coverage index for different Galileo constellation configurations
VI.2.
Green Index Optimisation
In this case, only the green index, which is the percentage of number of points on the grid representing the Earth’s surface covered by more than 4 satellites, has to be optimised. The reason for this choice is related to
M. Nugnes, C. Colombo, M. Tipaldi the GNSS accuracy requirements which are usually very tight. In Sec. II it has been explained that the minimum number of satellites to calculate the user’s position is equal to four. Of course, there is no possibility to improve the results in the position determination if only four measurements are provided for the trilateration process since they are associated to a unique solution. This is the reason why the green index has to be optimised in such a way that there are more satellites measurements. This gives the possibility of improving the accuracy of the determined position. The outcome of the green index optimisation process is shown in Fig. 21, which proves that a proper configuration can be achieved with 8 satellites per orbital plane. Indeed, the previous optimal number of 7 satellites per orbital plane does not guarantee the accuracy required from the Galileo constellation. An important remark is that Galileo original configuration involved 9 satellites per orbital plane, while the actual configuration relies on 8 operative satellites per plane. The G-CAT optimisation process outcome is in line with the actual Galileo constellation configuration, and this confirms the potentialities of using the G-CAT tool for the GNSS constellation preliminary design.
40 45 50 55 60 65 70 75 80 85 90020406080100120
Fig. 21. Mean green coverage index for different Galileo constellation configurations
VII.
Conclusion
In this paper, a system-level engineering approach for the coverage performance analysis of a generic GNSS constellation has been presented. The corresponding tool, named GNSS Coverage Analysis Tool (G-CAT), has been designed and implemented. Such tool provides different types of pointing scenarios and gives the possibility to analyse not only the nominal configurations of the GNSS constellation, but also to verify its robustness to transient or catastrophic failures. Each function, such as the orbit propagator, can be configured via specific panels of the G-CAT graphical user interface. Space system engineers can select the type of modelling for the Earth’s surface, the pointing mode, and can add the presence of catastrophic or transient failures before performing the GNSS coverage analysis. Moreover, such tool can also be used to determine the number of satellites per orbital plane given the corresponding orbital parameters and coverage requirements. The results have been validated by considering existing GNSS constellations, such as GPS and Galileo. Considering its current implementation status, the G-CAT tool can only be used during the initial phases of the GNSS constellation design. As for future work, the G-CAT models will be refined in order to provide space system engineers with a tool suitable for the design of future satellite constellations during more advanced project stages. In this regard, it is important to increase the total number of satellites handled by the tool, since future LEO constellations are expected to include hundreds of satellites. Another important aspect is to improve the modelling of both the ionospheric and tropospheric effects on the navigation signal propagation, such that the resulting coverage performance analysis becomes more realistic. A further refinement for both the orbit and attitude propagators is planned. The orbit propagator can provide more accurate state vectors for each satellite of the constellation, while the current ideal control law of the attitude propagator can be enhanced by modelling the real control torques provided by the spacecraft actuators (e.g., reaction wheels). More complex failure scenarios can also be implemented in order to include more difficult recovery procedures at constellation level and assess their impact on the coverage performance. Finally, it is planned to extend the G-CAT capabilities for other relevant GNSS applications, such as the ones related to the GNSS Safety-of-Life services.
Acknowledgements
References [1]
G. A. Orndorff and B. F. Zink, A Constellation Architecture for National Security Space Systems,
John Hopkins APL Technical Digest 29 (2010), 273-282. [2] [3]
C. Henry, OneWeb asks FCC to authorize 1200 more satellites (March 2018). Available: http://spacenews.com/oneweb-asks-fcc-to-authorize-1200-more-satellites/. [4]
C. Henry, Telesat says ideal LEO Constellation is 292 Satellites, but could be 512 (September 2012). Available: http://spacenews.c
M. Nugnes, C. Colombo, M. Tipaldi om/telesat-says-ideal-leo-constellation-is-292-satellites-but-could-be-512/. [5]
C. Henry, Jsat hopes its LeoSat stake will ward off cannibalization (June 2017). Available: http://spacenews.com/jsat-hopes-its-leosat-stake-will-ward-off-cannibalization/. [6]
Space Exploration Holdings, Application for Approval for Orbital Deployment and Operating Authority for the SpaceX NGSO Satellite System, Memorandum Opinion, Order and Authorization,
FCC 18-38 (2018) . [7] E. Music et al, The Total Electron Content from InSAR and GNSS: a midlatitude study,
IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 11 (2018), 1725-1733 . [8] J. Mashburn et al, Global Ocean Altimetry with GNSS Reflections from TechDemoSat-1,
IEEE Transactions on Geoscience and Remote Sensing 56 (2018), 4088-4097 . [9] V. Tyapkin et al, Characteristics of the Ground-Based User Equipment of Satellite Radio Navigation System Employing Pseudolites in Noisy Environments,
International Review of Aerospace Engineering (2018), 11 (2), 58-65. [10]
N. Sünderhauf et al, Multipath mitigation in GNSS-based localization using robust optimization,
IEEE Intelligent Vehicles Symposium (2012), 784-789 . [11] I. Kurniawati et al., Spatial Correlation and Diversity of Low-Latitude Multipath Multimode HF Skywave Radio Channels,
International Journal on Communications Antenna and Propagation (2019), 9 (2), 110-116. [12]
D. Coulot et al, Global Optimization of GNSS station reference networks,
GPS Solutions Journal 19 (2015), 569-577 . [13] S. Pan et al, New Approach for Optimising GNSS Positioning Performance in harsh observation environments,
Journal of Navigation 67 (2014), 1029-1048 . [14] S. Tadic et al, Design of GNSS Performance Analysis and Simulation Tools in the GAPFILLER Web Portal, st Telecommunication Forum (2013), 236-239 . [15] European Cooperation for Space Standardization, Space Product Assurance – Failure modes, effect (and criticality) analysis (FMEA/FMECA),
ECSS-Q-ST-30-02C (2009) . [16] M. Nugnes, C. Colombo, and M. Tipaldi, A System Engineering Tool for the Optimisation of a GNSS Constellation Design,
Proc. 5 th IEEE International workshop on Metrology for Aerospace (2018), 300-304 . [17] M. Nugnes, C. Colombo, and M. Tipaldi, Coverage Area Determination for Conical Fields of view considering an Oblate Earth,
Journal of Guidance, Control, and Dynamics 29 (2019) . [18] B. Soewito et al., Bluetooth Low Energy: Comparing 4 Trilateration Models in Indoor Positioning System,
International Journal on Communications Antenna and Propagation (2018), 8 (6), 500-509 . [19] D. Vallado,
Fundamentals of Astrodynamics and Applications , fourth ed. (Microcosm, Hawthorne, 2013). [20]
T. M. Habib, In-Orbit Spacecraft Inertia, Attitude, and Orbit Estimation Based on Measurements of Magnetometer, Gyro, Star Sensor and GPS Through Extended Kalman Filter,
International Review of Aerospace Engineering (2018), 11 (6), 247-251. [21]
J. R. Wertz,
Spacecraft Attitude Determination and Control , first ed. (Kluwer Academic Publishers, Dordrecht, 1978). [22]
Y. Meddahi A Backstepping Approach for Airship Autonomous Robust Control,
International Review of Aerospace Engineering (2013), 6 (6), 284-289 . [23] A. Quarteroni et al,
Matematica Numerica , fourth ed. (Springer, Milan, 2000). [24]
J. D. Biggs, A quaternion-based Attitude tracking controller for Robotic Systems,
Institute for Mathematics and Applications: Mathematics for Robotics (2015) . [25] V. Janssen, Understanding Coordinate Reference Systems, Datums and Transformations,
International Journal of Geoinformatics (2009) . [26] J. R. Wertz and W. J. Larson,
Space Mission Analysis and Design, third ed. (Microcosm, Hawthorne, 1999). [27]
K. Bennett,
MATLAB Applications for the Practical Engineer , first ed. (InTech, London, 2014). [28]
J. G. Proakis and M. Salehi,
Digital Communications, fifth ed. (McGraw-Hill, New York, 2008). [29]
D. Odijk, Fast Precise GPS Positioning in the Presence of Ionospheric Delays, Ph.D. dissertation, TU Delft, Delft, 2002. [30]
J. Haase, Accuracy and Variability of GPS Tropospheric Delay Measurements of Water Vapor in the Western Mediterranean,
Journal of Applied Meteorology (2003), 1547-1568. [31]
Department of Defense, Global Positioning System Standard Positioning Service Performance Standard, fourth ed., 2008. [32][32]