Cet ouvrage fait partie de la bibliothèque YouScribe
Obtenez un accès à la bibliothèque pour le lire en ligne
En savoir plus

Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores

De
9 pages
Resumen
Se presenta un modelo automático de corrección basado en la resolución orbital del satélite y la aplicación de un algoritmo que utiliza los ángulos de orientación para el cálculo de la posición geográfica del píxel. El modelo concluye con la implementación de remuestreo por el procedimiento del vecino más próximo. El modelo ha sido validado en un escenario completo de imágenes del sensor NOAA-AVHRR, especificando el análisis de los errores encontrados en unidades de píxel.
Abstract
A procedure of orbital resolution of the spacecraft position has been developed, together with the subsequent application of an algorithm using the orientation angles for the calculation of the geo-graphical pixel position. The model concludes with the implementation of the resampling algorithm of the images corrected through the nearest neighbour method. The model applied has been validated on a complete scenario of the NOAA-AVHRR sensor imagery, specifying the errors made in pixel units.
Voir plus Voir moins

Revista de Teledetección. 2000
Algoritmo automático de corrección geográfica
para imágenes NOAA-AVHRR: desarrollo y
análisis de errores
A. Calle, J. L. Casanova y J. Sanz
Correo electrónico: abel@latuv.uva.es
Laboratorio de Teledetección de la Universidad de Valladolid


RESUMEN ABSTRACT
Se presenta un modelo automático de corrección A procedure of orbital resolution of the spacecraft
basado en la resolución orbital del satélite y la position has been developed, together with the
aplicación de un algoritmo que utiliza los ángulos subsequent application of an algorithm using the
de orientación para el cálculo de la posición geo- orientation angles for the calculation of the geo-
gráfica del píxel. El modelo concluye con la im- graphical pixel position. The model concludes with
plementación de remuestreo por el procedimiento the implementation of the resampling algorithm of
del vecino más próximo. El modelo ha sido valida- the images corrected through the nearest neighbour
do en un escenario completo de imágenes del method. The model applied has been validated on a
sensor NOAA-AVHRR, especificando el análisis complete scenario of the NOAA-AVHRR sensor
de los errores encontrados en unidades de píxel. imagery, specifying the errors made in pixel units.


PALABRAS CLAVE: NOAA-AVHRR, correc- KEY WORDS: NOAA-AVHRR, automatic geo-
ción geográfica automática. graphical correction.


automáticamente para generar polinomios de ajus-INTRODUCCIÓN
te (Illera et al., 1996). La técnica de determinación
La corrección geográfica de las imágenes pre- de los parámetros de cruce con el ecuador ha sido
senta, en todos los casos estudiados en teledetec- empleada, también, por (Moreno and Meliá, 1993)
ción, grandes dificultades de cara a la obtención de para generar un modelo de corrección geométrica
resultados precisos. Existen múltiples fuentes de teniendo en cuenta varios factores de barrido del
error difícilmente tenidas en cuenta en la aplica- sensor que llevan a la obtención de resultados
ción de modelos teóricos funcionando de forma precisos; sin embargo, los métodos mencionados
totalmente automática, sin la interacción del usua- que aplican esta técnica emplean puntos de control
rio. Las imágenes corregidas, generadas por estos obtenidos de forma interactiva para el cálculo de
modelos, deben ser ajustadas finalmente mediante estas magnitudes, lo que demora de forma impor-
la identificación de puntos de control y la aplica- tante el procesado de la imagen.
ción de polinomios de ajuste, con la finalidad de En lo que respecta a los modelos automáticos,
obtener la máxima precisión. Sin embargo, el au- puede verse un planteamiento muy general en
tomatismo en todas las tareas de procesamiento de (Rosborough et al., 1994) En el que desarrolla una
imágenes debe ser considerado como un objetivo metodología para la corrección conociendo los
primordial en aquéllos casos que requieran la ge- ángulos estándar roll, pich y yaw, de orientación
neración de productos en tiempo real. Esto conlle- del satélite; además estos autores incluyen un pro-
va la búsqueda de un compromiso entre precisión cedimiento de control de altitud, a partir de un
exigida y tiempo de proceso. punto de control, que mejora en gran medida los
En el caso de la corrección del sensor AVHRR, resultados procedentes de la metodología automá-
los modelos más ampliamente utilizados toman al tica
menos un punto de control para resolver las condi- En lo que sigue, presentaremos un modelo basa-
ciones de cruce con el ecuador de la órbita ascen- do en Brunel and Marsouin, 1986, que aplica la
dente (Ho and Asem, 1986, Bachman and Bendix, metodología de Rosborough et al.,1994 calculando
1992). Estos procedimientos parten del conoci- los ángulos de orientación a partir de los paráme-
miento de los parámetros de cruce con el ecuador tros orbitales. Finalmente, el modelo concluye con
de la última órbita ascendente previa a la captura la implementación del algoritmo de remuestreo de
de la imagen (longitud geográfica y tiempo de las imágenes corregidas por el procedimiento del
cruce). En otros casos, se ha utilizado como mode- vecino más próximo.
lo de base el mencionado anteriormente, incluyen- La validación del modelo aplicado se ha realiza-
do más puntos de control de tierra, para una vez do sobre un escenario completo de imágenes del
resuelta la órbita, utilizar dichos puntos de control sensor NOAA-AVHRR, especificando los errores
Nº 14 – Diciembre 2000 1 de 9 A. Calle, J. L. Casanova y J. Sanz
cometidos en unidades de pixeles. Para todo el 3. Argumento de Perigeo, ω : Es el ángulo, me-0
desarrollo del estudio se ha elaborado un software dido en la elipse orbital, entre el cruce con el
específico, en entorno Windows. Finalmente se ecuador terrestre y el perigeo de la órbita.
discuten algunas dificultades encontradas en la 4. Anomalía Media, M : Es el desplazamiento 0
aplicación de la metodología de verificación con- angular recorrido por el satélite, si éste
cernientes a errores técnicos del sensor utilizado hubiera estado moviéndose con velocidad
(Brunel and Marsouin, 2000). Este modelo, así angular constante en una órbita circular.
planteado, genera unos buenos resultados y es de 5. Ascensión Recta del Nodo Ascendente, Ω : Es 0
fácil programación para ser incluido en programas el ángulo de longitud entre el punto Aries o
de procesamiento general. punto vernal y el punto en que el satélite cru-
za el ecuador terrestre en una órbita ascen-
dente (sur-norte). DESCRIPCIÓN TEÓRICA DEL
6. Inclinación de la Órbita, i: Es el ángulo entre MODELO
el plano de la elipse y el plano ecuatorial te-
El modelo de corrección consta de tres fases rrestre.
principales, cuya descripción teórica se expone en 7. Movimiento Medio, n: Es el promedio de la
los siguientes sub-apartados. En primer lugar se velocidad angular del satélite en la órbita,
realizará el cálculo de la posición del satélite a expresado en vueltas/día.
partir de los parámetros orbitales; posteriormente
se aplicarán las funciones de transformación (fi- De los elementos descritos, tres de ellos sirven
la,columna) →(long,lat), basadas en matrices de para localizar el satélite en la elipse orbital (e, M,
transformación y finalmente se realizará el proceso n) y otros tres para encontrar la orientación de la
de remuestreo de la imagen. En esta última etapa elipse respecto a una posición terrestre ( ω, Ω, i).
se expondrá la metodología de la corrección inver-
sa. Para calcular la posición del satélite, pueden se-
guirse varias metodologías. Una de ellas consiste
Determinación de la posición del satélite en la aplicación del modelo de propagación
NORAD-SGP4, sobre todo si se utilizan como
A partir de los parámetros orbitales podrá deter- parámetros orbitales los aportados por NORAD-
minarse la posición del satélite en un instante ante- two-line-elements. SGP4 produce como resultado
rior o posterior aplicando las ecuaciones que rigen las coordenadas cartesianas del satélite en un sis-
el movimiento en la órbita. El tiempo de validez de tema de referencia inercial, con uno de los ejes del
dichas magnitudes no es muy elevado en el caso de sistema posicionado en el punto aries o equinocio
los satélites polares que tienen una órbita de baja vernal en la época 1950. Deberá aplicarse una
altitud (en torno a 850 km.), debido a que la órbita matriz de triple rotación para generar estas coorde-
está modificándose constantemente por causa de nadas sobre un sistema en la época actual de cálcu-
las fuerzas de perturbación (Brunel and Marsouin, lo o la época 2000 y finalmente realizar un cálculo
1986): fuerzas de gravitación causadas por irregu- del tiempo sideral del meridiano de Greenwich
laridades en la distribución de la masa de la tierra, para obtener XYZ en un sistema en rotación soli-
fricción atmosférica (hasta 1000 km. de altitud), dario con la tierra. Posteriormente, la latitud y
presión de la radiación solar, fuerzas de gravita- longitud geográficas serán obtenidas aplicando las
ción del Sol y la Luna, etc. De todas ellas, la más transformaciones geodésicas al modelo de elipsoi-
importante en un gran orden de magnitud, es la de terrestre vigente. La otra alternativa es la reso-
debida al achatamiento de la Tierra. lución de la ecuación de Kepler y la aplicación
Los datos orbitales, para un instante inicial de- posterior de las ecuaciones de perturbación para la
nominado Época, t , son proporcionados por varios 0 precesión del plano de la órbita y el argumento de
tipos de boletines electrónicos difundidos a través perigeo. Si existe dificultad para encontrar el códi-
de la red Internet. En nuestro modelo se ha utiliza- go del SGP4 (originalmente en FORTRAN) para
do el NORAD “Two Line Element” (TLE), que insertarlo en los programas de procesamiento,
recibe este nombre porque resume la información presentamos una metodología alternativa que apor-
de cada satélite en dos líneas. En este boletín figu- ta unos resultados perfectamente comparables a los
ran los elementos de localización necesarios, que del modelo mencionado, siempre que los datos
son los siguientes: orbitales tengan antigüedades no superiores a los
tres días. Esto no supone ningún problema dado
1. Época, t : La fecha (con exactitud de milise-0 que es muy fácil encontrar estos ficheros actuali-
gundos) para la cual son emitidos los ele- zados para cada día en Internet.
mentos orbitales pertenecientes a ese TLE.
2. Excentricidad, e: Es la excentricidad de la El cálculo de la posición del satélite se realizará
elipse que traza el satélite en su movimiento en dos procesos: en primer lugar se calculará la
alrededor de la Tierra. posición en el plano orbital y posteriormente se
2 de 9 Nº 14 – Diciembre 2000 Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores
realizará el cambio del sistema de referencia para Estas coordenadas cartesianas corresponden a un
la determinación de las coordenadas geográficas de sistema de referencia cuyo centro se encuentra en
la posición del satélite.Localizar el satélite en su el foco de la elipse; o lo que es lo mismo, en el
propia elipse orbital consiste en resolver la ecua- centro de la Tierra, y la parte positiva del eje X es
ción de Kepler. Dicha ecuación se presenta en la la línea que une el mencionado centro con el peri-
forma de la denominada ecuación de anomalías, geo del satélite.
cuya expresión es: Las coordenadas cartesianas calculadas en el
plano de la órbita deberán ser transformadas a un
sistema de ejes terrestre, en reposo respecto a la  
  1 +e E  −1   Tierra, para obtener la localización del satélite M = E − e ·sin (E) T = 2 · tan · tan  1 −e 2 respecto a un punto de la superficie de la misma.   
  Para ello serán utilizados los parámetros orbitales
de orientación de la elipse, definidos con anteriori-
en la que e es la excentricidad, M es la anomalía dad. Sin embargo dos de ellos, el argumento de
media en el instante de cálculo y E es la anomalía perigeo y la ascensión recta del nodo ascendente,
excéntrica; que es un término que relaciona la sufren una variación constante provocada por una
distancia radial entre el satélite y el foco. T es la velocidad de rotación (también denominada prece-
anomalía verdadera, T, es el ángulo real barrido sión), de manera que es necesario recalcular sus
por el satélite en el momento exacto del cálculo, valores para el instante de análisis a partir de los
medido desde el perigeo en el sentido antihorario. valores aportados por el TLE para el instante de la
La anomalía media , M, se obtendrá a partir de la época. El valor instantáneo de la ascensión recta
anomalía media inicial M mediante el movimiento 0 del nodo ascendente, Ω, y el valor instantáneo del
medio, n, de la forma: M=M +n(t-t ). 0 0 argumento de perigeo, ω, vienen dados por la
expresión:
La resolución de la ecuación de Kepler requiere
la utilización de técnicas de cálculo numérico 3.5
 R cos(i) E  degmediante aproximaciones sucesivas, ya que se trata Ω = Ω +P ·(t −t ) ; P = −9.97· ·0 Ω 0 Ω 2 2   daysmm (1 −e )de una ecuación implícita. El procedimiento utili-  
zado se basa en el método de Newton-Raphson 3.5
2 R 5·cos (i) −1 deg(Crawford, 1995). A partir de las funciones:  E ω = ω +P ·(t −t ) ; P =4.98· ·0 ω 0 ω 2 2  diasmm (1 −e )  
f (E) =M +e • sin( E) −E =0f (E )N E =E −N +1 N con t el instante de cálculo, t la fecha de la época 0f '(E ) f '(E) =e • cos( E) −1N
de los parámetros orbitales, y R el radio ecuatorial E
terrestre. En el caso del satélite NOAA, la veloci-
en la que E es la estimación de E en la iterac-N dad de precesión del nodo ascendente P =d Ω/dt Ω
ción N-ésima. La solución quedará determinada tiene un valor cercano a 0.99º/día, por tratarse de
cuando se cumpla la condición: E -E < δ, N N-1 una órbita heliosíncrona.
asignando a δ un valor muy pequeño, dependiendo
de la precisión exigida (p.e. del orden de magnitud Las coordenadas cartesianas, (X , Y , ZTierra Tierra Tie-
-1010 ). La obtención de la anomalía excéntrica ) respecto a un sistema terrestre con el eje z coin-rra
conduce al conocimiento de la anomalía verdadera cidente con el eje de rotación de la Tierra y en
en el instante de cálculo, tal como se ha menciona- reposo respecto de ésta, serán calculadas a partir
do. de las coordenadas del plano orbital (Xorbi-
,Y ,0), mediante una matriz gaussiana de ta orbita
Finalmente, las coordenadas de la posición del transformación de sistema de referencia, según:
satélite, X e Y, respecto de la elipse orbital se
     obtienen a partir del parámetro semieje mayor de x c( ω)c( Ω) −s( ω)c(i)s( Ω) −s( ω)c( Ω) −c( ω)c(i)s( Ω) s(i)s( Ω) x Tierra     Orbita 
     la elipse, smm, calculado mediante el movimiento y = c( ω)s( Ω) +s( ω)c(i)c( Ω) −s( ω)s( Ω) +c( ω)c(i)c( Ω) −s(i)c( Ω) • yTierra Orbitaa     
z s( ω)s(i) c( ω)s(i) c(i) 0     Tierramedio, n, a partir de la 3ª ley de Kepler, expresada      
2 3matemáticamente a partir de: µ=n ·smm , en la que
µ es la constante gravitacional terrestre cuyo valor
3 -2es 398601.3 km ·s . Las coordenadas cartesianas donde, por brevedad se ha utilizado la notación
X e Y del satélite, referidas a su órbita, pueden c( ϕ)=cos( ϕ) y s( ϕ)=sen( ϕ).
calcularse fácilmente mediante la anomalía excén-
trica: Para el cálculo final de la latitud geográfica, ϕ , sat
y longitud geográfica, θ , del satélite, es necesario sat
  realizar el paso de coordenadas geocéntricas a X = smm ·  cos(E) − e orbit   geográficas, y además, utilizar el tiempo sideral 2Y = smm · sin(E) · 1 −eorbit del meridiano de Greenwich, que es el ángulo
Z =0orbit
Nº 14 – Diciembre 2000 3 de 9 A. Calle, J. L. Casanova y J. Sanz
formado por dicho meridiano y el punto aries, en el
instante de cálculo, para introducir la rotación
terrestre:

 
   2   R z  y −1 E −1  ϕ = tan ; θ =tan −15·tsat   sat sideral,Greenwich2    2 2 xR P x +y        

Corrección directa basada en los ángulos de
orientación.

El proceso de geo-referenciación de un píxel M
de una imagen, consiste en encontrar la relación
existente entre las coordenadas de la imagen C , M
F (columna, fila) y las coordenadas geográficas M
θ , ϕ (longitud, latitud). Por supuesto, tanto el M M Figura 1: Esquema de los sistemas de coordenadas en la correc-
ción directa. problema directo como el inverso deben estar
resueltos, ya que en el proceso de remuestreo, el
algoritmo comienza por la imagen corregida o en la que γ es el ángulo del espejo del sensor para
mapa, completándose con los valores de los pixe- el píxel M, determinado a partir de la columna que
les que son buscados en la imagen original. Por M ocupa en la imagen (1024.5-C )* δ. I es el M i
ello, la representación matemática que nos caracte- ángulo de inclinación orbital, β es el ángulo for-
riza el problema directo y el inverso se reduce al mado, sobre la elipse orbital, entre el satélite y el
-1 -1conocimiento de las funciones f, g, f y g , de la plano ecuatorial terrestre y equivale a la suma del
forma: argumento de perigeo en el instante de cálculo
(obtenido a través de la expresión de propagación
−1θ = f()C , F C = f ( θ , ϕ ) mediante la fórmula de perturbación ya explicada M M M
−1 ω= ω +P (t-t )) y la anomalía excéntrica, E, prove-0 ω 0ϕ = gC , F F = g( θ , ϕ)MI M M
niente de la resolución de la ecuación de Kepler.

Para localizar el punto terrestre M, deberemos
λ es el ángulo formado por el nodo ascendente resolver la ecuación vectorial:
del satélite y el meridiano de Greenwich. Para
OM = OS + SM en la que OM es el vector calcularlo hemos introducido una modificación
que une el centro de la Tierra con la superficie en respecto a la metodología propuesta por Brunel
el píxel M, OS es el vector de posición del satélite and Marsouin, 1986.: Para ello utilizamos la longi-
en ese instante y SM es el vector con origen el tud del nodo ascendente instantánea (utilizando la
satélite y final en el punto M. Todos los vectores velocidad de precesión por perturbación) restando
referidos a un sistema de referencia solidario con el tiempo sideral del meridiano de Greenwich en el
la tierra en rotación. La posición del satélite queda instante de cálculo, de la forma:
determinada, como se ha descrito en el anterior
apartado, a partir de las coordenadas geográficas λ = ( Ω + P (t − t ) ) − t0 Ω 0 sideral,Greenwich
(Lat Lon ); sus coordenadas cartesianas, en sat, sat
forma vectorial quedan de la forma (ver figura 1): El parámetro D, es la distancia escalar entre el
satélite y el punto M sobre la superficie terrestre.
RPR = Se calcula imponiendo la condición de satisfacer la
2X R cos()ψ cos(Lon )     R sat sat 2 P  ecuación del elipsoide terrestre: cos ()ψ −1 +1    2 OS = Y = R cos()ψ sin(Lon ) ; consat sat R E     2   Z R sin()ψ R sat    Ptan( ψ ) = tan(Lat )sat X =X + D· VX2  M satR 2 2 2 E X Y Z M M M+ + =1; Y = Y + D·VY M sat 2 2 2R R RE E P  Z = Z + D·VZEn cuanto al vector de posición SM del píxel,  M sat
con origen en el satélite, su expresión se calcula a
partir de matrices de rotación que utilizan los án- Una vez obtenidas las coordenadas cartesianas
gulos de orientación orbitales del satélite; su ex- XYZ) , en el sistema solidario terrestre, las coor-M
presión matemática es de la forma: denadas geográficas Long,Lat) se obtienen fácil-M
mente a partir de las mismas expresiones utilizadas
VX   − cos( γ )[cos( β ) cos( λ) − cos(I)sin( β )sin( λ)] + sin( γ )sin(I)sin( λ)  para el caso del posicionamiento del satélite, des-
   
SM = D VY = D − cos( γ )[cos( β )sin( λ) + cos(I )sin( β ) cos( λ)] − sin( γ )sin(I ) cos( λ)    critas anteriormente. Las coordenadas UTM serán VZ   − cos( γ ) sin(I )sin( β ) + sin( γ ) cos(I)    

4 de 9 Nº 14 – Diciembre 2000 Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores
obtenidas utilizando las expresiones aportadas en proyección UTM huso 30 N, manteniendo cons-
el apéndice anexo. tante la resolución espacial de todos los pixeles a
Es evidente, en todos los procesos anteriores, 1.1 km, que es la mejor resolución obtenida en el
que el instante de cálculo está determinado a partir nadir de la órbita. Este proceso ha consistido en
del tiempo de la línea de imagen que se está corri- construir una nueva imagen, con los pixeles de la
giendo. imagen original, una vez que se ha realizado la
redistribución espacial de los mismos. Para ello, el
La corrección inversa. algoritmo computacional seguido, consiste, de
forma esquemática, en los siguientes pasos:
La corrección inversa, es decir, la obtención de
la fila y columna de imagen a partir de las coorde- 1. Situarse en una coordenada determinada del
nadas de un mapa, está basada en que el plano de mapa
barrido del sensor es perpendicular a la superficie 2. Calcular, mediante las funciones de transfor-
terrestre, lo que se traduce en resolver la ecuación mación obtenidas, la posición en la que se
siguiente (Brunel and Marsouin, 1986): encuentra el píxel de dicha coordenada en la
imagen original
sin(I)sin( ψ ) + cos(I) cos( ψ )sin(Lon − λ) r sin(V − E)sat
tan( β) − − = 0 3. Acceder al fichero original para la lectura del
cos( ψ ) cos(Lon − λ) R cos( β ) cos( ψ ) cos(Lon − λ)
sat sat mencionado píxel y escribirlo en la posición
de la coordenada donde los ángulos han sido ya explicados y V es la
anomalía verdadera calculada a partir de la ecua-
Sin embargo, la aplicación de este algoritmo no ción de Kepler, como se ha visto en anteriores
es una operación inmediata, debido a que la posi-apartados. R es el radio del elipsoide terrestre a la
ción en la que debe encontrarse el píxel es una latitud geográfica de cálculo y r es el módulo del
posición intermedia o mezclada de varios pixeles. radio vector que une el centro de la tierra con el
Para solucionar este problema se han tenido en satélite.
cuenta tres procedimientos: vecino más próximo, Para resolver esta ecuación hemos adoptado un
interpolación bilineal y convolución cúbica. Fi-procedimiento diferente al propuesto por Brunel
nalmente se ha optado por el procedimiento de and Marsouin. Para ello hemos aplicado un proce-
vecino más próximo siguiendo el criterio de no dimiento de bisección, tomando como instantes de
transformar la información radiométrica de origen. cálculo iterativos iniciales los extremos de la cap-
Las ecuaciones de cálculo para las coordenadas tura de la imagen. Dicha ecuación es siempre con-
UTM se encuentran en el apéndice al final del vergente entre los límites de tiempo mencionados.
trabajo. Para el instante de convergencia, serán calculados
todos los ángulos, y con ellos, las coordenadas del
RESULTADOS punto M referidas al sistema de situa-
do en el satélite, y que vienen dadas por: Se describen a continuación las fases de obten-
ción de resultados: descripción del escenario de
X = R [sin(I)sin( ψ )sin( β) + cos( ψ )cos( β)cos(Lon − λ) + cos(I)cos( ψ )sin( β)sin(Lon − λ) ] − r cos(V − E)M ,O sat sat prueba, resultados obtenidos y el análisis de error. 
Y = R [sin(I)sin( ψ )cos( β ) − cos( ψ )sin( β)cos(Lon − λ) + cos(I)cos( ψ )cos( β)sin(Lon − λ) ] − r sin(V − E) M ,O sat sat
 Z = R[]cos(I)sin( ψ ) − cos( ψ )sin(I)sin(Lon − λ)M ,O sat
