Vista previa del material en texto
UNIVERSIDAD POLITÉCNICA DE MADRID
ESCUELA TÉCNICA SUPERIOR DE INGENIEROS
DE TELECOMUNICACIÓN
LOCALIZACIÓN Y SEGUIMIENTO DE
AERONAVES MEDIANTE SISTEMAS DE
MULTILATERACIÓN DE ÁREA EXTENSA
TESIS DOCTORAL
Jorge José Abbud Momma
Ingeniero de Telecomunicación
Madrid, 2015
Dpto. de Señales, Sistemas y Radiocomunicaciones
Grupo de Procesado de Datos y Simulación
E.T.S.I. DE TELECOMUNICACIÓN
UNIVERSIDAD POLITÉCNICA DE MADRID
LOCALIZACIÓN Y SEGUIMIENTO DE
AERONAVES MEDIANTE SISTEMAS DE
MULTILATERACIÓN DE ÁREA EXTENSA
TESIS DOCTORAL
Autor: Jorge José Abbud Momma
Ingeniero de Telecomunicación
Director: Gonzalo de Miguel Vela
Doctor Ingeniero de Telecomunicación
Madrid, 2015
TRIBUNAL
Tribunal nombrado por el Magnífico y Excelentísimo Sr. Rector de la
Universidad Politécnica de Madrid, el día .
Presidente:
Vocal:
Vocal:
Vocal:
Secretario:
Suplente:
Suplente:
Realizado el acto de defensa y lectura de la Tesis el día en la Escuela Técnica
Superior de Ingenieros de Telecomunicación.
CALIFICACIÓN:
EL PRESIDENTE LOS VOCALES
EL SECRETARIO
i
A Carolina, Teruyo y José Luis.
ii
iii
AGRADECIMIENTOS
A mi director, Gonzalo de Miguel, por su apoyo y colaboración, sabiendo que
compagino esta Tesis con mi actividad profesional.
A él y a Jesús Besada por informarme de la existencia del ESAV2011. Allí
presenté mi primer artículo en un congreso.
A Hiromi Miyazaki por informarme de la existencia del congreso EIWAC2013. A
Shigeru Ozeki y el equipo de ENRI por la organización del mismo.
A mi padre, que fue el primero en mencionar el Doctorado. A mi madre, que
desde la paciencia ha estado apoyándome en este camino.
A mi familia, tanto los que están como los que estuvieron, que siempre ha
tenido la formación como un pilar básico para el crecimiento y el desarrollo del
ser humano.
A mis Sensei José Manuel, Juan y Juan Antonio, por enseñarme la disciplina,
la constancia y el esfuerzo desde otra perspectiva, a través del Karate.
A Carolina, que con su cariño y su amor me ha apoyado desde el primer día en
este camino, sobre todo en los momentos más difíciles, dejando que dedique
casi todo el tiempo posible a esta Tesis. Sin ella, seguramente no hubiese ni
empezado. Gracias por aguantarme.
iv
v
ÍNDICE
ÍNDICE DE TABLAS .............................................................................................................. xii
RESUMEN ............................................................................................................................ xiv
SUMMARY ............................................................................................................................ xv
CAPÍTULO 1. INTRODUCCIÓN ..................................................................................... 1
1.1. Necesidades presentes y futuras del espacio aéreo ................................................................... 1
1.1.1. Crecimiento del transporte de pasajeros y de carga ........................................................... 1
1.1.2. Usuarios civiles y militares .................................................................................................. 3
1.1.3. Los aviones no tripulados ................................................................................................... 3
1.2. La gestión de un bien escaso: el espacio aéreo .......................................................................... 4
1.2.1. La fragmentación conduce a la ineficiencia ........................................................................ 4
1.2.2. El objetivo: un único espacio aéreo europeo ...................................................................... 7
1.3. El concepto CNS/ATM ............................................................................................................... 13
1.3.1. Descripción de los sistemas ............................................................................................. 14
1.3.2. La convergencia de los sistemas ...................................................................................... 26
CAPÍTULO 2. FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM ............................ 28
2.1. Descripción del sistema WAM ................................................................................................... 28
2.1.1. Definición .......................................................................................................................... 28
2.1.2. Historia de los sistemas de localización hiperbólica ......................................................... 28
2.1.3. Principio de funcionamiento de la multilateración ............................................................. 30
2.1.4. Arquitecturas existentes ................................................................................................... 31
2.2. Planteamiento del caso ideal ..................................................................................................... 39
2.2.1. Definiciones ...................................................................................................................... 39
2.2.2. Caso ideal ......................................................................................................................... 40
2.3. Errores que afectan a la determinación de posición .................................................................. 42
2.3.1. Efecto de los errores aleatorios ........................................................................................ 42
2.3.2. Efecto de los errores sistemáticos .................................................................................... 53
2.3.3. Concepto de pseudodistancia........................................................................................... 60
2.4. Planteamiento del problema matemático .................................................................................. 62
vi
CAPÍTULO 3. SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA
DE ERRORES DE PROPAGACIÓN Y SINCRONISMO ..................................................... 63
3.1. Necesidades del sistema ........................................................................................................... 63
3.2. Análisis del sistema de ecuaciones para el sistema de localización hiperbólico ....................... 64
3.3. Solución al problema muestra a muestra .................................................................................. 64
3.3.1. Mediante calibración ......................................................................................................... 64
3.3.2. Mediante tráfico de oportunidad ....................................................................................... 65
3.4. Solución al problema de seguimiento ........................................................................................ 69
3.4.1. Técnicas de seguimiento de blancos ................................................................................ 69
3.4.2. Seguimiento de los errores sistemáticos .......................................................................... 76
3.4.3. Arquitectura del sistema ................................................................................................... 79
3.4.4. Implementación ................................................................................................................ 80
CAPÍTULO 4. EXPERIMENTOS .................................................................................. 86
4.1. Diseño de los experimentos ...................................................................................................... 86
4.2. Ajustes experimentales .............................................................................................................88
4.2.1. Configuración de las estaciones base .............................................................................. 88
4.2.2. Análisis del intervalo de confianza de la media ................................................................ 90
4.2.3. Configuración del filtro de seguimiento ............................................................................. 91
4.2.4. Configuraciones de las aeronaves .................................................................................... 92
4.3. Resultados experimentales ....................................................................................................... 95
4.3.1. Escenario A ...................................................................................................................... 95
4.3.2. Escenario B ...................................................................................................................... 98
4.4. Discusión ................................................................................................................................. 102
4.4.1. Acerca de las hipótesis ................................................................................................... 102
4.4.2. Acerca del cumplimiento de los requisitos ...................................................................... 103
4.4.3. Otras consideraciones .................................................................................................... 104
CAPÍTULO 5. CONCLUSIONES Y TRABAJOS FUTUROS ...................................... 105
5.1. Conclusiones ........................................................................................................................... 105
5.2. Trabajos futuros ...................................................................................................................... 107
5.2.1. Mejoras del sistema de seguimiento ............................................................................... 107
5.2.2. Siguiente fase de los experimentos ................................................................................ 107
5.2.3. Otras líneas de investigación .......................................................................................... 107
vii
CAPÍTULO 6. BIBLIOGRAFÍA .................................................................................... 108
ANEXO I. SISTEMAS DE COORDENADAS Y TRANSFORMACIONES ................... 116
I.1. Sistema geodésico mundial - WGS ......................................................................................... 116
I.1.1. Elipsoide WGS84 ............................................................................................................ 116
I.1.2. Sistemas de coordenadas .............................................................................................. 117
I.1.3. Relación entre el nivel del mar y el elipsoide WGS84 .................................................... 120
I.2. Transformaciones de coordenadas ......................................................................................... 121
I.2.1. Llh – ECEF ..................................................................................................................... 121
I.2.2. ECEF – ENU................................................................................................................... 123
ANEXO II. MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN
HIPERBÓLICA 125
II.1. Método iterativo ....................................................................................................................... 126
II.1.1. Procedimiento ................................................................................................................. 126
II.1.2. Ejemplos ......................................................................................................................... 127
II.2. Métodos algebraicos ............................................................................................................... 128
II.2.1. Algoritmo de Bancroft ..................................................................................................... 128
II.2.2. Algoritmo de Chan, Ho ................................................................................................... 131
ANEXO III. UNSCENTED KALMAN FILTER ................................................................. 136
III.1. Introducción ............................................................................................................................. 136
III.2. Planteamiento del problema .................................................................................................... 136
III.3. Unscented Transform .............................................................................................................. 137
III.4. Unscented Kalman Filter ......................................................................................................... 139
III.5. Ejemplo ................................................................................................................................... 141
III.5.1. Problema ........................................................................................................................ 141
III.5.2. Resolución mediante EKF: desarrollo ............................................................................. 142
III.5.3. Resolución mediante UKF .............................................................................................. 142
III.5.4. Experimento.................................................................................................................... 143
ANEXO IV. COTA DE CRAMÉR-RAO PARA LOCALIZACIÓN HIPERBÓLICA ....... 145
IV.1. Localización esférica ............................................................................................................... 145
IV.2. Localización hiperbólica .......................................................................................................... 146
viii
ÍNDICE DE FIGURAS
Figura 1.1. Tráfico anual de pasajeros a nivel mundial (en billones de PKT). ............................. 2
Figura 1.2. Efecto de las zonas militares en las rutas aéreas. .................................................... 3
Figura 1.3. Trayectoria de un vuelo regular Madrid – Múnich. .................................................... 5
Figura 1.4. Sistemas constituyentes del la red europea de gestión del tráfico aéreo. ................ 8
Figura 1.5. Esquema con los problemas existentes en la gestión del tráfico aéreo europeo. .. 11
Figura 1.6. Esquema (ideal) del espacio aéreo europeo tras la implementación de SESAR. .. 11
Figura 1.7. Mapa con los nueve bloques funcionales de espacio aéreo (FAB). ........................ 12
Figura 1.8. Precisión (2D) de los sistemas de navegación. ....................................................... 14
Figura 1.9. El ASR-23SS (Raytheon) está compuesto por un PSR (abajo, con sus iluminadores
y su reflector en rejilla) y un SSR (arriba). .................................................................................. 16
Figura 1.10. Estación SSR. ........................................................................................................ 17
Figura 1.11. Principio de funcionamiento de la multilateración. ................................................. 22
Figura 1.12. Papel del sistema de WAM como constituyente del sistema de vigilancia en el
esquema del control de tráfico aéreo. ......................................................................................... 23
Figura 1.13. Esquema de intercambio de datos de vigilancia. .................................................. 24
Figura 2.1. Geometría de multilateración en 2D con tres estaciones receptoras. ..................... 30
Figura 2.2. Obtención de TDOA a partir de las medidas de TOA [13].......................................32
Figura 2.3. Obtención de TDOA mediante correlación cruzada de las señales recibidas [13]. 33
Figura 2.4. Topologías existentes para la sincronización de los receptores [13, 18]. ............... 34
Figura 2.5. Arquitectura de reloj común [13, 18]. ....................................................................... 35
Figura 2.6. Esquema de arquitectura distribuida [13, 18]. ......................................................... 36
Figura 2.7. Arquitectura de sincronismo “stand-alone” mediante GNSS [13, 18]. ..................... 37
Figura 2.8. Arquitectura de reloj acoplado con GPS [20]. .......................................................... 38
Figura 2.9. Arquitectura de sincronismo mediante GNSS "common-view" [13, 18]. ................. 38
Figura 2.10. Ramas de una familia de hipérbolas definidas por múltiplos enteros de σi. .......... 44
Figura 2.11. Desviación típica del error horizontal de posicionamiento para un ruido de medida
de TDOA de 10ns (unidad: m). ................................................................................................... 44
Figura 2.12. Desviación típica del error horizontal de posicionamiento para un ruido de medida
de TDOA de 33ns (unidad: m). ................................................................................................... 45
Figura 2.13. Influencia de la geometría de la constelación en la localización. .......................... 45
ix
Figura 2.14. Mapa de HDOP para la constelación C. ................................................................ 47
Figura 2.15. Mapa de VDOP para la constelación C. ................................................................ 47
Figura 2.16. Mapa de HDOP para la constelación H. ................................................................ 47
Figura 2.17. Mapa de VDOP para la constelación H. ................................................................ 48
Figura 2.18. Cota inferior del error cuadrático medio – componente horizontal (m). ................ 51
Figura 2.19. Cota inferior del error cuadrático medio – componente vertical (m). ..................... 51
Figura 2.20. Método de Bancroft: cociente de la desviación típica de la componente horizontal
de error por Hmin. ......................................................................................................................... 52
Figura 2.21. Método de Bancroft: cociente de la desviación típica de la componente vertical de
error por Vmin. ............................................................................................................................... 52
Figura 2.22. Método de Chan-Ho: cociente de la desviación típica de la componente horizontal
de error por Hmin. ......................................................................................................................... 53
Figura 2.23. Método de Chan-Ho: cociente de la desviación típica de la componente vertical
de error por Vmin........................................................................................................................... 53
Figura 2.24. Temperatura y presión de la atmósfera en función de la altitud. ........................... 54
Figura 2.25. Valor del coíndice de refracción (N) en función de la altitud. ................................ 55
Figura 2.26. Ley de Snell entre dos medios con diferente índice de refracción. ....................... 56
Figura 2.27. Diferencia entre la trayectoria recorrida por una onda radioeléctrica y la distancia
entre transmisor y receptor, en función de la altitud del transmisor. .......................................... 58
Figura 2.28. Cociente de la división de ∆P para una aeronave a una altitud determinada entre
∆P para una aeronave a una altitud de 14000m, para una misma distancia euclídea............... 58
Figura 2.29. Cociente de la división de ∆P para una aeronave a una altitud determinada entre
∆P para una aeronave a una altitud de 14000m, para una misma distancia euclídea............... 59
Figura 3.1. Función densidad de probabilidad de la aceleración en el modelo de Singer. ....... 73
Figura 3.2. Concepto de propagación de la media y la varianza: Monte Carlo, EKF y UT [78]. 76
Figura 3.3. Arquitectura del sistema. ......................................................................................... 79
Figura 3.4. Diagrama de bloques de la etapa de predicción del blanco [78]. ............................ 81
Figura 4.1. Diagrama de bloques del protocolo experimental. .................................................. 87
Figura 4.2. HDOP de la constelación empleada en los experimentos. ..................................... 89
Figura 4.3. VDOP de la constelación empleada en los experimentos. ...................................... 90
Figura 4.4. Escenario en base a los datos de la Tabla 4.1 y de la Tabla 4.4. ........................... 93
Figura 4.5. Escenario en base a los datos de la Tabla 4.1 y de la Tabla 4.5. ........................... 94
Figura 4.6. Sesgo de localización en el escenario A. ................................................................ 96
x
Figura 4.7. Desviación típica del error de localización en el escenario A. ................................. 96
Figura 4.8. Media del error de localización horizontal (unidad: m). ........................................... 97
Figura 4.9. Desviación típica del error de localización horizontal (unidad: m). .......................... 97
Figura 4.10. Valor absoluto de la media del error de localización vertical (unidad: m). ............ 98
Figura 4.11. Desviación típica del error de localización vertical (unidad: m). ............................ 98
Figura 4.12. Sesgo de localización en el escenario B. .............................................................. 99
Figura 4.13. Desviación típica del error de localización en el escenario B. ............................... 99
Figura 4.14. Media del error de localización horizontal (unidad: m). ....................................... 100
Figura 4.15. Desviación típica del error de localización horizontal (unidad: m). ...................... 100
Figura 4.16. Media del error de localización vertical (unidad: m). ........................................... 101
Figura 4.17. Desviación típica del error de localización vertical (unidad: m). .......................... 101
Figura 4.18. Cociente de la desviación típica de la componente horizontal del error entre Hmin.
................................................................................................................................................... 102
Figura 4.19. Cociente de la desviación típica de la componente vertical del error entre Vmin. 102
Figura I.1. Definición del sistema de coordenadas WGS 84. .................................................. 117
Figura I.2. Sistema de coordenadas geodésicas. .................................................................... 118
Figura I.3. El punto P es el centro del sistema de referencia local ENU. ................................ 119
Figura I.4. Ilustración de los distintos conceptos. .................................................................... 120
Figura I.5. Relación entre la altitud ortométrica y la altitud elipsoidal. ..................................... 120
Figura I.6. Altitud del geoide EGM2008 respecto del elipsoide WGS84. ................................ 121
Figura III.1. Principio de la Unscented Transform [69]. ............................................................ 137
Figura III.2. Valor de las medidas en función del tiempo. ........................................................ 143
Figura III.3. Comparación de las soluciones EKF y UKF con el valor real. ............................. 144
Figura III.4. Error en la estimación de posición con EKF y UKF. ............................................. 144
xi
xii
ÍNDICEDE TABLAS
Tabla 1.1. Comparativa del tráfico aéreo en 2010 entre Europa y EEUU. ................................... 6
Tabla 1.2. Cuadro resumen de los sistemas de vigilancia. ......................................................... 22
Tabla 3.1. Relación entre N y el valor mínimo de M. .................................................................. 65
Tabla 4.1. Posición de las estaciones base. ............................................................................... 88
Tabla 4.2. Parámetros característicos de los errores de sincronismo. ....................................... 88
Tabla 4.3. Parámetros del modelo de Singer. ............................................................................. 92
Tabla 4.4. Trayectorias de las aeronaves en el escenario A. ..................................................... 92
Tabla 4.5. Trayectorias de las aeronaves en el escenario B. ..................................................... 93
Tabla II.1. Estimación de la posición en cada iteración – caso convergente. .......................... 128
Tabla II.2. Estimación de la posición en cada iteración – caso divergente .............................. 128
Tabla III.1. Condiciones experimentales. .................................................................................. 143
xiii
xiv
RESUMEN
El continuo crecimiento de la demanda del transporte aéreo, junto con los
nuevos escenarios de intervención militar, están obligando a una optimización
en el uso del espacio aéreo. De este modo, la UE y los EEUU (a través de
SESAR y NextGen respectivamente) han asentado las bases para una nueva
gestión del tráfico aéreo (ATM). Con ello, se pretende aumentar la capacidad
de aeropuertos y rutas aéreas, otorgando mayor flexibilidad al uso del espacio
aéreo sin comprometer la seguridad de los usuarios.
Desde un punto de vista puramente técnico, la clave de este cambio de modelo
está en el conocimiento de la posición de cada aeronave en cada instante. En
este sentido, la tendencia en ATM es el uso de ADS-B como fuente principal de
posicionamiento. Sin embargo, debido a que este sistema está basado en la
difusión de la posición obtenida a través de GPS, es necesario un sistema de
seguimiento independiente. Actualmente, la intención es migrar del radar
secundario de vigilancia (SSR) a la multilateración de área extensa (WAM), con
el fin de mejorar la integridad de la posición para aplicaciones en ruta.
Aprovechando el rápido despliegue de ADS-B, se pretende reutilizar sus
estaciones base para WAM. Cada estación base que recibe el mensaje ADS-B
de la aeronave envía conjuntamente la medida del tiempo de llegada (TOA) de
dicho mensaje al centro de tráfico aéreo. La posición de la aeronave se obtiene
mediante multilateración, cuya técnica consiste en utilizar las medidas de TOA
de un mismo mensaje ADS-B obtenidas en las distintas estaciones base.
El objetivo es estimar la posición de cada aeronave con la mayor precisión
posible.
Para poder diseñar el sistema que permite alcanzar este objetivo, son dos los
aspectos básicos a estudiar. Por una parte, la identificación y posterior
caracterización de los errores (tanto sistemáticos como aleatorios) que afectan
a la medida de TOA. Por otra parte, es necesario el estudio de los sistemas de
seguimiento, basados en versiones sofisticadas del filtro de Kalman (IMM,
UKF).
Una vez establecidos estos dos pilares, la presente tesis doctoral propone un
sistema que permite efectuar el seguimiento de las aeronaves, corrigiendo los
efectos de las principales distorsiones que afectan a la medida de TOA: la
refracción troposférica y el error de sincronismo. La mejora en la precisión de la
localización ha sido evaluada mediante simulación de escenarios hipotéticos.
xv
SUMMARY
The ever-growing demand in the air transportation and the new military
intervention scenarios, are generating a need to optimize the use of the
airspace. This way, the EU and the USA (through SESAR and NextGen
respectively) have set the ground to overhaul the current air traffic management.
The intention is to enhance the capacity of airports and air routes, providing
greater flexibility in the use of airspace without jeopardizing the security of the
end-users.
From a technical perspective, the key for this change lies in the knowledge of
the aircraft position. The trend in Air Traffic Management (ATM) is to rely on
ADS-B as the main source for aircraft positioning. However, this system is
based on the aircraft’s self-declaration of its own (often GPS-based) navigation
solution. It is therefore necessary to have an independent surveillance system.
Nowadays, the intention is to gradually migrate from Secondary Surveillance
Radar (SSR) towards Wide Area Multilateration (WAM) in order to enhance
surveillance integrity for en-route applications.
Given the fast deployment of ADS-B, the aim is to use its base stations for
WAM. Each station sends the Time of Arrival (TOA) of the received ADS-B
messages to the air traffic center (ATC). The aircraft position is obtained
through multilateration, using the TOA of the same message measured by each
station.
The aim is to accurately estimate the position of each aircraft.
Knowledge from two key areas has to be gathered prior to designing such a
system. It is necessary to identify and then characterize the errors (both
systematic and random) affecting the TOA measurements. The second element
is the study of tracking systems based on sophisticated versions of the Kalman
filtering (e.g. IMM, UKF).
Based on this knowledge, the main contribution of this Ph.D. is an aircraft
tracking system that corrects the effects of the main errors involved in the TOA
measurement: tropospheric refraction and synchronization issues. Performance
gains in positioning accuracy have been assessed by simulating hypothetical
WAM scenarios.
xvi
1
CAPÍTULO 1. INTRODUCCIÓN
La multilateración de área extensa (WAM) es una técnica de vigilancia
independiente y cooperativa basada en la diferencia en los tiempos de llegada
de la señal transmitida por la aeronave a una serie de estaciones base. Esta
técnica, que está empezando a ser empleada en distintos puntos del planeta,
está llamada a sustituir los sistemas de radar secundario. Este capítulo
presenta la multilateración de área extensa desde la perspectiva de la gestión
del tráfico aéreo.
1.1. Necesidades presentes y futuras del espacio aéreo
Las necesidades relacionadas con la gestión del espacio aéreo vienen
determinadas por dos fenómenos: el continuo crecimiento del tráfico aéreo civil
y las necesidades de los distintos usuarios.
1.1.1. Crecimiento del transporte de pasajeros y de carga
En las dos últimas décadas, se ha observado el desarrollo económico en
distintos países de nuestro planeta, por ejemplo:
• los países del sudeste asiático, una vez que Japón y Corea del Sur han
llegado a cierta madurez económica y demográfica.
CAPÍTULO 1
2
• los BRIC (Brasil, Rusia, India y China), países con un fuerte crecimiento
del PIB. El caso de China es notable, dado que su PIB es ya el segundo
a nivel mundial.
• los MINT (México, Indonesia, Nigeria y Turquía), que teniendo menor
peso en PIB que los mencionados anteriormente, están experimentando
un mayor crecimiento.
• ciertos países exportadores de petróleo que buscan diversificar su
industria, véanse Emiratos Árabes Unidos o Qatar.
• los países de Europa oriental, que han adaptado sus industrias tras la
caída de la URSS, e intentan tener cierto peso dentro de la Unión
Europea.
Durante las últimas cuatro décadas, la demanda del transporte aéreo ha
desafiado a los distintos escenarios adversos. Prueba de ello es la evolución
del tráfico aéreo, que se ha visto multiplicado por diez1 (Figura 1.1) entre 1970
y 2010, a pesar de los contratiempos en forma de crisis.
En este sentido, laOACI prevé un crecimiento medio anual del 4,6% para el
tráfico aéreo de pasajeros, y un 6,6% para el tráfico aéreo de carga para el
horizonte de 2025.
Figura 1.1. Tráfico anual de pasajeros a nivel mundial (en billones de PKT) 2.
1 PKT (Pasajero-kilómetros transportados), en inglés RPK (Revenue Passenger-Kilometre): es
el resultado de multiplicar los asientos ocupados por pasajeros de pago en un tramo por la
distancia en kilómetros de dicho tramo. Se trata de una medida comercial de la demanda, que
pondera el número de pasajeros con la distancia recorrida.
2 Fuentes: OACI/ICAO (International Civil Aviation Organization), Airbus
INTRODUCCIÓN
3
1.1.2. Usuarios civiles y militares
Otro aspecto a optimizar es la segregación del espacio aéreo según la
naturaleza del usuario: avión civil o avión militar. Actualmente, las aeronaves
no siguen rutas directas, debido - entre otras razones - a la existencia de zonas
reservadas de espacio aéreo, principalmente para fines militares. Estas
restricciones impiden la creación de trayectorias óptimas, disminuyendo la
capacidad de las rutas aéreas y aumentando los costes (tiempo y combustible).
Figura 1.2. Efecto de las zonas militares en las rutas aéreas. 3
1.1.3. Los aviones no tripulados
Asimismo, una nueva “especie” empieza a participar en la saturación del
espacio aéreo. Los aviones no tripulados han de compartir a su vez el espacio
aéreo. Si bien actualmente su uso es limitado (misiones de ataque, vigilancia
de carreteras, o supervisión de cosechas agrícolas), es de esperar que el
espectro de servicios a los que pueden dar soporte siga ampliándose en las
próximas décadas.
3 Fuente: Dirección General de la Aviación Civil (DGAC)
CAPÍTULO 1
4
1.2. La gestión de un bien escaso: el espacio aéreo
Del apartado anterior se puede concluir que el número de usuarios del espacio
aéreo seguirá creciendo en las próximas décadas. Este apartado empieza
describiendo la situación actual del espacio aéreo europeo, comparándolo con
el estadounidense. Ante la ineficiencia en la gestión de su espacio aéreo, la UE
lanzó una iniciativa para gestionar el tráfico aéreo a nivel europeo mediante
cuatro regulaciones [1-4]. Esta iniciativa se concretó con la creación del
proyecto SESAR[5, 6]: un proyecto tecnológico y operativo con el fin de
modernizar la infraestructura de control de tráfico aéreo.
1.2.1. La fragmentación conduce a la ineficiencia
1.2.1.1. Europa: un espacio aéreo fragmentado y saturado
Europa se encuentra actualmente con un espacio aéreo fragmentado,
implicando un uso ineficiente de la capacidad. Dicha fragmentación es el
resultado de la historia de Europa, en la que el control del tráfico aéreo ha
estado estrechamente asociado a la soberanía y por tanto confinado dentro de
las fronteras nacionales4.
Esta falta de eficiencia puede constatarse en distintos puntos del sistema:
• en la actual red de rutas, optimizadas para las necesidades domésticas
• en los propios centros de control, en su mayoría sobredimensionados
• en los programas sobre la gestión del tránsito aéreo
• en la duplicación de sistemas, debido a la adopción no sincronizada de
los avances tecnológicos y de la adquisición no sistemática
Para poder gestionar el tráfico aéreo de manera segura y eficiente, es
necesario que los controladores de tráfico aéreo dispongan de la información
relativa a las aeronaves:
• que se encuentren dentro del espacio aéreo bajo su responsabilidad.
• que tengan previsto cruzar el espacio aéreo bajo su responsabilidad.
4 Durante el siglo XX, Europa ha sido testigo de las dos Guerras Mundiales (1914-1918 y 1939-
1945) y de la Guerra Fría (1947-1989).
INTRODUCCIÓN
5
A modo de ejemplo (Figura 1.3), una aeronave que cubre el trayecto entre
Madrid y Múnich ha de cruzar el espacio aéreo de cuatro países (España,
Francia, Suiza y Alemania). Es decir, a lo largo de su trayecto, esta aeronave
será detectada por distintos sistemas de vigilancia, dependiendo de la posición
y del área de cobertura de los mismos (segmento aeronave - sistema de
seguimiento).
Para que la información suministrada por dichos sistemas sea de utilidad para
la gestión del tráfico aéreo a nivel europeo, los sistemas de seguimiento que
dan soporte a la vigilancia han de ser transparentes al ATC.
Figura 1.3. Trayectoria de un vuelo regular Madrid – Múnich.
1.2.1.2. Comparación con el espacio aéreo estadounidense
Al contrario que su homólogo europeo, el espacio aéreo de los EEUU es único
(debido también a razones históricas). En este sentido, no es sorprendente el
hecho de que el espacio aéreo esté gestionado de manera más eficiente.
La Tabla 1.1 compara el tráfico aéreo entre Europa y EEUU, de la que se
puede destacar lo siguiente:
• la saturación de la capacidad, que afecta a más de 90 aeropuertos en
Europa, frente a únicamente 3 aeropuertos en EEUU.
• en Europa se controla un 40% menos de horas de vuelo, aun teniendo
un 60% más de controladores de tráfico aéreo que en EEUU.
Esta falta de eficiencia en la gestión del espacio aéreo europeo obligó a los
países de la UE a actuar, tal y como se describe en el siguiente apartado.
CAPÍTULO 1
6
Europa EEUU
Extensión (millones de km2) 11,5 10,4
Número de ANSP5 38 1
Controladores (operadores) de tráfico aéreo 16700 14600
Empleados totales 57000 35200
Vuelos controlados - IFR6 (millones) 9,5 15,9
Proporción de vuelos de/a los 34 aeropuertos principales 66% 63%
Proporción de vuelos de Aviación General7 4% 23%
Horas de vuelo controladas (millones) 13,8 23,4
Densidad relativa (horas de vuelo / km2) 1,2 2,2
Distancia media de vuelo dentro del espacio aéreo (NM) 557 493
Centros de ruta 63 20
Aeropuertos con servicios de ATC ≥450 509
denominados “coordinados“8 >90 3
Tabla 1.1. Comparativa del tráfico aéreo en 2010 entre Europa y EEUU9.
5 ANSP (Air Navigation Service Providers): Proveedores de Servicio de Navegación Aérea
6 Es decir, vuelos donde no es necesario tener contacto visual con el terreno.
7 Se trata del conjunto de vuelos civiles, excluyendo los vuelos instrumentales (Instrumental
Flight Rules).
8 Es aquel aeropuerto en el que la demanda de operación por parte de las compañías aéreas
supera la capacidad de las infraestructuras del aeropuerto para satisfacerla, por lo que se hace
necesaria la aplicación de procesos que determinen la prioridad en la autorización de los
vuelos.
9 Fuentes: EUROCONTROL, FAA
INTRODUCCIÓN
7
1.2.2. El objetivo: un único espacio aéreo europeo
1.2.2.1. Un nuevo marco regulatorio: el “Cielo único europeo”
La respuesta a la falta de eficiencia mencionada en el apartado anterior fue la
creación del “Cielo único europeo” (Single European Sky - SES). Se trata de
una iniciativa comunitaria, encuadrada dentro de los tratados de la UE, que
pretende reformar la arquitectura del sistema de gestión de tránsito aéreo
europeo, gestionado en aras de la eficiencia global del sistema, empezando por
el espacio aéreo superior, y de forma coherente con el espacio aéreo inferior.
El marco legal de SES pretende garantizar el cumplimiento de las necesidades
futuras de capacidad y de seguridad, así como con unos requisitos
medioambientales cada vez más exigentes. Dicho marco se erige sobre cuatro
reglamentos básicos [1-4], abarcando la provisión de servicios de navegación
aérea (ANS), la organización y uso del espacio aéreo, y la interoperabilidad de
la red europea de gestión del tráfico aéreo (EATMN).
El reglamento EC 552/2004 [4] (actualizado mediante EC 1070/2009 [7]),
perteneciente a SES, trata de definir los requisitos comunes para garantizar la
interoperabilidad entre los distintos sistemas utilizados para la gestión del
tráfico aéreo, estableciendo así un sistemaarmonizado de certificación de sus
componentes y sistemas.
Este reglamento persigue dos objetivos de manera simultánea:
• por una parte, conseguir la interoperabilidad entre los distintos sistemas,
componentes y procedimientos correspondientes de la red europea de
gestión del tráfico aéreo
• y a su vez, garantizar la introducción de nuevos conceptos de
explotación, acordados y validados, y tecnologías en el ámbito de la
gestión del tráfico aéreo.
Este reglamento también incluye la división de la red europea de gestión del
tráfico aéreo en ocho sistemas (Figura 1.4), enunciados a continuación:
1. sistemas y procedimientos de gestión del espacio aéreo (ASM)
2. sistemas y procedimientos de gestión de la afluencia del tránsito
aéreo (ATFM)
3. sistemas y procedimientos para los servicios del tránsito aéreo (ATS),
en particular sistemas de tratamiento de datos de vuelo (FDPS),
sistemas de tratamiento de datos de vigilancia (SDPS) y sistemas de
interfaz hombre-máquina (HMIS)
4. sistemas y procedimientos de comunicaciones (COM):
CAPÍTULO 1
8
a. para comunicación tierra-tierra (G/G)
b. para comunicación aire-tierra (A/G)
c. para comunicación aire-aire (A/A)
5. sistemas y procedimientos de navegación (NAV)
6. sistemas y procedimientos de vigilancia (SUR)
7. sistemas y procedimientos de servicios de información aeronáutica
(AIS)
8. sistemas y procedimientos para la utilización de información
meteorológica. (MET)
Entre ellos, los sistemas que componen el núcleo de la gestión del tránsito
aéreo son: los sistemas de comunicación (COM), navegación (NAV), y
vigilancia (SUR), que abarcaremos en el apartado 1.3.
Figura 1.4. Sistemas constituyentes del la red europea de gestión del tráfico aéreo.
Por otra parte, cabe destacar que la UE aprobó el reglamento (CE) 2150/2005
del 23 de diciembre de 2005 [8], por el que se establecen las normas comunes
para la utilización flexible del espacio aéreo10, eliminando así una restricción
importante asociada a este problema.
10
"[…] En particular, el presente Reglamento establece normas destinadas a potenciar la
cooperación entre los entes civiles y militares encargados de la gestión del tránsito aéreo que
operan en el espacio aéreo bajo la jurisdicción de los Estados miembros. […]"
INTRODUCCIÓN
9
1.2.2.2. El proyecto tecnológico y operativo: SESAR
Con el fin de concretar la iniciativa del SES, se creó en 2004 el proyecto
SESAR 11 . El objetivo de este programa es la modernización de la
infraestructura de control de tráfico aéreo en Europa, con el fin de satisfacer las
necesidades futuras del espacio aéreo en términos de capacidad y seguridad.
El concepto de SESAR se basa en la divulgación de los mismos datos a todos
los entes involucrados (líneas aéreas, controladores, aeropuertos e incluso
centros meteorológicos) en el proceso de ATM, con el fin de mejorar
drásticamente la planificación y la toma de decisiones para cada vuelo.
De este modo, se espera conseguir una mayor flexibilidad sin provocar retrasos
o distorsiones. La planificación puede hacerse más rápidamente, y más
importante aún, puede ser adaptada a las necesidades en tiempo real. Esto
conlleva a una mayor flexibilidad en la elección de la ruta aérea.
Otro aspecto importante es la divulgación de esta información a los ejércitos,
con el fin de eliminar la segregación del espacio aéreo12.
Así, los cuatro objetivos a largo plazo de SESAR son:
• Triplicar la capacidad del espacio aéreo
• Reducir a la mitad los gastos en gestión de tráfico aéreo
• Mejorar la seguridad en un factor 10
• Reducir el impacto medioambiental de cada vuelo en un 10%
El programa SESAR consta de tres fases:
• Definición (2005-2008): esta etapa concluyó con la entrega del plan
maestro de ATM, donde se define el contenido, al igual que los planes
de desarrollo y despliegue de los futuros sistemas de gestión de tráfico
aéreo. Esta fase fue liderada por EUROCONTROL13, cofinanciada por la
11 Single European Sky ATM Research
12 Cabe destacar que EUROCONTROL ha elaborado la especificación para la aplicación del
uso flexible del espacio aéreo [91] Posteriormente, esta misma entidad ha elaborado una
introducción a la trayectoria de misión [92].
13 EUROCONTROL: nombre abreviado de la Organización Europea para la Seguridad de la
Navegación Aérea (European Organisation for the Safety of Air Navigation). Fundada en
Bruselas en 1960, se trata de una organización civil y militar integrada por 38 estados
miembros que tiene como objetivo el desarrollo de un sistema seguro, eficaz y coordinado del
tráfico aéreo europeo: “Un cielo único europeo”.
CAPÍTULO 1
10
Comisión Europea, y ejecutada por un consorcio formado por las partes
implicadas en el transporte aeronáutico:
o autoridades aeroportuarias
o ANSP (AENA, DSNA, DFS)
o aerolíneas (Iberia, Air France-KLM, Lufthansa)
o constructores aeronáuticos (AIRBUS)
o principales proveedores de equipos (Indra, Selex, Thales).
• Desarrollo (2008-2016): se espera que en esta fase se produzcan los
sistemas tecnológicos de nueva generación establecidos en la fase de
definición. La financiación de esta fase (2100 M€ de presupuesto total)
se ha repartido a partes iguales entre EUROCONTROL, la CE, y la
Industria. Un total de 70 organismos participan en esta fase.
• Despliegue (2014-2020): esta tercera y última fase consiste en la
construcción e instalación a gran escala de la infraestructura definida en
la fase de desarrollo. Por ello, la responsabilidad recaerá íntegramente
sobre la Industria, que financiará los 20000 M€ previstos.
Se espera que SESAR produzca un cambio radical en la gestión del tráfico
aéreo en Europa.
Las dos figuras siguientes ofrecen una comparativa del estado del espacio
aéreo a día de hoy (Figura 1.5), donde se destacan los cuellos de botella
existentes, frente al estado del esquema del espacio aéreo donde las
soluciones han sido implementadas (Figura 1.6).
A modo de ejemplo, aspectos como la congestión - tanto en tierra como en las
TMA – y la espera en pista, son fáciles de percibir para el usuario final
(pasajero).
INTRODUCCIÓN
11
Figura 1.5. Esquema con los problemas existentes en la gestión del tráfico aéreo europeo. 14
Figura 1.6. Esquema (ideal) del espacio aéreo europeo tras la implementación de SESAR.
14 Fuente: SESAR Joint Undertaking
CAPÍTULO 1
12
1.2.2.3. El FAB como paso intermedio
Más recientemente, con el fin de estrechar la cooperación entre los distintos
estados (y por ende, entre sus respectivas autoridades de navegación aérea),
la Comisión Europea definió el concepto de bloque funcional de espacio aéreo
(FAB15), obligando a todos los países miembros de la UE a pertenecer a un
FAB antes del final de 2012. Dicha cooperación afecta a las relaciones entre
los países pertenecientes al mismo FAB en tres niveles: a nivel de Estado
(incluyendo la cooperación civil-militar), a nivel de autoridades nacionales de
supervisión (NSA16) y a nivel de ANSP.
Figura 1.7. Mapa con los nueve bloques funcionales de espacio aéreo (FAB). 17
15 Functional Airspace Block. Este concepto fue definido en el paquete normativo de SES-II (se
denomina así a la revisión de las normativas de SES efectuada en junio de 2008) del siguiente
modo: “un bloque de espacio aéreo basado en exigencias operativas y establecido con
independencia de las fronteras existentes, donde la prestación de servicios de navegación
aérea y las funciones conexas estén basadas en exigencias de rendimiento y optimizadas con
vistas a introducir, en cada bloque funcional de espacio aéreo, una cooperación reforzada entre
proveedores de servicios de navegación aérea o, cuando proceda, un proveedor integrado.”
16 National Supervision Authorities.En la fecha de redacción de la presente Tesis Doctoral,
existen 3 NSAs en España: la AESA (Agencia Estatal de Seguridad Aérea, perteneciente al
Ministerio de Fomento) por la parte civil, el Estado Mayor del Aire por la parte militar, y la
Secretaría de Estado de Cambio Climático representando el ámbito medioambiental.
17 Fuente: EUROCONTROL
INTRODUCCIÓN
13
Cumpliendo con esta legislación, España y Portugal propusieron componer el
bloque funcional de espacio aéreo del sur-oeste (SW FAB) ante la CE.
De este modo, ambos países tienen asignadas las siguientes tareas[9]:
• aumentar la cooperación en los servicios de CNS18/ATM y otros servicios
de soporte
• rediseñar los límites del espacio aéreo en la frontera hispano-portuguesa
• rediseñar el espacio aéreo superior mediante el estudio de un conjunto
de rutas con el fin de aumentar la capacidad, mejorar la eficiencia de los
vuelos y reducir el impacto medioambiental
Acuerdos similares han sido firmados en cada bloque aéreo, con el fin de
optimizar el espacio aéreo confinado en cada uno de ellos (Figura 1.7). Se
trata de una etapa intermedia en el largo camino hacia el espacio único
europeo.
1.3. El concepto CNS/ATM
El concepto CNS/ATM consiste en los sistemas de comunicación, navegación y
vigilancia que emplean tecnología digital, incluyendo sistemas de satélites junto
con varios niveles de automatización, son empleados para dar soporte a un
sistema global de gestión de tráfico aéreo[10].
Los sistemas de comunicación, navegación, y vigilancia son la clave del control
de tráfico aéreo. La mejora en las técnicas (constituyentes), la convergencia de
estos tres sistemas, y la interoperabilidad, permiten aumentar la capacidad del
espacio aéreo, manteniendo (o incluso mejorando) los niveles de seguridad.
La gestión de las aeronaves será finalmente una actividad de colaboración
entre el controlador, el piloto, y el centro de operaciones de la línea aérea. La
aeronave tendrá mayor autonomía y las tareas del controlador de tránsito aéreo
se desarrollarán aprovechando al máximo la automatización.
En este apartado describiremos los sistemas de comunicación, navegación, y
vigilancia, haciendo hincapié en este último y sus constituyentes. Analizaremos
también la convergencia de estos tres sistemas, y enunciaremos los beneficios
esperados en la gestión del tránsito aéreo gracias a las mejoras en los
sistemas CNS.
18 Comunicación, Navegación, y Vigilancia (Communications, Navigation and Surveillance).
CAPÍTULO 1
14
1.3.1. Descripción de los sistemas
1.3.1.1. Comunicación
Los sistemas de comunicación ofrecen el intercambio de datos aeronáuticos y
transmisión de datos entre usuarios y/o sistemas automatizados. Estos
sistemas también se emplean para dar soporte a funciones específicas de
navegación y vigilancia.
Existen diferencias fundamentales entre los sistemas aeronáuticos
convencionales de comunicación y los que forman parte de los nuevos
sistemas CNS/ATM, como por ejemplo:
• la mayoría de las comunicaciones rutinarias se efectúan mediante
intercambio de datos
• las comunicaciones por voz se emplean mayormente en situaciones
excepcionales y de emergencia
• las metas son la conectividad y la operativa a nivel global
De este modo, se espera un aumento en el tráfico de datos y un declive en la
comunicación por voz. Además, en aspectos como la comunicación tierra-tierra,
se está efectuando la migración hacia VoIP, que surgió del sector de las
telecomunicaciones, reflejando la transición de las comunicaciones analógicas
a las digitales.
1.3.1.2. Navegación
El elemento de navegación en los sistemas CNS/ATM está pensado para
determinar la posición de la aeronave de manera precisa y fiable a nivel
mundial mediante la introducción de la navegación basada en satélite (GNSS).
Esta mejora en exactitud y precisión permite (en teoría) - a los usuarios
optimizar las trayectorias, sin necesidad de sobrevolar puntos fijos.
Figura 1.8. Precisión (2D) de los sistemas de navegación.
INTRODUCCIÓN
15
La Figura 1.8 muestra una comparación de la precisión de las distintas
técnicas de ayuda a la navegación (tanto inercial como radionavegación).
La tendencia actual es el uso de sistemas de navegación inercial (ya sea
mediante láser o fibra óptica) con corrección periódica mediante un sistema de
navegación vía satélite (por ejemplo, GPS).
1.3.1.3. Vigilancia
El desarrollo habitual de las operaciones de vuelo en un cierto espacio aéreo
requiere que el ATC (desde tierra) disponga de medios para vigilar que la
separación (horizontal o vertical) entre aeronaves próximas no se reduzca por
debajo de un cierto valor. Unido a este concepto está el de la comunicación
directa rápida entre el ATC y cualquier aeronave que pueda estar en riesgo por
haber reducido su separación mínima con alguna otra aeronave19.
1.3.1.3.1. Técnicas de vigilancia (constituyentes)
Los sistemas de vigilancia se caracterizan mediante dos propiedades: si se
trata de un sistema independiente o no, y si se trata de un sistema colaborativo.
En el caso de los sistemas de vigilancia independientes, la posición de la
aeronave se calcula a partir de las medidas de los sensores en tierra, para ser
transmitida a los centros de control de tráfico aéreo (ATC).
La otra característica es si necesita que la aeronave emita señales para que el
sistema de vigilancia pueda funcionar. En caso afirmativo, se denomina que el
sistema es colaborativo.
Los sistemas de localización son los encargados de localizar las aeronaves con
precisión y exactitud. Son la clave para poder aumentar la capacidad del
espacio aéreo sin poner en riesgo la seguridad. En esta sección se presentan
las fuentes de información de posicionamiento que dan soporte al control de
tráfico aéreo. Estos sistemas son denominados como los constituyentes del
sistema de vigilancia, aplicando la terminología definida en SES.
1.3.1.3.1.1. Radar Primario de vigilancia (PSR)
Utilizados por primera vez durante la Segunda Guerra Mundial, estos radares
emiten señales en el sector de cobertura que tienen asignado. Los ecos
(reflexiones generadas por los blancos) de dichas señales permiten detectar los
blancos que se encuentren dentro del volumen de cobertura, y localizarlos
horizontalmente (latitud y longitud) mediante la combinación de distancia y
azimut.
19 Este asunto ha sido resuelto en Europa y EEUU, mediante el seguimiento obligatorio por
parte de las tripulaciones de los mensajes de resolución de conflictos generados por el ACAS
instalado en cada aeronave.
CAPÍTULO 1
16
Figura 1.9. El ASR-23SS (Raytheon) está compuesto por un PSR (abajo, con sus iluminadores
y su reflector en rejilla) y un SSR (arriba).
Sin embargo, en el ámbito civil, los PSR no son utilizados para obtener la
altitud de los blancos, por dos motivos principales: una resolución en elevación
insuficiente (cuyo efecto en el cálculo de la altura es un error creciente con la
distancia), y el hecho de que el canal vertical empleado en el tráfico aéreo es la
altitud barométrica.
Por otra parte, los PSR no son capaces de identificar los blancos.
Desde la perspectiva del control de tráfico aéreo, al tratarse de un sistema no
colaborativo20, el PSR es utilizado principalmente para la detección de objetos
móviles que puedan afectar al ATC y no lleven transpondedor21 (idea reforzada
tras los atentados del 11 de septiembre de 2001), verificando así la integridad
de la información enviada por las aeronaves equipadas a través de sus
transpondedores (SSR o ADS-B).
20 Es decir, sin que el blanco (la aeronave) emita una señal de contestación a la interrogación al
sistema.
21 No hay que olvidar que los radares detectan y localizan cualquier objeto capaz de reflejar la
señal emitida. Dependiendo del diseñoy de las características atmosféricas, los objetos
candidatos pueden ser por ejemplo: aeronaves, pájaros, helicópteros, globos aerostáticos, etc.
INTRODUCCIÓN
17
1.3.1.3.1.2. Radar secundario de vigilancia (SSR)
La necesidad de identificar aeronaves de un modo sencillo y fiable hizo que el
sistema IFF22 viera la luz. Este sistema, rebautizado SSR23 para el uso civil,
está basado en una estación de tierra y un equipo integrado a bordo de la
aeronave, conocido como transpondedor (equipo que combina un transmisor y
un receptor).
Figura 1.10. Estación SSR.
La estación de tierra transmite (1030MHz) un tren de pulsos denominado
“interrogación” a medida que el haz de su antena está rotando en azimut. En el
momento en que la aeronave reciba esta interrogación, ésta enviará una
respuesta (1090MHz) a la estación de tierra.
El contenido de la información enviada a la estación de tierra depende del
modo de interrogación:
• en el modo A, la aeronave se identifica ante la estación de tierra con un
código de los 4096 posibles (12 bits). Es evidente que de cara a un
centro de control, cada aeronave ha de tener un código único.
Idealmente, una aeronave podía completar un vuelo utilizando el mismo
código a lo largo de la duración del mismo.
Si bien en el momento en que se introdujo esta tecnología dicha
cantidad parecía suficiente, el hecho de que muchos de los códigos
estuviesen reservados para identificar aeronaves en operaciones
22 Identification Friend or Foe: creado durante la 2ª G.M., este sistema se basa en la
identificación positiva de una aeronave, con el fin de evitar el “fuego amigo”. Si una aeronave
con un transpondedor IFF respondía (automáticamente) a las interrogaciones (cierto patrón de
señales) de una estación de control, se consideraba “amigo”. Sin embargo, en caso negativo,
no se puede determinar si la aeronave es “enemiga” o “desconocida”.
23 Siglas del término inglés: ”Secondary Surveillance Radar”
CAPÍTULO 1
18
específicas, junto con el continuo crecimiento del tráfico aéreo
internacional, han hecho que esta cantidad se convierta en un factor
limitante de la capacidad del espacio aéreo, obligando a que la aeronave
cambie de código durante el vuelo.
• el modo B fue utilizado durante un tiempo en Australia de un modo
similar al modo A, pero no se emplea en la actualidad.
• en el modo C, la aeronave transmite la altitud barométrica en que se
encuentra. Únicamente se emplea un rango de 1280 posibles valores
(entre -1200ft y +126700ft) con una granularidad de 100ft24.
Si bien esta distribución era suficiente para monitorizar aeronaves
separadas por más de 1000ft, la creciente congestión en el espacio
aéreo hizo que fuera más importante determinar si un avión dejaba de
volar al nivel de vuelo que tenía asignado. Este hecho creó la necesidad
de que la información de la altitud barométrica tuviera una granularidad
más pequeña (25ft), que se incorporó en el radar secundario modo S,
junto con la interrogación selectiva.
• el modo D nunca ha sido utilizado a nivel operativo.
• el modo S (denominado “selectivo”) fue introducido por la OACI en 1983,
con el fin de paliar las carencias de los modos A y C [11].
Las interrogaciones pueden ser cortas (short squitter, 56 bits) o largas
(extended squitter, 112 bits). Los 24 bits iniciales contienen la dirección
de la aeronave25. Este campo ha sido dimensionado de tal modo que
cada aeronave tenga una única identificación (resolviendo así el
problema originado con el modo A).
Este modo S parte del principio de que cada interrogación está dirigida a
una única aeronave. Sin embargo, para que la estación SSR pueda
interrogar de modo selectivo, primero debe conocer la presencia de la
aeronave. Esta situación se resuelve alternando los períodos de
interrogación general (all-call), con los de interrogación selectiva (roll-
call). En los períodos all-call, las interrogaciones se efectúan con el fin
24 ft (foot): pie, unidad de medida de distancia equivalente a 0,3048m (pie internacional)
25 La dirección de la aeronave viene determinada por la ICAO, y está compuesta por 24 bits (un
total de 16777216 códigos posibles). El valor está definido por la autoridad responsable de
registrar las aeronaves en un país. A cada país le ha sido otorgado un bloque de códigos,
basándose en su tamaño y el volumen del tráfico aéreo. Internamente, cada país puede decidir
la distribución de estos códigos entre aeronaves civiles, militares, y vehículos aeroportuarios
para multilateración. Por ejemplo, España tiene otorgado un bloque de 262144 códigos, de tal
modo que las aeronaves registradas tienen una dirección ICAO que comienza por ‘001101’.
INTRODUCCIÓN
19
de vigilar las aeronaves con transpondedores de modo A y C, así como
para descubrir aeronaves con transpondedores modo S que acaban de
entrar en la zona de cobertura del SSR, con el fin de adquirir su
identidad.
Cabe destacar del modo S que existe una modalidad en la que hasta 16
interrogaciones largas pueden ser encadenadas, con el fin de transmitir
un máximo de 1280bits. Las respuestas de la aeronave pueden ser
encadenadas del mismo modo. El contenido detallado de los mensajes
transmitidos en modo S está especificado en [12].
Así, el modo S es el canal de comunicación más reciente entre aeronave
y estación de tierra. Su aplicación es cada vez mayor, y la información
transmitida desde la aeronave es más detallada (velocidad, cambios de
rumbo, identificación, nivel de servicio, eventos) debido al continuo
aumento del tráfico aéreo.
Desde la perspectiva del control de tráfico aéreo, la información adicional
proporcionada por la propia aeronave mediante SSR complementa a la
obtenida mediante el PSR. Es por ello que durante décadas, la combinación
PSR/SSR ha proporcionado las soluciones de identificación y posicionamiento
al control del tráfico aéreo.
1.3.1.3.1.3. ADS-B
ADS-B26 es un sistema de transmisión de datos en el que la aeronave difunde
su identidad y vector de estado (coordenadas, velocidad, etc.) a las estaciones
de tierra situadas dentro del área de cobertura.
Utilizando la perspectiva de la aeronave, se denomina “ADS-B out” a la
dirección de salida del enlace de datos (posición y datos). La información del
vector de estado viene determinada por el propio sistema de navegación de la
aeronave.
El uso de “ADS-B out”, complementado con los sistemas de localización, otorga
a los ATC la capacidad de obtener en tiempo real las posiciones de las
aeronaves con una mayor precisión (derivadas de GNSS) que si empleara
tecnología radar. Por otra parte, se cataloga como “ADS-B in” a la recepción de
mensajes ADS-B (emitidos por otras aeronaves), permitiendo visualizar el
tráfico aéreo cercano a las aeronaves equipadas con esta tecnología y en el
que se basan los sistemas anticolisión.
A nivel operativo, cabe destacar que el programa CASCADE de
EUROCONTROL está impulsando el uso de ADS-B en Europa, utilizando
Modo S Extended Squitter (ES) como el canal de comunicación.
26 Siglas de Automatic Dependent Surveillance-Broadcast
CAPÍTULO 1
20
1.3.1.3.1.4. Los sistemas de multilateración (MLAT/WAM)
Los sistemas de multilateración (MLAT) están basados en la localización de
emisores mediante el procesado de las señales transmitidas. Algunas técnicas
emplean la diferencia del tiempo de llegada (TDOA27) de la señal emitida en las
diferentes estaciones de tierra situadas en el área de cobertura. Otras técnicas
están basadas en la diferencia de frecuencia de la señal recibida (FDOA28). Lo
más importante de estas técnicas es que son capaces de estimar la posición de
la aeronave independientemente del contenido de la transmisión.
Debido a las mejoras en las prestaciones de los sistemas de multilateración, su
espectro de aplicación se ha visto extendido,abarcando aplicaciones de corta
distancia como la vigilancia aeroportuaria, y de media distancia como el Área
Terminal de Maniobras (TMA29). Esta mejora permite extender el espectro de
aplicaciones al control del tráfico aéreo, a través de la multilateración de área
extensa (WAM30).
Esto es obvio en las zonas continentales. Si a esto se añade el hecho de que
WAM es más asequible (menor inversión inicial, y coste de mantenimiento
inferior) que los sistemas PSR y SSR, es lógico que la introducción de WAM
sea favorecida en detrimento de la renovación de los sistemas SSR
existentes31, como en el caso de la República Checa[13].
WAM es también la elección para la vigilancia aérea en los valles: una
inversión en PSR o SSR sería demasiado elevada y desperdiciada, debido a
las limitaciones en cobertura a baja cota por culpa del terreno. A modo de
ejemplo están Innsbruck (Austria)[14] o Colorado (EEUU)[15].
WAM se está empleando también para la vigilancia aérea costera, donde se
pueden aprovechar las plataformas petrolíferas para instalar receptores WAM;
es el caso de los sistemas desplegados en el Mar del Norte o en el Golfo de
México[16].
Sin embargo, aún no existen ejemplos de sistemas WAM desplegados en
litorales o archipiélagos, con el fin de localizar y vigilar las aeronaves situadas
fuera de la constelación. En este sentido, [17] se encargó de presentar un
27 Time Difference Of Arrival
28 Frequency Difference Of Arrival
29 Terminal Maneuvering Area: se trata del área que rodea a un grupo de aeropuertos. Los
límites de ésta son variables, en función de las necesidades del tráfico, y a su vez pueden estar
divididas en sectores de control. En este espacio todos los vuelos deben estar controlados
30 Wide Area Multilateration
31 Aun siendo una inversión mucho mayor, la adquisición y renovación de los PSR puede
justificarse mediante criterios de seguridad nacional.
INTRODUCCIÓN
21
escenario donde todas las aeronaves se encuentran fuera de la constelación.
Dicho escenario se describe en el apartado 4.2.4.2.
Los sistemas de multilateración pueden ser pasivos o activos. Los sistemas
pasivos emplean las emisiones de las aeronaves que han sido generadas bien
como respuestas a interrogaciones efectuadas por otros equipos, o bien como
emisiones de modo ADS-B (se puede denominar también ADS-B out). Esto
tiene dos ventajas en el uso del espectro: primero, no es necesaria una licencia
de transmisión para el uso de este sistema; segundo, no hay un aumento en el
número de interrogaciones (1030MHz), y por ende, de respuestas (1090MHz).
Según [13], estos sistemas son la mejor solución para ser empleados en los
escenarios con un intercambio importante de datos en la banda de SSR, como
por ejemplo:
• zonas con mucho tráfico y una cantidad importante de aeronaves
equipadas de TCAS32
• zonas con infraestructura existente de vigilancia mediante MSSR
• zonas donde el uso del modo S es obligatorio
Por otra parte, los sistemas activos pueden solicitar respuestas dedicadas a
sus interrogaciones, además de tener la misma funcionalidad que los sistemas
pasivos. El interrogador WAM puede estar constituido por una antena
omnidireccional o a base de antenas sectoriales.
Estos sistemas son muy útiles en las aplicaciones de TMA, o zonas de
aproximación, puesto que se pueden obtener datos con mayor frecuencia,
incrementando la probabilidad de detección y mejorando así la precisión.
Actualmente, los ATC requieren el identificador (obtenido mediante modo A o
modo S) y la altitud barométrica (Modo C).
Adicionalmente, los sistemas WAM activos pueden ser utilizados para obtener
datos específicos de una aeronave. Por ejemplo, estos sistemas pueden
efectuar el seguimiento pasivo mediante TDOA comprobando la dirección
modo S, mientras que las interrogaciones pueden emplearse para obtener el
identificador modo A y la altitud barométrica.
32 Siglas de los términos “Traffic Collision Avoidance System” y “Traffic alert and Collision
Avoidance System”, se trata de los sistemas anticolisión equipados en aeronaves.
CAPÍTULO 1
22
Figura 1.11. Principio de funcionamiento de la multilateración.33
Sistema de
vigilancia
¿Independiente? ¿Colaborativo?
PSR SÍ: datos obtenidos por el radar NO: no es necesario que la
aeronave emita señales
SSR Parcialmente: es independiente
para la localización horizontal,
pero la altitud barométrica
procede de la aeronave
SÍ: se requiere que la aeronave
utilice un transpondedor
SSR Modo S NO : datos proporcionados por la
aeronave
SÍ: se requiere que la aeronave
utilice un transpondedor
ADS-B NO : datos proporcionados por la
aeronave
SÍ: se requiere que la aeronave
tenga la función ADS-B out
MLAT/WAM SÍ : datos obtenidos a partir del
tiempo de llegada de la señal
SÍ: se requiere que la aeronave
emita señales compatibles con SSR
o ADS-B
Tabla 1.2. Cuadro resumen de los sistemas de vigilancia.
El sistema de vigilancia tiende a combinar ADS-B con un sistema
independiente colaborativo (siendo MLAT el más económico), más un radar
primario en las zonas cercanas a los aeropuertos.
33 Fuente: THALES (Global Surveillance Solution Booklet)
INTRODUCCIÓN
23
1.3.1.3.2. Papel del sistema de vigilancia en el ATC
La vigilancia es una función clave en el control del tráfico aéreo. Los sistemas
de vigilancia son los “ojos” de los controladores de tráfico, mostrándoles quién
(qué aeronave), dónde y cuándo está en el cielo. Estos sistemas están en el
origen del proceso del control de tráfico aéreo, detectando las aeronaves y
enviando información detallada al sistema de ATC, permitiendo a los
controladores guiar el avión de manera segura. De esta manera, el control del
tráfico aéreo no es posible sin los sistemas de vigilancia, especialmente en las
áreas con un tráfico intenso.
El control de tránsito aéreo tiene tres aspectos: la planificación, la
monitorización, y el guiado.
La planificación proporciona las hojas de ruta de los distintos usuarios, como
plantilla para comprobar las trayectorias de los blancos detectados
(monitorización), y en caso de resolución de conflictos.
La monitorización consiste en comprobar que el comportamiento de las
aeronaves es acorde a dicha planificación. Así, los sistemas de vigilancia
permiten detectar, identificar y localizar las aeronaves, con el fin de vigilar su
comportamiento.
En caso de conflicto, el ATC necesita resolverlo mediante indicaciones al piloto
para evitar colisiones (guiado y control). En algunos casos, es necesario
considerar la actualización, retroalimentando a la fase de planificación.
La Figura 1.12 muestra cómo encaja el sistema de vigilancia (con WAM como
ejemplo de sistema independiente) en el esquema del control del tránsito aéreo.
Figura 1.12. Papel del sistema de WAM como constituyente del sistema de vigilancia en el
esquema del control de tráfico aéreo34.
34 Fuente: ERA
CAPÍTULO 1
24
1.3.1.3.3. ASTERIX: el intercambio de datos de vigilancia a
nivel europeo
Con el fin de proporcionar la vigilancia de blancos de manera transparente al
usuario, ASTERIX 35 fue introducido en 1986 por EUROCONTROL para
establecer un único formato para el intercambio de datos de vigilancia.
La comunicación entre dos sistemas diferentes - situados incluso en países
diferentes (Figura 1.13) - es posible debido al empleo de un formato común de
datos de vigilancia.
Figura 1.13. Esquema de intercambio de datos de vigilancia.
ASTERIX es un estándar relativo a los niveles de Presentación y Aplicación
(niveles 6 y 7) del modelo OSI36. El medio de transmisión de la información
codificada mediante ASTERIX puede ser cualquier medio de comunicación:
redes de área extensa (WAN), redes de área local (LAN), IP, etc.
La documentaciónde ASTERIX está dividida en Partes, donde cada Parte
agrupa los datos empleados para un propósito específico. Cada Parte contiene
35 All purpose STructured Eurocontrol suRveillance Information eXchange
36 El modelo de interconexión de sistemas abiertos (ISO/IEC 7498-1), también llamado OSI (en
inglés, Open System Interconnection) es el modelo de red descriptivo, que fue creado por la
Organización Internacional para la Estandarización (ISO) en el año 1984. Es un marco de
referencia para la definición de arquitecturas en la interconexión de los sistemas de
comunicaciones.
INTRODUCCIÓN
25
una o más categorías de datos. La información contenida en una categoría está
dedicada a un área específica de aplicación, estableciendo el formato en el que
los datos han de ser transmitidos entre los usuarios de estas aplicaciones.
Existen un máximo de 256 categorías, donde:
• las categorías 000 a 127 se utilizan para aplicaciones civiles y militares
• las categorías 128 a 240 están reservadas para aplicaciones civiles y
militares
• las categorías 241 a 255 se utilizan para aplicaciones extraordinarias
(tanto civiles como militares)
La Figura 1.13 presenta un ejemplo de intercambio de datos de vigilancia,
donde:
• las categorías 019 (mensajes de servicio para configurar la estación) y
020 (mensajes con informes de posición de aeronaves) definen los datos
de salida de multilateración.
• la categoría 021 (mensaje de datos de aeronave) define los datos de
salida de ADS-B.
• la categoría 023 (mensajes de servicio de las estaciones de tierra) define
los mensajes de servicio enviados por una estación de tierra CNS/ATM.
Estos mensajes son emitidos por estaciones de tierra que ofrecen
diferentes tipos de servicio, como por ejemplo: ADS-B, TIS-B37, FIS-B38,
y MLAT.
• las categorías 034 (mensajes de servicio) y 048 (mensajes de
notificación de blancos) definen los formatos de los datos que se
obtienen desde los siguientes equipos: PSR, PSR con procesado MTD39,
SSR y MSSR (SSR modo S)40.
• la categoría 062 (mensajes de pistas) define el formato de los datos de
salida de los sistemas de procesado de datos de vigilancia.
37 Traffic Information Service-Broadcast: servicio ofrecido por una estación de tierra, que
consiste en ofrecer la información de tráfico aéreo a las aeronaves equipadas con ADS-B. Es
especialmente importante cuando no todas las aeronaves están equipadas con ADS-B.
38 Flight Information Services-Broadcast: se trata de la difusión (tierra-aire) de información
meteorológica y aeronáutica
39 Moving Target Detection: detección de blancos móviles
40 SSR Modo S: se trata de un SSR que emplea interrogación selectiva a cada aeronave
CAPÍTULO 1
26
1.3.2. La convergencia de los sistemas
1.3.2.1. CNS como sistema de sistemas
Históricamente – incluso hoy en día, en algunos aspectos – los sistemas de
comunicación (COM), navegación (NAV) y vigilancia (SUR) no han estado
operando como un solo ente, parcialmente por razones de seguridad. La
separación COM-NAV-SUR garantizaba que si el sistema de vigilancia fallaba,
los aviones podían emplear COM y NAV.
Sin embargo, la evolución en los constituyentes del sistema de vigilancia (ADS-
B, WAM) y un control de tráfico aéreo basado cada vez más en los datos
obtenidos a partir de la fusión de sensores (en vez de los datos proporcionados
por un único sensor), están consiguiendo que CNS se establezca como un
sistema de sistemas.
La introducción de los sistemas de navegación vía satélite al mundo civil
(constituyente del sistema de navegación), junto con la introducción de ADS-B
(sistema de comunicación mencionado como constituyente del sistema de
vigilancia) como enlace de datos, muestran un ejemplo de la ruptura de las
fronteras que separan a los tres sistemas, que permite a las aeronaves declarar
su propia posición (en 3 dimensiones) de una manera más precisa, enviando
estos datos de navegación a la función de vigilancia.
1.3.2.2. Expectativas a nivel operacional
Se espera que los sistemas CNS aplicados a la gestión del tránsito aéreo
mejoren el manejo y la transferencia de la información, amplíen la vigilancia, y
mejoren la precisión en la navegación. Esto conduce a la reducción en la
separación entre aeronaves, permitiendo un aumento de capacidad en el
espacio aéreo.
Se espera también que los sistemas en tierra de CNS/ATM intercambien datos
directamente con los sistemas de gestión de vuelo a bordo de la aeronave
mediante un radioenlace de datos, ofreciendo mejoras en la detección y
resolución de conflictos mediante procesado, al igual que pasillos para evitar
conflictos, adaptándose rápidamente a los cambios en los requisitos de tráfico.
De este modo, se espera que el sistema de gestión de tránsito aéreo tenga
más flexibilidad para alojar el perfil de vuelo preferido por cada usuario,
ayudando a los operadores a reducir retrasos y costes operacionales.
1.3.2.3. Impacto en los constituyentes del sistema de vigilancia
Debido a esta convergencia, el papel de los sistemas PSR se verá reducido al
de la detección de aeronaves de manera independiente (tarea que sigue siendo
fundamental), y en menor medida a la localización de las mismas, dado que su
precisión es considerablemente menor a largas distancias.
INTRODUCCIÓN
27
Por otra parte, debido a la introducción del modo S (que proporciona el soporte
físico a ADS-B), los sistemas SSR seguirán dando paso a ADS-B.
La ventaja consiste en que se obtienen las posiciones de las aeronaves en
tiempo real con la precisión que otorga la navegación vía satélite. Sin embargo,
no hay que olvidar que ADS-B no hace más que transmitir la posición
determinada por el propio sistema de navegación de la aeronave (aviónica). Es
por ello necesario comprobar la integridad de esta información en todo
momento.
Dicha comprobación puede conseguirse mediante redundancia, es decir,
empleando diferentes tecnologías de localización mediante satélite (GPS,
GLONASS, Galileo, etc.), o mediante una técnica independiente (WAM).
De este modo, a medida que la gestión del tráfico aéreo se vaya apoyando en
ADS-B, más importante va siendo el papel de las técnicas de multilateración.
Este capítulo ha presentado la multilateración de área extensa (WAM) desde la
perspectiva de la gestión del tráfico aéreo. En el siguiente capítulo se describe
el principio de funcionamiento de WAM, se plantea el caso ideal, y se abordan
los errores que afectan a la determinación de la posición de los blancos.
28
CAPÍTULO 2. FUNCIONAMIENTO BÁSICO
DEL SISTEMA WAM
2.1. Descripción del sistema WAM
2.1.1. Definición
WAM proviene de las siglas del término inglés Wide Area Multilateration
(multilateración de área extensa). Se trata de una técnica colaborativa de
vigilancia basada en la localización hiperbólica, que permite localizar las
aeronaves a partir de la diferencia de los tiempos de llegada (TDOA) de una
misma señal. El tiempo de llegada de la misma señal emitida por la aeronave
depende de la distancia entre la estación de tierra y la aeronave.
A diferencia de los sistemas MLAT, pensados para las aplicaciones
aeroportuarias (tales como la monitorización de aproximación, o la supervisión
de la superficie aeroportuaria), los sistemas WAM están concebidos para la
vigilancia y el seguimiento del tráfico aéreo, tanto en ruta, como en el Área
Terminal de Maniobras.
2.1.2. Historia de los sistemas de localización hiperbólica
Los sistemas de localización hiperbólica fueron introducidos como sistemas de
radionavegación basados en la diferencia entre los tiempos de recepción de
dos señales en una misma estación. Esta diferencia es proporcional (en el caso
ideal) a la diferencia entre las distancias del receptor (aeronave) a las
estaciones transmisoras.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
29El lugar geométrico donde esta diferencia es constante es una hipérbola. El
cruce entre dos hipérbolas ofrece dos soluciones, con lo que es necesario
eliminar esta ambigüedad para obtener la posición. La solución puede ser una
tercera hipérbola, o la correlación con otros sistemas de navegación.
El primer sistema de localización hiperbólica fue “Gee”, introducido por la Royal
Air Force41 en 1941 durante la 2ª Guerra Mundial, tanto para los bombardeos
sobre Alemania, como para la navegación en el Reino Unido, especialmente
para los aterrizajes nocturnos. La diferencia de tiempos se medía a partir de un
osciloscopio situado en la aeronave. A partir de la fundación de la ICAO en
1946, Gee fue considerado como base de la norma mundial para la navegación,
pero la combinación de sistemas VOR y DME fue elegida a cambio: Gee fue
desmantelado en 1970.
A continuación apareció el sistema DECCA (1944), empleado por la Royal
Navy 42 para que los dragaminas pudieran preparar el desembarco de
Normandía. Este sistema estaba basado en la diferencia de fases entre
señales en tiempo continuo, en vez de la diferencia de tiempos. Sin embargo,
era necesario emplear otro sistema complementario para eliminar la
ambigüedad en la posición creada por la ambigüedad en fase.
Por otra parte, el sistema LORAN, empleado por la US Navy para navegación
marítima a larga distancia fue una versión de Gee (diferencia de tiempos),
aprovechando la propagación por onda de superficie.
En las décadas de los años 50 y 60 aparecieron dos sistemas derivados del
DECCA – basados en la diferencia de fase de las señales – que
proporcionaron una solución de navegación sin ambigüedad:
• LORAN-C, gracias a la introducción de los PLL43
• Omega, donde la desambiguación se pudo efectuar gracias a la
introducción de los sistemas de navegación inercial
Por su parte, la URSS desarrolló el sistema CHAYKA (equivalente a LORAN-C)
y el Alpha (una versión de Omega).
Todos estos sistemas están siendo sustituidos por los sistemas de navegación
por satélite (GPS, GLONASS, GALILEO, etc.), también basados en localización
hiperbólica.
41 Ejército del Aire del Reino Unido
42 Marina del Reino Unido
43 Phase-locked loop (lazo de seguimiento de fase): se trata de un sistema de control que
genera una señal cuya fase está relacionada con la fase de la señal de entrada.
CAPÍTULO 2
30
2.1.3. Principio de funcionamiento de la multilateración
El principio de funcionamiento de la multilateración (constituyente del sistema
de vigilancia según la definición de la EC 552/2004[4]) es el mismo que el de
los sistemas de navegación mencionados en el apartado anterior, reforzando el
aspecto de la convergencia de CNS (descrito en el apartado 1.3.2).
Aun siendo la multilateración también un sistema de localización hiperbólica, la
configuración es opuesta a los sistemas descritos en el apartado anterior: en
este caso es la aeronave la que transmite las señales, y las estaciones de tierra
son receptoras.
De este modo, este sistema consiste en un número de antenas que reciben
una misma señal transmitida por la aeronave, y un centro de procesado que
calcula la posición de la aeronave basándose en la diferencia de tiempos de
llegada (TDOA) de dicha señal en cada antena.
Partamos de la base de que la posición del emisor es desconocida. La
obtención de la diferencia del tiempo de llegada de una misma señal en dos
estaciones diferentes permite restringir el espacio a una rama de hiperboloide.
De este modo, en un caso ideal, la intersección de tres ramas de hiperboloide
define de manera unívoca la posición del emisor.
La Figura 2.1 ilustra el principio de funcionamiento para la localización de la
aeronave en dos dimensiones, utilizando tres receptores. La diferencia del
tiempo de llegada entre dos receptores de la misma señal emitida (rojo) permite
obtener una rama de hipérbola (azul), denominada de TDOA constante.
Figura 2.1. Geometría de multilateración en 2D con tres estaciones receptoras44.
44 Fuente: Electronic Navigation Research Institute, Japón
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
31
De este modo, tres estaciones son necesarias para determinar la posición en
dos dimensiones (latitud, longitud). En este caso, dicha posición horizontal se
complementa con la altitud barométrica, que ofrece una medida de altitud que
puede diferir significativamente de la altitud geométrica. Sin embargo, cabe
destacar que la normativa vigente relativa al tráfico aéreo sigue empleando la
altitud barométrica45.
De manera equivalente, al menos cuatro receptores son necesarios para
obtener la posición del emisor en tres dimensiones. Conociendo la posición de
la estación y los tiempos de llegada, se puede estimar la posición de la
aeronave (hasta en tres dimensiones).
2.1.4. Arquitecturas existentes
Los sistemas WAM pueden clasificarse empleando dos criterios:
• en función del método empleado para calcular el TDOA.
• en función del método de sincronización de los receptores.
El apartado 2.1.4.1 describe los métodos de cálculo de TDOA (a partir del
tiempo de llegada, o mediante correlación de las señales recibidas). Por su
parte, el apartado 2.1.4.2 ahonda en los métodos de sincronización de los
receptores.
2.1.4.1. Implementaciones existentes para calcular el TDOA
Según [13], existen principalmente dos métodos para calcular el TDOA: a partir
del TOA, o mediante correlación de las señales recibidas.
2.1.4.1.1. A partir del tiempo de llegada (TOA) de cada señal
El primero (Figura 2.2) consiste en medir los tiempos de llegada de la misma
señal (TOA), y calcular sus diferencias. Este método se emplea cuando se
puede medir con facilidad un flanco determinado, como es el caso de las
señales enviadas por los transpondedores (SSR) de las aeronaves.
El proceso de obtención de TDOA se describe a continuación:
• El bloque demodulador (“Down converter”) recibe la señal de
radiofrecuencia a 1090MHz, y la convierte en una señal de banda base.
La salida puede ser una señal I-Q, o una señal de vídeo.
45 La normativa existente hace referencia a los niveles de vuelo (Flight Level), basados en
altitud de presión, donde cada nivel corresponde a 100ft. Así, FL300 significa 30000ft de altitud
barométrica.
CAPÍTULO 2
32
• El bloque de digitalización (“Digitisation”) emplea un conversor A/D, o
similar, con el fin de digitalizar la señal obtenida a la salida del bloque
anterior.
• Una vez digitalizada la señal, el sistema TOA asociado al receptor
calcula el tiempo de llegada de la señal. El contenido del mensaje es
extraído para poder efectuar la correlación.
• Una vez obtenidas las medidas de TOA para cada receptor, el bloque de
correlación de TDOA agrupa las señales con el mismo contenido, puesto
que estas pertenecen a una misma transmisión de una aeronave
determinada. Así pueden calcularse las diferencias entre los tiempos de
llegada (TDOA).
• El algoritmo de localización (“TDOA algorithm”) calcula la posición de la
aeronave (X, Y, Z) a partir de las medidas de TDOA.
• Finalmente, dichas medidas de posición alimentan un filtro de
seguimiento (“tracker”), encargado de crear una pista para cada
aeronave, mejorando la precisión mediante filtrado y desestimando los
datos erróneos.
Figura 2.2. Obtención de TDOA a partir de las medidas de TOA [13].
2.1.4.1.2. Mediante correlación cruzada
El segundo método (Figura 2.3) consiste en obtener directamente el TDOA
mediante procesado, sin medir el TOA de cada señal. Estos sistemas se
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
33
emplean mayormente en entornos militares (ESM46), o para la localización de
teléfonos móviles durante las llamadas de emergencia.
El proceso de obtención de TDOA se describe a continuación:
• Las etapas de conversión a banda base y digitalizaciónson análogas a
las descritas en 2.1.4.1.1, convirtiendo la señal de radiofrecuencia a
banda base, para luego digitalizarla.
• Tras la digitalización, el bloque de correlación cruzada (‘Cross-
Correlation’) combina las señales digitalizadas para cada pareja de
estaciones. Asumiendo que la misma señal ha sido recibida en ambas
estaciones, puede calcularse la diferencia de tiempos de llegada. La
precisión de este proceso está determinada por el tipo de señal
digitalizada y el multitrayecto, entre otros factores. Cabe destacar que
los algoritmos empleados deben garantizar que la correlación cruzada
de señales no produzca resultados ambiguos o incorrectos,
especialmente cuando las señales recibidas son respuestas SSR que
intrínsecamente no ofrecen buenas propiedades de autocorrelación y
correlación cruzada.
• El algoritmo de localización y el filtro de seguimiento son análogos a los
descritos en la sección 2.1.4.1.1.
Figura 2.3. Obtención de TDOA mediante correlación cruzada de las señales recibidas [13].
46 Electronic Support Measures - Medidas de vigilancia electrónica. Dícese del apartado de
guerra electrónica encargado de detectar, interceptar, localizar, registrar y/o analizar las
fuentes de radiación electromagnética con el fin de localizar una amenaza de manera
inmediata (corto plazo), o bien con el fin de planificar operaciones (largo plazo).
CAPÍTULO 2
34
2.1.4.2. Métodos de sincronización de receptores
El sincronismo es fundamental para los sistemas de multilateración descritos
en el apartado anterior. Con el fin de calcular la posición, es necesario conocer
la diferencia de tiempos de llegada de una determinada señal a dos estaciones
receptoras en el sistema.
La señal recibe la marca de tiempo durante el proceso de digitalización.
Desafortunadamente, esta marca de tiempo no está asociada al instante en
que la señal llega a la antena del receptor. Dicha marca se ve afectada por el
retardo de grupo inherente al proceso de conversión a banda base. Este
retardo ha de ser conocido y considerado a la hora de calcular el TDOA.
Adicionalmente, con el fin de calcular el TDOA correctamente, las marcas de
tiempo han de estar referidas a una base de tiempos común. Es por ello
necesario un mecanismo para sincronizar los relojes locales, o en su defecto,
estimar el error para poder corregirlo a la hora de estimar la posición de cada
aeronave.
2.1.4.2.1. Topologías existentes
El árbol de la Figura 2.4 a continuación clasifica los métodos existentes para la
sincronización de los receptores.
Figura 2.4. Topologías existentes para la sincronización de los receptores [13, 18].
2.1.4.2.1.1. Sistemas de reloj común
Por una parte, los sistemas de reloj común (Figura 2.5) emplean receptores
sencillos, dejando la parte compleja en el centro de procesado. Estos sistemas
reciben las señales de RF de la aeronave, convirtiéndolas a frecuencia
intermedia (FI). La señal de FI es transmitida de cada receptor a la central de
procesado (ubicada normalmente en el centro de la constelación) a través de
un enlace analógico, mediante radioenlace o fibra. La conversión a banda base
y la digitalización se producen en dicha central de procesado empleando un
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
35
único reloj. Con esta arquitectura, no es necesario sincronizar los receptores
situados en cada estación, debido a que la marca de tiempo asociada a cada
señal se genera en un mismo punto del proceso. Sin embargo, el retardo de
grupo existente entre la recepción de la señal en la antena y la central es
significativo, al incluir la conversión a frecuencia intermedia y el enlace
analógico. Según [13], la incertidumbre en la medida del retardo es
proporcional al propio retardo total, dado que éste depende de las condiciones
atmosféricas. Por ello, estos sistemas suelen emplearse en constelaciones
pequeñas.
Figura 2.5. Arquitectura de reloj común [13, 18].
2.1.4.2.1.2. Sistemas de relojes distribuidos
Los sistemas de relojes distribuidos (Figura 2.6) emplean receptores más
complejos con el fin de evitar las incertidumbres derivadas del enlace a
frecuencia intermedia. El receptor se encarga de procesar la señal RF con el fin
de obtener el tiempo de llegada de la señal, tal y como se enuncia en 2.1.4.1.1.
Esta filosofía puede emplearse reutilizando un conjunto de estaciones ADS-B,
dado que los propios receptores ADS-B se encargan de decodificar la señal RF
y establecer la marca de tiempos. Esto es posible gracias a que los datos de
salida de los receptores ADS-B contienen el tiempo de llegada de la señal
ADS-B enviada por la aeronave (ASTERIX CAT-021, Data Items I021/073 e
I021/074)[19]. La región de cobertura resultante sería la intersección de las
regiones de cobertura asociadas a cada estación.
CAPÍTULO 2
36
Figura 2.6. Esquema de arquitectura distribuida [13, 18].
Aun así, es necesario un mecanismo que sincronice los relojes ubicados en
cada estación (ya sea dedicada a WAM, o ADS-B). Existen dos métodos de
sincronismo: mediante transpondedor, o mediante la señal de un sistema de
navegación por satélite (GNSS).
Por una parte, los sistemas sincronizados mediante transpondedor utilizan las
emisiones de un transpondedor (con un reloj de referencia) para calibrar el
retardo existente en cada receptor, y sincronizar los relojes asociados. Esto
implica que el transpondedor de referencia ha de tener línea de vista con cada
uno de los receptores, siendo necesarias las instalaciones de mástiles y torres.
Además, cuanto más importantes sean las distancias entre las estaciones base,
peor es el calibrado, debido a los cambios en el retardo inducidos por la
propagación troposférica.
La alternativa es el sincronismo empleando fuentes externas de tiempo, como
los existentes en los sistemas de navegación por satélite (Figura 2.7). Se trata
de aprovechar los requisitos de precisión asociados a las marcas de tiempo de
las señales GNSS, dado que dichas señales son básicas para los sistemas de
navegación47. Por ejemplo, el sistema GPS permite obtener marcas de tiempo
con un error inferior a 100ns. Puesto que lo importante en los sistemas de
47 Los satélistes GPS son controlados y operados por el Departamento de Defensa de los
EEUU.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
37
multilateración es el TDOA y no el TOA de las señales, es posible sincronizar
los receptores de manera que el error en la medida del TDOA esté por debajo
de 20ns empleando osciladores acoplados con GPS.
Esta técnica emplea un oscilador de cuarzo termoestabilizado o de rubidio. Los
principales factores que afectan a la deriva de dicho oscilador son la
antigüedad y la temperatura. El mecanismo de acoplo se efectúa de manera
similar al bucle de enganche de fase (PLL)48.
En la mayoría de estos osciladores, el bucle de filtrado es sustituido por un
microcontrolador que analiza y compensa los cambios de fase y frecuencia del
oscilador local, los efectos del envejecimiento, temperatura y otros parámetros
ambientales que a la deriva de dicho reloj (Figura 2.8) [20].
Figura 2.7. Arquitectura de sincronismo “stand-alone” mediante GNSS [13, 18].
El sincronismo entre las estaciones puede refinarse obligando a que los
distintos receptores reciban la señal GNSS de un mismo satélite (Figura 2.9).
Esto permite eliminar muchas fuentes de error, obteniendo errores de
sincronismo inferiores a 1ns.
48 Siglas del término inglés, phase-locked loop.
CAPÍTULO 2
38
Figura 2.8. Arquitectura de reloj acoplado con GPS [20].
Figura 2.9. Arquitectura de sincronismo mediante GNSS "common-view" [13, 18].
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
39
2.1.4.2.2. Topologías aplicables
En relación a los mecanismos de sincronización mencionados en el apartado
anterior, hay querecordar que la presente tesis doctoral está enfocada hacia la
localización y el seguimiento de aeronaves en ruta, con lo que el tamaño de la
constelación de receptores se presume importante.
Quedan por ello descartados, el uso de un sistema de reloj común (2.1.4.2.1.1),
y el uso de un transpondedor de referencia en un sistema de relojes
distribuidos.
De este modo, podemos determinar que el sistema bajo estudio requiere un
sincronismo mediante GNSS. Debido a la extensión de la constelación, no se
puede garantizar que los distintos receptores puedan ver simultáneamente al
mismo conjunto de satélites. Por ello, lo más razonable es asumir que la
arquitectura de sincronismo es mediante GNSS “stand-alone”.
Puesto que se espera que este sistema funcione en ausencia de GNSS, es
necesario considerar aspectos como las derivas de los relojes en cada
receptor49, y la manera de mitigar sus efectos.
2.2. Planteamiento del caso ideal
2.2.1. Definiciones
2.2.1.1. Sistema de coordenadas de referencia
El sistema de coordenadas de referencia que se va a emplear a lo largo de
esta tesis es ENU (eje x: Este – eje y: Norte – eje z: Arriba), cuyo centro es un
punto situado en latitud (φ0) y longitud (λ0), en el baricentro de la constelación
de las estaciones base (tomando sus latitudes y longitudes respectivas), y a
nivel del mar.
2.2.1.2. Ubicación de los actores
Es necesario definir las siguientes magnitudes basándose en el sistema de
referencia ENU definido anteriormente:
• Una determinada estación (i-ésima) base está ubicada en (φi,λi,hi),
perfectamente conocida. Para obtener la posición relativa (xi,yi,zi) al
sistema de coordenadas de referencia ENU definido en el apartado
anterior, se remite al lector al ANEXO I.
• Una determinada aeronave (j-ésima) está ubicada en (xj,yj,zj). La
posición de dicha aeronave es a-priori desconocida para el sistema de
vigilancia.
49 Más concretamente, a todas las estaciones menos una, que se empleará como referencia
CAPÍTULO 2
40
2.2.2. Caso ideal
El problema de la localización mediante WAM puede ser considerado como un
simple problema geométrico: se trata de encontrar el punto de concurrencia
entre las ramas hiperboloide de la constelación, donde cada rama es generada
por una TDOA.
El caso ideal asume el cumplimiento de los siguientes supuestos:
• el escenario en cuestión carece de atmósfera (vacío)
• el sincronismo entre relojes es perfecto (sin errores, sin deriva)
• los receptores de las estaciones base son perfectos y el sistema no
introduce ruido (no existe ruido)
2.2.2.1. Sistema de ecuaciones
El primer supuesto enunciado en 2.2.2 permite establecer la relación de
proporcionalidad entre la distancia euclídea (Rij) que separa la estación i de la
aeronave j, y el tiempo de llegada (TOA) de la señal transmitida por la
aeronave j y recibida por la estación i:
ijij R
c
⋅=
0
1
τ
( 2.1 )
En el caso de una única aeronave y un número N de estaciones, la ecuación
anterior pasa a ser el siguiente sistema de ecuaciones:
⋅=
⋅=
⋅=
⋅=
NjNj
jj
jj
jj
R
c
R
c
R
c
R
c
0
3
0
3
2
0
2
1
0
1
1
1
1
1
τ
τ
τ
τ
M
( 2.2 )
Con el fin de obtener una expresión en función de la incógnita (xj,yj,zj), se
emplea la definición de la distancia euclídea entre cada estación y la aeronave,
obteniendo así el siguiente sistema de ecuaciones TOA, donde cada una de
ellas define una esfera:
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
41
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
−+−+−=
−+−+−=
−+−+−=
−+−+−=
222
0
2
3
2
3
2
3
0
3
2
2
2
2
2
2
0
2
2
1
2
1
2
1
0
1
1
1
1
1
NjNjNjNj
jjjj
jjjj
jjjj
zzyyxx
c
zzyyxx
c
zzyyxx
c
zzyyxx
c
τ
τ
τ
τ
M
( 2.3 )
Para obtener el sistema de ecuaciones TDOA, es necesario tomar una de las
medidas de TOA como referencia, obteniendo así el siguiente sistema
compuesto por N-1 ecuaciones, siendo N el número de estaciones. Cada
ecuación del sistema define una rama de hiperboloide:
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
−+−+−−−+−+−=−
−+−+−−−+−+−=−
−+−+−−−+−+−=−
2
1
2
11
2
1
222
0
1
2
1
2
11
2
1
2
3
2
3
2
3
0
13
2
1
2
11
2
1
2
2
2
2
2
2
0
12
1
1
1
zzyyxxzzyyxx
c
zzyyxxzzyyxx
c
zzyyxxzzyyxx
c
jjNjNjNjjNj
jjjjjjj
jjjjjjj
ττ
ττ
ττ
M
( 2.4 )
2.2.2.2. Solución para el sistema de ecuaciones determinado
La solución al sistema de ecuaciones anterior depende del número de
estaciones. Dado que la solución necesaria es tridimensional, el número de
estaciones ha de ser mayor o igual que 4.
Con el fin de aprovechar los recursos existentes, es importante evitar que tres o
más estaciones sean colineales. La colinealidad de tres estaciones supone que
una ecuación sería redundante, al poder definirse como combinación lineal de
las otras.
La solución analítica al caso ideal, definida en [21], se detalla a continuación:
+−
+−
+−
+
−=
−
14
2
1,4
13
2
1,3
12
2
1,2
1
1,4
1,3
1,2
1
1,41,41,4
1,31,31,3
1,21,21,2
2
1
KKr
KKr
KKr
r
r
r
r
zyx
zyx
zyx
z
y
x
j
j
j
( 2.5 )
donde:
−=
−=
−=
11,
11,
11,
zzz
yyy
xxx
ii
ii
ii
( 2.6 )
CAPÍTULO 2
42
y:
( ) ( )
++=
++−++=−=
222
2
1
2
1
2
1
222
11,
iiii
iiiii
zyxK
zyxzyxrrr
( 2.7 )
Si el sistema está sobredeterminado, la solución puede obtenerse empleando
cuatro de las ecuaciones, siendo las restantes empleadas para corroborar la
solución.
Si la solución obtenida empleando 4 de las ecuaciones no pertenece al espacio
de soluciones definido por alguna de las ecuaciones restantes, significa que no
estamos en el caso ideal.
Para resolver dichos sistemas, hay que estudiar los fenómenos físicos que nos
impiden considerar el caso ideal, y entender sus efectos. En ello consiste el
siguiente apartado (2.3).
2.3. Errores que afectan a la determinación de posición
Puesto que la principal aportación de esta tesis está concebida para ser
utilizada en sistemas reales, ninguno de los supuestos enumerados en el
apartado 2.2 se cumple. Es por ello necesario identificar el impacto de los
distintos fenómenos físicos, para poder mitigar su efecto.
Los errores pueden clasificarse en dos categorías: aleatorios, y sistemáticos.
2.3.1. Efecto de los errores aleatorios
2.3.1.1. Fuentes de error
Los errores aleatorios varían de manera aleatoria en función del tiempo. Si bien
no se puede estimar el valor de una muestra en concreto en un instante de
tiempo determinado, es posible caracterizar su comportamiento mediante una
distribución de probabilidad.
Las principales fuentes de errores aleatorios son: el ruido en la señal recibida
(debido al ruido del transmisor, del receptor y del canal), y el ruido de
cuantificación inherente a la medida del TDOA.
Según [13], los errores aleatorios inducen un error en la estimación del tiempo
de llegada que puede llegar a ser de 10ns (equivalente a una pseudodistancia
de 3m), valor que usaremos como referencia en este estudio.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
43
2.3.1.2. Influencia del error
2.3.1.2.1. Influencia según la magnitud del error
En un sistema de localización hiperbólica, la influencia del error aleatorio en la
medida de TDOA se refleja mediante el término adicional n1i, donde los
subíndices 1 e i definen los miembros de la pareja de estaciones base a la que
se asocia el error aleatorio:
( )
( )
( )
+−⋅=−
+−⋅=−
+−⋅=−
NjNjjNj
jjjj
jjjj
nRR
c
nRR
c
nRR
c
11
0
1
1313
0
13
1212
0
12
1
1
1
ττ
ττ
ττ
M
( 2.8 )
Asumiendo una determinada constelación de estaciones, cadaecuación de
este sistema, define una familia de hiperboloides, cuya posición en el espacio
depende del valor del parámetro n1i.
Se asume además, que el error en la medida de TOA (y por ende, de TDOA)
es debido al ruido térmico asociado al receptor50.
Definiendo el ruido aleatorio mediante una v.a. gaussiana de media nula y
desviación típica σi, la probabilidad de que el objeto esté situado en una de las
ramas de la familia de hiperboloides está definida mediante una función
densidad de probabilidad gaussiana.
De este modo, la Figura 2.10 representa la región definida por una de las
ecuaciones51. La presencia del ruido aleatorio nos impide afirmar que el blanco
se encuentra en la rama verde. Lo único que se puede afirmar es que con un
95% de probabilidad, el blanco está situado entre la rama definida por Ni=-2σi y
la rama definida por Ni =2σi. Así, cuanto mayor sea la magnitud del error,
mayor será la región de incertidumbre. Cuanto mayor sea el valor de σ, mayor
será la separación entre las diferentes ramas de hipérbola.
El ruido aleatorio asociado a la medida del TOA induce un error de carácter
aleatorio en la posición. La varianza de dicho error aumenta con la magnitud
del ruido asociado a la medida. La Figura 2.11 muestra la varianza del error de
posición horizontal (proyección del error en el plano xy) para un error de
medida de TOA con una desviación típica de 10ns (equivalente a 3m), mientras
50 Esto implica que el conversor A/D (empleado para obtener el flanco de subida de la señal,
con el fin de medir su tiempo de llegada) tiene una resolución suficiente para poder asumir que
el ruido de discretización es insignificante respecto a la magnitud del ruido término.
51 Con el fin de simplificar la visualización, se ha ejemplificado con un sistema bidimensional.
CAPÍTULO 2
44
que la Figura 2.12 muestra la varianza del error de posición horizontal para un
error de medida de TOA con desviación típica de 33ns (aproximadamente 10m).
Figura 2.10. Ramas de una familia de hipérbolas definidas por múltiplos enteros de σi.
Figura 2.11. Desviación típica del error horizontal de posicionamiento para un ruido de medida
de TDOA de 10ns (unidad: m).
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
45
Figura 2.12. Desviación típica del error horizontal de posicionamiento para un ruido de medida
de TDOA de 33ns (unidad: m).
2.3.1.2.2. Influencia según la geometría de las estaciones
El apartado anterior describió la influencia de la magnitud del error aleatorio en
la localización del blanco, para una determinada geometría de estaciones base.
Este apartado aborda la influencia que tiene una determinada distribución del
error aleatorio, en función de la geometría de la constelación.
La Figura 2.13 muestra la intersección de familias de ramas, es decir, la
solución a dos de las ecuaciones del sistema definido en (2.8). En ambos
casos (A y B), el error aleatorio asociado a cada ecuación tiene la misma
magnitud. Debido a la constelación de las estaciones del caso B, la
incertidumbre en la intersección (definida en verde) es mayor. De este modo,
se puede afirmar que el error de una determinada magnitud es más perjudicial
en el caso B que en el caso A.
Figura 2.13. Influencia de la geometría de la constelación en la localización.
CAPÍTULO 2
46
Para analizar este hecho de manera cuantitativa, nos apoyaremos en el
concepto de GDOP 52 [22, 23], que especifica el efecto multiplicativo de la
geometría de las estaciones en el error en la medida del TDOA. Dicho
concepto fue introducido por los usuarios del sistema de navegación LORAN53,
siendo de especial utilidad para los sistemas de posicionamiento basados en
multilateración (GPS, WAM, etc.). La clave reside en que este valor únicamente
depende de la posición relativa del blanco y las estaciones base.
Por definición, GDOP es el cociente de la división del error de localización entre
el error en la medida de TDOA. Para separar la influencia del error en la
medida en la localización horizontal y la determinación de la altura del blanco,
se han definido dos parámetros auxiliares: HDOP y VDOP.
Con el fin de ilustrar estos conceptos, se han efectuado pruebas con dos
constelaciones formadas por 6 estaciones situadas a nivel del mar:
• la primera (constelación C), con 4 de las estaciones formando un
cuadrado de 100km de lado, y las dos restantes, en el punto medio de
dos de los lados opuestos.
• la segunda (constelación H), en la que las estaciones forman un
hexágono inscrito en un círculo de radio 50km.
La Figura 2.14 y la Figura 2.15 muestran respectivamente el HDOP y el VDOP
para la constelación C, mientras que la Figura 2.16 y la Figura 2.17 muestran
ambos valores para la constelación H. En cada figura, una cruz marca la
ubicación de una estación.
Los valores de HDOP y VDOP han sido calculados, asumiendo que la posición
del blanco se encuentra en una rejilla de 300km por 300km, con una
granularidad de 1km, y una altitud de 10km sobre el nivel del mar.
52 GDOP (Geometric Dilution Of Precision)
53 LORAN (LOng RAnge Navigation) es un sistema de navegación terrestre que permite
mediante un receptor, que los barcos y las aeronaves estimen su propia posición y velocidad a
partir de señales de baja frecuencia transmitidas por estaciones terrestres.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
47
Figura 2.14. Mapa de HDOP para la constelación C.
Figura 2.15. Mapa de VDOP para la constelación C.
Figura 2.16. Mapa de HDOP para la constelación H.
CAPÍTULO 2
48
Figura 2.17. Mapa de VDOP para la constelación H.
De estas figuras, puede constatarse que un error en la medida de TDOA, es a-
priori más perjudicial en la constelación H que en la constelación C en la
mayoría de los puntos de la región bajo estudio.
Esta constatación puede explicarse desde un punto de vista matemático: en las
cercanías del centro de la constelación, los TOA de la señal emitida por la
aeronave son similares. El sistema de ecuaciones de los TDOA tiene menos
ecuaciones linealmente independientes. Cuando el blanco se encuentra en el
centro de la constelación, el sistema únicamente dispone de una ecuación.
En una situación real, la elección de la constelación C (en detrimento de la
constelación H) se efectúa en la etapa de planificación de un sistema WAM.
2.3.1.2.3. Métodos de resolución
[24] recopila los métodos de resolución del sistema de ecuaciones (2.8). Estos
métodos son de dos tipos:
• un método iterativo, presentado en [23, 25]
• métodos algebraicos, como los presentados en [21, 26]
El lector puede encontrar el desarrollo de los algoritmos en el ANEXO II.
2.3.1.2.3.1. Método iterativo
Cronológicamente, el primer método en aparecer fue el método iterativo de
resolución de ecuaciones no lineales [23, 25]. El método, basado en el
desarrollo en serie de Taylor de las ecuaciones del sistema (2.8), empieza con
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
49
una estimación inicial que se refina mediante la minimización del error de
mínimos cuadrados asociado a la aproximación lineal.
Este método asume que el error de linearización es limitado. Sin embargo, [26]
mostró que los términos de orden alto son significativos a la hora de determinar
la solución en escenarios con GDOP desfavorable. Así, aunque la relación
señal a ruido sea elevada, no hay garantías acerca de que la solución obtenida
sea precisa.
Resumiendo, las desventajas más importantes de este método son:
• la necesidad de una estimación inicial
• el hecho de que la convergencia no está garantizada
• el coste computacional del proceso
Por otra parte, este método presenta las siguientes ventajas:
• las medidas independientes a una estación pueden ser ponderadas.
• la información obtenida a través de las medidas no basadas en TOA -
sino en DOA (direcciónde llegada) – se combinan correctamente, es
decir, con el factor geométrico adecuado, al poder ser ponderadas con
sus precisiones a-priori.
• la distribución estadística de la solución (DOP) puede calcularse de
manera sencilla.
• es sencillo detectar la no convergencia.
2.3.1.2.3.2. Método algebraico (Bancroft)
[26] presenta un algoritmo empleado originalmente para el cálculo inicial de la
posición de un receptor del sistema GPS, basado en la localización hiperbólica,
y que por ello podemos reutilizar.
Este método ofrece dos posibles soluciones a la ecuación cuadrática, de la que
únicamente una es la correcta. La principal ventaja de este método es que es
algebraico y no iterativo, siendo numéricamente más eficiente.
Sin embargo, las simulaciones muestran que la precisión en el eje vertical es
paupérrima, siendo utilizable la solución obtenida únicamente en el plano
horizontal (latitud, longitud). Esto es debido a que las estaciones de WAM están
ubicadas en la superficie de la Tierra.
2.3.1.2.3.3. Método algebraico (Chan, Ho)
Por su parte, [21] aproxima el estimador de máxima verosimilitud cuando el
ruido asociado a la medida del TDOA es pequeño. Las ecuaciones TDOA se
CAPÍTULO 2
50
transforman en una serie de ecuaciones lineales, asumiendo que la posición de
la aeronave y la distancia de la aeronave a la estación de referencia son
magnitudes independientes. Esto es posible asumiendo que la aeronave se
encuentra lejos de la constelación. El método de mínimos cuadrados ofrece
una primera solución a este problema.
En el caso de que ambas magnitudes sean comparables, se introduce la
relación entre ellas, y se obtiene la solución al problema. Sin embargo, las
simulaciones han mostrado errores de gran magnitud empleando esta segunda
iteración (campo cercano). No en vano, [24] sugiere que se emplee únicamente
la solución de campo lejano, para luego refinarlo con el algoritmo iterativo.
Por estos motivos, el cálculo empleado para obtener la posición del blanco – y
así inicializar el sistema de seguimiento (3.4.2) – se efectúa en dos etapas: la
primera consiste en el cálculo de campo lejano, y la segunda, en un algoritmo
iterativo.
2.3.1.2.4. Cota de Cramér-Rao
En estadística, la cota inferior de Cramér-Rao (CRLB54) expresa una cota
inferior para la covarianza de un estimador insesgado. Según [27], la expresión
del CRLB para un sistema de localización hiperbólico es:
( )[ ] 11CRLB −−= MCMRMMC TTT ( 2.9 )
donde:
−−−
−−−
−−−
=
Nj
Nj
Nj
Nj
Nj
Nj
j
j
j
j
j
j
j
j
j
j
j
j
R
zz
R
yy
R
xx
R
zz
R
yy
R
xx
R
zz
R
yy
R
xx
c
MMM
2
2
2
2
2
2
1
1
1
1
1
1
0
1
C
−
−
−
=
10001
0
0101
0011
OM
L
M
( 2.10 )
y R es la matriz de covarianza de las medidas de TOA.
54 CRLB: siglas del término en inglés Cramér-Rao Lower Bound. Este término es llamado así
en honor a Harald Cramér y Calyampudi Radhakrishna Rao
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
51
A partir de dicha matriz pueden calcularse, la cota inferior asociada a la
desviación típica de la componente horizontal de error, y la cota inferior
asociada a la desviación típica de la componente vertical de error:
( )
=
+=
3,3min
2,21,1min
CRLBabsV
CRLBCRLBH
( 2.11 )
Las siguientes figuras muestran las magnitudes Hmin (Figura 2.18) y Vmin
(Figura 2.19), con la constelación H, para una aeronave ubicada a 10km de
altitud, y una desviación típica en la medida del TOA de 10ns.
Figura 2.18. Cota inferior del error cuadrático medio – componente horizontal (m).
Figura 2.19. Cota inferior del error cuadrático medio – componente vertical (m).
Con el fin de comparar la eficiencia de los métodos de resolución algebraicos,
se han evaluado la desviación típica del error horizontal y vertical, empleando
CAPÍTULO 2
52
las mismas condiciones experimentales: constelación H, una única aeronave
ubicada a 10km de altitud, y desviación típica en la medida del TOA de 10ns).
Se observa que ambos métodos no son óptimos desde la perspectiva de la
cota de Cramér-Rao.
Por una parte, el método presentado por Bancroft [26] (Figura 2.20 y Figura
2.21), tiene un mejor rendimiento dentro de la constelación y a lo largo de cada
eje, empeorando cuando la componente x de la posición de la aeronave se
encuentra fuera de la constelación. Esta observación es válida tanto para la
componente horizontal, como para la componente vertical.
Figura 2.20. Método de Bancroft: cociente de la desviación típica de la
componente horizontal de error por Hmin.
Figura 2.21. Método de Bancroft: cociente de la desviación típica de la
componente vertical de error por Vmin.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
53
Por otra parte, el algoritmo para campo lejano empleado en el método
presentado por Chan y Ho [21] (Figura 2.22 y Figura 2.23), ofrece un mejor
rendimiento fuera de la constelación, y lejos de los ejes x e y. Sin embargo, su
rendimiento es nefasto para la componente vertical.
Figura 2.22. Método de Chan-Ho: cociente de la desviación típica de la
componente horizontal de error por Hmin.
Figura 2.23. Método de Chan-Ho: cociente de la desviación típica de la
componente vertical de error por Vmin.
2.3.2. Efecto de los errores sistemáticos
Los errores sistemáticos son sesgos en la medida que conllevan a que la
media de las medidas difiera del valor real del elemento medido. Cualquier
medida es proclive a tener errores sistemáticos, incluso de distinta naturaleza.
CAPÍTULO 2
54
Los errores sistemáticos pueden estar presentes a su vez, en el resultado de
una estimación basada en un modelo matemático o una ley física.
El apartado siguiente repasa los conceptos básicos de la propagación
troposférica, para entender su impacto en la medida del tiempo de llegada.
A partir de ahí, estudiaremos cómo se puede modelar dicho impacto.
2.3.2.1. Propagación troposférica
2.3.2.1.1. Propiedades de la atmósfera terrestre
Un aspecto fundamental en la propagación es el índice de refracción
radioeléctrica (n). Si bien en la Tierra dicho valor es muy cercano a la unidad,
[28] establece el concepto de coíndice de refracción (N). La relación entre
ambos es la siguiente:
6101 −⋅+= Nn ( 2.12 )
donde N depende de la presión atmosférica (P, en hPa), la temperatura (T) y la
presión del vapor de agua (e, en hPa). El comportamiento de la temperatura, y
de la presión del vapor de agua en función de la altitud pueden apreciarse en la
Figura 2.24.
Figura 2.24. Temperatura y presión de la atmósfera en función de la altitud.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
55
Por ello, [28] recomienda utilizar la siguiente expresión para modelar N en
función de la altitud (Figura 2.25):
( )
−⋅⋅+= −
0
6
0 exp101 h
hNhn
( 2.13 )
siendo N0=315 y h0=7,35km.
Figura 2.25. Valor del coíndice de refracción (N) en función de la altitud.
2.3.2.1.2. Método de trazado de rayos
El problema básico de la refracción es la predicción de la trayectoria de un rayo
en un medio de propagación refractivo. Se puede probar que si la dirección del
rayo es conocida en un punto concreto del espacio, y el índice de refracción es
conocido en cualquier punto del espacio, se puede determinar la trayectoria
completa del rayo.
La ley de Snell relaciona la dirección del rayo en cualquier punto a lo largo de
su trayectoria, con la dirección del gradiente del índice de refracción (Figura
2.26). Si el medio tiene un gradiente con una dirección constante, dicha ley
enuncia que si ψ es el ángulo entre el rayo y la dirección del gradiente del
índice de refracción, se cumple:
constantesin =ψn ( 2.14 )
CAPÍTULO 2
56
Figura 2.26. Ley de Snell entre dos medios con diferente índice de refracción.
Según [29], el gradientedel índice de refracción de la atmósfera puede
considerarse vertical. Así, empleando un modelo esférico de la Tierra, cuyo
origen está situado en el centro de la misma, se puede reformular la ley de
Snell tal y como sigue:
Csin =ψnr ( 2.15 )
donde r es la distancia euclídea que separa un punto dado del rayo y el centro
de la Tierra.
Definiendo el ángulo θ – llamado ángulo de elevación local del rayo - como el
ángulo complementario de ψ, la expresión puede reescribirse:
000 coscos θθ rnnr = ( 2.16 )
donde los valores con subíndice son los respectivos valores iniciales.
Caracterizando la atmósfera mediante un modelo de estratificación por capas
(donde sus propiedades únicamente varían con la altitud), y según [28, 30], un
haz radioeléctrico que atraviesa la porción inferior (no ionizada) de la atmósfera
experimenta curvaturas debidas al gradiente del índice de refracción.
Así, la propagación de la onda radioeléctrica a través de la troposfera hace que
la distancia recorrida por la luz sea mayor que la distancia que separa el
transmisor del receptor. Este recorrido (s) se calcula mediante la integración
numérica de la siguiente expresión, descrita en [29, 31] :
( ) ( )
( ) ( )( )[ ]{ }
∫
+−
=
jz
rzhnn
dzzn
zs
0 2
000
0
1cos1
,
θ
θ
( 2.17 )
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
57
donde θ0 es el ángulo de elevación del rayo en el receptor, n0 el índice de
refracción en el receptor, r0 el radio de curvatura de la Tierra, y zj la altitud de la
aeronave.
Dicha integral se resuelve numéricamente, tal y como se describe en [29].
La Figura 2.27 muestra la distorsión atmosférica, definida como diferencia
entre la trayectoria recorrida por la onda (tiempo de vuelo dividido por la
velocidad de la luz en la atmósfera estándar a nivel del mar) y la distancia
euclídea (∆P) debido a la propagación troposférica55. De dicha figura, puede
observarse lo siguiente:
• para una determinada altitud del transmisor, la distorsión atmosférica
puede aproximarse mediante una función parabólica de la distancia
entre transmisor y el receptor.
• para una determinada distancia euclídea, la distorsión atmosférica es
mayor cuanto menor sea la altitud del transmisor.
2.3.2.1.3. Modelo de la distorsión atmosférica
Como continuación al modelo de primer orden empleado en [32], [33] sugiere el
uso de un modelo de orden 2, argumentando que captura mejor la dependencia
de ∆P con la distancia euclídea para distancias mayores de 100km.
Así, [33] establece la siguiente ecuación:
( ) 2
21
0
, ijij
zz
jij
j
i RRzRP
j
αα +≈∆
=
( 2.18 )
donde Rij es la distancia euclídea entre transmisor (aeronave) y receptor
(estación base), mientras que los parámetros α1 y α2 son los coeficientes de
orden uno y dos.
Un trabajo posterior [17] compara las diferentes ramas de la familia de
parábolas, efectuando el cociente (C) de la división de ∆P a una altitud
determinada entre ∆P para una altitud de 14000m, para una misma distancia
euclídea.
Dicho cociente está representado en función de la distancia euclídea en la
Figura 2.28, pudiendo apreciarse que para una determinada altitud de la
aeronave, el cociente es prácticamente constante. Cabe destacar que la línea
roja discontinua en las dos figuras siguientes muestra el horizonte
radioeléctrico de la señal transmitida, siendo éste el límite de aplicabilidad del
cálculo del sesgo.
55 En este análisis, se asume que el receptor se encuentra al nivel del mar
CAPÍTULO 2
58
0 50 100 150 200 250 300 350 400 450 500
0
10
20
30
40
50
60
70
80
90
100
Slant range (km)
R
an
g
e
b
ia
s,
m
3km
4km
5km
6km
7km
8km
9km
10km
11km
12km
13km
14km
Limit of applicability
(Radio−wave horizon)
Figura 2.27. Diferencia entre la trayectoria recorrida por una onda radioeléctrica y la distancia
entre transmisor y receptor, en función de la altitud del transmisor.
La Figura 2.29 es otra representación del mismo cociente, pero esta vez, en
función de la altitud de la aeronave.
15 50 100 150 200 250 300 350 400 450 500
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
R
at
io
b
et
w
ee
n
t
h
e
R
an
g
e
b
ia
s
at
a
g
iv
en
a
lt
it
u
d
e
co
m
p
ar
ed
t
o
t
h
e
R
an
g
e
b
ia
s
at
1
4
k
m
a
lt
it
u
d
e
(e
q
u
al
s
la
n
t
ra
n
g
e
co
n
si
d
er
ed
),
n
o
u
n
it
Slant range, km
3km
4km
5km
6km
7km
8km
9km
10km
11km
12km
13km
Limit of applicability
(Radio-wave horizon)
Figura 2.28. Cociente de la división de ∆P para una aeronave a una altitud determinada entre
∆P para una aeronave a una altitud de 14000m, para una misma distancia euclídea.
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
59
Figura 2.29. Cociente de la división de ∆P para una aeronave a una altitud determinada entre
∆P para una aeronave a una altitud de 14000m, para una misma distancia euclídea.
En la Figura 2.29, establecemos una interpolación lineal para los valores de
altitud entre 3000m y 14000m, obteniendo los siguientes parámetros:
• función: C = -0.075z + 2.000
• coeficiente de determinación56: R2=0.982
con lo que podemos modelar C mediante una función afín, tal y como se
enuncia en [17]:
( )
( )
( )
−+≈
∆
∆
=
=
=
0
3
0
11
,
,
z
z
zRP
zRP
zC
j
ij
j
i
jij
j
i
j
kij
kij α
ρρ
ρρ
( 2.19 )
Así, la contribución de la propagación troposférica a la pseudodistancia puede
modelarse mediante el producto de las ecuaciones (2.8) y (2.9), resultando:
( ) ( ) ( ) ( )
−+⋅+≈⋅∆≈∆
=
0
3
2
21 11,,
0 z
z
RRzCzRPzRP
j
ijijj
zz
jij
j
ijij
j
i
j
ααα
( 2.20 )
A través de esta expresión, podemos observar que el sesgo debido a la
propagación troposférica puede modelarse mediante el producto de dos
factores: el primero se trata de una corrección para una altura de referencia,
que depende únicamente de la distancia euclídea; el segundo factor
corresponde a la corrección de la altura, dependiendo α3 de la altura de
referencia.
56 Este coeficiente indica la calidad del ajuste de una serie de datos a un modelo estadístico. Si
el ajuste es perfecto, este valor es 1, mientras que si el valor es 0, la serie de datos es
independiente del modelo.
CAPÍTULO 2
60
2.3.2.2. Error de sincronismo
El error de sincronismo existente entre dos estaciones base en un determinado
instante - ∆T(t) - puede definirse como la suma de los siguientes términos:
• el error existente en una observación previa: ∆T(t0)
• la deriva del reloj entre la observación previa (t0) y el instante actual (t):
∆T(t)-∆T(t0)
Asumimos que las estaciones base emplean osciladores de cuarzo termo-
estabilizado (OCXO57), cuya deriva del reloj es de 100ps para un intervalo de
medida de un segundo [34].
De este modo, en un intervalo de un segundo, la deriva de un OCXO puede
inducir un error en la pseudodistancia de menos de 3cm (caso peor).
Así, si bien [33] presenta un polinomio de segundo orden para caracterizar ∆T
basado en las evaluaciones de [35], [17] sugiere que la magnitud de la deriva
(determinada por los coeficientes de orden uno y orden dos) durante el
intervalo de estimación de los errores sistemáticos sea limitada.
Por ello, esta tesis caracteriza el error de sincronismo entre dos estaciones
mediante un escalar, tal como se define en la siguiente expresión:
1,01, ii TcS ∆⋅≈∆
( 2.21 )
donde ∆Ti,1 es el error de sincronismo existente entre el reloj de una estación
determinada y el reloj de la estación de referencia.
2.3.3. Concepto de pseudodistancia
Con el fin de aglutinar los efectos de los errores definidos en la sección 2.3, se
define el concepto de pseudodistancia. La pseudodistancia es el producto de la
velocidad de la luz en el vacío por el tiempo de llegada de una señal, y su
unidad es el metro.
Dicho valor únicamentepuede ser igual a la distancia que separa al transmisor
del receptor en el caso ideal, mencionado en 2.2.2.
Matemáticamente, la pseudodistancia puede expresarse mediante la suma de
los siguientes términos:
• la distancia euclídea entre transmisor y receptor (Rij)
• el efecto de la incertidumbre en el tiempo de emisión de la señal (εe)
57 Del ingles, Oven-Controlled Crystal Oscillator
FUNCIONAMIENTO BÁSICO DEL SISTEMA WAM
61
• el efecto de los errores sistemáticos, compuestos por la propagación
troposférica (εp) y el error de sincronismo (εs)
• el efecto de los errores aleatorios (εr)
ij
r
i
s
ij
p
j
eij
j
i
j
iij RcTOAc εεεετρ ++++=⋅=⋅= 00
( 2.22 )
Dado que la medida principal es el TDOA (entre dos estaciones) – equivalente
a la diferencia de pseudodistancias - de la misma señal emitida por una
determinada aeronave, el término εe desaparece, tal como se aprecia en la
expresión (2.13):
( ) ( ) ( ) ( ) ( )jrijrsisjpijpjejejijjij RR 11111 εεεεεεεερρ −+−+−+−+−=− ( 2.23 )
que puede ser reescrita del siguiente modo:
( ) ij
r
i
s
ij
pjijjij RR '''11 εεερρ +++−=−
( 2.24 )
donde εr’ es una v.a. gaussiana con el doble de varianza que εr, debido a que
se trata de la diferencia entre dos v.a. gaussianas de misma varianza.
CAPÍTULO 2
62
2.4. Planteamiento del problema matemático
Incorporando (2.18) y (2.19) en los términos de error de (2.24), obtenemos la
siguiente expresión:
( ) ( ) ( )[ ] jiijjjijijjij
jij
j
i
nTc
z
z
RRRRRR
TDOAc
1,1,0
0
3
2
1211
2
211
11,0
11 +∆⋅+
−++−++−
=−=⋅
ααααα
ρρ
( 2.25 )
donde c0 es la velocidad de la luz en el vacío y z0 es una altitud de referencia.
Así, nuestro problema de localización se convierte en un problema de
estimación estadística de los siguientes parámetros:
• las posiciones de las aeronaves (xj,yj,zj)
• los parámetros empleados en el modelo de distorsión atmosférica (αk)
• las estimaciones de los errores de sincronismo (c·∆Ti).
donde la función de transferencia es una función no lineal.
Adicionalmente, es necesario que la solución obtenida minimice el efecto del
error aleatorio, modelado mediante una v.a. gaussiana de media nula.
Este capítulo ha presentado el sistema WAM, ha planteado el caso ideal, y ha
analizado los errores que afectan a la determinación de la posición,
desembocando en el planteamiento del problema matemático. El siguiente
capítulo describe la solución a este problema.
63
CAPÍTULO 3. SOLUCIÓN AL PROBLEMA
DE MULTILATERACIÓN EN PRESENCIA
DE ERRORES DE PROPAGACIÓN Y
SINCRONISMO
3.1. Necesidades del sistema
El sistema bajo estudio ha de ofrecer una solución al problema planteado en el
CAPÍTULO 2. De este modo:
• el sistema ha de ser capaz de localizar y seguir las pistas de las
aeronaves, basándose en las señales que éstas han emitido.
• el sistema ha de ser capaz de mitigar el efecto de la propagación
troposférica, modelado mediante una función no lineal.
• el sistema ha de ser capaz de corregir los errores de sincronismo entre
las estaciones base.
• el sistema ha de ser capaz de reducir la varianza del error de
posicionamiento integrando medidas de una misma aeronave en
distintos instantes (mediante el filtrado de trayectorias).
• las únicas entradas al sistema serán los TOA de las señales emitidas por
el tráfico de oportunidad (se supone que no existen emisores de
calibración).
CAPÍTULO 3
64
3.2. Análisis del sistema de ecuaciones para el sistema
de localización hiperbólico
En un escenario genérico con una cantidad M de aeronaves y una cantidad N
de estaciones base en el área de cobertura, el sistema de ecuaciones (2.23)
está compuesto por:
• M(N-1) ecuaciones, al haber (N-1) medidas de TDOA por cada aeronave
• 3M+N+2 incógnitas:
o 3M provenientes de la posición (xj,yj,zj) de cada aeronave
o 3 debido a los parámetros α1, α2 y α3.
o N-1 provenientes de los errores de sincronismo
El apartado siguiente muestra dos métodos de resolución del problema para las
muestras de TDOA obtenidas en un instante.
3.3. Solución al problema muestra a muestra
Este apartado describe las dos técnicas presentadas en [32] para resolver el
sistema de ecuaciones (2.23).
3.3.1. Mediante calibración
La primera técnica consiste en determinar primero los parámetros del modelo
de distorsión atmosférica (α1, α2 y α3) a partir de las medidas de calibración, es
decir, a partir del TDOA de las señales generadas por una serie de emisores
fijos en posiciones conocidas. Una vez determinados, estos parámetros se
aplican a las ecuaciones asociadas a las medidas de TDOA de las señales
generadas por cada aeronave, con el fin de calcular su posición.
Esta técnica presenta dos inconvenientes: por una parte, es necesario que
cada emisor de calibración tenga línea de vista con cada estación base de la
constelación WAM. Por otra parte, debido a que dichos emisores son fijos y se
encuentran a muy baja altura, los parámetros α1, α2 y α3 se estiman a partir de
una serie de medidas que no permiten caracterizar los errores de propagación
a alturas mayores que las de la superficie de la tierra.
Por ello, esta técnica no es útil para nuestro problema. Es necesario recurrir a
un método que emplee las señales enviadas por las propias aeronaves. Este
método se presenta a continuación.
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
65
3.3.2. Mediante tráfico de oportunidad
3.3.2.1. Descripción de la técnica
La segunda técnica consiste en aumentar el número de estaciones base, y por
ende, el número de medidas en un instante. Con este aumento, se pretende
estimar simultáneamente la posición de cada aeronave, los errores de
sincronismo asociados a cada estación, y los parámetros α1, α2 y α3.
Con el fin de evitar que dicho sistema sea indeterminado, es necesario que el
número de ecuaciones sea mayor o igual que el número de incógnitas. Esto es,
la siguiente desigualdad ha de cumplirse:
( ) 231 ++≥− NMNM
( 3.1 )
es decir:
−
+
≥
4
2
N
N
M
( 3.2 )
La expresión (3.2) introduce dos restricciones por si sola:
• el sistema ha de estar compuesto por al menos 5 estaciones base
• el número mínimo de aeronaves depende del número de estaciones
base, según la relación definida en la Tabla 3.1
N valor mínimo de M
5 7
6 4
7 3
8 3
Tabla 3.1. Relación entre N y el valor mínimo de M.
De este modo, no es suficiente con tener una única aeronave para poder
corregir los efectos de los errores sistemáticos. Afortunadamente, desde la
perspectiva de la vigilancia aérea para el tráfico en ruta (donde es primordial
garantizar la separación entre aeronaves), la realidad juega a nuestro favor: en
caso de tener un número muy reducido de aeronaves en la zona de cobertura,
el efecto del error sistemático será de poca importancia en el control y
monitorización de posición, ya que la distancia entre las propias aeronaves
será grande. En el caso opuesto, una región con tráfico denso cumple
sobradamente con la restricción.
CAPÍTULO 3
66
Apoyándonos en este razonamiento, el sistema a desarrollar deberá basarse
en las medidas de TDOA obtenidas a partir del tráfico de oportunidad.
3.3.2.2. Descripción de la solución
Dado que la función que relaciona la medida de TDOA con la posición de la
aeronave no es lineal, [32] abre una línea de investigación con el fin de obtener
la solución al sistema tal y como sugiere [24].
Esta propuesta consiste en obtener la solución para cada punto de la
trayectoria en dos tiempos. Primero, se trata de encontrar una solución inicial
basada en la localización hiperbólica [21], tal y como se menciona en 2.3.1. Si
bien este paso obvia la existencia de errores sistemáticos, la estimación inicial
que se obtiene es un buen punto de partida para el siguiente paso.El segundo paso consiste en el uso de un algoritmo iterativo basado en la
minimización del error cuadrático medio (MMSE58) [25]. Este método parte de
la solución inicial anterior, y a través de una aproximación de primer orden de la
expresión (2.23), trata de converger a la solución.
En un instante determinado, el vector de incógnitas es el siguiente:
[ ]TNMMM TcTczyxzyx 3211,1,2111 ,,,,...,,,,,...,,, ααα∆∆=x
( 3.3 )
donde (xj,yj,zj) es la posición de la aeronave j-ésima, c∆Ti,1 es el error de
sincronismo entre la estación i-ésima y la estación de referencia (estación 1) α1,
α2 y α3 son los parámetros del modelo del error de propagación.
Cada iteración actualiza el vector definido en (3.3) del siguiente modo:
kkk xδxx += −1 ( 3.4 )
En cada iteración, se estiman las TDOA a partir del vector 1−kx . El vector de
error, definido como la diferencia entre las medidas de TDOA y las
estimaciones, es:
( )1−−= kk xfε τ ( 3.5 )
donde ( )1−kxf es el vector compuesto por las estimaciones de TDOA, y τ el
vector compuestos por las medidas de TDOA.
El incremento diferencial ( kxδ ) se obtiene resolviendo los términos de primer
orden del desarrollo en serie de Taylor del sistema de ecuaciones (2.23). Así, la
expresión (3.5) puede reescribirse del siguiente modo:
58 Minimum Mean Square Error
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
67
kkk
δxJε = ( 3.6 )
donde la matriz Jk es el jacobiano del sistema de ecuaciones (2.23) en la
iteración k-ésima.
Para un instante determinado, la matriz J se obtiene empleando las siguientes
expresiones59:
∇
∇
∇
∇
=
,1sN
i,1
3,1
2,1
,1sN,1sN,1sN
i,1i,1i,1
3,13,13,1
2,12,12,1
Z
Z
Z
Z
BA100
0
BA000
BA010
BA001
J
M
L
MMLMMM
L
L
L
( 3.7 )
donde:
=∇
A/CN
i,1i,1i,1
i,1
j
i,1i,1
i,1i,1
1
i,1
i,1
∆00
0∆0
00∆
L
MLMM
L
L
[ ]000=i,10
∂
∂
−
∂
∂
∂
∂
−
∂
∂
∂
∂
−
∂
∂
=
j
j
j
ij
j
j
j
ij
j
j
j
ij
zzyyxx
111
,,
ρρρρρρj
i,1∆
( )
jiji zr 31, 1 α+=A
( )
jiji zr 3
2
1, 1 α+=B
( )2211, ijijji rrz αα +=Z
( 3.8 )
siendo:
( ) ( )
j
ij
ij
j
ij
z
r
xx
x
32
1 12
1
αα
αρ
+
+
+
−=
∂
∂
( 3.9 )
59 El índice i corresponde a la estación i-ésima, mientras que el índice j corresponde a la
aeronave j-ésima
CAPÍTULO 3
68
( ) ( )
j
ij
ij
j
ij
z
r
yy
y
32
1 12
1
αα
αρ
+
+
+
−=
∂
∂
( ) ( ) ( )2213321 12
1
ijijj
ij
ij
j
ij
rrz
r
zz
z
ααααα
αρ
+++
+
+
−=
∂
∂
( ) ( ) ( )222
ijijijij zzyyxxr −+−+−=
Dado que este sistema suele ser sobredeterminado, se utilizará una solución
de minimización del error cuadrático medio para su solución. El incremento
diferencial se calcula tal y como sigue:
[ ] [ ] [ ] [ ]( ) [ ] [ ] kTT kkkkkk εSJJSJδx ⋅⋅⋅⋅⋅= −−− 111 ( 3.10 )
donde S es la matriz de covarianzas asociada a las medidas de TDOA.
3.3.2.3. Ventajas e inconvenientes de la solución
Trabajos como [32, 33] han buscado la forma de encontrar la solución a dicho
problema centrándose únicamente en las medidas de TDOA en un instante de
tiempo. Esto tiene la ventaja de que no es necesario tener conocimiento de las
trayectorias de los blancos. Sin embargo, el uso de este método con el fin de
efectuar el seguimiento de aeronaves contiene una serie de problemas.
3.3.2.3.1. Problema de convergencia
Cuando la magnitud del error aleatorio es importante, el algoritmo iterativo
efectuado en los procesadores tiene dificultades a la hora de converger.
El cálculo del paso (diferencia entre el vector solución entre una iteración y la
siguiente) requiere la inversión de matrices, como es el caso de la matriz de
gradiente A, definida en la expresión (15) de [33].
Supongamos el siguiente sistema de ecuaciones, donde el vector x es la
incógnita:
yxH =
( 3.11 )
H, x, y son definidos de la siguiente manera:
=
3
2
1
2
1
3231
2221
1211
y
y
y
x
x
hh
hh
hh
( 3.12 )
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
69
Siendo conceptualmente diferentes las magnitudes descritas por cada
dimensión del vector de solución x60, es posible que dichos valores difieran en
órdenes de magnitud. En esta situación, es muy probable obtener una matriz
mal condicionada, dificultando así la convergencia del algoritmo de MMSE.
Un método sencillo para mitigar este problema consiste en ponderar una de las
columnas de H, y la dimensión del vector x asociada, con el fin de no alterar la
ecuación, tal como aparece en la siguiente expresión:
=
⋅
3
2
1
2
1
32
31
22
21
12
11
y
y
y
xK
x
K
h
h
K
h
h
K
h
h
( 3.13 )
3.3.2.3.2. Problema a nivel conceptual
Analizando esta solución a nivel conceptual, nos encontramos ante un método
que - en caso de converger - obtiene el valor óptimo del vector x para un
determinado vector de medidas de TDOA. Sin embargo, no existe relación
alguna entre las componentes del vector x entre un instante y el siguiente.
De este modo, si el nivel de ruido (error en la medida) es suficientemente
elevado, es posible que tanto la posición de cada aeronave, como los
parámetros que caracterizan el modelo de la propagación troposférica, o
incluso los errores de sincronismo, varíen significativamente de una serie de
medidas a la siguiente (menos de 5 segundos de tiempo de experimento).
Se pueden mejorar las prestaciones si se incorpora un filtro de seguimiento
para relacionar estas medidas y las de las siguientes exploraciones (o si el
sistema es pasivo, de las siguientes respuestas de cada aeronave).
3.4. Solución al problema de seguimiento
3.4.1. Técnicas de seguimiento de blancos
La mayoría de las técnicas de seguimiento de blancos están basadas en
derivados del filtro de Kalman[36]. Dicho filtro opera de manera recursiva,
empleando las medidas (formalmente denominadas “observaciones”) con el fin
60 Un ejemplo es que una dimensión del vector solución sea la posición del blanco en metros,
mientras que otra dimensión del mismo vector sea un parámetro α asociado al modelo de
propagación troposférica
CAPÍTULO 3
70
de producir una estimación óptima (desde el punto de vista estadístico) del
estado del blanco.
De esta forma, el filtro deberá estar compuesto por la combinación de un
modelo de la dinámica del blanco y de un modelo de observación del mismo.
Mientras que el primero permite predecir el estado del blanco en un instante a
partir del anterior, el segundo estima la observación (sin ruido) partiendo del
estado del blanco.
En nuestro problema, el vector de estado del blanco está compuesto por: su
posición, su velocidad, y su aceleración en un sistema ENU, donde la latitud y
la longitud corresponden al centro de la constelación WAM, y la altitud es nula
respecto al elipsoide WGS84.
La clave para efectuar un seguimiento correcto consiste en extraer la
información necesaria del vector de estado del blanco (posición, velocidad,
etc.) a partir de las observaciones. En este sentido, un buen modelo del blanco
facilita en gran medida el proceso de extracción de la información. La elección
de dicho modelo tiene más peso en el éxito (o fracaso) del sistema de
seguimiento, cuanto más limitadas sean las observaciones del blanco. Debido
a su importancia, la discusión acerca de los modelos de dinámica de blanco
acontecerá en 3.4.1.1.
Por otra parte, en el marco de esta tesis, el modelo de observación es
justamente la función que estima la medida de TDOA a partir de la posicióndel
blanco. Puesto que dicha función no es lineal, tal como sugiere la expresión
(2.23), el apartado 2.3 analiza los distintos métodos existentes para aproximar
el comportamiento de esta función.
En el marco de esta tesis, tanto el cálculo de la posición del blanco como las
observaciones se ejecutan con un intervalo aproximadamente constante. Por
ello, podemos emplear las expresiones que caracterizan la dinámica del blanco
como las observaciones, en tiempo discreto.
3.4.1.1. Modelos de dinámica del blanco
Siguiendo la literatura [37, 38], la expresión clásica del modelo de dinámica del
blanco en tiempo discreto es:
[ ] [ ] [ ]( ) [ ]kwkukxfkx k +=+ ,1 ( 3.14 )
donde [ ] [ ] [ ]kwkukx ,, son respectivamente, el vector de estado, el vector de
control del blanco, y el error de proceso en el instante de tiempo k. A nivel de
notación, hay que destacar que el vector de estado del blanco [ ]kx , está
compuesto por la posición, la velocidad, y la aceleración a lo largo de una
dirección concreta.
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
71
[ ] ( )Tkxkxkxkx ][],[],[ &&&= ( 3.15 )
Hay que especificar que el blanco se considera como un punto en el espacio.
Esto simplifica el enfoque acerca de las ecuaciones de movimiento, ya que
únicamente es necesario caracterizar la dinámica del punto, sin mencionar la
dinámica de momentos.
Tal y como menciona [38], el proceso de control definido mediante el vector
u[k] - que podría contener información como la posición del acelerador, de los
alerones o del timón - responsable de la maniobra del blanco es de naturaleza
determinista, pero es a menudo desconocido para el sistema de seguimiento.
Una consecuencia lógica es modelar los elementos de dicho vector de estado
mediante variables deterministas que serían estimadas a partir de las
observaciones efectuadas durante el seguimiento. Este razonamiento dio pie a
métodos de estimación de los parámetros de entrada [39-47].
Alternativamente, se puede modelar el proceso de control como un proceso
estocástico, idea que ha sido mayoritariamente adoptada en el campo del
seguimiento de blancos. Estos modelos pueden clasificarse en los siguientes
tres grupos: los modelos de ruido blanco, los modelos de Markov, y los
procesos semi-Markovianos.
3.4.1.1.1. Modelos de ruido blanco
Los modelos de ruido blanco, son aquellos donde el proceso de control está
modelado como ruido blanco. Esto abarca los modelos polinómicos (de grado n,
donde los coeficientes de grado n+1 son modelados mediante ruido).
Dos casos particulares importantes son los modelos de velocidad constante,
donde la aceleración se modela mediante ruido (n=1) y los modelos de
aceleración constante, donde la sobreaceleración61 se modela mediante ruido
(n=2). Si bien estos modelos son sencillos y de aplicación directa, los cambios
de rumbo, o de aceleración de una aeronave son debidos a las maniobras de
los pilotos (manuales o automáticas), y por ello, modelar la sobreaceleración
como ruido blanco no es lo más apropiado.
3.4.1.1.2. Modelos de Markov
Los modelos de Markov son aquellos donde el estado del proceso de control
estocástico depende del estado del mismo únicamente en m instantes
anteriores (sistemas de memoria m), siendo m=1 el caso más sencillo.
61 Conocida también como tirón, sacudida o pique, se trata de la tasa de cambio de la
aceleración, es decir, la derivada de la aceleración con respecto al tiempo.
CAPÍTULO 3
72
De entre ellos, destacamos el modelo de Singer [48], que asume que la
aceleración del blanco a(t) está correlada en el tiempo: es decir, si el blanco
está acelerando (es decir, en una maniobra) en un instante de tiempo t, es
probable que esté acelerando en un instante posterior (t+τ) para un valor
suficientemente pequeño de τ. Por ejemplo, un giro lento (típico de aviones
comerciales) dará lugar a valores correlados de aceleración de hasta un minuto,
una maniobra evasiva dará lugar a períodos de entre 10 y 30 segundos,
mientras que las turbulencias pueden dar lugar a períodos de 1 ó 2 segundos.
El modelo representativo de la función de autocorrelación de la aceleración es:
( ) ( ) ( )[ ] 0 , 2 >=+= − ασττ ταetataEr m ( 3.16 )
donde 2mσ es la varianza de la aceleración del blanco y α es la inversa de la
constante de tiempo de la maniobra.
[48] demuestra que la función definida en (3.16) corresponde al
comportamiento de la sobreaceleración modelada mediante el siguiente
sistema LTI62:
( ) ( ) ( ) 0 , >+−= αα twtata& ( 3.17 )
donde w(t) es un ruido blanco de media nula y densidad espectral de potencia
2
2 mασ . Así, la representación en tiempo continuo en modelo de variables de
estado del modelo de Singer es:
( ) ( ) ( )twtxtx
+
−
=
1
0
0
00
100
010
α
&
( 3.18 )
La expresión equivalente a (3.17) en tiempo discreto (siendo T el intervalo de
muestra) es la siguiente:
[ ] [ ] 0 , 1 >+=+ − αα ak
T
wkaeka
( 3.19 )
donde akw es un ruido blanco de media nula y varianza ( )Tm e ασ 22 1 −− , mientras
que la expresión equivalente a (3.18) en tiempo discreto es:
[ ] [ ] [ ]
( )
( ) [ ] [ ]kwkxkwkxFkx +
+
=+=+
T-
T-
2T-
e00
e-110
e1-T1
1
α
α
α
α α
αα
( 3.20 )
62 Linear Time-Invariant: lineal e invariante en el tiempo
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
73
La varianza 2mσ está distribuida mediante el modelo de densidad de
probabilidad de la Figura 3.1, tendiendo las siguientes características:
• el blanco acelera o decelera un valor Amax con una probabilidad Pmax.
• durante la maniobra, el blanco tiene aceleración nula con una
probabilidad P0.
• para los valores de aceleración definidos en el intervalo ( ) { }0, maxmax −− AA ,
la probabilidad viene definida por una densidad de probabilidad uniforme
de valor ( ) maxmax0 221 APP −− .
Figura 3.1. Función densidad de probabilidad de la aceleración en el modelo de Singer.
Siendo los valores Amax, Pmax y P0 valores de diseño del modelo, la varianza de
la aceleración 2mσ es:
( )0max
2
max2 41
3
PP
A
m −+=σ
( 3.21 )
En cuanto al comportamiento asintótico del modelo, la expresión (3.20) se
aproxima a la de un modelo de aceleración constante para valores grandes de
τ (αT<<1), mientras que para valores pequeños de τ (α>>T, αT>>1) se
aproxima a un modelo de velocidad constante.
De este modo, el modelo de Singer caracteriza la dinámica del blanco a caballo
entre un modelo de velocidad constante y un modelo de aceleración constante,
donde el producto αT determina la proximidad del comportamiento a una de
sus dos asíntotas.
Debido a esta versatilidad, el modelo de aceleración de Singer es un modelo
habitualmente empleado para el seguimiento de blancos, tanto para vigilancia
aérea como para el control de misiles aire-aire [17, 37, 49-55].
CAPÍTULO 3
74
3.4.1.1.3. Procesos semi-Markovianos
En la práctica, muchas maniobras de los blancos se deben a aceleraciones de
media no nula, que adoptan valores discretos en ciertos períodos de tiempo.
Sin embargo, la dificultad se encuentra en que el sistema de seguimiento no
conoce ni los instantes en los que la media de la aceleración no es nula, ni los
valores discretos que dicha aceleración adopta.
Uno de los procesos aleatorios constantes a trozos, es el llamado proceso de
salto semi-Markoviano. Si bien los instantes de los saltos son independientes
de los estados anteriores (propiedad de proceso Markoviano), el proceso no
será Markoviano durante los instantes de tiempo entre dichos saltos (tiempo de
permanencia). Así, su estado en un instante determinado depende de los
estados en los instantes anteriores durante el tiempo de permanencia.
Un ejemplo de dicho procesoes el modelo de aceleración que resulta de la
combinación de un modelo de salto y el modelo de Singer presentado en
3.4.1.1.2:
( ) ( ) ( ) ( )tatutvta ~++−= β ( 3.22 )
donde ( )ta~ es la aceleración del modelo de Singer, v es la velocidad, u(t) es la
media de la aceleración (desconocida) y β es el coeficiente de arrastre.
Así, la representación en tiempo continuo en modelo de variables de este
modelo es63:
( ) ( ) ( ) wtutxtx
+
+
−
−=
1
0
0
0
1
0
00
10
010
α
β&
( 3.23 )
De este modo, la aceleración del blanco se modela mediante un modelo de
aceleración de Singer con media no nula, donde la aceleración se modela
mediante un proceso semi-Markoviano.
En esta tesis no ahondaremos en estos procesos, remitiendo al lector a la
literatura [56-66].
3.4.1.2. Técnicas de aproximación para filtrado no lineal
El filtro de Kalman fue concebido asumiendo que la función que transforma el
estado del blanco con su observación y su modelo de planta (de movimiento)
son lineales. Sin embargo, la expresión (2.23) - que define la estimación de la
63 La expresión de x(t) viene definida en (3.15)
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
75
observación (vector de TDOA) a partir del estado (vector posición) de los
blancos - es no lineal.
Existen muchas técnicas de aproximación para la estimación del estado en
funciones no lineales. La mayoría de ellas puede agruparse en los tres grupos
siguientes:
• aproximación de la función de observación
• aproximación de los momentos mediante muestreo
• aproximación del modelo estocástico.
En esta tesis se muestran técnicas pertenecientes a los dos primeros grupos,
remitiendo al lector a [67] para el tercer grupo.
3.4.1.2.1. Aproximación de la función de observación
Las técnicas que aproximan la función de observación (determinista y no lineal)
están basadas en las aproximaciones mediante series de Taylor. El ejemplo
más conocido es el filtro de Kalman extendido (EKF64) de primer orden, que
aproxima las funciones no lineales (en nuestro caso, el interés está en la
función de observación) mediante funciones lineales.
Dicha aproximación es válida siempre que el error asociado a la estimación del
estado del blanco sea suficientemente pequeño. De no ser así, el filtro no
convergería, aumentando el error en la estimación.
Estos filtros y sus múltiples variantes han sido estudiados e implementados en
numerosas aplicaciones, existiendo una amplia literatura al respecto. Esta tesis
presenta el empleo de esta técnica para la estimación de los sesgos, en los
apartados 3.4.3.2, 3.4.4.2 y 3.4.4.4.
3.4.1.2.2. Aproximación de los momentos mediante muestreo
Los valores estadísticos más empleados para caracterizar una distribución de
probabilidad son la media y la covarianza. Las técnicas que aproximan los
momentos de una distribución intentan aproximar directamente la media y la
covarianza, en vez de la función de observación. Estas técnicas pueden ser
analíticas o numéricas (mediante muestreo).
La gran desventaja de las técnicas analíticas radica en que suelen depender
del problema en cuestión. La generalidad de las técnicas basadas en muestreo,
combinada con el aumento de la capacidad de procesado en las últimas
décadas ha propiciado que éstas hayan sido empleadas en sistemas de
seguimiento.
64 Extended Kalman Filter
CAPÍTULO 3
76
El método Unscented Transformation (UT) [68-70] aproxima los momentos (la
media y la covarianza) de un vector aleatorio (en nuestro caso, el TDOA) que
está caracterizado mediante una función no lineal del vector de estado,
empleando una serie de muestras deterministas, denominadas “puntos sigma”
(Χ), con sus pesos asociados. La Figura 3.2 compara los conceptos de
propagación de los dos primeros momentos de la distribución según el modelo
de Monte Carlo, la aproximación de la función de observación de primer orden
(base del en EKF) y la UT.
El filtro basado en esta técnica es el denominado Unscented Kalman Filter
(UKF) [71-77]. Dicho filtro, presentado en el ANEXO III, ofrece una solución
más robusta que el EKF cuando el modelo de observación es no lineal. Por
este motivo, se emplea el UKF para el seguimiento de blancos mediante TDOA
en [17]. La implementación se detalla en los apartados 3.4.4.1 y 3.4.4.3.
Figura 3.2. Concepto de propagación de la media y la varianza: Monte Carlo, EKF y UT [78].
3.4.2. Seguimiento de los errores sistemáticos
3.4.2.1. Modelo de dinámica de los errores sistemáticos
Los errores sistemáticos descritos en el apartado 2.3.2 pueden variar en el
tiempo debido a los cambios en las condiciones atmosféricas y a las derivas en
los relojes situados en las estaciones base. Es por ello necesario efectuar el
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
77
seguimiento de los errores sistemáticos de manera simultánea al seguimiento
de los blancos.
Con el fin de determinar el modelo de dinámica de los errores sistemáticos,
analizamos el efecto de los errores sistemáticos para el caso de una aeronave
volando a 250m·s-1 con altitud constante, a una distancia euclídea de 100km de
una estación base determinada. Asumiendo que la aeronave emite cada 4
segundos, y partiendo de los datos obtenidos a partir del análisis de errores
efectuado en el apartado 2.3, la contribución total de los sesgos es de 0,42m
(0,3m debido a la refracción65 y 0,12m debido al sincronismo66).
Debido a esta lenta variación que exhibe el sesgo, se emplea un filtro recursivo
de orden 1 para corregir el sesgo (deriva lineal).
3.4.2.2. Estimación de los errores sistemáticos
El problema con la estimación de dichos errores en tiempo real consiste en que
sus ecuaciones de filtrado están vinculadas a las ecuaciones de seguimiento
asociadas a cada blanco.
A-priori, esto implicaría el uso de un sistema de seguimiento donde el vector de
estado sería el resultante de la concatenación de los vectores de estado del
blanco y de los errores sistemáticos. Este vector de estado está basado en el
mencionado en (3.3), al que se añaden la velocidad y la aceleración:
[ ]TNMMM TcTcavxavx 3211,1,2111 ,,,,...,,,,,...,,, ααα∆∆=z ( 3.24 )
donde los primeros 3M elementos corresponden a las posiciones de cada
aeronave, los N-1 siguientes corresponden a los errores de sincronismo de
cada estación respecto a la primera, y los tres últimos elementos son los
parámetros del modelo de distorsión debido a la propagación troposférica.
Según [79], dicho vector puede expresarse como una concatenación del vector
de posiciones (s), y el vector de parámetros de error (∆):
[ ]TTT E, Xz ≡
[ ]TMMM avxavx ,,,...,,, 111≡X
( 3.25 )
65 En 4 segundos, la distancia euclídea varía de 100km a 101km, con lo que el sesgo en la
pseudodistancia debido a la refracción atmosférica varía en 0,3m.
66 La deriva de un OCXO es de 100ps para un intervalo de medida de 1s. De este modo, en 4
segundos, el sesgo en la pseudodistancia debido al error de sincronismo varía en 0,12m.
CAPÍTULO 3
78
[ ]TNTcTcE 3211,1,2 ,,,,..., ααα∆∆≡
Definimos la función de transferencia h, que estima las observaciones (TDOAs)
a partir de los valores del vector de estado z. Adicionalmente, definimos la
función de error aleatorio g. Así, a partir de (3.25), tenemos:
( ) ( ) ( ) iEE nX,gX,hzh ⋅+= ( 3.26 )
Empleando los términos de primer orden del desarrollo en serie de Taylor, se
obtiene:
( ) ( ) ( ) ( ) nEEE GFXXFzhzh X +−+−+≈ ˆˆˆ ( 3.27 )
donde FX es el Jacobiano de la función h respecto al vector X, mientras que FE
es el Jacobiano de la función h respecto al vector E. Por otra parte, G es el
Jacobiano de la matriz de errores aleatorios respectoa una variable aleatoria
gaussiana de media nula y covarianza unidad.
Esta aproximación permite desacoplar las ecuaciones de filtrado en dos pares
de ecuaciones: el primero para el seguimiento de los blancos, mientras que el
segundo para la estimación de los errores sistemáticos. Dichas etapas de
filtrado se efectúan mediante Extended Kalman Filter (EKF).
Aplicando este método a nuestro problema, el sistema de seguimiento puede
descomponerse en dos sistemas entrelazados: el primero, efectuando el
seguimiento de los blancos, mientras que el segundo efectúa el seguimiento de
los parámetros empleados para la estimación del sesgo – αi y ∆Ti en la
expresión (2.23).
En este caso, el conjunto de estaciones base proporciona las medidas TDOA
que permiten obtener la posición de cada blanco, dando lugar a una única
trayectoria. El sistema bajo estudio emplea así, la simplificación propuesta en
[80] para el escenario de un único sensor.
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
79
3.4.3. Arquitectura del sistema
Basándonos en las reflexiones efectuadas en 3.4.1 y 3.4.2, el sistema
(presentado en [17]) ha de efectuar el seguimiento tanto de las aeronaves,
como de los parámetros utilizados en la caracterización de los errores
sistemáticos.
Figura 3.3. Arquitectura del sistema.
Por ello, tal como se avanzó en 3.4.2 [79], el sistema está compuesto por dos
filtros de seguimiento entrelazados, ilustrados en la Figura 3.3:
• el primero (en verde), para los vectores de estado67 de las aeronaves
• el segundo (en naranja), para los parámetros empleados en la
caracterización de los errores sistemáticos
67 Compuesto por los vectores posición, velocidad y aceleración
CAPÍTULO 3
80
3.4.3.1. Seguimiento de las aeronaves
Dado que la función que caracteriza el TDOA en base a la posición de la
aeronave es no lineal, el sistema emplea el Unscented Kalman Filter (UKF) [81]
– presentado en 3.4.1.2.2 - para efectuar el seguimiento de las aeronaves.
La etapa de predicción de UKF – desarrollada en 3.4.4.1 - está basada en el
modelo de aceleración de Singer [48]. Las entradas a este modelo son el vector
de estado de la aeronave y su matriz de covarianza asociada. Por otra parte, la
etapa de filtrado – desarrollada en 3.4.4.3 - combina la salida de la etapa de
predicción (vector de estado y su matriz de covarianza), la estimación de los
parámetros característicos del error sistemático y las medidas de TDOA.
La inicialización de este filtro requiere una estimación inicial del vector de
estado. Este aspecto será abordado en el apartado 3.4.4.5.
3.4.3.2. Seguimiento de los errores sistemáticos
Por otra parte, el seguimiento del modelo de errores sistemáticos es efectuado
mediante Extended Kalman Filter (EKF). La fase de estimación del EKF -
3.4.4.2 - utiliza la combinación del vector de estado y su matriz de covarianza
resultantes de la muestra anterior. En este caso, el vector de estado contiene
los errores de sincronismo y los parámetros α característicos del error causado
por la propagación.
La inicialización de este filtro de seguimiento se efectúa mediante un vector de
error nulo a la entrada (es decir, una estimación inicial sin sesgo). Los valores
de la matriz de covarianza estarán definidos como el resultado del compromiso
entre un sistema reactivo y un sistema estable.
3.4.4. Implementación
Tal como se definió en el apartado 2.2.1.1, el sistema de referencia se
encuentra situado en el centro de la constelación WAM. El vector de estado (X)
está compuesto por la combinación de los vectores tridimensionales de
posición, velocidad y aceleración para cada aeronave (M en total). La matriz de
covarianza del blanco (PX) será inicialmente una matriz diagonal de bloques.
Estos son, respectivamente:
[ ]
( )
=
=
Mj XXXX
T
Mj
PPPdiagP
XXXX
.....
.....
1
1
( 3.28 )
siendo:
[ ]
( )
=
=
222222222
jjjjjjjjjj zzzyyyxxxX
T
jjjjjjjjjj
diagP
zzzyyyxxxX
&&&&&&&&&
&&&&&&&&&
σσσσσσσσσ
( 3.29 )
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
81
Del mismo modo, el vector de estado que contiene los parámetros del modelo
de error sistemático (E), está compuesto por los (3) parámetros empleados en
el modelo de distorsión troposférica y los errores de sincronismo entre los
relojes de la estación base (N-1, donde N es el número de estaciones base).
Así, el vector E y su matriz de covarianza asociada (PE) se definen a
continuación:
[ ]
( )
=
∆∆∆∆=
∆∆∆∆
2222222
1,1,1,31,2321
1,1,1,31,2321 Ni TcTcTcTcE
T
Ni
diagP
TcTcTcTcE
σσσσσσσ
ααα
ααα LL
LL
( 3.30 )
3.4.4.1. Etapa de predicción (seguimiento del blanco)
Tal y como se adelantó en el apartado 3.4.3.1, la etapa de predicción en el
seguimiento del blanco se implementa mediante un UKF[81]. Esta etapa se
efectúa en tres fases, tal como aparece en la Figura 3.4:
• la primera fase consiste en generar el vector de estado aumentado
• la segunda fase consiste en la predicción del vector de estado asociado
a cada “punto sigma”, empleando el modelo de aceleración de Singer.
• la tercera y última fase consiste en la generación del vector de estado y
de su matriz de covarianza asociada a la salida de la etapa de
predicción.
Figura 3.4. Diagrama de bloques de la etapa de predicción del blanco [78].
Dado que el vector de estado que caracteriza a cada aeronave está compuesto
por 9 elementos, la cantidad de puntos sigma es 18M+1. Estos puntos sigma χi
se obtienen mediante la siguiente expresión:
CAPÍTULO 3
82
( ) ( )
( ) ( ) ( )
( ) ( ) ( )
+=−−−=−−
=+−−=−−
−−=−−
−−
−−
MMqPkkXkk
MqPkkXkk
kkXkk
kkX
q
kkX
q
18,...,19,1111
9,...,1,1111
1111
11
0
11
0
00
0
0
γχ
γχ
χ
( 3.31 )
siendo:
( )καγ += M9
( 3.32 )
El parámetro de ajuste α determina la dispersión de los puntos sigma alrededor
del valor medio (10-3 en la mayoría de los casos), y suele ser un valor
ligeramente positivo. κ es un parámetro de escalado secundario que suele ser
0.
La segunda fase consiste en la predicción del vector de estado asociado a
cada punto sigma ( )1−kkqχ , a partir del modelo de maniobra. Como ya se
mencionó en 3.4.1.1.2, el sistema empleado para caracterizar el modelo de
maniobra es el modelo de aceleración de Singer.
La tercera fase predice el vector de estado y la matriz de covarianza mediante
la siguiente expresión:
( ) ( )
( ) ( ) ( )[ ] ( ) ( )[ ]
−−−−−−=−
−=−
∑
∑
=
=
Tq
M
q
qq
cX
M
q
qq
m
kkXkkkkXkkWkkP
kkWkkX
11111
11
18
0
18
0
χχ
χ
( 3.33 )
donde:
( )
≠
=−
=
0 si,
18
1
0 si, 19
2
2
q
M
qM
W
q
m
α
α
( 3.34 )
y:
( ) ( )
≠
=+−+−
=
0 si,
18
1
0 si,119
2
22
q
M
qM
W
q
c
α
βαα
( 3.35 )
El conocimiento a-priori de la distribución de X es incorporado en el parámetro
de ajuste β. Un valor de 2 es óptimo para las distribuciones Gaussianas.
3.4.4.2. Etapa de predicción (error sistemático)
Tal y como se adelantó en 3.4.3.2, el seguimiento de los errores sistemáticos
se efectúa mediante EKF. Debido a los cambios de las condiciones
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
83
atmosféricas y los errores de sincronismo, los errores sistemáticos varían con
el tiempo. En la presente tesis doctoral, se considera que el modelo dinámico
de error presentado en [79] es apropiado, siendo las ecuaciones de predicción
(vector de estado y matriz de covarianza respectivamente) las siguientes:
( ) ( ) ( ) ( )[ ]
( ) ( ) ( ) ( )[ ]
−−−+=−
−−−+=−
T
EEEE APkkPAPkkP
EkkEAEkkE
01101
01101
( 3.36 )
donde (k│k-1) corresponde a la predicción en la muestra k en base al instante
anterior k-1, y (k-1│k-1) corresponde al estadofiltrado en la muestra anterior.
El ajuste de los valores iniciales de E y PE - es decir, E(0) y PE(0) - viene
determinado por la evaluación de los errores sistemáticos off-line. En base a
los valores típicos de α1, α2, α3, y los errores de sincronismo obtenidos durante
los experimentos realizados en [17, 33], los valores de E(0) y PE(0) han sido
elegidos del siguiente modo:
• E(0) es un vector nulo, dado que no se conoce el valor de ninguno de los
parámetros a-priori
• PE(0) es una matriz diagonal, donde sus elementos tienen un valor dos
órdenes de magnitud mayor que el cuadrado de su parámetro asociado,
con el fin de que el sistema; si los valores son demasiado pequeños, el
filtro tardará en converger
Por su parte, la matriz de pesos A se emplea para ajustar la sensibilidad del
filtro en la etapa de predicción.
Tras esta etapa, se han obtenido las predicciones de los parámetros relativos a
los blancos (posición, velocidad y aceleración) y a los errores sistemáticos.
3.4.4.3. Etapa de filtrado (seguimiento del blanco)
La etapa de filtrado en el seguimiento del blanco (mediante UKF) se efectúa en
dos fases. La primera consiste en la estimación del TDOA reutilizando la
filosofía mostrada en la Figura 3.4. De este modo, el proceso de aumento
descrito en (3.31) es aplicado a X(k│k-1), obteniéndose χq(k│k-1). El modelo
de observación descrito a continuación es aplicado a χq(k│k-1), donde q
pertenece al conjunto de números enteros comprendidos entre 0 y 18M (ambos
inclusive):
( ) ( ) ( )[ ] ( ) ( )mimijmjijmjijmjijji nnTT
z
z
RRRR
c
RR
c
TDOA −+∆−∆+
−+−+−+−=
0
3
22
21 11
11
ααα
( 3.37 )
Denominamos el resultado de este proceso como ( )1−∆ kkqτ , con el fin de
enfatizar que el resultado obtenido es un vector de estimaciones de TDOA para
CAPÍTULO 3
84
cada punto sigma. Nótese que los valores asociados a este modelo de
observación están incluidos en ( )1−kkE . Finalmente, se efectúa la media
ponderada de ( )1−∆ kkqτ mediante la expresión (3.33), con el fin de obtener el
vector ( )1−∆ kkT , que contiene las estimaciones de TDOA.
La segunda fase consiste en el filtrado del vector de estado mediante la
incorporación de las medidas de TDOA, tal y como se define en la siguiente
expresión:
( ) ( ) ( )( )
( ) ( )
−−=
−∆−∆+−=
∆∆
T
TTXX KKPkkPkkP
kkTTKkkXkkX
,1
11
( 3.38 )
donde:
1
,,
−
∆∆∆= TTTX PPK
( ) ( )[ ] ( ) ( )[ ]111118
0
, −∆−−∆−∆−−∆= ∑
=
∆∆ kkTkkkkTkkWP
q
M
q
qq
cTT ττ
( ) ( )[ ] ( ) ( )[ ]111118
0
, −∆−−∆−−−= ∑
=
∆ kkTkkkkXkkWP
q
M
q
qq
cTX τχ
( 3.39 )
3.4.4.4. Etapa de filtrado (error sistemático)
Finalmente, se efectúa la etapa de filtrado (EKF) correspondiente al
seguimiento del error sistemático mediante el proceso descrito en [79]. Así, las
expresiones para dicha etapa son las siguientes:
( ) ( ) ( )
( ) ( ) ( ) ( )
−−−−≈
+−≈
111
1
kkVPkkPkkPkkP
vkkPkkEkkE
T
EEEE
E
( 3.40 )
donde el vector v es la suma por filas de la matriz vm. vm, a su vez, es obtenida
del siguiente modo:
[ ] ( )[ ]11 −∆−−= −∆ kkTTDOAFEFFFv lTlm ( 3.41 )
siendo:
( ) TssTsl FFFF
1−
=
[ ]T
SSSSSSSS
MSNSNMM
fffffffF
,1,,31,3,22,21,2
............=
( 3.42 )
SOLUCIÓN AL PROBLEMA DE MULTILATERACIÓN EN PRESENCIA DE
ERRORES DE PROPAGACIÓN Y SINCRONISMO
85
( ) ( ) ( )
∂
∂
∂
∂
∂
∂
=
j
j
i
j
j
i
j
j
i
S
z
TDOA
y
TDOA
x
TDOA
f
ji
,,
,
[ ]T
MSNSNMM
fffffffF
,1,,31,3,22,21,2
............ ∆∆∆∆∆∆∆∆ =
( ) ( ) ( ) ( )
( )
( )
( )
∆∂
∂
∆∂
∂
∂
∂
∂
∂
∂
∂
=∆
1,1,2321
,...,,,,
,
s
ji
N
j
i
j
i
j
i
j
i
j
i
T
TDOA
T
TDOATDOATDOATDOA
f
ααα
∆= FFF l
S
T
ll PRFFE +=
donde R es la matriz diagonal que contiene las covarianzas de los errores en la
medida de TDOA. PS es la sub-matriz compuesta por los elementos de ( )kkPX
asociados a las posiciones de los blancos. Recordamos al lector que ( )1−∆ kkT
es la estimación de los TDOA, basados en las posiciones de los blancos
estimadas, y las estimaciones de los errores sistemáticos.
Por otra parte, la matriz V se obtiene del siguiente modo:
[ ] [ ]FFFEFFFV l
T
l −−= ∆
−
∆
1
( 3.43 )
Las expresiones (3.40) a (3.43) estiman el vector de estado (E) cuya varianza
(PE) es mínima.
3.4.4.5. Cuestiones de implementación
Puesto que este filtro de seguimiento es un filtro recursivo, es necesario
emplear un método para obtener el paso inicial (arranque en frío). Siguiendo
las indicaciones de [21], esta tesis sugiere el empleo de la técnica descrita en
3.3.2 para obtener una serie de muestras de posición para cada aeronave.
Alcanzado un cierto número de muestras, se emplea una regresión cuadrática,
con el fin de estimar la posición, la velocidad y la aceleración de cada blanco.
Con esta información se confecciona el vector de estado de cada blanco a la
hora de inicializar el filtro de seguimiento.
Este capítulo muestra la solución al problema de seguimiento de blancos
mediante multilateración, donde las medidas están afectadas por la existencia
de errores sistemáticos (propagación y sincronismo) y errores aleatorios
asociados a la propia medida.
El siguiente capítulo aborda los experimentos realizados para evaluar el
sistema de seguimiento diseñado en este apartado.
86
CAPÍTULO 4. EXPERIMENTOS
4.1. Diseño de los experimentos
En primer lugar, es necesario comentar que si bien se ha intentado validar el
sistema descrito en el CAPÍTULO 3 mediante tráfico de oportunidad real,
únicamente se podía acceder a los datos correspondientes a la posición
estimada por los sistemas de vigilancia, y no a los datos de TDOA.
Así, la validación de dicho sistema consiste en experimentos empleando el
método de Monte Carlo en un escenario hipotético WAM, donde intervienen
una serie de aeronaves y un número de estaciones base. Dicha validación se
efectúa midiendo la exactitud y la precisión de las estimaciones de posición de
las aeronaves a lo largo del tiempo del experimento.
El diseño del protocolo experimental (Figura 4.1) está compuesto por tres
bloques:
• el bloque generador de TDOA, que está compuesto a su vez por:
o la actividad generadora del vector de tiempos, que construye el
vector con las muestras de tiempo empleadas en la simulación
o la actividad generadora de las trayectorias, que utiliza los
parámetros de las aeronaves (rumbo, velocidad, altitud, etc.), y
cuya salida son las posiciones reales de cada aeronave
o la actividad generadora de los tiempos de transmisión68
68 Se ha modelado el intervalo de transmisión para cada aeronave mediante una función
gaussiana cuya media es de 4 segundos y desviación típica de 1ms.
EXPERIMENTOS
87
o la interpolación de la posición de cada aeronave, con el fin de
obtener su posición en el momento de transmisión de la señal
o la actividad encargada de calcular los TOA de las señales (y por
ende, los TDOA), en función de la posición de cada aeronave en
su respectivo instante de transmisión, la posición de cada
estación base y los errores descritos en los apartados 2.3.1 y
2.3.2.
• el sistema de seguimiento, descrito en el apartado 3.4.4.
• los ajustes del experimento, que permiten al usuario configurar:
o las estaciones base, a través de su posición, el error de
sincronismo inicial y la deriva de los relojes
o las trayectorias de las aeronaves, mediante la posición inicial (x, y,
y altitud sobre el nivel del mar), la velocidad y el rumbo.
o el funcionamiento del sistema de seguimiento, especificando el
número de muestras antes de inicializar los filtros de seguimiento,
los parámetros del modelo de maniobra (Singer) y los parámetros
de los filtros UKF
Figura 4.1. Diagrama de bloques del protocolo experimental.
CAPÍTULO 4
88
4.2. Ajustes experimentales
El escenario de simulación empleadoen la presente Tesis Doctoral está
compuesto por 6 estaciones base y 8 aeronaves. Cada ensayo tiene una
duración de 15 minutos, donde cada aeronave transmite la señal WAM cada 4
segundos. El sistema de seguimiento emplea las medidas obtenidas cada 4
segundos para iterar los filtros de Kalman. Por otra parte, dicho sistema de
seguimiento estima la posición de cada aeronave tres de cada cuatro muestras
mediante extrapolación.
El experimento consta de 1000 ensayos. Esta sección enumera los ajustes de
las estaciones base, las aeronaves y el sistema de seguimiento.
4.2.1. Configuración de las estaciones base
4.2.1.1. Presentación de la configuración
En el apartado 2.3.1.2.2 se analizaron dos configuraciones de las estaciones
base, donde se concluyó que la geometría de la constelación C era más
favorable que la constelación H. La posición de las estaciones base - en base a
la configuración C - está definida en la Tabla 4.1 a continuación. La posición en
el eje z ha sido calculada empleando el modelo esférico de la Tierra, con radio
6371km.
Estación x (km) y (km) z (km) Altitud sobre el nivel del mar (km)
1 -50 50 -0.0392 0
2 0 50 0.0804 1
3 50 50 -0.0192 0.2
4 -50 -50 0.0108 0.5
5 0 -50 0.0204 0.4
6 50 -50 -0.0092 0.3
Tabla 4.1. Posición de las estaciones base.
Por otra parte, los ajustes empleados para modelar el error de sincronismo
inicial y la deriva de reloj para cada estación se encuentran en la Tabla 4.2.
Error de sincronismo inicial (s) Deriva de reloj (s/s)
N(10-10,10-10) N(3·10-11, 3·10-11)
Tabla 4.2. Parámetros característicos de los errores de sincronismo.
El efecto del ruido del receptor ha sido modelado mediante una v.a. gaussiana
de media nula y de desviación típica de 10ns. Esta aleatoriedad en el ruido del
receptor induce una aleatoriedad (ruido) en la medida de pseudodistancia de
3m. Los experimentos se efectúan mediante simulación, empleando el método
de Monte Carlo.
EXPERIMENTOS
89
4.2.1.2. Impacto de la configuración y expectativas
Este apartado muestra el impacto de la ubicación de las estaciones, a través
del concepto de DOP (introducido en 2.3.1.2.2).
La Figura 4.2 presenta la componente horizontal de dicha magnitud (HDOP),
mientras que la componente vertical (VDOP) se muestra en la Figura 4.3.
De estas figuras, se observa lo siguiente:
• el valor del HDOP es mayor fuera de la constelación, creciendo
rápidamente a lo largo del eje x, que es el eje de simetría de la
proyección de la constelación sobre el plano horizontal.
• el comportamiento del VDOP a lo largo del eje x es análogo al
comportamiento mostrado por el HDOP.
• lejos del eje x, el comportamiento del VDOP es tal que dicho valor es
mínimo cerca de cada estación base, siendo creciente con la distancia la
estación base más cercana.
• comparando los valores HDOP y VDOP, se observa que el efecto
multiplicativo del error es mayor en el eje vertical. De hecho, el VDOP
puede ser mayor que 100 en algunas zonas de la constelación. En
dichas zonas, la expectativa acerca de la desviación típica del error en la
estimación de altitud mediante localización hiperbólica es superior a
420m 69.
Figura 4.2. HDOP de la constelación empleada en los experimentos.
69 Se asume que el error en la estimación del tiempo de llegada puede modelarse mediante
una v.a. gaussiana de media nula y 10ns de desviación típica
CAPÍTULO 4
90
Figura 4.3. VDOP de la constelación empleada en los experimentos.
4.2.2. Análisis del intervalo de confianza de la media
Antes de continuar con la explicación de los parámetros experimentales
restantes, es necesario comprobar que las medidas de TDOA simuladas son
suficientemente insesgadas. Para ello, se efectúa el análisis del intervalo de
confianza de la media.
El efecto del ruido del receptor (en la medida de TOA) se ha caracterizado
mediante una variable aleatoria gaussiana de media nula y desviación típica
10ns. Así, dicho efecto en la medida de TDOA se caracteriza mediante una
variable aleatoria (también) gaussiana de media nula, con una desviación típica
de 14,142ns70.
Dada una variable aleatoria X, la relación entre la varianza de la media de una
v.a. gaussiana ( 2xσ ) y la varianza de una muestra (
2
xσ ) es:
N
x
x
2
2 σσ = ( 4.1 )
donde N es el número de ensayos.
70 Al tratarse de una diferencia de dos v.a. gaussianas de misma nula y misma varianza, la
función densidad de probabilidad resultante es una v.a. gaussiana de media nula y con una
desviación típica de 10√2 (ns).
EXPERIMENTOS
91
Partiendo de la base de que el ruido del receptor se modela mediante una v.a.
gaussiana de media nula (siendo así también para la medida de TDOA), el
intervalo de confianza de la media viene determinado por [82]:
+−
N
z
N
z xc
x
c
σσ
; ( 4.2 )
donde zc es un parámetro que depende del intervalo de confianza definido. La
relación entre zc y la probabilidad de que dicho valor esté comprendido entre
los dos extremos del intervalo (confianza) es:
( ) α−=≤≤− 1zZzP
( 4.3 )
( ) ( )
2
1
α
−=≤=Φ cc zZPz ( 4.4 )
donde Φ es la función de distribución normal.
De este modo, al definir un nivel de confianza matemática del 95% (α=0,05), el
valor de zc de 1,96.
( )( ) ( ) 96.1975.0
2
1 111 =Φ=
−Φ=ΦΦ= −−−
α
zz ( 4.5 )
siendo el intervalo de confianza de la media [-0.877ns; +0.877ns].
Así, tomando un determinado par de estaciones base en un instante
determinado del experimento, existe un 95% de probabilidad de que la media
de las medidas de TDOA, esté comprendida entre -0.877ns y +0.877ns.
El valor de 0.877ns es inferior a la décima parte de la desviación típica de la
deviación típica de una medida de TDOA. Por ello, esta cantidad de ensayos
permite afirmar que el efecto del ruido del receptor en la medida observada por
un par de estaciones en un momento determinado es insesgado.
4.2.3. Configuración del filtro de seguimiento
Las simulaciones contemplan el arranque en frío del sistema de seguimiento,
mediante el mecanismo descrito en 3.4.4.5. En los experimentos, la estimación
del vector de estado (del blanco) inicial se efectúa tras la estimación de la
novena muestra de posición mediante localización hiperbólica [21] (con el fin de
inicializar el filtro en la décima muestra).
CAPÍTULO 4
92
El modelo de maniobra empleado es el modelo de Singer [48]. Los parámetros
característicos empleados en este experimento se encuentran en la Tabla 4.3.
Aceleración máxima 19,62 m·s-2 (2G)
Constante de tiempo de la maniobra (s) 16
Probabilidad de aceleración máxima 0,1
Probabilidad de aceleración nula 0,6
Tabla 4.3. Parámetros del modelo de Singer.
Por otra parte, las constantes de ajuste del UKF han sido definidas siguiendo la
recomendación de [81], es decir: 10-3 para α, 2 para β y 0 para κ.
4.2.4. Configuraciones de las aeronaves
Dos escenarios diferentes han sido empleados para modelar las trayectorias de
las aeronaves: en el primer caso (A), varias aeronaves llegan a estar situadas
dentro de la constelación; por otra parte, en el segundo (B) escenario, todas las
aeronaves se encuentran fuera de la constelación en todo momento.
4.2.4.1. Escenario A
Las características del movimiento de cada aeronave están descritas en la
Tabla 4.5 a continuación. La trayectoria de cada aeronave es tal que la
velocidad sobre el suelo, el rumbo y la altitud son constantes. Debido a la
curvatura de la Tierra, esto no significa que el vector velocidad sea constante.
La Figura 4.4 a continuación ilustra el escenario definido por la posición de las
estaciones base (rojo) y las trayectorias de las aeronaves (negro).
Aeronave x inicial (km) y inicial (km) Altitud (km)
Velocidad sobre
el suelo (m/s)
Rumbo (º)
1 -150 60 10 250 90
2 100 20 9 200 270
3 80 -20 8 270 270
4 -120 -40 10 230 90
5 20 80 7 190270
6 -30 0 8 210 90
7 -30 40 8 200 90
8 -120 -60 9 260 90
Tabla 4.4. Trayectorias de las aeronaves en el escenario A.
EXPERIMENTOS
93
Figura 4.4. Escenario en base a los datos de la Tabla 4.1 y de la Tabla 4.4.
4.2.4.2. Escenario B
Del mismo modo, las características del movimiento de cada aeronave están
descritas en la Tabla 4.5 a continuación.
Aeronave x inicial (km) y inicial (km) Altitud (km)
Velocidad sobre
el suelo (m/s)
Rumbo (º)
1 -150 60 10 250 90
2 100 70 9 200 270
3 80 -80 8 270 270
4 -120 -70 10 230 90
5 20 80 7 190 270
6 -30 -90 8 210 90
7 -30 90 8 200 90
8 -120 -60 9 260 90
Tabla 4.5. Trayectorias de las aeronaves en el escenario B.
CAPÍTULO 4
94
La Figura 4.5 a continuación ilustra el escenario definido por la posición de las
estaciones base (rojo) y las trayectorias de las aeronaves (negro). Nótese que
ninguna aeronave entra en la constelación de las estaciones base. Tal como se
explicará en la discusión (apartado 4.4), este escenario ha sido elegido
deliberadamente con el fin de retar al sistema de seguimiento.
Figura 4.5. Escenario en base a los datos de la Tabla 4.1 y de la Tabla 4.5.
EXPERIMENTOS
95
4.3. Resultados experimentales
Las magnitudes medidas en los experimentos son el sesgo y la desviación
típica del error de la posición estimada (a una frecuencia de 1Hz) de las
aeronaves. El objetivo es comparar dichos valores empleando el sistema de
seguimiento, frente a los valores obtenidos mediante localización hiperbólica.
En un determinado instante de un determinado ensayo, el error en la
estimación de la posición de una aeronave determinada (j) es:
( )T
jjjjjjjjj zzyyxxrr −−−=−= ˆ,ˆ,ˆˆε
( 4.6 )
donde jx̂ , jŷ , jẑ son las componentes de la posición estimada, mientras que
jx , jy , jz son las componentes de la posición real.
Con el fin de obtener las componentes horizontal y vertical del error, es
necesario procesar los valores jε obtenidos mediante la ecuación (4.1).
Teniendo en cuenta las distancias empleadas en el estudio, el error horizontal
inducido por el error vertical debido a la curvatura de la Tierra es limitado. De
este modo, puede asumirse que:
( )
( )
−=
−−=
jjV
T
jjjjH
zz
yyxx
ˆ
ˆ,ˆ
ε
ε
( 4.7 )
Con el fin de evaluar el comportamiento del sistema a nivel global, en vez de
centrarnos en una aeronave en concreto, mostramos el mínimo (azul), la
mediana (negro) y el máximo (rojo) de las magnitudes medidas.
En los siguientes subapartados, se muestran dos gráficas por escenario: la
primera corresponde al sesgo del error de localización, mientras que la
segunda muestra la desviación típica. En cada figura, las líneas continuas
representan las magnitudes obtenidas mediante el sistema de seguimiento
presentado en esta tesis doctoral [17], mientras que las líneas discontinuas
representan los valores obtenidos mediante localización hiperbólica.
4.3.1. Escenario A
En base a la Figura 4.6, se observa que el sesgo obtenido mediante el método
presentado en esta tesis doctoral [17] es similar al obtenido mediante
localización hiperbólica. Por otra parte, el filtro de seguimiento contenido en el
sistema presentado reduce las oscilaciones en el sesgo, y disminuye la
desviación típica del error en un orden de magnitud (Figura 4.7).
CAPÍTULO 4
96
Figura 4.6. Sesgo de localización en el escenario A.
Figura 4.7. Desviación típica del error de localización en el escenario A.
La Figura 4.8 y la Figura 4.9 a continuación muestran respectivamente el
sesgo y la desviación típica de la componente horizontal del error de
localización obtenido mediante el sistema de seguimiento introducido en esta
tesis para cada aeronave, en función de su ubicación.
Se observa que la componente horizontal del sesgo de la posición crece a
medida que la aeronave se aleja de la constelación, especialmente a lo largo
del eje x. En cuanto a la exactitud de la medida, la desviación típica del error se
mantiene por debajo de los 30m para cada aeronave71, excepto en el caso de
la aeronave que sobrevuela el eje x.
71 A modo de referencia, el A320-200 (AIRBUS) mide 37,57m de largo y 34,1m de
envergadura.
EXPERIMENTOS
97
Figura 4.8. Media del error de localización horizontal (unidad: m).
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
10
20
30
40
50
60
70
80
Figura 4.9. Desviación típica del error de localización horizontal (unidad: m).
Por su parte, la Figura 4.10 y la Figura 4.11 muestran respectivamente el valor
absoluto del sesgo, y la desviación típica del error vertical de localización
(obtenidos mediante el sistema de seguimiento introducido en esta tesis) para
cada aeronave del escenario A, en función de su ubicación.
En la Figura 4.10 se observa que el sesgo es menor en los puntos cercanos a
cada estación base. Se observa además que el sesgo de la estimación de
altitud de la aeronave que sobrevuela el eje x es elevado a lo largo de su
trayectoria. Ambas observaciones coinciden con las características del
comportamiento del VDOP.
La Figura 4.11 muestra que la desviación típica del error de localización
vertical es significativamente inferior al valor esperado (mencionado en el
apartado 4.2.1.2).
CAPÍTULO 4
98
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
20
40
60
80
100
120
Figura 4.10. Valor absoluto de la media del error de localización vertical (unidad: m).
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
20
40
60
80
100
120
Figura 4.11. Desviación típica del error de localización vertical (unidad: m).
4.3.2. Escenario B
En la Figura 4.12, se observa que empleando el método presentado en esta
tesis doctoral, el sesgo del error se reduce hasta una cuarta parte del obtenido
mediante localización hiperbólica. Las mejoras significativas pueden
observarse en el caso peor (curva que representa el máximo del sesgo).
Al igual que en el escenario A, el efecto del filtro de seguimiento puede
observarse mediante la mejora (hasta en un orden de magnitud) de la exactitud
de las estimaciones de posición (Figura 4.13), limitando la desviación típica de
los errores de localización – tanto horizontal como vertical – a 100m a lo largo
del experimento.
EXPERIMENTOS
99
Figura 4.12. Sesgo de localización en el escenario B.
Figura 4.13. Desviación típica del error de localización en el escenario B.
De manera análoga al escenario A, la Figura 4.14 y la Figura 4.15 muestran
respectivamente el sesgo y la desviación típica de la componente horizontal del
error de localización obtenido mediante el sistema de seguimiento introducido
en esta tesis para cada aeronave, en función de su ubicación.
Se observa que el comportamiento de la componente horizontal del sesgo de la
posición es similar al caso del escenario A, siendo creciente a medida que la
aeronave se aleja de la constelación, especialmente a lo largo del eje x. En
cuanto a la exactitud de la medida, la desviación típica del error se mantiene
por debajo de los 30m para cada aeronave.
CAPÍTULO 4
100
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
10
20
30
40
50
60
70
80
Figura 4.14. Media del error de localización horizontal (unidad: m).
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.81
x 10
5
x (m)
y
(m
)
0
10
20
30
40
50
60
70
80
Figura 4.15. Desviación típica del error de localización horizontal (unidad: m).
Por su parte, la Figura 4.16 y la Figura 4.17 muestran respectivamente el valor
absoluto del sesgo, y la desviación típica del error vertical de localización
(obtenidos mediante el sistema de seguimiento introducido en esta tesis) para
cada aeronave del escenario B, en función de su ubicación.
En la Figura 4.16 se observa que el sesgo es menor en los puntos cercanos a
cada estación base. En cambio, el sesgo de la estimación de altitud para las
cuatro aeronaves cuyo valor absoluto en el eje y es de 80km o más, el sesgo
es elevado a lo largo de su trayectoria. Ambas observaciones coinciden con las
características del comportamiento del VDOP.
La Figura 4.17 muestra que la desviación típica del error de localización
vertical se encuentra dentro del valor esperado (mencionado en el apartado
4.2.1.2), excepto en el caso de la aeronave 1 (y=60km).
EXPERIMENTOS
101
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
20
40
60
80
100
120
Figura 4.16. Media del error de localización vertical (unidad: m).
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (m)
y
(m
)
0
20
40
60
80
100
120
Figura 4.17. Desviación típica del error de localización vertical (unidad: m).
De manera adicional, la Figura 4.18 y la Figura 4.19 muestran
respectivamente el cociente de la desviación típica de la componente horizontal
del error entre Hmin, y el cociente de la desviación típica de la componente
vertical del error entre Vmin. Los valores de Hmin y Vmin se han obtenido a partir
de la cota de Cramér-Rao, tal y como se efectuó en el apartado 2.3.1.2.472.
Se observa que ambos cocientes - calculados en cada punto de la trayectoria
de cada aeronave - son inferiores a los obtenidos mediante los algoritmos de
Bancroft y Chan-Ho (Figura 2.20, Figura 2.21, Figura 2.22 y Figura 2.23),
confirmándose que el estimador formado por el sistema de seguimiento es más
eficiente que los estimadores formados por los métodos de localización
hiperbólica de Bancroft (2.3.1.2.3.2) y Chan-Ho (2.3.1.2.3.3).
72 se asume que el estimador del sistema de seguimiento es un estimador insesgado.
CAPÍTULO 4
102
-1.5 -1 -0.5 0 0.5 1 1.5
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
σxy / Hmin
x (km)
y
(k
m
)
0
5
10
15
20
25
30
Figura 4.18. Cociente de la desviación típica de la componente horizontal del error entre Hmin.
-1.5 -1 -0.5 0 0.5 1 1.5 2
x 10
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x 10
5
x (km)
σz / Vmin
y
(k
m
)
0
5
10
15
20
25
30
Figura 4.19. Cociente de la desviación típica de la componente vertical del error entre Vmin.
4.4. Discusión
4.4.1. Acerca de las hipótesis
4.4.1.1. Modelo atmosférico
En estos experimentos se ha empleado la atmósfera de referencia para
modelar el efecto de la propagación troposférica. Este modelo no tiene en
cuenta las diferencias de presión y temperatura existentes en una determinada
altitud (como efecto de anticiclones y borrascas).
4.4.1.2. Probabilidad de detección
En la actualidad, los receptores de SSR y 1090ES garantizan una probabilidad
de detección del 98%. Partiendo de la base de que las aeronaves transmiten
EXPERIMENTOS
103
cada 4 segundos, esto significa que existe un 2% de probabilidad de que no se
reciba la señal durante 8 segundos, en los que una aeronave en ruta puede
recorrer hasta 1,2NM73. La separación horizontal mínima de 5NM reduce la
criticidad de la probabilidad de no detección, y por ello los experimentos han
asumido que todas las señales transmitidas por las aeronaves han sido
recibidas y procesadas correctamente.
4.4.1.3. Sistema de coordenadas
El sistema de coordenadas empleado en la navegación aérea está basado en
el elipsoide WGS84[83] para la latitud y la longitud, y la altitud de presión. Sin
embargo, el sistema de coordenadas de referencia empleado en la presente
tesis es ENU (Este, Norte, Arriba), basado en un modelo esférico de la Tierra, y
cuyo origen está situado en el centro de la constelación. Si bien este modelo
alternativo es suficiente para evaluar el sistema de seguimiento, es necesario
emplear el WGS84 si se pretende utilizar este sistema en entornos reales.
4.4.2. Acerca del cumplimiento de los requisitos
Los requisitos de precisión [84] referentes a la localización de aeronaves
mediante WAM son los siguientes:
• [R1 - obligatorio] El error de posicionamiento horizontal, incluyendo el
error en la medida y el error debido a la latencia de los datos, ha de ser
inferior o igual a 500m RMS globalmente, e inferior o igual a 550m RMS
por vuelo.
• [R2 - recomendado] El error de posicionamiento horizontal, incluyendo
el error en la medida y el error debido a la latencia de los datos, ha de
ser inferior o igual a 350m RMS globalmente, e inferior o igual a 385m
RMS por vuelo.
En base a los resultados mostrados en los apartados 4.3.1 y 4.3.2, puede
deducirse que:
• en el escenario A, ambos requisitos son cumplidos empleando ambos
métodos. Hay que recordar que las trayectorias de varias aeronaves se
encuentran en zonas con GDOP favorable.
• en el escenario B (GDOP desfavorable), el empleo de la localización
hiperbólica no garantiza el cumplimiento del requisito R1. Sin embargo,
el sistema de seguimiento presentado en esta tesis doctoral garantiza el
73 Se asume una velocidad sobre el suelo de 1000km/h, siendo 900km/h de TAS y 100km/h de
viento de cola.
CAPÍTULO 4
104
cumplimiento de R1 durante el tiempo del experimento, e incluso cumple
con R2 durante los 12 primeros minutos.
• el hecho de que el escenario B produzca peores resultados que en el
escenario A, concuerda con los valores de GDOP particularizados en los
puntos de las trayectorias de las aeronaves.
• el sistema de seguimiento presentado en esta tesis doctoral ofrece
mejores resultados que la localización hiperbólica, sobretodo en la
exactitud, gracias al filtro de seguimiento.
4.4.3. Otras consideraciones
En el caso de que la ICAO decidiese dar el salto para cambiar la altitud de
navegación de la altitud de presión a la altitud elipsoidal (poco factible a corto-
medio plazo), sería imperativo reducir el sesgo en la estimación de la altitud, si
se pretende que WAM mejore la integridad de la altitud elipsoidal (ofrecida por
los GNSS). Esto requeriría un refinamiento del modelo de la estimación de
pseudodistancia en función de la altitud de la estación base.
105
CAPÍTULO 5. CONCLUSIONES Y
TRABAJOS FUTUROS
5.1. Conclusiones
En el capítulo 1 se identifica WAM como uno de los constituyentes del concepto
CNS/ATM, cuyo objetivo es dar soporte a un sistema global de gestión de
tráfico aéreo. La naturaleza de esta técnica de localización permite su empleo
aun en caso de fallo o inhibición de los sistemas de posicionamiento mediante
satélite (GPS, Galileo, GLONASS, etc.).
El capítulo 2 analiza el funcionamiento básico del sistema WAM, estudiando los
errores – tanto aleatorios (ruido de receptor) como sistemáticos (propagación
troposférica y sincronismo) – asociados a las medidas de TDOA.
Con relación a los errores aleatorios, esta tesis demuestra que para una misma
caracterización del ruido de receptor, la precisión de la posición del blanco
depende de la geometría de la constelación WAM. Este aspecto es capital para
los ANSP, a la hora de planificar una instalación WAM.
En cuanto a los errores sistemáticos, esta tesis doctoralaporta un nuevo
modelo de sesgo no lineal, en función de la distancia euclídea, y de la altitud de
la aeronave sobre el nivel del mar. Este capítulo se cierra mediante el
planteamiento del problema matemático.
El capítulo 3 analiza el sistema de ecuaciones mencionado anteriormente, y
presenta la solución a los dos problemas de multilateración en presencia de los
CAPÍTULO 5
106
errores identificados anteriormente: la localización muestra a muestra, y el
seguimiento.
En cuanto a la localización muestra a muestra, esta tesis ha analizado el uso
de un método basado en dos pasos, donde el primero consiste en emplear un
método de localización hiperbólica con el fin de estimar una solución inicial,
mientras que el segundo refina dicha solución mediante un algoritmo. Además
de resaltar el problema de convergencia asociado a cualquier algoritmo
iterativo, esta tesis enuncia que este enfoque no permite obtener una serie de
soluciones que optimicen las estimaciones (tanto de posición, como de sesgo)
para el tiempo total de experimento.
Esta tesis aporta un nuevo sistema que efectúa el seguimiento conjunto de
blancos y de errores sistemáticos para WAM, empleando dos filtros de Kalman
entrelazados: uno para el seguimiento de blancos, y otro para el seguimiento
de los errores sistemáticos.
Por otra parte, hay que destacar que debido a la naturaleza de la localización
hiperbólica, la función que caracteriza la medida de TDOA (sin considerar los
errores sistemáticos) no es lineal. Por ello, esta tesis aporta además el empleo
del Unscented Kalman Filter para el seguimiento de los blancos.
En el capítulo 4, esta tesis presenta un entorno de simulación escalable que
permite gestionar las trayectorias de las aeronaves, las posiciones de las
estaciones base, y los principales parámetros del sistema de seguimiento.
Los experimentos mediante simulación muestran que este sistema ofrece
mejoras en la precisión y la exactitud de la posición estimada, respecto de la
localización hiperbólica. Además, dichos resultados experimentales muestran
que el sistema de seguimiento presentado en esta tesis cumple con el requisito
obligatorio de localización horizontal para garantizar la navegación en ruta en
Europa, incluso en escenarios adversos.
La mejora en el eje vertical es notable, teniendo en cuenta los valores de
VDOP asociados a los puntos de las trayectorias de cada aeronave.74
Así, se puede afirmar que el nuevo sistema de seguimiento que aporta esta
tesis doctoral permite cumplir con los requisitos de navegación en ruta en
Europa, en caso de fallo temporal de los GNSS.
74 El requisito obligatorio de localización vertical para garantizar la navegación en ruta en
Europa está basado en la medidas de altitud barométrica, con lo que no forma parte del tema
de la presente tesis doctoral.
CONCLUSIONES Y TRABAJOS FUTUROS
107
5.2. Trabajos futuros
5.2.1. Mejoras del sistema de seguimiento
La arquitectura del protocolo experimental está pensada para poder modificar
de manera sencilla las condiciones de ensayo, pudiendo alterar las trayectorias
de las aeronaves, los errores de sincronismo y el modelo de maniobra
empleado por el filtro de seguimiento.
Al mismo tiempo, la separación en tres bloques bien diferenciados permite que
futuras mejoras del sistema de seguimiento puedan ser empleadas reutilizando
la arquitectura existente. De este modo, las líneas de investigación propuestas
en esta tesis son:
• la mejora del modelo de estimación de la pseudodistancia en función de
la altitud de la estación base.
• el empleo de medidas reales de TDOA.
• el empleo de modelos de maniobra más complejos a la hora de simular
las trayectorias de las aeronaves.
• el empleo de IMM para detectar y aplicar el tipo de maniobra más
apropiado (cambios de rumbo, ascensos, descensos, etc.); el IMM actúa
de “mezclador” entre una serie de filtros de seguimiento con la
arquitectura del sistema presentado en esta tesis, pero con distintos
modelos de dinámica del blanco.
La arquitectura del protocolo experimental está pensada para poder modificar
de manera sencilla las condiciones de ensayo, pudiendo alterar las trayectorias
de las aeronaves, los errores de sincronismo y el modelo de maniobra
empleado por el filtro de seguimiento.
5.2.2. Siguiente fase de los experimentos
Con el fin de dar el salto de la simulación al mundo real, es necesario emplear
datos de medidas reales de TDOA para comprobar las prestaciones del
sistema de seguimiento, y compararlas con los sistemas ofrecidos por la
Industria.
5.2.3. Otras líneas de investigación
El sistema de seguimiento está basado en el empleo de UKF para el
seguimiento de blancos, y el EKF para el seguimiento de los parámetros
empleados en el modelo de estimación de pseudodistancia. Esta tesis sugiere
también el empleo de otras filosofías para la localización, basados en los filtros
de partículas [85], o en algoritmos iterativos aplicados al ajuste en secciones
cónicas [86].
108
CAPÍTULO 6. BIBLIOGRAFÍA
[1] Comisión Europea, "Reglamento (CE) Nº 549/2004 del Parlamento Europeo
y del Consejo de 10 de marzo de 2004, por el que se fija el marco para la
creación del cielo único europeo." 2004.
[2] Comisión Europea, "Reglamento (CE) Nº 550/2004 del Parlamento Europeo
y del Consejo de 10 de marzo de 2004, relativo a la prestación de servicios de
navegación aérea en el cielo único europeo," 2004.
[3] Comisión Europea, "Reglamento (CE) Nº 551/2004 del Parlamento Europeo
y del Consejo de 10 de marzo de 2004, relativo a la organización y utilización
del espacio aéreo en el cielo único europeo," 2004.
[4] Comisión Europea, "Reglamento (CE) Nº 552/2004 del Parlamento Europeo
y del Consejo de 10 de marzo de 2004, relativo a la interoperabilidad de la red
europea de gestión del tránsito aéreo," 2004.
[5] Comisión Europea, "Reglamento (CE) No 219/2007 del Consejo de 27 de
febrero de 2007, relativo a la constitución de una empresa común para la
realización del sistema europeo de nueva generación para la gestión del
tránsito aéreo (SESAR)," 2007.
[6] Comisión Europea, "Reglamento (CE) No 1361/2008 del consejo de 16 de
diciembre de 2008 por el que se modifica el Reglamento (CE) no 219/2007,
relativo a la constitución de una empresa común para la realización del sistema
europeo de nueva generación para la gestión del tránsito aéreo (SESAR),"
2008.
[7] Comisión Europea, "Reglamento (CE) nº 1070/2009 del Parlamento
Europeo y del Consejo, de 21 de octubre de 2009, por el que se modifican los
Reglamentos (CE) nº 549/2004, (CE) nº 550/2004, (CE) nº 551/2004 y (CE) nº
552/2004 con el fin de mejorar el rendimiento y la sostenibilidad del sistema
europeo de aviación." 2009.
BIBLIOGRAFÍA
109
[8] Comisión Europea, "Reglamento (CE) Nº 2150/2005 de la Comisión de 23
de diciembre de 2005, por el que se establecen normas comunes para la
utilización flexible del espacio aéreo," 2005.
[9] DGAC and INAC, "SW FAB - Demonstration of Compliance with
Commission Regulation (EU) No 176/2011," 21/06/2012, 2012.
[10] ICAO, Global Air Navigation Plan for CNS/ATM Systems, Doc 9750 -
AN/963. 2002.
[11] ICAO, "Secondary Surveillance Radar - Mode S. Advisory Circular 174-
AN/110," 1983.
[12] ICAO, "Technical Provisions for Mode S Services and Extended
Squitter, Doc 9871," 2012.
[13] NLR, "Wide Area Multilateration - Report on EATMP TRS 131/04," 2005.
[14] (May 16, 2011). Austro Control Selects Sensis Wide Area Multilateration
For Austria’s Nationwide Surveillance Network. Available:
http://www.saabsensis.com/austro-control-selects-sensis-wide-area-
multilateration-for-austrias-nationwide-surveillance-network/.
[15] Sensis Multilateration The First FAA-certified Multilateration System For
Separation Of En Route Aircraft. Available: http://www.saabsensis.com/sensis-
multilateration-the-first-faa-certified-multilateration-system-for-separation-of-en-route-aircraft/.
[16] C. Daskalakis and P. Martone, "Alternative surveillance technology for the
gulf of mexico," in Digital Avionics Systems Conference, 2004. DASC 04. the
23rd, Salt Lake City, UT, U.S.A., 2004, pp. 1.D.4-1.1-8.
[17] J. Abbud and G. De Miguel, "Joint target tracking and systematic error
correction for wide area multilateration," in Proceedings of the ENRI
International Workshop on ATM/CNS (EIWAC2013), Tokyo, Japan, 2013, .
[18] S. Atkinson and C. e. a. Heyes, "GPS synchronisation of WAM systems -
pros and cons," in Enhanced Surveillance of Aircraft and Vehicles (ESAVS
2010), Berlin, Germany, 2010, .
[19] EUROCONTROL, EUROCONTROL STANDARD DOCUMENT FOR
SURVEILLANCE DATA EXCHANGE, Part 12 : Category 021, ADS-B Reports.
2011.
[20] M. A. Lombardi, "The use of GPS disciplined oscillators as primary
frequency standards for calibration laboratories," in 2008 NCSL International
Workshop and Symposium, Orlando, FL (USA), 2008, .
CAPÍTULO 6
110
[21] Y. T. Chan and K. C. Ho, "A simple and efficient estimator for hyperbolic
location," IEEE Transactions on Signal Processing, vol. 42, pp. 1905-1915,
1994.
[22] E. Kaplan and C. Hegarty, Understanding GPS: Principles and Applications.
Artech House, 2005.
[23] D. J. Torrieri, "Statistical Theory of Passive Location Systems," IEEE
Transactions on Aerospace and Electronic Systems, vol. 20, pp. 183-198, 1984.
[24] G. Galati et al, "Wide Area Surveillance using SSR Mode S Multilateration:
advantages and limitations," European Radar Conference (EURAD), 2005.
[25] W. H. Foy, "Position-Location Solutions by Taylor-Series Estimation," IEEE
Transactions on Aerospace and Electronic Systems, vol. AES-12, pp. 187-194,
March 1976, 1976.
[26] S. Bancroft, "An algebraic solution of the GPS equations," IEEE
Transactions on Aerospace and Electronic Systems, vol. AES-21, pp. 56-59,
1985.
[27] M. Deffenbaugh, J. G. Bellingham and H. Schmidt, "The relationship
between spherical and hyperbolic positioning," in OCEANS '96. MTS/IEEE.
Prospects for the 21st Century. Conference Proceedings (Volume:2). Fort
Lauderdale, FL, 1996, pp. 590-595.
[28] Unión Internacional de Telecomunicaciones, "Recomendación UIT-R
P.453-6, Índice de refracción radioeléctrica: su fórmula y datos sobre la
refractividad." 1997.
[29] L. V. Blake, Radar Range-Performance Analysis. 1986.
[30] Unión Internacional de Telecomunicaciones, "Recomendación UIT-R
P.834-6, Efectos de la refracción troposférica sobre la propagación de las
ondas radioeléctricas," 2007.
[31] D. K. Barton, Radar System Analysis and Modeling. 2004.
[32] G. De Miguel, J. Besada and J. Garcia, "Correction of propagation errors in
Wide Area Multilateration systems," Proceedings of the 6th European Radar
Conference, pp. 81-84, 2009.
[33] J. Abbud, G. De Miguel and J. Besada. Correction of systematic errors in
wide area multilateration. Presented at 2011 Tyrrhenian International Workshop
on Digital Communications – Enhanced Surveillance of Aircraft and Vehicles.
2011, .
[34] OVEN CONTROLLED CRYSTAL OSCILLATOR, AOCJY5 Series,
Datasheet. Available: http://www.abracon.com/Precisiontiming/AOCJY5.pdf.
BIBLIOGRAFÍA
111
[35] Y. S. Shmaliy, L. Arceo-Miquel, J. Munoz-Diaz and O. Ibarra-Manzano,
"Unbiased FIR estimates vs the sawtooth-corrected GPS-based measurements:
experimental evaluaion," 38th Annual Precise Time and Time Interval (PTTI)
Meeting, pp. 347-360, 2007.
[36] R. E. Kalman, "A New Approach to Linear Filtering and Prediction
Problems," Transactions of the ASME - Journal of Basic Engineering, vol. 82,
pp. 35-45, 1960.
[37] Y. Bar-Shalom, X. R. Li and T. Kirubarajan, Estimation with Applications to
Tracking and Navigation: Theory, Algorithms, and Software. 2001.
[38] X. Rong Li and V. P. Jilkov, "Survey of Maneuvering Target Tracking. Part
I: Dynamic Models," IEEE Transactions on Aerospace and Electronic Systems,
vol. 39, pp. 1333-1364, 2003.
[39] X. Rong Li and V. P. Jilkov, "Survey of Maneuvering Target Tracking - Part
IV: Decision-based methods." Proceedings of the 2002 SPIE Conference on
Signal and Data Processing of Small Targets, vol. 4728, April 2002, 2002.
[40] Y. T. Chan, A. G. C. Hu and J. B. Plant, "A Kalman filter based tracking
scheme with input estimation." IEEE Transactions on Aerospace and Electronic
Systems, vol. 15, pp. 237-244, Mar. 1979, 1979.
[41] J. Korn, S. W. Gully and A. S. Willsky, "Application of the generalized
likelihood ratio algorithm to maneuver detection and estimation." Proceedings of
the 1982 American Control Conference, Arlington, VA, 1982.
[42] P. L. Bogler, "Tracking a maneuvering target using input estimation." IEEE
Transactions on Aerospace and Electronic Systems, vol. 23, pp. 298-310, May
1987, 1987.
[43] P. L. Bogler, Radar Principles with Applications to Tracking Systems. Wiley,
1990.
[44] M. Farooq and S. Bruder, "Information type filters for tracking a
maneuvering target." IEEE Transactions on Aerospace and Electronic Systems,
vol. 26, pp. 441-454, 1990.
[45] I. H. Whang, J. Lee and T. Sung, "Modified input estimation technique
using pseudoresiduals." IEEE Transactions on Aerospace and Electronic
Systems, vol. 30, pp. 220-228, 1994.
[46] T. K. Sung and J. G. Lee, "A decoupled adaptive tracking filter for real
applications." IEEE Transactions on Aerospace and Electronic Systems, vol. 33,
pp. 1025-1030, 1997.
[47] H. Lee and M. Tahk, "Generalized input-estimation technique for tracking
maneuvering targets." IEEE Transactions on Aerospace and Electronic
Systems, vol. 35, pp. 1388-1402, 1999.
CAPÍTULO 6
112
[48] R. A. Singer, "Estimating optimal tracking filter performance for manned
maneuvering targets," IEEE Transactions on Aerospace and Electronic
Systems, vol. 6, pp. 473-483, July 1970, 1970.
[49] R. J. McAulay and E. J. Denlinger, "A decision-directed adaptive tracker."
IEEE Transactions on Aerospace and Electronic Systems, vol. 9, pp. 229-236,
Mar. 1973, 1973.
[50] J. R. Cloutier, J. H. Evers and J. J. Feeley, "Assessment of air-to-air missile
guidance and control technology." IEEE Control Systems Magazine, vol. 9, pp.
27-34, Oct. 1989, 1989.
[51] P. F. Easthope and N. W. Heys, "Multiple-model target-oriented tracking
system." Proceedings of the 1994 SPIE Conference on Signal and Data
Processing of Small Targets, vol. 2235, Apr. 1994, 1994.
[52] P. J. Costa, "Adaptive model architecture and extended Kalman-Bucy
filters." IEEE Transactions on Aerospace and Electronic Systems, vol. 30, pp.
525-533, Apr. 1994, 1994.
[53] C. N. D’Souza, M. A. McClure and J. R. Cloutier, "Spherical target state
estimators," Proceedings of 1994 American Control Conference, Baltimore, MD,
pp. 1675-1679, June 1994, 1994.
[54] G. Pulford and B. La Scala, "A survey of manoeuvring target tracking
methods and their applicability to over-the-horizon radar." Cooperative
Research Center for Sensor Signal and Information Processing, vol. Technical
Report CSSIP 14, July 1996, 1996.
[55] S. S. Blackman and R. F. Popoli, Design and Analysis of Modern Tracking
Systems. Boston, MA: Artech House, 1999.
[56] R. A. Howard, "System analysis of semi-Markov processes." IEEE
Transactions Military Electronics, vol. 8, pp. 114-124, Apr. 1964, 1964.
[57] R. L. Moose and P. L. Wang, "An adaptive estimator with learning for a
plant containing semi-Markov switching parameters." IEEE Transactions on
Systems, Man, Cybernetics, vol. 3, pp. 277-281, May 1973, 1973.
[58] R. L. Moose, "An adaptive state estimator solution to the maneuvering
target problem." IEEE Transactions on Automatic Control, vol. 20, pp. 359-362,
June 1975, 1975.
[59] N. H. Gholson and R. L. Moose, "Maneuvering target tracking using
adaptive state estimation." IEEE Transactions on Aerospace and Electronic
Systems, vol. 13, pp. 310-317, May 1977, 1977.
[60] R. L. Moose, H. F. VanLandingham and D. H. McCabe, "Modeling and
estimation of trackingmaneuvering targets." IEEE Transactions on Aerospace
and Electronic Systems, vol. 15, pp. 448-456, May 1979, 1979.
BIBLIOGRAFÍA
113
[61] X. Rong Li, "Model-set design for multiple-model estimation—Part I." in
2002 International Conference on Information Fusion, Annapolis, MD, 2002, pp.
26-33.
[62] S. S. Lim and M. Farooq, "Maneuvering target tracking using jump
processes." in Proceedings of the 30th IEEE Conference on Decision and
Control, Brighton, UK, 1991, pp. 2049-54.
[63] L. Campo, P. Mookerjee and Y. Bar-Shalom, "State estimation for systems
with sojourn-time-dependent Markov model switching." IEEE Transactions on
Automatic Control, vol. 36, pp. 238-243, Feb. 1991, 1991.
[64] D. D. Sworder, M. Kent, R. Vojak and R. G. Hutchins, "Renewal models for
maneuvering targets." IEEE Transactions on Aerospace and Electronic
Systems, vol. 31, pp. 138-149, Jan. 1995, 1995.
[65] D. R. Cox, Renewal Theory. London: Methuen, 1962.
[66] D. D. Sworder, P. F. Singer and R. G. Hutchins, "Image-enhanced
estimation methods." Proceedings of the IEEE, vol. 81, pp. 797-812, June 1993,
1993.
[67] X. Rong Li and V. P. Jilkov, "Survey of maneuvering target tracking - part
VI: Approximation techniques for nonlinear filtering," in Proceedings of 2004
SPIE Conference on Signal and Data Processing of Small Targets, San Diego,
CA, USA, 2004, pp. 537-550.
[68] S. J. Julier and J. K. Uhlmann, "A new approach for filtering nonlinear
sytems," in Proceedings of the 1995 American Control Conference, Seattle, WA,
USA, 1995, pp. 1628-–1632.
[69] S. J. Julier and J. K. Uhlmann, "A new extension of the kalman filter to
nonlinear systems," in Proceedings of AeroSense: The 11th International
Symposium on Aerospace/Defence Sensing, Simulation and Controls, Session:
Multi Sensor Fusion and Resource Management II, Orlando, FL, USA, 1997, .
[70] S. J. Julier, J. K. Uhlmann and H. F. Durrant-Whyte, "A New Method for
Nonlinear Transformation of Means and Covariances in Filters and Estimators,"
IEEE Transactions on Automatic Control, vol. 45, pp. 477-482, Mar. 2000, 2000.
[71] J. R. Van Zandt, "Boost phase tracking with an unscented filter," in 2002
SPIE Conf. on Signal and Data Processing of Small Targets , Orlando, FL,
USA, 2002, pp. 263-274.
[72] A. Farina, B. Ristic and D. Benvenuti, "Tracking a Ballistic Target:
Comparison of Several Filters," IEEE Transactions on Aerospace and Electronic
Systems, vol. 38, pp. 1916-1924, 2002.
[73] B. Ristic, A. Farina, D. Benvenuti and M. Arulampalam, "Performance
Bounds and Comparison of Nonlinear Filters for Tracking a Ballistic Object on
CAPÍTULO 6
114
Re-entry. ," IEE Proceedings Radar Sonar Navigation, vol. 150, pp. 65-70, April
2003, 2003.
[74] E. Kraft, "A quaternion-based unscented kalman filter for orientation
tracking." in 2003 International Conference on Information Fusion, Cairns,
Australia, 2003, pp. 42-54.
[75] T. M. Nguyen, V. P. Jilkov and X. Rong Li, "Comparison of sampling-based
algorithms for multisensor distributed target tracking," in 2003 International
Conference on Information Fusion, Cairns, Australia, 2003, pp. 114-121.
[76] M. Briers, S. Maskel and R. Wright, "A rao-blackwellised unscented kalman
filter," in 2003 International Conference on Information Fusion, Cairns, Australia,
2003, pp. 55-61.
[77] B. Ristic, M. Arulampalam and N. Gordon, Beyond the Kalman Filter.
Particle Filters for Tracking Applications. Artech House, 2004.
[78] E. A. Wan and R. Van der Merwe, "The unscented kalman filter," in Kalman
Filtering and Neural Networks, S. Haykin, Ed. John Wiley & Sons, Inc., 2001, pp.
221-280.
[79] Y. Bar-Shalom, "Systematic error estimation," in Multitarget-Multisensor
Tracking: Applications and Advances, Apartado 2.4.Anonymous 1992, pp. 38-
41.
[80] Y. Bar-Shalom, "Joint tracking and sensors' systematic error estimation," in
Multitarget-Multisensor Tracking: Applications and Advances, Apartado 2.B.
Anonymous 1992, pp. 56-63.
[81] S. J. Julier and J. K. Uhlmann, "Unscented Filtering and Nonlinear
Estimation," Proceedings of the IEEE, vol. 92, 2004.
[82] A. M. Law and W. D. Kelton, Simulation Modeling and Analysis. New York:
McGraw-Hill Education, 1991.
[83] US Department of Defense, Department of Defense World Geodetic
System 1984, its Definition and Relationships with Local Geodetic Systems
. 1997.
[84] EUROCONTROL, "EUROCONTROL Specification for ATM Surveillance
System Performance," vol. 1, 2012.
[85] A. Doucet and A. M. Johansen, "A Tutorial on Particle Filtering and
Smoothing: Fifteen years later," 2008.
[86] Z. Zhang, "Parameter Estimation Techniques: A Tutorial with Application to
Conic Fitting," Image and Vision Computing, vol. 15, pp. 59-76, 1997.
BIBLIOGRAFÍA
115
[87] National Imagery and Mapping Agency (NIMA), Department of Defense
World Geodetic System 1984, its Definition and Relationships with Local
Geodetic Systems. 2000.
[88] (May 06, 2013). Earth Gravitational Model 2008 (EGM2008). Available:
http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/.
[89] G. e. a. Cai, "Coordinate systems and transformations," in Unmanned
Rotorcraft Systems, Advances in Industrial Control, London: Springer-Verlag,
2011, pp. 23-34.
[90] B. R. Bowring, "Transformation from spatial to geographical coordinates,"
vol. 23, pp. 323-327, 1976.
[91] EUROCONTROL. EUROCONTROL specification for the application of the
flexible use of airspace (FUA). [General Public]. (1.1), 2009.
[92] EUROCONTROL, "Introduction to the Mission Trajectory," 25/05/2010,
2010.
116
ANEXO I. SISTEMAS DE COORDENADAS Y
TRANSFORMACIONES
I.1. Sistema geodésico mundial - WGS
El World Geodetic System (WGS), es un sistema comúnmente empleado en
cartografía, geodesia y navegación. Está compuesto por un sistema de
coordenadas para la Tierra, una superficie elipsoidal de referencia como
aproximación de la altitud, y una superficie gravitacional equipotencial (geoide)
que define el nivel nominal del mar.
La última revisión del sistema WGS es WGS 84, establecida en 1984 y
revisada en 2004. Sus precursores son WGS 60, WGS 66 y WGS 72. WGS 84
es el sistema de coordenadas de referencia empleado en GPS.
I.1.1. Elipsoide WGS84
El elipsoide WGS 84 [87] es un elipsoide de revolución definido por los
siguientes parámetros:
• semieje mayor (a): 6378137,0 m
• inversa del achatamiento75 (f): 298,257223563
75 El achatamiento (o elipticidad) es la medida de compresión de un círculo o una esfera a lo
largo de su diámetro para formar una elipse o un elipsoide de revolución (esferoide). La
notación usual para el achatamiento es f y su definición en términos de los semiejes de la
elipse o elipsoide resultante es: (a-b)/a, donde a es la longitud del semieje mayor, y b la
longitud del semieje menor
SISTEMAS DE COORDENADAS Y TRANSFORMACIONES
117
• Producto de la Constante Gravitacional (G) y la Masa de la Tierra (M):
3,986004418x1014 m3/s2
• Velocidad Angular de la Tierra (ω): 7,292115x10-5 rad/s
I.1.2. Sistemas de coordenadas
I.1.2.1 Sistema de coordenadas ECEF
El sistema de coordenadas ECEF 76 es un sistema cartesiano ortonormal,
ilustrado en la Figura I.1., cuyas propiedades se definen a continuación:
• origen: centro de masas de la Tierra
• plano fundamental: plano ecuatorial medio correspondiente a la eclíptica
del año 2000
• eje Z: dirección del polo medio convencional terrestre definido por el
IERS 77 , perpendicular al plano fundamental. También es el eje de
rotación de esta elipse de revolución
• eje X: intersección del meridiano de referencia IERS con la superficie
normal al eje Z que contiene al origen
• eje Y: es el producto vectorial del vector unitario en la dirección y sentido
del eje Z y del vector unitario en la dirección y sentido del eje X,
completando así el sistema de coordenadas ECEF
Figura I.1.Definición del sistema de coordenadas WGS 84.
76 Earth Centered Earth Fixed
77 International Earth Rotation and Reference Systems Service: http://www.iers.org/
ANEXO I
118
I.1.2.2 Sistema de coordenadas geodésicas
El sistema de coordenadas geodésicas es un sistema de coordenadas en que
la posición viene determinada por la 3-tupla de los siguientes elementos:
• la latitud geodésica φ, formada por el ángulo que forma el plano
ecuatorial con la perpendicular al elipsoide desde un punto dado. Se
toma positiva hacia el norte.
• la longitud geodésica λ, definida por el ángulo que forma el plano
meridiano de un punto dado. El meridiano de longitud cero es el
meridiano de referencia IERS, situado a 5,31’’ al este del meridiano de
Greenwich a la altitud del Royal Observatory de Greenwich (102,5
metros). Se toma positiva hacia el Este.
• la altitud geodésica (o elipsoidal) h, que es la distancia a un punto desde
el elipsoide medida a lo largo de la normal al elipsoide por este punto.
Este valor es positivo si el punto se encuentra fuera del elipsoide.
Figura I.2. Sistema de coordenadas geodésicas.
SISTEMAS DE COORDENADAS Y TRANSFORMACIONES
119
I.1.2.3 Sistema de coordenadas local ENU
Este sistema de coordenadas está centrado en un punto de referencia P, de
coordenadas geodésicas {φ0, λ0, h0}, y sus ejes vienen definidos a
continuación:
• U (arriba, del inglés up), cuya dirección viene definida por la recta SP, y
su sentido es el que va de S hacia P
• N (norte, del inglés North), cuya dirección y sentido vienen definidos por
la proyección del meridiano en el punto P al plano ortogonal al eje U.
• E (este, del inglés East), definido por el producto vectorial NxU
La superficie definida por los ejes E y N se conoce como plano tangente local.
Figura I.3. El punto P es el centro del sistema de referencia local ENU.
ANEXO I
120
I.1.3. Relación entre el nivel del mar y el elipsoide
WGS84
I.1.3.1 Geoide
El geoide es la superficie equipotencial que coincidiría con la superficie media
de la Tierra si los océanos y atmósferas estuvieran en equilibrio (relajación) con
relación a la Tierra, y extendida a través de los continentes, ilustrado en la
Figura I.4. Por ello, este elemento es considerado habitualmente como la
verdadera forma física de la Tierra, y no la figura regular definida por el
elipsoide de referencia.
La superficie del geoide se encuentra por fuera del elipsoide de referencia
cuando hay una anomalía positiva en la gravedad (debido a un exceso de
masa), y por dentro en caso de ser la anomalía negativa (déficit de masa). Las
diferencias en la gravedad surgen de la distribución heterogénea de la masa en
la Tierra.
Figura I.4. Ilustración de los distintos conceptos.
1. Océano; 2. Elipsoide de referencia; 3. Línea de la plomada local; 4. Continente; 5. Geoide
I.1.3.2 Altura ortométrica y altura elipsoidal
La altura ortométrica (H) de un punto P es la distancia de P al geoide, que
comparte la misma latitud y longitud que P. Es lo que se conoce comúnmente
como altitud sobre el nivel del mar. Por otra parte, la altitud elipsoidal (h) del
punto P es la distancia de ese punto al elipsoide de referencia.
Figura I.5. Relación entre la altitud ortométrica y la altitud elipsoidal.
SISTEMAS DE COORDENADAS Y TRANSFORMACIONES
121
En la práctica, se considera que la altitud ortométrica y la altitud elipsoidal son
colineales, con lo que:
NhH −≈ ( I.1 )
donde N es la altura del geoide respecto de la elipse de referencia.
El modelo gravitacional terrestre EGM200878 es un modelo geopotencial de la
Tierra, completo hasta el harmónico esférico de grado y orden 2159, que
contiene coeficientes adicionales hasta grado 2190 y orden 2159 [88].
La Figura I.6 a continuación muestra la altitud del geoide (valor N en la figura
Figura I.5), obtenida calculando los coeficientes harmónicos esféricos hasta
grado y orden 500.
Figura I.6. Altitud del geoide EGM2008 respecto del elipsoide WGS84.
De este modo, pueden obtenerse las coordenadas geodésicas de un punto
definido por su latitud, longitud y altura sobre el nivel del mar. En la sección I.2
se incluyen las transformaciones entre coordenadas geodésicas y ECEF, y
entre ECEF y ENU local.
I.2. Transformaciones de coordenadas
Este apartado ofrece las herramientas necesarias para poder representar la
posición de un punto en los tres sistemas de coordenadas definidos en la
sección I.1.2, es decir: geocéntricas (ECEF), geodésicas (latitud, longitud y
altura sobre el elipsoide WGS84), y local ENU.
I.2.1. Llh – ECEF
Esta sección incluye las transformaciones necesarias para obtener las
coordenadas cartesianas (sistema ECEF) de la posición de un punto a partir de
sus coordenadas geodésicas {φ0, λ0, h0}, y viceversa.
78 Earth Gravitational Model 2008, cuyos precursores son las versiones de 1996 y 1984.
ANEXO I
122
I.2.1.1 Llh a ECEF
Esta transformación permite obtener la posición de un punto P en el sistema
cartesiano ECEF {x, y, z}, a partir de su definición en las coordenadas
geodésicas {φ, λ, h} [89]:
( )
( )
( )[ ]
+−=
+=
+=
ϕ
λϕ
λϕ
sine1
sincos
coscos
2
hRz
hRy
hRx
N
N
N
( I.2 )
donde el radio de curvatura normal a la dirección del primer vertical de P (RN en
la Figura I.3) se obtiene mediante:
ϕ22 sin1 e
a
RN
−
= ( I.3 )
Por otra parte, el cuadrado de la excentricidad (e2) se calcula:
2
22
2e
a
ba −
= ( I.4 )
siendo a el semieje mayor (ya definido en I.1.1), y b el semieje menor. El
semieje menor (b) se obtiene a partir del semieje mayor (a) y del achatamiento
(f):
( )fab −= 1 ( I.5 )
Los valores de a y f asociados al elipsoide de referencia definido por WGS84
se encuentran en la sección I.1.1.
I.2.1.2 ECEF a Llh
Esta transformación es la operación inversa a la definida en I.2.2.1, permitiendo
obtener las coordenadas geodésicas {φ, λ, h} de un punto P, a partir de su
definición en el sistema cartesiano ECEF {x, y, z}.
En este apéndice presentamos un algoritmo basado en el método de Newton-
Raphson [90], que emplea el parámetro iterativo θ.
Definimos el valor p del siguiente modo:
22 yxp += ( I.6 )
Por otra parte, inicializamos θ del siguiente modo:
SISTEMAS DE COORDENADAS Y TRANSFORMACIONES
123
pb
za
arctan0 =θ ( I.7 )
En cada iteración, se calcula la latitud en función de dicho parámetro:
i
i
i
aep
bez
θ
θ
ϕ
32
32
cos
sin
arctan
−
′+
= ( I.8 )
El paso de la iteración consiste en actualizar el valor de θ:
=+ ii
a
b
θθ tanarctan1 ( I.9 )
Este método suele converger tras dos o tres iteraciones. Las coordenadas
geodésicas se calculan del siguiente modo:
−=
=
−
′+
=
NR
p
h
x
y
aep
bez
ϕ
λ
θ
θ
ϕ
cos
arctan
cos
sin
arctan
32
32
( I.10 )
donde RN viene definida en (I.3), y el cuadrado de la segunda excentricidad
(e’2) se calcula mediante:
2
22
2e
b
ba −
=′ ( I.11 )
I.2.2. ECEF – ENU
Esta sección incluye las transformaciones necesarias para obtener la posición
de un punto P con coordenadas cartesianas {x, y, z} (sistema ECEF) en el
sistema de coordenadas local ENU, cuyo origen se encuentra en el punto Q (tal
y como viene ilustrado en la Figura I.3) de coordenadas cartesianas {x0, y0, z0},
y viceversa[89].
I.2.2.1 ECEF a ENU
Las componentes del vector QP definen la posición de P en el sistema de
coordenadas local ENU centrado en Q. De este modo, la transformación
consiste en expresar el vector QP, del sistema ECEF al sistema ENU local.
ANEXO I
124
Esto se obtiene mediante la proyección del vector QP sobre la base de
vectores E, N y U definidos en el apartado I.1.2.3.
De este modo:
−
−
−
⋅=
= →
0
0
0
zz
yy
xx
R
u
n
e
P ENUECEFenu ( I.12 )donde RECEF→ENU es la matriz de rotación de ECEF a ENU a continuación:
−−−
−−
−
=→
00000
00000
00
sinsincoscoscos
cossinsincossin
0cossin
ϕλϕλϕ
ϕλϕλϕ
λλ
ENUECEFR
( I.13 )
I.2.2.2 ENU a ECEF
Esta transformación es la inversa a la descrita anteriormente. Así, para obtener
la posición del punto P – ubicado en {e, n, u} en el sistema de coordenadas
ENU centrado en Q - en coordenadas ECEF, el cálculo necesario es:
⋅+
=
−
→
u
n
e
R
z
y
x
z
y
x
ENUECEF
1
0
0
0
( I.14 )
donde RECEF→ENU
-1 es la matriz inversa de RECEF→ENU.
125
ANEXO II. MÉTODOS DE RESOLUCIÓN DEL
PROBLEMA DE LOCALIZACIÓN
HIPERBÓLICA
Este anexo abarca los métodos de resolución del problema de localización
hiperbólica definido por el siguiente sistema de ecuaciones:
( )
( )
( )
+−⋅=−
+−⋅=−
+−⋅=−
NjNjjNj
jjjj
jjjj
nRR
c
nRR
c
nRR
c
11
0
1
1313
0
13
1212
0
12
1
1
1
ττ
ττ
ττ
M
(II.1)
donde:
( ) ( ) ( )222
ijijijij zzyyxxR −+−+−= (II.2)
siendo (xi, yi, zi) la posición de cada estación, y (xj, yj, zj) la posición de la
aeronave – que es lo que debemos determinar. Dichas posiciones son relativas
a un sistema de coordenadas local (ENU), cuyo origen ha de ser especificado
por el usuario.
Por otra parte,
11 nnn ii −= (II.3)
y a su vez, n viene definida mediante una variable aleatoria gaussiana de
media nula y desviación típica σ.
Dichos métodos, enumerados en [24], están compuestos por un método
iterativo [23, 25] y dos métodos algebraicos, [21, 26].
ANEXO II
126
II.1. Método iterativo
II.1.1. Procedimiento
El método iterativo [23, 25] emplea el método del gradiente, basado en el
desarrollo en serie de Taylor de la expresión (II.1).
De este modo, la relación entre la posición estimada en una determinada
iteración y la estimación en la iteración anterior es:
[ ] [ ] [ ]kkk xδxx +−= 1 (II.4)
Donde x es el vector posición de la aeronave, y ix es el vector posición de la
estación i-ésima:
( )T
jjj zyx ,,=x (II.5)
En cada iteración, se estiman las TDOA a partir de la posición estimada
(vector [ ]1−kx ). Definimos el vector de error como la diferencia entre las
medidas de TDOA y sus estimaciones:
[ ] [ ]( )xxf −−−= 1kk τε (II.6)
donde f viene definida del siguiente modo:
( ) ( )TNjjj fff ,...,, 32=− xxf
( )jijij RR
c
f 1
0
1
−=
(II.7)
El incremento diferencial [ ]kxδ se obtiene resolviendo los términos de primer
orden del desarrollo en serie de Taylor de la expresión (II.1).
Así, (II.6) puede reexpresarse del siguiente modo:
[ ] [ ] [ ]kkk δxAε ⋅= (II.8)
donde A[k] es el jacobiano del sistema de ecuaciones (II.1) en la k-ésima
iteración respecto a las componentes del vector x. Por tanto, la matriz A es:
−
−
−−
−
−−
−
−
−
−
−−
−
−−
−
−
−
−
−−
−
−−
−
−
=
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
=
j
j
Nj
Nj
j
j
Nj
Nj
j
j
Nj
Nj
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
j
Nj
j
Nj
j
Nj
j
j
j
j
j
j
j
j
j
j
j
j
R
zz
R
zz
R
yy
R
yy
R
xx
R
xx
R
zz
R
zz
R
yy
R
yy
R
xx
R
xx
R
zz
R
zz
R
yy
R
yy
R
xx
R
xx
z
f
y
f
x
f
z
f
y
f
x
f
z
f
y
f
x
f
1
1
1
1
1
1
1
1
3
3
1
1
3
3
1
1
3
3
1
1
2
2
1
1
2
2
1
1
2
2
333
222
MMMMMM
A
(II.9)
MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN HIPERBÓLICA
127
Si el número de estaciones es superior a 4, el sistema de ecuaciones se vuelve
sobredeterminado, con lo que éste ha de ser resuelto mediante la expresión del
mínimo del error cuadrático medio, utilizando la matriz pseudoinversa a
continuación[25]:
( ) εSAASAδx 111 −−−= TT (II.10)
donde la matriz S es la matriz de covarianza de los errores de medida de
TDOA. Considerando que dichos errores de medida tienen misma distribución
estadística, a partir de la expresión (II.3), podemos deducir que dicha matriz es:
=
2
2
2
200
020
002
σ
σ
σ
L
MOMM
L
L
S
(II.11)
Al tratarse de un algoritmo iterativo, el autor sugiere el fin del bucle:
• cuando la norma de las componentes posicionales de δx se reduce por
debajo de un determinado umbral (por ejemplo, 50 metros multiplicado
por el número de aeronaves).
• al rebasarse un número de iteraciones (por ejemplo, 10)
• cuando la norma del vector δx en una iteración es superior a la anterior.
II.1.2. Ejemplos
El siguiente ejemplo está destinado a ilustrar la capacidad y las limitaciones de
este método. Para ello, consideramos el siguiente escenario bidimensional con
una aeronave (emisor) y tres estaciones:
• las estaciones están ubicadas en (-10000; 0), (0; 0) y (10000; 0).
• la aeronave está ubicada en (0; 0), con lo que la posición estimada es
exactamente el error en la estimación.
• el error de la medida del tiempo de llegada se modela mediante una
variable aleatoria gaussiana, de media nula y desviación típica 10ns.
II.1.2.1 Caso de convergencia
Dado que el algoritmo iterativo requiere una estimación inicial, se parte de la
hipótesis inicial de que la aeronave está ubicada en (5000; 5000).
La Tabla II.1 a continuación muestra la posición estimada en cada iteración.
ANEXO II
128
k x[k] y[k]
1 5000.000 5000.000
2 -1908.239 -1908.200
3 75.688 75.748
4 0.994 1.050
5 0.998 1.054
Tabla II.1. Estimación de la posición en cada iteración – caso convergente.
II.1.2.2 Caso de no convergencia
En este caso empleamos otra hipótesis inicial. Asumimos que la aeronave está
ubicada en (10000; 10000), y ejecutamos el método iterativo. Las estimaciones
en cada paso se reflejan en la siguiente Tabla II.2.
k x[k] y[k]
1 10000.000 10000.000
2 -26180.914 -26171.157
3 669003.193 668834.673
4 -11.772·109 -11.769·109
5 -4.726·1015 -4.724·1015
Tabla II.2. Estimación de la posición en cada iteración – caso divergente
Se observa que la posición estimada diverge desde la primera iteración (k=2).
A partir de ambos ejemplos, se constata que este método de localización
hiperbólica no es infalible. En caso de éxito, la convergencia se alcanza
rápidamente. Por otra parte, si la hipótesis inicial es desfavorable, la
divergencia del algoritmo es fácil de detectar.
II.2. Métodos algebraicos
II.2.1. Algoritmo de Bancroft
El método propuesto por Bancroft [26] asume que las medidas obtenidas son
de tiempo de llegada de la señal (TOA).
MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN HIPERBÓLICA
129
Definiendo la pseudodistancia como el producto del TOA por la velocidad de la
luz (ρij = c0·τij) - siendo j el índice de la aeronave que emite la señal, e i es el
índice de la estación base - el sistema de ecuaciones es el siguiente:
iijij nR +=ρ (II.12)
Se define el vector a, compuesto por la posición de cada estación base y el
error en la medida del TOA:
[ ]
iiiii nzyx=a (II.13)
La matriz A, a su vez, está compuesta por los vectores ai mencionados
anteriormente:
=
NNNN zyx
zyx
zyx
zyx
ρ
ρ
ρ
ρ
MMMM
3333
2222
1111
A (II.14)
donde N es el número de estaciones base que ha recibido la misma señal.
Por otra parte, se define i0 como el vector columna unidad, cuyo número de
elementos es igual al número de estaciones:
=
1
1
1
1
0
M
i (II.15)
y el vector columna r, que se obtiene del siguiente modo:
−++
−++
−++
−++
=
2222
2
3
2
3
2
3
2
3
2
2
2
2
2
2
2
2
2
1
2
1
2
1
2
1
2
1
NjNNN
j
j
j
zyx
zyx
zyx
zyx
ρ
ρ
ρ
ρ
M
r (II.16)
La matriz B, es la inversa generalizada:
( ) WAWAAB T-T 1= (II.17)
ANEXO II130
donde W es la inversa de la matriz de covarianza del error de medida de TOA:
1
2
2
3
2
2
2
1
000
000
000
000
−
=
Nσ
σ
σ
σ
L
MOMMM
L
L
L
W (II.18)
Después, se obtienen los vectores u y v mediante:
0iBu ⋅= (II.19)
rBv ⋅= (II.20)
Estos vectores permiten obtener los coeficientes asociados a la siguiente
ecuación de segundo orden:
02 =++ GFE λλ (II.21)
donde:
2
4
2
3
22
21
uuuuE −++= (II.22)
144332211 −−++= vuvuvuvuF (II.23)
2
4
2
3
22
21
vvvvG −++= (II.24)
Definimos λ1 y λ2 a las raíces de la ecuación (II.21). Estas raíces nos permiten
obtener el vector aumentado – que contiene la posición y el error en la medida
de tiempo - asociado a la aeronave del siguiente modo:
vuy += 2,12,1 λ (II.25)
donde y1 e y2 son los dos vectores candidatos (nótese la ambigüedad), que son
vectores columna de 4 elementos. La posición estimada de la aeronave se
obtiene mediante identificación de las tres primeras componentes:
=
j
j
j
k
k
k
z
y
x
y
y
y
3,
2,
1,
(II.26)
La ambigüedad se elimina calculando las pseudodistancias reutilizando cada
posición estimada.
MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN HIPERBÓLICA
131
Este método permite obtener la solución de manera cerrada, con lo que el
riesgo de la no convergencia está mitigado. Sin embargo, la limitación asociada
a este método es que las raíces de (II.21) pueden ser complejas. En este caso,
se aconseja emplear la parte real de la solución, y comprobar la validez de
dicho resultado, bien mediante un filtro de seguimiento (si el resultado obtenido
es muy diferente de la estimación anterior, se rechaza), o bien comprobando el
resultado con el obtenido mediante otro método.
II.2.2. Algoritmo de Chan, Ho
Este método [21] permite encontrar la solución de forma cerrada, y es válida
tanto para los casos en los que la aeronave está cerca, como para los casos en
que se encuentra alejada de la constelación. Cuando los errores asociados a la
medida del TDOA son pequeños, este método es una aproximación al
estimador de máxima verosimilitud.
Introduciendo:
( )[ ]2112 jjijij RRRR +−= (II.27)
podemos reescribir (II.2) del siguiente modo:
( ) ( ) 2222222
111
2
1 2222 jjjjijijiiiijjjijjij zyxzzyyxxzyxRRRRRR +++−−−++=+−+− (II.28)
Restando la distancia de la aeronave a la estación de referencia a cada lado:
( ) ( ) ( ) ( )2121212221121 2222 zyxzyxzzyyxxRRRRR jjjjijijijjijjij ++−+++−−−=−+− (II.29)
En el caso de un sistema determinado ( { }4,3,2,1∈i ), la solución de (II.29) viene
dada por:
( )
( )
( )
+−−
+−−
+−−
+
−
−
−
−−−
−−−
−−−
−=
−
2
1
2
4
2
14
2
1
2
3
2
13
2
1
2
2
2
12
1
14
13
12
1
141414
131313
121212
2
1
jjjj
jjjj
jjjj
j
jj
jj
jj
j
j
j
RRRR
RRRR
RRRR
R
RR
RR
RR
zzyyxx
zzyyxx
zzyyxx
z
y
x
(II.30)
Sin embargo, en presencia de ruido en la medida de TDOA, la solución al
problema es aquel punto que mejor se ajuste a las medidas de TDOA. Para ello,
definimos la variable za del siguiente modo:
=
j
j
j
j
R
z
y
x
1
az
(II.31)
ANEXO II
132
El vector de errores ψ que se obtiene a partir de (II.29) es:
0
aa zGh −=ψ (II.32)
donde:
( )
( )
( )
+−−
+−−
+−−
=
2
1
22
1
2
1
2
3
2
13
2
1
2
2
2
12
2
1
jNjjNj
jjjj
jjjj
RRRR
RRRR
RRRR
M
h
−−−−
−−−−
−−−−
−−−−
−=
jNjNNN
jj
jj
jj
RRzzyyxx
RRzzyyxx
RRzzyyxx
RRzzyyxx
1111
14141414
13131313
12121212
MMMM
aG
(II.33)
siendo el superíndice 0 el valor sin ruido (caso ideal).
Considerando que:
( ) ( ) ( )10110 nncRRc ijijjij −⋅+−=−⋅ ττ
( ) jjijij RRRR 11 +−=
(II.34)
el vector de error ψ es:
nnBn ⊗+= 200
2
1
ccψ
{ }
Njjj RRRdiag ,...,, 32=B
(II.35)
donde ⊗ es el producto elemento a elemento.
Para valores elevados de relación señal a ruido, el error asociado a las
medidas de TDOA tiende a ser gaussiano. Además, en la práctica se cumple
que Rij>>c0·ni, siendo así ψ un vector de error gaussiano, cuya covarianza es:
[ ] TT cE BQB20== ψψΨ (II.36)
donde Q es la matriz de covarianzas asociada a las medidas de TDOA.
Los elementos de za están relacionados entre sí mediante (II.2), con lo que
(II.34) es un sistema de ecuaciones no lineales, cuyas variables son xj, yj y zj.
MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN HIPERBÓLICA
133
II.2.2.1 Primera etapa: los elementos de za son independientes.
Este problema se aborda asumiendo primero que las variables xj, yj, zj y R1j son
independientes. Este primer paso se resuelve mediante mínimos cuadrados. La
solución final al problema se obtiene imponiendo la relación (II.2) a la primera
aproximación obtenida, y resolviendo mediante mínimos cuadrados. Este
procedimiento en dos pasos es una aproximación al estimador de máxima
verosimilitud para la localización de la aeronave transmisora.
Considerando que los elementos de za son independientes, el estimador ML
79
del vector za (suponiendo ruido gaussiano) es:
( ) ( ){ } ( ) hGGGzGhzGhz aaaaaaaa 1
111minarg −
−−− =−−= ΨΨΨ TT
T
(II.37)
que es la solución generalizada de (II.31) mediante mínimos cuadrados.
Sin embargo, ψ no es conocida, porque B contiene la información de las
distancias entre la aeronave y las estaciones base, con lo que es necesario
efectuar otra aproximación para resolver el problema.
Esta segunda aproximación consiste en asumir que la aeronave está lejos de la
constelación. Esto significa que los valores de Rij son similares para cualquier
valor de i, pudiendo efectuarse la aproximación B≈r0I , donde r0 es la distancia
de la aeronave a la estación de referencia (estación 1), y la matriz I es la matriz
identidad que tiene N-1 elementos, donde N es el número de estaciones.
Puesto que los escalados de la matriz ψ no afectan a la solución, (II.37) se
puede aproximar mediante:
( ) hGGGz aaaa 1
11 −−−≈ QQ TT (II.38)
Si la aeronave se encuentra cerca de la solución, el procedimiento consiste en
emplear (II.38) para tener una estimación inicial de la matriz B, y
posteriormente emplear (II.37).
II.2.2.2 Segunda etapa: relación entre los elementos de za.
La primera etapa nos permite obtener el vector za, que es la solución al
problema asumiendo que sus elementos son independientes. Sin embargo, la
estimación debe mejorar al introducir la relación entre la posición de la
aeronave, y la distancia entre la aeronave y la estación base de referencia. Por
ello, el algoritmo propone la siguiente iteración.
79 ML son las siglas del término inglés “Maximum Likelihood”, que significa “máxima
verosimilitud”.
ANEXO II
134
Asumiendo que el ruido en las medidas de TDOA es despreciable, el vector za
es un vector aleatorio centrado en la solución, cuya matriz de covarianza
asociada se obtiene de la siguiente forma:
( ) ( ) 110cov −−≈ aaa GGz ΨT (II.39)
Así, za puede expresarse del siguiente modo:
+=
+=
+=
+=
4
0
14,
3
0
3,
2
0
2,
1
0
1,
eRz
ezz
eyz
exz
ja
ja
ja
ja
(II.40)
Restando la posición de la estación a las tres primeras componentes de za, y
elevando al cuadrado, obtenemos el siguiente sistema de ecuaciones:
0
aa zGh ′′−′=′ψ (II.41)
donde:
( )
( )
( )
−
−
−
=′
2
4,
2
13,
2
12,
2
11,
a
a
a
a
z
zz
yz
xz
h
=′
111
100
010
001
aG
( )
( )
( )
−
−
−
=′
2
1
2
1
2
1
zz
yy
xx
j
j
j
az
(II.42)
Sustituyendo (II.40) en (II.41), se obtiene80:
80Las aproximaciones son válidas puesto que se asume que la magnitud del error es
despreciable, en comparación a la distancia entre la aeronave y la estación de referencia.
Esta es otra aproximación al procedimiento de máxima verosimilitud.
MÉTODOS DE RESOLUCIÓN DEL PROBLEMA DE LOCALIZACIÓN HIPERBÓLICA
135
( ) ( )
( ) ( )
( ) ( )
≈+=′
−≈+−=′
−≈+−=′
−≈+−=′
4
0
1
2
44
0
14
31
02
331
0
3
21
02
221
0
2
11
02
111
0
1
22
22
22
22
eReeR
ezzeezz
eyyeeyy
exxeexx
jj
jj
jj
jj
ψ
ψ
ψ
ψ
(II.43)
Así, la matriz de covarianzas de ψ’ es:
[ ] ( ) TaTE BzB ′′=′′=′ cov4ψψΨ
{ }01101010 ,,,diag jjjj Rzzyyxx −−−=′B
(II.44)
De este modo, la estimación ML de za’ es:
( ) hGGGz aaaa ′′′′′′≈′ −
−− 111
ΨΨ
TT
(II.45)
Aunque ψ’ no sea conocido, B’ puede aproximarse a partir de (II.36) y (II.39).
De este modo, la estimación final de la posición ofrece dos candidatos:
( ) ( )TT
jjj zyxzyx 111+′= az
y
( ) ( )TTjjj zyxzyx 111+′−= az
(II.46)
Al igual que en el método presentado por Bancroft, es necesario comprobar
cuál de las dos posibles soluciones es válida. La solución final es aquella que
define la posición de la aeronave de tal modo que las TDOA estimadas se
acerquen más a las medidas.
Este algoritmo es especialmente útil cuando la aeronave se encuentra lejos de
la constelación. Sin embargo, es más costoso que el método presentado por
Bancroft. Por otra parte, al igual que en el método de Bancroft es necesario
comprobar cuál de las dos soluciones obtenidas es la correcta.
136
ANEXO III. UNSCENTED KALMAN FILTER
III.1. Introducción
El filtro de Kalman es la estimación óptima para modelos de sistemas lineales
con ruido blanco tanto la fase de transición (vector de estado del blanco de un
instante, a partir del anterior), como en las medidas. No en vano, se trata de
uno de los métodos más empleados en los problemas de seguimiento y
estimación. Sin embargo, su aplicación a sistemas no lineales puede resultar
difícil. En estos casos, se suele emplear el Extended Kalman Filter (EKF), que
consiste en aproximar cada modelo no lineal empleado en el filtro mediante una
aproximación de primer orden, con el fin de poder aplicar el filtro de Kalman.
III.2. Planteamiento del problema
El objetivo es aplicar un filtro de Kalman a un sistema discreto no-lineal del tipo:
( ) ( ) ( ) ( )[ ]
( ) ( ) ( )[ ] ( )kkkkk
kkkkk
wuxhz
vuxfx
+=
=+
,,
,,,1
(III.1)
donde:
• X(k) es el vector de estado (de dimensión n) del sistema en el instante k
• u(k) es el vector de control
• v(k) es el vector que caracteriza el proceso de ruido (dimensión q)
• z(k) es el vector de observaciones
• w(k) es el vector de ruido asociado a la medida
UNSCENTED KALMAN FILTER
137
Se asume que los vectores de ruido v y w están compuestos por variables
aleatorias gaussianas de media nula. Además, estos vectores son
independientes entre sí, y sus matrices de covarianzas son las matrices
diagonales Q y R respectivamente.
El filtro de Kalman propaga los dos primeros momentos (media y varianza) de
la distribución de X(k) de manera recursiva, en dos etapas: estimación y
actualización.
En la práctica, el uso del EKF en sistemas no lineales tiene dos inconvenientes
[69]:
• la linearización puede producir filtros altamente inestables si el sistema
no tiene un comportamiento cuasi-lineal cerca de la solución.
• la obtención de las matrices jacobianas no es trivial en la mayoría de las
aplicaciones, siendo habitualmente una fuente de problemas a la hora
de la implementación.
Una de las alternativas al EKF es el empleo del Unscented Kalman Filter [69],
que aplica la filosofía de la Unscented Transform (UT) al filtro de Kalman. Este
método se presenta en el siguiente apartado.
III.3. Unscented Transform
La Unscented Transform (UT) presenta un método para calcular las
estadísticas de una variable aleatoria procesada mediante una transformación
no lineal. Se basa en la idea de que es más fácil aproximar una distribución
gaussiana que una función no lineal.
Figura III.1. Principio de la Unscented Transform [69].
Se escoge una serie de puntos (denominados puntos sigma), de modo que x y
Pxx son su media y covarianza respectivamente. Se aplica la función no lineal a
cada punto, con el fin de obtener una nube de puntos a la salida. Definimos la
media y la varianza de dichos puntos mediante y y Pyy respectivamente.
ANEXO III
138
Si bien a simple vista esta técnica puede parecerse a un método de Monte
Carlo, existe una diferencia fundamental: las muestras no se han escogido de
manera aleatoria, sino mediante un algoritmo determinista.
La variable aleatoria X de dimensión n, con media x y covarianza Pxx se
aproxima mediante los 2n+1 puntos ponderados definidos a continuación:
( )( )
( )( )
=+−=
=++=
=
+ nqPn
nqPn
qXXqn
qXXq
,...,1,
,...,1,
0
κ
κ
xX
xX
xX
(III.2)
donde κ es un número real, ( )( )
qXX
Pn κ+ es la fila o columna q-ésima de la
matriz raíz cuadrada de ( ) XXPn κ+ . Por otra parte, las constantes w
corresponden a los pesos asociados a cada punto sigma:
( )
( )
+
=
+
=
+
=
+
κ
κ
κ
κ
n
w
n
w
n
w
qn
q
2
1
2
1
0
(III.3)
Así, el proceso de transformación se compone de tres fases:
• la primera consiste en aplicar la función de transformación a cada punto
sigma, con el fin de obtener los puntos sigma transformados.
( )ii Xfy = (III.4)
• la segunda fase consiste en obtener la media ponderada de los puntos a
la salida de la función de transformación:
∑
=
=
n
i
iiw
2
0
yy (III.5)
• finalmente, la covarianza se obtiene mediante el producto exterior:
{ }{ }Ti
n
i
iiYY w yyyyP −−=∑
=
2
0
(III.6)
Tal y como se adelantó al principio de este apartado, la clave de este algoritmo
está en aproximar la distribución de la variable compuesta X, y no la función f.
Una vez explicado este principio, abordaremos su aplicación al filtro de Kalman
en el siguiente apartado.
UNSCENTED KALMAN FILTER
139
III.4. Unscented Kalman Filter
El Unscented Kalman Filter (UKF) es una extensión de la UT a la estimación
recursiva. El vector de estado se define como la concatenación del vector de
estado y los vectores de ruido:
( ) ( ) ( ) ( )[ ]TTTTa nvxx kkkk = (III.7)
En el caso particular (pero habitual) en el que tanto el ruido de proceso como el
ruido de medida son puramente aditivos, la complejidad computacional del UKF
puede reducirse. En este caso, el vector de estado no necesita ser aumentado
con los vectores de ruido, reduciendo así el número de dimensiones, y por
ende el número de puntos sigma81.
A continuación se presentan los pasos a seguir para una iteración del UKF: los
pasos 1 a 4 son necesarios para la etapa de estimación, donde se obtienen
( )1ˆ −kkx y ( )1−kkP . A dicha etapa le sigue la etapa de filtrado, cuyos
resultados son el vector de estado y su matriz de covarianza asociada, una vez
incorporada la información proporcionada por las observaciones (pasos 5 a 8).
1. Inicialización (k=0):
( ) ( )[ ]
( ) ( ) ( )( ) ( ) ( )( )[ ]
−−=
=
T
E
E
0ˆ00ˆ00
00ˆ
xxxxP
xx
(III.8)
2. Para una muestra dada (k≥1), calcular los puntos sigma:
( ) ( ) ( ) ( ) ( ) ( )[ ]11ˆ11ˆ1ˆ1 −−−−+−−=− kkkkkk PxPxx γγχ
(III.9)
donde:
( )καγ += n
(III.10)
El parámetro de ajuste α determina la dispersión de los puntos sigma alrededor
del valor medio, y suele ser un valor ligeramente positivo (10-3). Por otra parte,
κ es un parámetro de escalado secundario que suele ser 0.
81 En el caso en el que el ruido no pueda considerarse puramente aditivo, se invita al lector a
referirse a [78]
ANEXO III
140
3. Se transforma cada punto sigma obtenido mediante la función de
actualización f:
( ) ( ) ( )[ ]kkkk uf ,1 χχ =− (III.11)4. Se calcula el vector de estado y la matriz de covarianza empleando las
medias ponderadas:
( ) ( )
( ) ( ) ( )[ ] ( ) ( )[ ]
+−−−−−−=−
−=−
∑
∑
=
=
QxxP
x
x
Tq
n
q
qq
c
n
q
qq
m
kkkkkkkkWkk
kkWkk
1ˆ11ˆ11
11ˆ
2
0
2
0
χχ
χ
(III.12)
donde:
( )
≠
+
=
+=
0 si,
n2
1
0 si,
n
q
q
W
q
m
λ
λ
λ
(III.13)
( )
( )
≠
+
=+−+
+=
0 si,
n2
1
0 si,1
n
2
q
q
W
q
c
λ
βα
λ
λ
(III.14)
siendo:
( ) nn −+= καλ 2 (III.15)
El parámetro β se emplea para incorporar la información previa acerca de la
distribución de x. En el caso de distribuciones gaussianas, el ajuste óptimo es
β=2.
5. Se estiman las observaciones correspondientes a los puntos sigma
generados a partir de ( )1ˆ −kkx y ( )1−kkxP . Para ello, se generan los puntos
sigma:
( ) ( ) ( ) ( ) ( ) ( )[ ]11ˆ11ˆ1ˆ −−−−+−−= kkkkkkkkkkk xx PxPxx γγχ
(III.16)
6. y se aplica el modelo de observación h a cada punto sigma:
( ) ( )[ ]11 −=− kkkky χh (III.17)
UNSCENTED KALMAN FILTER
141
7. En este paso, se calcula la media ponderada de la observación asociada a
cada punto sigma:
( ) ( )∑
=
−=−
n
q
qq
m kkyWkk
2
0
11ŷ (III.18)
8. Finalmente, se incorpora el vector de observaciones, y(k):
( ) ( )[ ] ( ) ( )[ ]
( ) ( )[ ] ( ) ( )[ ]
( ) ( ) ( ) ( )[ ]
( ) ( )
−−=
−−+−=
=
−−−−−−=
+−−−−−−=
−
=
=
∑
∑
T
yyxx
yyxy
xy
yy
KKPPP
yKxx
PPK
yxP
RyyP
1
1ˆ1ˆˆ
1ˆ11ˆ1
1ˆ11ˆ1
1
2
0
2
0
kkkk
kkkykkkk
kkkkykkkkW
kkkkykkkkyW
Tq
n
q
qq
c
Tq
n
q
qq
c
χ
(III.19)
Este apartado ha presentado los pasos del Unscented Kalman Filter. En el
apartado siguiente se comparan los resultados a la salida de un Extended
Kalman Filter (EKF) y un UKF, dada una función de observación no lineal.
III.5. Ejemplo
III.5.1. Problema
Este problema consiste en estimar la posición de un blanco de movimiento
rectilíneo uniforme, a partir de una función de observación que corresponde a
la función seno de la posición (unidimensional). Sin embargo, existe un ruido
asociado a la medida, que puede asumirse gaussiano y de media nula.
Plasmando el problema matemáticamente, se trata del siguiente par de
funciones:
( ) ( )[ ]
( ) ( )[ ] ( )kkk
kk
wxhz
xfx
+=
=+ 1
(III.20)
donde:
( ) ( ) ( )[ ]
( )
( ) ( )[ ] RkwkwE
vkv
kvkxk
T
xx
T
x
=
≡
=x
(III.21)
En el siguiente apartado se presenta la resolución al problema mediante EKF.
ANEXO III
142
III.5.2. Resolución mediante EKF: desarrollo
El primer paso consiste en la estimación del vector de estado y su matriz de
covarianzas asociada, en base a la función de transición y a la solución de la
muestra anterior.
Empleando un modelo de maniobra que caracteriza un movimiento rectilíneo
uniforme, la matriz de transición F equivalente a la función f es:
∆
=
10
1 T
F (III.22)
donde ∆T es el tiempo entre dos muestras consecutivas.
Así, el vector de estado se obtiene mediante:
( )
( )
( )
( )
( )
( )
−
∆
=
−
−
∆
=
=−
xxx v
kxT
kv
kxT
kv
kx
kk
1
10
1
1
1
10
1
1x (III.23)
Y la matriz de covarianzas asociadas:
( ) ( ) ( )
∆
−−
∆
=−−=−=
10
1
11
10
1
111
T
kk
T
kkkk xxxxxxxx PFFPPP (III.24)
El segundo paso consiste en incorporar la información de la medida:
( ) ( ) ( ) ( )[ ]
( ) [ ] ( )
−=
−−+−=
=
+=
−
1
11ˆˆ
1
kkJkk
kkJkykkkk
J
RJJ
T
T
xx
xx
xx
PK-IP
xKxx
SPK
PS
(III.25)
donde J es la matriz jacobiana de la función h particularizada en x(k):
( )
( )
∂
∂
∂
∂
=
−=
−=
1
1
kkxxx
kkxx
v
h
x
h
J (III.26)
III.5.3. Resolución mediante UKF
Por otra parte, la resolución del problema mediante UKF se obtiene
implementando los pasos descritos en la sección III.4.
UNSCENTED KALMAN FILTER
143
III.5.4. Experimento
Con el fin de comparar el rendimiento de ambos métodos (EKF y UKF), se ha
realizado un experimento en base a los datos del problema expuesto en III.5.1.
Las condiciones están detalladas en la siguiente tabla:
Condición Valor unidad
Duración 2 s
∆T 10-2 s
x inicial 0 m
vx π/2 m·s
-1
h(x) 0,5 sin(x) no aplica
R 10-4 no aplica 82
Arranque en frío 3 muestras
Tabla III.1. Condiciones experimentales.
Cabe destacar que ambos filtros necesitan al menos 2 observaciones con el fin
de estimar la velocidad inicial del blanco. Por ello, ambos filtros únicamente
pueden entrar en juego a partir de la tercera muestra de tiempo.
En base a estos datos, la Figura III.2 muestra las observaciones durante una
ejecución del experimento.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
t (s)
h(
x)
-
S
in
u
ni
da
d
Figura III.2. Valor de las medidas en función del tiempo.
Como puede observarse, la función de observación no solamente no es lineal,
sino que también consta de un máximo local.
82 La unidad asociada a R es el cuadrado de la unidad asociada a la medida
ANEXO III
144
La Figura III.3 a continuación compara la posición del blanco (x) estimada por
cada filtro, junto con el valor real.
Se observa que la solución EKF es ruidosa e inestable para valores de t entre
0,8s y 1,2s. Precisamente, se trata del intervalo de tiempo centrado en t=1s,
momento en el que la no linealidad es más fuerte (el valor absoluto de h’’(x) es
máximo).
Por otra parte, se puede apreciar la robustez de UKF durante el tiempo del
experimento.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.5
1
1.5
2
2.5
3
3.5
t (s)
V
al
or
e
st
im
ad
o
de
x
(
m
)
EKF
UKF
Real
Figura III.3. Comparación de las soluciones EKF y UKF con el valor real.
La Figura III.4 a continuación permite comparar directamente el error cometido
en la estimación de la posición.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
t (s)
E
rr
or
e
n
la
e
st
im
ac
io
n
de
x
(
m
)
EKF
UKF
Figura III.4. Error en la estimación de posición con EKF y UKF.
Este anexo presenta la transformada Unscented (base del UKF), el desarrollo
de UKF y un experimento ilustrativo que compara el rendimiento de UKF y EKF
con un modelo de observación no lineal.
145
ANEXO IV. COTA DE CRAMÉR-RAO PARA
LOCALIZACIÓN HIPERBÓLICA
Este anexo muestra el desarrollo empleado en [27] para obtener la expresión
de la cota de Cramér-Rao para localización hiperbólica. La cota (inferior) de
Cramér-Rao (abreviada CRB o CRLB83), expresa una cota inferior para la
varianza de un estimador de un parámetro determinista.
IV.1. Localización esférica
La localización esférica es un método de localización basado en el tiempo de
llegada de una misma señal (TOA) transmitida por un emisor, y recibida por
varias estaciones base. Matemáticamente, se trata de un problema más
sencillo que la localización hiperbólica (basada en TDOA).
En el caso de la localización esférica, la posición del emisor se encuentra (en
un caso ideal) en la intersección de las esferas centradas en cada estación
base (receptora), y cuyo radio es el producto de la velocidad de la luz en el
vacío (caso ideal) por el TOA de la señal a cada estación base.
De este modo, el tiempo de llegada de la señal (TOA) emitida por una
aeronave j, y recibida por la estación base i, viene definido mediante:
( ) ( ) ( )
iijijijiijij ezzyyxx
c
eR
c
+−+−+−⋅=+⋅=
222
00
11
τ (IV.1)
83 siglas de los términos ingleses: Cramér-Rao Bound y Cramér-Rao Lower Bound, dicha cota
es llamada así en honor a Harald Cramér y Calyampudi Radhakrishna Rao.
ANEXO IV
146
donde Rij es la distancia euclídea entre la estación base i y la aeronave j, y ei
es el error asociado a la medidadel TOA en la estación base i.
El receptor estima la posición de la aeronave (x) en base a las medidas τ,
pudiendo establecerse la siguiente relación:
( ) exhτ += (IV.2)
donde h es la función de observación, que permite estimar el TOA en función
de la posición del blanco y de las estaciones base, mientras que e es el vector
formado por los errores de medida de TOA.
Sea R la matriz de covarianzas asociada al vector e; se presupone que los
errores de medida son independientes entre sí y de media nula (dando lugar a
que R sea una matriz diagonal).
La cota de Cramér-Rao proporciona el límite inferior de la covarianza asociada
a la estimación de la posición. Para calcular este límite, es necesario linearizar
la expresión que relaciona el TOA con la posición del blanco, particularizada en
la posición estimada del blanco (x0, y0, z0):
eCxτ += (IV.3)
donde C es la matriz jacobiana de τ respecto de x:
o
o
o
o
o
o
zz
yy
xx
ij
ij
ij
ij
ij
ij
zz
yy
xxi R
zz
R
yy
R
xx
czyx
=
=
=
=
=
=
−−−
=
∂
∂
∂
∂
∂
∂
=
0
1τττ
C (IV.4)
En este caso, la CRLB asociada a la estimación de posición es:
[ ] 11CRLB −−= CRCT (IV.5)
IV.2. Localización hiperbólica
Por su parte, el sistema de localización hiperbólica es un método de
localización basado en la diferencia entre los tiempos de llegada de una misma
señal transmitida por un emisor, y recibida por varias estaciones base.
Definiendo ∆τ como el TDOA de una señal emitida por una aeronave j y
recibida en dos estaciones base i y 1.
jij
j
i 11 τττ −=∆ (IV.6)
COTA DE CRAMÉR-RAO
147
La ecuación (IV.2) puede modificarse con el fin de establecer una relación entre
la posición del blanco y el TDOA, multiplicando por la matriz M, definida a
continuación:
−
−
−
=
10001
0
0101
0011
OM
L
M (IV.7)
De este modo, la ecuación de posicionamiento hiperbólico es:
eMCxMτ += (IV.8)
Finalmente, estableciendo un paralelismo con la relación entre las expresiones
(IV.3) y (IV.5), la CRLB asociada a la estimación de la posición mediante
localización hiperbólica es:
( )[ ] 11CRLB −−= MCMRMMC TTT (IV.9)