Estudio de Herramientas para el Análisis de Incertidumbre de Satélites en Órbita Baja Terrestre sometidos a Perturbaciones

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 Estudio de Herramientas para el Análisis de Incertidumbre de Satélites en Órbita Baja Terrestre sometidos a Perturbaciones Autor: Jesús Urbano Navajas Tutores: Rafael Vázquez Valenzuela oooooo. José Manuel Montilla García Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio de Herramientas para el Análisis de Incertidumbre de Satélites en Órbita Baja Terrestre sometidos a Perturbaciones Autor: Jesús Urbano Navajas Tutores: Rafael Vázquez Valenzuela Catedrático de Universidad José Manuel Montilla García Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado: Estudio de Herramientas para el Análisis de Incertidumbre de Satélites en Órbita Baja Terrestre sometidos a Perturbaciones Autor: Jesús Urbano Navajas Tutores: Rafael Vázquez Valenzuela José Manuel Montilla García 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: Agradecimientos Estas líneas pese a ser las primeras resultan ser las últimas que escribo y quizás las más complicadas de escribir también. Con este proyecto finalizo mi carrera y por fin puedo decir que soy ingeniero aeroespacial. Pero este logro no es solo mío, resultaría muy egoista si lo pensara, este logro es de mucha gente, toda la que me ha acompañado durante estos años. Este breve espacio del proyecto va dedicado a todas esa personas. A mis padres, Abundio y Ana. Un ejemplo continuo de cariño, trabajo y esfuerzo por sus hijos, esta carrera es casi más vuestra que mía. Con estos padres uno no puede estar más orgulloso de quien es y de donde viene. Gracias por inculcar en mi valores como el del trabajo bien hecho, la bondad, el cariño y el respeto, valores que gracias a vosotros llevaré conmigo toda la vida. Sois el mayor referente que uno puede tener. A mis hermanos, Jose y Juan Antonio. No seremos los hermanos que más hablen pero sí los que más se quieren. Gracias Jose por actuar siempre como hermano mayor y gracias Juan Antonio por enseñarme a afrontar la vida siempre con una sonrisa pese a las dificultades. A mi familia en general, a los que están y los que ya no. Espero que podáis estar orgullosos de mi. A mis amigos, los de Castro y los de Sevilla. Gracias por hacer este camino más ameno. Por las charlas interminables, los viajes para el recuerdo e incontables risas. Porque sois amigos para toda la vida. Y por último gracias a mis tutores, Rafa y José Manuel, por aconsejarme y antederme en todo momento y por haberme elegido para terminar mi carrera con este proyecto. Jesús Urbano Navajas Sevilla, 2023 I Resumen De bido al efecto de las perturbaciones, el cual es difícil de modelar con precisión, siempre aparece un cierto grado de incertidumbre asociado a la órbita. Para modelar y trabajar con dicha incertidumbre se puede utilizar la teoría de la probabilidad. El objetivo de este proyecto es el de estudiar y analizar distintos métodos que permitan la propagación de la incertidumbre a lo largo de la órbita, tales como el método Montecarlo, Linealización, Unscented Transform (UT) o los Elementos Orbitales Equinocciales. Asimismo, las perturbaciones más importantes para los satélites en órbita baja, el J2 y la resistencia aerodinámica, son estudiadas también. Además, el efecto que tienen en los elementos orbitales de la órbita y en la propia evolución de la incertidumbre es analizado de igual forma. III XÍndice 3.2 Variable aleatoria n-dimensional 24 3.3 Distribución Normal 25 4 Método Montecarlo 29 4.1 Definición 29 4.2 Desarrollo del método 29 4.2.1 Garantías con método Montecarlo y tamaño muestral 30 Generación de muestras 31 Modelo del sistema 31 Criterio de parada 31 Aumento del tamaño muestral 31 Ejemplo de uso 32 4.3 Resultados 32 4.4 Análisis de sensibilidad 35 4.4.1 Variación en excentricidad 36 4.4.2 Variación en semieje mayor 36 4.4.3 Variación en inclinación 37 4.5 Test de normalidad de Henze-Zirkler 37 5 Método de Linealización 41 5.1 Definición 41 5.2 Desarrollo del método 41 5.3 Resultados 43 5.4 Análisis de sensibilidad 45 5.4.1 Variación en excentricidad 45 5.4.2 Variación en semieje mayor 45 5.4.3 Variación en inclinación 45 6 Unscented Transform (UT) 49 6.1 Definición 49 6.2 Desarrollo del método 49 6.3 Resultados 50 6.4 Análisis de sensibilidad 51 6.4.1 Variación en excentricidad 52 6.4.2 Variación en semieje mayor 52 6.4.3 Variación con inclinación 53 6.4.4 Variación de κ53 7 Método de los Elementos Orbitales Equinocciales 59 7.1 Definición 59 7.2 Desarrollo del método 59 7.3 Resultados 61 7.4 Análisis de sensibilidad 61 7.4.1 Variación en excentricidad 62 7.4.2 Variación en semieje mayor 63 7.4.3 Variación en inclinación 63 8 Comparativa de métodos 67 8.1 Ventajas e inconvenientes de los métodos 67 Índice XI 8.1.1 Método Montecarlo 67 8.1.2 Método de Linealización 67 8.1.3 Unscented Transform 68 8.1.4 Elementos Orbitales Equinocciales Modificados 68 8.2 Comparativa de resultados 68 8.3 Análisis de sensibilidad con la excentricidad 71 9 Efecto de las perturbaciones 75 9.1 Efecto de las perturbaciones en la órbita 75 9.1.1 Efecto del J275 Consecuencias en la órbita 75 Resultados numéricos 76 9.1.2 Efecto de la resistencia aerodinámica 77 Resultados numéricos 77 9.2 Efecto de las perturbaciones en la incertidumbre 80 9.2.1 Efecto del J281 9.2.2 Efecto de la resistencia aerodinámica 82 9.2.3 Análisis de sensibilidad 82 10 Conclusiones y líneas futuras de trabajo 87 10.1 Conclusiones sobre la evolución de la incertidumbre y los métodos empleados 87 10.2 Conclusiones sobre el efecto de las perturbaciones 88 10.3 Líneas futuras de trabajo 88 Índice de Figuras 91 Índice de Tablas 93 Bibliografía 95 Notación aSemieje mayor eExcentricidad iInclinación ωArgumento de perigeo ΩAscensión recta del nodo ascendente o RANN MAnomalía media θAnomalía verdadera φLatitud λLongitud AR Ascensión recta δDeclinación LSímbolo astronómico de la Tierra JSímbolo astronómico del Sol $Símbolo astronómico de la Luna Primer punto de Aries Nodo ascendente Nodo descendente tTiempo TPeriodo orbital mvMasa del vehículo ρDensidad AÁrea pPresión de radiación solar εCoeficiente de reflectividad cVelocidad de la luz IFlujo de radiación solar CbCoeficiente balístico CDCoeficiente de resistencia aerodinámica vrel Velocidad relativa µParámetro gravitacional tg Función tangente arctg Función arcotangente sin Función seno cos Función coseno J2Perturbación del J2 XIII XIV Notación ω⊕Velocidad angular de la Tierra nVelocidad orbital media rVector posición vVector velocidad  hVector de momento específico nVector nodal eVector excentricidad LEO Órbita terrestre baja MEO Órbita terrestre media GEO Órbita geoestacionaria σDesviación típica o estándar σ2Varianza Cov(X1,X2)Covarianza de las variables X1yX2 ΣMatriz de covarianza µMedia N= (µ,σ2) Distribución normal o gaussiana de media µ y varianza σ2 PDF Función de densidad de propabilidad MC Método Montecarlo LIN Método de Linealización UT Método de la Unscented Transform EOE Elementos Orbitales Equinocciales EOK Elementos Orbitales Keplerianos EOEM Elementos Orbitales Equinocciales Modificados SCC Sistema de Coordenadas Cartesianas p,f,g,h,k,LElementos orbitales equinocciales modificados 1 Introducción 1.1 Motivación Resulta innegable el aporte que ha tenido el espacio en el desarrollo de la humanidad en sí misma. Desde el primer satélite puesto en órbita en la historia, el Sputnik I en 1957 por la URSS, el cual solo emitía pulsos de radio a cierta frecuencia, hasta los lanzamientos cada vez más frecuentes de la constelación de satélites Starlink de SpaceX, han pasado apenas 66 años. En este tiempo se han lanzado alrededor de 8900 satélites por más de 40 países. Según una estimación hecha en 2018 por las Naciones Unidad [ 3 ], unos 5000 permanecen en órbita. De ellos, unos 1900 están operativos, mientras que el resto ha cumplido con su vida útil y se han convertido en basura espacial. Los satélites se utilizan para diversos propósitos. Entre las modalidades más frecuentes se incluyen a los satélites de observación terrestre, tanto civiles como militares, satélites de comunicaciones, satélites de navegación o satélites meteorológicos. Además, no solo hemos lanzado satélites para orbitar a la Tierra, más de una docena de sondas espaciales han sido situadas en órbita alrededor de otros cuerpos celestes y han pasado a ser satélites artificiales de la Luna, Mercurio, Venus, Marte, Júpiter, Saturno, un cometa, el Sol y unos cuantos asteroides. Las órbitas de los satélites varían considerablemente en función de la finalidad de este. Las más usadas son la órbita baja terrestre (LEO 1 ) donde se ubica el 63% de los satélites operativos y la órbita geoestacionaria (GEO2) donde está el 37% restante. Sin embargo, aunque resulte sorprendente, la posición de todos estos satélites no se conoce exactamente. Debido a numerosos factores como las perturbaciones a las que se encuentran sometidos o maniobras realizadas siempre existe un grado de incertidumbre. Es decir, la posición del satélite no es un punto sino una nube de puntos probabilística donde las empresas estiman que está la posición de su satelite. Esto hace evidente la necesidad de una sólida y exhaustiva estimación de la posicíón del satélite, es decir, de esta nube de puntos. Para ello las empresas deben de contar con métodos y herramientas para su estimación y su análisis en todo instante de tiempo. Deben conocer cómo evoluciona dicha incertidumbre. Pues solo así se evitarían desgracias como la ocurrida en febrero de 2009 entre un satélite comercial en uso de la empresa Iridium Satellite LLC con otro satélite ruso declarado fuera de operaciones. El choque resultó en unos 600 fragmentos que se unían a las aproximadamente 9000 toneladas de basura espacial que estima la ESA que se encuentran actualmente orbitando alrededor de nuestro planeta. 1Órbitas con altura inferior a 2000 km. 2Órbita circular en el plano ecuatorial con una altura de 35786 km. 1 2Capítulo 1. Introducción Otro incidente más reciente ocurrió el 31 de octubre de 2022. China lanzaba un cohete con el objetivo de completar la última fase de su estación espacial. Tras ponerlo en órbita, las autoridades chinas permitieron que uno de sus módulos, de alrededor de 20 toneladas de peso, se desprendiera y cayera sin control hacia la Tierra. Esta caída descontrolada, a pesar de estar monitorizada por diversos organismos, provocó que varios aeropuertos de España y parte de Europa, como Francia e Italia, decidieran parar la actividad aérea, lo cual generó importantes retrasos en el tráfico aéreo. Afortunadamente, el 4 de noviembre de 2022 se confirmaba la caída del objeto en el océano pacífico. Con unos métodos para estimar la posición, la nube de puntos, de manera más precisa se podrían predecir posibles colisiones como estas, y evitarlas con las maniobras pertinentes. Para que dicha estimación sea lo más exacta posible, en estos métodos y modelos habrá que incluir las perturbaciones, que propiamente son las que provocan dicha incertidumbre. Así, como nuestro estudio se va a centrar en órbitas LEO, pues es donde están la mayoría de satélites, resultaría de lo más interesante estudiar como las perturbaciones más importantes para este tipo de órbita, el J2y la resistencia aerodinámica, interactúan con la incertidumbre. Este proyecto además surge también como complementación del trabajo pasado realizado sobre este tema en [10]. 1.2 Objetivo y alcance del proyecto El objetivo principal de este proyecto es el de analizar cómo evoluciona la incertidumbre a partir de errores iniciales a lo largo de una órbita tipo LEO con el paso del tiempo. El error inicial que se usará en este proyecto para analizar la evolución de la incertidumbre será de 50 m en posición en las tres direcciones. Además, se estudiarán diferentes métodos para realizar dicho análisis. Se pretende también estudiar las ventajas y desventajas que tienen todos los métodos, los resultados obtenidos con ellos y la determinación de los mejores métodos para llevar a cabo esta tarea. Por otro lado, se llevará a cabo distintos análisis de sensibilidad con todos los métodos para comprobar cuál es el efecto que tendría en la evolución de la incertidumbre en la órbita un cambio en las propiedades de ella. Finalmente, se procede a estudiar el efecto de las perturbaciones más importantes que sufren los satélites en órbita LEO, así como los cambios que provocan en la órbita y sobre todo en la evolución de la incertidumbre. 1.3 Estructura del proyecto El documento se encuentra estructurado de la siguiente manera: Capítulo 1: El capítulo empieza presentado la motivación principal que lleva a la realización de este proyecto. Además de esto se explica también el objetivo, alcance y estructura del mismo. Capítulo 2: Los conceptos básicos de la mecánica orbital más necesarios para el seguimiento del proyecto son introducidos: sistemas de referencias, ecuaciones del movimiento, elementos orbitales, así como una introducción a las perturbaciones. Capítulo 3: En este capítulo se presentan algunos conceptos básicos de estadística para facilitar el seguimiento del proyecto: variables aleatorias unidimensionales, n-dimensionales, distribución normal, momentos estadísticos, distribuciones continuas, distribuciones discretas. Capítulo 4: A partir de este capítulo y hasta el Capítulo 7 se presenta los diferentes métodos que se van a usar para analizar la evolución de la incertidumbre en la órbita. En este capítulo en concreto se estudia el método Montecarlo, se explica sus fundamentos, los resultados que se obtienen y los diferentes análisis de sensibilidad realizados. Capítulo 5: El método de Linealización es presentado en este capítulo, se explica sus fundamentos, los resultados que se obtienen y los diferentes análisis de sensibilidad realizados. 1.3 Estructura del proyecto 3 Capítulo 6: Se desarrollan las bases del método de la Unscented Transform (UT), los resultados obtenidos y los análisis de sensibilidad realizados. Capítulo 7: Para concluir con los métodos de estudios se analiza el método basado en los Elementos Orbitales Equinocciales, sus resultados y los análisis de sensibilidad. Capítulo 8: En este capítulo se compara los resultados obtenidos con los cuatro métodos y se presentan las ventajas e inconvenientes de la utilización de cada uno de ellos. Capítulo 9: En esta parte del proyecto se desarrolla más el efecto de las perturbaciones comentado en el Capítulo 2. Se centra en el estudio de dos de ellas: el J2 y la resistencia aerodinámica, las más importantes para LEO. Se estudia tanto su efecto en la propia órbita como en la evolución de la incertidumbre. Capítulo 10: Para terminar, se presenta en este capítulo todas las conclusiones extraídas y las posibles líneas de trabajo futuras. 2 Conceptos básicos de mecánica orbital El objetivo de este capítulo es presentar los conceptos y herramientas matemáticas básicas de la mecánica orbital que utilizaremos a lo largo de este proyecto y explicar los diferentes sistemas de referencia y su uso. Además, se explicará la teoría de perturbaciones así como los diferentes efectos que pueden tener en nuestro cuerpo de estudio. 2.1 Sistemas de referencia En la mecánica orbital, entendida esta como la rama de la ciencia encargada de la aplicación de la matemática y la física para la resolución de problemas prácticos relacionados con el movimiento de los astros y/o vehículos espaciales, es vital hacer una correcta y meticulosa elección del sistema de referencia bajo el cual vamos a trabajar, pues de no ser así podría llevarnos a problemas con una dificultad añadida. El sistema de referencia elegido dependerá en muchas ocasiones de cuál sea el cuerpo celeste al que nuestro objeto de estudio, normalmente un satélite, orbite. Encontrándonos así por ejemplo con sistemas de referencia heliocéntricos cuyo centro será el Sol ( J ) o sistemas de referencia geocéntricos con el centro en la Tierra (L). Además, existen pues otros sistemas de referencia auxiliares que serán útiles para ciertas misiones/tareas como puede ser simplificar la visualización de una órbita, lanzamientos espaciales o la toma de medidas astronómicas entre otras. Existen muchas clasificaciones distintas para dichos sistemas de referencia pero una de las más extendidas es la de dividir los sistemas de referencia en inerciales y no inerciales. Hay que tener especial cuidado cuando hablamos de sistemas de referencia inerciales pues formalmente, dicho sistema es aquel el cual se encuentra en reposo o en movimiento con velocidad constante (carece de aceleración) y como sabemos en el contexto que nos ocupa el centro del sistema de referencia va a estar situado en un cuerpo del Sistema Solar (el Sol , la Tierra u otro planeta) con lo que no son realmente inerciales; no obstante, en la práctica se desprecian las fuerzas de inercia (aceleraciones de arrastre, fuerzas de inercia, fuerzas de Coriolis) ya que suelen ser muy pequeñas frente al resto de fuerzas que actúan en el sistema (fuerza gravitatoria). 2.1.1 Sistemas de referencia inerciales Sistema Heliocéntrico El sistema de referencia heliocéntrico, el cual denotaremos con el superíndice ( J ), tiene su origen en el Sol. El plano que queda definido por el origen y las direcciones X e Y es el de la eclíptica, 5 12 Capítulo 2. Conceptos básicos de mecánica orbital • Anomalía media ( M ): Este valor nos indica la posición de nuestro cuerpo en la órbita, también se podría usar la anomalía verdadera θ pero M presenta mayores facilidades matemáticas a la hora de analizar perturbaciones. M=rµ a3t(2.17) Donde trepresenta al tiempo desde el último paso por el periapsis. Figura 2.7 Elementos orbitales de una órbita. 2.3 Teoría de perturbaciones Como hemos dicho en la Sección 2.2 el modelo de los dos cuerpos es poco realista si no incluimos las perturbaciones debido a las hipótesis realizadas. Con la primera hipótesis considerábamos nuestro sistema aislado del universo obviando otras fuerzas como pueden ser la resistencia atmosférica, la presión de radiación solar o la fuerza gravitatoria ejercida por otros cuerpos que en ocasiones pueden resultar bastante importantes. Con la segunda hipótesis se consideraban a los cuerpos de estudio con forma perfectamente esférica y por ende con una distribución de masas constante, afirmación errónea pues la distribución de masa va a depender de la latitud (φ) y la longitud (λ). A continuación se presenta un estudio de las principales fuentes de perturbación en nuestras orbitas. 2.3.1 Efecto de un tercer cuerpo Al inicio del Punto 2.3 rebatimos ya la primera hipótesis realizada ya que nuestro sistema no se encuentra realmente aislado del resto del universo pues, como es el caso, se ve afectado por ejemplo por la presencia de otros cuerpos celestes y por ende por las fuerzas que dichos cuerpos ejercen sobre nuestro sistema. A modo explicativo se va a proceder al estudio de la fuerza de perturbación a la que se ve sometido un cuerpo poco masivo, como es el caso de un satélite, orbitando la Tierra por la presencia de la Luna ( $ ) ya que como veremos al final de este apartado esta será la perturbación más importante para un cuerpo orbitando a la Tierra provocada por un tercer cuerpo. Una representación de este problema se puede ver en la Figura 2.8 2.3 Teoría de perturbaciones 13 Figura 2.8 Objeto orbitando a la Tierra con el efecto de la Luna [11]. En este problema la ecuación de movimiento, que quedó definida por la ecuación (2.6), la reescribiremos en el sistema de referencia inercial centrado en el centro de masas del sistema Tierra-Luna, (CML$): ¨ rIN =−µL r r3+µ$r$−r |r$−r|3(2.18) Donde r$=r2−r1(2.19) y se verifica por estar centrado nuestro sistema en el (CML$) la siguiente expresión: mLr1+m$r2= 0(2.20) Además sabemos que la ecuación de movimiento de la Luna queda: ¨ r$=−G(mL+m$)r$ r3 $ (2.21) Por tanto combinando adecuadamente las ecuaciones (2.19), (2.20) y (2.21) llegamos a: ¨ r1=µ$r$ r3 $ (2.22) Finalmente como podemos deducir de la Figura 2.8, rIN =r+r1 luego la expresión para la ecuación de movimiento teniendo en cuenta la presencia de la Luna es: ¨ r=−µL r r3+µ$r$−r |r$−r|3−µ$r$ r3 $ = γk+ γp(2.23) En donde se ha definido por comodidad  γk=−µL r r3(2.24)  γp=µ$r$−r |r$−r|3−µ$r$ r3 $ (2.25) la fuerza kepleriana clásica (  γk ) y la fuerza de perturbación debida al tercer cuerpo (  γp ), en nuestro caso la Luna. 14 Capítulo 2. Conceptos básicos de mecánica orbital Linealizando ahora la fuerza de perturbación podemos llegar a una expresión más sencilla del siguiente tipo:  γp≈µ$ r$"3r$  r$·r r$−r#(2.26) Donde se puede deducir fácilmente que el orden de magnitud para estas perturbaciones es de  γp≈µ$ r r$3(2.27) Ahora para que resulte más fácil vislumbrar el efecto de estas perturbaciones podemos realizar el cociente entre  γpy γkquedando finalmente:  γp γk≈µ$ µLr r$3 (2.28) para el caso de la Luna, mientras que para el caso general del efecto causado por la presencia de cualquier cuerpo del Sistema Solar al que denotaremos por el subíndice Cquedaría:  γp γk≈µC µLr rC3 (2.29) Como ya avanzamos las perturbaciones más importantes que va a ver un cuerpo orbitando a la Tierra por la presencia de un tercer cuerpo son las causadas por la presencia de la Luna, pues su valor es el de 10−5 para GEO y 10−7 para LEO, mientras que las causadas por el Sol serán del orden de 5·10−6 para GEO y 5·10−8 para LEO. Las perturbaciones producidas por los demás cuerpos del Sistema Solar quedan despreciadas frente a estas, las cuales comúnmente son llamadas perturbaciones lunisolares. 2.3.2 Forma no esférica de la Tierra: Al principio de la Sección 2.3 se discutió como la segunda hipótesis no puede ser válida pues la distribución de masas para un cuerpo cualquiera dependerá de la latitud ( φ ) y la longitud ( λ ) es este hecho el que introduce una nueva perturbación que definimos como:  γ=∇U(2.30) Donde U es el potencial gravitatorio el cual definiremos en coordenadas esféricas U=U(r,φ,λ) y vendrá dado para un cuerpo con forma arbitraria por la solución de la Ecuación de Laplace: 1 r2 ∂ ∂rr2∂U ∂r+1 r2cosφ ∂ ∂φ cosφ∂U ∂φ (2.31) Cuya solución general para el caso que nos ocupa, la cual viene desarrollada en [12], resulta: U=µ r"1− ∞ ∑ n=2 n ∑ m=0 Jnm RL rn pnm(sinφ)cos(m(λ−λnm))#(2.32) Expresión en la cual el primer término representa el potencial de una esfera, mientras que el resto (la serie) representa la desviación del modelo esférico, o cuan diferente es el potencial respecto a un cuerpo esférico. Los coeficiente Jnm y λnm son los asociados al armónico nm mientras que pnm representa el polinomio de Legendre de grado ny orden m, definido como: 2.3 Teoría de perturbaciones 15 pnm(x) = (1−x2)n/2dmpn(x) dxm(2.33) El argumento trigonométrico de los polinomios de Legendre en la ecuación (2.32) constituye armónicos esféricos superficiales permitiendo diferenciar dichos armónicos en zonales, sectoriales y teserales. • Armónicos zonales: Estos armónicos están definidos para cuando m=0 lo cual lleva al hecho de que el potencial no sufra variaciones respecto a λ y por tanto los armónicos coinciden con los del modelo con simetría de revolución. Entre dichos armónicos el J2 sin duda es la perturbación más importante debida a la forma no esférica de la Tierra pues el J22 y el J3 , las siguientes perturbaciones más importantes, son del orden de 1000 veces más débiles. • Armónicos sectoriales: Para cuando n y m coinciden provocando que no se produzcan variaciones en el potencial respecto a φ, dividiendo la esfera en sectores. • Armónicos teserales: Se producen cuando el grado y el orden, n y m , son distintos entre sí y a su vez distintos de 0, en este caso el potencial varía tanto con φ como con λ dividiendo por ende la esfera en cuadrados esféricos. Para una mejor comprensión de dichos armónicos se representan en la Figura 2.9. Figura 2.9 Representación de los distintos tipos de armónicos [8]. En esta visualización, donde l representa nuestra n , las áreas rojas representan U> 0, áreas azules representan U< 0 y las áreas blancas representan U≈ 0. Las líneas nodales del armónico están entre las áreas rojas y azules y correspondan a U=0 [8]. 2.3.3 Resistencia atmosférica: Dicha perturbación resulta la más importante en órbitas bajas y órbitas muy excéntricas con perigeo bajo. El efecto de esta en la órbita es la de su circularización, haciendo disminuir el radio de apogeo hasta que coincide con el de perigeo y una vez la órbita es circular el efecto será el de una caída en espiral hasta la reentrada del vehículo, efecto que se puede apreciar en la Figura 2.10 La resistencia tiene la dirección de la velocidad relativa, vrel , entendida esta como la diferencia entre la velocidad del cuerpo y la velocidad de la atmósfera, pero sentido opuesto: 16 Capítulo 2. Conceptos básicos de mecánica orbital Figura 2.10 Efecto de la resistencia aerodinámica [11].  γp=−SCD 2mV ρvrel vrel (2.34) En donde encontramos que mV es la masa del vehículo, S la superficie frontal, CD el coeficiente de resistencia aerodinámico y ρ la densidad del medio, precisamente por dicha densidad esta perturbación resulta de lo más complicada para modelar puesto que depende de múltiples factores: la naturaleza de la estructura molecular de la atmósfera, el flujo solar incidente o las interacciones geomagnéticas. Ahora bien existen diversos modelos que intentan caracterizar el comportamiento de la densidad como son los modelos estáticos y los modelos variables con el tiempo. Los modelos estáticos son los más simples pero también los menos precisos pues se basan en asumir todos los parámetros atmosféricos constantes mientras que los modelos variables con el tiempo son más fieles a la realidad pero también mucho más complejos pues tienen en cuenta efectos como las fluctuaciones diurnas, el ciclo de rotación solar de 27 días, rotación de la atmósfera entre otros [12]. 2.3.4 Presión de radiación solar: Como la resistencia, la presión de radiación solar es una perturbación no conservativa que se convierte más importante para grande altitudes. Uno de los aspectos más difíciles a la hora de analizar la radiación solar es el de modelar y predecir los ciclos y variaciones solares. Durante periodos de tormentas solares, episodios de mayor emisión de radiación solar, este efecto puede ser mucho mayor que cualquier otra perturbación, dependiendo claro está de la altitud en la que nos encontremos. Sin embargo, para periodos de baja actividad solar este efecto puede ser incluso despreciable frente a los demás [12]. Dicha perturbación como hemos avanzado está producida por la incidencia de la luz solar (fotones) en una superficie produciendo en ella un efecto mecánico, la presión de radiación solar la cual produce una perturbación que puede ser modelada como:  γp≈−pA(1+ε) mV eJ(2.35) Donde A es el área, ε el coeficiente de reflectividad que dependerá si el cuerpo es reflectante, absorbente o transparente, mV la masa del vehículo, eJ apunta en la dirección del Sol desde el vehículo y p es la presión que viene dada por p=I c siendo c la velocidad de la luz e I el flujo medio de radiación solar que se puede calcular como: 2.4 Clasificación y metodología para los efectos de perturbación 17 I=ILLL L2 (2.36) Expresión en la cual LL es la distancia Sol-Tierra, L la distancia del Sol a nuestro cuerpo e IL el flujo medio de radiación solar a 1 AU (149597870 km) del Sol (IL=1367 W m2). Para el estudio de esta perturbación también se tienen que tener en cuenta periodos de eclipses, pues ahí la fuerza desaparece, la sombra de la Tierra o el albedo de la Tierra, reflexión de la luz del Sol en la Tierra. Tanto en la Figura 2.12 como en la Figura 2.13 se muestra el efecto de todas estas perturbaciones para órbitas bajas y elevadas. 2.4 Clasificación y metodología para los efectos de perturbación Los efectos que producen todas estas perturbaciones en los elementos orbitales explicados en el Punto 2.2.2 permiten clasificar dichos efectos en: • Seculares: Los cuales tienen un comportamiento monótono en las variaciones que provocan en los elementos orbitales. Su efecto por ende es acumulativo, y por ello son las más importantes cuando se hace un estudio de perturbaciones. • Periódicos: Son aquellos cuya variación se repite cada cierto periodo ya sea de largo periodo, mayor que el periodo orbital2o de corto periodo, del orden del periodo orbital. En la Figura 2.11 se presenta la clasificación de los efectos que producen cada perturbación estudiada para cada elemento orbital. Para terminar con el tema de las perturbaciones, es necesario comentar los métodos más importantes hoy en día que se tienen para afrontar el estudio de los efectos de las perturbaciones: • Métododo de perturbaciones generales: Emplea de base los elementos orbitales y estudia como varían con el tiempo. • Métododo de perturbaciones especiales: No emplea elementos orbitales sino directamente las ecuaciones del movimiento, (2.6). Este a su vez se puede dividir en dos tipos de tratamiento del problema: – Método de Encke: Resolver directamente las ecuaciones del movimiento incluyendo perturbaciones y tantos cuerpos como sea necesario. – Método de Cowell: Resuelve las ecuaciones de una corrección δr sobre una solución de referencia de los dos cuerpos ¨ r=−µr r3 . Si la corrección crece es necesario recalcular la solución de referencia. 2Tiempo que tarda el cuerpo en completar una revolución en su órbita. 18 Capítulo 2. Conceptos básicos de mecánica orbital Figura 2.11 Clasificación de los efectos de las perturbaciones [11]. Figura 2.12 Efecto de perturbaciones en órbita baja [11]. 2.4 Clasificación y metodología para los efectos de perturbación 19 Figura 2.13 Efecto de perturbaciones en órbita elevada [11]. 3 Conceptos básicos de estadística La estadística es una rama de las matemáticas que se enfoca en la recolección, análisis e interpretación de datos. Su objetivo es describir y resumir los datos, así como inferir y hacer predicciones sobre la población 1 a partir de una muestra de datos. También se utiliza para evaluar la fiabilidad de las conclusiones que se obtienen a partir de los datos y para tomar decisiones basadas en ellos utilizando técnicas y métodos matemáticos para modelar dichos datos y extraer información útil de ellos. En resumen, la estadística es una herramienta fundamental en la investigación científica, el análisis de datos, la toma de decisiones y la resolución de problemas en una amplia variedad de campos. En este capítulo se van a presentar dichas técnicas y métodos matemáticos cuya comprensión es clave para el seguimiento de un proyecto como este. 3.1 Variable aleatoria unidimensional Una variable aleatoria unidimensional es una variable numérica que puede tomar valores numéricos aleatorios en un conjunto dado con una probabilidad asociada a cada valor posible. En otras palabras, una variable aleatoria unidimensional es una función que asigna un valor numérico a cada posible resultado de un experimento aleatorio, entendido este como un evento o proceso que puede producir diferentes resultados cuando se realiza en condiciones similares. 3.1.1 Distribuciones discretas y continuas La variable aleatoria unidimensional puede seguir una distribución discreta o continua siendo la principal diferencia el conjunto de valores que puede tomar la variable aleatoria. Una variable aleatoria discreta toma valores que son separados e individuales, como números enteros o contables. Es decir, los valores que puede tomar la variable se pueden enumerar mientras que una variable aleatoria continua toma valores que están en un intervalo continuo. Distribución discreta Las distribuciones discretas, las que siguen las variables discretas aleatorias, quedan definidas por una función de probabilidad o función de cuantía la cual describe la probabilidad de que una variable aleatoria discreta tome cada uno de sus posibles valores. p(x) = P(X=x)(3.1) En donde p(x) es la función de probabilidad, X la variable aleatoria discreta y x los posibles valores que puede tomar X. 1Conjunto de elementos o eventos similares que son de interés para alguna pregunta o experimento. 21 4 Método Montecarlo 4.1 Definición El método Montecarlo (MC) es una técnica numérica utilizada para resolver problemas matemáticos complejos mediante el uso de simulaciones estadísticas y aleatorias. Este método se basa en la generación de números aleatorios para simular un gran número de casos posibles de un problema, lo que permite estimar probabilidades y distribuciones de resultados. Este método es particularmente útil cuando no es posible obtener una solución analítica exacta a un problema, o cuando se requiere tomar en cuenta múltiples variables y posibles resultados. Este método puede proporcionar una aproximación de la solución con una precisión y confiabilidad deseada, pero requiere una gran cantidad de cálculos y simulaciones para producir resultados precisos. El método Montecarlo se utiliza en una amplia variedad de campos como la física, la ingeniería, la economía, la finanzas, las ciencias sociales y la informática, entre otros. Una definición más exhaustiva del método la podemos encontrar en [6]. 4.2 Desarrollo del método Para el caso que nos ocupa, dentro del marco de la mecánica orbital, este método se va a utilizar para observar el efecto de la incertidumbre en posición y velocidad en órbitas geocéntricas y ver su evolución en el tiempo. Para ello se genera una primera nube de puntos, del orden de los 1000, que seguirá una distribución normal alrededor del punto de referencia, el cual será la posición y velocidad inicial en nuestra orbita, r0yv0. Así, la media inicial de dicha distribución normal será las condiciones iniciales de posición y velocidad. Por otro lado, la matriz de covarianza inicial será diagonal con una desviación típica inicial en posición de 50 m en cada eje y la correspondiente en velocidad 1 . Calculada esta tomando derivada en la fórmula de la velocidad orbital: v=rµ r(4.1) Una vez tenemos esta primera nube de puntos, con el método Montecarlo vamos a propagar cada uno de estos puntos. Esto no es más que resolver las ecuaciones del problema de los dos cuerpos definidas en (2.6) para cada punto. Así, una vez que tenemos la nube de puntos propagada para un tiempo determinado podemos calcular fácilmente su nueva media y su matriz de covarianza y con ello a través de la expresión 1Para todos los métodos siguientes se usará esta misma matriz de covarianza inicial. 29 30 Capítulo 4. Método Montecarlo (3.17) el nuevo elipsoide de confianza viendo con esto como evoluciona la incertidumbre en nuestra órbita. En el Punto 4.3 se va a presentar los resultados obtenidos con este método para una órbita concreta, así como la realización de un test de normalidad para ver cuando deja dicha distribución de comportarse como una distribución normal. Además, se procederá a realizar un análisis de sensibilidad cambiando los elementos orbitales iniciales, que no es más que cambiar las condiciones iniciales con las que se realiza el método. 4.2.1 Garantías con método Montecarlo y tamaño muestral Al principio, en el Punto 4.2, decidimos que el número de puntos para el uso de este método sería de 1000. Sin embargo, no se justificó esta elección. A continuación, se va a presentar un algoritmo para poder tener ciertas garantías con el método Montecarlo, así como poder elegir adecuadamente un tamaño muestral. El algoritmo se ha extraído íntegramente de [14]. Figura 4.1 Diagrama del algoritmo. Como vemos en la Figura 4.1 el algoritmo se divide en tres grandes bloques, la generación de las muestras, el modelo de estudio y el criterio de parada. A continuación, se explica más en detalle cada uno de ellos. 4.2 Desarrollo del método 31 Generación de muestras Una de las características esenciales del método de Montecarlo es el uso de números aleatorios, por lo tanto, la forma en la que estos se generan constituye un punto importante a la hora de usar el método. Para nuestro caso, como vimos en el Punto 4.2, la nube de puntos inicial seguirá una distribución normal. De tal forma que la generación de la muestra se ha realizado a través del comando mvnrnd de MATLAB © , el cual con la matriz de covarianza y la media de la distribución normal es capaz de generar una nube de puntos aleatorios del tamaño elegido. Modelo del sistema El modelado del sistema es un aspecto clave del proceso, ya que la precisión de los resultados obtenidos en las simulaciones dependerá de la precisión del método y de la exactitud con la que el modelo representa el sistema físico que se pretende estudiar. Por lo tanto, es necesario un modelo que emule de la manera más exacta posible el comportamiento del sistema real. Para el estudio de la evolución de la incertidumbre en las órbitas usando el método Montecarlo está claro el modelo a emplear. Este modelo es el que se adelantó en el Punto 4.2 y no son más que las ecuaciones del problema de los dos cuerpos definidas en la expresión (2.6). De tal forma que el objetivo del método es resolver para cada uno de los puntos de la nube las ecuaciones y observar así su evolución. Criterio de parada En [ 14 ] se discute el criterio de parada que se usará para determinar el fin de las simulaciones realizadas con Montecarlo y por ende el tamaño muestral. Para ello se presentan dos posibles criterios. El primero de ellos consiste en comparar los valores obtenidos en la iteración con los obtenidos en varias iteraciones anteriores. Esta comparación consiste en restar los valores obtenidos en las diferentes iteraciones, calculando así los errores absolutos y comprobando que en varias iteraciones sucesivas estos errores estén por debajo de la tolerancia fijada. El motivo por el cual es necesario comparar con varios valores anteriores es que al comparar con un único valor anterior puede suceder que el error sea pequeño por “casualidad” sin haber llegado a la convergencia del método. El segundo procedimiento consiste en realizar varias simulaciones para un tamaño de la muestra elevado. Representando posteriormente el error cometido con respecto al tamaño de la muestra para cada una de las simulaciones, se obtiene una envolvente de errores que nos permite determinar el tamaño de la muestra que nos garantiza la precisión fijada. En nuestro caso resultaría más interesante usar el segundo de los procedimientos debido a las limitaciones en la memoria del ordenador, que podría verse colapsada en el caso de utilizar el primero de los procedimientos. Aumento del tamaño muestral En las sucesivas iteraciones del método será necesario aumentar el tamaño de la muestra en caso de buscar unos resultados más precisos. A la hora de cuantificar el aumento del tamaño de la muestra en cada iteración existen diferentes opciones. En [ 14 ] se detallan dos de ellas: el procedimiento aditivo, el cual en cada iteración añade un número de nuevas muestras del mismo tamaño que la muestra inicial, y el procedimiento multiplicativo, el cual añade en cada iteración un número de nuevas muestras del mismo tamaño que las muestras acumuladas hasta ese momento. El segundo de ellos, el procedimiento multiplicativo, resultaría más adecuado ya que el procedimiento aditivo tiene como inconveniente que al aumentar el número de iteraciones las simulaciones antiguas tienen más peso que las nuevas. 32 Capítulo 4. Método Montecarlo Ejemplo de uso Con el algoritmo ya definido se presenta en la Figura 4.2 los resultados obtenidos en [ 14 ] para el ejemplo de la estimación del número πcon el método Montecarlo. Figura 4.2 Reducción de la envolvente de error con el tamaño de la muestra. En esta figura se muestra la relación entre los errores relativos cometidos y el tamaño de la muestra. De tal forma que con el error relativo que nosotros exijamos a nuestros resultados podemos obtener el tamaño de muestra necesario. La diferencia entre los errores obtenidos en azul y en rojo es que en azul se obtiene el error relativo a partir del valor de π teórico y el caso del rojo es para cuando no se conozca los resultados a obtener, de tal manera que este error relativo se obtiene a partir del resultado obtenido en la simulación anterior. Ahora, para nuestro caso de estudio, la elección ha sido de un tamaño muestral de 1000 puntos pues analizando los resultados obtenidos con este algoritmo y con los obtenidos también de [ 10 ] se determina que este tamaño es suficientemente grande como para dar resultados fiables. Se deja como línea futura de trabajo realizar un análisis como el de la Figura 4.2 para el caso de la evolución de la incertidumbre en órbitas. 4.3 Resultados Para hacer uso del método primero de todo debemos de definir cuál va a ser nuestra órbita de referencia y sus elementos orbitales iniciales, Tabla 4.1. Una vez que tenemos esto definido podemos hacer uso del método y observar como es la evolución de la incertidumbre, es decir, de los elipsoides de confianza tras una revolución. Como se observa en la Figura 4.3 vemos como al principio este elipsoide empieza prácticamente como una esfera, pues como dijimos en el Punto 4.2 las desviaciones típicas iniciales van a ser las 4.3 Resultados 33 Tabla 4.1 Elementos orbitales iniciales. a e i ωΩM 6878 km 0.01 20◦90◦45◦0◦ Figura 4.3 Elipsoides de confianza en una revolución (MC). mismas en los 3 ejes. Sin embargo, a medida que pasa el tiempo se observa que en el sistema local LHVH , definido este con el eje OX como el eje radial el OY perpendicular al OX en el mismo plano en el sentido del movimiento y OZ completando el triedro, hay un crecimiento de la incertidumbre sobre todo en dirección OY , convirtiéndose con el paso del tiempo en prácticamente una línea. Primero de todo, el crecimiento de la incertidumbre se debe a la naturaleza no lineal de las ecuaciones del problema de los dos cuerpos, que hace que pequeñas variaciones en las condiciones iniciales del problema se conviertan en grandes variaciones con el paso del tiempo, lo que se denomina ”efecto mariposa”. Por otro lado el crecimiento mayor en dirección OY local es debido simplemente a la 3ªLey de Kepler 2 . Las órbitas con menor semieje mayor tienen menos periodo orbital. Por ello se va a producir un desfase entre los diferentes puntos de la nube (los que tienen más semieje mayor se van quedando atrás, y al contrario). Este desfase se traduce en una mayor dispersión en la dirección OY local, mientras que en las otras direcciones locales no vemos tantos cambios. Para entender mejor estos resultados sería adecuado presentar como varía con el tiempo la matriz de covarianza, es decir, tanto las desviaciones típicas en los tres ejes como la suma de los autovalores, la cual nos da una idea del tamaño de la incertidumbre. Esto se puede ver en las Figuras 4.5 y 4.6. En estas figuras observamos el efecto visto en la Figura 4.3. Tanto la suma de autovalores como las desviaciones típicas tienen por un lado una componente secular creciente y por otro lado tienen una componente periódica de periodo igual al periodo orbital, pues, será en el perigeo cuando mayor sea dicha incertidumbre al haber aquí mayores cambios en velocidad y posición, lo que se traduce en 2T=2πqa3 µ 34 Capítulo 4. Método Montecarlo Figura 4.4 Representación de ejes locales. Figura 4.5 Comparación desviaciones típicas (MC). una mayor dispersión de la incertidumbre. Además, en cuanto a las desviaciones típicas se observa lo mencionado anteriormente, el claro predominio del crecimiento en dirección OY local. También resultaría interesante mostrar como sería este elipsoide de confianza si en vez de en una revolución lo propagamos durante mucho más tiempo, días por ejemplo. El resultado se muestra en la Figura 4.7, y es que al cabo de cierto tiempo el elipsoide de confianza deja de ser representativo pues, la distribución deja de comportarse como una gaussiana. De hecho, se muestra como la nube de puntos propagada empieza a curvarse por los extremos, cosa que el elipsoide no puede capturar. El tiempo concreto a partir del cual deja de comportarse como una distribución normal se presentará en el Punto 4.5, evidentemente cuanto más tiempo pase más evidente será esta diferencia. 4.4 Análisis de sensibilidad 35 Figura 4.6 Suma de autovalores (MC). Figura 4.7 Elipsoide de confianza versus nube de puntos propagada tras un par de días. 4.4 Análisis de sensibilidad Para facilitar la comprensión de los resultados obtenidos sería conveniente observar como dichos resultados se verían afectados al cambiar ciertos parámetros, en concreto se van a realizar análisis de sensibilidad con la excentricidad, el semieje mayor y la inclinación de la órbita. 36 Capítulo 4. Método Montecarlo 4.4.1 Variación en excentricidad Como observamos en la Tabla 4.1 la órbita utilizada inicialmente es prácticamente circular debida a la excentricidad tan pequeña que tiene. Para ver cuál es el efecto de esta en la evolución de la incertidumbre se va a probar con excentricidades iguales a 0.3 y 0.6 respectivamente, obteniendo los resultados que se muestran en la Figura 4.8. Con estas gráficas se pueden extraer cierta información realmente interesante. Lo primero que es fácil de ver a simple vista es el hecho de que al aumentar la excentricidad en la órbita aumenta el crecimiento de la incertidumbre, pues tanto la suma de autovalores como las desviaciones típicas por separado son mayores, lo que significa que el elipsoide de confianza es mayor. Esto es debido a que al aumentar la excentricidad los cambios en velocidad y por ende en posición son más abruptos y por tanto hay una mayor dispersión de la incertidumbre. Hecho que se magnifica en el perigeo de la órbita, punto donde el satélite pasa a mayor velocidad, como se puede observar en la periodicidad de las gráficas de periodo igual al periodo orbital. Para facilitar la comprensión del efecto de la excentricidad en la incertidumbre se presenta en la Tabla 4.2 los valores de desviaciones típicas y suma de autovalores para un tiempo concreto de 20000s en los tres casos. Tabla 4.2 Variación en excentricidad (MC). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.02 0.063 122.63 e0=0.30.29 16.39 0.091 232.57 e0=0.60.62 33.88 0.24 1031.8 Con esta tabla es más fácil observar el aumento de la incertidumbre con el aumento de excentricidad en la órbita de partida. 4.4.2 Variación en semieje mayor Para observar el efecto que tendría en la incertidumbre un cambio en el tamaño de la órbita se procede a variar el semieje mayor de los iniciales 6878 km a 8000 km y 9000 km respectivamente obteniendo los resultados de la Figura 4.9 Estas gráficas nos permiten ver lo primero de todo que ahora tenemos el efecto contrario que con la excentricidad, pues al aumentar el valor del semieje mayor ahora el crecimiento de la incertidumbre es más lento. Este hecho es debido a que al aumentar el tamaño de la órbita estamos aumentando el periodo orbital y la distancia de perigeo, lo cual va a hacer que haya menor variación de velocidad y por tanto también en posición, haciendo que las dispersión de la incertidumbre sea menor conforme aumentemos dicho tamaño orbital. Para facilitar la comprensión del efecto de aumentar el semieje mayor en la incertidumbre se presenta en la Tabla 4.3 los valores de desviaciones típicas y suma de autovalores para un tiempo concreto de 20000s en los tres casos. Tabla 4.3 Variación en semieje mayor (MC). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.02 0.063 122.63 a0=8000km 0.26 9.75 0.07 89.46 a0=9000km 0.12 7.72 0.11 52.35 4.5 Test de normalidad de Henze-Zirkler 37 Con esta tabla queda mucho más claro el efecto de la disminución de la incertidumbre para un tiempo concreto con el aumento del semieje mayor. 4.4.3 Variación en inclinación El último parámetro que vamos a modificar para realizar un análisis de sensibilidad en la evolución de la incertidumbre es la inclinación, la cual en un principio es de 20 ◦ y la vamos a modificar a 40 ◦ y 60◦respectivamente, obteniendo los resultados que se pueden ver en la Figura 4.10. Ahora resulta curioso observar como al contrario que pasaba en los casos 4.4.1 y 4.4.2 el efecto no está tan claro. Se puede concluir por ende que la inclinación en un principio no va a modificar la incertidumbre en la órbita, pues no modifica la velocidad de paso del satélite en la órbita y por tanto no aumenta ni disminuye la dispersión de la incertidumbre. Esto se puede observar de manera más clara en la Tabla 4.4 en la que se compara para los tres casos mencionados y para un tiempo concreto de 20000s la evolución de la incertidumbre. Tabla 4.4 Variación en inclinación (MC). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.02 0.06 121.8 i0=40◦0.16 11.03 0.06 122.63 i0=60◦0.16 11.11 0.06 123.63 4.5 Test de normalidad de Henze-Zirkler El test de normalidad multivariable de Henze-Zirkler es una prueba estadística utilizada para evaluar si un conjunto de datos multivariados se ajusta a una distribución normal. Esta prueba es una extensión del test de normalidad univariada de Anderson-Darling y se utiliza para comprobar si las variables aleatorias en un conjunto de datos tienen una distribución normal conjunta. El test se basa en una distancia de divergencia entre la distribución empírica y la distribución normal teórica. Esta distancia es medida mediante una función de densidad de probabilidad que se ajusta a los datos de entrada y se compara con una función de densidad normal teórica. La hipótesis nula del test es que los datos siguen una distribución normal multivariada. Si el valor asociado con el test es menor que un nivel de significancia previamente establecido, entonces se rechaza la hipótesis nula y se concluye que los datos no siguen una distribución normal multivariada. Una explicación más en detalle se puede encontrar en [13]. Si aplicamos dicho test a nuestra órbita inicial y posteriormente a las distintas modificaciones realizadas en el análisis de sensibilidad observamos: Tabla 4.5 Resultados test de normalidad. Órbita inicial Tiempo test Variación en eTiempo test Variación en aTiempo test Variación en iTiempo test ≈23000 s e0=0.3≈11000 s a0=8000 km ≈41000 s i0=40◦≈23000 s e0=0.6≈5200 s a0=9000 km ≈49000 s i0=60◦≈23000 s Tabla en la cual se recoge el tiempo en el cual la incertidumbre sigue una distribución gaussiana, observando como para el caso en el que hacemos crecer la excentricidad de nuestra órbita la distribución deja de ser gaussiana mucho antes, lo que está en concordancia con lo que vimos en el Punto 4.4.1 en el cual observábamos que la incertidumbre crecía más rápidamente con el crecimiento de la excentricidad. Por otro lado al crecer el semieje mayor de la órbita, es decir, crece el tamaño de la órbita y también su periodo, el tiempo en el cual deja de ser una distribución normal aumenta hecho que 44 Capítulo 5. Método de Linealización Como se puede observar en la Figura 5.1 los resultados, pese a usar distintos métodos, son coherentes con los obtenidos en el Capítulo 4. Los elipsoides de confianza empiezan inicialmente como una esfera debido a la igualdad de las desviaciones típicas iniciales en los tres ejes y con el paso del tiempo su tamaño va aumentando. Este efecto se ve magnificado sobre todo en la dirección OY local pues, como ya explicamos en el Punto 4.3, esto es debido al efecto de la variación del periodo orbital en los puntos de la nube. Todos estos efectos que se observan visualmente en la Figura 5.1 se pueden entender mejor a través de las Figuras 5.2 y 5.3 donde se muestra la evolución temporal de la matriz de covarianza. Figura 5.2 Comparación desviaciones típicas (LIN). Figura 5.3 Suma de autovalores (LIN). 5.4 Análisis de sensibilidad 45 5.4 Análisis de sensibilidad Al igual que se hizo para el método Montecarlo, se va a analizar como cambia la evolución de la incertidumbre si cambiamos la órbita inicial de referencia. En concreto se varía la excentricidad, el semieje mayor y la inclinación. 5.4.1 Variación en excentricidad Para poder observar los cambios que se producen en la evolución de la incertidumbre si cambiamos la excentricidad inicial de nuestra órbita de referencia, se procede a hacer un análisis con excentricidades iguales a 0.3 y 0.6, obteniendo los resultados que se muestran en la Figura 5.5. De nuevo, estos resultados son coherentes con los obtenidos en el Capítulo 4. La incertidumbre crece con la excentricidad debido a que al aumentar esta los cambios en velocidad y por ende en posición son mayores. Estos efectos se acentúan por tanto en el perigeo de la órbita, punto de mayor velocidad de paso, como se observa en la periodicidad de la Figura 5.5 de periodo igual al periodo orbital. Para facilitar la comprensión de dichos efectos se presenta en la Tabla 5.1 los valores de desviaciones típicas y suma de autovalores para un tiempo concreto de 20000s en los tres casos. Tabla 5.1 Variación en excentricidad (LIN). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 11.02 0.061 123.23 e0=0.30.29 16.19 0.084 232.57 e0=0.60.63 33.88 0.17 1133.6 5.4.2 Variación en semieje mayor En el caso de una variación en el semieje mayor de la órbita y por ende en su tamaño, se procede a hacer un estudio variando el semieje mayor de los iniciales 6878 km a 8000 km y 9000 km respectivamente obteniendo los resultados mostrados en la Figura 5.6. Al igual que en el Punto 4.4.2 observamos que, al aumentar el tamaño de la órbita, la incertidumbre tarda más en crecer, pues, ahora los cambios en velocidad y posición son menores y en consecuencia la dispersión de la incertidumbre es menor. Como en el Punto 5.4.1, a continuación se presenta la Tabla 5.2 con los valores de las desviaciones típicas y suma de autovalores para un tiempo de 20000s en los tres casos. Tabla 5.2 Variación en semieje mayor (LIN). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 11.02 0.06 123.23 a0=8000km 0.16 9.45 0.07 89.47 a0=9000km 0.11 7.43 0.11 55.25 5.4.3 Variación en inclinación El último análisis de sensibilidad que vamos a hacer es al cambiar la inclinación de los iniciales 20 ◦ a 40◦y 60◦respectivamente, obteniendo los resultados de la Figura 5.4. En estas gráficas se observa el mismo efecto que se observó en 4.4.3 y es que la inclinación no va a tener apenas efecto en la incertidumbre. La diferencia radica es que con el método de Linealización no se captura ese pequeño efecto que si se veía en la Figura 4.10 con el método Montecarlo. 46 Capítulo 5. Método de Linealización A continuación se presenta en la Tabla 5.3 la comparativa de la incertidumbre en la órbita para los tres casos comentados y para un tiempo concreto de 20000s. Tabla 5.3 Variación en inclinación (LIN). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 11.02 0.06 123.23 i0=40◦0.15 11.02 0.06 123.23 i0=60◦0.15 11.02 0.06 123.23 Como vemos no hay diferencias en la incertidumbre con la inclinación si hacemos el análisis con el método de Linealización. Figura 5.4 Variación de σy suma de autovalores con i0(LIN). 5.4 Análisis de sensibilidad 47 Figura 5.5 Variación de σy suma de autovalores con e0(LIN). 48 Capítulo 5. Método de Linealización Figura 5.6 Variación de σy suma de autovalores con a0(LIN). 6 Unscented Transform (UT) 6.1 Definición El método de Unscented Transform (UT) es una técnica utilizada en la teoría de estimación estadística para aproximar la propagación de incertidumbre a través de funciones no lineales. Este método se utiliza principalmente en la estimación de estados en sistemas dinámicos no lineales. Su objetivo es obtener una representación precisa de la distribución de probabilidad de un vector de estados después de que ha pasado a través de una función no lineal. Esto es especialmente útil cuando la función no lineal no puede ser aproximada de manera precisa mediante una expansión de primer orden como la serie de Taylor. El método de UT utiliza un conjunto de puntos llamados sigma points para representar la distribución de probabilidad del estado. Estos puntos se eligen de manera estratégica para capturar los momentos de la distribución, como la media y la covarianza. La idea clave detrás de este método es propagar estos puntos a través de la función no lineal para obtener una estimación de la distribución resultante. La principal ventaja de la UT es que puede capturar de manera más precisa la propagación de la incertidumbre a través de funciones no lineales en comparación con los métodos de aproximación lineal tradicionales. Sin embargo, también tiene algunas limitaciones, como la necesidad de seleccionar adecuadamente los sigma points y la complejidad computacional asociada con el cálculo de las transformaciones no lineales. 6.2 Desarrollo del método Una vez definidos los fundamentos del método vamos a presentar la matemática que rodea a la UT. Expresiones que se han extraído de [2]. Para la aplicación del método el primer paso es definir cuantos sigma points necesitamos y cuales son. Con respecto a la cantidad de sigma points en [ 2 ] se justifica que el número más aceptado en la literatura es 2n+1 , siendo n las dimensiones de nuestro problema. Así un número más elevado requeriría un excesivo coste computacional. Por otro lado, para la definición de estos sigma points ( S ) hay distintos métodos, sin embargo nosotros nos quedaremos con el que se define en [2]: 49 50 Capítulo 6. Unscented Transform (UT) S[0]=µ0(6.1) S[i]=µ0+p(n+κ)Σ0[i](6.2) S[i+n]=µ0−p(n+κ)Σ0[i](6.3) para i∈{1,···,n}. En estas expresiones µ0 representa la media y Σ0 la matriz de covarianza de nuestra distribución inicial. Recordemos que esta distribución inicial como dijimos en el Punto 4.2 es la misma que para el método Montecarlo y de Linealización. Ahora que tenemos los 2n+1sigma points de la UT, se procede a propagarlos para un instante t cualquiera a través de la resolución de las ecuaciones de los dos cuerpos definidas en (2.6). Una vez propagados los puntos simplemente queda definir la nueva media y covarianza en el instante t: µ(t) = 2n ∑ i=0 wiSi(6.4) Σ(t) = 2n ∑ i=0 wi(Si−µ(t))(Si−µ(t))T(6.5) donde wi son los pesos que se usan en la UT asociados cada uno con un sigma points y definidos como: w[0]=κ n+κ(6.6) w[i]=1 2(n+κ)(6.7) w[i+n]=1 2(n+κ)(6.8) con i∈{1,···,n} y con κ un parámetro libre de la UT que en la literatura típicamente se escoge como κ=n−3. Con estos pesos definidos ya tenemos la covarianza en cualquier instante de tiempo, expresión (6.5), y por ende la evolución de la incertidumbre en la órbita. 6.3 Resultados Una vez que conocemos ya tanto los fundamentos del método como la matemática tras él, pasamos a su implementación numérica para la misma órbita de referencia usada para los otros dos métodos, Tabla 4.1. Así, podemos obtener de nuevo la evolución de la incertidumbre tras una revolución completa en la órbita pero con otro método. En la Figura 6.1 se puede ver como se obtienen de nuevo resultados coherentes con los métodos anteriores, y es que al inicio el elipsoide de confianza vuelve a ser prácticamente una esfera y con el paso de tamaño aumenta su tamaño. Alargándose sobre todo en la dirección OY local debido a ese efecto de la disminución y aumento del periodo orbital en función del semieje como ya explicamos en el Punto 4.3. Este efecto se puede vislumbrar más claramente en las Figuras 6.2 y 6.3 donde se muestra la evolución temporal de la matriz de covarianza, es decir, de la incertidumbre en la órbita. 6.4 Análisis de sensibilidad 51 Figura 6.1 Elipsoides de confianza en una revolución (UT). Figura 6.2 Comparación desviaciones típicas (UT). 6.4 Análisis de sensibilidad En este punto se va a proceder a examinar como varía la evolución de la incertidumbre si cambiamos ciertos parámetros de la órbita inicial de referencia como son la excentricidad, el semieje mayor o la inclinación. Además, se realizará un análisis para observar cuan relevante resulta el parámetro κ en la UT al ser este un parámetro libre a elegir. 52 Capítulo 6. Unscented Transform (UT) Figura 6.3 Suma de autovalores (UT). 6.4.1 Variación en excentricidad Para observar cual es el efecto de la excentricidad en la evolución de la incertidumbre para este nuevo método, se procede a variar la excentricidad de la inicial 0.01 a 0.3 y 0.6 respectivamente con lo que obtenemos los resultados de la Figura 6.5. Al igual que en los anteriores métodos se observa una tendencia creciente con la excentricidad debido al efecto que tiene esta en la velocidad de paso y por tanto en la dispersión de la incertidumbre. Efectos explicados ya en los anteriores métodos, Punto 4.4.1. En la Tabla 6.1 se puede ver la evolución de la incertidumbre para un tiempo de 20000s en los tres casos estudiados. Observándose la tendencia claramente creciente mencionada anteriormente. Tabla 6.1 Variación en excentricidad (UT). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.1 0.064 123.22 e0=0.30.29 16.19 0.09 262.57 e0=0.60.63 33.66 0.24 1133.7 6.4.2 Variación en semieje mayor Para el caso del semieje mayor se procede a hacer un análisis de los cambios en la incertidumbre variando el semieje mayor a 8000 km y 9000 km respectivamente. El resultado de este se puede ver en la Figura 6.6. Ahora la tendencia es la contraria que para la excentricidad, al aumentar el semieje mayor la incertidumbre crece más lentamente. Esto es debido al aumento del tamaño de la órbita y por ende a la disminución de la velocidad de paso y de los cambios en posición. Explicado este efecto ya en el Punto 4.4.2. Para clarificar el efecto mencionado se presenta en la Tabla 6.2 los datos numéricos del análisis realizado para los tres casos en un tiempo de 20000s. 6.4 Análisis de sensibilidad 53 Tabla 6.2 Variación en semieje mayor (UT). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.1 0.06 123.22 a0=8000km 0.16 9.45 0.07 89.47 a0=9000km 0.11 7.43 0.11 55.25 6.4.3 Variación con inclinación El último parámetro que se va a modificar es la inclinación, la cual se va a variar de los iniciales 20◦a 40◦y 60◦, obteniendo los resultados mostrados en la Figura 6.7. Con la UT se constata que, al igual que con Montecarlo y la Linealización, el efecto de la inclinación es limitado en la evolución de la incertidumbre. En la Tabla 6.3 se presenta la comparativa de la incertidumbre en la órbita para los tres casos comentados y para un tiempo concreto de 20000s. Tabla 6.3 Variación en inclinación (UT). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.16 11.1 0.06 123.23 i0=40◦0.16 11.1 0.06 123.22 i0=60◦0.16 11.1 0.06 123.22 6.4.4 Variación de κ Al inicio del capítulo, en concreto en el Punto 6.2, se mencionó de la existencia en la UT del parámetro κ , un parámetro libre a elegir por el usuario que afecta tanto a los pesos como a la definición de los sigma points. En la literatura normalmente κ se toma como n -3 pues, como se explica en [ 2 ], este valor minimiza el desajuste de momento de cuarto orden. Sin embargo, este parámetro puede tomar el valor que nosotros consideremos y es por ello que vamos a analizar su efecto en la UT. En nuestro caso como tenemos n =6, κ según la literatura valdría 3 pero para la realización de este análisis se ha probado con valores de κde 1 a 6 obteniendo los resultados de la Figura 6.4. En la Figura 6.4 se observa como el valor de κ=3 se sitúa en medio mientras que los valores menores se adelantan y los mayores se retrasan. De todas formas, el efecto a gran escala de κ parece que no es demasiado relevante para los resultados obtenidos con el método, así que se opta por un κ=3siguiendo así con lo utilizado en la literatura. 60 Capítulo 7. Método de los Elementos Orbitales Equinocciales donde a,e,ω,Ω,i,θ son los Elementos Orbitales Keplerianos ( EOK ) ya estudiados en el Capítulo 2. A su vez, para el uso del método también va a ser preciso obtener los EOK a partir de los EOEM para con ellos pasar de nuevo al SCC: a=p 1−f2−g2(7.7) e=pf2+g2(7.8) i=2arctanph2+k2(7.9) Ω=arctank h(7.10) ω=arctang f−Ω(7.11) θ=L−ω−Ω(7.12) (7.13) Una vez definidas estas transformaciones ya estamos en disposición de hacer uso del método. Lo primero es que, al igual que con el resto de métodos, inicialmente tenemos nuestra Σ y µ de la distribución normal inicial de nuestra órbita de referencia que ya quedaron definidas en el Punto 4.2. Ahora, y haciendo uso de lo aprendido en el Capítulo 6 con la UT, vamos a pasar esta Σ y µ al espacio de los EOEM . Para ello obtenemos los sigma points de la distribución inicial en cartesianas y lo transformamos con las expresiones (7.1) a los EOEM . Con las fórmulas (6.4) de la UT calculamos la nueva Σyµen el espacio de los EOEM. Una vez tenemos nuestra Σ y µ iniciales en EOEM y los sigma points pasamos a la propagación. Como adelantamos en el Punto 7.1 todos los elementos van a permanecer constantes exceptuando L. Así, la propagación queda: p(t) = p0(7.14) f(t) = f0(7.15) g(t) = g0(7.16) h(t) = h0(7.17) k(t) = k0(7.18) L(t) = L0+n0△t(7.19) donde n0=qµ a3 0 está referido al istante inicial t0. Con los sigma points ya propagados el cálculo de la nueva µ y Σ en el espacio de los EOEM es inmediato pues, simplemente hay que volver a usar las expresiones (6.4) de la UT para obtener esta nueva media y matriz de covarianza propagadas. Ahora, para poder comparar este método con los otros procedemos a volver al SCC . Para ello, se vuelve a proceder de la misma forma en la que procedió incialmente para pasar al espacio de los EOEM pero al revés. Es decir, los sigma points propagados en el espacio de los EOEM se pasan al SCC a través de las expresiones (7.7) y con las mimas expresiones de la UT que antes, las definidas en (6.4), se obtiene la nueva media y matriz de covarianza pero ahora en el SCC. 7.3 Resultados 61 7.3 Resultados Con este último método que vamos a usar para el estudio de la incertidumbre en la órbita se va a proceder de manera completamente análoga que con los anteriores métodos. Primero se implementa numéricamente el método, como quedó ya explicado en el Punto 7.2, y posteriormente se extraen los resultados para una revolución completa de la órbita. Esto se muestra en la Figura 7.1. Figura 7.1 Elipsoides de confianza en una revolución (EOEM). Así, como se puede observar en la figura el resultado es igual que los métodos usados previamente. De nuevo, se observa un crecimiento de la incertidumbre, el cual se acentúa en dirección OY local. Este efecto quedó ya explicado en el Punto 4.3. Ahora, para que quede más clara esta evolución, se va a proceder de nuevo a mostrar la evolución de la matriz de covarianza, esto es, de la suma de autovalores y desviaciones típicas. Estos resultados se muestran en las Figuras 7.2 y 7.3. La tendencia creciente está en concordancia con los resultados obtenidos con los anteriores métodos. Ahora la diferencia es que esta componente periódica, al menos para la órbita inicial, no queda tan clara. Es decir, los EOEM parecen no captar del todo bien estos aumentos de incertidumbre al pasar por el perigeo de la órbita, que recordemos, es debido al aumento de los cambios en velocidad y posición en dicho punto. Lo que provoca un aumento de la dispersión de la incertidumbre. Esta diferencia entre métodos quedará más clara en el posterior análisis de sensibilidad y en el Capítulo 8 donde se compararán todos ellos. 7.4 Análisis de sensibilidad Al hilo de las diferencias observadas en el Punto 7.3 en la suma de autovalores y desviaciones típicas con respecto a los demás métodos, resulta de lo más interesante realizar el análisis de sensibilidad que se hizo en los anteriores métodos en este. Este análisis consiste, como en los demás, en variar la excentricidad, semieje mayor e inclinación y observar el cambio en la evolución de la incertidumbre en la órbita. 62 Capítulo 7. Método de los Elementos Orbitales Equinocciales Figura 7.2 Comparación desviaciones típicas (EOEM). Figura 7.3 Suma de autovalores (EOEM). 7.4.1 Variación en excentricidad Para el caso del efecto de la excentricidad en la evolución de la incertidumbre se va a variar la excentricidad de la inicial a 0.3 y 0.6 respectivamente. Con estos cambios se obtienen los resultados de la Figura 7.4. En esta figura se muestra que, como cabía esperar, el efecto sigue en concordancia con los anteriores métodos, y es que la tendencia es claramente creciente con la excentricidad debido al aumento de la dispersión de la incertidumbre por el aumento de la velocidad de paso con el aumento de la excentricidad. El cambio más significativo es que con el aumento de la excentricidad si se percibe el efecto periódico de los resultados, que comentamos en el Punto 7.3, de manera más clara. Esto es debido a 7.4 Análisis de sensibilidad 63 que al aumentar la excentricidad la velocidad de paso por el perigeo lo hace también, provocando un aumento mayor en la incertidumbre. En la Tabla 7.1 se muestra dicha evolución de la incertidumbre para el caso concreto de una propagación de 20000s. Tabla 7.1 Variación en excentricidad (EOEM). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 12.68 0.069 160.91 e0=0.30.15 17.81 0.07 317.4 e0=0.60.17 31.14 0.15 969.9 7.4.2 Variación en semieje mayor Para el caso del análisis de sensibilidad con el semieje mayor se procede a variar este de los iniciales 6878 km a 8000 km y 9000 km respectivamente, obteniendo los resultados de la Figura 7.5. En este caso la tendencia principal de decrecimiento de la incertidumbre con el crecimiento del semieje mayor vuelve a coincidir con los anteriores métodos. Pero ahora, al contrario que pasaba con el caso de la excentricidad, esta componente periódica no se acentúa sino que es prácticamente igual que la inicial. Para hacer más claro el efecto del semieje mayor en la incertidumbre para el caso de los EOEM se presenta la Tabla 7.2 donde se muestra los resultados del análisis para un tiempo concreto de 20000s. Tabla 7.2 Variación en semieje mayor (EOEM). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 12.68 0.069 160.91 a0=8000km 0.17 10.41 0.059 108.42 a0=9000km 0.16 8.97 0.077 80.61 7.4.3 Variación en inclinación Ahora, para el caso del análisis de sensibilidad con la inclinación de la órbita se procede a variar la inclinación de la inicial a 40 ◦ y 60 ◦ respectivamente, con lo que se obtienen los resultados mostrados en la Figura 7.6. Al igual que pasaba con los anteriores métodos, debido a que el cambio de inclinación no provoca un cambio en la velocidad de paso en la órbita y por ende en la dispersión de la incertidumbre, los resultados para los tres casos son prácticamente iguales. En la Tabla 7.3 se muestran los resultados concretos para el caso de una propagación de 20000s. Observando como, efectivamente, los cambios son mínimos para los tres casos. Tabla 7.3 Variación en inclinación (EOEM). σx(km) σy(km) σz(km) Suma autovalores Órbita inicial 0.15 12.68 0.069 160.91 i0=40◦0.15 12.68 0.069 160.91 i0=60◦0.15 12.68 0.069 160.91 64 Capítulo 7. Método de los Elementos Orbitales Equinocciales Figura 7.4 Variación de σy suma de autovalores con e0(EOEM). 7.4 Análisis de sensibilidad 65 Figura 7.5 Variación de σy suma de autovalores con a0(EOEM). 66 Capítulo 7. Método de los Elementos Orbitales Equinocciales Figura 7.6 Variación de σy suma de autovalores con i0(EOEM). 8 Comparativa de métodos Este capítulo viene a cerrar la parte del estudio de métodos para el análisis de la incertidumbre en órbitas. Así, se procederá a estudiar y analizar las ventajas e inconvenientes que tienen la utilización de cada uno de los métodos y se realizará comparativas de los resultados obtenidos con ellos. Para hacer esta comparativa, se tomará el método Montecarlo como “realidad”, pues, con suficientes puntos para la simulación es el que mejor la captura, ya que resuelve directamente las ecuaciones. 8.1 Ventajas e inconvenientes de los métodos Como avanzamos con anterioridad, lo primero a realizar para la elección de los mejores métodos es, sin duda, conocer en detalle cuáles son las ventajas e inconvenientes principales que tienen cada uno de estos métodos. 8.1.1 Método Montecarlo En cuanto a las ventajas que presenta este método destacamos: •Capacidad para modelar y propagar incertidumbres complejas y no lineales. •Puede considerar fuentes de incertidumbre tanto aleatorias como sistemáticas. • Permite obtener distribuciones de probabilidad para los resultados, lo que facilita el análisis estadístico. •Puede adaptarse a diferentes niveles de detalle y complejidad del modelo orbital. En contraposición sus desventajas principales son: • Requiere un gran número de muestras para obtener resultados precisos, lo que puede llevar a un alto coste computacional. • La precisión y calidad de los resultados dependen en gran medida de la representatividad de las muestras generadas. • Puede resultar difícil de implementar y validar correctamente debido a la complejidad asociada a la generación de muestras aleatorias. 8.1.2 Método de Linealización Las ventajas principales que presenta el uso de la Linealización son: •Método apropiado para sistemas donde las no linealidades son pequeñas. • Es computacionalmente eficiente en comparación con los métodos no lineales más complejos. 67 68 Capítulo 8. Comparativa de métodos • Puede proporcionar una solución analítica aproximada, lo que facilita el análisis matemático. Los inconvenientes de su uso por otro lado resultan ser: • No tiene en cuenta las no linealidades importantes, lo que puede introducir errores significativos con el paso del tiempo. •No proporciona una medida directa de la distribución de probabilidad de los resultados. •Puede requerir conocimiento detallado del modelo linealizado y su derivación. 8.1.3 Unscented Transform Las principales ventajas de este método son: • Aproxima la propagación de incertidumbre en sistemas no lineales sin requerir una Linealización explícita. •Proporciona una estimación precisa de la media y covarianza de la distribución resultante. •Menos susceptible a errores de muestreo en comparación con el método Montecarlo. Por otro lado, los inconvenientes que presenta son: • Requiere una mayor cantidad de evaluaciones del modelo en comparación con el método de Linealización, lo que puede aumentar la carga computacional. •No proporciona una descripción completa de la distribución de probabilidad resultante. • Puede ser más complejo de implementar y entender en comparación con el método de Linealización. 8.1.4 Elementos Orbitales Equinocciales Modificados En cuanto a las ventajas de trabajar en el espacio de los EOEM destacamos: • Con los EOEM se puede llegar a propagar la PDF del estado completo usando una suma gaussiana unidimensional. • Cinco de las seis coordenadas evolucionan linealmente bajo el problema de los dos cuerpos. •La incertidumbre crece principalmente solo a lo largo de una de las seis coordenadas. Las desventajas que presenta en contraposición son: • Debido a la forma de propagación del método puede no capturar tan bien la incertidumbre como los demás métodos. 8.2 Comparativa de resultados Una vez hemos analizado cuáles son las ventajas e inconvenientes de cada método y podemos tener una idea más clara de ellos a la hora de hacer una elección, vamos a pasar a comparar los resultados que se obtienen con cada uno de ellos. 8.2 Comparativa de resultados 69 Para ello lo primero que es importante definir es la órbita de referencia, pues, para que la comparativa sea coherente esta debe ser la misma. En concreto, usaremos la que hemos estado utilizando hasta ahora y es la definida por sus elementos orbitales en la Tabla 4.1, pero que volvemos a indicar en la Tabla 8.1 para que quede más claro en este capítulo. Tabla 8.1 Elementos orbitales iniciales. a e i ωΩM 6878 km 0.01 20◦90◦45◦0◦ Como la órbita es la misma que hemos estado usando hasta ahora para estudiar los métodos individualmente, solo queda superponer los resultados para ver cuan de semejantes o diferentes son del método Montecarlo. Estos resultados se muestran en la Figura 8.1 donde se pueden observar para una simulación de medio día la evolución de la incertidumbre, esto es, de la suma de autovalores y desviaciones típicas con los cuatro métodos. Con estos resultados se pueden extraer cierta información interesante sobre los métodos. Lo primero que resulta claro es que el método de los EOEM resulta ser el peor debido a las diferencias con los demás métodos. Esto se puede explicar por la forma en la que se propaga con los EOEM , y es que, como se dijo en el Punto 8.1.4, la incertidumbre crece principalmente solo a lo largo de una de las seis coordenadas mientras que las otras cinco permanecen constantes. Esto hace que el método no capture bien la incertidumbre. Con los otros métodos, la Linealización y la UT, no resulta tan claro pues, se puede observar como al menos para cortos periodos de tiempos ambos métodos se difieren del Montecarlo de una manera prácticamente similar. Para poder ver mejor las diferencias entre estos dos métodos se procede a analizar la evolución para mayor tiempo, en concreto se muestra en las Tablas 8.2, 8.3 y 8.4 los resultados para 20000s, 7 días y 30 días. Con estos resultados se puede concluir que el método de Linealización, como se comentó en el Punto 8.1.2, pierde exactitud con el tiempo al no tener en cuenta la naturaleza no lineal de las ecuaciones del movimiento. Tabla 8.2 Evolución de la incertidumbre en 20000s con tres métodos. σx(km) σy(km) σz(km) Suma autovalores Tiempo (s) MC 0.16 11.02 0.063 122.63 3.82 LIN 0.15 11.02 0.062 123.23 0.034 UT 0.16 11.1 0.064 123.23 0.062 Tabla 8.3 Evolución de la incertidumbre en 7 días con tres métodos. σx(km) σy(km) σz(km) Suma autovalores Tiempo (s) MC 10.33 324.69 0.15 105530 107.81 LIN 0.15 333.51 0.058 111230 0.42 UT 11.81 333.11 0.15 111100 1.43 Tabla 8.4 Evolución de la incertidumbre en 30 días con tres métodos. σx(km) σy(km) σz(km) Suma autovalores Tiempo (s) MC 199.69 1393.7 0.13 1982200 433 LIN 0.13 1431.2 0.09 2048300 1.63 UT 214.87 1399.3 0.13 2004200 5.38 76 Capítulo 9. Efecto de las perturbaciones En cambio si la órbita es retrógrada ( i∈[90◦,180◦) ) entonces ˙ Ω>0 por lo que la línea de nodos avanzaría hacia el Este ganando ascensión recta. Para órbitas polares ( i=90◦ ) la línea nodal se mantendrá fija. •Retroceso/avance del perigeo: En la ecuación (9.2) podemos ver que la variación del argumento del perigeo se anula para determinados valores de la inclinación: ˙ ω=0⇔5 2sin2(i)−2=0→sin(i) = r4 5→i=arcsin r4 5!=63.4◦ 116.4◦(9.4) Por lo tanto ˙ ω>0 , cuando i∈[0◦,63.4◦) o i∈(116.6◦,180◦] , indica que el perigeo avanza en la dirección del movimiento, fenómeno conocido como Avance del perigeo. El efecto contrario, Retroceso del perigeo, se presenta cuando ˙ ω<0 , es decir cuando i∈(63.4◦,116.6◦) . Notemos que hay dos inclinaciones críticas en las que la línea de ápsides no gira: i=63.4◦ei=116.6◦. Resultados numéricos Llegados a este punto tenemos que aclarar que la implementación computacional del efecto del J2 en la órbita no se hará siguiendo las expresiones (9.1), (9.2) y (9.3) sino que directamente se resolverán numéricamente las ecuaciones del problema de los dos cuerpos añadiendo el efecto del J2 en ellas. Esto es debido a que dichas expresiones solo son válidas para elementos medios y nosotros trabajaremos con elementos orbitales osculantes 1 . Así, las ecuaciones del movimiento en cartesianas quedan: ¨ r=−µr r3+aJ2 pert (9.5) con aJ2 pert =       µx r3h3 2J2RT r25z r2−1i µy r3h3 2J2RT r25z r2−1i µz r3h3 2J2RT r25z r2−3i        (9.6) donde J2es 0.00108248, valor extraido de [11] y r=px2+y2+z2. Expresiones del efecto del J2en cartesianas extraídas de [1]. Con la ecuación (9.5) ya completa con el efecto del J2 se procede a analizar su efecto en nuestra órbita de referencia que recordemos viene dada por sus elementos orbitales iniciales en la Tabla 4.1. Así, los resultados obtenidos son los que se muestran en las Figuras 9.1 , 9.2 y 9.3. En las Figuras 9.2 y 9.3 se muestra por un lado el efecto mencionado en el Punto 9.1.1, el efecto secular decreciente de Ω por ser una órbita directa y a su vez un efecto periódico de corto periodo. Por otro lado, en el caso de ω se muestra un comportamiento secular creciente debido a que la inclinación de la órbita es menor que la crítica de 63.4 ◦ y también la componente periódica de corto periodo. El hecho de que el efecto sea tanto secular como periódico está en concordancia con lo mostrado en la Figura 2.11. En la Figura 9.1 sin embargo, se muestra de una manera más visual el efecto que tiene estos cambios en los elementos orbitales en la propia órbita, observando un efecto oscilante en ella. 1 Elementos correspondientes a la órbita definida por el vector posición y velocidad en un determinado instante de tiempo. 9.1 Efecto de las perturbaciones en la órbita 77 Figura 9.1 Efecto del J2en la órbita de referencia. 9.1.2 Efecto de la resistencia aerodinámica Con respecto a la resistencia aerodinámica, las variaciones que provoca en los elementos orbitales corresponden con su efecto de circularizar la órbita. Esto es, la disminución de la excentricidad y del semieje mayor. Efectos explicados ya en el Punto 2.3.3. Resultados numéricos Para la implementación del efecto de la resistencia aerodinámica, al igual que se hizo con el J2 , se va a resolver numéricamente las ecuaciones de los dos cuerpos en cartesianas introduciendo en ellas los términos de la resistencia aerodinámica quedando: ¨ r=−µr r3+aDrag pert (9.7) con 78 Capítulo 9. Efecto de las perturbaciones Figura 9.2 Evolución de Ωcon el tiempo. Figura 9.3 Evolcución de ωcon el tiempo. aDrag pert =   −1 2 1 Cbρvr(vx+ω⊕y) −1 2 1 Cbρvr(vy+ω⊕x) −1 2 1 Cbρvrvz   (9.8) donde ω⊕ es la velocidad angular de la Tierra , Cb el coeficiente balístico que escogeremos como 9.1 Efecto de las perturbaciones en la órbita 79 6.69·10−8 por semejanza con otros satélites LEO, vx,vy,vz las componentes de la velocidad, vr el módulo de la velocidad relativa, definida esta como: vrel =v−vatm (9.9) con vatm =ω⊕×r(9.10) Mientras que ρ es la densidad, calculada esta a través de un modelo exponencial como el siguiente: ρ=ρ0exp−hellp −h0 H(9.11) donde hellp es la altura del satélite y ρ0 , h0 y H son tres parámetros tabulados en función de hellp como se muestra en la Figura 9.4. Tabla extraida integramente de [12]. Figura 9.4 Tabla para la obtención de los parámetros del modelo exponencial de densidad [12]. Ahora, con las ecuaciones del movimiento ya perfectamente completas y definidas se procede a analizar su efecto en nuestra órbita de referencia que, al igual que en 9.1.1, queda definida por sus elementos orbitales iniciales en la Tabla 4.1. Los resultados obtenidos se muestran en las Figuras 9.5, 9.6 y 9.7. En las Figuras 9.6 y 9.7 se muestra la evolución de la excentricidad y del semieje mayor que como avanzamos son ambos decrecientes. Tanto el semieje mayor como la excentricidad tienen por un lado una componente secular decreciente y por otro lado una componente de corto periodo. Efectos en concordancia de nuevo con lo mostrado en la Figura 2.11. En cuanto a la Figura 9.5 se muestra de una manera más visual el efecto de la resistencia aerodinámica, la circularización de la órbita. 80 Capítulo 9. Efecto de las perturbaciones Figura 9.5 Efecto de la resistencia aerodinámica en la órbita de referencia. Figura 9.6 Evolución de acon el tiempo. 9.2 Efecto de las perturbaciones en la incertidumbre En este apartado vamos a proceder a analizar cual es el efecto de las perturbaciones, tanto del J2 como de la resistencia aerodinámica, en la evolución de la incertidumbre en la órbita. Hasta ahora, hemos analizado por varios métodos dicha evolución pero no hemos tenido en cuenta perturbaciones 9.2 Efecto de las perturbaciones en la incertidumbre 81 Figura 9.7 Evolución de econ el tiempo. de ningún tipo. Para este análisis, en vista de los resultados obtenidos en el Capítulo 8, se va a utilizar el método Montecarlo y la UT para la propagación de la incertidumbre. Esto es debido a que la UT es el método más interesante para su uso, mientras que el Montecarlo lo usaremos como reflejo de la realidad. 9.2.1 Efecto del J2 Para estudiar el efecto que tiene la perturbación del J2 en la evolución de la incertidumbre con Montecarlo o con la UT, simplemente hay que introducir en las ecuaciones del movimiento las expresiones en cartesianas de la aceleración producida por el efecto del J2 . Expresiones que encontramos en (9.6). De nuevo se utilizará como referencia la órbita de la Tabla 4.1 y la media y matriz de covarianza inicial del Capítulo 2. Así, haciendo uso de ambos métodos obtenemos los resultados de la Figura 9.9. En esta figura, como ya avanzamos en el Punto 9.2, se ha tomado como “realidad” los resultados obtenidos en el método Montecarlo pues, como se explicó en el Capítulo 8, es este método el que más fiel refleja la realidad a resolver directamente las ecuaciones del movimiento para una nube de puntos inicial. Considerando esto, se observa como al principio el efecto es despreciable, pero con el paso del tiempo el Montecarlo captura un aumento de la incertidumbre, mientras que la UT no captura tan bien este efecto debido a que no deja de ser una aproximación a través de los sigma points. La razón de este aumento en la incertidumbre es debido al efecto que produce el J2 en la órbita, generando cambios en la forma y la orientación. Esto, como vimos en los análisis de sensibilidad realizados para todos los métodos, provocará cambios en la evolución de la incertidumbre de las órbitas, en este caso el aumento de ella. 82 Capítulo 9. Efecto de las perturbaciones 9.2.2 Efecto de la resistencia aerodinámica Para el análisis de la resistencia aerodinámica en la evolución de la incertidumbre en la órbita se procede de manera similar que para el caso del J2 . Esto es, se procede a añadir, como se explicó en el Punto 9.1.2, las expresiones en cartesianas de la aceleración producida por la perturbación de la resistencia aerodinámica a las ecuaciones del movimiento. Para el análisis, se utiliza una órbita de referencia con los mismos elementos orbitales iniciales que los de la Tabla 4.1 exceptuando la excentricidad que elegiremos de 0.008 y un semieje mayor de 6600 km. Esta elección no es arbitraria, sino que se elige una altura de órbita muy baja, en torno a 200 km, para que el efecto de esta perturbación sea mayor. Así como una excentricidad que hace prácticamente circular a la órbita para que el efecto se pueda observar también durante más tiempo. Con la órbita ya definida y las ecuaciones del movimiento completas con el término de la resistencia aerodinámica, se procede a realizar el análisis con ambos métodos obteniendo los resusltados de la Figura 9.10. En esta figura se muestra como la resistencia provoca ahora una disminución de la incertidumbre, con un efecto mucho mayor al que tenía el J2 . Esto es debido a esta altura tan baja elegida de la órbita. Así pues, la diferencia entre Montecarlo y la UT sigue siendo que la UT al ser una aproximación no captura del todo bien este efecto, pero se ve también ese comportamiento decreciente de la incertidumbre. La razón por la que disminuye la incertidumbre de nuevo está en concordancia con los análisis de sensibilidad realizados hasta ahora. Y es que al disminuir la excentricidad la incertidumbre baja más rápidamente que aumenta por la disminución del semieje mayor. 9.2.3 Análisis de sensibilidad Para que quede más claro el efecto de la resistencia aerodinámica se va a hacer un estudio de sensibilidad donde se va a proceder a variar la altura de la órbita definida en el Punto 9.2.2. En la Figura 9.8 se muestran los resultados para una h =200 km, 300 km y 400 km respectivamente. En esta figura se puede observar el efecto claro de disminución de la incertidumbre, tanto suma de autovalores como desviaciones típicas, con la altura. Este efecto es debido a que al aumentar la altura, el efecto de la resistencia aerodinámica es menor y por tanto no produce tantos cambios en posición ni velocidad, no hay tanta dispersión de la incertidumbre. También se puede observar como al cabo de un tiempo, unas horas, el satélite a una altura de 200 km empezaría a caer, pues hay un crecimiento brusco de la suma de autovalores y desviaciones que lo indica. 9.2 Efecto de las perturbaciones en la incertidumbre 83 Figura 9.8 Efecto de la altura en σy suma autovalores. 84 Capítulo 9. Efecto de las perturbaciones Figura 9.9 Efecto del J2 en σy suma autovalores en MC y UT. 9.2 Efecto de las perturbaciones en la incertidumbre 85 Figura 9.10 Efecto de la resistencia aerodinámica en σy suma autovalores en MC y UT. 92 Índice de Figuras 6.1 Elipsoides de confianza en una revolución (UT) 51 6.2 Comparación desviaciones típicas (UT) 51 6.3 Suma de autovalores (UT) 52 6.4 Variación de σy suma de autovalores con κ54 6.5 Variación de σy suma de autovalores con e0(UT) 55 6.6 Variación de σy suma de autovalores con a0(UT) 56 6.7 Variación de σy suma de autovalores con i0(UT) 57 7.1 Elipsoides de confianza en una revolución (EOEM) 61 7.2 Comparación desviaciones típicas (EOEM) 62 7.3 Suma de autovalores (EOEM) 62 7.4 Variación de σy suma de autovalores con e0(EOEM) 64 7.5 Variación de σy suma de autovalores con a0(EOEM) 65 7.6 Variación de σy suma de autovalores con i0(EOEM) 66 8.1 Comparación de suma de autovalores y σcon los distintos métodos estudiados 70 8.2 Comparación de suma de autovalores y σcon los distintos métodos estudiados (e0=0.3) 72 8.3 Comparación de suma de autovalores y σcon los distintos métodos estudiados (e0=0.6) 73 9.1 Efecto del J2en la órbita de referencia 77 9.2 Evolución de Ωcon el tiempo 78 9.3 Evolcución de ωcon el tiempo 78 9.4 Tabla para la obtención de los parámetros del modelo exponencial de densidad [12] 79 9.5 Efecto de la resistencia aerodinámica en la órbita de referencia 80 9.6 Evolución de acon el tiempo 80 9.7 Evolución de econ el tiempo 81 9.8 Efecto de la altura en σy suma autovalores 83 9.9 Efecto del J2 en σy suma autovalores en MC y UT 84 9.10 Efecto de la resistencia aerodinámica en σy suma autovalores en MC y UT 85 Índice de Tablas 4.1 Elementos orbitales iniciales 33 4.2 Variación en excentricidad (MC) 36 4.3 Variación en semieje mayor (MC) 36 4.4 Variación en inclinación (MC) 37 4.5 Resultados test de normalidad 37 5.1 Variación en excentricidad (LIN) 45 5.2 Variación en semieje mayor (LIN) 45 5.3 Variación en inclinación (LIN) 46 6.1 Variación en excentricidad (UT) 52 6.2 Variación en semieje mayor (UT) 53 6.3 Variación en inclinación (UT) 53 7.1 Variación en excentricidad (EOEM) 63 7.2 Variación en semieje mayor (EOEM) 63 7.3 Variación en inclinación (EOEM) 63 8.1 Elementos orbitales iniciales 69 8.2 Evolución de la incertidumbre en 20000s con tres métodos 69 8.3 Evolución de la incertidumbre en 7 días con tres métodos 69 8.4 Evolución de la incertidumbre en 30 días con tres métodos 69 93 Bibliografía [1] Paloma Belén Mangas Carbajo, Desarrollo de un simulador de vuelo orbital. estudio de perturbaciones, Trabajo fin de carrera, Departamento de Bioingeniería e Ingeniería Aeroespacial. Universidad Carlos III de Madrid, 2013. [2] Donald Ebeigbe, Tyrus Berry, Michael M. Norton, Andrew J. Whalen, Dan Simon, Timothy Sauer, and Steven J. Schiff, A generalized unscented transformation for probability distributions, ArXiv (2021). [3] United Nations Office for Outer Space Affairs, Number of satellites in low earth orbits, extraído el 25-06-2023 de http:// www.unoosa.org/ oosa/ en/ ourwork/ topics/ space-debris/ index.html. [4] Juan Luis Gonzalo Gómez, Propagation and Optimal Control of Space Trajectories Using Perturbation Methods, Ph.D. thesis, Departamento de Física Aplicada a las Ingenierías Aeronáutica y Naval. Universidad Politécnica de Madrid, 2017. [5] Nathan D. Aragon Joshua T. Horwood and Aubrey B. Poore, Gaussian sum filters for space surveillance: Theory and simulations, Journal of guidance, control, and dynamics 34 (2011), no. 6. [6] Jan Krejčí, Uncertainty Propagation for Tracking of Moving Objects, Ph.D. thesis, Universidad de Bohemia del Oeste, 2018. [7] Pierluigi Di Lizia, Uncertainty estimation and propagation with space applications [sesión de conferencia], Politécnico de Milán, (2021). [8] Matt Miller, Geofísica de la tierra sólida, Apuntes de clase (Universidad de Concepción), extraídos de https:// www.mttmllr.com/ geoTS_files/ GTS_S07.pdf . [9] Ester Gutierrez Moya, Estadística e investigación operativa, Apuntes de clase (ETSI Sevilla). [10] Manuel Ruiz Pérez, Uncertainty propagation in low earth orbits, Trabajo fin de grado, Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos. ETSI Sevilla, 2019. [11] Rafael Vázquez Valenzuela, Mecánica orbital y vehículos espaciales, Apuntes de clase (ETSI Sevilla). [12] David A. Vallado, Fundamentals of astrodynamics and applications, 4 ed., Springer, 2001. [13] Ming Zhou and Yongzhao Shao, A powerful test for multivariate normality, Journal of Applied Statistics 41 (2014), no. 2, 351–363. 95 96 Bibliografía [14] Manuel Ángel Zapata Habas, Análisis de trayectorias de crucero de aviones comerciales sujetas a incertidumbre en los datos, Trabajo fin de grado, Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos. ETSI Sevilla, 2015.