Source: http://www.youscribe.com/catalogue/documents/autres/algoritmo-automatico-de-correccion-geografica-para-imagenes-2012763
Timestamp: 2017-07-26 10:34:10+00:00

Document:
Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores - A. Calle - Autres
Algoritmo automático de corrección geográfica para imágenes NOAA-AVHRR: desarrollo y análisis de errores
A. Calle (Auteur)
ResumenSe 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.AbstractA 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.
Revista de Teledetección 1988-8740 2000 número 14
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-0desarrollo 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 0aplicació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 0fá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.52 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) −1Nde 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 orbitaFinalmente, 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 =0orbitNº 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 , MF (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 iello, 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 Mniente 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 satcon 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-Mpresió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)sattan( β) − − = 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 10Figura 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 ( δ ) iNú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 350L55000T Las conclusiones más destacables a partir de la 300 EE 54000 figura 4 son las siguientes: El error encontrado 250 QQ 53000XX 200 para los valores de fila y columna están compren-52000150 didos en la mayoría de los casos en el intervalo ±2 5100010050000 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 1024a estar situado en un área extrema del 23% de la Día de la secuenciaTEQX 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 = 0Dí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 2048pixeles.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 iaprecia, salvo excepciones, siempre sufre retraso, de esperar que al cometer un error ∆δ , el error en ique 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 imiliradianes, 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 UTMUTMde 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 = α + L0do 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 2X − 5e5  a  sinh( ψ )UTM 2ψ = 1 − tan( α) = 500a  3  cos( η)1Y − AUXLAT()ζ0 UTMη = ()1 − a + ζ tan( τ ) =cos( α) tan( η)2a1-500-1000 donde las constantes tienen los valores: -1500L =(-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 3HARROD, 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 ; β = α ; γ = α24 3 27 ILLERA, P., DELGADO, J.A. and CALLE, A. C1 = 6397376.634, y E2 = 0.00676817031996. 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

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución