Notas de investigación

 

Desarrollo de curvas espectrales del crecimiento anual de la vegetación, usando sensores remotos

 

Development of remote sensing vegetation spectral growth curves

 

Fernando Paz-Pellat^a, Martín Alejandro Bolaños-González^a, Luis Alberto Palacios-Sánchez^a, Fermín Pascual-Ramírez^a, Edgardo Medrano-Ruedaflores^a, Fernando Ibarra-Hernández^a

 

^a Colegio de Postgraduados, Campus Montecillo. pellat@colpos.mx. Correspondencia al primer autor.

 

Recibido el 23 de marzo de 2010.
Aceptado el 18 de noviembre de 2010.

 

RESUMEN

El desarrollo de curvas del crecimiento de la vegetación (pastizales y matorrales), a nivel histórico, plantea retos importantes por falta de estadísticas en México. La tecnología de sensores remotos en plataformas espaciales, por medio de índices de vegetación o IV, permite realizar una aproximación en esta dirección. Este trabajo discute un esquema de modelación paramétrica (indirecta) de la curva de crecimiento de la biomasa en píxeles individuales de sensores Advanced Very High Resolution Radiometer o AVHRR y Moderate Resolution Imaging Spectroradiometer o MODIS, permitiendo así su caracterización usando algoritmos de pre y post-proceso. Los resultados obtenidos de la aplicación de los algoritmos discutidos en este trabajo, permiten la obtención de curvas anuales del crecimiento a la escala de 100 ha (AVHRR) y 6.25 ha (MODIS), que pueden ser usadas en estudios de cambio climático y del manejo de agostaderos en México.

Palabras clave: Curvas de crecimiento, Agostadero, índices NDVIcp e IVCP, Modelo paramétrico, AVHRR, MRIS, Cambio climático.

 

ABSTRACT

Development of historic vegetation growth curves for both grasslands and shrublands in Mexico becomes a real challenge due to lack of statistics. Remote sensing technology using satellites, allows obtaining a first approximation through vegetation indices (VI). The present study discusses an indirect parametric modeling of the biomass growth curve in individual pixels of both Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS), thus allowing their characterization using pre and post process algorithms. Results obtained through application of the algorithms discussed in the present paper, allow obtaining annual growth curves at 100 ha (AHVRR) and 6.25 ha (MODIS) scales that can be used in studies on climate change and on grassland/shrubland management in Mexico.

Key words: Growth curve, Gassland/shrubland, NDVIcp index, IVCP index, Parametric model, AVHRR, MODIS, Climate change.

 

El contar con curvas de crecimiento de la vegetación (pastizales naturales y matorrales, principalmente) a escala de país, permite analizar los efectos del cambio climático, así como el estado de este recurso en relación a las actividades ganaderas, entre otras. En especial, una base histórica de curvas de crecimiento a escala de celdas de una malla nacional (celdas de dimensiones de 1,000, 500 y 250 m) permitiría analizar los procesos de herbivoría que impactan la condición del agostadero en México, además de generar información para la estimación de las capacidades de carga animal locales en forma dinámica.

Los sensores remotos en plataformas satelitales, disponibles desde los 70s, son una tecnología de carácter espacialmente exhaustivo, y con resolución temporal diaria de revisita al mismo sitio (escalas de 1,000 a 250 m), donde sobresale el sensor AVHRR (Advanced Very High Resolution Radiometer; resolución espacial de 1,100 m a nadir; con información de cobertura de México de 1986 a la fecha) y MODIS (Moderate Resolution Imaging Spectroradiometer; resolución espacial de 1,000 a 250 m, con información de 2000 a la fecha). Así, los sensores remotos brindan una oportunidad única de reconstruir la historia de manejo y efectos climáticos de la vegetación en México.

Los sensores remotos caracterizan la vegetación usando el alto contraste entre una banda fotosintética como el rojo (R) y una no fotosintética como el infrarrojo cercano (IRC)⁽¹⁾. Este par de bandas están disponibles en los sensores AVHRR y MODIS, por lo que pueden usarse para el desarrollo de índices espectrales de la vegetación o IV, de los cuales se ha publicado gran cantidad^(2,3,4). En el crecimiento de la vegetación, en las fases iniciales, los sensores remotos observan una mezcla suelo-vegetación (vegetación con cobertura parcial del suelo; caso típico en pastizales y matorrales en la ganadería extensiva del país), por lo que los IV generalmente son diseñados para minimizar el efecto del suelo y obtener información únicamente de la vegetación⁽⁵⁾.

El objetivo de este trabajo fue el desarrollo de una metodología para parametrizar la curva anual de crecimiento espectral de pastizales y matorrales a nivel nacional, y poder realizar análisis de cambios temporales a escala de pixeles individuales en imágenes satelitales. En este trabajo se utilizó el índice NDVIcp⁽⁴⁾ o Índice de Vegetación de Diferencias Normalizadas del Colegio de Postgraduados, para estimar el crecimiento de la vegetación usando un modelo paramétrico. Este modelo, usando una transformación del NDVIcp, es la base para el desarrollo de un algoritmo de pre-proceso y post-proceso de las bandas espectrales del sensor AVHRR (extensible al MODIS) para generar curvas anuales de crecimiento de la vegetación a escala píxel por píxel de una imagen satelital.

El índice NDVIcp fue diseñado para tener una relación lineal con el índice de área foliar o IAF en la etapa vegetativa de la vegetación (hasta un IAF aproximado de 3), minimizando el efecto del suelo debajo de ésta. Considerando que el IAF tiene una relación aproximadamente lineal con la biomasa aérea o Bm⁽⁶⁾ podemos hacer estimaciones robustas de Bm; aun para un umbral del nivel de estrés, donde la vegetación cambia su estrategia de acumulación de IAF y Bm, para incrementar las reservas de fotosintatos en la parte radical⁽⁷⁾.

La Figura 1a muestra el patrón temporal generalizado del crecimiento anual de la vegetación (pastizales) para IAF. En caso de matorrales, el IAF no inicia en cero.

La Figura 1b muestra el patrón descrito en la Figura 1a, pero para el espacio espectral del R-IRC, suelo debajo de la vegetación variando desde suelo seco (S1) a húmedo (S2), representando las condiciones esperadas durante el monitoreo por sensores remotos de un píxel.

Los patrones mostrados en la Figura 1, curvas iso-suelo (la vegetación varía y las propiedades ópticas del suelo son constantes), tienen los siguientes patrones:

Patrón V. Esta etapa va desde la emergencia (línea del suelo o IAF = 0; en matorrales es una línea para un IAF dado específico) hasta la antesis. En el espacio temporal, se caracteriza por un crecimiento exponencial del follaje (fase V1), hasta lograr el cubrimiento máximo del suelo (fotosintéticamente) y después hay una transición a un crecimiento tipo lineal (fase V2). En alguna vegetación se presenta otra fase exponencial (fase V3) después de la fase lineal, y se caracteriza por cambios en las propiedades ópticas de la hoja y nuevos órganos con propiedades ópticas diferentes a la hoja (fase reproductiva). La fase exponencial inicial, en el espacio espectral, se caracteriza por la evolución de las curvas iso-IAF (igual vegetación con propiedades ópticas diferentes del suelo) hasta alcanzar el punto de saturación del rojo (reflectancia de medio denso o R¥). La fase lineal está representada por una línea recta vertical, donde el valor del R permanece constante, hasta alcanzar el valor máximo de reflectancia del infinito de la banda del IRC (IRC¥). La fase exponencial 2 presenta un patrón donde el R se incrementa (cambios de color) y el IRC permanece constante o decrece.

Patrón M. Va desde la antesis hasta antes del inicio de la senescencia y se caracteriza por un IAF constante (nacimiento y muerte equilibrado de hojas). En el espacio del R, este patrón se caracteriza por un punto donde R y el IRC permanecen constantes o por una reducción del IRC.

Patrón S. Va desde el inicio de la senescencia hasta donde la Bm está totalmente muerta (línea del residuo o IAF = 0; en matorrales es una línea para un IAF dado específico). El espacio temporal se caracteriza por un patrón exponencial. El patrón es similar en el espacio del R-IRC y se caracteriza por el incremento del R debido a la degradación de los pigmentos de las hojas y la reducción del IRC por el deterioro de la estructura del mesófilo.

Curva de crecimiento de la vegetación usando el índice NDVIcp

Para analizar el patrón de crecimiento asociado al NDVIcp, se utilizan los datos experimentales de un pastizal natural (Bouteloua spp.) en la cuenca de Walnut Gulch en Arizona, Estados Unidos. La documentación de este experimento puede encontrarse en estudios publicados^(8,9,10). Bolaños et al⁽¹¹⁾ analizan este experimento para el desarrollo de un modelo de la geometría sol-sensor en las reflectancias. Para contar con mediciones simultaneas del NDVIcp y la Bm, se interpolaron algunos datos (polinomios) para hacer coincidir las mediciones disponibles.

Puesto que las reflectancias de la vegetación son dependientes de la cantidad de hoja y sus propiedades ópticas⁽¹²⁾, la comparación de índices espectrales debe ser realizada en términos de biomasa verde o viva. La relación entre la Bm y el NDVIcp es de tipo no lineal para la etapa vegetativa y de senescencia. Para aproximar linealmente una relación del NDVIcp con la Bm, se usó la transformación 1/NDVIcp, que muestra buen comportamiento. La ventaja de la transformación 1/NDVIcp, es que permite establecer modelos lineales de crecimiento de la vegetación en el espacio del tiempo (Figura 2). Esto es importante ya que simplifica los algoritmos operacionales para la estimación (parametrización) de los patrones anuales del crecimiento de los pastizales y matorrales a una escala píxel por píxel. El modelo paramétrico (funciones lineales) puede ser ajustado (definición de sus parámetros) usando mediciones del índice NDVIcp. De esta manera, el modelo de crecimiento de la vegetación es usado y no los valores directos de las mediciones del NDVIcp en imágenes satelitales.

Modelo temporal general para el índice NDVIcp

De acuerdo a los patrones mostrados en la Figura 2, la Figura 3a muestra el modelo desarrollado para la evolución temporal (X= día juliano) de 1/NDVIcp (=Y). En la Figura 3b se muestra el patrón temporal del NDVIcp, el cual no es lineal.

En la Figura 3a, al inicio de la etapa vegetativa, se presenta una meseta inicial (suelo desnudo o vegetación perenne). La etapa vegetativa inicia en (X1, Y1) y termina en (X2i, Y2). Después de esta etapa se considera una meseta intermedia entre (X2i, Y2) y (X2f, Y2) y al final, una etapa de senescencia entre (X2f, Y2) y (X3, Y3). El modelo propuesto no considera en forma completa la etapa inicial de la fase exponencial inicial (V1), ni la exponencial final (V3), de la etapa vegetativa-reproductiva, por la aproximación lineal utilizada. Esto no tiene mayores implicaciones en la práctica, ya que estas transiciones no lineales son aproximadas por la intersección de dos líneas rectas (ver más adelante).

Las dos etapas generales del crecimiento de la vegetación, utilizando 1/NDVIcp, pueden estimarse de:

Equivalencia entre el NDVIcp y la biomasa aérea de la vegetación

Considerando un espacio geográfico dado (píxel de una imagen satelital) y un tiempo específico (por ejemplo, el tiempo donde Bm la alcanza su valor máximo en el año), la capacidad de carga animal está relacionada con la Bm seca total (kg M.S.) por:

donde FA es el factor de aprovechamiento de la biomasa (0.5 u otro valor), FAG es el factor de biomasa aprovechable por el ganado (dependiente de la composición florística y disponibilidad forrajera en cada sitio; de 0 a 1), FAA es un factor de ajuste por condiciones abióticas (cercanía a cuerpos de agua y pendiente del terreno⁽¹³⁾); y 4,928 kg MS es el requerimiento de alimento para mantener una unidad animal todo el año (ingesta diaria del 3 % del peso corporal de un bovino con un peso de 450 kg). La Bm está dada en términos absolutos y es el resultado de multiplicar la densidad de biomasa (kg m^-2) por la superficie (m²) que representa esta estimación, para un tiempo específico o ciclo. En el caso de usar la biomasa viva, se puede usar otro factor de ajuste por contenido de humedad.

Si consideramos un par de años consecutivos (k y k-1), suponiendo que los factores de ajuste no cambian apreciablemente, podemos establecer entre ellos la relación:

La ventaja de la relación de razones dada por la ecuación (3) es que filtra los factores de ajuste, evitando así tener que calibrarlos localmente. Este tipo de razonamiento puede extenderse para una serie de tiempo de la Bm, desde el año mal año n:

donde las r son razones entre años consecutivos, que reflejan los cambios de la biomasa. Así, es posible construir indicadores del estado de la vegetación en un píxel en relación a un año base. Las razones r pueden analizarse en relación a condiciones de manejo (nivel de pastoreo) o climáticas (precipitación).

De la relación (3) y (4), para evitar el uso de factores de ajuste locales, es necesario que la relación entre la Bm y un índice de vegetación IV cumpla:

Dada la relación entre 1/NDVIcp y la Bm, relación lineal con una constante aditiva (Figura 2), es necesario hacer una transformación para obtener una relación similar a la mostrada en (6):

El nuevo índice IVCP (índice de Vegetación del Colegio de Postgraduados), permite obtener una relación similar a (3), ya que la constante aditiva entre 1/NDVIcp y Bm es filtrada. El valor de (1/ NDVIcp)_min se puede obtener de la evolución temporal del índice 1/NDVIcp. Así, considerando a las etapas vegetativa y de senescencia en forma independiente, en la etapa vegetativa se puede usar Y1 = (1/NDVIcp)_min y en la de senescencia Y3 = (1/NDVIcp)_min. La Figura 3c muestra el modelo del IVCP para esta estrategia de parametrizar (1/NDVIcp_min.

Para poder contar con un registro histórico de curvas de crecimiento para valorar la condición del agostadero, se utilizaron imágenes satelitales del sensor AVHRR^(14,15), disponibles para México desde 1986 (cobertura nacional). El software utilizado denominado SPIA (Sistema para Proceso de Imágenes AVHRR) consiste de dos algoritmos generales: pre-proceso y post-proceso.

Algoritmo de pre-proceso

El procesamiento de imágenes satelitales del sensor AVHRR, algoritmo de pre-proceso, consistió de los siguientes pasos:

Referenciación geográfica. Las imágenes satelitales del sensor AVHRR, después de su lectura, son referenciadas geográficamente usando un modelo orbital⁽¹⁶⁾ para su re-proyección en coordenadas Cónicas Conforme de Lambert en una malla nacional con celdas de 1 km x 1 km, usando el método del Vecino más Cercano. El error de referenciación fue de 1 a 3 píxeles, que no tiene repercusión importante, para la implementación, usando polígonos de uso del suelo 1:250,000 del INEGI.

Discriminación de nubes y sombras. El software SPIA realiza un proceso de generación de códigos de calidad de la información contenida en cada píxel. E n especi al , este p roceso sirve p ara discriminar píxeles con nubes, los cuales son sacados del proceso del NDVIcp. La base del algoritmo de discriminación de nubes el algoritmo CLAVR-1(Clouds from AVHRR-Phase I) usado por la NASA⁽¹⁷⁾ (North American Space Agency).

Al aplicar el algoritmo CLAVR-1, modificado, las posibles clases de objetos son: claro, sombra, nieve/ hielo. Siguiendo el esquema propuesto en el algoritmo CLAVR-1, se utiliza la clase "mixto" (p= parcial) para denotar la situación donde existe ambigüedad (una o más pruebas fallan para etiquetarlo como nube o claro) en la clasificación:

En la porción solar del espectro electromagnético (0.3 a 3 μm), la radiación solar incidente domina sobre la radiación térmica emitida. La cantidad de radiación reflejada depende de la composición química y física de la superficie terrestre y la atmósfera. Esta composición es lo que distingue a las clases de superficies terrestres y a constituyentes atmosféricos. Las diferencias en su composición, generalmente, hacen a algunos altamente reflectivos (nubes, nieve y hielo), a otros moderadamente reflectivos (desiertos, suelo desnudo) y, todavía, a otros, los hacen poco reflectivos (océano, vegetación densa). Esta propiedad puede ser usada para detectar la presencia de estos componentes, utilizando las diferencias en sus reflectancias.

Las diferencias en composición también pueden causar que estas componentes reflejen la radiación solar en forma más eficiente en ciertas partes del espectro, que en otras. Esta propiedad puede ser usada para discriminar la presencia de ciertos componentes superficiales/atmosféricos. Por ejemplo, la nieve refleja altamente en longitudes de onda menores a 1 μm y pobremente, más allá de 1.5 μm, mientras que la vegetación verde refleja pobremente debajo de 0.7 μm, pero lo hace en forma moderada más allá de 0.8 μm. Las nubes tienden a dispersar la radiación en forma uniforme con las longitudes de onda, mientras que la bruma tiende a dispersar más en longitudes de onda menores.

En el infrarrojo térmico (3 a 20 μm), los componentes superficiales/atmosféricos pueden ser discriminados, pero como resultado de procesos físicos diferentes: absorción y re-emisión de radiación, en vez de reflexión y dispersión. La radiación infrarroja es proporcional a la temperatura de la escena, por lo que es posible separar superficies calientes y frías (por ejemplo, nubes sobre desiertos en el día). Las firmas espectrales son aquellas cuya emisividad varía con la longitud de onda, tal como entre las nubes cirrus delgadas y el vapor de agua atmosférico en la región de 10 a 12 μm del espectro electromagnético.

Las técnicas para discriminar pfxeles claros usando el sensor AVHRR, utilizan las cinco bandas (b) disponibles: rojo o R (b1, centrado en 0.63 μm) que mide la radiación visible reflejada, la cual es atenuada débilmente por absorción de ozono y por la dispersión molecular y de partfculas; infrarrojo cercano o IRC (b2, centrado en 0.83 μm), el cual es moderadamente atenuado por vapor de agua y poco afectado por la dispersión molecular y de partfculas; infrarrojo térmico reflectivo o IRTR (b3, centrado en 3.7 μm); infrarrojo térmico 1 o IRT1 (b4, centrado en 10.8 μm) e infrarrojo térmico 2 o IRT2 (b5, centrado en 11.9 μm), los cuales detectan la radiación en el infrarrojo térmico, donde la atenuación es debida principalmente al vapor de agua, incrementando en resistencia al aumentarse la longitud de onda. El IRTR es también sensible a la radiación solar reflejada por la superficie terrestre y las nubes, y es sólo débilmente atenuado por la dispersión de partículas.

En el algoritmo modificado, todos los valores de reflectancias (R e IRC) son considerados como exo-atmosféricos y están dados en porcentajes. Los valores de temperaturas de brillantez (ecuación de Planck) IRT1 e IRT2 son también exo-atmosféricos y están dados en grados Kelvin (°K). La banda IRTR no es usada en la construcción del algoritmo para detectar píxeles claros, ya que los errores asociados a su calibración, en la generación de los sensores AVHRR de los ochentas y noventas, son muy altos.

La Figura 4 muestra el diagrama de flujo del algoritmo implementado en el sistema SPIA, donde la variable U representa valores umbrales para las pruebas implementadas.

El proceso dinámico de imágenes consistió en la asignación de códigos numéricos de la calidad de información disponible. La imagen reproyectada y referenciada es analizada píxel por píxel, para asignarle un código numérico. Los primeros códigos están asociados a la calidad de la información, por lo que es necesario revisar los indicadores de calidad contenidos en las imágenes AVHRR^(14,15), asf como verificarlos al procesar la información extraída de las imágenes.

Estimación del índice NDVIcp. A nivel de cada pfxel de la imagen satelital, se estimó el índice NDVIcp⁽⁴⁾, después de proyectar las reflectancias a nadir usando un modelo de geometría sol-sensor⁽¹¹⁾. Los parámetros empíricos del índice NDVIcp fueron optimizados para reducir efectos atmosféricos⁽¹⁸⁾. Por razones de eficiencia numérica, el sistema SPIA multiplica los valores del NDVIcp por 1000, por lo que en lo siguiente se considera esta situación.

Algoritmo de post-proceso

Primeramente, los datos generados por el sistema SPIA (a nivel pfxel por pfxel) día Juliano y NDVIcp son analizados si cumplen el criterio temporal de aceptación (X= día juliano) de XII < X < XFF (valores de defecto XII= 90 y XFF= 340; donde II significa inicio del ciclo de crecimiento y FF el final). Los datos que pasan el criterio, son analizados para determinar si su número N > NPZ (valor de defecto: NPZ= 10) es suficiente para continuar los análisis. Los datos que cumplen esta condición son ordenados de menor a mayor (base día juliano). Esto permite analizar la estructura temporal de los NDVIcp (denominados RY) en forma temporalmente ordenada. Se permiten días juliano con más de un dato (varias imágenes en el mismo día).

El siguiente paso es, determinar el valor máximo de RY (RYmax), asociado a Xmax, utilizando el criterio de búsqueda temporal XIM < X < XFM (valores de defecto: XIM= 180 y XFM= 334). Estos valores utilizan un conocimiento previo del pico de crecimiento de la vegetación. Con los datos ordenados (RY) en tiempo, i = 1,2,...,n, se analizan pares de datos contiguos para definir tres patrones básicos: decreciente (patrón 1), creciente (patrón 2) y plano (patrón 3), (Figura 5).

Para definir el criterio de aceptación de parejas contiguas de datos, se calcula el índice:

Así, la reglas de aceptación a los patrones básicos son:

Donde E es un valor decisión (valor de defecto: E= 0.2).

Como resultado final de este primer proceso de detección de patrones, todas las parejas contiguas de valores son etiquetadas con un patrón específico.

Análisis de patrones planos de la etapa inicial del crecimiento

El proceso inicial (primera revisión) general usado en todas las etapas de crecimiento de la vegetación (RY), consiste en el análisis de pares continuos de patrones, PP, para determinar si ambos pertenecen al patrón general buscado (patrón plano en este caso). Para esto, el análisis es realizado de los patrones etiquetados anteriormente (tres pares de datos a la vez). Los posibles patrones en PP(j-1) y PP(j) continuos (3²) están mostrados en la Figura 5.

De los patrones de la Figura 5, el Cuadro 1, muestra los valores de RY (y X) que son considerados como pertenecientes al patrón plano inicial de la curva de crecimiento (primera revisión), analizando los datos entre X = 2 y X = Xmax +1. En el caso de los patrones 1-2 y 2-1, Figura 5, se utiliza como criterio de aceptación:

De los valores aceptados en esta primera revisión, se verifica si el primer dato (primer día juliano de la serie original) está entre las parejas aceptadas. Si no está este dato, entonces se revisa si cumple la condición abs{RY(inicial de la serie aceptada)-RY(1)} < VCI (valor de defecto: VCI = 10). Si es así, este valor es agregado a la lista de parejas aceptadas. Con esto, se termina la primera revisión de posibles parejas de datos para caracterizar la etapa inicial (plana) del crecimiento de la vegetación.

La segunda revisión consiste en determinar qué valores usar de RY para calcular el promedio de estos como una aproximación al valor de RY1 (Y1 = 1/RY1). Para esto se revisa que no existan valores temporales mayores que Xmax + 1 y si es así, estos se rechazan. El criterio para considerar que un valor de RY contribuye al promedio es usar:

Donde: PRY representa un promedio calculado en forma progresiva en función del número de datos analizados y E1 es un valor que define el criterio de aceptación (valor de defecto: E1 = 0.1). Este proceso sesga el promedio (ecuación 10), a favor de los primeros valores de la serie temporal analizada. El último valor de RY aceptado de la serie, RYimi, es guardado para ser utilizado más adelante. Esto termina la segunda revisión de las posibles parejas, para caracterizar la etapa inicial (plana) del crecimiento de la vegetación.

La tercera revisión considera el promedio de RY, PRYI, calculado de los valores aceptados en la segunda revisión y utiliza el criterio:

Donde: E2 es el valor que define el criterio de aceptación (valor defecto: E2 = 0.08). Con esta última revisión se corrige el sesgo introducido en la segunda revisión. Los valores aceptados de RY son usados para calcular el promedio y con esto, estimar Y1 {= 1/promedio(RY)}.

En cualquiera de las revisiones, si no hay datos aceptados, el proceso define Y1 = -999 (dato faltante). Al existir un solo dato, se acepta éste para estimar Y1. La primera (Xic1, RYic1) y última (Xicn, RYicn) pareja de datos aceptados es guardada para análisis posteriores.

Análisis de patrones planos de la etapa final del crecimiento

La selección de parejas de datos para caracterizar esta etapa final plana del crecimiento de la vegetación, es similar al caso de la etapa inicial descrita ya. Los rasgos que la diferencian son: a) los valores analizados en esta etapa son los que van desde Xn-1 (Xn es el último valor de la serie original) hasta Xmax -1, en orden decreciente; b) se verifica que el último valor de la serie esté incluido en la primera revisión, usando como valor de aceptación a VCF (valor de defecto: VCF = 10), en vez de VCI; c) se guarda RYimf, último valor (primer valor), antes de la segunda revisión; d) la primera (Xfd1, RYfd1) y última (Xfdn, RYfdn) pareja de los datos aceptados es guardada para análisis posteriores; e) se calcula Y3 {= 1/ promedio (RY)}, de los valores aceptados para esta etapa.

Análisis de patrones de la etapa creciente

Se analizan las parejas de datos entre Xicn y Xmax + 1 para revisar los pares continuos de patrones, PP, para definir qué valores aceptar en la primera revisión del proceso. Esto es realizado de acuerdo al Cuadro 1. Para el par PP 2-3, se revisa que la ecuación 10 cumpla dRYP > E para aceptar el valor RY(i+1) en forma adicional a los inicialmente aceptados.

Posterior a la selección inicial de valores de patrón creciente, se revisa si el último valor del patrón plano inicial (Xicn, RYicn) está en la lista de pares aceptados. Si no es así, se agrega a la lista. Los valores aceptados, en pares contiguos, son revisados para ver si su patrón es creciente.

Adicionalmente, se revisa que:

como criterio de aceptación (que no estén muy debajo del valor del patrón plano inicial), usando la proporción DVD como valor de decisión (valor de defecto: DVD = 0.90).

Los valores de RY aceptados para la etapa creciente son transformados a valores Y (Y = 1/RY = 1/ NDVIcp) y analizados usando un proceso iterativo, regresión lineal, de eliminación de valores que no se ajusten al modelo lineal dado en la Figura 5a:

a) Se utiliza un criterio general de aceptación de valores que se ajustan a una línea recta (decreciente) por regresión lineal. El criterio de aceptación es que el coeficiente de correlación R² > R2U (valor de defecto: R2U = 0.8).

