Predicción de distancia mínima entre aeronaves a partir de datos ADS-B

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster Universitario en Ingeniería Aeronáutica Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Autor: Antonio Rivas Rincón Tutor: Alfonso Valenzuela Romero Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 i Trabajo Fin de Máster Máster Universitario en Ingeniería Aeronáutica Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Autor: Antonio Rivas Rincón Tutor: Alfonso Valenzuela Romero Profesor Contratado Doctor Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 ii iii Trabajo Fin de Máster: Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Autor: Antonio Rivas Rincón Tutor: Alfonso Valenzuela Romero El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal iv v A mi eterna amiga María. Por todos estos momentos que nos arrebataron. vi vii Agradecimientos Una vez llegado a este punto en la vida es inevitable echar la vista atrás y recordar, con cierta nostalgia, el camino seguido y todas las personas y vivencias que han formado parte de él y que, por tanto, forman parte de mí. De todos ellos he aprendido algo, de todos ellos me llevo algo y, sin duda, sin todos ellos hoy no sería la persona que soy ni estaría afrontando este reto al que hoy me enfrento. Es entonces, durante este recuerdo, cuando soy consciente de la cantidad de amigos, compañeros y profesores que deberían tener su reconocimiento en este humilde y tímido agradecimiento. Volviendo de nuevo al pasado, allá por febrero de 2018 comienza mi último período como estudiante universitario y salen a relucir las ganas por acabar junto con el nerviosismo por hacer frente al trabajo de fin de máster. Y, en este revuelo compuesto también por las horas en el trabajo y la falta de tiempo, recuerdo lo oscuro y confuso que se me hacía este proyecto. No sabía por donde empezar ni mucho menos cómo acabar. Sin embargo, a medida que avanzaba en este proceso y me adentraba en el túnel la luz cada vez se hacía más clara debido, principalmente, a la ayuda, guía y colaboración que en todo momento me ha prestado Alfonso Valenzuela Romero. Ha sido enorme la paciencia y desenfado con la que siempre recibía cualquiera de mis dudas o dificultades, a pesar de haber sido un proceso muy largo en los que hemos tenido momentos en los que el avance era muy lento. Sin embargo, y a pesar de estos momentos de agobio, siempre me hacía ver una idea tan clara como necesaria en la vida: “si vas a hacer algo, hazlo bien”, idea que transmite a todos sus alumnos. Gracias a mis padres por la educación que me han dado desde que nací, por enseñarme a mirar el mundo con una sonrisa, sea cual sea la sitacuón, y por el permanente apoyo que me han mostrado desde que tengo uso de razón. Fueron ellos los que me abrieron las puertas del mundo de la ingeniería aeroespacial, un mundo en el que una vez que entras, te atrapa y te cambia la vida. Es imposible no acordarme de mis compañeros durante estos cuatro años, mis fieles compañeros de batalla sin los cuales estos años hubieran sido traumáticos. Son ellos las personas capaces de transformar el día anterior a un examen en bromas y risas, esas personas en las que, aunque el agobio esté presente, siempre van a hacerte más leve cualquier trabajo interminable o cualquier práctica ininteligible, y quienes nunca dudan en ayudarte. Sin duda, su amistad es uno de los mayores regalos que me llevo. Y por último, pero no por ello menos importante, agradecer a todas esas personas que durante el desarrollo de este trabajo han tenido que soportar, involuntariamente, la descripción cada paso que daba, transformando cada uno de nuestros encuentros en una pequeña sesión de desahogo de mi agobio. Antonio Rivas Rincón Sevilla, 2020 xiv 8 Bibliografía 53 Apéndice A - Derivadas Parciales 55 A.1. Matriz Ak-1 del Filtro de Kalman extendido 55 A.2. Matriz Wk-1 del Filtro de Kalman extendido 55 A.3. Gradiente del curso 56 A.4. Gradiente de la velocidad respecto a tierra 57 A.5. Gradiente de la separación entre dos aeronaves 57 xv ÍNDICE DE TABLAS Tabla 2.1. Descripción de la estructura seguida en el archivo “.bst” para almacenar la información recibida 5 Tabla 4.1. Categoría para la precisión de la posición horizontal 16 xvi xvii ÍNDICE DE FIGURAS Figura 2.1. Receptor SBS-3 4 Figura 2.2. Ejemplo para ilustrar la estructura de la información almacenada en el archivo “.bst” 4 Figura 3.1. Resumen de las etapas del Filtro de Kalman discreto 10 Figura 3.2. Resumen de las etapas del Filtro de Kalman extendido 12 Figura 4.1. Resumen del Filtro de Kalman aplicado 17 Figura 4.2. Trayectoria rectilínea completa y vista en detalle. 22 Figura 4.3. Evolución de la desviación típica de la posición para una trayectoria rectilínea 23 Figura 4.4. Comparación entre el curso y rumbo geográfico para una trayectoria rectilínea 24 Figura 4.5. Desviación típica del curso geográfico para una trayectoria rectilínea 25 Figura 4.6. Comparación entre la velocidad respecto a tierra dada por ADS-B y la calculada a partir de la posición para una trayectoria rectilínea 25 Figura 4.7. Desviación típica de la velocidad respecto a tierra para una trayectoria rectilínea 26 Figura 4.8. Ejemplo de trayectoria con viraje 27 Figura 4.9. Ejemplo de trayectoria con viraje sin aplicar la corrección del intervalo temporal de evaluación 27 Figura 4.10. Evolución de la desviación típica de la posición para una trayectoria con viraje 28 Figura 4.11. Representación ampliada del efecto del viraje sobre la desviación típica de la posición 29 Figura 4.12. Comparación entre el curso y rumbo geográfico para una trayectoria con viraje 29 Figura 4.13. Desviación típica del curso geográfico para una trayectoria con viraje 30 Figura 4.14. Comparación entre la velocidad respecto a tierra dada por ADS-B y la calculada a partir de la posición para una trayectoria con viraje 31 Figura 4.15. Desviación típica de la velocidad respecto a tierra para una trayectoria con viraje 31 Figura 5.1. Ejemplo de trayectoria predicha 36 Figura 5.2. Evolución de la desviación típica de la posición predicha 37 Figura 5.3. Evolución de la desviación típica para una predicción de 20 minutos 38 Figura 5.4. Figura ilustrativa de cómo se obtiene la distancia mínima predicha 38 Figura 5.5. Ejemplo representativo de la evolución de la separación entre aeronaves con el tiempo 39 Figura 5.6. Ejemplo en el que una pareja de aeronave presenta dos distancias mínimas relativas 39 Figura 5.7. Ejemplo representativo de la predicción de la distancia mínima 40 Figura 5.8. Ejemplo representativo de la predicción de la distancia mínima cuando se tienen diferentes tramos de aproximación 40 Figura 5.9. Representación de la trayectoria para una pareja de aeronaves con diferentes tramos de aproximación 41 Figura 6.1. Error de predicción según la antelación 44 Figura 6.2. Desviación estándar del error de predicción según la antelación 44 Figura 6.3. Representación de las trayectorias para la pareja 4CA7B4 – 4CA8EA 45 xviii Figura 6.4. Representación de los puntos de mínima separación con error normalizado alto 46 Figura 6.5. Evolución de la separación para la pareja 4CA7B4 – 4CA8EA 46 Figura 6.6. Error de predicción para antelación menor a 600 segundos 47 Figura 6.7. Desviación estándar del error de predicción para antelación menor a 600 segundos 47 Figura 6.8. Trayectorias de la pareja de aeronaves 020098 - 020073 48 Figura 6.9. Evolución de la separación para la pareja 020098 – 020073 48 Figura 6.10. Error de predicción para antelación menor a 600 segundos tras corrección 49 Figura 6.11. Región de confianza del 95% para el error de predicción normalizado 50 xix Glosario/Terminología Aeronave rastreada Aeronave cuyos mensajes ADS-B son recibidos por la estación receptora en tierra. El tiempo transcurrido desde la recepción del último mensaje debe ser inferior al tiempo límite de rastreo. Tiempo límite de rastreo Tiempo transcurrido desde el instante en el que se recibe el último mensaje de una aeronave rastreada a partir del cual, si no se recibe ningún mensaje nuevo, la aeronave pasa de “rastreada” a “perdida”. Aeronave perdida Aeronave para la cual el tiempo transcurrido desde la recepción de su último mensaje ADS-B es superior al tiempo límite de rastreo. Posición mejorada Posición geográfica obtenida como resultado de la aplicación del Filtro de Kalman. Intervalo temporal de evaluación Diferencia de tiempo entre dos posiciones de una aeronave a partir de las cuales se estima su curso geográfico y velocidad respecto a tierra. Instante de predicción Instante de tiempo en el que se realiza la predicción de la distancia mínima para una pareja de aeronaves. Instante de distancia mínima Instante de tiempo en el que se alcanza la separación mínima entre una pareja de aeronaves. Esta separación entre dos aeronaves puede ser de dos tipos: - Predicha: pronosticada por el modelo de predicción. - Real: obtenida como resultado del análisis a posteriori de la separación entre la pareja de aeronaves a lo largo del tiempo. Error de predicción Diferencia entre la distancia mínima predicha y la distancia mínima real. Error de predicción normalizado Error de predicción normalizado respecto a la distancia mínima real. Se da en porcentaje. 1 1 INTRODUCCIÓN 1.1 Motivación y Objetivo Aunque antiguamente la información relacionada con la navegación aérea y las aerolíneas se consideraba de uso privado y era sólo accesible para aquellos implicados, con el paso del tiempo se ha experimentado una cierta liberalización de dicha información. Además, como consecuencia del desarrollo de la navegación aérea, la información disponible en abierto es cada vez es más extensa. De hecho, hoy en día es posible obtener una gran cantidad de información relacionada con el estado de una aeronave en vuelo, así como con el tráfico aéreo, ya sean trayectorias realizadas por aeronaves en el pasado, planes de vuelo planificados o, incluso, información del estado de la aeronave en tiempo real. Por otra parte, es un hecho contrastado el aumento de la demanda del espacio aéreo por parte de aerolíneas, por lo que los proveedores de servicios se esfuerzan en mejorar la oferta y, así, intentar satisfacer estas necesidades. Algunas de las medidas que los proveedores de servicios de navegación aérea (ANSPs, “Air Navigation Service Providers”) están desarrollando se centran en el concepto de “Free Route”. Al aumentar la complejidad de la red aérea se hace más complejo adelantarse a una posible violación de las distancias mínimas de separación. En este sentido, una mejora sería tener en cuenta el efecto del viento a la hora de predecir las trayectorias de las aeronaves. Recientemente se ha desarrollado una mejora de la comunicación Mode-S en la que las aeronaves adecuadamente equipadas pueden transmitir información meteorológica, denominada “Meteorogical Routine Air Report” (MRAR). Poder disponer de información meteorológica captada por las aeronaves implica la posibilidad de crear una gran red de estaciones meteorológicas distribuidas por todo el territorio y a diferentes altitudes. De esta manera, se abre una vía para poder conocer mejor el estado y la evolución de la atmósfera, uno de los factores más importantes en la navegación aérea debido a sus efectos sobre la aeronave y su carácter impredecible. Este proyecto surge inicialmente con la intención de incorporar los datos de vientos obtenidos gracias al MRAR para mejorar la predicción de la posición de las aeronaves y deducir la importancia de este factor, pues, entre otros factores, el viento es responsable de la degradación de la predicción de la separación con otras aeronaves conforme el horizonte temporal aumenta. Para ello se dispone del receptor SBS-3, capaz de recibir señales Mode-S/ADS-B con un alcance de hasta 400 km. Sin embargo, durante la investigación inicial se encontraron 3 problemas para alcanzar este objetivo: 1. Para recibir la información MRAR se debe interrogar al transpondedor a bordo solicitando dicha información. Sin embargo, el receptor SBS-3 no tiene capacidad de interrogación, sólo de escucha. 2. Como continuación del inconveniente anterior, en el momento en el que se inició este proyecto España no disponía de la infraestructura suficiente para interrogar mediante Mode-S EHS en la zona sur. Sin embargo, se espera que para 2020 se tenga capacidad para recibir datos Mode-S EHS en todo el país [10]. 3. Por tanto, ante esta situación ha sido imposible conocer qué información Mode-S es capaz de decodificar el receptor SBS-3. Ante esta problemática se optó por darle un giro al objetivo del proyecto. De esta manera, el nuevo y definitivo objetivo consiste en desarrollar un predictor de trayectorias para el corto plazo (hasta un horizonte de 20 minutos) basado en medidas ADS-B. Con este predictor se pretende emular las capacidades de predicción que tienen los proveedores de servicios de navegación aérea y analizar su degradación con el horizonte de predicción, especialmente en cuanto a la predicción de la distancia mínima entre aeronaves. 2 2 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 1.2 Metodología Para alcanzar el objetivo mencionado anteriormente se ha desarrollado un proceso en el que, partiendo únicamente de la posición de la aeronave transmitida mediante ADS-B, sea capaz de obtener la predicción de la separación mínima entre las aeronaves que sobrevuelan el espacio aéreo analizado. Este proceso consta de dos fases claramente diferencias: 1. Filtro de Kalman. Se diseña un Filtro de Kalman extendido con el objetivo de mejorar la posición recibida inicialmente mediante ADS-B. 2. Predicción de trayectorias futuras. A partir de la posición mejorada se estima la trayectoria futura para cada una de las aeronaves. De esta forma, se puede obtener fácilmente la distancia mínima entre aquellas aeronaves que se estén aproximando entre sí. Con la intención de poder evaluar la calidad de las predicciones y las separaciones obtenidas se realizará un análisis a posteriori en el que se comparen estos resultados con los valores obtenidos en la realidad. 1.3 Estructura del documento Este documento consta de 8 capítulos. El primero de ellos pretende presentar el objetivo de este proyecto. En el capítulo 2 se hace una breve introducción a la situación de partida y al receptor usado para este proyecto. Por otro lado, en el capítulo 3 se describe de manera teórica las dos modalidades del Filtro de Kalman, el discreto y el extendido para, en el capítulo 4, dar paso a la caracterización de las ecuaciones que definen el Filtro de Kalman extendido aplicado en este proyecto. El siguiente punto, el capítulo 5, consiste en la descripción del proceso seguido para estimar la trayectoria futura de las aeronaves y, a partir de ella, la separación entre las aeronaves circundantes. Finalmente, en el capítulo 6 se muestran los resultados. En ellos se realiza un análisis a posteriori en el que se pretende ilustrar la calidad de las predicciones tomando como referencia la distancia real entre las aeronaves. Para ello, se representa la evolución de la diferencia entre predicción y realidad según la antelación de las predicciones. 3 2 DATOS DE VIGILANCIA Desde el inicio de la aviación comercial, ésta ha experimentado un crecimiento prácticamente constante. Este crecimiento ha obligado al desarrollo de tecnología que pueda garantizar una gestión del tránsito aéreo de manera segura, económica, y eficiente. Como consecuencia de esta necesidad, se han desarrollado paulatinamente tres tipos de vigilancia dependiendo de la tecnología disponible: • Vigilancia Independiente No Cooperativa: provista por los radares primarios (PSR), los cuales sólo pueden determinar la posición de las aeronaves circundantes a partir del reflejo de la señal emitida. • Vigilancia Independiente Cooperativa: asociada a los radares secundarios (SSR). Éstos hacen posible la comunicación gracias a unos elementos a bordo denominados “transpondedores” capaces de responder a interrogaciones desde tierra. La posición horizontal de la aeronave es determinada por los sistemas en tierra mientras que la aeronave proporciona en su respuesta la altitud y su identidad. • Vigilancia Dependiente Cooperativa: basado en la tecnología ADS-B (“Automatic Dependent Surveillance – Broadcast”), mediante la cual la aeronave difunde su posición, altitud, velocidad, identidad, así como otros parámetros adicionales, a todos los usuarios que estén escuchando. Esta difusión se realiza automáticamente, a una tasa aproximada de 1 Hz, y no como respuesta a una interrogación previa. La Vigilancia Independiente Cooperativa se fundamenta en el uso de la técnica de interrogación “Mode-S”. Esta técnica de interrogación pretende evitar errores en la identificación de las aeronaves en espacios muy congestionados mediante la asignación de una dirección única a la aeronave, lo cual permite el establecimiento de una comunicación selectiva. A su vez, el Mode-S EHS (“Enhanced Surveillance”) es capaz de transmitir otros parámetros de la aeronave, como por ejemplo la altitud, rumbo magnético, velocidad respecto a tierra, IAS o Mach, entre otros. Aquellos radares equipados con Mode-S EHS pueden, a su vez, interrogar a las aeronaves solicitando información adicional. Esta información depende de dos factores, denominados “Downlink Format” (DF) y “Binary Data Store” (BDS). El BDS se trata de un mensaje incluido en la interrogación que realiza el radar secundario al transpondedor y que define el tipo de información que quiere recibir. Dentro de los diferentes códigos BDS se ha introducido recientemente el MRAR, del que se había hablado en el capítulo anterior y que corresponde al BDS 4,4. Este código permite obtener observaciones directas del viento y la temperatura que experimenta la aeronave en cada instante. Debido a la imposibilidad tecnológica para disponer de esta información, este proyecto se basa en los datos obtenidos mediante ADS-B, aunque se plantea como posible desarrollo futuro continuar con esta línea de investigación. 2.1 Datos de vuelo Para que la señal transmitida por la aeronave se traduzca en información útil para el resto de usuarios son necesarios dos factores adicionales: una antena receptora y un software decodificador. Receptor Como se comentaba anteriormente, para este proyecto se ha hecho uso del receptor diseñado por la compañía KINETIC denominado “SBS-3”. Este dispositivo es especialmente interesante por tratarse de un receptor “stand-alone”, es decir, es capaz de realizar el procesamiento y demodulación de numerosas señales de audio y datos de manera simultánea e interna. Sus principales características son: 10 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝐾𝑘=𝑃𝑘−𝐻𝑘𝑇(𝐻𝑘𝑃𝑘−𝐻𝑘𝑇+𝑅𝑘)−1 (3-15) 𝑥𝑘=𝑥𝑘−+𝐾𝑘(𝑧𝑘−𝐻𝑘𝑥𝑘−) (3-16) 𝑃𝑘=(𝐼−𝐾𝑘𝐻𝑘)𝑃𝑘− (3-17) A la vista de estas ecuaciones, una vez obtenida la Ganancia de Kalman, 𝐾𝑘 se incorpora el valor de la medida recibida, 𝑧𝑘, para poder obtener la estimación del estado a posteriori. Finalmente, se computa la estimación de la covarianza a posteriori. Una vez se completan las dos etapas el proceso se repite tomando la estimación a posteriori para predecir la estimación a priori del siguiente instante. Este aspecto es el que le confiere un carácter más práctico que otros filtros. En la Figura 3.1 se muestra un resumen del proceso a seguir por el Filtro de Kalman discreto para cada instante de evaluación. 3.2 Filtro de Kalman extendido Cuando el proceso viene descrito por un sistema de ecuaciones diferenciales no lineales se emplea el Filtro de Kalman Extendido. Este filtro consiste, simplemente, en un Filtro de Kalman que se linealiza en torno al valor del estado y del control en el instante 𝑡𝑘−1. En este caso, el vector de estado se modela como: 𝑥𝑘=𝑓(𝑥𝑘−1,𝑢𝑘−1,𝑤𝑘−1) (3-18) donde 𝑓(𝑥𝑘−1,𝑢𝑘−1,𝑤𝑘−1)∈ℜ𝑛 corresponde a un vector de funciones no lineales que relaciona el estado en el instante 𝑡𝑘−1 con el estado en el instante 𝑡𝑘. Por otro lado, la medida viene dada por: Etapa de predicción 1. Proyección del estado: 𝑥𝑘−=𝐴𝑥𝑘−1+𝐵𝑢𝑘−1 2. Proyección de la covarianza a priori: 𝑃𝑘−=𝐴𝑃𝑘−1𝐴𝑇+𝑄𝑘 Etapa de actualización 1. Cálculo de la Ganancia de Kalman: 𝐾𝑘=𝑃𝑘−𝐻𝑘𝑇(𝐻𝑘𝑃𝑘−𝐻𝑘𝑇+𝑅𝑘)−1 2. Actualización de la estimación: 𝑥𝑘=𝑥𝑘−+𝐾𝑘(𝑧𝑘−𝐻𝑘𝑥𝑘−) 3. Actualización de la covarianza: 𝑃𝑘=(𝐼−𝐾𝑘𝐻𝑘)𝑃𝑘− Estimación inicial para 𝑥0− y 𝑃0 Figura 3.1. Resumen de las etapas del Filtro de Kalman discreto 11 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝑧𝑘=ℎ(𝑥𝑘,𝑣𝑘) (3-19) Donde ℎ(𝑥𝑘,𝑣𝑘)∈ℜ𝑚 es definido por el vector de funciones no lineales que relaciona el estado, 𝑥𝑘, con la medida, 𝑧𝑘. En la práctica no es posible saber los valores del ruido 𝑤𝑘 y 𝑣𝑘. Sin embargo, al igual que se hizo para el Filtro de Kalman discreto, es usual suponer que estas variables siguen una distribución normal con media nula y covarianzas 𝑄𝑘 y 𝑅𝑘, respectivamente: 𝑤𝑘 ~𝑁(0,𝑄𝑘) (3-20) 𝑣𝑘 ~ 𝑁(0,𝑅𝑘) (3-21) Reescribiendo las ecuaciones para el sistema y la medida una vez linealizadas se tiene: 𝑥𝑘≈𝑥𝑘+𝐴𝑘−1(𝑥𝑘−1−𝑥𝑘−1)+𝑊𝑘−1𝑤𝑘−1 (3-22) 𝑧𝑘≈𝑧𝑘+𝐻𝑘(𝑥𝑘−𝑥𝑘)+𝑉𝑘𝑣𝑘 (3-23) donde: - 𝑧𝑘: estimación de la medida en el instante 𝑡𝑘. - 𝐴𝑘−1: jacobiano de la función 𝑓 respecto a 𝑥. La componente 𝑖,𝑗 de la matriz viene dada por: 𝐴[𝑖,𝑗]𝑘−1=𝜕𝑓[𝑖] 𝜕𝑥[𝑗](𝑥𝑘−1,𝑢𝑘−1,0) (3-24) - 𝑊𝑘−1: jacobiano de la función 𝑓 respecto a 𝑤𝑘−1. La componente 𝑖,𝑗 de la matriz viene dada por: 𝑊[𝑖,𝑗]𝑘−1=𝜕𝑓[𝑖] 𝜕𝑤[𝑗](𝑥𝑘−1,𝑢𝑘−1,0) (3-25) - 𝐻𝑘: jacobiano de la función ℎ respecto a 𝑥. La componente 𝑖,𝑗 de la matriz viene dada por: 𝐻[𝑖,𝑗]𝑘=𝜕ℎ[𝑖] 𝜕𝑥[𝑗](𝑥𝑘) (3-26) - 𝑉: jacobiano de la función ℎ respecto a 𝑣. La componente 𝑖,𝑗 de la matriz viene dada por: 𝑉[𝑖,𝑗]𝑘=𝜕ℎ[𝑖] 𝜕𝑣[𝑗](𝑥𝑘,0) (3-27) Es importante observar que, tras linealizar el sistema, el error entre el estado real y el linealizado varía. Se tiene que: - 𝜀𝑘: error entre el estado real del sistema y el estado linealizado. - 𝜂𝑘: error entre la observación real y la observación linealizada. Es decir, antes de linealizar el ruido venía dado por 𝑤𝑘 y 𝑣𝑘. Sin embargo, como resultado de linealizar, el ruido asociado al sistema y a la observación corresponden a 𝜀𝑘 y 𝜂𝑘, las cuales siguen una distribución normal: 𝜀𝑘~ 𝑁(0,𝑊𝑘𝑄𝑘𝑊𝑘𝑇) (3-28) 12 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝜂𝑘~ 𝑁(0,𝑉𝑘𝑄𝑘𝑉𝑘𝑇) (3-29) A partir de este punto, una vez se ha linealizado el sistema, el procedimiento a seguir es similar al aplicado para el Filtro de Kalman discreto. Por tanto, las etapas que componen el Filtro de Kalman extendido son: ▪ Etapa de predicción. Responsable de la proyección del estado y de la covarianza en el instante 𝑡𝑘−1 para obtener el estado a priori para el siguiente instante, 𝑡𝑘. En algunos textos también se denomina “actualización de tiempo”. Las ecuaciones que forman parte de esta etapa son: 𝑥𝑘−=𝑓(𝑥𝑘−1,𝑢𝑘−1,0) (3-30) 𝑃𝑘−=𝐴𝑘−1𝑃𝑘−1𝐴𝑘−1 𝑇+𝑊𝑘−1𝑄𝑘−1𝑊𝑘−1 𝑇 (3-31) ▪ Etapa de actualización. Responsable de la retroalimentación o corrección del Filtro de Kalman. Las ecuaciones que forman parte de esta etapa son: 𝐾𝑘=𝑃𝑘−𝐻𝑘𝑇(𝐻𝑘𝑃𝑘−𝐻𝑘𝑇+𝑉𝑘𝑅𝑘𝑉𝑘𝑇)−1 (3-32) 𝑥𝑘=𝑥𝑘−+𝐾𝑘(𝑧𝑘−ℎ(𝑥𝑘−,0)) (3-33) 𝑃𝑘=(𝐼−𝐾𝑘𝐻𝑘)𝑃𝑘− (3-34) Una vez se completan las dos etapas el proceso se repite tomando la estimación a posteriori para predecir la estimación a priori del siguiente instante. Este aspecto es el que le confiere un carácter más práctico que otros filtros. Al igual que se hizo en el caso del Filtro de Kalman discreto en la Figura 3.2 se muestra un resumen del proceso a seguir por el Filtro de Kalman extendido para cada instante de evaluación. Etapa de predicción 1. Proyección del estado: 𝑥𝑘−=𝑓(𝑥𝑘−1,𝑢𝑘−1)+𝑤𝑘−1 2. Proyección de la covarianza a priori: 𝑃𝑘−=𝐴𝑘−1𝑃𝑘−1𝐴𝑘−1 𝑇+𝑊𝑘−1𝑄𝑘−1𝑊𝑘−1 𝑇 Etapa de actualización 1. Cálculo de la Ganancia de Kalman: 𝐾𝑘=𝑃𝑘−𝐻𝑘𝑇(𝐻𝑘𝑃𝑘−𝐻𝑘𝑇+𝑉𝑘𝑅𝑘𝑉𝑘𝑇)−1 2. Actualización de la estimación: 𝑥𝑘=𝑥𝑘−+𝐾𝑘(𝑧𝑘−ℎ(𝑥𝑘−)) 3. Actualización de la covarianza: 𝑃𝑘=(𝐼−𝐾𝑘𝐻𝑘)𝑃𝑘− Estimación inicial para 𝑥0− y 𝑃0 Figura 3.2. Resumen de las etapas del Filtro de Kalman extendido 13 4 APLICACIÓN DEL FILTRO DE KALMAN A TRAYECTORIAS DE AERONAVES Tras esta introducción al Filtro de Kalman faltaría por conocer el sistema matricial resultante de la adaptación de dicho filtro al caso bajo estudio. Para ello se han tenido en cuenta las siguientes consideraciones: - Se considera el plano horizontal de la trayectoria. - Para una mayor robustez y capacidad de aplicación en diferentes situaciones se va a tener en cuenta sólo la información básica transmitida por un sistema de vigilancia independiente cooperativa. Por tanto, de los 17 parámetros de los que está compuesto cada línea de mensaje ADS-B registrado, como se puede ver en la Tabla 2.1, sólo se tendrá en cuenta la latitud, longitud, el instante en el que se recibe y el valor de “Is On Ground”. Respecto a este último parámetro se tienen 2 posibles situaciones: 1. Aeronave en el aire (IoG = 0): se realiza el estudio general. 2. Aeronave en tierra (IoG = -1): se considera como aeronave perdida. Por tanto, para aplicar el Filtro de Kalman extendido es necesario determinar los siguientes elementos: - Variables de estado. - Variables de control. - Predicción del estado y su covarianza. - Medidas y su covarianza. En este capítulo, en primer lugar, se definen las variables a las que se aplicará el Filtro de Kalman extendido. La determinación de la velocidad con respecto a tierra y del curso a partir de las posiciones se muestra en el apartado 4.2. Por último, en el apartado 4.3, se muestran y se comentan dos casos representativos que permiten ilustrar comportamiento del Filtro de Kalman desarrollado. 4.1 Definiciones 4.1.1 Variables de estado Como ya se ha comentado en varias ocasiones, el objetivo en este de trabajo de aplicar un Filtro de Kalman es obtener la mejor estimación posible de la posición geográfica de una aeronave. Esta posición viene definida por su proyección sobre el plano horizontal. Por tanto, el vector de estado corresponde a: 𝑥𝑘=[𝜙𝑘 𝜆𝑘] (4-1) donde: - 𝜙𝑘: latitud de la aeronave en el instante 𝑡𝑘. - 𝜆𝑘: longitud de la aeronave en el instante 𝑡𝑘. Es necesario destacar que la posición transmitida por ADS-B toma como referencia el modelo terrestre WGS- 84, por tanto, estas variables son referidas a este modelo. 14 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 4.1.2 Variables de control Las variables que intervienen en la estimación de la posición corresponden al curso geográfico, 𝜓, y a la velocidad respecto a tierra, 𝑉𝑔. 𝑢𝑘=[𝜓𝑘 𝑉𝑔𝑘] (4-2) Estas dos magnitudes se obtienen a partir de las posiciones anteriores, por lo tanto, se considera que su error no depende de la posición en el instante de evaluación. Más adelante, en el apartado 4.2 se describirá cómo se obtiene, así como la estimación de los errores asociados. 4.1.3 Predicción del estado y su covarianza Se considera que la aeronave se desplaza siguiendo un movimiento rectilíneo y uniforme. Esta suposición podría entrar en conflicto cuando se analicen trayectorias con virajes, sin embargo, debido a que el tiempo entre mensajes es pequeño (aproximadamente de 1 segundo), se observará en los resultados que es posible aproximar dicho viraje mediante desplazamientos rectilíneos en intervalos de tiempo reducidos. El sistema que determina la evolución de la posición viene dado por [7]: 𝑑𝜙 𝑑𝑡=𝑉𝑔cos(𝜓) 𝑅𝑡+ℎ (4-3) 𝑑𝜆 𝑑𝑡=𝑉𝑔sin(𝜓) (𝑅𝑡+ℎ)cos (𝜙) (4-4) donde: - 𝑅𝑡: valor medio del radio terrestre (6371 m). - ℎ: altitud de la aeronave. En este trabajo se ha supuesto que durante todo el intervalo entre el instante inicial y el de evaluación, definido como Δ𝑡, la velocidad respecto a tierra, 𝑉𝑔, el curso, 𝜓, y la altitud, ℎ se mantienen constantes. Como se mencionó anteriormente, estas variables se obtendrán como función de las posiciones anteriores. En el apartado 4.2 se desarrollará el proceso seguido para obtener estas variables, así como sus características. Es importante aclarar que las ecuaciones (4-3) y (4-4) se basan en un modelo de tierra esférico. Por ello, la diferencia entre este modelo y el WGS-84, en el que se basan las variables de estado, será absorbida por el Filtro de Kalman. Integrando las ecuaciones diferenciales anteriores y discretizándolas se llega a las ecuaciones que permiten estimar la proyección del estado a priori: 𝜙𝑘−=𝑉𝑔𝑘−1cos(𝜓𝑘−1) 𝑅ℎ𝑘−1 (𝑡𝑘−𝑡𝑘−1)+𝜙𝑘−1 (4-5) 𝜆󰆹𝑘−=tan(𝜓𝑘−1)𝑙𝑛|sec(𝜙𝑘)+tan(𝜙𝑘) sec(𝜙𝑘−1)+tan(𝜙𝑘−1)|+𝜆󰆹𝑘−1 (4-6) Para cerrar la etapa de predicción faltaría por conocer cómo obtener la covarianza del estado a priori. Esto es posible aplicando directamente la expresión (3-31), la cual se vuelve a mostrar a continuación: 𝑃𝑘−=𝐴𝑘−1𝑃𝑘−1𝐴𝑘−1 𝑇+𝑊𝑘−1𝑄𝑘−1𝑊𝑘−1 𝑇 (4-7) Como es fácil de imaginar, una vez se han definido las variables de estado, la matriz de covarianza del estado adopta la forma: 15 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝑃𝑘=[(𝜎𝜙𝑘)2(𝜎𝜙𝜆𝑘)2 (𝜎𝜙𝜆𝑘)2(𝜎𝜆𝑘)2], (4-8) donde: - 𝜎𝜙𝑘: desviación típica asociada a la latitud en el instante 𝑡𝑘. - 𝜎𝜆𝑘: desviación típica asociada a la longitud en el instante 𝑡𝑘. - 𝜎𝜙𝜆𝑘: indicador del grado de correlación entre la latitud y la longitud en el instante 𝑡𝑘. En cuanto al error asociado al sistema linealizado se tiene que éste depende de ambas variables de control. Como ya se ha comentado, tanto el curso como la velocidad respecto a tierra se obtienen a partir de las posiciones anteriores, por lo que llevan asociado un error que, a su vez, afecta a la predicción de la posición. Por tanto, considerando que el curso y la velocidad son independientes entre sí se tiene que: 𝑤𝑘=[𝜎𝜓𝑘 𝜎𝑉𝑔𝑘] (4-9) 𝑄𝑘=[𝜎𝜓𝑘 20 0 𝜎𝑉𝑔𝑘 2] (4-10) En el apartado 4.2 se describe con más detalle cómo se obtienen las dos variables de control, así como el proceso seguido para estimar la desviación típica asociada a cada variable, 𝜎𝜓𝑘 y 𝜎𝑉𝑔𝑘. Por último, las matrices 𝐴𝑘−1 y 𝑊𝑘−1 son las responsables de poder estimar cómo aumenta la varianza con el tiempo. Las expresiones concretas que definen cada una de las componentes de ambas matrices se pueden consultar en el Anexo A.1 y A.2. Aun así, en el apartado 4.3 se mostará gráficamente la propagación de la desviación típica asociada a la posición para un par de ejemplos reales. 4.1.4 Medidas y su covarianza Para el caso estudiado la medida de la posición se extrae directamente de los datos almacenados en el archivo “.bst”, es decir, no es necesario ninguna transformación entre la información recibida y el vector de estado. Esto significa que el modelo de medición discretizado toma la forma mostrada a continuación: 𝜙𝑘𝐴𝐷𝑆−𝐵=𝜙𝑘+𝑣1𝑘 (4-11) 𝜆𝑘𝐴𝐷𝑆−𝐵=𝜆𝑘+𝑣2𝑘 (4-12) Siendo 𝑣1 y 𝑣2 el ruido asociado a la latitud y longitud recibida, respectivamente. Cada uno de ellos se asume que es blanco y que su función de densidad sigue una distribución normal. Es decir: 𝑣1𝑘 ~ 𝑁(0,(𝜎𝜙𝑘 𝐴𝐷𝑆−𝐵)2) (4-13) 𝑣2𝑘 ~ 𝑁(0,(𝜎𝜆𝑘 𝐴𝐷𝑆−𝐵)2) (4-14) Por tanto, el ruido asociado a la medida, 𝑣𝑘, está compuesto por la incertidumbre asociada a la posición transmitida por ADS-B. 𝑣𝑘=[𝑣1𝑘 𝑣2𝑘] (4-15) Al estar compuesto por dos normales independientes y asumiendo que los errores en las medidas no están correlacionados se tiene que: 16 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝑉𝑘=[1 0 0 1] (4-16) En este proyecto, a pesar de aplicar un Filtro de Kalman extendido, la medida viene dada por un sistema lineal en el que la medida se obtiene directamente. Esto implica que la matriz 𝐻𝑘 se corresponde con la matriz identidad: 𝐻𝑘=[1 0 0 1] (4-17) Casi todas las variables que definen las etapas del Filtro de Kalman están determinadas, sólo faltaría por conocer cómo obtener la covarianza de la medida, 𝑅𝑘. Para determinar la covarianza haría falta conocer el error asociado a la posición transmitida por ADS-B, 𝑣𝑘. Este error varía según la aeronave, pues depende de los sistemas a bordo que use para determinar la posición, y, a su vez, para una misma aeronave su valor varía con el tiempo. En este sentido la tecnología ADS-B presenta una interesante ventaja, pues incluye en la señal transmitida la calidad de la velocidad y posición transmitidas. Concretamente, la precisión de la posición horizontal viene indicada por un grupo de 4 bits en el mensaje codificado. Como se observa en la Tabla 4.1 este grupo de bits puede tomar 12 valores o categorías diferentes, cada uno de ellos relacionado con un valor concreto de la incertidumbre asociada a la posición horizontal transmitida por ADS-B [1]. Esta incertidumbre define el radio de un círculo centrado en la posición horizontal real de la aeronave de manera que se puede asegurar que la posición transmitida por ADS-B se encuentra dentro de dicho círculo con un 95% de probabilidad. Considerando que sigue una distribución normal es fácil ver que este valor de probabilidad define un intervalo de confianza del doble de la desviación típica. Es decir, la desviación típica asociada a la medida corresponde a la mitad del valor indicado en la Tabla 4.1. Categoría de la precisión (NACp) Incertidumbre estimada (EPU) 0 ≥ 10 NM 1 < 10 NM 2 < 4 NM 3 < 2 NM 4 < 1 NM 5 < 0.5 NM 6 < 0.3 NM 7 < 0.1 NM 8 < 0.05 NM 9 < 30 m 10 < 10 m 11 < 3 m Tabla 4.1. Categoría para la precisión de la posición horizontal Sin embargo, el receptor SBS-3 usado no permite conocer la incertidumbre transmitida por ADS-B ni acceder 17 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B al mensaje codificado. Por lo tanto, ante esta situación se ha elegido como error para la posición transmitida el valor mínimo permitido por EUROCONTROL según el documento “EUROCONTROL Specification for ATM Surveillance System Performance (Volume 2 Appendices)” [5]. Según EUROCONTROL, la precisión de la posición horizontal debe ser igual o superior a la Categoría 8. Esta categoría se aplica a zonas donde la señal de GPS incluye una degradación intencionada normalmente por tratarse de territorios en conflictos. Como este no es el caso de la zona cubierta por el receptor SBS-3 se toma como valor de incertidumbre el definido por la Categoría 9, es decir, 30 m. Este valor para el error de la posición horizontal se va a considerar constante en el tiempo e igual para todas las aeronaves analizadas. Sin embargo, se debe transformar esta desviación típica dada en metros a grados geográficos, siendo esta transformación dependiente del punto en el que se encuentre la aeronave. Por lo tanto, la covarianza de la medida se obtiene como sigue: 𝐶𝑎𝑡𝑒𝑔𝑜𝑟í𝑎 9 →2𝜎𝐴𝐷𝑆−𝐵=30 𝑚 𝜎𝐴𝐷𝑆−𝐵=𝜎𝑥=𝜎𝑦=15 𝑚 → {𝜎𝜙𝑘 𝐴𝐷𝑆−𝐵=𝑝(𝜎𝐴𝐷𝑆−𝐵,𝜙𝑘𝐴𝐷𝑆−𝐵,𝜆𝑘𝐴𝐷𝑆−𝐵) 𝜎𝜆𝑘 𝐴𝐷𝑆−𝐵=𝑞(𝜎𝐴𝐷𝑆−𝐵,𝜙𝑘𝐴𝐷𝑆−𝐵,𝜆𝑘𝐴𝐷𝑆−𝐵) (4-18) 𝑅𝑘=[(𝜎𝜙𝑘 𝐴𝐷𝑆−𝐵)20 0(𝜎𝜆𝑘 𝐴𝐷𝑆−𝐵)2] (4-19) 4.1.5 Resumen del Filtro de Kalman extendido aplicado al cálculo de trayectorias Finalmente, como consecuencia de los resultados anteriores las etapas que componen el Filtro de Kalman extendido se resumen en las ecuaciones mostradas en la Figura 4.1: La estimación inicial necesaria para iniciar el filtro, 𝑥0− y 𝑃0, se obtiene directamente de la medida proporcionada por ADS-B. A partir de este instante tanto la variable de estado como su covarianza vienen dadas por el Filtro de Kalman extendido. Etapa de predicción 1. Proyección del estado: { 𝜙𝑘−=𝑉𝑔𝑘−1cos(𝜓𝑘−1) 𝑅ℎ𝑘−1 (𝑡𝑘−𝑡𝑘−1)+𝜙𝑘−1 𝜆󰆹𝑘−=tan(𝜓𝑘−1)𝑙𝑛|sec(𝜙𝑘)+tan(𝜙𝑘) sec(𝜙𝑘−1)+tan(𝜙𝑘−1)|+𝜆󰆹𝑘−1 2. Proyección de la covarianza a priori: 𝑃𝑘−=𝐴𝑘−1𝑃𝑘−1𝐴𝑘−1 𝑇+𝑊𝑘−1𝑄𝑘−1𝑊𝑘−1 𝑇 Etapa de actualización 1. Cálculo de la Ganancia de Kalman: 𝐾𝑘=𝑃𝑘−(𝑃𝑘−+𝑅𝑘)−1 2. Actualización de la estimación: 𝑥𝑘=𝑥𝑘−+𝐾𝑘(𝑧𝑘−𝑥𝑘−) 3. Actualización de la covarianza: 𝑃𝑘=(𝐼−𝐾𝑘)𝑃𝑘− Estimación inicial para 𝑥0− y 𝑃0 Figura 4.1. Resumen del Filtro de Kalman aplicado 18 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 4.2 Cálculo del curso y de la velocidad respecto a tierra Para cerrar el problema del Filtro de Kalman sólo faltaría conocer el curso, 𝜓, y la velocidad respecto a tierra, 𝑉𝑔, junto con sus desviaciones típicas. Como ya se ha mencionado, ambas variables se obtienen a partir de las posiciones mejoradas de la aeronave bajo estudio. En este apartado se detallará el procedimiento seguido para obtener tanto el valor las variables de control en cada instante de evaluación como la covarianza asociada. 4.2.1 Curso Antes de poder definir las ecuaciones que permiten estimar el curso es necesario realizar algunas consideraciones: ▪ Se asume el modelo de Tierra esférico. ▪ Se considera que las aeronaves navegan según una ruta loxodrómica en la que el rumbo (o curso para el caso de estudio) se mantiene constante. ▪ Se toma como referencia el Norte Geográfico. Por lo tanto, la ecuación que permite estimar el curso teniendo en cuenta estas consideraciones es [6]: 𝜓=atan ( 𝜆𝐵−𝜆𝐴 𝑙𝑛(𝑡𝑎𝑛(𝜋4−𝜙𝐴 2) 𝑡𝑎𝑛(𝜋4−𝜙𝐵 2)) ) (4-20) Donde: - 𝜆𝐵: longitud del punto final. Medida en radianes [rad]. - 𝜆𝐴: longitud del punto inicial. Medida en radianes [rad]. - 𝜙𝐵: latitud del punto final. Medida en radianes [rad]. - 𝜙𝐴: latitud del punto inicial. Medida en radianes [rad]. A la vista de esta expresión se deduce que existe una singularidad cuando la aeronave vuela siguiendo un paralelo (𝜙𝐴=𝜙𝐵). En ese caso, el curso corresponde a: - 𝜓=90º: si 𝜆𝐵>𝜆𝐴 y 𝜆𝐵−𝜆𝐴<180º. - 𝜓=270º: si 𝜆𝐵<𝜆𝐴 y 𝜆𝐴−𝜆𝐵<180º. Para este proyecto el curso se va a obtener a partir de dos posiciones que deben cumplir: - La diferencia temporal entre ellas debe ser aproximadamente igual al intervalo temporal de evaluación. - El punto final del tramo corresponde a la última posición mejorada de la aeronave. Este intervalo temporal de evaluación se establece en 10 segundos. Sin embargo, aunque la tasa de transmisión para la tecnología ADS-B se establece en 1 Hz, en la práctica se observa que esta frecuencia puede variar con el tiempo e incluso se pueden producir pérdidas temporales de la señal. Por tanto, el intervalo temporal variará según el instante en el que se recibe el mensaje más próximo al intervalo de 10 segundos respecto al instante de evaluación. De ahora en adelante, este instante que define el límite inferior del intervalo se define con el índice “j”. De esta manera, 𝑡𝑗 es el mayor tiempo que hace que se cumpla que 𝑡𝑘−1−𝑡𝑗<10 segundos. Además, el valor establecido para el intervalo temporal se puede ver afectado en las siguientes situaciones: 1. Cuando se rastrea una aeronave por primera vez o tras recuperar la señal. En estos casos, al disponer de menos de 10 segundos de información, el intervalo de tiempo es inferior al intervalo temporal de evaluación. Una vez se alcanzan los 10 segundos de diferencia entre el instante de evaluación y el primer mensaje recibido se sigue el proceso normal. 19 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 2. Cuando se detecte un viraje. En primer lugar, para detectar si una aeronave está virando el criterio que se ha tomado es el siguiente: en cada instante de evaluación se comprueba si la posición mejorada y la obtenida a partir de ADS-B difieren más de 3 veces de la desviación típica de la posición ADS-B, 𝜎𝐴𝐷𝑆−𝐵 (correspondiente a un intervalo de confianza del 99.73%). Es decir, si: 𝐷𝑖𝑠𝑡(𝑃𝑜𝑠𝐹𝐾,𝑃𝑜𝑠𝐴𝐷𝑆−𝐵)>3𝜎𝐴𝐷𝑆−𝐵 (4-21) O, lo que es lo mismo: 𝑅𝑡acos(sin(𝜙𝑘)sin(𝜙𝑘𝐴𝐷𝑆−𝐵)+cos(𝜙𝑘)cos(𝜙𝑘𝐴𝐷𝑆−𝐵)cos(𝜆𝑘𝐴𝐷𝑆−𝐵−𝜆󰆹𝑘))>3𝜎𝐴𝐷𝑆−𝐵 (4-22) Cuando se cumpla esta condición el intervalo temporal se reduce en 1 segundo para el siguiente instante de evaluación, hasta un mínimo de 3 segundos. De esta manera se consigue mejorar notablemente el seguimiento del viraje al tomar posiciones más cercanas entre sí, aunque, por el contrario, se obtiene un aumento del error asociado. Cuando deje de cumplirse la condición dada por la expresión 4-21, es decir, una vez finalice el viraje, el intervalo temporal de evaluación aumentará en 1 segundo para cada instante de evaluación hasta alcanzar de nuevo los 10 segundos preestablecidos. Por lo tanto, para estimar el curso en el instante de evaluación, 𝜓𝑘−1, se hace uso de las posiciones mejoradas en 𝑡𝑘−1 y 𝑡𝑗. Es decir, se está considerando que el curso en 𝑡𝑘−1 se mantiene constante en el intervalo [𝑡𝑘−1,𝑡𝑘]. Esta consideración es consecuencia de la suposición tomada según la cual la aeronave sigue un movimiento rectilíneo y uniforme. Teniendo en cuenta estos comentarios y aplicando la notación seguida durante este proyecto se obtiene que la ecuación que define el curso en 𝑡𝑘−1 es: 𝜓𝑘−1=atan ( 𝜆𝑘−1−𝜆𝑗 𝑙𝑛(𝑡𝑎𝑛(𝜋4−𝜙𝑗 2) 𝑡𝑎𝑛(𝜋4−𝜙𝑘−1 2)) ) (4-23) Faltaría por conocer la estimación del error asociado al curso. Una manera de estimar el error asociado a una variable definida por una función conocida se expone en [8], según el cual: si se considera 𝑦1 una función tal que 𝑦1=𝑔1(𝑢1,𝑢2) se puede sustituir por una aproximación lineal centrada en el valor medio (𝑢1,𝑢2) de las variables aleatorias 𝑢1 y 𝑢2. Es decir: 𝑦1≈𝑔1(𝑢1,𝑢2)+∇𝑔1(𝑢1,𝑢2)[𝑢1−𝑢1 𝑢2−𝑢2] (4-24) Donde ∇𝑔1(𝑢1,𝑢2) corresponde al gradiente de la función 𝑔1 evaluado en 𝑢1 y 𝑢2. Por tanto, se tiene que: 𝐸[𝑦1]≈𝑔1(𝑢1,𝑢2) (4-25) 𝜎2[𝑦1]≈∇𝑔1(𝑢1,𝑢2) ∙𝑉𝑔1(𝑢1,𝑢2)∙ ∇𝑔1(𝑢1,𝑢2)𝑇 (4-26) Siendo 𝑉𝑔1(𝑢1,𝑢2) la matriz de covarianza de las variables 𝑢1 y 𝑢2 evaluada en 𝑢1 y 𝑢2. Aplicando este razonamiento a las ecuaciones usadas para calcular el curso se tiene que: 26 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 4.7. Desviación típica de la velocidad respecto a tierra para una trayectoria rectilínea 4.3.2 Trayectoria con viraje Para ilustrar el comportamiento del Filtro de Kalman durante un viraje se ha tomado la trayectoria seguida por la aeronave con distintivo de llamada “34110C” el día 1 de julio de 2013 durante la fase de aproximación. Al igual que se realizó para la trayectoria rectilínea, en la Figura 4.8 se muestra la trayectoria definida tras aplicar el Filtro de Kalman junto con las posiciones predicha y proporcionada por ADS-B para cada instante. 27 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 4.8. Ejemplo de trayectoria con viraje Como ya se mencionó en el apartado 4.2.1 los virajes presentan la peculiaridad de que el intervalo temporal de evaluación se reduce para mejorar la predicción de dicho viraje a partir de segmentos rectilíneos. Para ilustrar gráficamente la mejora en el seguimiento del viraje aplicando esta corrección se propone comparar los resultados obtenidos en la Figura 4.8 con los que se obtendrían si no se corrigiera el intervalo temporal de evaluación entre posiciones y se mantuviera constantemente en 10 segundos, representado en la Figura 4.9. Figura 4.9. Ejemplo de trayectoria con viraje sin aplicar la corrección del intervalo temporal de evaluación A la vista de estos resultados se demuestra una mejora notable en el seguimiento del viraje como resultado de un filtro más ágil. Sin embargo, como consecuencia de reducir el intervalo de evaluación y tomar posiciones más próximas entre sí se obtiene un aumento inevitable de la desviación típica asociada a la posición. 28 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B A continuación, en la Figura 4.10 se muestra la evolución de la desviación asociada a la latitud y longitud cuando. En esta representación se observan diferentes intervalos en los que aumenta la desviación típica. Según sus características estos incrementos se asocian a dos fenómenos diferentes: 1. Durante 5 intervalos (aproximadamente, entre 10 y 50 segundos; 95 y 115 segundos; 150 y 200 segundos; 250 y 270 segundos; 300 y 315 segundos después de iniciarse el análisis de la trayectoria) se pierde la señal, lo que conlleva un aumento del error de predicción y, por tanto, también de la desviación típica asociada a la posición mejorada. Aún así, siempre se obtiene un valor menor que el que se tendría teniendo únicamente en cuenta la posición predicha o la dada por ADS-B. 2. El viraje se inicia a partir de los 380 segundos desde el inicio de la trayectoria. Este hecho lleva asociado un incremento de la desviación típica como consecuencia de reducir el intervalo de evaluación. En la Figura 4.11 se muestra con más detalle el efecto del viraje sobre la desviación típica asociada a la posición. En esta representación se observa un aumento de la desviación típica, especialmente rápido al inicio del viraje. Figura 4.10. Evolución de la desviación típica de la posición para una trayectoria con viraje 29 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 4.11. Representación ampliada del efecto del viraje sobre la desviación típica de la posición Para el caso del viraje es especialmente interesante la comparación entre el curso y el rumbo. Esta comparación se muestra en la Figura 4.12 de manera similar a la seguida para la trayectoria rectilínea. Figura 4.12. Comparación entre el curso y rumbo geográfico para una trayectoria con viraje Al igual que sucedía para el caso anterior, mientras la aeronave se desplaza en una dirección constante tanto el curso como el rumbo presentan un comportamiento similar, siendo la diferencia entre ellas de, 30 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B aproximadamente, 2º. Sin embargo, cuando se produce el viraje se produce un retraso en la evolución del curso respecto al rumbo debido al intervalo temporal de evaluación, como era de esperar. Aproximadamente, el retraso en la evolución del curso es de 8 segundos. En cuanto a la desviación típica asociada al curso, ésta se representa en la Figura 4.13. Observando estos resultados llama la atención la existencia de 4 intervalos en los que la desviación típica aumenta, llegando a valores entre los 4º y 14º de desviación. Estos 4 intervalos se pueden diferenciar según el origen de dicho aumento: - Los tres primeros intervalos (aproximadamente, entre 10 y 50 segundos; 170 y 200 segundos; y 315 y 320 segundos después de iniciarse el análisis de la trayectoria) este incremento de la desviación típica se debe a la pérdida de señal, lo que provoca que al volver a producirse la recepción de datos ADS-B el intervalo de evaluación disminuye. - El cuarto intervalo (aproximadamente, entre 380 a 425 segundos después de iniciarse el análisis de la trayectoria) se debe tanto a la reducción del intervalo de evaluación durante el viraje como a la recuperación de dicho intervalo de evaluación una vez finaliza el viraje. Es decir, durante el viraje se reduce la diferencia de tiempo entre las posiciones tomadas para calcular el curso, pero una vez se detecta que el viraje finaliza se produce el proceso contrario, poco a poco se aumenta la diferencia temporal entre las posiciones hasta alcanzar los 10 segundos. Este proceso de corrección del intervalo de evaluación introduce cierta inestabilidad, efecto que se observa en la Figura 4.13. Figura 4.13. Desviación típica del curso geográfico para una trayectoria con viraje En relación a la velocidad, al tratarse de una maniobra de aproximación, va disminuyendo con el tiempo, como se observa en la Figura 4.14. Esta continua disminución hace que la velocidad respecto a tierra calculada a partir de la posición sea más inestable que para el caso de la trayectoria rectilínea. Además, en relación a lo mencionado para el curso durante el viraje, la mayor diferencia entre la velocidad proporcionada por ADS-B y la obtenida a partir de las posiciones anteriores se produce durante el tramo de viraje y de su posterior corrección, lo que refuerza la importancia de tomar un intervalo de evaluación adecuado. 31 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 4.14. Comparación entre la velocidad respecto a tierra dada por ADS-B y la calculada a partir de la posición para una trayectoria con viraje Observando la evolución del error de la velocidad en la Figura 4.15 se tiene que la desviación típica presenta una evolución similar a la del curso, en el que se distinguían 4 intervalos en los que se producía un aumento de dicha desviación típica. Por lo tanto, como cabe de esperar, el origen de dicho incremento de la desviación típica es el mismo que el descrito para el caso del curso. Figura 4.15. Desviación típica de la velocidad respecto a tierra para una trayectoria con viraje 33 5 DISTANCIA MÍNIMA PREDICHA Y REAL En este punto se dispone de la siguiente información para cada una de las aeronaves rastreadas: - Trayectoria mejorada: conjunto de puntos obtenidos tras aplicar el Filtro de Kalman, así como la desviación típica asociada a cada posición. - Instante de tiempo para cada posición mejorada. Es importante tener en cuenta que la última posición de cada aeronave rastreada puede corresponder a instantes diferentes. - Curso y velocidad respecto a tierra de la aeronave a lo largo de su trayectoria. Por lo tanto, partiendo de esta información, el objetivo de este capítulo es doble: 1. Predecir la distancia mínima entre aquellas aeronaves que se aproximan, así como el instante en el que se produce. 2. Determinar la separación mínima real entre las aeronaves mediante un análisis a posteriori. Esta separación mínima se obtiene a partir de la posición mejorada. Para ello es necesario tener en cuenta las siguientes consideraciones tomadas: - La distancia mínima y el instante en el que se producirá se obtienen por parejas. - Las dos aeronaves que forman la pareja deben estar aproximándose entre ellas. En caso contrario no tiene sentido intentar predecir cuál será el punto de mínima distancia. 5.1 Predicción de la distancia mínima El proceso que permite obtener la separación mínima entre las aeronaves rastreadas consta de los siguientes pasos: 1. Actualizar la posición de la aeronave cuya posición se haya recibido antes. 2. Determinar si las dos aeronaves se aproximan. 3. Propagar las trayectorias. 4. Obtener la distancia mínima. A continuación, se describe cada uno de estos pasos. 5.1.1 Actualización de la posición de la aeronave cuya posición se haya recibido antes Por lo general, la última posición de cada una de las aeronaves de la pareja se tiene en instantes diferentes. Por esta razón, es necesario propagar tanto la posición como la desviación típica de aquella aeronave cuya última posición se haya obtenido antes. Para ello se vuelve a hacer uso de las ecuaciones (4-5) y (4-6). De esta manera, es posible conocer la posición de cada aeronave en el mismo instante de tiempo (se tendrá para una aeronave una posición estimada mientras que para la otra aeronave su posición es predicha) y, por tanto, poder compararlas para determinar si las aeronaves se aproximan o se alejan. 5.1.2 Aproximación entre la pareja Como se había mencionado anteriormente sólo tiene sentido predecir la distancia mínima entre dos aeronaves cuando éstas se estén acercando. Para determinar si una pareja de aeronaves se está aproximando el criterio 34 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B tomado se basa exclusivamente en las posiciones de cada una de ellas, al igual que se viene haciendo a lo largo del proyecto. Como ya se sabe, la posición se trata de una variable estocástica de la cual se conoce su estimación del error. Por tanto, teniendo en cuenta este carácter estocástico, se considerará que dos aeronaves se aproximan o se alejan según términos probabilísticos. Es decir, si se considera que la separación entre dos aeronaves sigue una distribución normal y, además, es independiente de sus valores anteriores se tiene que: 𝐷(𝑡𝑘)~ 𝑁(𝐷𝑘,𝜎𝐷𝑘 2) (5-1) 𝐷(𝑡𝑙) ~ 𝑁(𝐷𝑙,𝜎𝐷𝑙 2) (5-2) Donde 𝐷𝑘 se refiere a la estimación de la distancia entre dos aeronaves en el instante de evaluación, 𝑡𝑘, mientras que 𝐷𝑙 se corresponde a la distancia entre las mismas aeronaves en un instante, 𝑡𝑙 , anterior al instante de evaluación un tiempo Δ𝑡. Debido a que una gran diferencia temporal entre los dos instantes considerados podría dar lugar a una mala interpretación del comportamiento de las aeronaves se considera que la diferencia entre los instantes correspondientes a las dos distancias a comparar será lo más próxima posible a 30 segundos de diferencia, pero siempre que sobrepase los 15 segundos. La separación entre las aeronaves A y B en un instante cualquiera, 𝑘, viene dada por la expresión: 𝐷𝑘=𝑅𝑡acos(sin(𝜙𝐴𝑘)sin(𝜙𝐵𝑘)+cos(𝜙𝐴𝑘)cos(𝜙𝐵𝑘)cos(𝜆𝐵𝑘−𝜆𝐴𝑘)) (5-3) Al ser conocida la expresión que permite estimar la separación entre la pareja de aeronaves es posible conocer su incertidumbre, 𝜎𝐷, aplicando el mismo procedimiento que el ya empleado en el capítulo 3 para la posición, el curso y la velocidad respecto a tierra. De esta manera, se tiene que: 𝑦3=𝐷𝑘≈𝑔3(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘)+∇𝑔3(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘) [ 𝜙𝐴𝑘−𝜙𝐴𝑘 𝜆𝐴𝑘−𝜆󰆹𝐴𝑘 𝜙𝐵𝑘−𝜙𝐵𝑘 𝜆𝐵𝑘−𝜆󰆹𝐵𝑘 ] , (5-4) donde el gradiente se define como: ∇𝑔3= [𝜕𝐷𝑘 𝜕𝜙𝐴𝑘𝜕𝐷𝑘 𝜕𝜆𝐴𝑘𝜕𝐷𝑘 𝜕𝜙𝐵𝑘𝜕𝐷𝑘 𝜕𝜆𝐵𝑘], (5-5) cuyas expresiones de las derivadas parciales se pueden consultar en Apéndice A. Ya que se supone que los errores entre las aeronaves son independientes, y también entre las propias coordenadas la varianza corresponde a: 𝑉𝑔3= [ 𝜎𝜙𝐴𝑘 20 0 0 0 𝜎𝜆𝐴𝑘 20 0 0 0 𝜎𝜙𝐵𝑘 20 0 0 0 𝜎𝜆𝐵𝑘 2 ] . (5-6) Por tanto, la desviación típica quedaría definida como: 𝜎𝐷𝑘=√∇𝑔3(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘)𝑉𝑔3(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘)∇𝑔3(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘)𝑇 (5-7) 35 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝜎𝐷𝑘 2=(𝜕𝐷𝑘 𝜕𝜙𝐴𝑘)2𝜎𝜙𝐴𝑘 2+(𝜕𝐷𝑘 𝜕𝜆𝐴𝑘)2𝜎𝜆𝐴𝑘 2+(𝜕𝐷𝑘 𝜕𝜙𝐵𝑘)2𝜎𝜙𝐵𝑘 2+(𝜕𝐷𝑘 𝜕𝜆𝐵𝑘)2𝜎𝜆𝐵𝑘 2. (5-8) Una vez se conoce la distribución probabilística de la separación entre dos aeronaves para cualquier instante evaluado queda definir el criterio que permite determinar si esta pareja de aeronaves se aproxima o se aleja. Lógicamente, desde un punto de vista determinista dos aeronaves se aproximan cuando la separación entre ellas disminuye con el tiempo. Trasladando esta idea a un proceso estocástico, como el seguido en este proyecto, se considera que dos aeronaves se aproximan cuando la probabilidad de que la separación en el instante de evaluación sea menor que la separación en un instante anterior es, como mínimo, del 90%. Este criterio se puede resumir en: 𝑃(𝐷𝑗>𝐷𝑘)≥0.9. (5-9) Considerando que 𝐷𝑗 y 𝐷𝑘 son independientes entre sí y teniendo en cuenta las propiedades de la distribución normal para la diferencia entre variables aleatorias normales independientes se sabe que: 𝑃(𝐷𝑗>𝐷𝑘) ~ 𝑁(𝐷𝑘−𝐷𝑗, 𝜎𝐷𝑘 2+𝜎𝐷𝑗 2) (5-10) 5.1.3 Propagación de la posición Cuando, siguiendo el criterio anterior, se determina que la pareja de aeronaves se aproxima, se procede a la propagación de la posición de cada una de ellas a partir del instante de evaluación y hasta un máximo de 20 minutos después. Para ello se vuelve a hacer uso de las ecuaciones (4-5) y (4-6). Teniendo en cuenta que la predicción de la trayectoria futura se compone de predicciones discretas (trayectoria definida por posiciones en instantes concretos) es importante determinar el tiempo de paso de esta propagación, es decir, la diferencia temporal entre cada par de instantes consecutivos en los que se predice la posición. Por ejemplo, si dos aeronaves se desplazan a 240 m/s y en sentidos opuestos, al coger un tiempo de paso de predicción de 1 segundo se estaría cometiendo un error en la medida de 480 m. A su vez, un tiempo de paso pequeño supondría un elevado número de predicciones, las cuales podrían ralentizar demasiado la predicción. Por esta razón, como parámetro de diseño se ha decidido tomar 0.1 segundos como tiempo de paso, ya que conlleva un error asumible y el programa es capaz de procesarlo con fluidez. Cuando se realiza la propagación de la posición, cuanto más se aleje la posición predicha mayor será la desviación típica asociada a dicha posición. Para poder conocer cómo se propaga la incertidumbre asociada a la posición se presenta en la Figura 5.1 la trayectoria predicha, en color rojo, para la aeronave con matrícula “4CA8D5” una vez se estabiliza el comportamiento del Filtro de Kalman. 43 6 RESULTADOS Los ejemplos mostrados anteriormente son una muestra de la diversidad de los casos analizados. Por esta razón se propone analizar esta pluralidad de casos de forma global, siendo el objetivo de este análisis conocer la precisión de la predicción según su antelación respecto al instante de distancia mínima real. Por ello, tras aplicar el Filtro de Kalman a cada uno de los mensajes recibidos y realizar posteriormente las predicciones de las separaciones entre las diferentes aeronaves se llega al análisis a posteriori. Para este estudio, a partir del instante en el que se realiza la predicción y aquel en el que se produce realmente la separación mínima, se obtiene la antelación con la que se realiza la predicción. De esta forma, las predicciones se ordenan según su antelación y se agrupan cada 5 segundos. Es decir, el primer grupo reúne aquellas predicciones realizadas entre 0 y 5 segundos antes de alcanzar el instante de mínima separación, el segundo grupo reúne aquellas predicciones con una antelación entre 5 y 10 segundos, y así sucesivamente. Para cuantificar cómo de buena es cada predicción se ha decidido tomar como parámetros de evaluación los dos siguientes: - Valor medio del error de predicción, 𝜇𝑒𝑝. Para cada una de las predicciones que componen cada grupo se calcula el error de predicción, es decir, la diferencia entre la distancia predicha y la real. A partir de esta información es fácil obtener el valor medio de la diferencia entre distancias cada 5 segundos. - Desviación típica del error de predicción, 𝜎𝑒𝑝. El procedimiento es similar al seguido para el valor medio, salvo que en este caso interesa la dispersión de las medidas que conforman cada grupo o, en otras palabras, la desviación existente para cada grupo entre el error de cada una de las predicciones y su valor medio. Con esta definición de los parámetros de evaluación se estaría estimando la precisión de la predicción sólo prestando atención al error de predicción absoluto, sin tener en cuenta el valor de la distancia mínima. Es decir, se estaría dando el mismo peso o relevancia si se obtiene un error de predicción de 50 m cuando la separación mínima entre las aeronaves es de 10 km que para otro caso en el que el error es también de 50 m, pero la separación mínima es de 100 km, cuando realmente la predicción es mucho más precisa para el segundo caso. Debido a ello se decide normalizar el error de predicción respecto a la distancia mínima real, quedando definido cada parámetro de evaluación como se muestra a continuación: 𝜇𝑒𝑝=1𝑁∑𝐷𝑖min 𝑟𝑒𝑎𝑙−𝐷𝑖𝑚𝑖𝑛 𝑝𝑟𝑒𝑑 𝐷𝑖min 𝑟𝑒𝑎𝑙 𝑁 𝑖=1 (6-1) 𝜎𝑒𝑝=√1𝑁∑(𝐷𝑖min 𝑟𝑒𝑎𝑙−𝐷𝑖𝑚𝑖𝑛 𝑝𝑟𝑒𝑑 𝐷𝑖min 𝑟𝑒𝑎𝑙 −𝜇𝑒𝑝)2 𝑁 𝑖=1 , (6-2) donde 𝑁 corresponde al número de predicciones que componen el grupo en cuestión. Con estas ideas aclaradas se está en disposición de pasar al plano final de los resultados. Para este análisis global se ha aplicado el proceso desarrollado en este proyecto a los mensajes recibidos entre las 16:54:45 y 23:54:31 del día 15 de mayo de 2013. En resumen, durante estas 7 horas se han recibido y, posteriormente, analizado 372387 mensajes. Tras someter cada uno de estos mensajes al proceso descrito el análisis global arroja los resultados mostrados en las Figuras 6.1 y 6.2. En el eje vertical izquierdo de la Figura 6.1 se muestra la evolución del error medio de la predicción según la antelación con la que se realiza, mientras que en la Figura 6.2 se representa la dispersión del error de predicción. A su vez, en el eje vertical derecho se indica el número de predicciones que componen cada uno de los grupos de 5 segundos en los que se ha dividido el eje horizontal en el que se refleja la antelación de la predicción. Como no puede ser de otra forma el eje vertical derecho es el mismo para ambas figuras. 44 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Para ambas representaciones se observa cómo a medida que aumenta la antelación respecto al instante de mínima distancia el número de predicciones decrece. Teniendo en cuenta que la predicción se inicia cuando una pareja de aeronaves comienza a aproximarse, al aumentar la antelación es menos probable que dicho fenómeno se haya producido, ya sea porque aún no se ha empezado a rastrear alguna de las aeronaves de la pareja o porque en instantes posteriores se aproximarán debido a que una de ellas realizará un viraje. Por ejemplo, hasta 12 minutos antes del punto de mínima separación la pareja de aeronaves puede estar alejándose o, incluso, que aún no se haya recibido un mensaje de alguna de ellas, por lo que sólo se tendrán predicciones a partir de dicho instante. En la Figura 6.1 se observa que a partir de una antelación de 600 segundos el número de predicciones es inferior a 150. Figura 6.1. Error de predicción según la antelación Figura 6.2. Desviación estándar del error de predicción según la antelación Observando ambas figuras y atendiendo tanto a los resultados obtenidos como al número de predicciones según su antelación se pueden diferenciar dos regiones. La primera abarca las predicciones realizadas con una antelación inferior a 600 segundos, mientras que la segunda se extiende desde los 600 segundos de antelación hasta el límite de 1200 segundos. 45 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B En esta segunda región llama especialmente la atención el tramo entre 12 y 15 minutos de antelación, en el que se observa un error normalizado medio superior a 10, es decir, el error entre la distancia mínima predicha y la real es mayor a 10 veces el valor de esta última. A su vez, la dispersión del error normalizado este tramo también alcanza valores muy elevados, siendo generalmente superior a 10. En relación a la desviación típica para algunos valores de antelación comprendidos en este tramo se obtienen desviaciones muy bajas, pero esto es debido a que el número de predicciones con dicha antelación es inferior a 5. A primera vista, estos resultados son muy llamativos, especialmente para la única predicción que se realiza con 890 segundos de antelación respecto al punto de mínima distancia, en el que se alcanza un error normalizado de -72.29. El valor negativo indica que la distancia mínima predicha es mayor que la distancia mínima real. Aunque en esta segunda región el número de predicciones es reducido y no supone una muestra representativa del comportamiento del modelo aplicado, a continuación se muestra una explicación de por qué se obtienen estos valores. Todas las predicciones realizadas en este intervalo (entre 12 y 15 minutos de antelación) con un error normalizado superior a 5 pertenecen a un único caso, correspondiente a la pareja de aeronaves 4CA7B4 – 4CA8EA. En la Figura 6.3 se representa la trayectoria realizada por cada aeronave, en la que la aeronave 4CA7B4 (en color azul) se comienza a rastrear a, aproximadamente, 40º de latitud norte y 3.5º longitud oeste. A partir de este instante sigue una trayectoria en dirección suroeste y, tras recuperar la señal, se desplaza en dirección noreste. Es durante este tramo cuando se recibe el primer mensaje de la aeronave 4CA8EA (en color negro) a, aproximadamente, 40.3º de latitud norte y 4.3º de longitud oeste, cuya trayectoria se inicia en dirección suroeste para luego virar en dirección norte-noreste, describiendo lo que se podría asemejar a un lazo. Figura 6.3. Representación de las trayectorias para la pareja 4CA7B4 – 4CA8EA Una vez aclarada las trayectorias de la pareja en la Figura 6.4 se representan aquellas predicciones cuyo error normalizado es superior a 20 (en color verde) junto con el punto de mínima separación real (en color rojo). En esta figura se observa la gran diferencia entre la posición en la que se produce la mínima separación y las predichas para la aeronave 4CA8EA como consecuencia del primer viraje realizado por esta aeronave. A esta gran diferencia se le suma el hecho de que la posición en la que se produce la mínima separación real es casi coincidente (en el plano horizontal), lo que conlleva un valor de distancia muy pequeño. Esto hace que al calcular el error entre predicción y realidad y normarlizarlo respecto a esta última se obtengan valores tan elevados. 46 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 6.4. Representación de los puntos de mínima separación con error normalizado alto Para más detalle en la Figura 6.5 se muestran la evolución de la separación real entre la pareja 4CA7B4 – 4CA8EA respecto al tiempo transcurrido desde que se empieza a rastrear a la primera aeronave de la pareja junto con el valor de la separación mínima real y los predichos. En esta representación se tiene que la distancia mínima predicha está entre 40 y 100 km, mientras que la separación mínima real obtenida es de 993.8 m. Además, otro hecho a destacar es que, aunque el punto predicho es claramente diferente al punto de separación mínima real, el instante en el que se alcanzan unos y otros es similar. Figura 6.5. Evolución de la separación para la pareja 4CA7B4 – 4CA8EA Centrando ahora la atención en la primera región, es decir, en aquellas predicciones realizadas con una antelación inferior a 10 minutos. En este tramo se dispone de un número de predicciones suficientemente significativo como para sacar conclusiones más determinantes respecto al comportamiento del modelo aplicado. Para una mejor visualización de los resultados obtenidos en esta primera región se muestra en las Figuras 6.6 y 6.7 una vista ampliada de la evolución del error normalizado medio y su desviación típica, respectivamente. 47 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 6.6. Error de predicción para antelación menor a 600 segundos Figura 6.7. Desviación estándar del error de predicción para antelación menor a 600 segundos En esta región destacan los resultados obtenidos cuando la antelación de la predicción está comprendida entre 220 segundos y 325 segundos. En este tramo se dispara tanto el error normalizado medio como la desviación típica. Esto se debe a que en este intervalo se producen numerosos casos similares al obtenido para la pareja 4CA7B4 – 4CA8EA mostrado anteriormente. Estos casos se caracterizan por: 1. Realización de un viraje por parte de una de las aeronaves que compone la pareja, lo que provoca que el punto de mínima separación real se distancie de los puntos de mínima separación predichos antes de realizar el viraje. 2. La separación mínima obtenida en la realidad suele ser pequeña, del orden de 1 km. A modo de ejemplo, en la Figura 6.8 se muestra uno de los casos que componen este tipo de predicciones. En esta figura se representan las trayectorias para la pareja 020098 - 020073. La aeronave 020073 se desplaza en dirección norte hasta que realiza un viraje en dirección noroeste. Por otro lado, la aeronave 020098 describe la misma trayectoria, pero en sentido inverso. Es por ello que el punto de mínima separación prácticamente coincide en el plano horizontal. A su vez, se representan las posiciones en las que se predice que se alcanzará la 48 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B mínima separación y cuyo error normalizado es superior a 10. Figura 6.8. Trayectorias de la pareja de aeronaves 020098 - 020073 En esta representación se demuestra la primera de las características comentadas anteriormente, es decir, la presencia de un viraje. Para la demostración de la segunda característica de estos casos se añade la Figura 6.9, en la que se representa la evolución de la distancia para la pareja 020098 - 020073 desde que se recibe el primer mensaje de una de las aeronaves de la pareja. El valor de las distancias mínimas predichas está comprendido entre los 2 km y los 30 km, mientras que el valor obtenido como separación mínima en la realidad es de 46.51 m. Además, la mayoría de los puntos en los que se predice que se alcanzará la distancia mínima no distan más de 15 segundos desde el instante en el que posteriormente se producirá el punto de mínima separación real. Figura 6.9. Evolución de la separación para la pareja 020098 – 020073 El número de predicciones cuyo error normalizado es superior a 5 para antelaciones comprendidas entre 220 segundos y 325 segundos es de 75, es decir, el 0.26% de las predicciones consideradas en este intervalo. Sin embargo, su influencia es tan grande que si no se tuvieran en cuenta la evolución del valor medio del error normalizado que se obtendría sería la mostrada en la Figura 6.10. 49 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 6.10. Error de predicción para antelación menor a 600 segundos tras corrección Como era de esperar, al eliminar aquellos casos cuyo error normalizado es mayor que 5 para las predicciones realizadas con una antelación entre 220 segundos y 325 segundos los resultados obtenidos en este intervalo mejoran en cuanto a precisión de la predicción y mantienen la tendencia mostrada en esta primera región. Considerando el error normalizado en valor absoluto se obtiene una tendencia ascendente cuanto mayor es la antelación con la que se realiza la predicción. Este es el comportamiento esperado, pues cuanto mayor es la antelación de la predicción existe una mayor probabilidad de que se produzca algún fenómeno (el más común suele ser un viraje) que provoque un desvío de la trayectoria respecto a la predicha. Realizando un análisis cuantitativo del valor medio del error normalizado se obtiene que, por lo general, éste no sobrepasa el límite del 10% hasta que las predicciones se realizan con una antelación superior a 330 segundos. Especialmente llamativo son los resultados obtenidos cuando la antelación es inferior a 150 segundos, en los que el error normalizado medio es inferior, en valor absoluto, al 2%. Es necesario destacar que el valor medio del error normalizado toma mayoritariamente (salvo para antelaciones pequeñas) valores negativos. Atendiendo a la expresión 6-1 el valor negativo se traduce en una sobreestimación de la separación mínima predicha, es decir, la separación mínima que se alcanza en la realidad es inferior a las separaciones predichas. Para analizar la evolución de la desviación típica en esta primera región en la Figura 6.11 se muestra la evolución de la región delimitada por el doble de la desviación típica del error normalizado respecto a su valor medio. Esta área, o región de confianza, representa para cada una de las antelaciones consideradas el intervalo dentro del cual se estima que estará el error normalizado con una probabilidad del 95%. En la Figura 6.11 se observa que cuando la antelación es inferior a 40 segundos la región de confianza se mantiene por debajo del 20%. Sin embargo, a medida que aumenta la antelación de la predicción la disparidad de casos aumenta, lo que conlleva un aumento de la desviación típica. Este hecho se traduce en un ensanchamiento de la región de confianza a medida que crece la antelación, salvo cuando la antelación es superior a 520 segundos (en este tramo se tiene una menor dispersión en los errores de predicción como consecuencia de una reducción en el número de predicciones analizadas). 50 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B Figura 6.11. Región de confianza del 95% para el error de predicción normalizado Como ya se ha comentado anteriormente, al considerar que las aeronaves se desplazan siguiendos trayectorias rectilíneas y uniformes la realización de un viraje conlleva que la predicción se desvíe de la realidad. Este fenómeno es la principal causa del incremento del error de predicción y de la desviación típica conforme aumenta la antelación de la predicción. Del conjunto de predicciones realizadas con una antelación inferior a 600 segundos, se tiene un total de 12881 de predicciones afectadas por virajes, lo que supone el 7.87% de las predicciones en este intervalo. Sin embargo, a pesar de ser una cantidad relativamente baja su influencia tienen una gran relevancia en los resultados obtenidos y son los responsables principales de que la región de confianza se ensanche al aumentar la antelación. 51 7 CONCLUSIONES Y LÍNEAS FUTURAS Este proyecto nace con el objetivo de predecir la separación mínima entre aeronaves a partir de la posición recibida mediante ADS-B. Para ello, cada vez que se lee uno de los mensajes recibidos se le aplica el proceso siguiente, el cual se puede dividir en dos partes: 1. En primer lugar, se ha desarrollado un Filtro de Kalman extendido en el que el sistema ha sido modelado aplicando las ecuaciones (4-3) y (4-4), mientras que la observación necesaria para actualizar el estado del filtro se ha tomado directamente de la posición proporcionada por ADS-B. A su vez, para estimar los errores asociados a la predicción y a la medida se han seguido dos procedimientos diferentes: - Predicción: la estimación del error asociado al vector de estado se ha obtenido a partir de la covarianza asociada al sistema linealizado. - Observación: en este caso, para la desviación típica de las observaciones se ha considerado el límite establecido por EUROCONTROL para la posición transmitida por ADS-B. Además, como parte del Filtro de Kalman es necesario conocer el curso y la velocidad de las aeronaves. Estas magnitudes se han conseguido estimar a partir de las posiciones mejoradas anteriores. De esta manera, se consigue una mejora de la posición de las aeronaves respecto a la obtenida mediante ADS-B. 2. Posteriormente, se determina aquellas aeronaves que se están aproximando y, tras propagar sus trayectorias, se estima la distancia mínima de separación junto con el instante en el que se alcanzará. Finalmente, para poder evaluar la precisión del modelo de predicción se ha realizado un análisis a posteriori. Este análisis se basa en la distancia mínima real, la cual se obtiene una vez se han leído todos los mensajes y se dispone de las trayectorias finales. Por tanto, se ha tomado como parámetro de evaluación el error entre la distancia mínima predicha y la real normalizado respecto a la distancia mínima real. Tras analizar los datos recogidos durante 7 horas y determinar el error normalizado del conjunto de predicciones realizadas se extraen dos resultados principales: por un lado, el valor medio del error normalizado según la antelación con la que se realiza la predicción; por otro lado, la desviación del error normalizado respecto al valor medio según la antelación con la que se realiza la predicción. A la vista de la evolución de ambos resultados se observa, como era de esperar, un aumento de ambos parámetros cuanto mayor es la antelación con la que se realiza la predicción, tanto porque la precisión de la predicción decae como por una mayor presencia de virajes. Sin embargo, existe un par de tramos en los que se obtienen resultados realmente llamativos por su elevado valor tanto para el error normalizado como para la desviación típica. Estos valores son consecuencia de dos fenómenos. Por un lado, la realización de un viraje posterior a las predicciones cuando el horizonte temporal es elevado (para horizontes temporales pequeños se obtiene un buen comportamiento). Por otro lado, una separación mínima en la realidad reducida, lo que conlleva que las predicciones afectadas lleven asociadas un error normalizado muy elevado. Aún así, es necesario tener en cuenta que el modelo desarrollado se basa únicamente en la posición transmitida por la aeronave, lo que lo hace más robusto e independiente. Por tanto, a partir de una información muy básica es capaz de obtener una buena estimación de las trayectorias futuras cuando las aeronaves se desplazan con rumbo constante y con aceleraciones suaves. En este sentido, a lo largo del proyecto se han tomado diferentes suposiciones o simplificaciones que han permitido adaptar el alcance de este trabajo a los medios disponibles. Este hecho implica la existencia de diferentes mejoras o desarrollos futuro que permitan ampliar el campo de investigación iniciado a partir de este proyecto. A continuación, se proponen algunas de las líneas de mejora que se podrían seguir: 1. Retomar el planteamiento inicial para el proyecto, es decir, obtener información meteorológica a partir del registro MRAR del Mode-S EHS. Aunque al inicio del proyecto no se disponían de medios para 58 Predicción de distancia mínima entre aeronaves a partir de datos ADS-B 𝜕𝐷𝑘 𝜕𝜆𝐵𝑘(𝜙𝐴𝑘,𝜆󰆹𝐴𝑘,𝜙𝐵𝑘,𝜆󰆹𝐵𝑘)=𝑅ℎ(cos𝜙𝐴𝑘cos𝜙𝐵𝑘sin(𝜆󰆹𝐵𝑘−𝜆󰆹𝐴𝑘)) √1−(sin𝜙𝐴𝑘sin𝜙𝐵𝑘+cos𝜙𝐴𝑘cos𝜙𝐵𝑘cos(𝜆󰆹𝐵𝑘−𝜆󰆹𝐴𝑘))2 (A-20)