Analysis and Optimization of Earth Observation Micro-Constellations

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Grado Ingeniería Aeroespacial Analysis and Optimization of Earth Observation Micro-Constellations Autor: Francisco Pallarés Chamorro Tutor: Rafael Vázquez Valenzuela Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Trabajo de Fin de Grado Ingeniería Aeroespacial Analysis and Optimization of Earth Observation Micro-Constellations Autor: Francisco Pallarés Chamorro Tutor: Rafael Vázquez Valenzuela Profesor Titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Trabajo de Fin de Grado: Analysis and Optimization of Earth Observation Micro-Constellations Autor: Francisco Pallarés Chamorro Tutor: Rafael Vázquez Valenzuela El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Contents 1 Introduction 1 1.1 Scope of the project 2 1.2 Structure of the document 3 2 Preliminary considerations 5 2.1 Set up 5 2.1.1 Geocentric equatorial frame 5 2.1.2 Perifocal frame 6 2.1.3 Geographic equatorial frame 7 2.1.4 Vehicle Velocity, Local Horizontal frame 8 2.1.5 Spacecraft frame 9 2.1.6 Normalised frame 10 2.2 Orbital Mechanics 11 2.2.1 Two-body problem 11 2.2.2 Orbit perturbations 14 2.2.3 Orbit propagation 17 2.3 Objective Regions 18 2.4 Satellites 20 2.5 Ground stations 21 2.6 Observations 22 2.6.1 Swaths 23 2.6.2 Modes 33 2.7 Satellite-Antenna transmissions 33 3 Computational geometry modelling 35 3.1 Geometry 35 3.2 Observation windows 37 3.3 Transmission windows 40 3.4 Windows incompatibility 41 3.5 Subregions 42 4 Problem optimization 45 4.1 Mathematical description 45 4.1.1 Mathematical entities 45 4.1.2 Mathematical relations 48 4.1.3 Variables 48 4.1.4 Constraints 49 4.1.5 Objective function 50 4.2 Mathematical resolution 50 I II Contents 4.2.1 Observations scheduling sub-problem 51 4.2.2 Observations and transmissions scheduling problem 52 5 Results 55 5.1 Resources 55 5.1.1 Satellites 55 5.1.2 Ground stations 56 5.2 Scenarios 56 5.2.1 SPAIN 57 5.2.2 SPAIN AND ITALY 64 5.3 Final analysis 69 6 Conclusions and future work 71 List of Figures 73 List of Tables 75 Bibliography 77 1 Introduction Remote sensing techniques provide the capacity of collecting information of spots from afar. This tool has an important application on Earth’s monitoring from space, a practice which help to gather multiple data from any point on Earth’s surface in short periods of time: vegetation biomass, water quality, surface profile, surface temperature, human infrastructures spread... This knowledge is used for a diversity of purposes of human interest, such as military, agricultural, meteorological or ecological ends. In order to acquire this information by remote sensing practice an entire infrastructure has been set, divided into space and ground segments. The space segment comprises the satellite or satellites, also called Earth Observation Satellites (EOS), which are equipped with sensors for the reception of data, usually in the form of electromagnetic signals. These satellites establish transmissions with the ground segment by means of antennas placed in certain spots, in order to sent the information collected. Then the ground segment is responsible for the processing of the data acquired for its use on a particular application. Earth observation satellites have been used for decades now. Historically, large satellites were required in order to do the observation tasks. These satellites were manually coordinated and their time resolution, linked with the time needed to revisit and acquire data from the exact location [ 11 ], was very low. Currently, the miniaturization tendencies have also arrived to the space sector and smaller satellites are now developmed. The emergence of the CubeSats, included in the nanosatellites group, and the growth of the number of their launches make clear the future trend: Earth observation with small and numerous satellites. Figure 1.1 Amount and forecast of CubeSats’ launches per year. The increasing amount of EOS will allow a coordinated solution for Earth’s surface monitoring with a higher time resolution. However, this concept inevitably requires a more efficient and automatised scheduling of the observations and transmission times for each of the satellites of the constellation involved in the task. In addition, the raise of time resolution together with the specific demands of clients calls for a real-time analysis, in which the computing time for the resolution plays an important roll. 1 8Chapter 2. Preliminary considerations The relationship between this frame and the geocentric equatorial frame is given from the Greenwich mean sidereal time at the instant of time (GMST), which states the angle between the Vernal equinox ( X⊕ -axis) and the Greenwich meridian (Xgph-axis). θg=GMST ·15 o hr (2.11) cos(θg) =~ i⊕·~ igph (2.12) Hence, taking into account the right ascension ( α ) and declination ( δ ) angles definition on the Geocentric frame, the transformation between the frames is immediate: λ=α+θg(2.13) φ=δ(2.14) This is the key frame for the representation and computation of geographical zones and allows to relate the model of visibility used with the topography on Earth’s surface. It is also appropriate to indicate that for this project, a simplification of ideal spheric Earth has been used. 2.1.4 Vehicle Velocity, Local Horizontal frame In order to analyse the attitude of the satellite a local reference frame is used. On the Vehicle Velocity, Local Horizontal frame (VV LH ), ZVV LH -axis is along the negative position vector (nadir): ~ kVV LH k −~ Rsat (2.15) while YVVLH -axis is defined along the negative orbit normal: ~ jVV LH k −~ Rsat ×~ Vsat (2.16) Hence, the XVV LH -axis is towards velocity, but not necessary parallel to it. Figure 2.4 VVLH reference frame. The rotation from VVLH frame to Geographical frame can be described by means of the following equivalences: ~ iVV LH cgph =~ jVV LH cgph ×~ kVV LH cgph (2.17) ~ jVV LH cgph =−~ Rsat cgph ×~ Vsat cgph |~ Rsat cgph ×~ Vsat cgph|(2.18) ~ kVV LH cgph =−~ Rsat |~ Rsat |,(2.19) 2.1 Set up 9 therefore, the rotation matrix between the two frames is: TVV LH cgph = ~ iVV LH cgph ~ jVV LH cgph ~ kVV LH cgph  (2.20) It is important to note that VVLH frame is a local frame while the Geographic frame is an Earth-centred; thus, not only a rotation, but also a translation defined by the vector ~ Rearth , which represents Earth’s center position from the satellite, is needed for changing between these two frames. ~ RearthcVV LH =Rsat   0 0 1 (2.21) 2.1.5 Spacecraft frame The spacecraft reference frame ( sc ) is bounded to the satellite and rotates with it. At an initial instant before the manoeuvring, when the sensor( Zsc -axis) points towards the nadir, the Spacecraft frame matches the VVLH frame. The manoeuvre achieves the desired aim direction by rotations along the Xsc , Ysc and Zsc -axes, directing Zsc-axes towards the objective: ~ ksc =~uaim (2.22) Once the definitions between these two last frame is made , it is important to identify the attitude angles of a satellite. •Roll angle (ω): Angle defined along the positive Xsc-axis direction. •Pitch angle (ϕ): Angle defined along the positive Ysc-axis direction. •Yaw angle (χ): Angle defined along the positive Zsc-axis direction. In order to aim at a certain objective on Earth’s surface, only two manoeuvres along their axes are needed. In our case, first rolling ( ω ) followed by pitching ( ϕ ) manoeuvres have been chosen. These rotations can be described by their rotation matrices: Troll =  1 0 0 0 cos(ω)sin(ω) 0−sin(ω)cos(ω) ,(2.23) Tpitch =  cos(ϕ)0−sin(ϕ) 0 1 0 sin(ϕ)0 cos(ϕ) .(2.24) Therefore, the matrix of aiming, which corresponds to the rotation matrix from VVLH frame to Spacecraft frame, is: Taim =TpitchTroll =  cos(ϕ)sin(ω)sin(ϕ)−cos(ω)sin(ϕ) 0 cos(ω)sin(ω) sin(ϕ)−sin(ω)cos(ϕ)cos(ω)cos(ϕ) .(2.25) The rotation matrix from the Spacecraft frame to the VVLH frame, as well as the expression of the Zsc -axis in the VVLH frame are as follows: 10 Chapter 2. Preliminary considerations TsccVV LH =T0 aim =  cos(ϕ)0 sin(ϕ) sin(ω)sin(ϕ)cos(ω)−sin(ω)cos(ϕ) −cos(ω)sin(ϕ)sin(ω)cos(ω)cos(ϕ) ,(2.26) ~ ksccVV LH =~uaimcVV LH =  sin(ϕ) −sin(ω)cos(ϕ) cos(ω)cos(ϕ) (2.27) 2.1.6 Normalised frame The normalised frame ( norm ) represents a Non-inertial Earth-centred frame and is used for the computation of the limits of the visibility swath or visible region. It is defined by the position of the satellite ( ~ Rsat ), which belongs to Xnorm -axis; and the aim direction, contained in the XZnorm -plane and pointing to the positive values of the Znorm -axis. For the computation mentioned, the radius of Earth ( R⊕ ) is used to nondimensionalise the distances; therefore, the points obtained from the algorithm are multiplied by this magnitude to get the real position on this frame. The Normalised frame is described by means of its relation with the local frames, VVLH frame and the Spacecraft frame, as follows: ~ inorm k − ~ kVV LH ,(2.28) ~ jnorm k~ inorm ×~uaim ,(2.29) ~ knorm k~ inorm ×~ jnorm.(2.30) It is important to define the off-nadir angle ( Γ ) which measures the angle between the aim direction and the nadir, strongly related with the quality of the images taken by the satellite. For larger angles of aim the satellite capture Earth’s surface farther from the normal, worsening the quality of the taken image. cos(Γ) =~uaim ·− ~ inorm=~uaim ·~ kVV LH =cos(ω)cos(ϕ).(2.31) The rotation from Normalised frame to VVLH frame can be expressed by means of the previous equivalences: ~ inormcVV LH =  0 0 −1 ,(2.32) ~ jnormcVV LH =−~ inormcVV LH ×~uaimcVV LH | ~ inormcVV LH ×~uaimcVV LH |(2.33) ~ knormcVV LH =~ inormcVV LH ×~ jnormcVV LH .(2.34) therefore, the rotation matrix is: TnormcVV LH = ~ inormcVV LH ~ jnormcVV LH ~ knormcVV LH  (2.35) As in the transformation between VVLH and Geographic frames, the Normalised frame is Earth-centred and a translation given by the vector ~ Rsat is needed for the change of frame. ~ Rsat cnorm =Rsat   1 0 0 (2.36) 2.2 Orbital Mechanics 11 2.2 Orbital Mechanics On Earth observation missions, the satellite’s position determination along time is an essential task. Once the attitude and specifications of the sensor are fixed, the relative position between this one and Earth’s surface defines which area is within its visibility limits. For this reason an accurate position determination is required in this project. With this in mind, an introduction to orbital mechanics is presented in this section, in order to determine the orbits and positions of the satellites considered, taking also into account the most representative perturbation involved in the problem. 2.2.1 Two-body problem The position of the satellites on this problem is ruled by Newton’s laws of motion, by which is stated the relation between the forces and the accelerations governing the movement of any body. In particular, Newton’s second law establishes the relation between the acceleration of a body, ~ai , described on an inertial frame, and the sum of the forces, ~ Fi j, acting on it by introducing the gravitational mass of this body, mi, ∑ ∀j ~ Fi j =mi~ai. By the integrations of the body’s acceleration, its velocity and position are determined along time. Hence, the motion of the satellite is uniquely defined by their masses and the forces acting on them once an initial state has been set. Newton’s gravitational law The gravitational force is the main factor governing the orbital motion of bodies. This force, F , describes the attraction between two particles with masses, m1 and m2 , as proportional to each of their masses and inversely proportional to the square of the distance, r, between them, F=Gm1m2 r2 where Gis established as the gravitational constant, G=6.6742 ×10−11 m3 Kgs2 Figure 2.5 Gravitational forces [1]. Therefore, including the vectorial analysis of the attraction between the two point of masses described by their positions, ~ R1and ~ R2, the force, ~ F12 acting on 1due to 2attraction, ~ F12 =Gm1m2 |~ R2−~ R1|3~ R2−~ R1, or defining the relative distance, r=|~ R2−~ R1|, and the versor pointing 2from 1,~u12 =~ R2−~ R1 |~ R2−~ R1|, ~ F12 =Gm1m2 r2~u12. 12 Chapter 2. Preliminary considerations The force acting in the other body is defined by Newton’s third law, ~ F21 =−~ F12. Equations of motion Given an inertial reference frame, and the position vectors of both point of masses relative to this frame, ~ R1 and ~ R2, the motion of each one is defined by Newton’s second law, m1 ¨ ~ R1=Gm1m2 r2~u12,(2.37) m2 ¨ ~ R2=−Gm1m2 r2~u12.(2.38) The potential energy, V, of the gravitational force, ~ F, between these two masses is established by, V=−Gm1m2 r,(2.39) from which the gravitational force can be derived thanks to its conservative character, ~ F=−∇V.(2.40) It is important to describe thoroughly the relative motion of one point, 2, around the other, 1, due to the non-inertial character of this approximation. Consequently, the relative acceleration, ¨ ~r=¨ ~ R2−¨ ~ R1,(2.41) ¨ ~r=−G(m1+m2) r2~ur,(2.42) and defining the gravitational parameter, µ, µ=G(m1+m2),(2.43) the equation of relative motion can be rewritten as, ¨ ~r=−µ r3~r.(2.44) The Eq. (2.44) consists on a non-linear second-order differential equation which governs the motion of the point of mass 2, relative to 1. The double integration of this equation together with its six constraints of integration, such as initial relative velocity and position, lead to the description of the relative position as function of time. This equations describe the ideal situation with two point masses, while the bodies involved in this project are the satellites and the planet Earth. However, by considering the Earth as a perfect sphere, its gravitational field created outside its interior region is identical to the created one by a point mass in the center of this sphere with the Earth’s mass, M⊕ . In addition, due to the difference of masses, M⊕msat , and sizes, the satellite can be as well considered as a point mass and the previous description for the relative motion can be used to describe the position in space between these two entities. The gravitational parameter, once again due to the difference between masses, is practically function of the greater mass only, and it is defined for the planet Earth, µ≈GM⊕=3.986004 ×105Km3 s2. Position in the orbit Once the differential equation which describes the relative motion of the satellite around the Earth is established, by introducing the angular momentum of the orbit, h , which is a constant, together with other 2.2 Orbital Mechanics 13 necessary of parameters, the relative position of the orbiting body can be defined directly on its orbital plane, r=h2 µ 1 1+ecosθ.(2.45) This one represent a parametric equation of a conic curve, where e is the eccentricity of the orbit, related to the type of conic curve, and θ is the true anomaly, the angle between the position vector of the orbiting body and the position vector of the periapsis point, i.e. the point of minimum distance between the satellite and the central body. Figure 2.6 Orbital plane description [1]. The artificial satellites orbits around Earth are elliptical, hence the eccentricity, e<1, and in particular, with regard to Earth observation satellites, their orbits are practically circular, e≃0. This description defines the position of the orbiting body in the perifocal reference frame ( p f ). From the coordinates on this frame, the position in the geocentric equatorial ( ⊕ ) frame can be established by a rotation transformation (Subsec. 2.1.2). Position as function of time In order to describe the variation of the position in time two additional angles are established: •Eccentric anomaly: The definition of the eccentric anomaly, E , is shown in Fig. 2.7 using the auxiliary circle which cirumscribes the ellipse. Figure 2.7 Eccentric anomaly definition [1]. 14 Chapter 2. Preliminary considerations The eccentric anomaly, E, is related to the true anomaly, θ, through the next equation: tan θ 2=r1+e 1−etan E 2.(2.46) •Mean anomaly: In this case, the mean anomaly, M , is the angle between the vector of a fictitious body moving around the ellipse at a constant angular speed, n=qµ a3 , and the periapsis vector. This angular speed defines an orbital period, T, equal to the one involving the real object. M=nt =2π Tt.(2.47) The mean anomaly, M, is also described by the eccentric anomaly by the relation, M=E−esinE,(2.48) also known as the Kepler’s equation. Therefore, the true anomaly, θ , of the object and the elapsed time, t , are related using the previous expressions. 2.2.2 Orbit perturbations The previous model for the determination of the position of an orbiting object lies in the assumption of a situation where the gravitational field created by a perfect homogeneous sphere of matter is the only agent with which the body interacts. However, in the real situation other factors participate in satellite’s motion, perturbing the ideal one. This factors, also called perturbations, need to be analysed in order to get a precise position of the orbiting bodies. A model used in order to address the problem of perturbation might be rewriting Eq. (2.44) adding a term, ~p, representing the acceleration due to the perturbations, ¨ ~r=−µ r3~r+~p.(2.49) The effects of some of the perturbation which affect a generic satellite orbiting Earth at a height of 1000 km can be estimated by stating the orders of magnitude of the perturbing acceleration: a0=µ r2 •Earth’s oblateness: p∼10−2a0 •Lunar gravity: p∼10−7a0 •Solar gravity: p∼10−7a0 •Solar radiation pressure: p∼10−9a0 •Atmospheric drag: p∼10−10a0 In case of the atmospheric drag, its perturbation becomes rapidly negligible with altitudes. Nevertheless, for altitudes of the Earth observation satellites orbits, its presence manifests and it should be considered. The other predominant perturbation considered in this analysis is the effect of the gravitational field created by a non-spherical Earth. Due to the heterogeneity of Earth, its gravitational field does not perfectly matches the one created by a perfect homogeneous sphere and an analysis of its potential should be carried out. The effects of the perturbations mentioned are divided into three groups: •Secular effects: represent linear variations in the orbital elements. •Long-period effects: represent periodic variations in the orbital elements whose periods is greater than the orbital one. •Short-peeriod effects: represent periodic variations in the orbital elements whose periods is less than the orbital one. 2.2 Orbital Mechanics 15 Atmospheric Drag The atmospheric drag is caused by the interaction between the spacecraft and the particles of the atmosphere. The dependence of atmospheric air’s density on the altitude has been an important matter of study for the aerospace sector, and as a result several models has been developed, most of them by interpolating values of density of the air at certain altitudes. In Fig. 2.8 is presented the US Standard Atmosphere 1976. Figure 2.8 US Standard Atmosphere 1976 [1]. As shown, the density at the altitudes where orbiting bodies exist is low, but their effects still can be noticed eventally deorbiting the satellites. The drag is modelled as a force which oppose to the relative velocity of an object moving on a fluid, in this case the atmosphere, therefore the relative velocity between the satellite and an atmosphere which rotates with Earth must be established, ~vrel =~v−~vatm =~v−~ Ωe×~r, with ~ Ωethe angular velocity vector of Earth, ~ Ωec⊕=  0 0 7.2921159 ×10−5  rad s. Then the drag can be expressed as the vector, ~ D=−D~urel, defining the versor, ~urel,with the direction of the relative velocity,~vrel, vrel =|~vrel | ~urel =~vrel vrel . The magnitude of the drag, D , can be estimated thanks to the drag coefficient of the body, CD , the density at the altitude of the orbit, ρand the frontal area of the body A, D=1 2ρv2 relCDA, 16 Chapter 2. Preliminary considerations from which the perturbing acceleration can be completely defined, ~p=−1 2ρvrel CDA m~vrel. Usually, the model uses the ballistic coefficient, B, in stead of the drag coefficient, B=CDA m, or the analogous BSTAR coefficient, B∗, which employs a reference value of the air density, ρ0, B∗=ρ0B 2=ρ0CDA 2m. Gravitational perturbations The gravitational potential per unit of mass given a perfect and homogeneous sphere is given by Eq. (2.39) , V=−µ r, and the acceleration is given by its gradient, ~a=−∇V=−µ~r r3. Nevertheless, the Earth’s form is not a perfect sphere, but more similar to an oblate spheroid. This fact has consequences on the gravitational field which is generated. In order to model the gravitational perturbation, a term representing the perturbation of the gravitational potential, Φ(r,φ), is used, V=−µ r+Φ(r,φ). This term only considers the zonal harmonics, this is the differences between a perfect sphere and Earth along the latitudes; for this reason it is called the rotationally symmetric perturbation, and it is only function of the distance, r, and the latitude, φ. The rotationally symmetric perturbation potential can be expressed as infinite series: Φ(r,φ) = µ r ∞ ∑ k=2 JkR rk Pk(cosφ) where Jk are the zonal harmonics, R is the equatorial radius and Pk are the Legrende polynomials. Some of the first zonal harmonics are, J2=0.00108263 J3=−2.33936 ×10−3J2 J4=−1.49601 ×10−3J2 ... while the Legrende polynomials can be obtained from Rodrigues’ formula, Pk(x) = 1 2kk! d dxkx2−1k Finally, the acceleration due to gravitational perturbations can be derived from the gradient of the rotationally symmetrical perturbation potential, ~p=−∇Φ 2.2 Orbital Mechanics 17 Sun-synchronous orbits: An important secular effect of the gravitational perturbations strongly related with Earth observation satellites is the nodal regression, whereby the right ascension of the ascending node of the orbit, Ω , varies along time due to the effect of the first zonal harmonic, J2. This variation can be expressed as, ˙ Ω=−3 2 nJ2R2 a2(1−e2)2cosi where ais the semimajor axis of the elliptical orbit, and nrepresents the mean motion of the orbit, a=h2 µ 1 1−e2, n=rµ a3. For a given combination of semimajor axis, a , and eccentricity, e , the inclination of the orbit, i , fixes the variation of the right ascension, ˙ Ω . This fact can be used to achieve a desired rate of change. In case of Earth observation, the variation matches the variation of the sun position during the Earth translation movement, hence the illumination of the spot captured by the sensor at a given latitude remains approximately equal during the year. For the considered LEO satellites on this project: n≃15.26 rev day , a≃6866km, e≃0, i≃97.4o, establishing a value for the nodal regression of ˙ Ω∼0.99 o day , which matches the sun position variation around the geocentric frame of reference. 2.2.3 Orbit propagation At this point, the necessary definitions for the determination of the position of the satellites at every instant of time have been set out. However, regarding orbital perturbations, in stead of numerically integrate the equation of motion (2.49) , it is usually preferred to describe the orbital elements at a given instant of time, and analyse their evolution by using a propagation. Here is explained the procedure to establish an initial state of the satellite and propagate it along time. Two Line Element set The state of a satellite at a certain instant of time, or epoch, is described in its Two Line Element (TLE) set. This data format gather the orbital elements of the object into two lines of 69 characters which are used for the propagation (osculating parameters). Figure 2.9 Two Line Element set format [12]. 24 Chapter 2. Preliminary considerations •Cone of vision: The vision of the satellite’s sensor will be modelled as the points from the space inside a cone. This cone will be defined by its vertex, Psat , its axis vector, ~uaim , and its aperture angle, FOV (Field of vision). In order to define in a better way the surface points of this cone, an auxiliary vector perpendicular to the axis, ~w , will be used. The sum of these two vectors is equal to a third vector, ~v=~uaim +~w, which defines the generatrix lines of the cone’s surface through it vertex, Psat . Figure 2.15 Cone of Vision. The vector ~w depends on the parameter θ , defined as the angle between the vector and the plane y=0 . For each value of θthe points define a surface line of the cone which contains the vertex. ~rcone(θ) = Psat +σ[~uaim +~w(θ)].(2.58) Hence, the vector is defined as follows, ~uaim ·~w=0,(2.59) |~w|=|~uaim|tanFOV 2,(2.60) ~w·~ j=|~w|sin(θ).(2.61) Defining ~uaim as a versor (unity vector), |~uaim|=1, therefore, |~w|=tanFOV 2,(2.62) ~w·~ j=tanFOV 2sin(θ);(2.63) finally, taking into account Eq. (2.57) ~w=  −tanFOV 2sin(Γ)cos(θ) tanFOV 2sin(θ) −tanFOV 2cos(Γ)cos(θ) .(2.64) •Scan sector: In case the observation region is shifted creating a swath, the vision entity which describes the gathered area is simplified from the cone of vision to a scan sector. This sector sweeps the surface under the satellite while is moving acquiring the desired information. The scan sector is defined by 2.6 Observations 25 its vertex or center, which corresponds to the satellite’s position ( Psat ); its angle, equal to the angle of vision of the sensor ( FOV ); and the plane in which it lies, containing the aim line. The orientation of this plane will be set relying on the definition of vector ~w(2.64) , which is established as perpendicular to the satellite’s ground track at every time step. Once, known the velocity of the satellite on the simplified situation frame ( Psat ∈X ), it is useful to determine the versor which defines the trajectory parallel to Earth’s surface, to wit, the satellite’s ground track. ~ N=~ Vsat −~ Vsat ·~ i~ i,(2.65) ~n=~ N |~ N|;(2.66) therefore, the expression for vector ~nis: ~n=  0 ny nz  ;(2.67) The computation of the edges of the sector can be deduced from the algorithm of computation of the generatrix lines of the cone of vision (Eq. (2.58) to (2.64) ). Nevertheless, as mentioned before, the vector ~wsatisfies the equation: ~w·~n=0.(2.68) Taking into account the expression of ~w described on Eq. (2.64) , the values of θ for the sector edges fulfill: tan(θ) = nz ny cos(Γ);(2.69) which give a pair of values for θ: θ1=arctan nz ny cos(Γ)!,(2.70) θ2=arctan nz ny cos(Γ)!+180o.(2.71) Figure 2.16 Scan Sector. 26 Chapter 2. Preliminary considerations Instantaneous visible region For the first situation, in which a single shot of Earth’s surface is taken on a certain moment from the satellite’s position, the computation of every point of intersection between the cone of vision, which models the visible space region for the satellite’s sensor, and the sphere is needed to define the visible surface region at this instant. This intersection between the cone and the sphere is computed for each line as follows: For a given θ, P?=Psat +σ?[~uaim +~w(θ)],(2.72) σ?:|OP?|=1.(2.73) Covering all the values for θ,0<θ<2π, all the intersection points are calculated. In case the cone intersects the sphere only partially, the rest of the points needed to close the visibility region are computed considering tangent lines to the sphere through the point Psat . The points selected are the points of tangency between these lines and the sphere which are inside the cone of vision, and represent the limit of visibility not because of the sensor’s aperture (FOV ), but due to Earth’s curvature. Figure 2.17 Limit of possible visibility. These points fulfill the equations deducted from the tangency condition: OP?⊥Psat P?,(2.74) |OP?|=1,(2.75) which implies, x?=1 A,(2.76) y?2+z?2=1−1 A2.(2.77) In order to visualise the results, three examples are shown (Fig. 2.18a, and 2.18b). First two of them have the same position of the satellite ( A=1.3 ) and aperture ( FOV =10o ), however in the second (Fig. 2.18b) a bigger off-nadir angle ( Γ=45.3o ) is considered, obtaining a much more elongated shape of the visible region. On the third example (Fig. 2.18c), a greater off-nadir angle ( Γ=50o ) produces a partial intersection and the limits for the visible region are computed as discussed before. 2.6 Observations 27 (a) Intersection Cone-Sphere, Example 1 . (b) Intersection Cone-Sphere, Example 2 . (c) Intersection Cone-Sphere, Example 3 . Figure 2.18 Visibility regions. Shifted visible region In case the satellite carries out a "scanning" of Earth’s surface while moving, the resolution is similar to the one already explained. However, instead of modelling the vision of the satellite as the region inside a cone of visibility, it is defined as the surface scanned by a sector whose vertex corresponds to the position of the satellite. The intersection between the edges of the scan sector and Earth’s surface at each time step produces both limits of the swath of visibility while the satellite moves. The rest of the algorithm used to compute the intersection points between these two lines and the sphere is identical to the one used previously to compute the intersection of the cone (Eq. (2.72) to (2.77) ), with 28 Chapter 2. Preliminary considerations the exception that now the angle θ only takes values of θ1 and θ2 , generating the pair of points pursued P? 1=P?(θ1)and P? 2=P?(θ2). FOV =30o,Γ=0o Figure 2.19 Visibility swath and regions . As shown on Fig. 2.19, it can be deduced that the limits of the visibility swath are tangent on every point of a certain region of visibility produced by the visibility cone model. Nevertheless, due to the complexity of tangency computation on every single visibility region, a simpler model has been chosen and the definition used for the scan sector, although providing good results for combinations of small sensor aperture ( FOV ) and aiming angle ( Γ ), does not fulfill the tangency condition between regions and swath under certain extreme conditions. (a) FOV =10o,Γ=50o.(b) FOV =30o,Γ=50o. Figure 2.20 Swath’s limits tangency error. It is important to mention that this algorithm is not prepared to compute the swath on situations where the cone of visibility intersects only partially Earth’s sphere. These situations often involve large aperture angles of the sensor and/or angles of aim not small enough to generate a proper surface image, hence they are not the scope of this project. Geographical frame description Once the visibility limits are described in the normalised frame ( norm ), the coordinates of their points are changed to the geographical frame ( gph ) in order to locate the visibility regions with respect to other geographical data present on Earth’s surface. 2.6 Observations 29 These changes of reference frames are guided by algebraic transformation which apply the translation vectors and rotation matrices already described on this document (Sec. 2.1). Some results obtained for different latitudes are shown below. The orbital parameters for the EOS Satellite no.1 from the Planet Labs’ constellation Flock 2P have been used, Dateo f E poch : 22/08/2019 14 : 00 : 00.001 B?=0.22087 ×10−3R−1 ⊕ i=97.3677o Ω=299.1570o e=0.0010827 ω=152.8799o M=308.5990o n=15.26243907 rev day as well as a hypothetical value for the aperture of the sensor FOV =10o. In stead of a defined objective to observe on Earth’s surface, different fixed values for roll and pitch angles have been used. All simulations correspond to time windows defined by their Universal Time on 24 August, 2019. •Western Europe 10:30:00 - 10:40:00 To begin with, the results for the mentioned observation satellite as it flies over western Europe are presented. As can be seen, at medium latitudes not so far from the equator, the swath’s limits are practically parallel to the satellite’s ground-track, whereas approaching the pole they start to be distorted. ω=0o,ϕ=0o Figure 2.21 Nadir aiming observation. 30 Chapter 2. Preliminary considerations (a) ω=30o,ϕ=0o.(b) ω=−30o,ϕ=0o. Figure 2.22 Rolling observation. (a) ω=0o,ϕ=30o.(b) ω=0o,ϕ=−30o. Figure 2.23 Pitching observations. 2.6 Observations 31 (a) ω=30o,ϕ=30o.(b) ω=−30o,ϕ=30o. (c) ω=−30o,ϕ=−30o.(d) ω=30o,ϕ=−30o. Figure 2.24 Rolling and pitching observations. 32 Chapter 2. Preliminary considerations •North Pole 11:55:00 - 12:05:00 In order to study the behaviour of the algorithm around the poles, this extreme situation has been simulated. ω=−55o,ϕ=0o Figure 2.25 Pole coverage. As it is shown, the procedure taken to define the limits of the visible regions and swath has the capability to properly perform around the singularity of the poles •South Pole 11:05:00 - 11:10:00 An additional exposition of observation results is done for extreme latitudes, in this case, over the south pole by analysing different satellite rolling attitudes. (a) ω=0o,ϕ=0o. (b) ω=30o,ϕ=0o.(c) ω=−30o,ϕ=0o. Figure 2.26 South pole observations. 2.7 Satellite-Antenna transmissions 33 2.6.2 Modes Concerning the attitude configuration for the analysis of the scheduling project under the scope of this project, once the observation window of an objective region from a satellite begins the attitude is fixed until it finishes the capture. After the observation concludes, the satellite is allowed to accomplish a slew manoeuvre in order to change its attitude. Nevertheless, while the available attitudes configuration are not continuous, in order to address the problem by computational methods a discretisation is needed. Consequently, each attitude configuration involving specific values for roll, ω , and pitch, ϕ , angles is identified as a mode of observation of the satellite; while a finite number of modes are established to create the observation swaths of every satellite. As a consequence of considering different modes of observation for the same satellite, the algorithm should consider the presence of observations which are nor compatible between each other. The aspect of incompatibility between windows will be discussed in following sections. 2.7 Satellite-Antenna transmissions While there still is free space, the acquired observation by the satellites are stored in their RAM. These memories have a maximum storage capacity, which implies that the data must be transmitted to the ground stations, not only for making it accessible, but also for allowing the satellite to continue with its observations. Therefore, the transmissions consist on radio communication between a satellite and a ground station, both provided with the necessary instrumentation, such as the antennas. To establish the communications certain conditions must be fulfilled, mostly related with the visibility conditions. As the satellite orbits, a transmission with certain ground station is available while it flies over a zone where there is good visibility between both agents, i.e. when the satellite is within the ground station’s line of sight. The requirement compliance might be affected not only by the relative position between the ground station, the satellite and the surface of Earth, but also due to other factors such as meteorological phenomena. Elevation for transmission Although the meteorological phenomena remains out of the scope of this project, the atmospheric interference is taken into account. Therefore, a transmission will be available granted that the elevation, γ , of the satellite from the ground station location exceeds a minimum value, so that the microwaves travel along a shorter path inside the atmosphere. Figure 2.27 Minimum elevation for transmission. 40 Chapter 3. Computational geometry modelling observation windows are merged into just one. OWnew ←/0 ; OW0←OWold ; while OW06=/0 do j0←j∈OW0; OW00 ←OW0−j0; OWnew ←OWnew ∪j0; while OW00 6=/0 do j00 ←j∈OW00 ; ttot ←max{tf j0,tf j00 }− min{ts j0,ts j00 }; tuse ←(tf j0−ts j0)+(tf j00 −ts j00 ); if ttot −tuse ≥tmin then OW00 ←OW00 −j00; else j?←j0∪j00; ts j?←min{ts j0,ts j00 }; tf j?←max{ts j0,ts j00 }; OW00 ←/0; OWnew ←OWnew −j0; OWnew ←OWnew ∪j?; end end OW0←OW0−j0 end Algorithm 1: Observation windows joining algorithm 3.3 Transmission windows As the observation windows, the transmission ones are as well defined by their start an finish points, governed by the only considered restriction of satellite minimum elevation (Sec. 2.7). In order to identify each one of the edges and define the satellite’s attitude for each window, the next procedure has been followed: 1. For each satellite, its positions during time are analysed and checked whether the elevation relative to the any of ground stations fulfills the restriction established or not. 2. Once the transmission availability with each ground station is characterised for every point of the satellite’s orbit, the first and last points where the transmission is available for each interval are established as the window’s start and finish points. 3. For the start and finish instants of each window the relative positions of the ground station and satellite are used to calculate its pointing attitudes. 4. Each transmission window is described by the satellite and ground station involved and the start and finish times, ts and tf , while the start and finish attitudes will be employed to determine the slew times to accomplish the transmission. Idle time The transmission window is defined in terms of transmission availability regarding the relative positions between receiver and transmitter. Nonetheless, as a result of the transmission data-rate and the data hypothetically stored in the satellite’s memory, all data might be transmitted before the transmission window closes. In this case, the time between the end of the data transmission, t , and the finish of transmission window, tf , is defined as the idle time of the transmission, tidle. tidle =tf−t;ts≤t≤tf. During this intervals of time, the transmissions are still available but no needed. In order to implement the foresee-download policy in the scheduling algorithm [ 14 ], the idle periods of time must be take into account to state which transmission windows are not completely exploited. 3.4 Windows incompatibility 41 Figure 3.7 Transmission windows between a satellite and a ground station given the minimum elevation, γmin . 3.4 Windows incompatibility The established windows, observation or transmission ones, are defined for any available occasion between satellites and observation regions or ground stations. Nevertheless, not every combination of selected windows for the scheduling solution is feasible. As an example a satellite cannot acquire two different observations happening simultaneously. The incompatibility of the windows is strongly related with the time management on every source, ground station or satellite. Not only one observation or transmission operation is allowed at each instant, but also the auxiliary operations must be considered, such as the satellite slew operations or the ground stations resetting. Hence the incompatibility must be identified for every combination of windows. Observations incompatibility For two observations, {1,2} , performed by the same satellite a time is required for the slew operation, tslew 1,2 (2.53). Therefore the compatibility of both observation windows is defined by the following inequation: max{tf 1,tf 2}− min{ts 1,ts 2} ≥ (tf 1−ts 1)+(tf 2−ts 2)+tslew 1,2.(3.1) Transmission incompatibility In case of two transmissions, {1,2} , received by the same ground station, the resetting time of the ground station must be obeyed (Sec. 2.5), treset . The compatibility condition between both transmission windows is established as: max{tf 1,tf 2}− min{ts 1,ts 2} ≥ (tf 1−ts 1)+(tf 2−ts 2)+treset .(3.2) Mixed incompatibility Not only the slew operations are needed in order to accomplish two different observations with different attitudes; when a transmission is started or finished the attitude of the satellite is such that its antenna points the ground station involved. Hence a slew manoeuvre will be carried out between observations and transmission windows. Three situations may be presented: 42 Chapter 3. Computational geometry modelling •Observation before transmission (tf obs ≤ts trs) : In this case the attitudes considered for the slew operation are the attitude for the observation, {1} , and the attitude at the transmission start, {2s} . The compatibility restriction is described as follows: max{tf 1,tf 2}− min{ts 1,ts 2} ≥ (tf 1−ts 1)+(tf 2−ts 2)+tslew 1,2s.(3.3) •Transmission before observation (tf trs ≤ts obs) : For this situation the attitudes considered to establish the slew operation are the attitude at the transmission finish, {1f} , and the attitude for the observation, {2}. The compatibility restriction is described as follows: max{tf 1,tf 2}− min{ts 1,ts 2} ≥ (tf 1−ts 1)+(tf 2−ts 2)+tslew 1f,2.(3.4) •Simultaneous observation and transmission : Under this situation the windows are inherently incompatible. The simultaneity of both windows is described by the expression max{tf 1,tf 2}− min{ts 1,ts 2} ≤ (tf 1−ts 1)+(tf 2−ts 2)(3.5) 3.5 Subregions Considering one observation window related to an objective region (Fig. 3.8), this is divided into sub-parts which can be classified into two groups: •Covered zones : The zones born from the intersection between the objective region and the observation are covered by the observation itself. •Uncovered zones: The zones defined by the subtraction of the observation from the objective region are not covered by the considered observation. Figure 3.8 Objective region divided by observation . Based on the previous idea, examining the interactions between all the objective regions and the entire number of observation windows generated, the objective region is divided into a series of distinct zones. This regions will be named subregions (SR), and each one of them will be covered by certain satellites with certain modes after certain number of orbit revolutions. The objective of the optimization algorithm is to determine which subregion is covered by which observation window in order to fulfill some specific constraints and optimize a certain feature of this coverage. 3.5 Subregions 43 Subregions algorithm In order to define the geometry of the subregion a recursive algorithm has been used in which, starting from the original objective region, one at a time, every observation window is considered, and the new subregions generated are included to the previous subregion set. (a) First considered observation window. (b) Second considered observation window. Figure 3.9 Subregions computation algorithm. As seen in Fig. 3.9, the algorithm takes the existing subregions set at the previous step and define the new one by carrying out boolean operations between the new observation window and the subregions. This algorithm is better explained in algo. 2. Meanwhile the algorithm build the new subregions, these ones are characterized with information about the objective region to which they belong as well as the satellites an attitudes used for covering them. In addition the coverage matrix, Qi j , defined later in this document, containing information about the coverage relation between subregions and observation windows is set up. SRnew ←/0 ; while SRold 6=/0 do i?←i∈SRold ; Qnew i j ←Qold i j ∀i6=i?,j; Inters ←i?∩j?; while Inters 6=/0 do i0←i∈Inters ; SRnew ←SRnew ∪i0; Qnew i0j←Qold i?j∀j; Qnew i0j?←1; Inters ←Inters −i0; end Subtr ←i?−j?; while Subtr 6=/0 do i0←i∈Subtr ; SRnew ←SRnew ∪i0; Qnew i0j←Qold i?j∀j; Qnew i0j?←0; Subtr ←Subtr −i0; end end Algorithm 2: Subregion set, SR , and coverage matrix, Q , update given a new considered observation window, j?. Since the boolean operation must be executed between the new observation window and the whole number of pre-existing subregions, the complexity of this algorithm increases exponentially with the number of 44 Chapter 3. Computational geometry modelling observation windows already considered. This might raise the computation time for large objective regions in which a great number of observation windows are involved. With the objective of avoiding this problem, a divide-compute-assemble process is implemented for each objective region: 1. Maximum values for latitude height, ∆φmax , and longitude width, ∆λmax , are set depending on the orientation of the observation windows. The objective region is checked to fit into the geographical rectangle established by these values and, in case not, it is divided into parts which fulfill the maximum width and height guidelines. 2. The algorithm described (Fig. 2) is executed for each one of the parts in which the objective region has been divided. (a) Division of the regions into smaller auxiliary ones. (b) Subregion computation for each auxiliary region. Figure 3.10 Subregions computation process. 3. Once the subregions are defined for each part, an assemblage of the parts occurs. The contiguous subregions of different parts are merged into one subregion by the union boolean operation. This process is repeated for every objective region participating in the problem. In Fig. 3.11 a final result of the algorithm is presented for the Iberian Spanish territory and Balearic Islands Figure 3.11 Subregions for certain objective regions and satellites. 4 Problem optimization A final mathematical description about the model used to address the observations and transmissions scheduling problem as well as the algorithms and tools used to solve it are presented in this chapter. 4.1 Mathematical description For the resolution of the problem is convenient to thoroughly describe it before. The scheduling problem, as an optimization problem, is governed by a number of constraints which must be obeyed and an objective to optimize. The constraints define the set of feasible solution: schedules which fulfill a series of previously stated requirements, while the objective determines the best schedule (or schedules) from the feasible solution set. For the observation and transmission scheduling problem constraints might take into account limitations such as the available energy of the satellites to perform slew and observation operations, whereas the objective might consist in minimizing the transmission of data or the total cost of the operations. In order to address the problem with the optimization algorithms, a mathematical model of these constraints and objective must be done. Once the parameters and variables are precisely described, the mathematical form of the problem goal consist on a objective function of the variables governed by the parameters which has to be maximized; and the constraints are expressed as equalities or inequalities related to the variables and governed by the parameters as well. In particular, the observation and transmission scheduling problem is typically modelled in terms of integer variables, more specifically, most of them binary data, indicating dual information about the entities (covered or not, scheduled or not, compatible or not...). A description of this model used for the particular problem of observation-transmission problem under the scope of the project is made as follows: 4.1.1 Mathematical entities Time interval The scheduling problem is considered to be solved during a time interval characterized by initial and final instants, t∈[t1,tend ] . The magnitude of the time span, tspan =tend −t1 , characterises the scale of the scheduling problem, allowing a greater number of objective regions observed but less rescheduling responsiveness for longer time spans and vice versa. A discretisation of the continuous time interval is made to address the problem. The time step, tstep , establishes the distance between the discrete time instants and, as a consequence, the number of subintervals between each discrete instants given a time span, Nsubi =&tspan tstep ', as well as the number of time instant when the state is defined, N=Nsubi +1, 45 46 Chapter 4. Problem optimization and therefore the discrete instants, t1≤tτ≤tend :τ=1,2, ... , N. Hence, an integer time parameter, τ, is defined for the complete time characterisation of the problem: τ=&tτ−t1 tstep +1':τ∈Z,0≤τ≤N.(4.1) Objective regions set The union of the R disconnected regions, r , whose observation is required establishes the objective regions set, Ob jReg: [ ∀r r=Ob jReg. Ground stations set The union of the G ground stations, g , considered for the data transmission forms the ground stations set, GrdSt: [ ∀g g=GrdSt. Observation satellites set The union of the S satellites, s , participating in the capture and transmission of data defines the observation satellites set, ObsSat: [ ∀s s=ObsSat. Modes set The Mfixed modes, m, used by the satellite to capture Earth’s surface form the modes set, Modes: [ ∀m m=Modes. Observation configurations set Based on the previous definitions, each combination of satellite, s , and mode of observation, m , generates an specific observation configuration, k . The K observation configurations, k , generate the Observation configurations set, ObsCon f : [ ∀k k=ObsCon f . Hence, a direct relation between the observation configurations and the observation satellites and modes of observation is established: k⇔s,m. Observation windows set The observation windows set, OW , is composed by the J observation windows, j , produced thanks to the combinations of the previous parameters: [ ∀j j=OW. A given observation window is determined by the e -th observation chance of the r objective region by the k observation configuration, involving the satellite s=s(k) and observation mode m=m(k) . Given the previous entities the start, τs rke , and finish, τf rke , times are extracted for each observation window j . Therefore: j={[τs rke,τf rke]:r∈Ob jReg,k∈ObsCon f ,e∈Z}.(4.2) 4.1 Mathematical description 47 Therefore, a direct relation for each observation window with the objective region, observation configuration and observation chance is established: j⇔r,k,e. The observation windows subset generated by a certain satellite, s , definition will be useful for the resolution: [ j=j(s) j=OWs Transmission windows set The transmission windows set, TW, is composed by the Dtransmission windows, d: [ ∀d d=TW. A certain transmission window is established by the e -th transmission chance between the g ground station and the s satellite. Given the previous entities the start, τs gse , and finish, τf gse , times are extracted for each transmission window d. Therefore: d={[τs gse,τf gse]:g∈GrdSt,s∈ObsSat,e∈Z}.(4.3) Similar to the observation windows case, a direct relation for each transmission window with the ground station, observation satellite and transmission chance is established: d⇔g,s,e. Other characteristic transmission windows subsets are defined: •Satellite transmission windows : formed by the transmission windows involving a certain satellite, s : [ d=d(s) d=TWs. •Ground station transmission windows : formed by the transmission windows involving a certain ground station, g: [ d=d(g) d=TWg. Subregions set The subregion set, SR , is composed by the I subregions, i , created by the intersection of all the observation window, j, with the objective regions, r: [ ∀i i=SR. It is useful to define additional subregions subsets: •Covered subregions: formed by the subregions observed at an observation window, j: [ ∀i∈j =SRj. •Satellite subregions: formed by the subregions observed by a satellite, s: [ ∀i∈OWs =SRs. 48 Chapter 4. Problem optimization Acquisitions set The acquisitions set, Acq, is composed by the scheduled observation windows, j?: [ ∀j=j? j=Acq;j∈OW. Once the problem is finished, the acquisition set represent the solution for the observation scheduling problem. Downloads set The downloads set, Dwn, is composed by the scheduled transmission windows, d?: [ ∀d=d? d=Dwn;d∈TW. Once the problem is finished, the downloads set represent the solution for the transmission scheduling problem. 4.1.2 Mathematical relations Observation windows coverage The observation windows coverage is modelled by the coverage matrix, Qi j , which states whether the i subregion is covered by the jobservation window or not. Qi j =1if i∈j 0if i/∈j;i∈SR,j∈OW.(4.4) Windows incompatibility The time incompatibility between windows is modelled using a series of matrices. These matrices states whether two generic windows are incompatible, i.e. if both windows are allowed to be scheduled together (Sec. 3.4). •Observation windows incompatibility (Robs j1,j2): establishes whether the observation windows j1 and j1are incompatible or not. Robs j1j2=1if j1,j2incompatible 0if j1,j2compatible ;j1,j2∈OW.(4.5) •Transmission windows incompatibility(Rtrs d1,d2): establishes whether the transmission windows d1 and d2are incompatible or not. Rtrs d1d2=1if d1,d2incompatible 0if d1,d2compatible ;d1,d2∈TW.(4.6) •Mixed windows incompatibility(Rmix j,d): establishes whether the observation window j and the transmission window dare incompatible or not. Rmix jd =1if j,dincompatible 0if j,dcompatible ;j∈OW,d∈TW.(4.7) 4.1.3 Variables Observation scheduling The observation scheduling variable, xj , states whether an observation window j is adopted as an acquisition for the solution or not: xj=1if j∈Acq 0if j6∈ Acq ;j∈OW.(4.8) 4.1 Mathematical description 49 Transmission scheduling The transmission scheduling variable, yd , states whether an transmission window d is adopted as a download for the solution or not: yd=1if d∈Dwn 0if d6∈ Dwn ;d∈TW.(4.9) Storage The storage variable, Lsτ , describes the storage level of the satellite s at the time τ . The storage is measured by the time of observation still not transmitted gathered in memory, T, at the specific time. Lsτ=Tsτ;s∈ObsSat,τ∈Z,0≤τ≤N.(4.10) 4.1.4 Constraints Complete coverage The complete coverage of the objective regions must be fulfilled. To address this the constraint which states that each subregion, i, must be covered at least once by the acquisition set, Acq, is defined: ∑ ∀j∈OW Qi jxj≥1; ∀i∈SR.(4.11) This constraint in some cases might be very restrictive, typically on analysis where the time span, tspan , is not big enough, the area of observation regions, Area(S∀r∈ObsReg r) , is too big and/or the observation parameters do not allow to get a feasible solution for this constraint. In these situations and when the resolution procedure allows it, the problem description will hold a relaxation of this constraint in which the coverage of the maximum number of subregions is pursued. Schedule compatibility The windows chosen for the schedule solution must not interfere between them, every window incorporated to the solution must be compatible with the other ones. This constraint is defined for each scheduled observation and transmission windows as follows: •Acquisitions: Every scheduled observation, j? , must not be incompatible with the rest of the scheduled observations and transmissions: ∑ ∀j∈OW Robs j?jxj<1; ∀j?∈Acq.(4.12) ∑ ∀d∈TW Rmix j?dyd<1; ∀j?∈Acq.(4.13) •Downloads: Every scheduled transmission, d? , must not be incompatible with the rest of the scheduled transmissions and observations: ∑ ∀d∈TW Rtrs d?dyd<1; ∀d?∈Dwn.(4.14) ∑ ∀j∈OW Rmix jd?xj<1; ∀d?∈Dwn.(4.15) Storage range The storage of each satellite, s , must remain within the safe range at every instant of time, τ . The safe range is defined by it minimum value, set to zero, and a maximum value, MaxStg. 0≤Lsτ≤MaxStg;∀s∈Sat,∀τ∈Z,0≤τ≤N.(4.16) 56 Chapter 5. Results 5.1.2 Ground stations Regarding the ground stations, the locations shown on Fig. 2.13 and described in Table 2.2 are used, with the exception of no. 6 in Antarctica, which is no longer available. In addition,a total of 4 antennas are located on each site, establishing a total number of 28 antennas considered for every scenario. Moreover, the parameters involved in the communications with the ground stations are established as well thanks to the indications of Planet Labs. The reset time of each antenna, treset and the satellite’s minimum elevation required, γmin, are defined bellow: treset =4min γmin =10 o For the transmission-observation speed ratio value, δtrs (2.78) , some estimations have been made about the size of the observation area [5] and the downlink speeds [2]. Observed area ratio ≃250 km2 s Data size Observed area ≃15 kB km2 ⇒DRobs ≃3750 kB s DRtrs ≃220Mbps ⇒δtrs ≃7.5 5.2 Scenarios Several scenarios are considered for the previous resources varying in the objective regions considered. For each of them the computational demands for the instances definition and the final resolution are exhibit. Instances The following entities are computed: •Observation windows: After the definition of the swaths and the objective region, a total number of NOW observation windows are computed. From all of them, the only ones occurring at daytime, NOWD , are selected. The computing time needed for this algorithm is provided, TOW . •Subregions: Next, the objective region is divided into NSR subregions and the coverage matrix is calculated. The computing time needed for this process is provided, TSR. The computational costs of this process strongly depends on the morphology of the observation windows considered. The number of subregion increases rapidly when working with observation windows that cross between them, while if they remain parallel less subregions are created and the cost of the computation decrease. In the particular case of the sun-synchronous orbits with which the scenarios are solved, the pair of different directions of the observation windows for each satellite belong to the two observation chances during the whole day, of which only one occurs at day time. The constellations selected have parallel directions of their observation for the daytime passes (Fig. 2.12), hence the number of subregions, NSRD , and the computing time, TSRD are much lower than the values for the consideration of both passes, NSRtot and TSRtot . •Transmission windows: Afterwards, the number of transmission windows between each satellite and each ground station antenna, NTW , is computed. The computing time needed for this process is provided, TTW . •Incompatibility matrices: Finally, an algorithm computes the three incompatibility matrices for each observations and/or transmissions combinations. The computing time needed for this algorithm is provided, TR. Resolution Once the Instances for each scenario are set, the resolutions are carried out for both observations sub-problem and observations and transmissions complete problem. Different methods are employed for the resolution depending on which problem is addressed: 5.2 Scenarios 57 Observations sub-problem •Matlab ILP solver: The intlinprog.m function of the software an be used thanks to the linearity of this sub-problem. The algorithm first look for the optimal solution of the continuous linear problem (LP) and provides the value of its objective function, f? LP . Afterwards, it applies different kind of cuts in order to obtain a lower bound of the objective function value for the integer linear problem (ILP), flb . Ultimately, it performs the computation of the optimal solution for the integer linear problem thanks to aBranch and Bound method and the optimal objective function value is obtained, f? . The computing time taken for this resolution is given, T?. •Observations heuristic algorithm: The previously mentioned heuristic algorithm (Algo. 3) is used as an intermediate step between the optimal solution of the observations sub-problem and the heuristic one of the complete problem. The algorithm is executed for different sizes for the candidates subset, p , and number of repetitions, Nsol , in order to obtain different heuristic solutions, fheu , to compare. As previously mentioned (Subsec. 4.2.1) the Algo. 3 might not fulfill the full-coverage constraint (Eq. (4.11) ). The decision made in this project is to accept the results as feasible even though portions of the objective regions is not covered, and to evaluate the size of the area from the objective regions which remains uncovered, Auncov rel =Area not covered Area o f ob jective region. In case a number of solution is generated, as in GRASP algorithm, the process also provides the best solution which covers the entire regions if it is found. On each case the computing times for the resolutions are evaluated, Theu. Observations and transmissions sub-problem •Observations and transmissions heuristic algorithm: On the complete problem, the only method used is the observations and transmissions heuristic algorithm designed (Algo. 4). As done for the heuristic resolution of the observations sub-problem, the algorithm is executed for different sizes for the candidates subset, p , and number of repetitions, Nsol , in order to obtain different heuristic solutions, fheu , to compare. In addition, index Auncov rel is presented as for the observations sub-problem algorithms. The computing time for each case is evaluated, Theu. As introduced in Sec. 4.2 and in the previous paragraphs, the way the heuristic algorithms are designed for this project also requires the previous computation of the subregions, NSR , in which the objective regions are divided. This fact substantially increase the overall necessary computing time with respect to other heuristic resolutions which avoid the computation of the mentioned subregions; specially on scenarios where a variety of constellations is used and observation windows which cross between them are generated, rising the number of the created subregions. Nevertheless, the objective of the designed heuristic algorithms, and in particular the of Algo. 4, is to provide the closest solution to the optimal one for the observations and transmissions scheduling problem, since the exact algorithms are not able to address this complete problem. Therefore, in order to contextualise the results obtained from the heuristic algorithm used for the resolution of the complete problem and link them with the exact resolution of the observations scheduling sub-problem, the results provided by the heuristic algorithm solving the observations sub-problem (Algo. 3) are as well presented, representing a middle step between both main resolutions of sub-problem and complete problem. 5.2.1 SPAIN In this case, the political territories of Spain are considered as objective regions for the resolution of the scheduling problems. Aob j =5.1602 ×105km2. 58 Chapter 5. Results Figure 5.1 Spanish territories. Instances The parameters evaluating the problem instances computation performance for this particular scenario are shown in Table 5.1. (a) Subregions created by all the observation windows. 5.2 Scenarios 59 (b) Subregions created by daytime observation windows . Figure 5.1 Subregions created . Table 5.1 Instances computation for Spanish territories coverage . Observation windows Subregions Transmission windows Incompatibility matrices Entities NOW =976 NOWD=488 NSRtot =59611 NSRD=1359 NTW =12668 - Computing time [s]TOW =51.246 TSRtot =2165.456 TSRD=25.412 TTW =67.781 TR=0.741 60 Chapter 5. Results Resolution The resolutions for both sub-problem and complete problem and the performance of the algorithms used are finally analysed for the scenario with the Spanish territories as objective regions. •Observations scheduling sub-problem (c) Scheduled observations. (d) Gantt chart for acquisitions. Figure 5.2 Optimal solution provided by MATLAB ILP solver for the observations scheduling sub-problem concerning the coverage of the Spanish territories. 5.2 Scenarios 61 Table 5.2 Spain coverage observations sub-problem: Algorithms performances. MATLAB ILP solver f? LP [acq.]62.000010 flb [acq.]63 f?[acq.]63 T?[s]0.151 GREEDY algo. Partial coverage allowed fheu [acq.]72 Auncov rel [-] 0 Theu [s]0.377 GRASP algo. Partial coverage allowed (p=2) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Nsol =500 Full coverage No Yes No Yes No Yes No Yes fheu [acq.] 70 - 68 73 66 73 66 71 Auncov rel [-] 0.465% -0.491% 01.082% 03.237% 0 Theu [s] 1.904 5.786 11.677 56.922 GRASP algo. Partial coverage allowed (p=4) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Nsol =500 Full coverage No Yes No Yes No Yes No Yes fheu [acq.] 69 74 67 70 66 70 67 69 Auncov rel [-] 3.910% 06.348% 06.818% 01.084% 0 Theu [s] 1.290 6.020 11.561 57.557 GRASP algo. Partial coverage allowed (p=6) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Nsol =500 Full coverage No Yes No Yes No Yes No Yes fheu [acq.] 68 70 68 72 65 71 65 70 Auncov rel [-] 2.837% 00.644% 01.698% 03.047% 0 Theu [s] 1.795 6.023 12.047 57.478 GRASP algo. Partial coverage allowed (p=8) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Nsol =500 Full coverage No Yes No Yes No Yes No Yes fheu [acq.] 67 - 67 - 67 71 65 70 Auncov rel [-] 1.175% -2.625% -3.459% 04.739% 0 Theu [s] 2.502 6.171 11.849 57.952 GRASP algo. Partial coverage allowed (p=10) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Nsol =500 Full coverage No Yes No Yes No Yes No Yes fheu [acq.] 67 - 67 73 67 70 65 71 Auncov rel [-] 1.606% -5.352% 03.578% 02.605% 0 Theu [s] 1.809 6.198 11.856 58.351 62 Chapter 5. Results •Observations and transmissions scheduling problem Table 5.3 Spain coverage observations and transmissions problem: Algorithms performances. GREEDY algo. Partial coverage allowed fheu [acq.]72 Auncov rel [-] 0 Theu [s]20.753 GRASP algo. Partial coverage allowed (p=4) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Full coverage No Yes No Yes No Yes fheu [acq.] 68 72 67 71 67 71 Auncov rel [-] 1.829% 02.740% 02.305% 0 Theu [s] 108.021 514.922 987.790 GRASP algo. Partial coverage allowed (p=6) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Full coverage No Yes No Yes No Yes fheu [acq.] 70 - 67 71 68 72 Auncov rel [-] 0.832% -4.758% 01.019% 0 Theu [s] 105.097 510.122 984.468 GRASP algo. Partial coverage allowed (p=8) Solutions obtained Nsol =10 Nsol =50 Nsol =100 Full coverage No Yes No Yes No Yes fheu [acq.] 67 - 67 72 66 71 Auncov rel [-] 0.791% -3.069% 02.843% 0 Theu [s] 104.390 104.390 986.837 The solution provided by the GRASP algo. which randomly chooses between p=8 candidates and provides Nsol =10 has been selected to illustrate the results obtained. In this case, 99.2% of the total area is covered by 67 acquisitions. In this solution, 62 downloads are scheduled to gather the data acquired before the end of the day. (a) Scheduled observations. 5.2 Scenarios 63 (b) Gantt chart for acquisitions . (c) Gantt chart for acquisitions and downloads . (d) Storage evolution on satellites. Figure 5.2 Solution provided by GRASP heuristic algo. (p=8,Nsol =10) for the observations and transmissions scheduling problem concerning the coverage of the Spanish territories. 64 Chapter 5. Results 5.2.2 SPAIN AND ITALY In the next scenario, the Italian political territories, together with the Spanish ones, are considered in order to increase the size of the problem. Aob j =8.3654 ×105km2. Figure 5.3 Spanish and Italian territories. Instances The parameters evaluating the problem instances computation performance for this particular scenario are shown in Table 5.4. Figure 5.4 Subregions created by daytime observation windows. Table 5.4 Instances computation for Spanish and Italian territories coverage . Observation windows Subregions Transmission windows Incompatibility matrices Entities NOW =1968 NOW D=981 NSRD=2579 NTW =12668 - Computing time [s]TOW =85.821 TSRD=53.780 TTW =76.340 TR=0.826 5.2 Scenarios 65 Resolution The resolutions for both sub-problem and complete problem and the performance of the algorithms used are finally analysed for the scenario with the Spanish territories as objective regions. •Observations scheduling sub-problem (a) Scheduled observations. (b) Gantt chart for acquisitions. Figure 5.5 Optimal solution provided by MATLAB ILP solver for the observations scheduling sub-problem concerning the coverage of the Spanish and Italian territories. 72 Chapter 6. Conclusions and future work scheduling of constellations with transmit-while-observe capabilities or with large storage capacities which relax the downloading demands. Finally, the resolution of the scheduling problems is carried out using different methodologies. For the observations scheduling sub-problem the performance of the MATLAB ILP solver make the designed heuristic algorithm uninteresting, since it is capable to reach the exact solution in less computing time. Regarding the complete problem, the heuristic algorithm designed successfully accomplishes its objective and places the required downloads in order to allow the scheduling of more observations when the storage capacity constraint block them. In addition, although the algorithm occasionally is not able to find solutions which full-cover the objective regions, the areas remaining uncovered are small in comparison with the total area of the objective. An analysis for the different combinations of the number of candidates, p , to join the solution on each iteration of the algorithm, and the number of solutions generated, Nsol , has been established for each scenario. This analysis shows that when more candidates are considered, more number of solutions must be generated to get better solutions, substantially increasing the computing time of the resolution. Furthermore, the rise of the parameter p complicates the detection of full-coverage solutions. Therefore, small number of candidates, p, are preferred to decrease the time needed to obtain a proper schedule. Improvements and future work Here are now presented a series of considerations which have not been addressed in this project but might be interesting for future studies: • As mentioned, in this project, the Earth’s surface is modelled as a perfect sphere in order to determine the swaths of the satellites. This fact might produce inaccurate coordinates for the limits of the visibility region. Hence, a better modelling of the Earth is proposed using an ellipsoid in stead of a sphere, e.g. the WGS-84 standard ellipsoid. • In the approach taken in this document, the available attitudes take discrete values and the slew operations of the satellites are only allowed while they are not observing. The consideration of continuous values of attitude and their capability to slew and observe simultaneously might be interesting for full-coverage better solutions. • As previously introduced, the resolution carried out address the problem concerning specific regions to be covered. A bigger scale analysis might be performed for the scheduling pursuing the complete coverage of the Earth land mass. • For future projects, other considerations constraining the scheduling problem could be analysed. Some of them might include the duty cycle constraints for satellites and ground stations or possible failed observations and transmissions due to meteorological or misalignment uncertainties. • The procedure taken on the project detects the transmission windows and their incompatibility between the rest of the operation windows. Then, once one transmission is scheduled in the resolution, all their incompatible windows are discarded even though only a brief downloading lapse is used from the entire transmission window. A more optimised solution might be achieved considering this fact and positioning the download interval at the best moment within the transmission window. • In stead of analysing a transmit-or-observe policy in which no simultaneous observations and transmissions are allowed, future projects can address problems concerning transmit-while-observe constellations, removing the observation-transmission incompatibility constraint. • Due to the growth in the amount of satellites in LEO orbit, space agencies are developing constellations of satellites whose task is to relay the information and data between the entities involved, such as the European Data Relay System (EDRS). For future concerns, this constellations might substitute the ground stations for better transmission availabilities. • In terms of the algorithms employed for the resolution, more efficient heuristic algorithms, which do not need the computation of elements needed for the exact resolution, could be used in order to improve their time performances. In addition, commercial integer programming solvers (Gurobi, SCIP, ...) might be also used to solve the coupled observations and transmissions scheduling problem. List of Figures 1.1 Amount and forecast of CubeSats’ launches per year 1 1.2 Planet Labs’s Dove 2 2.1 Geocentric reference frame [1] 5 2.2 Perifocal and geocentric equatorial frames [1] 6 2.3 Geographic reference frame [1] 7 2.4 VVLH reference frame 8 2.5 Gravitational forces [1] 11 2.6 Orbital plane description [1] 13 2.7 Eccentric anomaly definition [1] 13 2.8 US Standard Atmosphere 1976 [1] 15 2.9 Two Line Element set format [12] 17 2.10 Region of interest modification 19 2.11 Region overlapping solution 19 2.12 Flock 2P and Flock 3P constellations 21 2.13 PlanetLabs ground stations spots and ground tracks of representative Doves satellites [2] 22 2.14 Simplified geometrical situation 23 2.15 Cone of Vision 24 2.16 Scan Sector 25 2.17 Limit of possible visibility 26 2.18 Visibility regions 27 2.19 Visibility swath and regions 28 2.20 Swath’s limits tangency error 28 2.21 Nadir aiming observation 29 2.22 Rolling observation 30 2.23 Pitching observations 30 2.24 Rolling and pitching observations 31 2.25 Pole coverage 32 2.26 South pole observations 32 2.27 Minimum elevation for transmission 33 3.1 Example 1 of polyshape object 35 3.2 Example 2 of polyshape object 36 3.3 Polyshape with outlier 36 3.4 Boolean operations 37 3.5 Swath polyshapes for a certain satellite 38 3.6 Observation windows for Iberian and Balearic Spanish regions, created by Flock 3P constellation 39 3.7 Transmission windows between a satellite and a ground station given the minimum elevation, γmin 41 3.8 Objective region divided by observation 42 3.9 Subregions computation algorithm 43 3.10 Subregions computation process 44 73 74 List of Figures 3.11 Subregions for certain objective regions and satellites 44 5.1 Spanish territories 58 5.1 Subregions created 59 5.2 Optimal solution provided by MATLAB ILP solver for the observations scheduling sub-problem concerning the coverage of the Spanish territories 60 5.2 Solution provided by GRASP heuristic algo. (p=8,Nsol =10)for the observations and transmissions scheduling problem concerning the coverage of the Spanish territories 63 5.3 Spanish and Italian territories 64 5.4 Subregions created by daytime observation windows 64 5.5 Optimal solution provided by MATLAB ILP solver for the observations scheduling sub-problem concerning the coverage of the Spanish and Italian territories 65 5.5 Solution provided by GRASP heuristic algo. (p=4,Nsol =10)for the observations and transmissions scheduling problem concerning the coverage of the Spanish and Italian territories 68 List of Tables 2.1 Planet Labs’Flock 2P and 3P constellations 20 2.2 Planet Labs ground stations (2017) 22 5.1 Instances computation for Spanish territories coverage 59 5.2 Spain coverage observations sub-problem: Algorithms performances 61 5.3 Spain coverage observations and transmissions problem: Algorithms performances 62 5.4 Instances computation for Spanish and Italian territories coverage 64 5.5 Spain and Italy coverage observations sub-problem: Algorithms performances 66 5.6 Spain and Italy coverage observations and transmissions problem: Algorithms performances 67 75 Bibliography [1] Howard D. Curtis, Orbital mechanics for engineering students, Butterworth-Heinemann, 2013. [2] Kiruthika Devaraj, Ryan Kingsbury, Matt Ligon, Joseph Breu, Vivek Vittaldev, Bryan Klofas, Patrick Yeon, and Kyle Colton, Dove high speed downlink system, (2017). [3] Felix R Hoots and Ronald L Roehrich, Models for propagation of norad element sets, Tech. report, Aerospace Defense Command, Peterson AFB (CO): Office of Astrodynamics, 1980. [4] The MathWorks Inc., Computational geometry: Elementary polygons. [5] Bryan Klofas, Planet labs ground station network, Tech. report, 13th Annual CubeSat Developers Workshop, Cal Poly SLO, 2016, 13th Annual CubeSat Developers Workshop, Cal Poly SLO. [6] Ana Laura Zabala López, Análisis y desarrollo de algoritmos para la planificación optimizada de adquisición de imágenes terrestres, Master’s thesis, Escuela Técnica Superior de Ingeniería (Universidad de Sevilla), 2013. [7] Sergio Rubio Madroñal, Análisis y desarrollo de algoritmos para la planificación optimizada de misiones espaciales de observacion de la tierra, Master’s thesis, Escuela Técnica Superior de Ingeniería (Universidad de Sevilla), -. [8] NASA, Remote sensors | earthdata, Aug 2019. [9] Federico Perea, Rafael Vazquez, and Jorge Galan-Viogue, Swath-acquisition planning in multiplesatellite missions: an exact and heuristic approach, IEEE Transactions on Aerospace and Electronic Systems 51 (2015), no. 3, 1717–1725. [10] Planet Labs, Planet labs specifications: spacecraft operations and ground systems, 1.0 ed., 2015. [11] Jérôme Théau, Temporal resolution, pp. 1150–1151, Springer US, Boston, MA, 2008. [12] David Vallado, Paul Crawford, Ricahrd Hujsak, and TS Kelso, Revisiting spacetrack report# 3, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2006, p. 6753. [13] Rafael Vazquez, Federico Perea, and Jorge Galán Vioque, Resolution of an antenna–satellite assignment problem by means of integer linear programming, Aerospace Science and Technology 39 (2014), 567–574. [14] Pei Wang, Gerhard Reinelt, Peng Gao, and Yuejin Tan, A model, a heuristic and a decision support system to solve the scheduling problem of an earth observing satellite constellation, Computers & Industrial Engineering 61 (2011), no. 2, 322–335. 77