Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite

i Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite m Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Autor: Fernando Soler Lanagrán Tutor: Rafael Vázquez Valenzuela Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2015 ii iii Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Autor: Fernando Soler Lanagrán Tutor: Rafael Vázquez Valenzuela Profesor titular Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2015 iv v Proyecto Fin de Carrera: Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Autor: Fernando Soler Lanagrán Tutor: Rafael Vázquez Valenzuela El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2015 El Secretario del Tribunal vi vii A mi familia A mis maestros viii Índice Índice viii Notación x 1 Introducción 11 1.1 Contextualización 11 1.2 Justificación del proyecto 12 1.3 Objetivos 12 1.4 Estructura del trabajo 13 2 Sistemas de referencia 14 2.1 Conceptos previos y observaciones 14 2.2 Sistema Geocéntrico Ecuatorial 15 2.3 Sistema Geográfico 15 2.4 Sistema Topocéntrico 16 2.5 Transformaciones 17 2.5.1 Transformación de a 17 2.5.2 Transformación de a 17 3 Elementos orbitales 19 4 Obtención de las efemérides de los satélites 22 4.1 Propagación de órbitas 22 4.1.1 Conversiones entre distintas representaciones matemáticas de una órbita. 22 4.1.2 Creación de un algoritmo de propagación 27 4.2 Obtención de las efemérides reales de la constelación GPS 36 4.2.1 GPS Week (“Semana GPS”) 36 4.2.2 Acceder a datos IGS con el número GPSweek 38 4.2.3 Descomprimir archivo .Z 38 4.2.4 Lectura e interpretación de los archivos .sp3 39 4.2.5 Procesado de los datos procedentes de los archivos .sp3 40 4.3 Interpolación de los datos de las órbitas 42 4.3.1 Análisis del interpolador 43 5 Calidad del servicio de posicionamiento 48 5.1 Visibilidad de los satélites. Ángulo de máscara. 48 5.1.1 Cálculo del número de satélites visibles 49 5.2 Cálculo de la posición con cuatro satélites. 50 5.2.1 Linealización de las ecuaciones 51 5.2.2 Inicio del sistema de posicionamiento. Tipos de arranque. 52 5.3 Algoritmo de mínimos cuadrados para GPS. 52 5.4 Factores PDOP y cálculo de errores. 54 5.5 Parámetros de medida de la calidad del servicio 57 6 Cálculo de rutas 59 6.1 Ruta ortodrómica. 59 6.2 Ruta loxodrómica. 60 ix 7 Cálculo de la calidad del servicio 62 7.1 Implementación del cálculo de los factores DOP y de la disponibilidad 62 7.1.1 Estación en tierra 62 7.1.2 Vuelo 63 7.1.3 Región 64 8 Interfaz de la herramienta 66 8.1 Representaciones 3D y paquete M_map 66 8.2 Inicio de la herramienta 68 8.3 Análisis de estación en tierra 69 8.4 Análisis de continuidad en vuelo 84 8.5 Análisis de una región 88 9 Práctica de Fundamentos de Navegación Aérea 94 9.1 Observaciones previas 94 9.2 Boletín de la práctica 95 9.3 Cuestionario a los alumnos 97 9.3.1 Resultados de la valoración de los alumnos 98 10 Análisis con la herramienta 100 10.1 Comparación del análisis de estación 100 10.2 Comparación de análisis de disponibilidad manipulando satélites en la constelación 103 10.3 Disponibilidad global 105 11 Constelación Galileo 110 12 Conclusiones 119 12.1 Trabajo futuro 120 Apéndice A. Boletín de prácticas 121 Referencias 127 Sistemas de referencia 16 16 2.4 Sistema Topocéntrico El sistema Topocéntrico está también ligado a la Tierra, y tiene su origen en el punto donde se encuentra el observador ( ). El plano EXY es tangente a la superficie de la Tierra, la dirección EX apunta al Norte, la EY al Este, y la EZ hacia el centro de la Tierra [5]. Las coordenadas angulares asociadas a este sistema son: elevación , azimut y la distancia al objeto . Figura 2.3. Sistema de referencia Geográfico. X Y Z Ecuador D E N Figura 2.4. Sistema de referencia Topocéntrico. 17 17 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 2.5 Transformaciones Para pasar los vectores de un sistema de posición a otro emplearemos matrices de transformación. Para simplificar la notación, denotaremos al sistema de referencia Geográfico con un subíndice , el sistema Geocéntrico Ecuatorial con un subídice , y el Topocéntrico con una . Equation Chapter (Next) Section 1 2.5.1 Transformación de a La diferencia entre ambos sistemas es un giro en torno al eje OZ de ángulo el valor de GST en un determinado instante. cos(GST) sin(GST) 0 sin( ) cos( ) 0 0 0 1 G I C GST GST       (2.1) La transformación inversa consiste en multiplicar por la matriz inversa a la anterior, que en este caso coincide con su transpuesta:     IT I G G G I I C C C (2.2) 2.5.2 Transformación de a Para realizar dicha transformación primero realizamos un giro en torno al eje de valor el ángulo para pasar a un sistema intermedio . De esta forma alinearemos los ejes y OE. A continuación, para hacer coincidir el eje X con ON y el eje Z con OD realizamos un giro negativo en torno al eje de valor . Figura 2.5. Transformación G-T. Sistemas de referencia 18 18 De esta forma, la matriz de transformación queda: sin( ) 0 cos( ) cos( ) sin( ) 0 sin( )cos( ) sin( )sin( ) cos( ) 0 1 0 sin( ) cos( ) 0 sin( ) cos( ) 0 cos( ) 0 sin( ) 0 0 1 cos( )cos( ) cos( )sin( ) sin( ) T G C                                                          (2.3) Por último, si quisiéramos una transformación entre I y T o viceversa, podemos hacer la transformación multiplicando de forma sucesiva las matrices anteriores, por ejemplo: T G T I I G C C C (2.4) Equation Chapter (Next) Section 1 19 19 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 3 ELEMENTOS ORBITALES na órbita genérica puede definirse mediante 6 parámetros que determinan de forma total y unívoca las características, posición y orientaciones de la órbita, así como la posición que ocupa el cuerpo dentro de esta en un instante determinado. La revisión de conceptos relacionados hecha en este capítulo está basada en el material correspondiente a la referencia número [5]. Una órbita es una trayectoria que sigue un cuerpo en torno a otro como consecuencia de la atracción gravitatoria. Dicha trayectoria es siempre una cónica. En este proyecto nos centraremos en casos de órbitas cerradas elípticas. El cuerpo celeste en torno al cual orbita el satélite ocupa el foco de la elipse. El punto de la órbita más cercano al foco se denomina periapsis, y el punto más lejano es el apoapsis. Particularizando para el caso de órbitas en torno a la Tierra, se les llama respectivamente perigeo y apogeo. Figura 3.1. Parámetros geométricos de una elipse. La línea que une ambos puntos se denomina línea de ápsides. Por otro lado necesitamos un plano de referencia para ubicar el de la órbita. En el caso de órbitas en torno a la Tierra dicho plano de referencia es el plano del Ecuador. Si la órbita está inclinada respecto al plano del Ecuador, a lo largo de su movimiento el satélite ocupará tanto posiciones por debajo como por encima de dicho plano. El punto en que el satélite cruza el plano ecuatorial en dirección Sur-Norte se denomina nodo ascendente , mientras que si lo hace en dirección Norte-Sur se denomina nodo descendente . Se define el vector en la dirección que une el foco con el nodo ascendente, y el vector excentricidad tiene de módulo la excentricidad de la órbita y de dirección la que une foco con periapsis. U apoapsis periapsis Elementos orbitales 20 20 Una vez vistos todos los elementos explicados anteriormente, podemos definir los elementos orbitales. Dichos parámetros son (considerando exclusivamente el caso de órbita elíptica que nos atañe): · El parámetro (semieje mayor) determina su tamaño. · La excentricidad , determina la forma. · El parámetro , llamado el argumento del perigeo (o de periápsis) es el ángulo medido en el plano orbital y en la dirección del movimiento entre y . · El parámetro es la ascensión recta del nodo ascendente, también llamado RAAN por sus siglas en inglés. Es el ángulo medido en el sentido contrario a las agujas del reloj entre y la dirección . · La inclinación mide el ángulo entre el plano de la órbita y el plano de referencia. · Un parámetro que permita localizar el cuerpo dentro de la propia órbita en un instante dado. Se suelen utilizar la anomalía verdadera o un desde periapsis, pero también otros como (anomalía media) o (anomalía excéntrica). Como se ha explicado, los 6 parámetros anteriores determinan de forma unívoca una órbita cualquiera, pero existen casos en que dichos elementos orbitales no están bien definidos. Estos casos son importantes porque se dan en la práctica y, de hecho, algunos deberemos considerarlos en este proyecto (tal es el caso de las órbitas circulares). - Órbitas elípticas ecuatoriales: no están definidos ni , ya que al estar contenidas en el plano ecuatorial no existe línea de nodos. En dicho caso se sustituyen por el ángulo (longitud del perigeo) entre y , medido en el sentido contrario de las agujas del reloj, de manera que . - Órbitas circulares no ecuatoriales: y no están bien definidos (no existe línea de ápsides). Se sustituyen por un ángulo llamado argumento de la latitud, que mide la posición del cuerpo desde el nodo ascendente en el sentido del movimiento. Plano del Ecuador Plano de la órbita Figura 3.2. Órbita genérica. 21 21 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite - Órbitas circulares ecuatoriales: no existe ni línea de nodos ni línea de ápsides. Por lo tanto ni , ni ni están bien definidos. Se sustituyen por , la longitud verdadera, que es la que forma el cuerpo con , medido en el sentido contrario de las agujas del reloj. Equation Chapter (Next) Section 1 Obtención de las efemérides de los satélites 22 22 4 OBTENCIÓN DE LAS EFEMÉRIDES DE LOS SATÉLITES a creación de la herramienta software para el análisis de una constelación por satélite necesita una gran cantidad de información y de cálculos intermedios a procesar entre las órdenes del usuario y los resultados que se le proporcionan. Para lograr lo anterior hemos de alcanzar una serie de objetivos básicos: 1) La obtención de la posición de los satélites, ya sea a partir de una fuente externa de datos reales, o mediante la simulación de los mismos (propagación de órbitas). 2) La creación de algoritmos y funciones que permitan calcular todos los parámetros que se espera que la aplicación sea capaz de obtener. 3) La implementación de lo anterior en una aplicación, más concretamente a través de GUIDE en MATLAB, que ofrezca al usuario un entorno de trabajo lo más sencillo, claro y eficiente posible. En este y los siguientes capítulos se expondrán los métodos desarrollados para alcanzar las metas anteriormente planteadas. En este capítulo vamos a describir los dos métodos empleados para la obtención de las efemérides de los satélites: la propagación de órbitas a partir de sus elementos orbitales y la obtención de efemérides reales de la constelación GPS a través del servicio web del IGS (International GNSS Service). 4.1 Propagación de órbitas La propagación de órbitas consiste en predecir o estimar la posición futura de los satélites a partir de efemérides conocidas en una determinada época. Como se explicó anteriormente, si conocemos todos los elementos orbitales de un satélite (6 parámetros) conocemos su órbita de forma precisa y podemos realizar la propagación empleando las leyes de la mecánica orbital. 4.1.1 Conversiones entre distintas representaciones matemáticas de una órbita. 4.1.1.1 Función para transformar vectores posición y velocidad en elementos orbitales. Para el cálculo de los elementos orbitales a partir de los vectores velocidad y posición podemos emplear las fórmulas que nos proporciona la mecánica orbital [5] y la definición de los parámetros geométricos de las órbitas que se explicaron anteriormente: En primer lugar determinamos el parámetro vectorial , que es el momento angular específico: h r v (4.1) Y con él obtenemos el vector excentricidad:   /v h r r e     (4.2) Donde el parámetro es el Parámetro Gravitacional de las Masas Combinadas, que en nuestro caso puede aproximarse por con la constante de gravitación universal y la masa de la Tierra. El valor de para la Tierra es de aproximadamente . L 23 23 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Calculando el módulo de obtenemos el valor de la excentricidad: ee (4.3) A continuación calculamos la energía específica : 2 2 v r    (4.4) Y con ello el semieje mayor : 2 a    (4.5) Para obtener empleamos: 1 cos cos e r e r er er         (4.6) Para resolver la ambigüedad: si , el satélite está viajando de apoapsis a periapsis. Entonces . Sabiendo que es perpendicular al plano orbital, se obtiene que:   cos(i) , 0, Iz h hk i hh      (4.7) Por otro lado: 2 2 2 2 , ,0 I y x y x x y I kh n h h h h h h kh         (4.8) Para obtener empleamos: cos( ) nx  (4.9) Resolviendo la ambigüedad sabiendo que: si el vector está a la izquierda de , luego . Por último, es el ángulo entre y , luego: cos( ) e n e   (4.10) Corregimos el cuadrante con: si el vector excentricidad está bajo el plano de referencia y . Ya hemos determinado, con las fórmulas (4.3), (4.5), (4.6), (4.7), (4.9) y (4.10) los elementos , , , , , y . Nótese que para poder aplicar las fórmulas anteriores es preciso que los vectores estén expresados en el sistema de referencia Geocéntrico Ecuatorial. Las ecuaciones anteriores, pueden ser implementadas en MATLAB para crear una función, llamada “rvael” que proporcione: [a, e, inc, w, Omg, theta]=rvael(r, v) Sin embargo el algoritmo anterior no está completo si no se adapta para el cálculo de los casos límite mencionados en el capítulo 3. Para ello hacemos lo siguiente: Obtención de las efemérides de los satélites 24 24 Si la órbita es circular no ecuatorial, cambiamos por y dejamos indefinido como NaN. Si la órbita es ecuatorial con excentricidad, cambiamos por y dejamos indefinido como NaN. Si la órbita es ecuatorial circular, cambiamos por y dejamos y indefinidos como NaN. En caso de que la órbita sea parabólica, se sustituye por . Aunque no será necesario en este proyecto, dado que los satélites de navegación siguen órbitas cerradas. 4.1.1.2 Función para transformar elementos orbitales en vectores posición y velocidad. Para este caso es conveniente definir un sistema de referencia auxiliar, denominado perifocal y que identificaremos con el superíndice , cuyas características son las siguientes: Está centrado en el foco, con apuntando hacia periapsis, y el eje hacia “arriba”, paralelo al vector . Este sistema de referencia se utiliza porque en él es más sencillo calcular los vectores y : Se puede demostrar [5] que, definiendo un sistema de referencia polar sobre el perifocal, con parámetro y empleando como ángulo , podemos expresar: cos sin 1 cos 0 Fp re         (4.11) sin cos 0 F ve p          (4.12) 2 h p   (4.13) Plano del Ecuador Plano de la órbita Figura 4.1. Sistema de referencia perifocal. 25 25 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Una vez obtenemos los vectores y sólo nos falta hacer una transformación para expresarlos en el sistema de referencia Geocéntrico Ecuatorial, para lo cual necesitamos la matriz . Obtendremos primero la matriz : Realizamos un giro en torno al eje con un ángulo y pasamos al sistema intermedio : Giramos en torno al eje un ángulo para pasar al sistema : Plano del Ecuador Plano de la órbita Obtención de las efemérides de los satélites 32 32 Esto significa que nuestro propagador nos proporcionaría una información pobre en torno al perigeo. Figura 4.5. Propagación de una órbita con alta excentricidad con intervalos de tiempo equiespaciados. La Figura 4.5 nos permite ilustrar este efecto. Dicha figura se ha obtenido con una órbita de elementos orbitales . 2) Obtención de la anomalía excéntrica por métodos numéricos. Otro aspecto a considerar es que la ecuación (4.23) es fácil de resolver para obtener conociendo la anomalía excéntrica , ya que es conocido. Por lo tanto es inmediato obtener la anomalía media a partir de . Pero si queremos hallar el valor de a partir de es necesario recurrir a métodos numéricos para obtener la solución. En general los métodos numéricos empleados para resolver este tipo de ecuaciones requieren un valor inicial para iterar. El valor de es adecuado para ello, y la solución converge. El problema es que la resolución numérica incrementa notablemente el tiempo de cálculo, lo cual debemos tener en cuenta a la hora de crear un algoritmo de cálculo ya que con mayor número de puntos calculados el efecto de esta ralentización se hace más notable. Repitiendo la simulación de la Figura 4.5 en un script “Untitled” empleando el temporizador, obtenemos los datos de la siguiente tabla: 33 33 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 4.6. Temporizador aplicado a la función de propagación con intervalos de tiempo equiespaciados. De la tabla anterior nos interesa el tiempo de “fzero”, que resuelve la ecuación (4.23), respecto al tiempo de ejecución de “propagador_n_dt” que es la función principal. Se observa que el tiempo consumido por la resolución numérica de es el del total, confirmando así el importante peso de este paso del algoritmo respecto al tiempo de cálculo total. 4.1.2.4 Optimización de la función “propagador” Vamos a estudiar de qué forma se puede optimizar la respuesta de la función ante las dos observaciones realizadas en el apartado anterior. 1) La forma más evidente de solventar este problema es realizar la discretización de los puntos en que se calculan los elementos con otra variable que dé lugar a una discretización más uniforme en el espacio. Dado que la relación entre el tiempo y la anomalía media es lineal, emplear nos proporcionaría el mismo resultado que teníamos anteriormente. Las otras dos opciones son la anomalía verdadera y la anomalía excéntrica. Sin embargo el uso de cualquiera de estas dos variables supone igualmente un problema: ahora tenemos intervalos de tiempo irregulares, y las épocas pasan a ser una variable que calcular. 2) Dado que calcular a partir de (o ) requiere el empleo de métodos numéricos computacionalmente y temporalmente más costosos, realizar la discretización en base a o resuelve este problema. Aunque de nuevo el tiempo sería una variable que calcular (pero ya con ecuaciones que sí tienen una resolución directa sin necesidad de métodos numéricos). De esto se deduce que, en un principio, emplear discretización en o reduce el tiempo de cálculo (resolvemos el punto 2) y, además, permite calcular puntos de la órbita en el espacio con otra distribución. Obtención de las efemérides de los satélites 34 34 Probemos primero empleando una distribución equiespaciada de la anomalía verdadera . Repetimos la simulación de la Figura 4.5. y obtenemos lo siguiente: Figura 4.7. Propagación de una órbita con alta excentricidad con intervalos de equiespaciados. Curiosamente sucede el fenómeno contrario. Sin embargo era de esperar, ya que si empleamos la anomalía verdadera para el cálculo con intervalos equiespaciados, el arco que se recorre sobre la órbita es mayor cuanto mayor sea el radio. Es decir, en el perigeo, donde el radio es menor, hay una mayor concentración de puntos que en el apogeo, donde el radio es máximo. Respecto al tiempo de cálculo: Figura 4.8. Temporizador aplicado a la función de propagación con intervalos de equiespaciados. Efectivamente el tiempo de cálculo se ha reducido drásticamente: frente a . Y evidentemente esta diferencia es mayor cuantos más puntos se calculen. Esta solución resuelve el segundo de los problemas, pero no el primero. 35 35 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite A continuación probamos con una distribución equiespaciada de . Repetimos el análisis de los casos anteriores, obteniendo lo siguiente: Figura 4.9. Propagación de una órbita con alta excentricidad con intervalos de equiespaciados. Como vemos, la anomalía excéntrica proporciona puntos equiespaciados. Además también soluciona el problema del tiempo de cómputo: Figura 4.10. Temporizador aplicado a la función de propagación con intervalos de equiespaciados. Por lo tanto se llega fácilmente a la conclusión de que si pretendemos optimizar la concentración de puntos a lo largo de la órbita y reducir el tiempo de cálculo, la solución es utilizar una discretización en los valores de la anomalía excéntrica . Sin embargo este proceso ha sido aplicado para una única órbita. Ahora toca plantearse qué sucede si queremos que el algoritmo calcule varias órbitas a la vez. Obtención de las efemérides de los satélites 36 36 La ventaja de utilizar el interpolador original con un vector de tiempo reside en que proporciona valores para todas las órbitas en unas épocas determinadas y conocidas a priori. Es decir, para cada instante de tiempo del vector, tenemos las posiciones correspondientes de todos los satélites. Si propagamos con valores de equiespaciados obtendremos puntos de las órbitas en instantes de tiempo arbitrarios unos respecto a otros. Además, mientras que con el vector de tiempo las órbitas se propagan hasta donde lo harían en el tiempo especificado por el vector (es decir, si una órbita es más rápida que otra, es posible que esta haya dado una vuelta entera en su órbita mientras la otra aún no ha completado una vuelta a la suya), si el vector de valores de llega hasta, por ejemplo , la propagación muestra todas las órbitas con una vuelta completa. Por lo tanto, este tipo de interpolación tiene también sus desventajas. Dado que para este proyecto interesa propagar varios satélites a la vez y tener dichos datos para un conjunto de épocas conocidas, conviene conservar el propagador original. Aún así, si la órbita es circular, podemos simplificar los cálculos ya que la ecuación (4.23) queda: ME 4.2 Obtención de las efemérides reales de la constelación GPS El IGS (International GNSS Service) es un servicio que proporciona, entre otras cosas, información de distintos sistemas de navegación por satélites a través de internet. Entre sus bases de datos es posible encontrar por ejemplo las efemérides reales de la constelación GPS (o GLONASS) en cualquier época. Existen tres tipos de datos: “UltraRapid”, proporcionados cuatro veces al día y que ofrecen información basada en observación y predicción de las órbitas; “Rapid”, disponibles con 17 horas de retraso respecto al la época correspondiente, y “Final”, disponibles con 12 días de retraso [6]. Los últimos son los datos más precisos. El interés de esta base de datos para el presente proyecto es evidente: obteniendo de ella los datos de las órbitas y procesando dicha información, podemos realizar los análisis con la constelación real sin la necesidad de realizar una propagación. Nos centraremos en obtener los datos tipo “Final” para una fecha deseada. 4.2.1 GPS Week (“Semana GPS”) Entre las particularidades del sistema GPS está el hecho de que cuenta con un calendario propio, que organiza el tiempo en días y semanas con lo que se conoce como GPS Week. A cada semana le corresponde un número, siendo por ejemplo la primera semana de 1994 la número 729. Esta semana es, de hecho, la primera semana de la que existen datos de la constelación. Dado que no fue declarada con capacidad operacional inicial hasta diciembre de 1993, no es de extrañar que la base de datos comience poco después, el 1 de enero de 1994, que dio comienzo a su primer año con capacidad operativa reconocida (más adelante en abril de 1995 fue declarada con capacidad operacional total)[2]. La segunda semana de 1994 es por tanto la 730, y así sucesivamente. Hay además otra peculiaridad, y es que el conteo de días de una semana no se realiza como se acostumbra en la mayoría de países de habla hispana, que es como establece el estándar internacional ISO 8601 [7] de representación de fechas y horas: el lunes es considerado primer día de la semana y esta transcurre hasta llegar al domingo. Una semana GPS comienza en un día domingo y llega al próximo sábado, que termina dicha semana GPS para dar comienzo a la siguiente en el domingo posterior. De acuerdo con lo anterior, los días se numeran del 0 al 6, siendo 0 el domingo, 1 el lunes, y así hasta el 6 sábado. Este número se suele poner a continuación del número de la semana a la que corresponde. Por ejemplo, el sábado 1 de enero de 1994 tendría el número: 7296. El domingo 2 de enero de 1994 es el 7300. 37 37 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 4.2.1.1 Convertir fecha a GPS Week Los datos de la base de datos están organizados por días siguiendo este formato, por lo que necesitamos ser capaces de hacer la conversión entre el formato usual de fecha “día/mes/año” y GPS Week. La forma más sencilla es tomar de referencia un día del cual conozcamos su número GPS Week, contar los días a partir de él, y con ello, obtener el número GPS Week de la fecha que queremos. Para realizar la cuenta de días transcurridos lo idóneo es recurrir a la fecha juliana o días julianos JD, que consisten en una cuenta de días y fracción transcurridos desde el mediodía del 1 de enero del año 4713 a.C. Esta escala temporal fue creada en el año 1582 por Joseph Scaliger de Leyden para poder fechar fenómenos astronómicos o históricos lejanos en el tiempo evitando ambigüedades como años bisiestos o cambios de calendarios [9]. Para calcular el día juliano (JD) correspondiente a una fecha, podemos emplear la fórmula [5]: 9 77 275 12 367 1721013.5 49 M A INT M JD A INT INT D                     (4.30) Donde la función denota la parte entera de un número real, representa el año, el mes y el día. La fórmula anterior proporciona el dato a las 00:00 del día correspondiente. Para obtenerlo a otra hora basta con sumar la hora UT dividida por 24. A nosotros nos basta con la ecuación (4.30). Tomaremos de referencia el día 2 de enero de 1994 (02/01/1994), al que corresponde la semana GPS 730, y el día 0. Denotamos su día juliano como . El día juliano de la fecha que deseamos convertir es . El número de semanas transcurridas se puede hallar como: 0 7 e w JD JD N INT     (4.31) Y el día como: 00 7 77 ee w JD JD JD JD D INT                    (4.32) Por tanto la semana GPS correspondiente es: 730 w GPSweek N (4.33) Por lo que el número GPSweek podemos formarlo como: 10 w NúmeroGPSweek GPSweek D   (4.34) Obtención de las efemérides de los satélites 38 38 4.2.2 Acceder a datos IGS con el número GPSweek Si conocemos el número GPSweek de la fecha que queremos simular, podemos generar un enlace web para descargar el archivo. El formato general de dicho enlace es el siguiente: ftp://igscb.jpl.nasa.gov/igscb/product/*GPSWeek*/igs*númeroGPSweek*.sp3.Z Donde los elementos entre asteriscos deben ser sustituidos por el correspondiente número. Si la semana GPS correspondiente es menor de 1000, debe ir precedida por un 0, de modo que siempre sea un número de 4 dígitos. Por ejemplo, para el 2 de enero de 1994 el enlace quedaría: ftp://igscb.jpl.nasa.gov/igscb/product/0730/igs07300.sp3.Z Con los algoritmos desarrollados hasta ahora ya nos es posible hallar la semana GPS de una fecha concreta y descargar el archivo correspondiente de la base de datos IGS. 4.2.3 Descomprimir archivo .Z Los archivos de la base de datos IGS están comprimidos en un formato de compresión “.Z”. La herramienta será implementada en MATLAB, y dicho formato de compresión no es reconocido por el programa, luego aunque podamos crear líneas de código para descargar los archivos, aun necesitamos ser capaces de realizar la descompresión de forma automática, para lo que necesitamos poder hacerlo mediante ejecución de comandos externos a MATLAB a través de instrucciones al propio MATLAB en el workspace. Lo primero que necesitamos es encontrar un programa que funcione en línea de comandos. En este proyecto se empleará “Unarchiver” [9]. Una vez tenemos el programa, necesitamos que se encuentre en el mismo directorio que el archivo a descomprimir. Si esto se cumple podemos dar a MATLAB la orden de ejecutar el comando externo y habremos descomprimido con éxito el archivo. Los archivos de la base de datos del IGS se encuentran en formato “.sp3”, que puede ser leído con visores de archivos de texto. 39 39 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 4.2.4 Lectura e interpretación de los archivos .sp3 El archivo .sp3 original tiene la siguiente apariencia (tomando las primeras 30 filas): Figura 4.11. Datos contenidos en un archivo .sp3. Las primeras dos filas muestran información relativa a la fecha, (en la primera fila tenemos la fecha, en el ejemplo es el 17 de junio de 1997, en la segunda fila aparecen la semana GPS que es 910 y el día juliano correspondiente 50616). En la tercera y cuarta fila tenemos una lista de los satélites disponibles. En la fecha del ejemplo hay disponibles 25 satélites y se detalla el número correspondiente a cada uno de ellos. En la fila 23 podemos observar de nuevo la fecha en formato año/mes/día seguido de dos ceros, que representan el primero la hora, y el segundo los minutos. Si avanzamos más en el documento, estas épocas van avanzando en saltos de cuartos de hora: 0 15, 0 30, 0 45, 1 0, 1 15, etc, hasta llegar a 23 45. Figura 4.12. Cabecera de los datos correspondientes a las 01:15 am del 17/06/1997. Obtención de las efemérides de los satélites 40 40 La información correspondiente a las 00:00 del día siguiente estaría en su correspondiente archivo .sp3. Esto significa que los datos que obtenemos del IGS están espaciados temporalmente 15 minutos unos de otros. Más adelante abordaremos cómo conseguir información con una discretización en el tiempo más fina. A partir de la línea 23 (y las análogas a lo largo del documento, tal y como la 153 en la Figura 4.12) comienza la lista de datos de los satélites para cada hora. Los números de la segunda columna identifican a qué satélite corresponden los datos. La tercera columna es la coordenada en el sistema de referencia Geográfico, según nuestra notación: . La cuarta y la quinta columna son respectivamente y . Estos datos están dados en km. La última columna hace referencia a la desviación del reloj del satélite en microsegundos [11]. Para poder trabajar con los datos necesitamos convertirlos a formato matricial a través de un script que tome como entrada el archivo sp3 [12]. Los datos que obtenemos con ello aparecen en una matriz de siete columnas y una cantidad variable de filas (depende del número de satélites disponibles en la fecha): Figura 4.13. Datos “en bruto” en formato matricial obtenidos del archivo sp3. En la primera columna tenemos la GPSweek, la segunda son los segundos transcurridos desde el inicio de dicha semana GPS hasta el momento al que corresponden los datos, la columna 3 contiene el número del satélite, las columnas de 4 a 6 las componentes del vector posición y la columna 7 la desviación del reloj. 4.2.5 Procesado de los datos procedentes de los archivos .sp3 Llegados a este punto ya hemos obtenido satisfactoriamente la información de la constelación para la fecha deseada. Sin embargo lo ideal es procesar estos datos para que coincidan con el formato que utilizamos en la propagación, es decir: [xI, yI, zI, xG, yG, zG, GST, lambda, phi] 41 41 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Recordemos que el formato que hemos decidido utilizar en la herramienta es de la forma: 1 0 1 1 1 1 1 2 0 2 1 2 2 1 0 1 1 0 0 0 1 (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) (t ) jM jM i i i j i M N N N N M xG xG xG xG xG xG xG xG xG xG xG xG xG xG xG xG xG               Para crear esta matriz debemos recorrer la columna 4 de la matriz de datos del IGS atendiendo a la numeración de satélites de la columna 1, de modo que podamos reordenar los datos en la nueva matriz que queremos generar. Por ejemplo, en la Figura 4.11 la cabecera de la línea 23 nos dice que esos datos corresponden a las 00:00, que tomamos como . De modo que esos datos los asignamos al elemento , donde es o según corresponda. Para obtener y seguimos el mismo procedimiento que empleamos en la propagación con las ecuaciones (4.28) y (4.29). Para obtener los vectores posición en el sistema de referencia Geocéntrico Ecuatorial necesitamos conocer GST, que debemos calcularlo. En primer lugar necesitamos conocer . Existe una fórmula que permite calcular con bastante precisión el para un día dado en días julianos a las 00:00, que es precisamente el que necesitamos: Primero calculamos el “Siglo Juliano” : 0 2451545 36525 JD T  (4.35) A partir de ello se obtiene , en grados, como [5]: 2 8 3 0 0 0 0 100.4606184 36000.77004 0.000387933 2.583 10GST T T T       (4.36) A continuación obtenemos el con (4.26) y calculamos el vector posición en el sistema de referencia Geocéntrico Ecuatorial de la siguiente manera:     1 I IG IG Ij I x r y C GST r z        (4.37) Con esto hemos obtenido los datos que nos faltaban. Sin embargo en esta ocasión proporcionaremos una salida más, que es la lista de cuáles de los satélites de la constelación están disponibles, ya que nos será de utilidad más adelante. [lista_sat, rgx, rgy, rgz, xI, yI, zI, lambda, phi] = Datos_GPS_IGS(Dia, Mes, Año, GPS_week) Calidad del servicio de posicionamiento 48 48 5 CALIDAD DEL SERVICIO DE POSICIONAMIENTO n este capítulo se van a explicar algunos conceptos relacionados con la calidad del servicio de posicionamiento que ofrece una constelación como la GPS, así como los procedimientos para calcular algunos parámetros importantes relacionados. Varias partes de este capítulo se basan en la referencia [3]. 5.1 Visibilidad de los satélites. Ángulo de máscara. El empleo de la constelación GPS como sistema de posicionamiento se basa en la posibilidad de determinar la posición de un observador conocidas las posiciones de los satélites y empleando señales de radio. Por lo tanto, para que el sistema funcione correctamente necesitamos un número suficiente de satélites y que estos sean visibles por el observador. En general, si un satélite está bajo el horizonte implica que el observador no puede recibir su correspondiente señal. Pero considerar un satélite visible desde el momento en que cruza el horizonte es poco realista (salvo casos como un vuelo o un barco en alta mar) ya que en general los elementos del entorno (edificios, accidentes geográficos…) impiden ver los satélites hasta que estos se eleven un poco más en el cielo. Al ángulo mínimo que deben elevarse los satélites desde el horizonte hacia el cenit se le llama ángulo de máscara. Generalmente este ángulo es función del azimut, ya que los mencionados objetos que obstaculizan la recepción de la señal pueden tener cualquier tipo de disposición alrededor del observador. Sin embargo esto complica los cálculos y necesitaríamos conocer la distribución de estos elementos para cada caso concreto. En este estudio asumiremos un ángulo de máscara constante en todas las direcciones. Este parámetro nos determinará el número de satélites visibles para el observador en una época determinada. No obstante, sería posible tener en cuenta un ángulo de máscara variable dado que los algoritmos que desarrollaremos en este capítulo nos proporcionarán para cada satélite y época el correspondiente azimut y elevación. Bastaría con proporcionar a la herramienta la función que relacione el ángulo de máscara para cada azimut, de modo que en cada época la herramienta evalúe la visibilidad de los satélites según si la elevación que tienen es o no superior al ángulo de máscara correspondiente a su azimut. Dicha implementación es sencilla, pero no se ha realizado porque solicitar al usuario la función de ángulo de máscara en función del azimut es más propio de un análisis avanzado, y el objetivo de esta herramienta es ser de empleo sencillo e ilustrativo. Con un ángulo de máscara constante es suficiente en una primera versión didáctica. Sin embargo más adelante se explicará cómo sería el proceso para ello. E 49 49 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 5.1.1 Cálculo del número de satélites visibles La forma de determinar si un satélite es visible o no en una determinada época desde una posición sobre la Tierra consiste en verificar si su elevación respecto al horizonte es mayor que el ángulo de máscara. Para hacer esto primero necesitamos conocer la posición del satélite respecto a la estación en Tierra situada en . A dicho vector lo llamaremos . El vector posición del satélite lo llamaremos . Si conocemos las coordenadas y de la estación que queremos analizar, podemos aproximar el vector posición de la estación como: ( R el radio de la Tierra) cos cos cos sin sin ee G e e e e rR          (5.1) El vector lo obtuvimos en el capítulo anterior, luego sólo nos queda hacer: G G G e s r r (5.2) Y transformamos el vector al sistema de referencia Topográfico centrado en la estación con la matriz (2.3): TG T G s C s (5.3) El vector se puede expresar como: cos(h)cos(Az) cos(h)sin(Az) sin(h) N TE D s s s s s             (5.4) Figura 5.1. Vector en sistemas de referencia Geográfico y Topocéntrico. Calidad del servicio de posicionamiento 50 50 Por lo que las coordenadas angulares y podemos obtenerlas haciendo:  1 sin D s hs   (5.5)  1 tan E N s Az s   (5.6) Con y . Ya podemos imponer que la condición para que un satélite sea visible en una época es que cumpla: min hh (5.7) Donde es el ángulo de máscara. Esta información hemos de tenerla en cuenta ya que aquellos satélites no visibles en una época no aportarán nada al servicio de posicionamiento y por lo tanto no deben ser tenidos en cuenta. 5.2 Cálculo de la posición con cuatro satélites. Las medidas del receptor GPS se denominan observables, y con ellas se puede determinar la posición del observador. A partir de la señal enviada por un satélite es posible determinar el tiempo en que se enviaron. Comparando con el tiempo de recepción, obtenemos la diferencia de tiempos Sin embargo los relojes de cuarzo empleados por los receptores tienen errores, que llamaremos , de modo que . Llamando a la distancia entre el observador y el satélite, y siendo la velocidad de la luz, podemos decir que: uu r c t ct ct       (5.8) Donde es la pseudodistancia. Si llamamos al vector posición del satélite y al vector posición del receptor, se cumple que: r s u (5.9) Por lo tanto, para cada satélite visible en una época tenemos una ecucación del tipo: i u i ct s u     (5.10) Donde se considera conocido ya que los satélites transmiten sus efemérides al receptor con suficiente precisión. De esta forma, para cada satélite tenemos una esfera en la que puede encontrarse el observador, determinada por la distancia al satélite. Dado que la intersección de dos esferas determina una circunferencia, y la intersección de dos circunferencias determina en general dos puntos, necesitamos al menos tres satélites para determinar la posición. Además el valor de es desconocido, por lo que el número mínimo de satélites necesarios es de cuatro. Tenemos las siguientes ecuaciones: 11 22 33 44 u u u u ct s u ct s u ct s u ct s u                 (5.11) 51 51 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 5.2.1 Linealización de las ecuaciones Si definimos ui u i i ui xx u y s y zz                      (5.12) tenemos que 2 2 2 (x x ) (y y ) (z z ) (x ,y ,z ,t ) i i u i u i u u i u u u u ct f          (5.13) Supongamos que conocemos una estimación inicial de u y u t , denotadas por ˆ u y ˆu t , que deben estar suficientemente cerca de la posición real. Definimos: ˆˆ ˆ ˆˆ ,, u u u u u i i u u x u u u y t t t s u ct z                    (5.14) Si linealizamos las funciones en torno a la estimación inicial tendremos que:               ˆ ˆˆˆ , , , , , , ˆˆ ˆ ˆ ˆ ˆ ˆˆ , , , , , , ˆ ˆˆˆ , , , ˆˆ ˆ ˆ ˆ ˆ ˆˆ , , , , , , i i u u u u i u u u u u u u u i u u u u i u u u u i u u u u u u uu i u u u u i u u u u uu uu f x y z t f x x y y z z t t f x y z t f x y z t f x y z t x y xy f x y z t f x y z t zt zt                      (5.15) Por tanto, la linealización queda: ˆHx         (5.16) Donde: 1 1 1 11 2 2 2 22 3 3 3 33 4 4 4 44 ˆ1 ˆ1 ,, ˆ1 ˆ1 x y z u x y z u x y z u x y z u a a a x a a a y xH a a a z a a a ct                              (5.17) El valor es la estimación de la pseudodistancia y los valores precedidos de son desviaciones respecto a un valor de referencia estimado de . Los coeficientes , y son de la forma: 2 2 2 ˆ (x x ) ˆˆ ˆ ˆ (x x ) (y y ) (z z ) i i u xi ui u i u i u f ax          (5.18) Y análogamente con y . Calidad del servicio de posicionamiento 52 52 Por tanto para determinar x simplemente basta con hacer 1 xH      (5.19) y se obtienen los errores respecto a la estimación inicial. 5.2.2 Inicio del sistema de posicionamiento. Tipos de arranque. Generalmente para la estimación inicial que se emplea en la linealización se puede tomar la posición calculada inmediatamente antes ya que será próxima a la posición real en la época. Si no disponemos de un valor inicial no tenemos en torno a qué linealizar. Esto sucede si el receptor, durante el tiempo que ha estado inactivo, ha cambiado su localización significativamente respecto a la última posición que calculó, o si dicha información se pierde por ejemplo por falta de una fuente de alimentación. En dicho caso se necesita obtener de algún modo información suficiente para resolver las ecuaciones (5.11) y comenzar a proporcionar la posición del usuario. A este proceso se le llama “GPS lock”. Existen varios tipos de arranque, que serán más o menos rápidos en función de la información de la que disponemos. El más rápido (“hot start”) se realiza cuando el receptor conserva la información de la última posición calculada y los satélites visibles, el almanaque empleado (información aproximada de las órbitas y el estado de los satélites de la constelación), el tiempo UTC y basándose en esta información se intenta obtener la posición con los mismos satélites. Este método sólo funciona si la posición es aproximadamente la misma que se tenía la última vez que se desactivó el dispositivo. El proceso es algo más lento (“warm start”) cuando el dispositivo conserva la última posición calculada, los datos del almanaque y el tiempo UTC pero no sabe cuáles son los satélites visibles desde su posición (por ejemplo, porque esta ha cambiado). En dicho caso realiza un reset de la posición y trata de calcular qué satélites son visibles en base a los datos del almanaque y obtener las señales de los satélites para calcular una nueva posición. Finalmente el proceso más lento (“cold start”) se da cuando el dispositivo GPS no dispone de ninguna información y por lo tanto tiene que obtener toda la información de la que carece y procesarla, situar los satélites y comenzar a estimar la posición sin información o cálculos previos de referencia [14]. 5.3 Algoritmo de mínimos cuadrados para GPS. El proceso detallado en el apartado 5.2 es aplicable cuando disponemos de información de cuatro satélites. En la práctica se emplean los datos de más de cuatro satélites, ya que si contamos con más información podemos conseguir resultados más precisos. Entonces hay que resolver un problema sobredeterminado, para lo que se emplea normalmente el algoritmo de mínimos cuadrados. Además, dado que las medidas contienen error, se podría usar un modelo del error del tipo Hx      (5.20) Donde para el caso general de satélites se tiene que   es , H es , x es y  es un modelo de error para la pseudodistancia que cumple que sigue una distribución normal (0, ) n N . Este modelo es apto para ser resuelto con el algoritmo de mínimos cuadrados, que además permitirá estimar el error que se está cometiendo en x a partir del valor de H y de  . Una interpretación geométrica del método de mínimos cuadrados es la búsqueda de la mejor aproximación al vector   mediante vectores generados por las columnas de la matriz H . Esta mejor aproximación es la proyección ortogonal de   sobre el espacio columna de H . 53 53 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Para ilustrarlo gráficamente: Como vemos, consiste en tomar la solución disponible lo más cercana a   , y el vector es el que nos da tal solución. Para hallarlo basta con tener en cuenta que es perpendicular al espacio columna de , ya que es la proyección ortogonal sobre dicho espacio. Por lo tanto:   ˆˆ 0, T T T H H x H H H x         (5.21) Lo que nos lleva a [15]:   1 ˆTT x H H H      (5.22) Respecto al modelo de error  , es razonable suponer que la matriz de covarianza  tiene la siguiente forma: 2 1 2 2 2 n           (5.23) Donde 2 i  es la varianza del error de cada pseudodistancia. No obstante es usual despreciar la diferencia de errores entre satélites. Asumimos que el error sigue una distribución con la misma varianza en todos ellos, de modo que se toma 22 i UERE   . El término “UERE” procede de “User Equivalent Range Error”, una estimación cuyo valor típico es 27 1.5 UERE m   según si el sistema de posicionamiento es el standard (SPS) o el sistema de posicionamiento preciso (PPS) y procede de las siguientes fuentes de error, sumadas con RSS: 2 2 2 12UERE causa causa       : - Segmento espacial: error del reloj (1.1 m), cálculo de la órbita (0.8 m). - Segmento usuario: efectos atmosféricos, ruido del receptor y resolución, efectos multicamino: 7-1.4 m (SPS - PPS). Espacio columna de la matriz Figura 5.2. Ilustración geométrica simple para caso del método de mínimos cuadrados. Calidad del servicio de posicionamiento 54 54 De modo que en primera aproximación, la ecuación (5.23) queda: 2 UERE n I   (5.24) Con la matriz identidad de orden . Por lo tanto la covarianza del resultado será:   1 2 1 2 1 ˆ(H H) H (H H) (H H) T T T T UERE n UERE Cov x I H        (5.25) Definimos: 1 (H H) T G  (5.26) Y llegamos a:   2 ˆUERE Cov x G   (5.27) El significado físico de este término es el que sigue:                           2 2 2 2 ˆ u u u u u u u u u u u u u u u u u u u u u u u u u u u u Var x Cov x y Cov x z Cov x t Cov y x Var y Cov y z Cov y t Cov x Cov z x Cov z y Var z Cov z t Cov t x Cov t y Cov t z Var t                                              (5.28) Estas varianzas son el producto de dos factores: 2 UERE  , que depende de la señal, y , que depende de , la cual es a su vez función de la geometría, y nos permiten estimar el error cometido. 5.4 Factores PDOP y cálculo de errores. De la anterior matriz nos interesan los valores de la diagonal, ya que nos dicen la varianza en las diferentes direcciones y el tiempo. De esto se saca la conclusión de que los errores cometidos por la señal ( UERE  ) son multiplicados a su vez por factores que dependen de la geometría ( ). Estos factores son llamados factores DOP (dilution of precision). Los valores usados típicamente son: GDOP (Geometric Dilution of Precision)= 11 22 33 44 G G G G   . PDOP (Position Dilution of Precision)= 11 22 33 G G G . HDOP (Horizontal Dilution of Precision)= 11 22 GG . VDOP (Vertical Dilution of Precision)= 33 G TDOP (Time Dilution of Precision)= 44 G Mediante estos factores DOP podemos hallar rápidamente una estimación de la precisión de nuestro GPS: / Z UERE t UERE VDOP TDOP c     (5.29) 55 55 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Otra medida comúnmente usada es el 2DRMS: círculo que contiene del 95% al 98% de los puntos de la distribución estadística que sigue la posición estimada en el plano horizontal. Se obtiene como un círculo de radio 2DRMS: 22 UERE DRMS HDOP   (5.30) Veamos cómo obtener también intervalos de confianza para los valores en la dirección vertical. Dado que estamos asumiendo que la distribución del error sigue una normal, la probabilidad de intervalos para esta distribución puede consultarse en cualquier libro de estadística, ya que aparecen tabuladas. Podemos asumir que en torno al 68% de los valores están a una distancia UERE VDOP   del valor estimado, alrededor del 95.4% están a dos desviaciones típicas 2UERE VDOP   de la media y alrededor del 99.7% están a 3UERE VDOP   . Los factores PDOP son factores que multiplican el error, por lo que es conveniente que sean lo más pequeños posible. Estos factores son función de la geometría, cuyo efecto depende de las posiciones de los satélites en la época. Desde un punto de vista matemático, si los satélites no cubren bien el cielo la matriz H (cuyas componentes son función de las posiciones de los satélites, como veremos a continuación) se vuelve una matriz mal condicionada al tener vectores en su interior que son prácticamente unos combinaciones lineales de otros. Esto implica que al hacer la inversión de H necesaria para obtener la matriz G estaremos dividiendo por un determinante próximo a 0, con lo que la matriz G tendrá componentes con un valor elevado, lo que significa que se disparan los errores. Si todos los satélites se encuentran próximos unos a otros en el cielo la información que proporcionan es pobre y redundante, pero si están situados de forma espaciada proporcionan información suficiente para determinar la posición del receptor con buena precisión. Figura 5.3. Ilustración acerca del factor PDOP. Calidad del servicio de posicionamiento 56 56 Hasta ahora se ha hablado de direcciones para el cálculo de la posición, pero sin especificar en qué sistema de referencia. Los coeficientes , y pueden calcularse en cualquiera de ellos. Sin embargo si se calculan en, por ejemplo, el sistema Geocéntrico Ecuatorial, obtendríamos unas matrices y que nos darían un "HDOP" para el plano ecuatorial (que es el que habríamos usado). Sin embargo, si usamos el sistema de referencia Topocéntrico nos aseguramos de que el HDOP obtenido es para un plano que coincide con la superficie de la tierra para el observador. En algunas referencias [4], no hacen esto último, sino que calculan las matrices H y G directamente con los datos de las efemérides en el sistema Geocéntrico Ecuatorial y luego las transforman con matrices de rotación para poder obtener los factores HDOP y VDOP. En este proyecto se ha decidido emplear directamente el sistema de referencia Topocéntrico, ya que ahorra la necesidad de transformar matrices y facilita el cálculo de los coeficientes , y , como veremos a continuación. Si esto es así y llamamos a su tal y como en los apartados anteriores, tenemos que: iu i i u iu xx r y y zz         (5.31) Si el vector está expresado en el sistema de referencia Topográfico, mirando la Figura 5.1: cos(h)cos(Az) cos(h)sin(Az) sin(h) iu i i u iu xx r y y r zz                          (5.32) Por lo que los elementos , y se pueden calcular fácilmente como: 222 222 222 ( ) ( ) cos(h)cos(Az) ( ) ( ) ( ) ( ) ( ) cos(h)sin(Az) ( ) ( ) ( ) ( ) ( ) sin(h) ( ) ( ) ( ) i u i u xi i i u i u i u i u i u yi i i u i u i u i u i u zi i i u i u i u x x x x ar x x y y y y y y y y ar x x y y y y z z z z ar x x y y y y                             (5.33) Conocidos los valores de y (cuya obtención describimos en el apartado 5.1.1) de los N satélites que sean visibles en una época podemos calcular la matriz correspondiente como: 1 1 1 1 1 2 2 2 2 2 cos(h )cos(Az ) cos(h )sin(Az ) sin(h ) 1 cos(h )cos(Az ) cos(h )sin(Az ) sin(h ) 1 cos(h )cos(Az ) cos(h )sin(Az ) sin(h ) 1 cos(h )cos(Az ) cos(h )sin(Az ) sin(h ) 1 i i i i i N N N N N H               (5.34) Y de ahí mediante la ecuación (5.26). 57 57 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Efectivamente la matriz H es función de la geometría, y vemos que sus componentes son función de las posiciones de los satélites, por lo que como se comentó anteriormente si estos no se encuentran separados cubriendo bien el cielo, sus filas son prácticamente combinaciones lineales y la matriz está mal condicionada. 5.5 Parámetros de medida de la calidad del servicio La calidad del servicio que ofrece un sistema de posicionamiento por satélite no solo depende de la precisión que puede alcanzar, sino que se define en base a una serie de parámetros que en conjunto dan una idea de la calidad global que posee. Algunos de ellos son la disponibilidad, la capacidad, la continuidad, la cobertura y la integridad. La capacidad es el número de usuarios que pueden usar el sistema simultáneamente. En el caso de la constelación GPS es en principio ilimitada, ya que los satélites son pasivos, es decir, simplemente emiten sus correspondientes señales de radio, por lo que cualquier receptor que las reciba puede utilizar el sistema. La cobertura es la superficie o el espacio en general donde la calidad del servicio es suficiente para que un usuario pueda determinar su posición con una precisión determinada. En el caso de la constelación GPS la cobertura es idealmente global. La disponibilidad es el porcentaje de tiempo que el sistema es “utilizable”, dentro de su área especificada de cobertura. Esto significa que debe cumplir unos requisitos mínimos previamente especificados. Una definición típica de estos requisitos para la constelación GPS es que el usuario tenga un . Un ejemplo a nivel global en la Figura 5.3. Este concepto también puede definirse a nivel local, de modo que para una estación en tierra con unas coordenadas concretas se considere en qué momentos el servicio está o no disponible. La continuidad de un sistema de navegación respecto a una misión u operación se define como la probabilidad de que dicho sistema sea “utilizable” de forma continua durante toda la duración de dicha misión u operación. Esta definición llevada al caso que nos atañe equivaldría a por ejemplo imponer un límite máximo en el PDOP durante la misión. Por último, se define integridad de un sistema de navegación como la capacidad de dicho sistema para advertir que el sistema no debe ser utilizado (debido a que no está operativo o a que contiene errores). El sistema GPS no proporciona por sí mismo ningún tipo de mecanismo de integridad. Debido a que el sistema GPS es muy susceptible a fallos no esperados, y a la variabilidad de la precisión por los factores DOP, este sistema no es admisible para aplicaciones de navegación aérea, pero sí puede emplearse hibridado con otros sistemas como las IMU, además es ampliamente empleado en otros sectores como los vehículos no tripulados. Cálculo de la calidad del servicio 64 64 También podemos combinar (7.4) con la ecuación (6.15) para obtener la distancia angular en la ruta loxodrómica, que queda de la misma forma: j lox Vt R     (7.6) De nuevo los puntos de partida son conocidos: , , y . Con ello las ecuaciones (6.10) a (6.13), según corresponda, nos proporcionan el curso inicial. Ya podemos emplear las ecuaciones de (6.16) y (6.17) para calcular la posición de la aeronave en el instante , teniendo en cuenta que para cada época debemos resolverlas con el correspondiente valor de obtenido con (7.6). 7.1.3 Región Por último cabe plantearse cómo ampliar estos algoritmos al caso de una región. En primer lugar se discretiza una zona delimitada por unos valores mínimos y máximos de latitud y longitud definidos por el usuario. Figura 7.1. Discretización de Europa en 81 rectángulos. Dicha región se discretiza con un nivel de discretización seleccionable por el usuario. Para cada rectángulo se selecciona su punto medio, de modo que tenemos una lista de valores de latitud y una lista de valores de longitud. Si D es el número de filas y columnas: 1,2...,D; 1,2..., ij i j D   De modo que cada recuadro tiene definidos sus valores y . Con ello podemos realizar los análisis de número de satélites visibles y PDOP para cada punto de la discretización del mismo modo que lo realizamos para una estación fija en la superficie. Sin embargo, al tener dos índices para identificar a cada recuadro (i y j) seguimos necesitando otro índice para determinar a qué época pertenece la información. Es decir, el PDOP del recuadro el la época es: (t ) i j k PDOP 65 65 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Con esta notación no podemos almacenar la información en formato matricial ya que necesitamos tres índices para denotar un elemento. Para solucionarlo cambiaremos la notación a una notación de un índice : (i 1) D jI    (7.7) De esta manera el formato de los datos sería: 1 0 1 1 1 2 0 2 1 2 1 0 1 1 (t ) (t1) (t ) (t ) (t ) (t ) (t ) (t ) (t ) M M D D D D D D M PDOP PDOP PDOP PDOP PDOP PDOP PDOP PDOP PDOP PDOP             Donde el número de fila corresponde al índice . Para pasar a notación de nuevo, podemos hacer: I iD    (7.8) (i 1) DjI    (7.9) Donde la función   denota división con redondeo hacia el número entero mayor. Si ya tenemos los valores de PDOP para los puntos de la región en todas las épocas, podemos hacer cálculos como los mayores valores de PDOP en la región, mayores lapsos de tiempo sin disponibilidad sabiendo dónde se producen, calcular la disponibilidad para distintos valores máximos de PDOP, etc. Todas estas opciones se mostrarán en el capítulo siguiente. Equation Chapter (Next) Section 1 Interfaz de la herramienta 66 66 8 INTERFAZ DE LA HERRAMIENTA n el presente capítulo se va a realizar una explicación de la herramienta implementada con una interfaz en GUIDE, MATLAB. Veremos tanto la interfaz (desarrollada con ayudas de mathworks [17]) como las opciones que ofrece para modificar parámetros. Así mismo se explicará la información que proporciona al usuario como resultado de los análisis. 8.1 Representaciones 3D y paquete M_map Para que la herramienta quede completa es necesario disponer de medios para representar gráficamente las posiciones de los satélites, de modo que sea fácil apreciar las características de las órbitas y que el usuario disponga de un medio intuitivo para apreciar por dónde pasa cada satélite y sus posiciones, ya que esto está estrechamente ligado con la visibilidad y la calidad del servicio en un momento dado. Lo primero que necesitamos es la capacidad de representar las trazas de las órbitas, que son los puntos de la superficie sobre los cuales pasan los satélites. Para ello se ha empleado un paquete libre de mapeado, llamado M_map. Este paquete ofrece múltiples opciones, entre ellas la posibilidad de realizar mapas terrestres con diferentes tipos de proyecciones. A nosotros nos interesa la proyección de tipo cilíndrica equidistante, ya que permite representar toda la tierra en una superficie finita. Sin embargo este tipo de proyección deforma la imagen de la superficie, más notablemente cuanto más nos acercamos a los polos, de tal modo que estos puntos se transforman en rectas. Nótese que en la Figura 8.1 toda la línea con latitud 90º Norte corresponde al Polo Norte, mientras que el Polo Sur es la línea de latitud 90º Sur. Figura 8.1. Proyección Cilíndrica Equidistante realizada con M_Map. Sin embargo latitud y longitud siguen una escala lineal y uniforme en la figura. Una vez disponemos de una representación de la superficie, con las trazas empleamos las matrices y obtenidas en el capítulo 4 para representar las posiciones de los satélites sobre la Tierra. Dado que conocemos los valores de todos los satélites para cada época podemos representar tanto las curvas completas como las posiciones en una época dada. E 67 67 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Por otro lado, necesitamos una visualización 3D de las órbitas. En primer lugar realizamos una visualización de la Tierra en 3D a la que añadimos un mapa de los continentes [18]. A continuación, empleando las matrices también obtenidas en el capítulo 4, hacemos la representación en sistema Geocéntrico Ecuatorial de las órbitas. De nuevo al disponer de los datos para todas las épocas podemos hacer representación de las órbitas completas y de la posición de cada satélite en una época dada. Figura 8.2. Representación de dos órbitas circulares en el sistema de referencia Geocéntrico Ecuatorial. Otro objetivo interesante es la realización de animaciones, de modo que podamos observar la evolución de las posiciones de los satélites en sus órbitas y en sus trazas, aportando al usuaro información para conocer qué satélite es visible y en qué momento, ya que las animaciones incluyen información de la hora a la que corresponden las imágenes. El usuario podrá imponer un ángulo de máscara y la herramienta calcula la elevación de cada satélite, así que basta con comprobar la visibilidad de cada satélite durante la animación, de modo que se señale con algún tipo de signo aquellos satélites que sean visibles por el observador. Este signo será una línea que una al observador con aquellos satélites que están siendo visibles en un momento determinado. A cada satélite le corresponderá un color aleatorio, de modo que puedan distinguirse unos de otros, y al mismo tiempo cada uno llevará una etiqueta con su correspondiente número, que lo identifica y permite también saber cuál es cuál y diferenciar uno del resto. La línea de unión que empleamos para señalar la visibilidad también será del color correspondiente al satélite visible. Por otro lado, dado que la discretización en el tiempo puede ser del orden de un minuto para resultados de buena resolución, y teniendo en cuenta que el tiempo de la simulación es de un día, es recomendable poder acelerar la animación. Esta opción está disponible para el usuario como se explicará más adelante. Interfaz de la herramienta 68 68 8.2 Inicio de la herramienta La herramienta pretende ofrecer al usuario un análisis de la constelación GPS para la fecha que se desee. De modo que lo primero que necesita la herramienta es conocer la fecha que el usuario desea analizar. Una vez se conoce, mediante los procedimientos explicados en el apartado 4.2, la aplicación descarga y procesa de forma automática los datos correspondientes a dicha fecha. Para ello habría que pulsar (Figura 8.3) el botón “Descargar datos e iniciar”. Figura 8.3. Ventana de Inicio de la herramienta. También debemos considerar la posibilidad de que el usuario no disponga de acceso a internet pero ya posea los datos de la constelación para una fecha. Este caso se puede dar de hecho en el aula, si el equipo de alguno de los alumnos no dispone de conexión a internet. En dicho caso basta con que un compañero le proporcione el archivo .sp3 de la fecha correspondiente. Se indica al programa la fecha de la simulación y se pulsa la opción “Iniciar sin descarga”. Cuando la herramienta termina de obtener los datos de la constelación, abre automáticamente el módulo de análisis. Figura 8.4. Ventana del menú del módulo de análisis. 69 69 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Esta ventana da paso a los diferentes tipos de análisis que puede realizar la herramienta. En caso de querer seleccionar una fecha nueva, basta con pulsar la opción “Volver”. En este caso se abre la ventana de inicio y el usuario puede seleccionar una nueva fecha para el estudio. 8.3 Análisis de estación en tierra Figura 8.5. Ventana de análisis de una estación en tierra. En primer lugar contamos con una visualización 3D de las órbitas: Figura 8.6. Visualización 3D de las órbitas. Interfaz de la herramienta 70 70 Esta visualización nos permite por ejemplo observar los planos orbitales y qué distribución tienen los satélites en ellos. Dicha visualización puede rotarse para facilitar la visión de los números de los satélites y para poder ver las órbitas desde diferentes perspectivas. Al lado aparece una visualización de las trazas de la constelación, que representan los puntos de la superficie sobre los que pasan los satélites. Figura 8.7. Traza de la constelación GPS. Esta imagen nos permite observar sobre qué puntos pasa cada satélite, así como la posición exacta de cada uno de ellos a la hora que desee el usuario. Para esto existe una barra deslizante que permite al usuario modificar la hora. Esta se encuentra bajo la visualización 3D. Figura 8.8. Barra deslizante para modificar la hora de la visualización de satélites en las trazas. Del mismo modo, podemos observar en una imagen 3D la disposición de los satélites en sus órbitas para la hora que desee el usuario, así como una indicación de si son o no visibles con una línea de unión entre la estación y los satélites visibles. Para ello se fija la hora en la barra deslizante y se pulsa la opción “3D”. 71 71 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 8.9. Satélites visibles desde Sevilla el 12/12/2012 con ángulo de máscara de 10º a las 06:08. Dado que el objetivo de este apartado es realizar el análisis para una estación en tierra, se debe permitir al usuario la opción de introducir las coordenadas de la localización que desee. Figura 8.10. Personalización de parámetros de la estación y de la simulación. Como se ve en la Figura 8.10 el usuario puede introducir las coordenadas de la estación (por defecto aparecen las correspondientes a Sevilla, visibles en la imagen). Existe una lista desplegable para indicar las direcciones norte, sur, este y oeste. Interfaz de la herramienta 72 72 La localización de la estación en el mapa se indica con una cruz roja sobre un rombo negro. Figura 8.11. Señalización de la estación. Aunque se de la opción de introducir grados, minutos y segundos para dar mayor precisión, es perfectamente posible introducir valores como y dejar en blanco las demás casillas, lo cual sería equivalente, en ese ejemplo, a haber introducido y . Otro parámetro importante es el ángulo de máscara. Este parámetro se puede personalizar, como se ve en la Figura 8.10. En el capítulo 4 se explicó el interpolador programado en la herramienta. Como se dijo, tanto la discretización del tiempo como el grado de interpolación son seleccionables por el usuario. La discretización en el tiempo permite obtener mayor resolución en los resultados obtenidos. Para el grado de interpolación se recomiendan valores de 3 o 5. Para más información ver Apartado 4.3.1. Otra opción interesante es la información del número de satélites disponibles en una fecha. En el capítulo 4 explicamos cómo obtener la información de los satélites de la constelación. De entre los datos que la herramienta obtiene a partir de ellos está la lista de los satélites disponibles, que permite a la aplicación mostrar una lista de “checkbox” correspondientes a dichos satélites: Figura 8.12. Panel de satélites disponibles en la fecha. 73 73 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite En caso de que algún satélite no estuviera disponible no aparecería en la lista. Un ejemplo: Figura 8.13. Satélites de la constelación nominal. Sólo 24 de 32 disponibles. Este ejemplo ha sido obtenido pulsando la opción de “Cargar constelación nominal”, que te permite calcular la constelación con los 24 satélites nominales para poder realizar comparaciones entre la constelación real en una fecha y la nominal. En fechas actuales en general hay en torno a 32 satélites disponibles. Sin embargo si el usuario selecciona fechas años atrás sucede que hay menos satélites disponibles (porque estos no habían sido aún puestos en órbita). También puede ser que en una determinada fecha no haya algún satélite disponible por otras razones. En dicho caso el satélite desaparecerá de la lista. Para volver a cargar la constelación real en la fecha después de haber cargado la nominal, se puede pulsar el botón “Cargar constelación real”. Además de cargar los satélites disponibles, la herramienta ofrece la opción de desactivarlos, y realizar los análisis con dichos satélites fuera de servicio. De esta manera podemos observar de qué modo empeora por ejemplo el PDOP para una estación determinada según qué satélites estén o no activos. Otro parámetro importante incluido en este recuadro (Figura 8.12-13) es el de PDOP máximo. Más adelante se harán observaciones al respecto. Por último el panel de resultados da acceso a la representación del resultado de los análisis realizados: Figura 8.14. Panel de resultados para el análisis de una estación. El primer botón sirve para hacer efectivas las modificaciones impuestas en la configuración del análisis por el usuario (cambios en la localización de la estación, en el ángulo de máscara, satélites desactivados, PDOP máximo…). Interfaz de la herramienta 80 80 Si animamos las trazas de este análisis sucede lo siguiente: Figura 8.27. Visibilidad de los satélites en la animación de las trazas para Sydney. La figura pareciera dar información errónea. Hay líneas que enlazan la estación con satélites que parecen estar alejados, y la supuesta línea de visión atraviesa el globo de lado a lado para llegar de la estación al satélite. Sin embargo los satélites que muestra como visibles efectivamente lo son, ya que aunque parezcan estar “lejos”, realmente están cerca de la estación, pues los meridianos 180ºE y 180ºW son el mismo. El problema es que si únicamente indicamos a MATLAB que cree una línea de unión, esta cruza todo el mapa, ya que conecta puntos con longitudes negativas y puntos con longitudes positivas, con lo que las rectas de unión atraviesan el meridiano de Greenwich. Evidentemente MATLAB desconoce la continuidad que existe en los meridianos 180º, y que ese es el camino más corto. 81 81 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Para arreglarlo simplemente tenemos que hacer una corrección si el satélite visible en la época cumple que (t ) e i j     . El resultado obtenido entonces es el siguiente: Figura 8.28. Corrección de la animación para el caso de la Figura 8.23. Ahora las rectas representadas son más fieles a línea visual entre cada satélite y la estación. En estas figuras se puede observar también la indicación del tiempo, que aparece en la parte superior. En la Figura 8.28, la imagen corresponde a las 00:34. Si en la Figura 8.14 pulsamos la opción “Animar 3D” se obtiene una imagen del globo terráqueo y las órbitas de los satélites. En dicha visión 3D la Tierra rota, tal y como debe, respecto al sistema de referencia Geocéntrico Ecuatorial, que es en el que está hecha la animación, y se observan los satélites orbitando. En caso de que el usuario seleccione un ángulo de máscara, aparecerá un cono ubicado con su vértice en la estación, que marca la región del espacio sobre esta de la que tiene visibilidad. Así mismo, aquellos satélites que en una determinada época sean visibles por la estación estarán unidos a ella con una línea que en este caso sí que representa la línea real de visión. Es por esto que esta línea aparece para aquellos satélites que entren en el cono de visibilidad (región con elevación mayor que el ángulo de máscara) y desaparece cuando un satélite abandona dicha región y deja de ser visible. Interfaz de la herramienta 82 82 Figura 8.29. Animación 3D de la visibilidad de los satélites. El cono de visibilidad por supuesto ha de estar siempre situado con su vértice en la localización de la estación, por lo que se mueve solidario a la rotación de la Tierra. Además, como se aprecia en la Figura 8.29, los satélites aparecen numerados para permitir su identificación. Por último, hay una etiqueta en la que se muestra la hora correspondiente durante la animación, del mismo modo que sucede con la animación de la traza. Esta representación puede rotarse, ampliarse o reducirse como el usuario desee para facilitar una visualización lo más cómoda posible. 83 83 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 8.30. Animación 3D de las órbitas, con indicación de la visibilidad y hora. De esta manera, empleando las animaciones podemos observar de forma muy intuitiva qué satélites son visibles en cada época. Así mismo la representación gráfica del historial de visibilidad de cada satélite (Figura 8.21) permite ver de forma precisa e inmediata qué satélites son visibles y cuáles no, en qué épocas, etc. Además la barra deslizante (Figura 8.9) nos permite observar las posiciones de los satélites en sus trazas de forma estática para la época que desee el usuario, así como su disposición en el espacio pulsando la visualización 3D. Interfaz de la herramienta 84 84 8.4 Análisis de continuidad en vuelo Figura 8.31. Ventana de análisis de continuidad en vuelo. En primer lugar tenemos una visualización de las trazas de la constelación, así como la representación de la trayectoria que seguirá la posición de la aeronave. Figura 8.32. Traza de la constelación GPS y trayectoria del avión. La trayectoria se representa con una línea más gruesa de color azul. Los puntos de origen y destino se marcan con cruces azul claro y rojo respectivamente. La posición de la aeronave se representa con un rombo de color negro. 85 85 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 8.33. Visualización del vuelo. En la figura, entre Sevilla y Hamburgo. Esta imagen nos permite observar sobre qué puntos pasa cada satélite, así como la posición exacta de cada uno de ellos a la hora que desee el usuario. Lo mismo sucede con la trayectoria de la aeronave. Para esto existe una barra deslizante que permite al usuario modificar la hora, de modo que se puedan visualizar las posiciones de todos los elementos de forma estática en la época que se desee. Esta barra se encuentra a la derecha de la imagen de las trazas (Figura 8.31). Los puntos de origen y destino son modificables, así como la velocidad de vuelo. Figura 8.34. Modificaciones de los parámetros del vuelo. El tipo de ruta (Capítulo 6) es seleccionable con una lista desplegable: ortodrómica o loxodrómica. Por último, del mismo modo que se explicó en el apartado anterior, se puede seleccionar la constelación nominal y la constelación real. Un panel similar al de las Figuras 8.12-13 aparece en esta ventana, que nos permite desactivar o activar los satélites disponibles para personalizar el análisis. Interfaz de la herramienta 86 86 El valor de la discretización en el tiempo así como el grado de la interpolación se pueden editar al igual que en el análisis de la estación en tierra. Figura 8.35. Edición de parámetros para el interpolador. Editor de hora de partida para el vuelo. En esta figura observamos un nuevo panel, llamado “hora de partida”. Dado que simularemos un vuelo, la herramienta nos va a proporcionar la información de la disponibilidad durante el tiempo que dure, por lo que podemos seleccionar también la hora a la que este vuelo comienza. Las posiciones de los satélites durante el vuelo dependerán de la franja horaria durante la que transcurra. Fijada la hora de salida, el tiempo de llegada viene determinado por la velocidad del vuelo y la trayectoria a seguir. El panel de resultados tiene un aspecto similar al de la Figura 8.13. Figura 8.36. Panel de resultados para el análisis en vuelo. En la parte superior apreciamos la casilla de edición del PDOP máximo deseado por el usuario para los análisis de disponibilidad. La función de los botones “Número de satélites disponibles”, “PDOP” y “Visibilidad de los 87 87 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite satélites” es la misma que en el análisis para una estación fija. Sin embargo ahora los análisis se muestran para el tiempo del vuelo, durante el cual queremos conocer la disponibilidad. Por ejemplo, para el análisis Sevilla-Hamburgo de las figuras anteriores, con salida a las 04:27, usando la constelación nominal con 4 satélites menos e imponiendo un PDOP máximo de 2, obtenemos lo siguiente: Figura 8.37. Evolución del PDOP para un caso de vuelo. Figura 8.38. Análisis de la disponibilidad para el caso de la Figura 8.27 con PDOP máximo de valor 2. Interfaz de la herramienta 88 88 Por último, la opción “Animar trazas” de la figura 8.36 nos permite animar las trazas. La diferencia con el análisis en tierra es que la posición de la aeronave cambia con el tiempo. De nuevo, se nos da la opción de acelerar la velocidad de la animación. Figura 8.39. Animación de las trazas para el análisis en vuelo. 8.5 Análisis de una región Figura 8.40. Ventana de análisis de una región. 89 89 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite En esta ocasión no disponemos de un mapa con las trazas, ya que este puede consultarse fácilmente en los dos módulos anteriores, además, si se quiere saber con exactitud que satélites son más o menos tiempo visibles en la región, más adelante veremos que existe la opción de ver el historial de visibilidad de cada uno de ellos, como en los análisis en tierra y vuelo. En este modulo el mapa se muestra para permitir al usuario visualizar la región seleccionada. Figura 8.41. Mapa con la región de Europa seleccionada. Los límites de la región son seleccionables, así como el nivel de discretización del mismo. Como venimos comentando en los anteriores capítulos, la proyección utilizada para el mapa deforma la imagen. Esto significa que puntos con mayores latitudes aparecen ampliados en cuanto a longitud se refiere. Por lo tanto dichas zonas tienen una discretización más fina en cuanto a puntos sobre la superficie real, ya que la discretización se está realizando en valores de latitud y longitud, que sí guardan una escala lineal y uniforme en la proyección cilíndrica equidistante. Práctica de Fundamentos de Navegación Aérea 96 96 ser tedioso y confuso. Esto podía dificultar la práctica, alargarla, y desviar la concentración de los alumnos en los resultados debido a las configuraciones que de debían realizar entre cada análisis. A continuación un ejemplo de la configuración de STK para analizar acceso y cobertura para una región: Figura 9.3. Párrafos extraídos del boletín de la práctica para STK [19]. En la figura vemos una explicación clara y paso a paso del proceso a seguir para realizar los análisis necesarios para alcanzar los objetivos de la práctica con STK. Sin embargo se aprecia a simple vista la gran cantidad de manipulación de opciones y pasos necesarios. El boletín de la herramienta desarrollada en este proyecto busca ser más sencillo (dado que la herramienta está más directamente enfocada a cumplir los objetivos de la práctica). Se puede consultar en el Apéndice A. 97 97 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 9.3 Cuestionario a los alumnos Para valorar la opinión de los alumnos, comprobar la facilidad de uso y eficacia de la aplicación, y por supuesto recaudar propuestas de mejoras, recomendaciones, detección de fallos, y cualquier otra valiosa aportación por parte de los alumnos (que son el usuario objetivo de la herramienta) se elaboró un cuestionario para que los alumnos lo rellenaran al finalizar la práctica. El formato era el siguiente: Figura 9.4. Formato de cuestionario de opinión para los alumnos. Práctica de Fundamentos de Navegación Aérea 98 98 9.3.1 Resultados de la valoración de los alumnos A continuación vamos a analizar los resultados del cuestionario y mostrar los comentarios más importantes de los alumnos. El cuestionario fue respondido por 45 alumnos. Facilidad de uso de la aplicación Interés de la práctica Claridad del boletín Media 4,62 4,28 4,40 Desviación típica 0,61 0,82 0,65 En general los resultados obtenidos son bastante optimistas. Los alumnos están satisfechos con la facilidad de uso de la aplicación y la claridad del boletín. Lo mismo sucede con el interés de la práctica, que es sin embargo el punto más dado a variedad de opiniones, como se ve en su mayor valor de la desviación, ya que el interés en el tema de la práctica es algo más subjetivo. El punto que en general ha recibido mayor puntuación es la facilidad de uso. Teniendo en cuenta que los alumnos se enfrentaron a la práctica y la herramienta en un estado aun sin pulir (Apartado 9.1), los resultados obtenidos son satisfactorios. Algunos de los comentarios más importantes (por ser los más repetidos o por el contenido en sí) son: 1) Falta de espacio en el boletín para escribir las soluciones. Un alumno propone una tabla. 2) Falta de reloj que indique la hora en la animación. Ausencia de una opción para pausar la animación y ver la disposición de los satélites en un momento dado, o la posibilidad de manipular el tiempo de la animación. 3) En el módulo de región, no hay opción de hacer zoom para ver mejor los límites de la región. Además no se puede saber qué satélites sobrevuelan la región ya que no hay historial de visibilidad para cada satélite (aún no se había implementado para la región). 4) Al ampliar la visualización 3D en el panel de parámetros del análisis de estación en tierra, la imagen sale de su recuadro. Como se ve, gran parte de los comentarios hacen referencia a las limitaciones comentadas en el Apartado 9.1. Aparte de estos comentarios, había muchos valorando positivamente la herramienta, que en general se considera fácil e intuitiva. Consideran la práctica interesante y afirman que les ha ayudado a comprender los conceptos vistos en la asignatura, y que los resultados gráficos son fáciles de interpretar 99 99 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite A continuación se va a explicar cómo se han resuelto las observaciones hechas por los alumnos: 1) Este problema es de fácil solución. Por ejemplo elaborar una tabla como propone uno de los alumnos es una buena opción. En cualquier caso tan solo se debe editar el boletín. 2) Como se ha explicado en el capítulo 8, todas las visualizaciones de trazas u órbitas cuentan ya con un reloj, así como una numeración de los satélites para facilitar su identificación. Además, se ha implementado la opción de manipular el tiempo de las visualizaciones con la barra deslizante de tiempo, que permite ver trazas y órbitas en 3D en la época que desee el usuario de forma estática. Además se ha mejorado la animación 3D, que ahora también cuenta con numeración de los satélites e indicación de la visibilidad. 3) En el módulo de región se ha añadido la opción de hacer zoom para ver de cerca el mapa y los límites de la región. Además como se explicó en el capítulo 8, ahora se puede acceder fácilmente a la información de qué satélites son visibles desde la región. 4) Respecto a este punto, es una limitación propia de la interfaz GUIDE de MATLAB. Las observaciones de los alumnos han sido tenidas en cuenta en la mejora de la herramienta hasta su estado final, y se han implementado más mejoras explicadas en el Capítulo 8. Se espera por tanto que en el futuro se alcancen mayores niveles de satisfacción por parte de los alumnos. Análisis con la herramienta 100 100 10 ANÁLISIS CON LA HERRAMIENTA n este capítulo se van a realizar varios análisis con la herramienta, para intentar reproducir gráficas de interés de algunas referencias [3] y comparar resultados para probar la exactitud, lógica y coherencia de los resultados que proporciona la herramienta. 10.1 Comparación del análisis de estación La siguiente figura, obtenida con el software comercial STK, muestra un análisis local para Boston, correspondiente a la fecha 23/12/1994. Se ha empleado un ángulo de máscara de 7º, y se ha impuesto un PDOP máximo de 6. Estos datos y las coordenadas de la ciudad pueden consultarse en la esquina superior derecha de la figura. Figura 10.1 Análisis de la constelación nominal para Boston empleando STK. Vamos a intentar reproducir dichos análisis empleando la herramienta. Hacerlo es muy sencillo. En primer lugar acudimos al inicio de la herramienta y fijamos la fecha a 23/12/1994, y entramos al módulo de análisis de estación en tierra. Fijamos las coordenadas a 42.35ºN, 71.08ºW. Fijamos el ángulo de máscara a 7º, y establecemos un PDOP máximo de 6. E 101 101 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 10.2. Configuración en la herramienta para reproducir el análisis de la Figura 10.1. Solo queda pulsar “Actualizar datos”. Veamos los resultados obtenidos. Visibilidad de los satélites: Figura 10.3. Visibilidad de los satélites para el análisis de la Figura 10.2. Comparemos satélite a satélite con los resultados de la referencia. Se ha girado respecto a un eje horizontal la Figura 10.8 para que coincida el orden de los satélite, ya que nuestra gráfica los ordena en orden ascendente y la que proporciona el STK los ordena en orden descendente. Análisis con la herramienta 102 102 Figura 10.4. Comparativa de visibilidad de satélites entre STK y la herramienta desarrollada. Observamos que no hay apenas diferencia entre ambas gráficas, más allá de las separaciones verticales entre los datos de cada satélite ya que nuestra gráfica deja un espacio en blanco y la obtenida con STK elimina el espacio. Pero la información de disponibilidad es exactamente igual para todos los satélites, lo cual nos sirve para confirmar la veracidad de los análisis realizados por la herramienta desarrollada en el proyecto en cuanto a este tipo de análisis se refiere. Puede comprobarse también mirando los datos de cada satélite en las Figuras 10.6 y 10.8. Si esta gráfica es idéntica, también lo será la de “Número de satélites visibles”, pues es la suma para cada época. 103 103 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 10.2 Comparación de análisis de disponibilidad manipulando satélites en la constelación Ahora vamos a intentar hacer una comparación más interesante: la disponibilidad en una situación más desfavorable (desactivando satélites). Para ello emplearemos la siguiente figura: Figura 10.5. Análisis de la constelación nominal para Boston, con los satélites 16, 25 y 26 fuera de servicio. Ponemos a prueba una vez más la herramienta, no solo en cuanto al número de satélites visibles, sino en cuanto a las simulaciones de PDOP y disponibilidad con satélites desactivados (16, 25 y 26 fuera de servicio). Además en esta ocación el software STK detecta dos gaps en los que el servicio no está disponible: uno a las 12:00 y otro poco antes de las 21:00. Trataremos también de comparar la forma de PDOP en función del tiempo. Primero configuramos nuestra herramienta para la situación especificada: Figura 10.6. Configuración de la herramienta con los satélites 16, 25 y 26 desactivados. Análisis con la herramienta 104 104 Veamos los resultados obtenidos: De nuevo la gráfica de la visibilidad de los satélites es idéntica: se han eliminado los correspondientes 3 satélites. En cuando al PDOP y los gaps: Figura 10.7. Comparación de espacios sin disponibilidad según análisis con la herramienta y STK. La herramienta detecta también dos gaps. Uno centrado en las 12:00 y otro que finaliza a las 20:50 (justo antes de las 21:00). De nuevo la herramienta nos proporciona unos resultados idénticos a los obtenidos con el software comercial STK. Respecto al PDOP: Figura 10.8. Comparación de la evolución del PDOP con STK y con la herramienta. Las diferencias, de nuevo son mínimas. Los análisis y comparaciones realizados en este apartado nos dejan ver, en primer lugar, la precisión de los cálculos y análisis que realiza la herramienta. Además nos permiten confiar en las configuraciones realizadas con la herramienta: los ajustes realizados en las coordenadas de la estación, ángulo de máscara, satélites fuera de servicio, o límite de PDOP máximo nos han 105 105 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite dado unos resultados idénticos a los obtenidos con STK, a pesar del número de restricciones que estábamos imponiendo a la simulación. Además, se observa que reproducir este tipo de análisis es sumamente sencillo, ya que la herramienta proporciona una interfaz intuitiva para realizar las modificaciones. La lista con los satélites disponibles indicando el número correspondiente a cada uno ha permitido seleccionar correctamente aquellos que se debían desactivar. 10.3 Disponibilidad global Para continuar probando la calidad de los análisis de la herramienta vamos a intentar reproducir una nueva gráfica. Para ello se va a emplear la herramienta, más concretamente el módulo de análisis de región, para analizar la disponibilidad del servicio de posicionamiento a nivel global con la constelación nominal. El objetivo es recrear la siguiente gráfica: Figura 10.9. Disponibilidad global para distintod valores de PDOP. Para ello empleamos el módulo de análisis de región y discretizamos el mapa completo. Empleamos una discretización de 50 (2401 puntos), un “dt” de 5 minutos y grado de interpolaión 5 con la constelación nominal. Realizamos los cálculos para ángulo de máscara 0º, 2.5º, 5º y 7.5º. Las gráficas las obtenemos con la opción “Disponibilidad Región” en el panel de resultados, ya que proporciona los valores de la disponibilidad para diferentes PDOP. Constelación Galileo 112 112 Empleando el propagador para reproducir las órbitas y el interpolador para manipular la discretización del tiempo como hicimos con la aplicación GPS, obtenemos las efemérides de la constelación y las usamos para recrear una aplicación de análisis de la constelación nominal Galileo. Vamos a realizar varios análisis y comparar con la constelación GPS. Figura 11.1. A la izquierda la constelación Galileo. A la derecha la constelación GPS. Mientras que los datos obtenidos con el propagador nos dan planos exactamente coincidentes al haber sido obtenidos matemáticamente, los datos de la constelación GPS, obtenidos de la base de datos de IGS muestran órbitas aproximadamente en el mismo plano, ya que en la realidad hay pequeñas desviaciones y es imposible posicionar los satélites en los planos orbitales de forma exacta, además a lo largo del tiempo las perturbaciones modifican las órbitas, que se tienen que corregir. Las efemérides de las órbitas GPS en 1994 se encuentran mejor contenidas en los 6 planos correspondientes: Figura 111.2. Constelación GPS nominal en 1994. 113 113 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Figura 11.3. Arriba constelación Galileo. Abajo constelación GPS. Respecto a las trazas, vemos que tanto la constelación Galileo como la GPS consiguen una distribución que cubre toda la superficie (hasta las latitudes correspondientes de 55º y 56º). Sin embargo la constelación Galileo tiene una distribución mucho más uniforme que la GPS, debido a que esta última sigue la distribución no uniforme comentada anteriormente. Constelación Galileo 114 114 Para comparar los valores de PDOP empleamos, por ejemplo, un análisis de estación para Sevilla: Figura 11.4. PDOP con Galileo (arriba), GPS nominal (centro) y GPS real con 32 satélites (abajo). Los valores de PDOP son muy similares en ambas constelaciones. Tanto la constelación GPS nominal como la Galileo tienen valores de PDOP en torno a 1.5, aunque la Galileo es más “estable”, ya que al ser regular se repite el mismo patrón, mientras que la GPS muestra algunos picos con máximos más altos. Sin embargo teniendo en cuenta que la GPS nominal tiene 3 satélites menos, logra un rendimiento bueno en este caso. La constelación GPS real es la que tiene menores PDOP, lo cual es razonable ya que tiene mayor número de satélites. 115 115 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Para el número de satélites visibles: Figura 11.5. Satélites visibles con Galileo (arriba), GPS nominal (centro) y GPS real con 32 satélites (abajo). De nuevo tanto la constea lación Galileo como la GPS tienen valores similares para la nominal, aunque la constelación Galileo alcanza valores de 12 y no baja de 7, mientras que la GPS alcanza 6 en un determinado momento y no alcanza valores superiores a 11. Ambas constelaciones funcionan de forma similar, destacando la GPS si tenemos en cuenta que son 24 satélites. La constelación GPS real logra visibilidad de más satélites debido al mayor número de ellos en esta constelación. Por lo tanto, a nivel local en Sevilla ambas tienen un funcionamiento muy parecido. Pero este tipo de análisis puede cambiar notablemente el resultado dependiendo de la localización. Constelación Galileo 116 116 Para comparar con un enfoque más general, vamos a analizar globalmente empleando el módulo de análisis de región. Para un ángulo de máscara de 7.5º: Figura 11.6. Disponibilidad global en función de PDOP para Galileo (izquierda) y GPS (derecha) nominales. GPS nominal Galileo nominal Disponibilidad PDOP ~ 0 1.8 1.5 ~ 50 1.88 1.55 ~70 2 1.66 ~ 90 2.5 1.85 100 - 2.9 La constelación GPS nominal no alcanza disponibilidad del 100% (anque se queda en torno al 99%) para los valores considerados de PDOP, además vemos que tiene mayores valores de PDOP para la misma disponibilidad. De este análisis se desprende que la constelación Galileo tiene mejor calidad que la GPS nominal a nivel global. Figura 11.7. Gaps para GPS nominal, PDOP máximo 5. Para Galileo la disponibilidad es 100% (sin gaps). 117 117 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite Para un ángulo de máscara de 10º: Figura 11.8. Disponibilidad global en función de PDOP para Galileo (izquierda) y GPS (derecha) nominales. GPS nominal Galileo nominal Disponibilidad PDOP ~ 0 1.9 1.8 ~ 50 2 1.85 ~70 2.1 1.85 ~ 90 2.65 1.95 100 - 3.3 Evidentemente al aumentar el ángulo de máscara la constelación GPS nominal no alcanza disponibilidad del 100% para los valores considerados de PDOP, además vemos que sigue teniendo mayores valores de PDOP para la misma disponibilidad. Figura 11.9. Gaps para GPS nominal, PDOP máximo 5. Para Galileo la disponibilidad es 100% (sin gaps). Constelación Galileo 118 118 Con estos análisis podemos comprobar que la constelación Galileo nominal ofrece (idealmente) un sistema de posicionamiento con mayor precisión y disponibilidad, a nivel global, que el GPS nominal. Dicho esto, un análisis más interesante sería observar cómo se ve afectada la calidad del servicio en cada constelación ante fallo de satélites, sin embargo este tipo de estudio no es trivial, ya que depende del satélite que se desee desactivar, (sobre todo en la constelación GPS, que tiene una distribución irregular) y el efecto de estas modificaciones es más notable en regiones concretas que a nivel global. Pero los resultados obtenidos aquí nos permiten esperar buena calidad del servicio de Galileo en el momento en que esté disponible. Es importante también comentar que una de las principales características de Galileo es que será capaz de proporcionar información en tiempo real sobre la integridad del sistema [4], lo que permitirá tener más confianza en él. 119 119 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite 12 CONCLUSIONES n este el último capítulo se van a mostrar las conclusiones extraídas de la realización del trabajo: - Empleando los conocimientos adecuados, es posible crear herramientas de software libre capaces de proporcionar la misma calidad en los resultados que la obtenida por software comercial. La programación en MATLAB y sus paquetes de creación de interfaz (GUIDE) ofrecen buenas posibilidades para ello. - La principal ventaja de las herramientas de este tipo es que están totalmente abiertas a modificaciones por parte del usuario. Se pueden añadir funcionalidades en un futuro y mejorar o modificar las actuales, así como cambiar la forma de mostrar los resultados. - Una vez desarrollados los algoritmos, también es crucial trabajar la implementación, teniendo en cuenta las limitaciones de la interfaz y de la plataforma en que se esté programando (en este caso MATLAB). - Con el empleo de la herramienta se pueden observar claramente conceptos relacionados con el servicio de posicionamiento, como la dependencia del PDOP con el número de satélites: Cuando aumenta el número de satélites disminuye el PDOP, si el número de satélites visibles disminuye, aumenta el PDOP. Como se aprecia en la figura, los “escalones” coinciden. E Conclusiones 120 120 - A la hora de obtener información fundamental, como en este caso han sido las efemérides de la constelación, podemos optar por la simulación o por tomarla de una base de datos, lo que nos garantiza unos datos más cercanos a la realidad. Sin embargo con esto último no basta. En ocasiones debemos combinar las fuentes con otras herramientas para obtener información más útil para la aplicación que deseamos alcanzar, como fue en este caso la interpolación para aumentar la resolución en el tiempo de los datos. - De nuevo, relacionado con los datos obtenidos de las efemérides desde el servicio IGS, es importante interpretar correctamente la información de la fuente y adaptarla a los formatos con los que se ha decidido trabajar, creando los algoritmos que nos permitan pasar de unos formatos a otros. - La geometría de las órbitas y la distribución de los satélites en ellas es realmente un aspecto fundamental a la hora de diseñar una constelación de navegación por satélites. Dicha geometría afecta tanto a la calidad del servicio como a la vulnerabilidad ante fallos de satélites. Para esto es importante tener en cuenta tanto el servicio a nivel local (estación) como global (región). - El empleo de este tipo de herramientas puede ser muy útil para ayudar a los alumnos en su aprendizaje, como confirman los resultados de las encuestas realizadas. Además la evaluación y crítica de la herramienta que realizaron los alumnos tras práctica ha sido fundamental a la hora de mejorarla y de buscar una interfaz lo más útil e intuitiva posible. - Actualmente se sigue avanzando en los sistemas de navegación por satélite. Finalizada la constelación Galileo, si se garantiza la interoperabilidad con otras constelaciones (capacidad de emplear varios sistemas juntos para alcanzar mejores resultados) mejorará el servicio y quizá su uso pueda extenderse a aplicaciones que requieran mayor precisión. Además, la capacidad de proporcionar información de la integridad de Galileo permitirá aumentar la confianza en este sistema. 12.1 Trabajo futuro En este proyecto se han tratado conceptos básicos relacionados con constelaciones de navegación por satélite y se ha implementado una herramienta que permita cumplir los objetivos definidos en el inicio de este proyecto. Dichos objetivos se han alcanzado satisfactoriamente, pero cabe plantearse nuevas posibilidades para mejorar y continuar la herramienta con trabajo futuro. - Implementar ángulo de máscara variable con el azimut para permitir análisis más avanzados, como se explica en el apartado 7.1.1. - Añadir selección de otros factores DOP como HDOP o VDOP. - Mejorar la propagación considerando las perturbaciones. - Eliminar la hipótesis de Tierra esférica implementando un modelo como puede ser el WGS-84. - Mejorar el análisis de vuelo eliminando las hipótesis aplicadas (velocidad constante, altitud nula…). - Del mismo modo que se ha implementado la constelación Galileo, extender la herramienta para analizar más casos, como el aumento de constelaciones con satélites extra geoestacionarios, análisis de otras constelaciones como GLONASS, etc. - En relación con el punto anterior, añadir análisis de la interoperabilidad, pudiendo analizar combinaciones de varias constelaciones. - Añadir nuevas gráficas y opciones de muestra de resultados de los análisis. Un ejemplo pueden ser los “sky plots”, que muestran la posición de los satélites en el cielo desde el punto de vista del observador, o la señalización en el mapa de la cobertura de cada satélite. - Estudiar la aplicación a sistemas como el GPS diferencial. 121 121 Herramenta para el cálculo de precisión y coberturas de un sistema de navegación por satélite APÉNDICE A. BOLETÍN DE PRÁCTICAS FUNDAMENTOS DE NAVEGACIÓN AÉREA Práctica 5: Disponibilidad global y regional del sistema GPS. Continuidad. Nombre del alumno: En esta práctica se pretende estudiar los conceptos de disponibilidad y continuidad del sistema GPS. Evaluación: Se pide rellenar el boletín de prácticas con ciertos resultados observados u obtenidos de la aplicación de análisis en MATLAB. Introducción: Conceptos de PDOP, disponibilidad y continuidad del sistema GPS. El empleo de la constelación GPS en la navegación aérea se basa en la posibilidad de determinar la posición de un observador conocidas las posiciones de los satélites y empleando señales de radio. Por lo tanto, para que el sistema funcione correctamente necesitamos un número suficiente de satélites y que estos sean visibles por el observador. En general, si un satélite está bajo el horizonte implica que el observador no puede recibir su correspondiente señal. Pero considerar un satélite visible desde que se encuentra en el horizonte es poco realista, ya que en general los elementos del entorno (edificios, accidentes geográficos) impiden ver los satélites hasta que estos se eleven un poco más en el cielo (concepto de ángulo de máscara). Esto nos determina el número de satélites visibles por el observador en un determinado momento en el tiempo. Pero aparte de tener visibles un número suficiente de satélites, necesitamos que la información que estos proporcionan sea útil para determinar nuestra posición con suficiente precisión.