Desarrollo de algoritmos usando datos GNSS

¡Bienvenidos una vez más a este rinconcito dedicado a la geomática!

En primer lugar deciros que he estado un poco ausente las últimas semanas (exámenes, trabajos y demás cosas con las que no os quiero aburrir) pero ya estoy aquí de nuevo con las pilas cargadas.

Este pequeño parón me ha permitido organizar algunos aspectos del blog, la semana que viene os contaré todas las novedades (no os lo perdáis 😉) y como la semana que viene (último domingo del mes) va a estar ocupada vamos a adelantar a hoy nuestra sección especial de #nuestros trabajos, así que, ahora si, ¡empezamos!

Esta vez vamos a presentar un trabajo realizado por Jorge Hernández Olcina, ingeniero en Geomática y Topografía y estudiante del máster en ingeniería Geomática y Geoinformación en la UPV. Os dejo con el interesante trabajo.

Título: Desarrollo de algoritmos usando datos GNSS

Autor: Jorge Hernández Olcina

Tutor: Ana Belén Anquela Julián

Tipo de trabajo: Trabajo fin de grado

Una de las tecnologías más utilizadas en la actualidad son los sistemas de navegación global por satélite (GNSS), tanto a nivel de usuario como a nivel profesional, ya que permiten desarrollar una gran cantidad de aplicaciones, como pueden ser sistemas de navegación y guiado, localización, cartografiado, gestión del territorio etc. Éstos son ejemplos que cada vez están más presentes en las vidas cotidianas de las personas, pues las infraestructuras globales creadas para el uso del GNSS permiten la aparición de multitud de explotaciones profesionales.

La motivación de este trabajo es profundizar en el funcionamiento del GNSS (dentro de muy poquito este contenido también aparecerá en el blog), además de comprender y adquirir los conocimientos básicos de cálculo que llevan integrados algunos dispositivos que usan esta tecnología.

En el presente trabajo se profundiza en la metodología de cálculo a partir de algoritmos matemáticos que permiten obtener las coordenadas de un dispositivo que incorpore la tecnología GNSS.

Datos de partida

Para ello se va a partir de una rutina de algoritmos existente, a partir de la cual se van a modificar dichos programas además de desarrollar nuevas rutinas de cálculo que se adapten a las necesidades deseadas.

Resumiendo, el trabajo consta de las siguientes partes:

- Modificar la estructura de algoritmos y scripts existentes con el fin de adaptarlos a los datos actuales. La rutina dada fue programada con objeto de cálculo de la posición de un receptor con efemérides transmitidas, en el año 2001. Actualmente los ficheros de observables y navegación han cambiado, por consiguiente, será necesario adaptar la rutina a los nuevos ficheros.

The GPS Easy Suite–Matlab code for the GPS newcomer, Kai Borre: [fuente: Artano]
The Easy Suite es una colección de scripts de Matlab que tienen como finalidad comprender  los aspectos básicos del sistema GPS. Vienen dados en dos colecciones:

• GPS Easy Suite I (2003) (10 scripts): Permiten realizar operaciones básicas como calcula posición de los satélites con efemérides transmitidas, cambios de formato de fechas, cálculo de coordenadas a partir de pseudodistancias etc. En ellos se encuentra la rutina que se va a emplear en el presente trabajo:

• “easy3.m”. Cálculo de la posición del receptor mediante posicionamiento absoluto con código (SPP - standard point positioning).

- Generar nuevos scripts y algoritmos que realicen el mismo procedimiento de cálculo que la rutina dada solo que en este caso, los datos de partida son efemérides precisas.

• En este trabajo uno de los objetivos es obtener las coordenadas de los satélites cada 30 segundos. Debido a que los ficheros de efemérides precisas solo facilitan las posiciones de los satélites cada 15 minutos, es preciso realizar una interpolación entre épocas para poder obtener la posición de los mismos cada 30 segundos. Para ello, la mejor forma es aplicando el polinomio de Lagrange, el cual nos dará una solución muy aproximada de la órbita del satélite, donde la incógnita del polinomio será la variable tiempo, la cual se podrá sustituir para obtener la posición del satélite en el momento deseado.

          Por tanto, asumimos valores funcionales (tj) que definen las épocas:

          tj,j=0,1,…,n  Entonces:

Es la definición de las funciones de correspondencia base (l_j (t)) de grado “n” relacionado con una época “t” arbitraria.

El valor funcional interpolado en la época “t” sigue el sumatorio:

Ésta es la definición matemática del polinomio. En la siguiente figura se puede apreciar un ejemplo del comportamiento del polinomio. En las partes extremas el polinomio aparece sesgado, pero en la parte central el polinomio se aproxima bastante a la función dada. También se puede observar que el polinomio siempre coincide en los puntos utilizados para la definición  del mismo f(tj).

Con todo esto, se va a generar un código que calcule la posición de un receptor con ambas soluciones (efemérides transmitidas y precisas), y que, a su vez, calcule la posición de los satélites que intervienen en la observación.

Metodología y resultados

Este apartado se va a dividir en dos grandes bloques:

1. Desarrollo de la interpolación de Lagrange: Se desarrollan tres scripts con el fin de  conseguir los resultados deseados, al final, se consiguió generar una función que calculara la interpolación para los satélites en una época determinada.

La función queda adaptada de forma que cada vez que se llama a ésta, devolverá solo la interpolación de ese satélite en la época indicada. Queda, por tanto, definida como:

function [pos_int] = Inter_Lgrge_Easy3(satelite,t_GPS,coorsat,tmp)

Donde las variables de entrada son:

