1. Educación

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