Descripción del escenario de prueba y resulta-
dos Finalmente, la columna y la fila ocupadas en la
imagen por el píxel M serán:
Para llevar a cabo un análisis de verificación del
modelo expuesto, se han analizado los resultados  Z−1 M ,0 tan obtenidos tras la aplicación del modelo sobre imá- X t − t M ,0 converg prim _ linea    C = +1024.5 ; F = *6 genes del sensor NOAA-AVHRR. Con la finalidad M M  δ 1000i   de abarcar un análisis más global, se ha utilizado
una configuración completa de todas las situacio-
siendo t el tiempo para el que converge la converg nes geométricas posibles, que en el caso del sensor
ecuación resuelta por el método de bisección, AVHRR se compone de 10 imágenes consecutivas
t el tiempo de la captura de la primera línea prim_linea correspondientes a la imagen diaria de las 13:00
de imagen y δ el paso del espejo (para el AVHRR i GMT (aproximadamente), en trayectorias ascen-
con un valor cercano a 0.95 miliradianes). Si la dentes sur-norte, que cubren la totalidad de la
trayectoria es ascendente la fila será la comple- Península Ibérica. Esta configuración está mostra-
mentaria al tamaño en filas de la imagen. da en la figura 2, en la que puede verse desde la
primera imagen de la serie, con una longitud de
Proceso de remuestreo de la imagen cruce con el ecuador terrestre, en sentido ascen-
dente, de 353.1945º, hasta la última imagen de la
Para crear las imágenes corregidas geográficamen- serie con longitud de cruce de 18.7047º. Ambas
te se ha aplicado un proceso de remuestreo en situaciones son posiciones extremas que provocan
Nº 14 – Diciembre 2000 5 de 9 A. Calle, J. L. Casanova y J. Sanz
1 2
3 4
5 6
7 8
9 10
Figura 2: Situaciones geométricas analizadas para el sensor AVHRR.
una gran deformación geométrica con la trayecto- dientes a la secuencia descrita, seleccionando zo-
ria sub-satélite (línea vertical en el centro de la nas distantes del nadir. De esta manera se logra
imagen). Una descripción más detallada de la una descripción completa del procedimiento. En
configuración de los parámetros de cruce con el cada escena se ha especificado la siguiente infor-
ecuador puede obtenerse a través del análisis de la mación: Día de la secuencia (véase la figura 2),
figura 3 donde se especifican los parámetros distancia al nadir promedio en unidades de pixeles
TEQX (tiempo de cruce con el Ecuador) y LEQX (ND) que da una idea clara del alejamiento de la
de todas las imágenes en estudio. En lo que respec- escena y por consiguiente de los errores de pers-
ta al sensor AVHRR, las características que afec- pectiva que pueden afectar. Esta distancia está
tan a la geometría son los datos que siguen: tomada con el origen en la traza del satélite y va-
riando en el intervalo [0,-1024] para las zonas
FOV = ±55.4º; Barrido del sensor= 2500 km; situadas a la derecha de la imagen y [0,+1024]
Velocidad de captura de lineas: 166 millisegundos por linea
para las zonas situadas a la izda. Obviamente y Ángulo de paso del espejo del sensor: 0.95 mrad ( δ ) i
Número de pixeles por línea: 2048 como se desprende del modelo, los errores aumen-
tan a medida que aumenta la distancia a la traza del
Parámetros de cruce con el Ecuador satélite. X-err y Y-err son los errores promedios
encontrados en la escena. 56000 350
L55000T Las conclusiones más destacables a partir de la 300 EE 54000 figura 4 son las siguientes: El error encontrado 250 QQ 53000
XX 200 para los valores de fila y columna están compren-
52000
150 didos en la mayoría de los casos en el intervalo ±2 51000
10050000 pixeles con la excepción de situaciones dotadas de
5049000 gran deformación; téngase en cuenta que la escena
48000 0 correspondiente al día 8, con un X-err=+3 equivale
151 16 2 3 17 184 195 620 7 21 8 22 9 23 1024
a estar situado en un área extrema del 23% de la Día de la secuencia
TEQX LEQX imagen, y en el caso aún más extremo de deforma-
Figura 3: Parámetros de cruce con el Ecuador en la secuencia ción (véase la figura 2) de la escena correspondien-
analizada
te al día 10, con un X-err=+4, la distancia al nadir
es de 880, lo que representa a un área extrema de
la imagen del 14%. Otra importante conclusión es Los resultados gráficos obtenidos tras el análisis
individual de cada una de las imágenes puede que el error encontrado en las filas es independien-
apreciarse en la figura 4, donde, por brevedad, se te de la distancia al nadir, producido muy proba-
muestran las escenas más representativas de los blemente por errores en la determinación orbital
(un error de 167 milisegundos equivale a una línea resultados de la serie completa. Se ha representado
una escena de una imagen de cada 2, correspon- de imagen); sin embargo, el error en las columnas
6 de 9 Nº 14 – Diciembre 2000 Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores
Día 4
ND=152
X-err = 0
Y-err = +1

Día 2 Día 8
ND=670 ND=-784
X-err = -1 X-err = +3
Y-err = +2 Y-err = 0
Día 6
ND=544 Día 10 ; ND=-880
X-err = -1 X-err = +4
Y-err = +2 Y-err = 0
Figura 4: Resultados obtenidos
si presenta dependencia con la distancia al nadir
(como era esperable)