b) Cuando las parejas (X, Y) iniciales no cumplen el criterio de aceptación, se realiza un proceso iterativo no considerando una pareja a la vez y calculando el R² de este ejercicio. La pareja con R² mayor, después de no considerarla, es eliminada y el proceso es repetido hasta que el R² general (después de la eliminación) >R2U o queden sólo dos valores. Los parámetros A1 y B1 de la línea recta ajustada, son guardados para representar el patrón creciente (decreciente en Y), de acuerdo a la Ecuación (6).

c) Cuando sólo hay un par de valores, se asigna A1 = -999 y B1 = -999.

d) Se revisa que B1 < 0 para aceptar a B1.

El proceso definido determina la caracterización de la etapa creciente (en NDVIcp) de la curva de crecimiento de la vegetación.

Análisis de patrones de la etapa decreciente

La caracterización del patrón de la etapa decreciente de la vegetación sigue las mismas consideraciones que la etapa creciente, con las siguientes modificaciones:

• Se analiza la serie temporal que va de Xfd1 a Xmax -1, en forma decreciente

• La aceptación de los pares contiguos de patrones está definida en el Cuadro 1, sustituyendo 1 por 2.

• Se revisa si está el primer valor (en tiempo) del patrón plano final, Xfd1, y si no está, lo agrega a la serie aceptada del criterio del Cuadro 1 modificado. Algo similar hace con el último valor (en tiempo) del patrón creciente, Xmpc.