• satelite: Es el número del satélite a interpolar.
• t_GPS: Es el tiempo de emisión de la señal por el satélite, es decir, la época en la que se desea obtener la posición del satélite. 
• coorsat: Matriz de coordenadas del satélite. Está compuesta por las coordenadas del satélite del cual se desea obtener la interpolación, que se almacenan fila por fila por cada una de las épocas. Proceden del archivo “MatrizCoorSat2.txt”.
• mp: Matriz de las épocas de interpolación. Proceden del archivo “Epocas2.txt”.

Devuelve:

• pos_int: Vector columna con las coordenadas del satélite interpolado en la época definida.

Con esta función se podrá obtener la posición de cada uno de los satélites dependiendo de la observación dada en el fichero de observables, en la época de salida de la señal del satélite.

2. Modificación de la rutina de “easy3”: Se va a explicar cómo se han modificado los distintos scripts que forman parte de la rutina para que funcionen con ficheros RINEX de la actualidad. También se ha modificado el programa para que calcule las coordenadas de una estación en distintas épocas y en distinto periodo de tiempo al programado inicialmente por Kai Borre, con ambas soluciones (efemérides transmitidas y precisas, estas últimas necesitan de la interpolación).

Se va a modificar toda la rutina, tanto el programa principal como las funciones asociadas al mismo, con el fin de poder calcular las coordenadas del receptor con las dos soluciones en un intervalo de observación de 2 horas, con observables cada 30 segundos, un total de 240 épocas.

El procedimiento de cálculo es el mismo en ambas soluciones, tanto en efemérides transmitidas como en efemérides precisas, cuya única diferencia es la obtención de la posición de los satélites, la cual se ha explicado en apartados anteriores. El proceso es el siguiente.

1. Preparación de los datos: Se preparan los datos necesarios para el cálculo. Para ello se emplean las funciones mencionadas en el apartado anterior.
2. Una vez preparados los datos éstos entran en la función que realiza el cálculo, “recpo_lsa.m” para efemérides transmitidas y “recpo_efp.m” para efemérides transmitidas. Como la obtención de la posición de los satélites ya se ha explicado, se va a generalizar el cálculo que realizan estas dos funciones, pues es el mismo procedimiento.
3. La forma de cálculo es simple, por cada época de observación se va a generar un sistema de ecuaciones que tenga como incógnitas las coordenadas de la estación y el estado del reloj. En cada una de las épocas aparece un número aleatorio de satélites observados, los cuales van a generar las ecuaciones del sistema, es decir, por cada satélite se añadirá una ecuación nueva.El sistema de ecuaciones se resolverá siguiendo un ajuste pos mínimos cuadrados, el cual devuelve una posición aproximada del receptor y el estado del reloj del mismo. Para la obtención de la matriz de coeficientes (matriz A) se necesitan los valores de la posición del satélite, el observable de la pseudodistancia y una posición previa del receptor. Lo mismo se necesita para obtener el vector de observables (vector K), solo que en este caso se añadirán las correcciones del error troposférico e ionosférico.Como se ha comentado, es necesario disponer de una posición previa del receptor, pero como al principio ésta se desconoce el procedimiento seguido es simple: Se calcula una primera iteración donde se parte que la posición del receptor es 0, y se calculan unas coordenadas previas, las cuales tendrán mucho error. Las nuevas coordenadas se almacenan en el vector de posición, y se vuelve a iterar realizan el mismo cálculo pero en este caso se empleará la nueva posición obtenida. Este proceso se repite tantas veces hasta que el error se reduce, es decir, se itera el proceso hasta obtener unas matrices estables y por tanto una posición del receptor buena. En la rutina de esay3 se itera hasta seis veces, pues son suficientes para obtener una posición buena.
4.  Con todas las posiciones obtenidas por cada época observada (concretamente 240 posiciones calculadas), se realiza una media aritmética de todas ellas obteniendo unas coordenadas promediadas de la estación.
5. Las coordenadas obtenidas se encuentran el sistema de coordenadas ECEF. Para una mejor comparación con la solución oficial, se ha realizado una transformación de las mismas a coordenadas geográficas y posteriormente a coordenadas UTM
6. El programa no solo calcula la posición del receptor, sino que inicialmente ya venía programado un ploteo de la variación de las coordenadas en las distintas épocas de cálculo. Además se ha añadido un ploteo de las órbitas de los satélites y el cálculo de la diferencia que se produce entre ambas soluciones, en las épocas de trabajo, con el fin de poder ver la diferencia entre calcular las coordenadas de los satélites de una forma o de otra.

Coordenadas

Posición de los satélites

Conclusiones

Con todo ello, se puede concluir que:

- Para la interpolación de Lagrange es necesario el uso de algoritmos y scripts de programación, ya que la complejidad del cálculo solo permite realizarlo por medio de la computación.

- La diferencia de cálculo de la posición de los satélites con ambas soluciones es de unos metros, por tanto, se puede concluir que ambas metodologías de cálculo son correctas, pues respetan el margen de error que las instituciones indican.

- Las coordenadas del receptor obtenidas, con ambas soluciones, se han mejorado con la aplicación de esta metodología, con respecto al cálculo que realizan algunos dispositivos a partir de la pseudodistancia, además destacar la mejora de unos 30-40 cm que las efemérides precisas generan con respecto a las efemérides transmitidas.

                                                                           ------

¿Que os ha parecido el trabajo de esta semana? Si tenéis cualquier duda o algún trabajo que queráis compartir no dudes en ponerte en contacto conmigo.

Dentro de poco tendréis en el blog los enlaces a los trabajos completos que vamos compartiendo, por si queréis profundizar un poco más.

Y nos vemos el domingo que viene con una entrada especial y algunas novedades (¡espero que tengáis tantas ganas como yo de conocerlas!)

¡Hasta la próxima!