Análisis de errores
66%
El análisis de los errores encontrados debe reali-
zarse en función y localización dentro de la ima-
gen. Por ello, a partir de los resultados obtenidos
en el anterior apartado podemos realizar un es-
16% 16%quema de zonas de riesgo (o zonas de validez, si se
quiere). De forma aproximada y para ofrecer una
12% 12%visión cualitativa, podemos concluir en que en el
7% 7% 66% central de la imagen, los errores en X son del
orden de ±1; en el 8% siguiente hacia los extre- Figura 5: Niveles de X-err en una imagen completa de 2048
pixeles.X-err= ±1y ±2, ±3, ±4 mos, los errores son del orden de ±2 (estas áreas
están solapadas y no existe un límite claro); el
las calculadas a través del modelo, debido a laárea del 10% siguiendo hacia los extremos, los
limitación de la resolución espacial. Además, en errores aproximados son ±3 y finalmente en el área
las zonas de máxima distorsión, que coinciden con
de los extremos equivalente al 14%, el error es ±4.
la resolución espacial más baja y las máximas Este es el mejor resultado para esta zona (léanse
distancias al nadir, se repiten los mismos valores
las otras) ya que si dividiéramos esta zona extrema
de un píxel para ocupar distintas posiciones (ob-encontraríamos mayores errores en los puntos más
sérvense distintas situaciones de las imágenes
extremos; por esa razón esta clasificación debe
analizadas en el apartado anterior). Estos proble-entenderse de forma orientativa. Véase la figura 5.
mas podrían eliminarse si se utiliza el modelo de
Respecto a la fuente de errores, en primer lugar
corrección de forma inversa; es decir, deformando no debe olvidarse que el proceso de verificación se
el mapa a superponer en vez de corregir la imagen.
está realizando sobre imágenes remuestreadas,
Dado que el mapa es un vector, carece de la limi-Debe tenerse en cuenta, como luego se menciona-
tación de la resolución espacial
rá, que los errores de un proceso de corrección
remuestreando la imagen son mayores que los
En lo que respecta al error en X, dado que es
producidos al deformar un mapa. El remuestreo,
creciente con la distancia al nadir, puede estar sin embargo es indispensable para el posterior
provocado por dos causas: en primer lugar el mo-
tratamiento de la imagen. Este procedimiento pue-
delo de elipsoide elegido, o incluso por la exacti-de dar lugar a georreferenciar pixeles cuyas coor-
tud con la que la superficie de la tierra puede estar
denadas geográficas no coincidan exactamente con
perfectamente modelizada por esa superficie. La
Nº 14 – Diciembre 2000 7 de 9 A. Calle, J. L. Casanova y J. Sanz
otra causa de propagación de errores apuntada es la La figura 7 representa el error de tiempo del re-
referida a imprecisiones cometidas al cuantificar el loj del satélite NOAA-14 durante el período
ángulo de paso del espejo del sensor, al que nos 01/01/01/01/94 94 (a (a pesar de no estpesar de no estaar en vueloo, se represen-, se represen-
ta para establecer el cero) hasta 29/08/00. Como se hemos referido en anteriores apartados como δ . Es i
aprecia, salvo excepciones, siempre sufre retraso, de esperar que al cometer un error ∆δ , el error en i
que ha llegado a alcanzar 1.5 segundos de magni-la localización de los pixeles se vea incrementado
tud. a medida que aumenta la distancia al nadir del
satélite. En este sentido, hemos simulado un error
-3en el paso del espejo del sensor ( ∆δ) de 6·10 CONCLUSIONES i
miliradianes, de forma arbitraria. Este resultado
La decisión de adoptar un modelo automático de puede consultarse en la figura 6 que representa la
corrección debe estar analizada sopesando la nece-dependencia del error en la geo-referenciación, con
sidad de obtener resultados automáticos sacrifi-la distancia al nadir. La ecuación del ajuste que
cando la precisión. Aún en ese caso, como hemos aporta la expresión del error en función de la dis-
visto, si las imágenes utilizadas son de buena cali-tancia al nadir es:
dad en cuanto a la resolución espacial, el proceso
autautoommááttiicco puede ser adopto puede ser adoptado con garantado con garantííaas. s.
−8 3 −6 2 X −error =2·10 *ND +1·10 ·ND +0.0145·ND +7.15 Hemos establecido una caracterización del error
en la dimensión X; en el caso de la Y, mejores
resultados serán encontrados conociendo con exac-siendo ND la distancia del píxel al nadir. Este
polinomio de grado 3 es el que aporta el mejor titud el retraso del reloj del satélite, y mejorando el
2 modelo orbital o utilizando otros más exactos coeficiente de correlación (R =1). La magnitud del
polinomio cúbico muestra la incidencia de la im- como el SGP4.
Por otPor otra partra parte, heme, hemoos apls apliicado elcado el m modelodeloo a a iimmá-á-precisión en el espejo del sensor sobre el proceso
de geo-referenciación. genes de muuyy di diferentferentes es zonas zonas geográfigeográficas cas (Perú, (Perú,
con órbitas ascendentes) encontrando resultados
Una de las principales dificultades encontradas similares a los aquí expuestos, por lo que se puede
concluir en la generalidad del modelo. en la aplicación del proceso de validación ha sido
el posicionamiento exacto del satélite a partir de
los parámetros orbitales, a partir del tiempo de APÉNDICE: ECUACIONES DE
capcaptutura dra dee la la pprrimimera era línlínea de imagen, en unidades agen, en unidades TRANSFORMACIÓN DE
de mde milisegundos. En el proceso de validación con con COORDENADAS. el sensor AVHRR, se obtenía este dato a través de
la propia transmisión del satélite. Se comprobó la
existencia de un error sistemático de desplaza- Cálculo de las coordenadas UTM-30N
miento en las imágenes georeferenciadas y se
pensó en un posible error en el reloj del satélite (en Las coordenadas geográficas, longitud y latitud,
este caso el NOAA-14). Esta imprecisión fue con- se obtienen a partir de las correspondientes coor-
firmfirmada a través del Departamento ento de de Navegación denadas UTM referi referidas aldas al huso 30 nort huso 30 nortee, X, X y y UTMUTM
de NESDIS/NOAA (E. Harrod, 2000). En el caso Y mediante las expresiones siguientes: UTM,
concreto del NOAA-14, el error más elevado que
hemos encontrado ha sido del 1500 milisegundos  3 2Lat = 1 + E cos ()ζ − E sin()ζ cos()ζ ( τ − ζ) ( τ − ζ) + ζ 2 2el 30 de junio de 1997 (fuente: NOAA Navigation 2 
Department), lo que equivale a un error aproxima- Lon = α + L0
do a nivel de la superficie terrestre, de 10 kilóme-
ttrros. Sios. Si se di se dispone de estee dat datoo, no es di, no es difífícicill iinnttrro-o- habiendo calculado previamente las funciones
ducirlo ducirlo en en el mel modelo de corrección. En el caso de matemáticas que se refieren a continuación:
las escenas analizadas en el presente trabajo, la
corrección ha sido de 1100 milisegundos. Y C  X − 5e5 1UTM 1 UTM 2 ζ = ; a = ; a = E cos ()ζ 1 2 22  CE 2 a1 + E cos ()ζ2  1 2
X − 5e5  a  sinh( ψ )UTM 2ψ = 1 − tan( α) = 500
a  3  cos( η)1
Y − AUXLAT()ζ0 UTMη = ()1 − a + ζ tan( τ ) =cos( α) tan( η)2
a1
-500
-1000 donde las constantes tienen los valores:

-1500
L =(-3* π/180.) 0
-2000 CE2=(6366197.724/0.9996)
Días transcurridos (01/01/94-29/08/00) E2=(0.0067681703)
Figura 7: Error de tiempo (milisegundos) del reloj del NOAA-14
8 de 9 Nº 14 – Diciembre 2000
Error del reloj
(milisegundos)Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores
de Meteorologie Spatiale (CMS) B.P. 147-Lannion Función auxiliar de cálculo AUXLAT
Cedex. SATMOS notes No 2.
BRUNEL, P. and MARSOUIN, 2000. Operational La función auxiliar de cálculo AUXLAT, para el
AVHRR navigation results. Int. Journal of Remote parámetro, ϕ, se calcula de acuerdo a la siguiente
Sensing, 2000, Vol. 21 No.5, pp 951-972.
ecuación: CRAWFORD, P.S. 1995. Kepler’s equations in
‘C’. Int. Journal of Remote Sensing, 1995, Vol.
AUXLAT( ϕ) =C *( ϕ − b α +b β −b γ ) 16 No.3, pp 549-557. 1 1 2 3
HARROD, EMILY D. NOAA/NESDIS/OSDPD/IPD.
E/SP13 Room 0308 FB#4 5200 Auth Road habiendo calculado las funciones matemáticas:
Suitland, Md. 20746-4304. Comunicación personal.
2 4 HO, D. and ASSEM, A. 1986. NOAA AVHRR sin(2 ϕ) 3b +sin(2 ϕ)cos ( ϕ) 5b +sin(2 ϕ)cos ( ϕ)1 2b = ϕ + ; b = ; b =1 2 3 image referencing. Int. Journal of Remote Sen-2 4 3
3 5 352 3 sing, 1986, Vol. 7, pp 895-904. α = E ; β = α ; γ = α2
4 3 27 ILLERA, P., DELGADO, J.A. and CALLE, A.
C1 = 6397376.634, y E2 = 0.0067681703
1996. A navigation algorithm for satellite ima-
ges. Int. Journal of Remote Sensing, 1996, Vol.
REFERENCIAS 17 No.3, pp 577-588.
MORENO, J AND MELIÁ, J. 1993. A Method for BACHMANN, M. and BENDIX, J. 1992. An
Accurate Geometric Correction of NOAA AVHRR
improved algorithm for NOAA-AVHRR image HRPT Data. IEEE Transactions on Geoscience and
referencing. Int. Journal of Remote Sensing, Rem. Sen., Vol 31, No.1, 204-226
1992, Vol. 13 No.16, pp 3205-3215. ROSBOROUGH, G.,W., BALDWIN, D.G. AND
BRUNEL, P. and MARSOUIN, 1986. Geographical EMERY, W.J. 1994. Precise AVHRR Image Naviga-
navigation of NOAA AVHRR series imagery. Centre tion. IEEE Transactions on Geoscience and Rem.
Sen., Vol 32, No.3, 644-657

Nº 14 – Diciembre 2000 9 de 9

Un pour Un
Permettre à tous d'accéder à la lecture
Pour chaque accès à la bibliothèque, YouScribe donne un accès à une personne dans le besoin