Interpolación de las coordenadas de los satélites GPS para el posicionamiento geodésico I.

Resumen. Los datos de las efemérides GPS  que contienen las posiciones de los satélites G.P.S. (coordenadas x, y, z) , se obtienen en intervalos de 5 ó 15 minutos, sin embargo, cuando se realizan observaciones satelitales se requiere que los datos de las efemérides estén menos espaciados, por lo tanto, es necesario realizar interpolaciones para obtener las coordenadas de los satélites, para la hora de observación. En esta disertación, se discutirán los métodos de interpolación conocidos como el polinomio de Lagrange, y el método de Neville.

II.

Introducción.

Antes de iniciar con el tema central de esta ponencia, se recordaran algunos conceptos relacionados con el movimiento de los satélites de acuerdo al problema de los dos cuerpos en Mecánica Celeste. El problema de los dos cuerpos establece que: " Dados dos puntos de masas 'M' y 'm' separados una distancia 'r' y que se atraen según la ley de la Gravitación Universal:

donde G es la constante de Gravitación Universal, encontrar la trayectoria del desplazamiento de uno de los puntos con respecto al otro ". De acuerdo al estudio del potencial de un cuerpo, se tiene que el potencial de un cuerpo con simetría esférica es equivalente al potencial de un punto situado en su

centro y de masa igual a la masa del cuerpo considerado, por lo tanto, si se asume en primera aproximación que la tierra es una esfera, entonces, se puede considerar a la tierra como un punto. Sabemos que la componente principal del potencial terrestre en adición al potencial en 1/r es del orden de una milésima del término principal. Por lo tanto, si se desprecian estos términos suplementarios, se puede reducir el movimiento de un satélite a la solución del movimiento del problema de los dos cuerpos, que constituye una primera aproximación del movimiento real. Despreciando la masa m del satélite con respecto a la masa M de la tierra, se tiene que el centro de gravedad de la tierra es el centro de gravedad del sistema tierra – satélite. Despreciando el movimiento de traslación y suponiendo que el sistema de coordenadas inercial tiene su centro en el centro de gravedad de la tierra (el error que se tiene por esta hipótesis, se corrige considerando las perturbaciones lunisolares), entonces, esto nos conduce a estudiar el movimiento de un satélite "s", que es atraído por un punto fijo "o", situado en el origen de un sistema de coordenadas inerciales oxyz (fig. 1), que produce sobre el satélite "s" una aceleración:

Fig.1 Sistema de coordenadas inerciales o, x, y, z

Donde,  = G m es la constante geocéntrica de la gravitación,  es la ascención recta geocéntrica de s, y  es la declinación geocéntrica de s. Relacionando las expresiones (1) y (2), se tiene la siguiente ecuación de movimiento del satélite en forma vectorial:

La ecuación (3) es la forma vectorial de una ecuación diferencial de segundo orden con seis constantes de integración, que son los seis parámetros orbitales de Kepler. La posición de los satélites en el espacio, se pueden analizar en primera aproximación, considerando las orbitas normales, en un plano orbital que permanece fijo en el espacio, donde los seis elementos de Kepler juegan un papel importante.

Fig. 2. Plano orbital en la esfera de direcciones

Los seis elementos de Kepler son:

a: semieje mayor de la elipse orbital e: excentricidad de la elipse orbital

i: Inclinación de la orbita : Ascención recta del nodo ascendente : Argumento del perigeo v: Anomalía verdadera

La anomalía verdadera "v" es el único elemento de Kepler que es función del tiempo para orbitas normales, los cinco elementos restantes permanecen constantes. La ecuación de movimiento del satélite, se puede expresar como:

Realizando el producto escalar de la ecuación (4) con el vector de la velocidad:

se



tiene

que:

pero

y 

La integración de esta última expresión, es la integral de la energía: El primer

término del lado derecho de

la

ecuación (8) es la energía

cinética del satélite con masa unitaria (m= 1), y el segundo término su energía potencial, por lo que la energía mecánica total del movimiento del satélite permanece constante, esto implica que la energía se conserva cuando no existen fuerzas externas. Ahora consideremos un sistema de coordenadas polares en el plano de la órbita, el polo es el origen de coordenadas del sistema tridimensional y el eje polar una dirección arbitraria en el plano de la órbita.

Fig. 3. Componentes de la velocidad

La velocidad "V" del satélite, se descompone en una componente radial "Vs", y en una componente normal al radio vector r, "Vn" (Fig. 3):

con

Como la fuerza aplicada pasa por el origen (fuerza central), su momento con respecto al origen es cero (constante).

donde "C" es la constante de las áreas. 

El área (dA) descrita por el radio vector "r" en el intervalo de tiempo dt es:

Aplicando el teorema de la conservación de la energía en su forma diferencial:

y ya que OS y F son colineales,

de donde integrando, se tiene que:

 Donde "h" es la constante de la ecuación de la energía. Considérese el siguiente diagrama:

 

que es la ecuación diferencial de las trayectorias posibles. Reescribiendo la expresión anterior:

sea:

La ecuación (17) se transforma en:

Esta ecuación diferencial, se resuelve con respecto a "" y se integra: La ecuación se resuelve por variables separables

Integrando:

Esta es la ecuación de una cónica de foco "0" Haciendo:

Por lo tanto la ecuación (18), resulta:

"e" es la excentricidad, "p" es el parámetro, y "a" es el semieje mayor, 0 es el ángulo polar en la dirección del vértice más próximo al foco "0" (punto llamado el perigeo), el ángulo "v" contado a partir de esta dirección es llamado la anomalía verdadera. De la ecuación (20), se deduce que:

De donde:

Y por lo tanto, la integral de la energía (ecuación (16)), se expresa por:

Analizando la ecuación (22), se puede observar que depende del signo de

el género de la cónica

h.

a) Si

h = 0, la cónica es una parábola

b) Si

h  0, la cónica es una hipérbola

c) Si

h  0, la cónica es una elipse

Estudio del movimiento elíptico. El satélite se desplaza alrededor de la tierra en una órbita elíptica, con la tierra en uno de los focos, la elipse se determina por "a", "e", y "S" según la fig. 4, los primeros dos elementos "a", "e", definen el tamaño de la elipse, mientras que "S" define la posición del satélite en el plano de la órbita. El movimiento medio y el tiempo definen a "S". M = n (t – t0)

Fig. 4 Sistemas de coordenadas en el plano orbital Donde; a: semieje mayor de la órbita satelital b: semieje menor de la órbita satelital F: foco de la órbita satelital donde se localiza la Tierra (geocentro) ae: distancia focal de la órbita satelital

q1, q2: sistema de coordenadas geocéntricas x1, x2: sistema de coordenadas cartesianas de la órbita v: anomalía verdadera E: anomalía excéntrica r: radio vector del satélite medido desde el geocentro perigeo: punto más cercano a la tierra del satélite apogeo: punto más alejado a la tierra del satélite s: posición del satélite en su órbita s': proyección de la posición del satélite en un círculo de radio "a" que es tangente a la órbita del satélite en el perigeo y el apogeo. órbita: trayectoria que describe el satélite en el espacio

La posición del satélite en cualquier tiempo "t", se define por sus coordenadas polares (r, v). El radio vector de posición del satélite "r" puede calcularse utilizando las coordenadas (q1, q2), que como se puede observar en la fig.4, se describen por las siguientes ecuaciones:

El paso siguiente es la transformación de las coordenadas

q1, q2

a las

coordenadas inerciales X, Y, Z, mediante tres rotaciones al dextrorso. La primera rotación se hace en el plano de la elipse girando al dextrorso el sistema

q1, q2

desde el perigeo hasta el nodo ascendente (fig. 5), esto se hace alrededor del eje Z, y se representa por R3 (– )

Fig. 5. Movimiento orbital del satélite

q q q

La siguiente rotación involucra el giro de las coordenadas 1, 2, 3 por un ángulo

i

de la inclinación de la órbita en una dirección al dextrorso alrededor del eje X, y por

i

lo tanto se indica por R1 (– ). La tercera rotación es del ángulo  en la dirección al dextrorso alrededor del eje Z, por lo que se indica por R3 (– ). 

Esto define las coordenadas cartesianas inerciales del satélite en el espacio en términos de los seis elementos Keplerianos.

III.

Desarrollo