• Se revisa que no haya valores X < Xmpc y que sean decrecientes

• Se revisa si los valores cumplen el criterio:

para aceptarlos (valor defecto: DVC = 1.10).

• Los valores RY aceptados, son convertidos a Y, y analizados, usando una regresión lineal iterativa de los pares (X, Y), similar a la usada en la etapa creciente.

• Los parámetros A2 y B2 de la línea recta ajustada por regresión lineal son guardados y caracterizan el patrón decreciente (en NDVIcp) de la vegetación.

Análisis de patrones de la etapa de meseta del crecimiento

Las parejas entre Xicn y Xfd1 son analizadas para caracterizar el patrón plano de la meseta inicial, siguiendo un proceso similar al discutido para la etapa inicial y final del crecimiento de la vegetación. El promedio de los valores aceptados es usado para calcular Y2 {= 1/promedio (RY)}, si éste existe.

El valor estimado de la meseta, en caso de existir, es revisado en comparación de la intersección de las dos líneas rectas etapa decreciente y creciente en Y (Figura 5a).

Cuando los parámetros de las lfneas rectas (A1, B1) y (A2, B2) existen, la intersección es calculada de:

Si Y2 y Y2int existen (y son mayores que cero), entonces se revisa Y2 < Y2int y si esto sucede, Y2 = Y2int. De otra manera, Y2 es usada.

Aplicación del algoritmo de detección de patrones

El proceso de caracterización de los patrones de la curva de crecimiento de la vegetación, en base anual píxel por píxel, se aplicó a todos los píxeles en los tipos de vegetación de pastizales y matorrales del país y, en cada uno, se obtienen ocho valores que definen los parámetros de la curva de crecimiento anual: Y1, Y2, Y2int, Y3, A1, B1, A2, B2. Los parámetros X1, X2i, X2f y X3 se obtienen de los parámetros evaluados.

En la Figura 6a se muestra un caso extremo, típico de valores temporales de 1/NDVIcp, usando varios satélites del sensor AVHRR. Los valores altos de 1/NDVIcp corresponden a situaciones con efectos atmosféricos residuales (píxeles contaminados con nubes). Los datos de la Figura 6a fueron analizados con el algoritmo desarrollado y la Figura 6b muestra los resultados obtenidos. Se observa en la Figura 6b que los valores aceptados en cada etapa, caracterizan en forma adecuada (para los valores de defecto usados), la curva de crecimiento de la vegetación.

La Figura 7a muestra los valores de 1/NDVIcp para un píxel, caso extremo típico, del municipio de Hidalgo, Durango. Aplicando el algoritmo desarrollado a la Figura 7a, se obtienen los resultados mostrados en la Figura 7b. El modelo ajustado en la figura mencionada muestra el caso donde sólo hay dos datos aceptados para caracterizar el patrón decreciente de la curva de crecimiento de la vegetación. Asimismo, esta figura muestra que el patrón inicial del crecimiento resulta más difícil de caracterizar por falta de datos (se eliminan muchos de ellos en el algoritmo de filtrado de patrones).

La utilización de la tecnología de sensores remotos para la estimación de las curvas de crecimiento de la vegetación usando índices espectrales, permite tener una aproximación al estudio de los procesos de cambio, climáticos o de manejo, a una escala de píxeles individuales de imágenes satelitales.

