Tema 7. Corrección atmosférica de imágenes de satélite

1. Fundamentos teóricos

Una imagen de satélite en bruto contiene unos valores numéricos denominados niveles digitales (ND) que el satélite obtiene a partir de la energía recibida mediante una ecuación lineal. Para recuperar los valores de energía recibida es necesario aplicar la inversa de esa ecuación lineal.

La energía recibida por el sensor es por su parte la suma de la energía reflejada por el terreno y la energía dispersada y reflejada por la atmósfera. Por otra parte la energía que llega al suelo es la energía solar por la transmisividad de la atmósfera.

Esuelo = t·E0

Esatelite = r·Esuelo + Eatmosfera

r = Esatelite - Eatmosfera / t·E0

Este esquema se complica aún más ya que sólo es válido en terreno perpendicular a los rayos solares. En caso contrario:

Esuelo = t·sin(a)·E0

Donde a es el ángulo que forman los rayos solares respecto a la superficie del terreno. Ya que la cantidad de radiación que llega a un punto del terreno depende de este ángulo.

Es necesario corregir todos estos problemas ya que de no ser así:

Aparecerían efectos debidos a la mayor o menor iluminación
Los resultados aparecerían difuminados por el efecto de la atmósfera
Sería difícil comparar dos imágenes de dos fechas distintas con diferente influencia de la atmósfera

2 Obtención de parámetros

Existen métodos sofisticados y complejos para obtener los parámetros relacionados con la corrección atmosférica sin embargo requieren gran abundancia de datos acerca del estado de la atmósfera en un momento determinado. Sin embargo existen otros métodos más sencillos y que funcionan relativamente bien:

2.1 Mínimo del histograma

Sirve para eliminar la energía reflejada por la atmósfera hacia el satélite. Se basa en la asunción de que, puesto que la reflectividad en el oceano es prácticamente nula, toda la energía recibida desde los pixels de oceano procede de la atmósfera. Estos pixels de oceano coincidirán con los valores mínimos de los histogramas de las imágenes de ahí el nombre del método.

En aquellos casos en que no se dispone en la imagen de pixels de oceano se podrían utilizar embalses aunque los resultados son evidentemente peores ya que un embalse pude ganar reflectividad debido a la turbidez, baja profundidad, desarrollo de algas, etc.

2.2 Cálculo de iluminación

Sabiendo la pendiente y orientación de una porción del terreno y la posición del sol puede calcularse fácilmente el ángulo que forma el sol con la superficie del terreno.

La posición del sol se obtiene a partir de la latitud del lugar el día y la hora y la pendiente y orientación del terreno a partir de un modelo digital de terreno.

3 Aplicación en GRASS

En primer lugar vamos a calcular sin(a). Para ello es necesario disponer de mapas de pendiente y orientación que pueden obtenerse apartir del mapa de elevaciones con el módulo r.slope.aspect

r.slope.aspect elevation=name slope=name format=degrees aspect=name min_slp=5

Una vez obtenidas pendientes y orientaciones se utiliza el módulo r.sun para obtener un mapa del ángulo que formaban los rayos solares con la superficie del terreno en el día y la hora del paso del satélite

r.sun elevin=mde50M aspin=orientacion50M slopein=pendiente50M incidout=angulo latitude=38 dej=60 lum_time=10.7

Para obtener los senos de los ángulos se utiliza r.mapcalc:

r.mapcalc 'sinang=eval(pordef=sin(41.17),if(isnull(angulo),pordef,sin(angulo)))'

Los valores de radiación recibida en el sensor, descontada la radiación procedente de la atmósfera, pueden obtenerse con una simple operación de r.mapcalc. Se trata de un programa que actúa como una calculadora de mapas de propósito general. La forma general sería

r.mapcalc 'Esensorb=(+NDb-minimob+(-a0b/a1b))*a0b+a1b';

donde NDb es el mapa de niveles digitales para la banda b, minimob, a0b y a1b los valores de mínimo del histograma, a0 y a1 respectivamente para la banda b. Estos valores aparecen en la tabla 1

Tabla 1: Coeficientes necesarios para obtener valores de energía a partir de los niveles digitales mediante la ecuación E=a0+a1ND

Banda	a0	a1	E0	Mínimo
1	-6.2	0.786	1970	45
2	-6	0.817	1843	25
3	-4.5	0.64	1555	17
4	-4.5	0.635	1047	14
5	-1	0.128	227.1	7
7	-0.35	0.0424	80	78

CALCULO DE REFLECTIVIDADES

r.mapcalc 'albedob=eval(t=Esensorb*3.141592/(sinang*E0b),if(t>1,null(),t))'

Con esta última operación se obtiene una buena aproximación al albedo.