Analysis and optimization of snapshot-mode earth observation 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 Fin de Grado Grado en Ingeniería Aeroespacial Analysis and Optimization of SnapshotMode Earth Observation Constellations Autor: Julián Mandrión Ramírez 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, 2020 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Analysis and Optimization of Snapshot-Mode Earth Observation Constellations Autor: Julián Mandrión Ramírez 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, 2020 Trabajo Fin de Grado: Analysis and Optimization of Snapshot-Mode Earth Observation Constellations Autor: Julián Mandrión Ramírez 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: Resumen Análisis y Optimización de Constelaciones de Observación Terrestre en Modo Snapshot Palabras clave Satélites de observación terrestre, optimización, constelaciones, toma de imágenes, algoritmos heurísticos, planificación, satélites ágiles, camino de observación. Resumen El problema propuesto en este documento es la realización de una planificación optimizada para la toma de imágenes de una constelación formada por satélites ágiles. Este tipo de satélites presenta una mayor libertad para modificar su actitud con respecto a sus predecesores. Los objetivos de este proyecto fueron principalmente: minimizar el tiempo de observación de un cierto escenario, dada una nube de puntos objetivos en un área, y maximizar la calidad de las imágenes tomadas por la constelación. Para poder realizar simulaciones y obtener resultados se crearon dos modelos de realidad. El primero consiste en un cierto vehículo sobrevolando un plano, que serviría para ensayar los algoritmos creados, y el segundo es un modelo de una constelación de satélites heliosíncronos orbitando la Tierra. Se estudió el movimiento óptimo de giro de la cámara para el apuntamiento entre dos puntos. Se llegó a la conclusión de que este movimiento era una rotación alrededor del eje de Euler y se analizó como resolver este movimiento de forma numérica. Dada la imposibilidad de obtener la secuencia óptima de observación de un escenario de manera exacta, se desarrolló un Algoritmo de Selección de Objetivos (TSA) para buscar una secuencia con un tiempo de observación aceptable. Además un algoritmo de optimización ha sido probado, obteniendo de manera exitosa caminos de observación alternativos para el problema plano. Se han realizado simulaciones en varios escenarios globales y locales. En los escenarios globales se han estudiado las tendencias que presenta el número de puntos objetivo observados durante 24 horas de simulación. Por otra parte, se han simulado la observación de escenarios locales sobre la península ibérica. El objetivo de estas simulaciones era observar las tendencias del tiempo de observación total al variar ciertos parámetros. Además se ha medido el efecto de ciertos acontecimientos aleatorios como la presencia de nubes en el tiempo de observación. Conclusiones Las conclusiones que pueden extraerse de los resultados obtenidos son las siguientes: •Escenarios globales I II Chapter 0. Resumen – En este tipo de escenarios la manera de seleccionar los puntos no es relevante, ya que los satélites pueden observar todos los puntos que entran en cobertura. – Al aumentar la altitud ( h ) o el máximo ángulo de rotación de la cámara ( Ψ ) de los satélites , el número de objetivos observados sube hasta un máximo (el número de objetivos generados). Este incremento en el numero de puntos observado conlleva una reducción de la calidad de las imágenes tomadas. – Aumentar el número de satélites de la constelación siempre aumenta el número de puntos objetivo observados siempre que no se llegue al máximo. – Cuando una constelación observa todos los puntos generados y cubre toda la tierra, puede ocurrir en ciertas ocasiones que al incrementar el número de satélites de la constelación se reduzca de forma significativa el tiempo requerido para observar esos puntos. •Escenarios locales –Al aumentar Ψohde los satélites, el tiempo de observación del escenario disminuye. Además, este tiempo tobs toma valores que son aproximadamente múltiplos de 24 h . Este descenso en el tiempo de observación conlleva una pérdida de calidad de las imágenes tomadas. – Al aumentar el número de satélites presentes en una constelación, el tiempo de observación de los escenarios generalmente se reduce. – Ninguno de los dos métodos de selección de puntos estudiados parece tener una especial ventaja sobre el otro. En ciertas ocasiones "Sort by distance" funciona mejor, mientras que en otras "Sort by topt " resulta más adecuado. – El tiempo de observación de los escenarios se incrementa de forma significativa si se tienen en cuenta puntos objetivos de observación repetida. Por esta razón es de suma importancia tener en cuenta los fenómenos aleatorios que fuerzan a observar los objetivos en varias pasadas. Si no se tienen en cuenta será difícil obtener una estimación razonable del tiempo de observación. •Algoritmo Iterated Greedy Search –A pesar de ser un algoritmo efectivo para el problema de planificación del Flow Shop, caminos de observación alternativos mejores que la solución inicial sólo han sido obtenidas en el caso de escenarios del problema plano. Contents Resumen I 1 Introduction 1 1.1 Earth Observation Satellites 1 1.1.1 Types of Earth observation satellites 2 1.2 Planning of agile satellites constellations 3 1.3 Objectives and scope of the project 4 1.4 Document structure 4 2 Modelling of the problem 5 2.1 Flat problem 5 2.1.1 Position of the vehicle as a function of time 5 2.1.2 Position of the points to observe 5 2.1.3 Logic of analysing the flat problem 6 2.2 Orbital problem 6 2.2.1 Hypothesis of the model 6 2.2.2 Reference frames 6 2.2.3 Position of the vehicle as a function of time 7 2.2.4 Position of the targets to observe 7 2.3 Variables related to the vehicle 8 2.4 Orbital elements 8 2.4.1 Definition 8 2.5 Satellite constellations 9 2.6 Pointing process 10 2.6.1 Static case 10 2.6.2 Dynamic case 12 2.6.3 Optimum rotation time in motion 12 2.7 Kinematics of the problem 13 2.7.1 Pointing to a target 13 2.7.2 Pointing maintenance 14 2.7.3 Examples 15 2.8 Observation restrictions 16 2.9 Groundtracks 19 2.10 Graphical representation 19 2.10.1 Flat problem 20 2.10.2 Orbital problem 21 3 Target Selection Algorithm (TSA) 25 3.1 General search algorithm for a constellation 25 Data input 25 III 6Chapter 2. Modelling of the problem 2.1.3 Logic of analysing the flat problem The results obtained of analysing this problem are not useful because the reality of the Earth observation satellites is an orbital problem. However, several algorithms as the Target Selection Algorithm or optimization ones could be tested analysing this problem. Adapting the algorithms from the flat problem is easy and the trends of the results in both problems could be similar. Furthermore, thanks to the different simplifications involved in this model compared to the orbital problem, it will be possible to better understand the concepts involved in the different equations and constraints involved in both problems. 2.2 Orbital problem The vehicle is an Agile Earth Observation Satellite, which is located in Low Earth Orbit (LEO). The units are km for distance, sfor time and rad for angles. In this model some simplifications have been made: 2.2.1 Hypothesis of the model • Spherical Earth: the geometry of a sphere is enough to approximate the Earth. However, an ellipsoidal model would be more realistic. •Perturbations: in this document no perturbations are taken into account. • Circular orbit: Thank to this hypothesis, an analytical solution of the position of the vehicle (as a function of time) can be obtained. The altitude of the satellite is constant and this has an advantage: the optical sensors of the satellite could be designed to operate in that height. In reality, the orbits of observation satellites are ellipses but they are very similar to circumferences. • Sun-synchronous orbit: It is convenient that the lighting conditions are the same at a given latitude after any number of revolutions, and this can be achievable with this type of orbit. 2.2.2 Reference frames Two reference frames are used in the definition of the trajectory of the vehicle. The first one is the geographical frame (R) , which rotates with the Earth. The xG and yG axis are in the equatorial plane and the zG axis points to the North Pole, the origin OG is the center of the Earth, and the plane OGxGyG is coincident with the Greenwich Meridian. The second reference frame is fixed and does not rotate. It is called Earth-centered inertial frame (I) . The xI and yI axis are in the equatorial plane and the zI axis points to the North Pole. The origin OG is the center of the Earth like the geographical frame. However, the xIaxis always points to the First Point of Aries (). We can change between the reference frames with the following rotation matrix: TG I=  cosβsinβ0 −sinβcosβ0 0 0 1  (2.3) where β=β(t) = w⊕·t+GST0 , w⊕=7.29 ·10−5rad/s , t is the elapsed time since 00:00 UT, in seconds, and GTS0 is the Greenwich sidereal time at the beginning of the day (00 : 00 UT) . It depends on the day and the year, but it is insignificant for the analysis so its value will be considered as 0. Spherical coordinates are used for representing the groundtrack of the satellite in the geographical frame. These coordinates are: φ latitude, λ longitude and r which is the radius. In this particular case r is the radius of Earth, denoted as R⊕. Then, the relations between spherical and Cartesian coordinates are the following: xG=R⊕cosφcosλ yG=R⊕cosφsinλ zG=R⊕sinφ    (2.4) 2.2 Orbital problem 7 Figure 2.1 Spherical coordinates (Adapted from [5]). φ=arctan zG √(xG)2+(yG)2 λ=arctan yG xG    (2.5) It must be considered that there are some ambiguities in the inverse trigonometric functions when transforming from Cartesian coordinates to longitude λ, which can be resolved with the signs of xGand yG. 2.2.3 Position of the vehicle as a function of time Using rotating matrices and defining a set of 6 parameters called "orbital elements", described in Section 2.4, the position of the vehicle, referred to Earth-centered inertial frame, can be obtained. ~rv I(t) =   cosΩcosu−sinΩcosisinu sinΩcosu+cosΩcosisinu sinisinu  (2.6) where r=p=a=R⊕+h,his the altitude of the satellite, and if we have a circular orbit uis: u=u0+n·t n =rµ⊕ (R⊕+h)3(2.7) where µ⊕=398600.442 km3/s2 is the standard gravitational parameter of the Earth and u0 depends on the type of simulation. The inclination i can be obtained from h because the orbit is sun-synchronous, using the following formula obtained from [10]. i=arccos −0.0989 (R⊕ R⊕+h)7/2(2.8) 2.2.4 Position of the targets to observe The target points are contained in the surface of the sphere, and they are referred to the geographical reference frame, so only 2 spherical coordinates are needed (φp,λp) . The position vector of a point is the following one: ~rp G=  R⊕cosφpcosλp R⊕cosφpsinλp R⊕sinφp  (2.9) 8Chapter 2. Modelling of the problem 2.3 Variables related to the vehicle There are a few of parameters which depends on the vehicle itself, which are described below: •Ωmax . This is the maximum angular velocity achievable by the vehicle. Is is supposed to be constant in all directions so the inertia of the satellite is not considered. •tAF . Time which the camera needs to take a photo. •Ψ. This is the maximum angle of pitch and roll. It is considered the same value for both angles. 2.4 Orbital elements 2.4.1 Definition The orbital elements are the parameters required to uniquely identify a specific orbit and, given an arbitrary epoch (initial time t0 ), it is the minimum set of data which allows to determine the position of a body in an orbit at any instant of time. For an arbitrary orbit in space, it is necessary to locate the orbital plane with respect to a reference plane. To achieve that, two parameters are needed ( i , Ω ). Given the orbital plane, and the Earth-centered reference inertial frame OxOyOzO such that the plane OxOyO is coincident with the orbital plane, 3 parameters are needed to determine the orbit. Given a fourth parameter ( θ1 , E , M , or ∆t ), the position of the body can be determined. These six parameters will be the so-called orbital elements. Some conditions must be considered: Figure 2.2 Orbital elements of a 2D orbit (Obtained from [10]). •The reference plane will be the equatorial plane for geocentric orbits •The intersection between the orbital plane and the reference plane determines the line of nodes • The orbit intersects the line of nodes at two points called the orbital nodes. The ascending node (  ) is where the orbiting object moves north through the plane of reference, and the descending node (  ) is where it moves south through the plane. •The node vector~nis an unit vector whose direction is the line of nodes and points to . Then, the six minimum parameters are the following: 2.5 Satellite constellations 9 •Ω : Right ascension of the ascending node (RAAN). It is the angle, measured counter-clockwise, between  and ~n , where  is an unit vector which points to the First Point of Aries, a reference point in celestial coordinate systems. •a , e : Semimajor axis and eccentricity. They are typical parameters of an ellipse and they determine the orbit type and its shape and size. Sometimes they can be replaced by p , T or n , where n is the mean orbital speed, T is the period of the orbit and p is a parameter p=a(1−e2) . Moreover, the eccentricity vector~eis a vector with module ewhich points from origin to the perigee. •ω : Argument of the perigee. (or periapsis). It is the angle, measured in the orbital plane and in the movement direction, between~nand~e •i : Inclination. It is the angle between the reference plane and the orbital plane. Its sense is indicated by ~nand it takes values between 0and πrad. •θ1 : True anomaly. Given an epoch, it determines the position of the body in the orbit. It can be replaced by ∆t , M or E(H) , where ∆t is time increment between t and t0 , M is the mean anomaly and E(H) is the eccentric anomaly. Figure 2.3 Orbital elements of a 3D orbit (Obtained from [10]). In this analysis, some special orbits, where the number of orbital elements is less than six, might appear: • Non-equatorial circular orbit: ω y θ1 are not defined (the line of apsides does not exist). They are replaced by the argument of latitude u=ω+θ1 and it measures the position of the body from the ascending node in movement direction. In this document, non-equatorial circular orbits are analysed with the parameters [Ω,u,i,a,e]. 2.5 Satellite constellations To achieve global coverage so all regions of Earth can be observed, it is necessary to use more than one satellite. 10 Chapter 2. Modelling of the problem A constellation is a set of satellites located in coordinated orbits. In Earth observation constellations, more than one satellites are used in order to offer full or near full coverage of the Earth. The constellations commonly used for Earth observation consists of satellites which share their orbital plane. Examples of this type of constellation are the constellations of cube-sats used by the enterprise Planet, and the satellites Pleiades, which are a set of two agile satellites. The orbital elements of the constellations studied in this document are detailed in Section 3.1. Finally, an example of the Pleiades satellites is shown in Figure 2.4. This constellation consists of two satellites in the same orbital plane which are equidistant in the parameter u. Figure 2.4 Scheme of the organization of the Pleiades constellation. 2.6 Pointing process This section describes the procedure of pointing to a certain target on the surface of the Earth with an agile satellite. First, the variables involved in the rotation of the vehicle are defined. Two types of reference frames are used: the fixed one (R) , defined in the modelling of the flat problem, and the geographical one (G) , which is defined in the modelling of the orbital problem and it is considered as a fixed frame although this is not completely accurate. The characteristics of the observation vehicle ( Ωmax , tAF ,Ψ) are also considered in this section. Then, the superscripts (R) and (G) are going to be omitted as they change if the flat or orbital problem is solved. 2.6.1 Static case In this case the vehicle is stationary in every moment, and the camera is moving from a target a to a target b . 2.6 Pointing process 11 The pointing unit vector of the camera to a generic point pis defined as: ~ Cp(t) = ~rp−~rv(t) |~rp−~rv(t)|(2.10) Where ~rp is the position of the target point respect the reference frame expressed in Cartesian coordinates, and ~rv(t)is the position of the satellite for a given time t. Once the position of the 2 points involved in the turn and the position of the vehicle in a time t are known, the rotation which minimizes the time is a turn around an axis perpendicular to both pointing vectors. This axis is called as Euler Axis [10], and can be calculated as: ~e0= ~ C0 a∧~ C0 b |~ C0 a∧~ C0 b| (2.11) where: ~ C0 a=~ Ca(t0)~ C0 b=~ Cb(t0)(2.12) are the pointing vectors to each point in a time t0 . The vector ~e0 is an unit vector, and the angular rotation speed in the pointing process between aand bcan be obtained from this vector: ~ ω0=Ω~e0(2.13) where Ω is the rotation speed of the satellite. The rotation angle performed by the camera in the pointing process between the two targets in the static case can be obtained as: θ(t0) = θ0=arccos[~ Ca 0·~ Cb 0] = arccos[~ Ca(t0)·~ Cb(t0)] (2.14) If the angle of rotation and the maximum rotation speed of the satellite Ωmax are known, the minimum rotation time is: t∗ 0=θ0 Ωmax (2.15) Figure 2.5 Static pointing process between two targets (Adapted from [5]). 12 Chapter 2. Modelling of the problem 2.6.2 Dynamic case In this case, the vehicle is in motion and the camera is pointing from an initial target a at a time t0 , to the target b during an unknown time t . The pointing vectors are defined differently than the static case. The initial pointing vector to the target ais: ~ Ca=~ra−~rv(t0) |~ra−~rv(t0)|=~ Ca 0(2.16) The pointing vector to the target bafter a time tfrom t0is: ~ Cb=~ Cb(t) = ~rb−~rv(t0+t) |~ra−~rb(t0+t)|(2.17) This vector is a function exclusively of the pointing time t . The angle of rotation performed by the camera is: θ(t) = arccos[~ Ca·~ Cb(t)] (2.18) This angle is exclusively determined by the pointing time t , since the positions of the targets a and b are known and the initial pointing vector is known given the starting time t0. The Euler rotation axis as a function of tcan be obtained as: ~e=~ Ca∧~ Cb(t) |~ Ca∧~ Cb(t)|(2.19) For any value of t , we will have a specific value to ~ Cb(t) and therefore a value of θ(t) . Knowing a value of the rotation angle, the associated rotation speed can be obtained as: Ω(t) = θ(t) (t+t0)−t0 =θ(t) t(2.20) However, the rotation speed is limited by Ωmax: Ω(t)∈[0,Ωmax](2.21) 2.6.3 Optimum rotation time in motion As exposed in the document by Meléndez [ 5 ] there is a minimum value of the pointing time t which optimizes the rotation. Thus, 2 different methods to obtain this value of t are described in that document. However, in this project we use the "Method 2" which is faster for tolerance greater than 10−12 , which is sufficient for the intended accuracy of the program. The convergence of this method has also been demonstrated in the project of Meléndez. To obtain the optimum rotation time in motion t∗, the following iterative method must be performed: •For all n≥0we calculate iteratively: t∗ n+1=θ(t∗ n) Ωmax =arccos[~ Ca·~ Cb(t∗ n)] Ωmax (2.22) when the initial value for n=0 correspond to the optimum rotation time of the static case t∗ 0 . The stop condition will be: pn+1≤tol (2.23) where tol is the desired tolerance, and pn+1is: pn+1=|t∗ n+1−t∗ n|(2.24) In this project a tolerance of tol =10−6is taken. 2.7 Kinematics of the problem 13 Figure 2.6 Dynamic pointing process between two targets (Adapted from [5]). 2.7 Kinematics of the problem Once the optimum rotation time in motion of the vehicle is known, the pointing process between two known target points can be defined. The following data are required to define the pointing process: •~ Ca: Pointing vector to the initial point aat an initial time t0: ~ Ca=~ Ca 0=~ra−~rv(t0) |~ra−~rv(t0)|(2.25) •~ Cb : Pointing vector to the target point b after the optimal rotation time t∗ has elapsed from the initial time t0. ~ Cb=~ Cb(t∗) = ~rb−~rv(t0+t∗) |~rb−~rv(t0+t∗)|(2.26) •θ: Rotation angle: θ=arccos[~ Ca·~ Cb](2.27) •~e: Euler rotation axis: ~e=~ Ca∧~ Cb |~ Ca∧~ Cb|(2.28) The process is divided into two main movements: The first movement consists in the rotation of the vehicle to point the camera from the starting target a to the target point b ; and the second movement consists in keeping the camera directed at point b during the aiming pointing time tAF , which the camera needs to stabilize and take the photography. 2.7.1 Pointing to a target In the pointing of the camera from a target point a to a target b , the angular rotation speed vector is defined as: ~ ω=Ωmax ·~e(2.29) 14 Chapter 2. Modelling of the problem The differential equation which governs the motion of the pointing vector is: d~ C(t) dt =d~ CM(t) dt +~ ω∧~ C(2.30) where ~ C is the pointing vector of the camera of the vehicle in the reference frame used in the problem (flat problem or orbital problem), ~ CM is the pointing vector with respect to a reference frame centered in the vehicle. This reference frame rotates with the vehicle with the angular rotation speed defined in (2.29) so the pointing vector of the camera is always contained in the xaxis of the frame. Thus, the vector ~ CM remains constant in time with respect to that frame, so is true that d~ CM(t) dt =0 and we have: d~ C(t) dt =~ ω∧~ C(2.31) Since the initial value of the pointing vector is known: ~ C(t=t0) = ~ Ca(2.32) we can solve the differential equation (2.31) using a numerical method and obtain the evolution of the camera vector in each instant ~ C=~ C(t)from the point ato b, where: ~ C(t=t0+t∗) = ~ Cb(2.33) 2.7.2 Pointing maintenance This is the process which keeps the camera pointing to the target b during the time tAF . The differential equation is the same which we used in the previous section: d~ Cb(t) dt =~ ω∧~ Cb(t)(2.34) However, in this case ~ ω is not given by the optimum rotation of the vehicle, so a new expression is deduced to obtain ~ ω. Once ~rb , position vector of point b , is known, the pointing vector to target b in a certain instant of time t is defined as: ~ Cb(t) = ~rb−~rv(t) |~rb−~rv(t)|(2.35) Although the vector ~ Cb is known at every time, the differential equation is solved, so external disturbances can be taken into account in the future. The expression of the derivative of ~ Cb can be obtained, since all the variables involved in this expression are already known. d~ Cb(t) dt desired =~ Cb(t)0=1 |~rb−~rv(t)| d dt [~rb−~rv(t)]+ d dt 1 |~rb−~rv(t)|[~rb−~rv(t)] (2.36) The development of the equation (2.36) is detailed in the annex of the document. Since ~ Cbis an unit vector, the angular rotation speed is defined as: ~ ω(t) = ~ Cb(t)∧~ Cb(t)0(2.37) 2.7 Kinematics of the problem 15 To show that the angular speed of (2.37) guarantees the pointing maintenance to target b , this expression is replaced in the differential equation (2.34). d~ Cb(t) dt =~ ω∧~ Cb(t) = [~ Cb(t)∧~ Cb(t)0]∧~ Cb(t) = =−~ Cb(t)∧[~ Cb(t)∧~ Cb(t)0] = − ~ Cb(t)~ Cb(t)0 ~ Cb(t)·~ Cb(t)~ Cb(t)·~ Cb(t)0 =(2.38) =~ Cb(t)0[~ Cb(t)·~ Cb(t)]−~ Cb(t) [~ Cb(t)·~ Cb(t)0] where ~ Cb(t)·~ Cb(t) = |~ Cb|2=1, and 1 2·d dt |~ Cb(t)|2=~ Cb(t)·~ Cb(t)0=0since |~ Cb|=1. Finally we obtain: d~ Cb(t) dt =~ Cb(t)0=d~ Cb(t) dt desired (2.39) Thus, with the angular rotation speed defined in (2.37) and an initial condition of the pointing vector ~ Cb=~ Cb(t0+t∗) , the differential equation (2.34) can be solved using a numerical method, obtaining the evolution of the pointing vector in the pointing maintenance process. The equation must be solved until the time t=t0+t∗+tAF is reached, where: ~ Cb=~ Cb(t0+t∗+tAF )(2.40) 2.7.3 Examples In Figures 2.7 and 2.8 two examples of the pointing process is represented for the flat and orbital problem. Figure 2.7 Pointing process in flat problem (Adapted from [5]). 22 Chapter 2. Modelling of the problem Hence, the parameter νs(t)can be calculated as: xl(t)2+yl(t)2+zl(t)2=|~rl(t)|2=R2 ⊕(2.58) Where: |~rl(t)|2=|~rv(t)+νs(t)~ C(t)|2=|~rv(t)|2+νs(t)2|~ C(t)|2+2νs(t) [~rv(t)·~ C(t)] (2.59) Thus, combining (2.58) and (2.59) we have: |~rv(t)|2+νs(t)2|~ C(t)|2+2νs(t) [~rv(t)·~ C(t)]−R2 ⊕=0(2.60) which is a quadratic equation with νs(t)as variable. Then, this variable is going to be named as τ(t). a(t)τ(t)2+b(t)τ(t)+ c(t) = 0(2.61) where: a(t) = |~ C(t)|2;b(t) = 2~rv(t)·~ C(t);c(t) = |~rv(t)|2−R2 ⊕(2.62) Now the parameter τ(t)can be obtained: −b(t)±pb(t)2−4a(t)c(t) 2a(t)=τ1(t) τ2(t)(2.63) When both solutions are positives, the solution which has less absolute value is taken since the other solution corresponds to the intersection of the line with the opposite side of the Earth. However, when one solution is positive and the other is negative, the positive one should be taken as the pointing vector points to the Earth. Once τ(t)is known, replacing in (2.53) we obtain: ~rt(t) = ~rv(t)+τ(t)~ C(t) =   xv(t)+τ(t)C1(t) yv(t)+τ(t)C2(t) zv(t)+τ(t)C3(t)  (2.64) On the other hand, as the Earth is considered as a sphere of radius R⊕ , the groundtrack of the position of the satellite is: ~rtv(t) = R⊕ ~rv(t) |~rv(t)|=R⊕ pxv(t)2+yv(t)2+zv(t)2  xv(t) yv(t) zv(t)  (2.65) We can observe an example of graphical representation of the camera and satellite groundtrack in Figure 2.14 for the case of orbital problem. In this figure it is verified that the groundtrack of the camera are not straight lines. 2.10 Graphical representation 23 Figure 2.14 Graphical representation of vehicle and camera in orbital problem. 3 Target Selection Algorithm (TSA) In this section we explain the target selection algorithm of the constellation which is used in the simulations. This algorithm also allows to obtain an initial sequence for the optimization algorithm. 3.1 General search algorithm for a constellation The general algorithm, given a set of target points in a certain scenario, selects the observation path for all the satellites of the constellation. In Figure 3.5 an schematic flow chart of the algorithm is shown. In order to start the algorithm, some parameters are assumed known. Then the program will go into a loop, and it will evaluate the situation of each satellite step by step, until all the targets have been observed or a certain time has elapsed in the simulation. The various modules which constitutes this algorithm are detailed hereunder. Data input The algorithm receives the following data as input: • Parameters of each vehicle of the constellation: for example, focusing time tAF , maximum rotation speed Ωmax , maximum deviation of nadir Ψ , etc. Some of these parameters depend on the type of problem studied, either the flat problem or the orbital one. –Flat problem: Height and speed of the satellite must be provided: hi,Vi –Orbital problem: The following data of each satellite is required by the algorithm: Ωi,u0i,hi. The constellations studied in this documents consists of sun-synchronous satellites in the same orbital plane which are equidistant in the parameter u . Thus, the parameters ( Ω , h ) are the same for all satellites, and u0iis a vector with the following structure: [u01,u02,...u0i,...u0n] Where nis the number os satellites of the constellation, and u0iis: u0i=u01 +2π n(i−1) • Set of target points: It consists in a matrix rp with dimensions nx2 which will contain the position of the target points. First two columns of rp indicate the Cartesian coordinates in the plane z=0 ( rxi , ryi ) for the flat problem, while these columns represent the spherical coordinates of latitude and longitude on the sphere of Earth ( rλi , rφi ) in the case of orbital problem. This matrix has the following structure: rp|f lat =     rx1ry1 rx2ry2 . . .. . . rxn ryn      rp|orbit =     rλ1rφ1 rλ2rφ2 . . .. . . rλnrφn      (3.1) 25 26 Chapter 3. Target Selection Algorithm (TSA) Main loop In this module, the algorithm will check if all the target points has been observed or if the time in the simulation has reached a certain value. If the condition is not met, the loop is maintained; if the condition is met, the loop is terminated. Constellation loop In this loop, it is checked that each satellite of the constellation has accomplished the action that correspond to the current step of the algorithm. Once every satellite has been checked the loop is ended. Satellite function Is the main function and is in charge of generating the observation path for each satellite. It will be detailed later. This function must receive the data of the vehicle i and the set of target points as input, while it return the value of pointing vector ~ Ci associated with the satellite in next step of the algorithm, as well as information about the targets: when they are observed, which satellite observes each target, etc. Results and graphical representation When the main loop ends, the algorithm shows results of the simulation as observation time of the scenario, order in which the targets of rp are observed, etc. A graphical representation of the observation path is also obtained if necessary. 3.2 Satellite function: Target selection algorithm This is the function which generates the path of the groundtrack of the camera in each step of the simulation. Then, this part of the algorithm is detailed in Figures 3.6 and 3.7 for the orbital problem. This algorithm is based on the evaluation of the state of each satellite. This is represented in the state variable which may take various values. For each satellite and simulation step, state is stored in a matrix with dimensions sxn where s is the number of steps of the algorithm and n is the number of satellites of the constellation. This matrix is both input and output of the function and has the following structure:           state(1,1)state(1,2)··· state(1,j)··· state(1,n) state(2,1)state(2,2)··· state(2,j)··· state(2,n) . . .. . .. . .. . .. . .. . . state(i,1)state(i,2)··· state(i,j)··· state(i,n) . . .. . .. . .. . .. . .. . . state(s,1)state(s,2)··· state(s,j)··· state(s,n)           (3.2) The state variable may take the following values: •state =0 : the satellite searches for possible observable targets in the current step and chooses the one with the highest priority. To do this the search of targets algorithm is executed, which is composed by three subfunctions that we will explain later: "Selection of points under coverage", "Sort by distance or by topt algorithm" and "Point check algorithm". Once a target has been selected by the algorithm, state changes to value 1 and the chosen target point is stored in the memory. Furthermore, the optimum rotation axis and the needed time for turning are calculated. In case of orbital problem, if no target points are encountered, the waiting algorithm is run. For flat problem, state takes the value 3. •state =1 : The satellite performs the pointing to a target process described in Section 2.7.1. Once this process ends and the camera aims to the target, the state variable takes the value 2. •state =2 : The satellite performs the camera maintenance process on a point described in Section 2.7.2. Once a time TAF has elapsed and the photography has been taken, the state variable takes the value 0 . 3.2 Satellite function: Target selection algorithm 27 •state =3 : This state is only valid for flat problem, and it is part of the waiting algorithm. Its similar to state =0 but target points are searched as if the satellite has double Rci than it has reality. Hence, the selected target is not under coverage at that time, but the camera can prepare the observation of that target. These targets are called as "far targets". If a target is selected, the state variable takes the value 4, if not, the waiting algorithm is run. •state =4 : It is similar to state =1 but the satellite points to a far target. If the point enters under coverage before the process has finished, the state variable takes the value 1 . If the pointing process has finished, the state variable takes the value 5. •state =5 : It is similar to state =2 but the satellite points to a far target. If the point enters under coverage before the process has finished, the state variable takes the value 2 . If the pointing process has finished, the state variable takes the value 0. Next, 4 subfunctions related to state =0are explained. 3.2.1 Points under coverage selection function This function selects the points within the satellite coverage circumference in the current step of the simulation. This points are stored in a vector which is input of the next subfunctions. Thus, a number of targets Ncob potentially observable will be chosen. Figure 3.1 Targets selected after points under coverage selection function for flat problem. An example of the operation of this algorithm is shown in Figure 3.1, where a case of the flat problem is represented because it is more understandable than the orbital problem. 3.2.2 Points selection algorithm Two points selection algorithm are used in this document: "Sort by distance" and "Sort by minimum rotation time". These two methods are used in this analysis, however other types of selection algorithms can be used, such as ordering the points based on a priority function, or simply have them sorted randomly. 28 Chapter 3. Target Selection Algorithm (TSA) Sort by distance This function sorts all the target points under coverage based on the distance from the target where the camera is pointing. This points are output of the points under coverage selection function, they are sorted and stored in a vector so the closest points the camera groundtrack has index i=1 and the farthest point has index i=Ncob. One example of selecting the closest point is shown in Figure 3.2. Figure 3.2 Target selected after points selection algorithm for flat problem. Sort by rotation time This function sorts all the target points under coverage based on minimum rotation time topt from the target where the camera is pointing. 3.2.3 Point check algorithm This algorithm checks that all target points obtained after the points selection algorithm are, indeed, observable. Although a point is under coverage at a time t0 , after the rotation time t0+t∗+tAF this point might no longer be covered, so this function is necessary to obtain true observable targets. For a target point to be observable, the nadir deviation restrictions described in Section 2.10 must be met both at the initial time t0and after the rotation time has elapsed: φp=φ(t0) = arccos[~ Cp·~n] = arccos[~ C(t0)·~n(t0)] (3.3) φf p=φ(t0+t∗+tAF ) = arccos[~ Cf p·~nf] = arccos[~ C(t0+t∗+tAF )·~n(t0+t∗+tAF )] φp≤Ψand φf p≤Ψ(3.4) Thus, the condition is applied to the first target of the sorted by distance points vector. If this target does not fulfil the condition, it is applied to the next point and so on until an observable point is obtained. An example of the operation of this function is shown in Figure 3.3. Finally, that checked target is stored in the memory to be used by the rest of modules of the general algorithm. It is stored in a matrix in a way which allows to distinguish between targets chosen by different satellites of the constellation. If a chosen point is a vector ~pe, this matriz has the following structure: ~pe1~pe2··· ~pei··· ~pen(3.5) 3.2 Satellite function: Target selection algorithm 29 Figure 3.3 Target selected after point check algorithm for flat problem. Where nis the number of satellites of the constellation. 3.2.4 Waiting algorithm This algorithm is executed when state =0 and no observable points is found. There are two posible actions to take depending of the problem: • Flat problem. In this case, the pointing vector ~ C is maintained as constant for the next step of the algorithm. ~ C(t+∆t) = ~ C(t) • Orbital problem. In this problem, keeping the pointing vector ~ C constant might be problematic. The selection algorithm works choosing targets based on the camera groundtrack; therefore, since the target surface is a sphere, it could happen that the satellite points towards space, where there is no groundtrack. This problem is illustrated in Figure 3.4. To avoid errors, if no observable points are found, the camera points to the centre of Earth in the direction of nadir ~n . To do this, the process of pointing to a point described in Section 2.7.1 is performed, having: ~ Cb=~n Where ~n is the nadir. It would be good to prepare the camera to observe the next point before it is under coverage. This is possible in the case of the flat problem where all targets are observable at any pass. However, in the orbital problem, it is difficult to select the targets and prepare the camera as not all points are observable at any pass. The only targets which can be selected are the ones under coverage and have been checked with "Point check algorithm". 30 Chapter 3. Target Selection Algorithm (TSA) Figure 3.4 Problem of the waiting algorithm in the orbital problem. 3.2 Satellite function: Target selection algorithm 31 Figure 3.5 General algorithm of the constellation. 38 Chapter 3. Target Selection Algorithm (TSA) Figure 3.11 Iterated Greedy Search algorithm flow chart2. 3.3 Optimizations Algorithms: Iterated Greedy Search (IGS) 39 Figure 3.12 Local Search Algorithm flow chart. 4 Results and analysis 4.1 Target points generation In this section we explain how the cloud of target points is generated in every scenario of the simulations. The objective points are located on the observation surface, either the plane z=0 for the flat problem, or the terrestrial sphere in the orbital problem. We proceed to detail each of the scenarios to analyse that have been created: 4.1.1 Flat Problem In this case, the target points are generated randomly in a rectangle with length L=200 UD , and width w=Rci UD , where Rci is the instrumental coverage radius of the satellite. In this way all targets are observable Then, random values of xp∈[0,200]UD e yp∈[−Rci/2,Rci/2]UD will be generated following an uniform distribution. An example of a scenario for the flat problem is the one shown in Figure 4.1. Figure 4.1 Example of escenario in flat problem. 4.1.2 Orbital Problem In the case of the orbital problem, two types of scenarios are analysed in this document: a global scenario of the complete Earth and a local scenario on the Iberian Peninsula. 41 42 Chapter 4. Results and analysis Global scenario For this type of scenario, points are generated on the entire earth at random. Random length values λp∈[−180o,180o] are generated following an uniform distribution, while latitude values φp∈[−90o,90o] are generated following a normal distribution with mean µ=0o and standard deviation σ=3 8i , where i is the inclination of satellites orbits. Hence it is achieved that there are more target points at the equator than at the poles. Some clarifications should be made about the way the target points are generated: • If a value generated for φp , is greater than 90o or lower than −90o , the value is eliminated and a new one is obtained. • The target points do not follow a uniform distribution because the distances are smaller next to the poles than next to the equator in a cylindrical equidistant projection. Moreover few people live next to the poles so the fact that there are more targets next to the equator makes sense. An example of this type of scenario is shown in Figure 4.2. Figure 4.2 Example of global scenario in orbital problem. Local scenario For this case, the points are generated in a uniform way in a spherical rectangle located on the Iberian Peninsula. Random values of λp∈[−9.2o,3.2o] y φp∈[−35.83o,43.71o] following a uniform distribution. In this case they are uniform in latitude because the zone is small and the distances do not vary too much as in the global scenario. They are used as coordinates of each target point. An example of this type of scenario is shown in Figure 4.3. 4.2 Repeated observation targets In reality, there are target points that must be observed several times, since during one pass there may be clouds over that area, or the camera could have obtained a poor quality photograph. The aim of this document is to observe the sensitivity of the scenario observation time to the presence of several target points that must be observed during several consecutive passes. To determine how many times each objective point must be observed, random integer numbers greater than 1, which are assigned to each objective point, are generated. A number is generated following a normal distribution with mean µ=1 and standard deviation σ=1 . If this number is lower than 1, if will be 4.2 Repeated observation targets 43 Figure 4.3 Example of local scenario in orbital problem. eliminated and a new number will be created. If the number is greater than 1, it will be rounded so that it is integer. This number is denoted as nmax Using this method, it will be a great probability that nmax =1 , so only a few target points must be observed in more than one pass. At the beginning of each pass, the list of target points is checked: If the number of times n a target has been seen is less than nmax , that target can be observed once during this pass. On the other hand, if n=nmax the selection point algorithm will not take that target into account. 44 Chapter 4. Results and analysis 4.3 Global scenario simulations If agile satellites are equipped with an infrared camera, it is possible to obtain images of points from any latitude. Therefore, it might be interesting to study the behaviour of the constellations in case the target points are distributed throughout the globe. The behaviour of a constellation of satellites, located in the same orbital plane and equispaced with respect to the u0 variable, is simulated. This way of organizing the constellation is typical of Earth observation satellites such as the Pleiades satellites or the constellations of the enterprise Planet. For the simulations, the two scenarios of the Figure 4.4 are generated, one with 300 points and the other with 600. It is necessary to comment that the map which is reflected in these images is not exact and is merely indicative. The targets are generated in such a way that the majority are close to the equator and in less concentration near the pole, and they are uniformly distributed in λ. Figure 4.4 Global simulation scenarios. Once the scenarios are generated, the behaviour of the constellation is simulated during 24 hours and the number of targets observed is obtained. The Target Selection Algorithm described in Chapter 3 is used in these simulations. The number of target points is not high enough to saturate the satellites, in other words, when a satellite passes through an area, all the points which are within the coverage can be observed without having to pass again. However, the number of target points is high enough to say that if a constellation observes all of them, then it covers the entire earth. Furthermore, since the satellite is not saturated, in these simulation it does not seem so important to know at all time the evolution of the pointing vector ~ C. Parameters of the simulations We proceed to define the fixed parameters of the simulation, common to all satellites of the constellation: •In these simulations the RAAN does not affect the results and is given a null value: Ω=0. •The epoch does not affect the results in these simulations: GST0=0 •Maximum angular speed of the satellite: Ωmax =1.5o/s •Time needed to take a photography: Ta f =10 s As the target are generated following an uniform distribution in λ and the satellite must cover the entire Earth, the point where the satellite starts in t0 does not matter. This is why the RAAN and the epoch does not affect the results. However, these parameters will be required to obtain a real observation path. The parameters that are variable between simulations are the altitude h , the maximum angle of rotation of the camera Ψand the number of satellites (between 1 and 5). The altitudes studied in this document are: 4.3 Global scenario simulations 45 • Altitude of a satellite constellation of the enterprise Planet (Data obtained from [ 1 ]): h=487.86 Km •Altitude of the Pleiades satellite constellation (Data obtained from [4]): h=694 Km •Average altitude between both constellations: h=590 Km Furthermore, the parameter Ψwill take the following values: •Ψ=15o •Ψ=30o •Ψ=45o Figure 4.5 Example of simulation with two satellites. The black groundtrack corresponds to the first satellite and the grey groundtrack corresponds to the second one. Varying Ψ has some interest since, although the maximum value of the angle of rotation for a certain satellite has a value Ψsat that depends on its technical characteristics, if the angle is forced to be lower, photos will be taken of higher quality as there is less deviation from the nadir. In the same way, by reducing the altitude of the constellation, photos with more resolution can be obtained as the camera is closer to the ground. The parameter u0 must be taken in a way that the satellites are equispaced in the orbital plane. An explanation of how u0 changes if the number of satellites increases is shown in Section 3.1. In these simulations u01 =0o. In short, the orbital elements of the satellites at t0are: •The RAAN is common for all satellites: Ω=0 • Inclination is common for all satellites and is calculated using the altitude as the satellites are sunsynchronous. It is calculated using the following formula described in Section 2.2.3. i=arccos −0.0989 (R⊕ R⊕+h)7/2(4.1) •The orbits are circular: a=h,e=0 •The parameter u0changes between satellites and is calculated using the following formula: u0i=u01 +2π n(i−1) Where nis the total number of satellites of the constellation. 46 Chapter 4. Results and analysis Results The results of simulations are shown in Table 4.1. In Figure 4.5 the path of a constellation of two satellites with altitude h=694 Km is shown. This case is a simulation of the Pleiades satellite constellation. It is observed that the groundtracks of both satellites surround the Earth in 24 h , and for the value of Ψ=45o that these satellites have, all the points can be observed, so the Earth is perfectly covered. Therefore, it seems that these simulations imitate the behaviour of a real constellation [ 4 ] and they allow us to obtain some conclusions. To better understand the trends in the results, h−Ψ plots are made for various constellation examples and for both scenarios. These plots are shown in Figure 4.6. Figure 4.6 h−Ψplots of various constellations in global scenario. The h−Ψ plots have interest since they allow knowing when a constellation can cover the earth in one day, and give an idea of the quality of the images: the closer to the origin (0,0) , the higher quality the photographs obtained will have. The variation of the points observed with the parameter Ψ is much greater than with the altitude. On the other hand, the plots are similar in shape for the 300 and 600 target points scenarios. This means that if the satellite is not saturated and if the constellation is able to cover the entire earth all points will be observed. A simple way to optimize the organization of these satellites is to find a constellation that observes all the points, and whose parameters ( Ψ , h ) are the closest to ( 0 , 0 ). This ensures that each satellite has a "band" where it will obtain the images, and where the nadir deviation is the smallest possible, assuring a good quality of the photographs. 4.3 Global scenario simulations 47 Table 4.1 Results of the simulations in the global scenario. Number of satellites h(Km)12345 487.86 57 110 157 196 209 590 73 138 190 218 258 694 79 146 199 240 280 Ψ=15o300 points scenario Number of satellites h(Km)12345 487.86 135 210 273 297 300 590 151 254 290 300 300 694 168 261 298 300 300 Ψ=30o300 points scenario Number of satellites h(Km)12345 487.86 211 290 300 300 300 590 238 300 300 300 300 694 254 300 300 300 300 Ψ=45o300 points scenario Number of satellites h(Km)12345 487.86 116 208 289 356 430 590 155 270 372 448 514 694 179 309 438 506 551 Ψ=15o600 points scenario Number of satellites h(Km)12345 487.86 268 426 543 586 600 590 313 489 573 599 600 694 330 524 594 600 600 Ψ=30o600 points scenario Number of satellites h(Km)12345 487.86 423 575 600 600 600 590 466 598 600 600 600 694 507 600 600 600 600 Ψ=45o600 points scenario On the other hand, it would also be interesting to study the targets observed as a function of the number of satellites available in the constellation. This trend is shown in the plots of the Figure 4.7. Figure 4.7 Targets observed as a function of the number of satellites for a given altitude. For small values of Ψ or low number of satellites, it is not possible to obtain full coverage. In this case for 54 Chapter 4. Results and analysis • As h increases, the speed of the satellite when flying over the zone is lower and it spends more time over the area. This parameter also modifies the coverage radius but in a smaller way. Table 4.4 Number of points per satellite pass in each simulation using "Sort by distance". Number of satellites h(Km)12345 487.86 – 4.17 – 1.92 2.63 590 – 4.55 – 4.17 4.17 694 5 5 4.17 5 5 Ψ=15o50 points scenario Number of satellites h(Km)12345 487.86 5 4.55 5.56 6.25 5 590 6.25 5.56 8.33 7.14 7.14 694 7.14 6.25 8.33 7.14 8.33 Ψ=30o50 points scenario Number of satellites h(Km)1 2345 487.86 7.14 8.33 7.14 7.14 8.33 590 6,25 10 10 10 10 694 16,67 12,5 10 10 10 Ψ=45o50 points scenario Number of satellites h(Km)12345 487.86 – 4.17 – 2.44 2.86 590 – 5.88 – 3.03 5 694 5.88 4.35 5.56 4.18 5.88 Ψ=15o100 points scenario Number of satellites h(Km)12345 487.86 6.25 5.88 6.67 5.56 5.56 590 6.67 6.67 8.33 7.69 7.14 694 9.09 8.33 10 9.09 8.33 Ψ=30o100 points scenario Number of satellites h(Km)12345 487.86 9.09 12.5 10 10 10 590 10 11.1 11.1 11.1 10 694 14.3 14.3 12.5 16.7 16.7 Ψ=45o100 points scenario A compromise between quality and observation time must be valued by the operator of the satellite constellation. Furthermore, a study of how much the variation of parameters h and Ψ affects the quality of the photo can be quite interesting in the valuation, because it might be that these problems are solved by improving the technical capacity of the satellite and the camera with low cost. 4.4.2 Point selection by minimum rotation time In this section, the "Sort by topt " algorithm is used to create an observation path. The trends when varying the simulation parameters are similar to those obtained for the "Sort by distance" algorithm. The results of the simulations are shown in table 4.5. It is interesting to compare the results of both methods to check if one of them is better than the other. Hence, the results of the simulations obtained are represented for both algorithms in Figure 4.14. The values of Ψ=15o for h=487.86 , 590 Km are not represented because this situations have not full coverage of the Earth and they are problematic. In the graphs of the comparison it is observed that for certain occasions the "Sort by distance" algorithm is better, while in others the "Sort by topt " algorithm is the most appropriate. However, there are cases where both algorithms generate the same observation path. This fact is quite surprising since in the experiments 4.4 Local scenario simulations 55 Figure 4.14 Comparison of the observation time tobs for both point selection methods. carried out by Meléndez [ 5 ], it seems that the selection algorithm by minimum rotation time was better for any scenario, as illustrated in Figure 4.15. However, this experiment consisted of observing the greatest number of points during one pass, and it is a different situation from the problem of this document. It is also interesting that the lower Ψ is, as the problem is more restricted, both algorithms find the same path in most cases. Therefore, as the optimal path depends on the scenario, the altitude of the satellites and the parameter Ψ , it seems that the best option would be to obtain two alternative paths with both methods and study which of them is the most appropriate. In this way, calculating another alternative path, the observation time of the scenario could be reduced by several days. Finally, the number of average points observed per pass is shown in the table 4.6. The trends and values of this result remain similar to the simulations with the algorithm "Sort by distance". 56 Chapter 4. Results and analysis Table 4.5 Observarion time of each scenario tobs in hours using "Sort by topt ". Number of satellites h(Km)1 2 3 4 5 487.86 – 192 – 168 121 590 – 192 – 96.1 72.4 694 336 168 144 72.4 71.7 Ψ=15o50 points scenario Number of satellites h(Km)12345 487.86 431 168 144 72.0 48.2 590 240 96.3 48.3 24.3 24.7 694 216 95.7 71.9 47.9 24.5 Ψ=30o50 points scenario Number of satellites h(Km)12345 487.86 168 71.6 47.8 48.1 47.9 590 191 72,5 47.8 24.2 24.2 694 95.5 47.8 24.2 23.9 24.1 Ψ=45o50 points scenario Number of satellites h(Km)12345 487.86 – 407 – 383 263 590 – 311 – 311 144 694 671 407 215 215 120 Ψ=15o100 points scenario Number of satellites h(Km)1 2 3 4 5 487.86 647 383 216 144 72.4 590 288 191 96.6 96.1 48.6 694 288 144 72.4 72.0 71.8 Ψ=30o100 points scenario Number of satellites h(Km)12345 487.86 288 120 72.4 71.8 47.9 590 216 95.6 48.3 25.0 47.7 694 120 71.6 47.8 24.8 24.1 Ψ=45o100 points scenario Figure 4.15 Number of targets observed as a function of the number of points in the scenario for Meléndez’s experiment. Figure obtained from [5]. 4.4 Local scenario simulations 57 Table 4.6 Number of points per satellite pass in each simulation using "Sort by topt ". Number of satellites h(Km)12345 487.86 – 4.17 – 2.63 2.63 590 – 4.55 – 4.17 4.17 694 5 5 4.17 4.17 5 Ψ=15o50 points scenario Number of satellites h(Km)12345 487.86 3.85 4.55 3.85 5 5 590 6,25 7.14 6.25 7.14 7.14 694 5 6.25 8.33 10 7.14 Ψ=30o50 points scenario Number of satellites h(Km)12345 487.86 8.33 8.33 8.33 7.14 5.56 590 6.25 6.25 8.33 8.33 8.33 694 10 10 10 10 8.33 Ψ=45o50 points scenario Number of satellites h(Km)12345 487.86 – 4.17 – 2.44 2.86 590 – 5.88 – 3.03 5 694 5.27 4.35 5.56 4.18 5.88 Ψ=15o100 points scenario Number of satellites h(Km)12345 487.86 5.26 4.17 5.26 5.56 6.67 590 9.09 7.69 9.09 7.14 7.14 694 9.09 10 11.1 9.09 8.33 Ψ=30o100 points scenario Number of satellites h(Km)12345 487.86 9.09 9.09 10 10 11.1 590 11.1 10 14.3 12.5 10 694 14.3 14.3 12.5 14.3 16.7 Ψ=45o100 points scenario 4.4.3 Simulations with repeated observation targets In this part the results obtained in Section 4.4.1 are compared with simulations where some points have to be observed several times. This pretends to recreate reality, where there are sometimes clouds over the target point, or the camera has had an error and has not been able to take the picture correctly. The 50 targets scenario in Figure 4.10 is analysed and the results of these simulations are shown in table 4.7. Furthermore, an example of an observation path is shown in Figure 4.17, where it can be seen that some objective points are observed several times. 58 Chapter 4. Results and analysis The comparison between simulations is shown in Figure 4.16. The result is quite interesting because the observation time increases significantly on the stage since some points must be observed several times. This means that if these random phenomena such as the presence of clouds are not taken into account, the actual observation time of the points can be much longer than expected. However, the trends remain similar in both cases. Figure 4.16 Comparison between simulations with targets observed after several attempts, and simulations with target points observed at once. 4.4 Local scenario simulations 59 Figure 4.17 Observation path in a simulation with repeated observation targets. 60 Chapter 4. Results and analysis Table 4.7 Observation time tobs in hours and number of points per pass in the 50 targets scenario. Number of satellites h(Km)12345 487.86 – 671 – 264 216 590 – 432 – 192 168 694 886 479 216 192 192 Ψ=15oObservation time Number of satellites h(Km)1 2 3 4 5 487.86 647 264 216 144 96.6 590 407 216 144 95.7 95.9 694 408 192 121 72.8 72.1 Ψ=30oObservation time Number of satellites h(Km)12345 487.86 312 168 120 72.4 72.1 590 240 120 72.5 48.7 48.0 694 168 96.3 72.4 47.8 25.1 Ψ=45oObservation time Number of satellites h(Km)12345 487.86 – 1.32 – 1.61 1.72 590 – 2.08 – 2.27 1.79 694 2.08 1.85 2.63 2.27 1.79 Ψ=15oTargets per pass Number of satellites h(Km)12345 487.86 2.63 2.94 2.63 2.50 2.50 590 3.57 3.33 3.13 3.85 2.94 694 3.33 3.33 3.57 3.85 3.85 Ψ=30oTargets per pass Number of satellites h(Km)12345 487.86 3.85 3.57 3.33 4.17 3.57 590 5 4.55 5 4.55 4.55 694 5.56 5 4.17 5.56 5.56 Ψ=45oTargets per pass 4.5 Iterated Greedy Search algorithm results 61 4.5 Iterated Greedy Search algorithm results Finally, the iterated greedy search algorithm has been adapted to the satellite observation scheduling problem. This algorithm takes as its initial solution the observation sequence obtained by applying the "Sort by distance" algorithm, and tries other solutions until it finds a better one. For the flat problem, this algorithm finds better solutions than the initial one. Three scenarios have been generated with 20, 30 and 40 objective points and the algorithm has been iterating during a period of 8 hours in each case. The simulations of the flat problem have been performed using the following parameters: •V=1.2UD/UT •h=100 UD •Ta f =10 UT •Ωmax =1.5o/UT •Ψ=30o In addition, the iterated greedy search algorithm uses the following destruction and temperature parameters: •d=10 •T=3 The results of applying the IGS algorithm to the initial sequence obtained by searching the points with minimum distance are shown in table 4.8. Table 4.8 Observation time of the initial sequence ti and the optimized one topt for several scenarios of the flat problem. Number of targets tobs (UT)20 30 40 ti663 910 459 topt 639 889 434 This table shows that the improvement in the total observation time is small and this might be due to the following reasons: • The initial solution is very good, close to optimal and therefore the improvement of this solution is very difficult. •The iterated greedy search algorithm is stuck at local minimums and needs more diversification. • The problem is so restricted that the set of good solutions is very small and it is difficult to find better solutions. The observation paths obtained are shown in Figures 4.18, 4.19 and 4.20. It has been tried to optimize the orbital scenario over the Iberian Peninsula with 20, 30 and 40 points and with 1 and 2 satellites during 8 hours of calculation. However, it has not been possible to improve the initial solution obtained with the point selection algorithm by minimum distance. Moreover, the scenarios in Figure 4.10 could not be analysed since the calculation time per iteration exceeds 8 hours. Possible reasons why better observation paths have not been found are the following ones: 62 Chapter 4. Results and analysis • By increasing the computation time per iteration in the orbital case, fewer iterations are performed than in the flat problem. • The orbital problem is more restricted than the plane and therefore it is even more difficult to optimize the initial solution. Figure 4.18 Initial observation path on the left and optimized with Iterated Greedy Search on the right in a 20 targets scenario. 4.5 Iterated Greedy Search algorithm results 63 Figure 4.19 Initial observation path on the left and optimized with Iterated Greedy Search on the right in a 30 targets scenario. 70 List of Figures 4.9 h−Ψplots of various constellations in global scenario 49 4.10 Local simulations scenarios 50 4.11 Simulation example with a constellation of 3 satellites 51 4.12 h−Ψplots of various constellations in local scenario 52 4.13 Observation time tobs as a function of the number of satellites for a given altitude 53 4.14 Comparison of the observation time tobs for both point selection methods 55 4.15 Number of targets observed as a function of the number of points in the scenario for Meléndez’s experiment. Figure obtained from [5] 56 4.16 Comparison between simulations with targets observed after several attempts, and simulations with target points observed at once 58 4.17 Observation path in a simulation with repeated observation targets 59 4.18 Initial observation path on the left and optimized with Iterated Greedy Search on the right in a 20 targets scenario 62 4.19 Initial observation path on the left and optimized with Iterated Greedy Search on the right in a 30 targets scenario 63 4.20 Initial observation path on the left and optimized with Iterated Greedy Search on the right in a 40 targets scenario 64 List of Tables 4.1 Results of the simulations in the global scenario 47 4.2 Observation time of a constellation of 3, 4 and 5 satellites 48 4.3 Observation time of each scenario tobs in hours using "Sort by distance" 53 4.4 Number of points per satellite pass in each simulation using "Sort by distance" 54 4.5 Observarion time of each scenario tobs in hours using "Sort by topt " 56 4.6 Number of points per satellite pass in each simulation using "Sort by topt " 57 4.7 Observation time tobs in hours and number of points per pass in the 50 targets scenario 60 4.8 Observation time of the initial sequence tiand the optimized one topt for several scenarios of the flat problem 61 71 Bibliography [1] F. Pallarés Chamorro, Analysis and Optimization of Earth-Observation Micro-Constellations, Proyecto Fin de Carrera (Ingeniería Aeroespacial). Dep. Ingeniería Aeroespacial y Mecánica de Fluidos. E.T.S.I. Sevilla, 2019. [2] Earth Observation Portal, Mini-satellite deimos-2 mission data and specifications, extracted on 07/11/2020 from https://directory.eoportal.org/web/eoportal/satellite-missions/d/deimos-2. [3] ,Pleiades-HR (High-Resolution Optical Imaging Constellation of CNES), extracted on 25/10/2020 from https://earth.esa.int/web/eoportal/satellite-missions/p/pleiades. [4] Earth Observing System, PLÉIADES 1 CHARACTERISTICS, extracted on 25/10/2020 from https://eos.com/pleiades-1/. [5] J. Meléndez López, Algoritmos Heurísticos para la Adquisición de Múltiples Imágenes por Satélites Ágiles, Proyecto Fin de Carrera (Ingeniería Aeronáutica). Dep. Ingeniería Aeroespacial y Mecánica de Fluidos. E.T.S.I. Sevilla, 2015. [6] National Aeronautics and Space Administration, Image of the Earth from TIROS 1, extracted on 25/10/2020 from https://www.nasa.gov/vision/earth/features/bm_gallery_3.html. [7] ,TIROS, the Nation’s First Weather Satellite, extracted on 25/10/2020 from https://www.nasa.gov/multimedia/imagegallery/image_feature_1627.html. [8] T.Stützle R.Ruiz, "A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem", European Journal of Operational Research, vol 177, 2033–2049, 2007. [9] F. Montoro Sánchez, Análisis y Desarrollo de Algoritmos para la Planificación Optimizada de Adquisición de Imágenes por Satélites Ágiles, Proyecto Fin de Carrera (Ingeniería Aeronáutica). Dep. Ingeniería Aeroespacial y Mecánica de Fluidos. E.T.S.I. Sevilla, 2014. [10] Rafael Vázquez Valenzuela, Mecánica Orbital y Vehículos Espaciales. Apuntes de clase, E.T.S.I. Sevilla, 2020. Extracted from "http://www.aero.us.es/move/desc.html". [11] X.Wang, G.Wu, L.Xing, W.Pedrycz, Agile Earth observation satellite scheduling over 20 years: formulations, methods and future directions, extracted on 25/10/2020 from https://arxiv.org/pdf/2003.06169.pdf. 73