La obtención de curvas paramétricas espectrales del crecimiento de la vegetación en términos de archivos históricos, permite contar con información para reconstruir el pasado del manejo del agostadero (pastizales y matorrales) del país, al menos desde 1986 a la fecha (sensor AVHRR en una malla nacional con celdas de 100 ha) o del 2000 a la fecha (sensor MODIS, con celdas de 6.25 ha). El uso de estas series de tiempo de los patrones anuales del crecimiento del agostadero, es analizado en otros trabajos de los autores, para estudiar los efectos del cambio climático y las prácticas de manejo asociadas a las actividades ganaderas en el país.

 

LITERATURA CITADA

1. Tucker CJ. Red and photographics infrared linear combination for monitoring vegetation. Remote Sens Environ 1979;8(2):127-150.         [ Links ]

2. Verstraete MM, Pinty B. Designing optical spectral indexes for remote sensing applications. IEEE Trans Geosci Remote Sens 1996;34(5):1254-1265.         [ Links ]

3. Gilabert MA, González J, García FJ, Meliá J. A generalized soil-adjusted vegetation index. Remote Sens Environ 2002;82(2-3):303-310.         [ Links ]

4. Paz F, Palacios E, Bolaños M, Palacios LA, Martínez M, Mejía E, Huete AR. Diseño de un índice espectral de la vegetación: NDVIcp. Agrociencia 2007;41(5):539-554.         [ Links ]

5. Huete AR. A Soil-Adjusted Vegetation Index (SAVI). Remote Sens Environ 1988;25(3):295-309.         [ Links ]

6. Goudriaan J, van Laar HM. Modelling potential crop growth processes, textbook with exercises. Current issues in production ecology. Dordrecht: Kluwer Academic Publishers; 1994.         [ Links ]

7. Stockle CO, Donatelli M, Nelson R. CropSyst, a cropping systems simulation model. Eur J Agron 2003;18(3-4):289-307.         [ Links ]

8. Huete AR, Hua G, Qi J, Chehbouni A. Normalization of multidirectional red and nir reflectances with SAVI. Remote Sens Environ 1992;41 (2-3):143-154.         [ Links ]

9. Chehbouni A, Kerr YH, Qi J, Huete AR, Sorooshian S. Toward the development of a multidirectional vegetation index. Water Resour Res 1994;30(5):1281-1286.         [ Links ]

10. Qi J, Huete AR, Cabot F, Chehbouni, A. Bidirectional properties and utilization of high-resolution spectra from a semiarid watershed. Water Resour Res 1994;30(5):1271-1279.         [ Links ]

11. Bolaños M, Paz F, Palacios E, Mejía E, Huete A. Modelación de los efectos de la geometría sol-sensor en la reflectancia de la vegetación. Agrociencia 2007;41(5):527-537.         [ Links ]

12. Jacquemoud S, Bacour C, Poilve H, Frangi JP. Comparison of four radiative transfer models to simulate plant canopy reflectance: direct and inverse mode. Remote Sens Environ 2000;74(3):471-481.         [ Links ]

13. Holechek JL, Pieper RD, Herbel CH. Range management, principles and practices. Englewood Cliffs. NJ: Prentice Hall; 1989.         [ Links ]

14. Kidwell KB editor. NOAA polar orbiter data users guide (TIROS-N, NOAA-6, NOAA-7, NOAA-8, NOAA-9, NOAA-10, NOAA-11, NOAA-12, NOAA-13 and NOAA-14). NOAA-NESDIS-NCDC. Suitland, MD. 1998.         [ Links ]

15. Goodrum G, Kidwell KB, Winston W, editors. NOAA KLM user's guide. NOAA-NESDIS-NCDC. Suitland, MD. 2001.         [ Links ]

16. Palacios LA, Pascual F, Paz F. Modelo de navegación satelital de imágenes de satélite multiespectrales para el monitoreo de zonas de riego. Congreso Nacional de Irrigación (CD-ROM). Asociación Nacional de Especialistas en Irrigación. Acapulco, Guerrero. 2005:9.         [ Links ]

17. Stowe LL, Davis PA, McClain EP. Scientific basis and initial evaluation of the CLAVR-1 global clearlcloud classification algorithm for the Advanced Very High Resolution Radiometer. J Atmos Oceanic Technol 1999;16(6):656-681        [ Links ]

18. Paz F, Bolaños M, Palacios E, Palacios LA, Martinez M, Huete AR. Optimización del índice espectral de la vegetación NDVIcp. Agrociencia 2008;42(8):925-937.         [ Links ]