La interpolación polinomial es uno de los problemas fundamentales en análisis numérico, y consiste básicamente en seleccionar una función p(x), a partir de diversas clases de funciones, de tal manera que la gráfica de y = p(x), pase a través de los datos conocidos (xi, yi), i = 1,…….., n, los puntos pueden representar mediciones de un problema físico, o pueden ser obtenidos de una función conocida, en nuestro caso, se desconoce la función, y se pretende obtener un

polinomio que nos permita interpolar datos conocidos, que son las coordenadas precisas de los satélites G.P.S.. Como se comentó en un principio, las efemérides precisas de los satélites, son proporcionadas de manera gratuita por diversas instituciones como por ejemplo el N.G.S. de los Estados Unidos de Norteamérica. En esta presentación se consideran dos métodos de interpolación, el Polinomio de Lagrange y el Método de Neville. 

Polinomio de Lagrange

La expresión general para el polinomio de interpolación de Lagrange está dada por:

con

La variable x, se encuentra solo en el numerador de cada término y los denominadores son números, entonces:

i

n

i

i

Así L es un polinomio de grado , nótese que cuando L (x) es evaluada en x = x ,

i

cada factor en la ecuación precedente es 1, pero cuando L (x) es evaluada en

j

i j

algún otro nodo x , uno de los factores es cero, y así, L (x ) = 0 para

i

j

 ,

entonces es posible interpolar cualquier función "f " por el polinomio de interpolación de Lagrange:

Entonces, un polinomio de Lagrange de orden 1, se representa por la siguiente expresión:

Un polinomio de Lagrange de orden 2, se representa por la siguiente expresión:

y así sucesivamente. Con estos dos ejemplos, podemos observar que el orden del polinomio es de (n – 1), donde n es el número de datos. Analizando la forma del polinomio de Lagrange, se observa que los datos que se van a interpolar, se deben considerar de manera separada, es decir, el valor de x

en el polinomio p(x), representa la coordenada x, o la coordenada y o la coordenada z del satélite. Así que, los polinomios de interpolación de Lagrange, estarían dados por: p(x), p(y) y p(z). Un ejemplo de la aplicación del polinomio de Lagrange para la interpolación de las coordenadas de los satélites G.P.S., es el siguiente: Sean los siguientes datos para las coordenadas precisas de los satélites: * P

2014 5 6 0 0 0.00000000 1 13355.037945 -18951.987446

12803.387070

11.208738

* P

2014 5 6 0 5 0.00000000 1 13320.680631 -18442.871889

13562.550679

11.209122

* 2014 5 6 0 10 0.00000000 P 1 13289.875785 -17905.549153

14295.528916

11.209572

* 2014 5 6 0 15 0.00000000 P 1 13263.885400 -17340.988109

15000.911043

11.209979

* 2014 5 6 0 20 0.00000000 P 1 13243.920415 -16750.251476

15677.340498

11.210358

* 2014 5 6 0 25 0.00000000 P 1 13231.134938 -16134.491820

16323.517545

11.210780

* 2014 5 6 0 30 0.00000000 P 1 13226.620727 -15494.947179

16938.201788

11.211149

En la primera línea de la tabla anterior está el año, mes, día, horas, minutos y segundos de inicio. En la segunda línea de la tabla anterior esta una letra "p" que indica posición, el número del satélite, las coordenadas del satélite en kilómetros, y el registro del reloj en microsegundos. En el polinomio de Lagrange, se consideran los datos de la tabla anterior, excepto los datos correspondientes a las 0

h

15

min

0

seg

, ya que este dato se va a calcular

con el polinomio y se comparará con el dato correspondiente a la tabla, con el propósito de observar las diferencias obtenidas al utilizar el polinomio desarrollado. Entonces el orden del polinomio es de orden 5, ya que se consideran 6 series de datos. La tabla de datos es la siguiente:

i

xi

fi

0

0

13355.037945

1

5

13320.680631

2

10

13289.875785

3

20

13243.920415

4

25

13231.134938

5

30

13226.620727

La columna "i" son los datos consecutivos. La columna "xi" son los tiempos de observación en minutos La columna "fi" son las coordenadas "x" para cada hora de observación en kilómetros. Para la interpolación de las coordenadas "y" y "z", se procede de manera similar que para la interpolación en la coordenada "x". Entonces el polinomio de interpolación de Lagrange es el siguiente:





Método de Neville.

El método de Neville está basado en el siguiente teorema: Teorema. Sea "f" definida en los (k+1) puntos x0, x1,…., xk y sean xi y xj dos puntos distintos de este conjunto. Sean P0, 1,…, coinciden con "f" en x0, x1,…, xi

i – 1, i+1,…, k(x)

– 1,

el polinomio de Lagrange que

xi+1,…, xk (el punto xi es el único que no se

encuentra en esta lista). Similarmente, sea P0, 1,…,

j – 1, j + 1, k(x)

el polinomio de

Lagrange que coincide con "f" en x0, x1,…, xj – 1, xj+1,…, xk. Por supuesto que P0, 1,…, i – 1, i+1,…, k(x)

y P0, 1,…,

j – 1, j + 1, k(x)

son polinomios de grado (k – 1).

Entonces el polinomio de Lagrange P0, 1,…, k(x) a través de los (k + 1) puntos: x0, x1,…, xk puede calcularse con la siguiente expresión: La idea del Método de Neville es usar recursivamente los polinomios de Lagrange de potencias bajas, a fin de calcular polinomios de Lagrange de potencias más altas. Entonces: P0,1 es el polinomio interpolante que pasa por (x0, f(x0)), (x1, f(x1) P1,2 es el polinomio interpolante que pasa por (x1, f(x1)), (x2, f(x2)) . . . P0,1,2 es el polinomio interpolante que pasa por (x0, f(x0)), (x1, f(x1)), (x2, f(x2)) y así sucesivamente. Una de las características del método de Neville, es que se basa en la relación de recurrencia derivada de las diferencias divididas para obtener los polinomios de interpolación. El método de Neville se aplica cuando se quiere interpolar f(x) en un punto x = p para polinomios de interpolación de Lagrange de muy alto orden.

Por ejemplo, sean tres puntos distintos x0, x1, x2 en los que se puede evaluar f(x), es decir, f(x0), f(x1), f(x2), a partir de estos tres puntos, podemos construir un polinomio de orden cero (constante) que se aproxime a f(p). Los polinomios de Lagrange de primer orden son:

Observamos que f(xi) = Pi(x) 

Y similarmente

Entonces el polinomio de tercer orden está dado por:

Que es justamente el polinomio de Lagrange de tercer orden que interpola los puntos x0, x1, y x2. Ahora se aplica el método de Neville para resolver el mismo problema del inciso anterior.

Se busca calcular el valor de la coordenada x del satélite, que corresponde a las 0 horas 15 min 0 seg, a través del método de Neville.

Resumen de resultados:

Método de interpolación

Coordenada

x (en metros), para la

observación de las 0 hrs. 15 min 0 seg.

IV.

Efemérides precisas

13263885.4000

Polinomio de Lagrange

13263885.4129

Método de Neville

13263885.4122

Conclusiones.

Se puede observar que la diferencia para la coordenada x del satélite, entre las efemérides precisas publicadas por el National Geodetic Survey y las encontradas por medio de los dos métodos de interpolación, es de aproximadamente 1 centímetro, lo que nos indica que ambos métodos de interpolación dan resultados confiables, y las diferencias podrían minimizarse, aumentando el número de datos, lo que nos resultaría en polinomios de mayor grado. V.

Bibliografía



Fowles, G.R., Analytical Mechanics, Saunders College Publishing, Usa, 1993.



Torge, w., Geodesy, Walter De Gruyter Inc., U.S.A., 1991



Rapp, R., Advanced Theoretical Geodesy, Notas de Clase de posgrado, Fort Clayton, Panama.



Mueller, Ivan, Satellite Geodesy, Notas de Clase de posgrado, Fort Clayton, Panama, 1985.



Soler, T., Teoría Y Aplicaciones de Levantamientos con Clase, Aguascalientes, Ags., 1994.



Levallois, Mécanique Céleste



Chen Berlin, Polynomial Interpolation



U.A.A., Notas de clase de Análisis Numérico, 1999



Leykekhman, Dmitriy, Polynomial Interpolation



Karris, Steven T., Numerical Analysis

G.P.S., Notas de