Introducción
El índice de área foliar (IAF) proporciona información acerca de la cantidad de superficie fotosintética del bosque en relación con la superficie total (Weiss, Barret, Smith, Jonckheere y Coppin, 2004), la cual afecta las condiciones micrometeorológicas como la luz, agua, temperatura, humedad, entre otras y regula los intercambios de energía y de materia entre la vegetación y la atmósfera (Bréda, 2003; Chason, Baldocchi y Huston, 1991). Por lo tanto, el IAF tiene una relación directa con procesos vitales como la fotosíntesis, la respiración, la traspiración, la productividad y la intercepción de la precipitación (Schlerf, Atzberger y Hill, 2005). Cualquier cambio en el IAF, ya sea por heladas, tormentas, defoliación, sequía, u otros, ocasiona modificaciones en la productividad del rodal (Bréda, 2003). Consecuentemente, mapas precisos con la distribución espacial del IAF son necesarios para un mejor entendimiento de la estructura y función de los bosques y para cuantificar diferentes variables biofísicas y/o biológicas en ellos. Entre estas variables están la estructura, la composición vegetal, la producción primaria neta (PPN), el aumento de la biomasa epigea y sus implicaciones en los ciclos del carbono y otros nutrientes (Jarlan et al., 2008).
Existen distintas técnicas para la medición del IAF en campo (Weiss, Barret, Smith, Jonckheere y Coppin, 2004). Entre ellas se encuentran las técnicas directas, que generalmente implican la cosecha de las hojas o de hojarasca; técnicas mixtas que combinan información de campo con ecuaciones alométricas (Gower, Kucharik y Norman, 1999); y las técnicas indirectas, que utilizan equipos comerciales tales como LiCor LAI-2000 o TRAC, o bien usan fotografías hemisféricas para proporcionar información sobre la geometría del dosel que se asocia con el IAF (Rich, 1990). Sin embargo, los métodos de medición del IAF en campo no son prácticos para la estimación a escala de paisaje o a escalas espaciales grandes. Por lo que se han utilizado sensores remotos calibrados con una muestra de mediciones puntuales en campo para poder obtener así estimaciones continuas del IAF a escala de paisaje y/o región, mediante interpolación.
Uno de los métodos más utilizados para estimar el IAF y la fenología del dosel de los bosques, a partir de imágenes de satélite, consiste en asociar el IAF con características espectrales de la imagen y sus valores trasformados; particularmente, las bandas roja e infra-roja y/o los índices de vegetación, tales como el NDVI (normalized difference vegetation index, por sus siglas en inglés) o el EVI (enhanced vegetation index) (Tillack, Clasen, Kleinschmit y Förster, 2014, Aguirre-Salado et al., 2011, Morisette et al., 2006). Sin embargo, una fuerte limitación de este procedimiento es que los índices de vegetación y los valores espectrales de las bandas se “saturan” con altos niveles de biomasa, es decir, no varían aun cuando la biomasa continúe aumentando (Song, 2013). Esta situación limita fuertemente la posibilidad de generar mapas de IAF particularmente en bosques con altos niveles de biomasa, como los del área de estudio que oscilan entre 113 t ha-1 y 230 t ha-1 (Dai et al., 2014), los cuales son útiles para proponer acciones de manejo y conservación.
Un método potencial para superar los problemas de saturación de los índices de vegetación es utilizar imágenes satelitales de alta resolución, que permitan relacionar los datos de campo con información que incluya no solo datos espectrales de las imágenes, sino también información espacial de estas en diferentes formas, tales como análisis de textura, medidas de dependencia espacial (parámetros geoestadísticos de semivariogramas de las imágenes), entre otras (Gray y Song, 2012). Las imágenes con alta resolución espacial proporcionan potencialmente más información de textura que datos de mediana resolución porque permiten distinguir la estructura más fina en el bosque, por ejemplo, la estructura de la vegetación, la estructura de las copas, la altura de los árboles y el espaciamiento entre ellos. Esta estructura fina se manifiesta en imágenes con alta resolución espacial como variaciones en la tonalidad de la imagen (Zhou et al., 2014). Durante los últimos años, algunos estudios han demostrado el potencial que tienen las imágenes de alta resolución para estimar y mapear el IAF, particularmente en bosques de coníferas y selvas tropicales húmedas (Gray y Song, 2012; Beckschäfer, Fehrmann, Harrison, Xu y Kleinn, 2014; Zhou et al., 2014). Sin embargo, existen relativamente pocos estudios de este tipo en selvas secas, las cuales difieren en su fenología de las selvas húmedas y los bosques de coníferas.
Las selvas tropicales secas presentan una marcada fenología estacional como respuesta a la duración y severidad de la sequía. La pérdida de las hojas es una estrategia común para reducir la pérdida de agua (Poorter y Markesteijn, 2008). Muchas especies logran sobrevivir a largos periodos de estiaje gracias a la pérdida de hojas (Gao, Giese, Brueck, Yang y Li, 2013). En la mayoría de las especies vegetales de los bosques tropicales secos, el crecimiento, la reproducción, la germinación de las semillas y el establecimiento de las plántulas se encuentran limitados a la corta temporada de lluvias (Sánchez-Azofeifa et al., 2005). Por lo tanto, la cuantificación de los cambios temporales del IAF en las selvas tropicales secas resulta de importancia para lograr entender los procesos ecosistémicos (Maass, Vose, Swank y Martínez-Yrízar, 1995). Esto acentúa la necesidad de realizar análisis a partir de imágenes multi-temporales como en el presente trabajo.
A pesar de que los bosques tropicales caducifolios presentan una amplia distribución (42% de la extensión total de los bosques tropicales), los datos de IAF para selvas tropicales deciduas son escasos en comparación con los disponibles para selvas tropicales húmedas (Kalácska, Calvo-Alvarado y Sánchez-Azofeifa, 2005; Quesada et al., 2009). Hasta el año 2003 solo 8% de los valores de IAF incluidos en la base de datos mundial correspondían a regiones de bosques tropicales secos (Asner, Scurlock y Hicke, 2003). El estudio de los cambios estacionales en el IAF en las selvas tropicales secas es aún más escaso, a pesar de ser sumamente importante para el entendimiento de los procesos de intercambio de energía y materia en ese tipo de ecosistemas (Liu, Chen, Jin y Qi, 2015).
Objetivos
El estudio tuvo dos objetivos 1) cuantificar el IAF tanto para la estación de lluvias (mayor cantidad de tejido foliar) como para la de estiaje (menor cantidad de tejido foliar); y 2) Ajustar modelos que permitan mapear de manera continua el IAF en ambas estaciones (lluvia y estiaje) a partir de imágenes de alta resolución y que puedan ser utilizados en selvas tropicales secas con altos niveles de biomasa y características similares a las del sitio de estudio.
Materiales y métodos
Área de estudio
El área de estudio se localiza en el SO del estado de Yucatán, dentro de la Reserva Biocultural Kaxil Kiuic, A. C. (20°04’-20°08’ N, 89°32’-89°5’ W). El clima es tropical, cálido-subhúmedo (Aw) con lluvias en verano y sequía de noviembre a abril. La temperatura media anual es de 26 °C. La precipitación media anual se encuentra en el intervalo entre 1000 mm y 1200 mm, presentándose principalmente en los meses de junio a octubre (Hernández-Stefanoni et al., 2015). Topográficamente, se encuentra en una zona de lomeríos bajos con suelos de piedra caliza, pendiente moderada (10° - 25°) y baja elevación (60 m snm -160 m snm) que se intercalan con planicies de suelos profundos (Flores y Espejel, 1994). La vegetación presente es de selva mediana subcaducifolia (Flores y Espejel, 1994), en la que entre 50% y 75% de las especies pierden sus hojas durante la temporada de estiaje. La mayor parte del arbolado tiene una altura de entre 8 m y 13 m y una edad de 20 años a 90 años, aproximadamente a partir del abandono, tras ser sometida a la práctica agrícola tradicional de roza, tumba y quema. Las especies más abundantes en esta selva son Neumillspaughia emarginata (H. Gross) S.F. Blake, Gymnopodium floribundum Rolfe., Bursera simaruba (L.) Sarg., Piscidia piscipula (L.) Sarg. y Lysiloma latisiliquum (L.) Benth (Fig. 1).
Muestreo de datos en campo
Como parte de un estudio más amplio de dinámica del carbono forestal (Comisión Nacional Forestal [ Conafor], 2017), se establecieron 32 parcelas circulares de 400 m2, distribuidas de manera sistemática en un área de 3 km × 3 km alrededor de una torre de flujos que mide el intercambio neto de CO2 y agua entre la vegetación y la atmósfera (Fig. 1). La mayoría de las parcelas se encuentra dentro de la Reserva Biocultural Kaxil Kiuic. En cada parcela, con ayuda de un GPS, se tomaron las coordenadas centrales como referencia para el trazo y establecimiento de las parcelas. Se adquirieron fotografías hemisféricas en el centro de cada parcela, una vez durante el mes de octubre de 2014 y otra en abril de 2015, con la intención de obtener estimaciones del IAF en las épocas de mayor y de menor tejido foliar. Las fotografías hemisféricas se tomaron siguiendo el protocolo establecido por la Comisión Nacional Forestal y la Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (Conafor-Conabio) dentro del Inventario Nacional Forestal y de Suelos (Infys) (Comisión Nacional Forestal [Conafor], 2010). Una de las especificaciones más importantes es la que indica que las fotografías deben tomarse exclusivamente en horarios cercanos al amanecer o atardecer, cuando el sol se encuentre cercano al horizonte y no afecte las mediciones creando reflejos que subestimen o sobrestimen el IAF (Weiss, Barret, Smith, Jonckheere y Coppin, 2004; Beckschäfer, Seidel, Kleinn y Xu, 2013). Se optó por la opción de “bracketing”, consistente en capturar imágenes en campo con diferentes exposiciones, dado que esto permite seleccionar aquellas fotografías con el contraste suficiente para distinguir claramente entre cielo y vegetación (Sasaki, Imanishi, Ioki, Song y Morimoto, 2013; Frazer, Lertzman y Trofymow, 1997).
Cálculo del IAF a partir de fotografías hemisféricas
Las fotografías hemisféricas se procesaron en 2 fases: La primera consistió en un control de calidad para seleccionar la imagen con el mejor contraste entre el cielo abierto y la vegetación. En este estudio, tras una inspección visual, se seleccionaron las fotografías con exposición -1, las cuales proporcionaron el valor de IAF estimado más cercano al obtenido de manera directa e independiente a través de la medición de la producción anual de hojarasca (Nafarrate-Hecht; 2017). Posteriormente se convirtieron en imágenes en blanco y negro (proceso de binarización) con ayuda del software SideLook 1.1.01 (Nobis y Hunziker, 2005). Se eligió este software para el preprocesamiento de las fotografías hemisféricas debido a que permite la interacción de un usuario experimentado al momento de la asignación del umbral, lo cual hace posible obtener una imagen binaria fiel a la original, cosa que no sucede con el software Hemisfer 2.0.1.0, que asigna el umbral automáticamente, generando errores en un número importante de imágenes, debido a la presencia de nubes que son clasificadas como tejido foliar (Schleppi et al., 2014). La transformación de la imagen a color en una en blanco y negro es la parte más crítica del proceso de análisis (Jonckheere et al., 2004; Beckschäfer, Seidel, Kleinn y Xu, 2013) y de ella dependerá la precisión del IAF calculado. Algunos autores aseveran que debe ser el algoritmo del software especializado el que determine si los pixeles serán blancos (cielo) o negros (vegetación) (Jonckheere et al., 2004; Weiss, Baret, Smith, Jonckheere y Coppin, 2004; Beckschäfer, Seidel, Kleinn y Xu, 2013; Schleppi et al., 2014) para así asegurar un proceso repetible y sin sesgo. Sin embargo, esta selección automática no siempre genera una contraparte binaria fiel (Frazer, Lertzman y Trofymow, 1997; Nobis y Hunziker, 2005), sobre todo cuando se encuentran nubes o sitios con poca vegetación. En este estudio, se optó por una selección manual a criterio de una experta en el procesamiento de fotografías hemisféricas del Infys (la primera autora). En la segunda fase se analizaron estas imágenes binarias con el software Hemisfer 2.0.1.0 para obtener el IAF, encontrando que la ecuación propuesta por Chen y Cihlar (1995) generó los mejores valores.
Procesamiento de las imágenes de satélite
Se adquirieron dos escenas de imágenes de satélite RapidEye con una resolución espacial de 5 m, una del 1 de febrero de 2013 (la fecha más cercana a la época de lluvias que se pudo obtener con una cobertura libre de nubes y cuando la mayor parte de la vegetación aún estaba verde) y otra del 31 de marzo del 2013, durante la estación seca. Las imágenes fueron corregidas radiométrica y atmosféricamente. Los valores de los pixeles en las 32 parcelas, que correspondían a la banda roja (630 nm - 685 nm) y borde rojo (690 nm - 730 nm) se identificaron y fueron extraídos. Las bandas: roja, infrarrojo-cercano y/o infrarrojo-medio han sido recomendadas para discriminar la vegetación y se asocian con el IAF (Nagendra, 2001). En este estudio se incluyó la banda del borde rojo, en lugar de la banda infrarroja-cercana dado que tiene mayor potencial de discriminar la vegetación (Kross, McNairn, Lapen, Sunohara y Champagne, 2015).
El NDVI del borde rojo y el EVI se calcularon y extrajeron en las 32 parcelas de muestreo de cada estación del año. El segundo índice presenta menor saturación en condiciones de dosel denso, debido a que es más sensible a las características estructurales de la vegetación (Gao, Huete, Ni y Miura, 2000), además, incluye las constantes C1 = 6 y C2 = 7.5, que son coeficientes asociados a la resistencia de aerosoles y a la banda azul (440 nm - 510 nm) utilizada para la corrección atmosférica. También incorpora L, la minimización de la sensibilidad de los valores de reflectancia a variaciones del suelo, así como un factor de ganancia (G = 2.5) (Gao, Huete, Ni y Miura, 2000; Jiang, Huete, Didan y Miura, 2008). Estos dos índices se calcularon con GRASS GIS 7.0 (GRASS Development Team, 2015), utilizando las siguientes fórmulas:
Las medidas de textura se definen como la cuantificación de la variabilidad en los valores de reflectancia de los píxeles vecinos y su disposición espacial en un área determinada (Haralick y Shanmugam, 1973). Se calcularon dos métricas de textura de primer orden y doce de segundo orden para ambas bandas (borde rojo y rojo) y para los índices de vegetación (NDVI y EVI) usando GRASS GIS 7.0 (GRASS Development Team, 2015). Las medidas de textura de primer orden se calcularon a partir de los valores de imágenes originales que se obtuvieron dentro de un determinado tamaño de ventana y no consideran las relaciones de píxeles vecinos. En este grupo se encuentran el promedio de los pixeles en la ventana (mean) y la varianza de los mismos (nvar). Por otro lado, las medidas de textura de segundo orden se calculan con base en una matriz de coocurrencia de niveles grises (también llamada matriz de dependencia espacial de tonos de grises) y consideran las relaciones entre píxeles vecinos (Haralick y Shanmugam, 1973). Estas medidas se pueden agrupar en tres categorías, la primera se basa en el grado de contraste entre píxeles e incluye el momento u homogeneidad de la diferencia inversa (IDM) y el contraste (CONTR); el segundo se basa en la organización de píxeles dentro de una ventana (segundo momento angular -ASM), mientras que la tercera categoría está basada en datos estadísticos como el promedio de la suma (SA), varianza (VAR), varianza de la suma (SV), varianza de la diferencia (DV), entropía (ENTR), entropía de la suma (SE) y entropía de la diferencia (DE), correlación ( CORR) y medidas de correlación (MOC 1). Se calcularon las medidas de segundo orden a 0, 45, 90 y 135 grados y se promediaron para obtener un valor único de textura sin dirección (Haralick y Shanmugam, 1973). En la tabla 1 se muestra una lista de las métricas de textura empleadas. Las métricas de textura de primer y segundo orden se calcularon usando dos tamaños de ventana de 3 × 3 y 5 × 5 píxeles dado que las imágenes RapidEye tienen una resolución espacial de 5 m, por lo que con dichos tamaños de ventana se representarían 225 m2 y 625 m2, respectivamente, que están apenas por abajo y por arriba del tamaño de la parcela: 400 m2.
Medida (Abreviación) | Fórmula | Descripción |
---|---|---|
Medidas de primer orden | ||
Mean (mean) | Valor medio de los pixeles en la ventana. | |
Variance (nvar) | Valor de varianzade los piexeles en la ventana. | |
Medidas de segundo orden | ||
Sum Average (SA) | Σ i=2 2Ng i px+ y (i) | Valor medio de la suma de los valores de pixeles vecinos. |
Entropy (ENTR) | Σ i Σ j p(i , j) log( p(i , j)) | Medida de desorden o aleatoriedad. Valores altos corresponden a mayor complejidad en la imagen. |
Sum Entropy (SE) | Σ i=2 2 Ng p x+ y (i) log {p x+ y (i) } | Medida de entropía de la suma de los valores de pixeles vecinos. Valores altos corresponden a mayor complejidad en la imagen. |
Difference Entropy (DE) | −Σ i=0 Ng−1 p x−y (i) log {p x−y (i) } | Medida de entropía de la diferencia de los valores de pixeles vecinos. Valores altos corresponden a mayor complejidad en la imagen. |
Sum of Squares: Variance (VAR) | Σ j Σ j (i−μ) 2 p(i , j) | Variabilidad o "Varianza" en niveles de grises en la imagen. |
Sum Variance (SV) | Σ i=2 2 Ng (i− var) 2 p x+ y 2 | Varianza en niveles de grises para la suma de los valores de pixeles vecinos. |
Difference Variance (DV) | varianza of p x - y | Varianza en niveles de grises para la suma de los valores de pixeles vecinos. |
Angular Second Moment (ASM) | Σ i Σ j [ p(i , j)] 2 | Segundo momento angular o “Uniformidad” medida de la homogeneidad local, opuesto de Entropía. |
Contrast (CONTR) | Σ n=0 N g−1 n 2 [Σ i=1 Ng Σ j=1 Ng p(i , j)] | Variaciones locales de niveles de grises (medidas del contraste de la imagen). |
Inverse Difference Moment (IDM) | Σ i Σ j (1/1+(i+ j) 2 )p(i , j) | “Homogeneidad” medida directa de homogeneidad local. |
Correlation (CORR) | (Σ i Σ j (ij) p(i , j)−μ x−μ y )/(σ x σ y) | Mide la dependencia lineal en valores de niveles de grises. |
Information Measures of Correlation (MOC) | HXY−HXY 1 / max {HX ,HY } | Mide la dependencia lineal en valores de niveles de grises en los valore de pixeles vecinos. |
Notación: HXY= −Σi Σj p(i , j)log( p(i , j)); HXY1= −Σi Σj p(i , j) log {px(i) py(j)}; HXY2 = −ΣiΣj px(i) py(j) log {px(i) py(j); pi : frecuencia del nivel de gris i; pj : frecuencia del nivel j; p(i,j): matriz de frecuencias relativas de la entrada (i,j) en una matriz co-ocurrencia; p(i,j): matriz de frecuencias relativas (i,j) entrada en una matriz co-ocurrencia; px(i): i-ésima entrada en una matriz de probabilidades marginales obtenido sumando las filas (i,j); Ng: niveles de gris.
Análisis de datos
Dado que uno de los objetivos de esta investigación fue ajustar modelos para predecir el IAF usando como variables explicativas los valores de reflectancia, índices de vegetación y métricas de textura, se utilizó un modelo que combina la regresión lineal múltiple con kriging ordinario de los residuales (obtenidos de la regresión). Este método considera tanto la dependencia espacial de las observaciones (residuales de la regresión), como las relaciones entre la variable dependiente (IAF) y las variables explicativas (obtenidas de las imágenes) (Webster and Oliver, 2001).
El estimador de regresión con kriging del IAF Zrk(x) se define como la suma de la estimación de la regresión Zr(x) obtenida como una función lineal entre el IAF y las variables explicativas (refectancia de las banda, índices de vegetación y textura de la imagen), más los residuales de este modelo de regresión que se estiman utilizando el kriging ordinario Єok(x), según la siguiente ecuación (Hernández-Stefanoni, Gallardo-Cruz, Meave y Dupuy, 2011):
Para hacer la estimación de kriging con regresión, primero se utilizó un modelo de regresión lineal múltiple, en donde los valores la refectancia de la banda roja, la banda del borde rojo, el NDVI, el EVI y 52 medidas de textura (14 medidas de textura de la imagen calculadas para cada una de las siguientes capas: la banda roja, la banda del borde rojo, el NDVI y el EVI) se usaron como variables explicativas y el IAF como variable dependiente. Las variables explicativas se transformaron según fue necesario con 1 / x, log10 (x), log10 (x + 1) o sqrt (x) para cumplir con los supuestos de linealidad. Se utilizó el procedimiento llamado bestsubset, particularmente la función regsubsets en el paquete ‘leaps’ del software R (R Development core team, 2012) para seleccionar el mejor modelo de todos los subconjuntos posibles de variables explicativas. En este caso, se seleccionaron cinco modelos candidatos para los datos de cada temporada, que incluía una, dos, tres, cuatro y cinco variables explicativas. Posteriormente, se utilizó el menor valor de akaike AIC para seleccionar el mejor modelo. Dado que la multicolinearidad entre las variables predictoras puede causar problemas en el modelado, se verificó que las variables explicativas consideradas para el análisis no estuvieran correlacionadas o expresaran una colinealidad pequeña, con un factor de inflación de varianza menor a 2. Posteriormente, los residuales del modelo de regresión lineal anterior, se utilizaron para calcular un semivariograma experimental, utilizando la función gstat de R para ajustar modelos esféricos, gaussianos, lineales y exponenciales a los semiovariogramas experimentales (Pebesma, 2004). El coeficiente de determinación (R2) resultante del ajuste de los modelos fue el criterio utilizado para seleccionar el mejor modelo de semivariograma.
Para validar los modelos de predicción del IAF por medio del método de regresión con kriging, y evaluar su poder predictivo se utilizó una validación cruzada. En este procedimiento, una observación se elimina temporalmente del conjunto de datos y los valores de muestreo restantes se utilizan para predecir la variable. Esta validación cruzada arroja una lista de valores estimados de IAF, que se pueden relacionar con los obtenidos en el campo. El poder de predicción de los modelos se evaluó a través del coeficiente de determinación (R2), obtenido de una regresión simple de los valores observados en campo y los estimados, y del error cuadrático medio (RMSE), que es la raíz cuadrada del promedio de las diferencias al cuadrado entre los valores estimados y predichos.
Finalmente, utilizando las ecuaciones de regresión ajustadas parada cada temporada, así como las capas con las bandas de reflectancia, los índices de vegetación, la textura de las imágenes y los residuales estimados con kringing , se generaron los mapas con la distribución espacial del IAF para las temporadas de lluvia y estiaje usando algebra de mapas. Estos mapas fueron clasificados en categorías de IAF para su visualización.
Resultados
Índice de área foliar obtenido en campo
En la temporada de lluvias el IAF tuvo una media de 3.37 y una desviación estándar de 0.44, mientras que para el estiaje la media fue de 2.49 con 0.38 de desviación estándar. El IAF fue significativamente menor durante la época de sequía que durante la de lluvias (p < 0.05) (Fig. 2) y la magnitud de las diferencias varió entre parcelas entre 4.8% y 48.7%.
Modelización de la distribución espacial del IAF
En la tabla 2 se muestran los modelos candidatos de las relaciones entre el IAF y las variables explicativas, así como el mejor modelo para cada temporada basado en el AIC y la multicolinearidad, mientras que en la tabla 3 se muestra el resumen de las estadísticas de la regresión entre el IAF y las variables explicativas, solamente para el tamaño de ventana de 5 × 5, que presentó los mejores resultados. Los resultados de la regresión múltiple indicaron relaciones estadísticamente significativas entre el IAF y la reflectancia de las bandas y las métricas de textura de la imagen, tanto para la temporada de lluvias, como para la de estiaje. En el primer modelo de regresión (temporada de lluvias), el IAF se relacionó negativamente con las medidas de textura y de manera positiva con la media de la banda 4 (borde del rojo). Por otro lado, en la temporada de estiaje, el IAF se asoció negativamente con la textura y con la banda 3 (banda roja). El porcentaje de variación explicado por los modelos de regresión fue ligeramente mayor en la temporada de estiaje R2 aj = 0.63 comparado con la temporada de lluvias R2 aj = 0.58 (Tabla 3).
Temporada | Núm de variables | Variables explicativas | AIC |
---|---|---|---|
Lluvias | 1 | MOC1_evi | -63.11 |
2 | MOC1_evi, Corr_band4 | -67.83 | |
3 | MOC1_evi, nvar_band4, Corr_band4 | -70.32 | |
4 | MOC1_evi , band4, mean_band4, Corr_band4 | -77.86 | |
5* | ndvi, Var_ndvi, Var_evi, MOC1_evi, Corr_band4 | -82.88 | |
Estiaje | 1 | MOC1_band4 | -65.61 |
2 | DV_ndvi, SE_evi | -69.84 | |
3 | DV_ndvi, Contr_ndvi, SE_evi | -73.75 | |
4 | DV_ndvi, Contr_ndvi, SE_evi, MOC1_band4 | -82.67 | |
5 | DV_ndvi, Contr_ndvi, SE_evi, Corr_band4, band3 | -89.53 |
* indica que existe Multicolinearidad entre variables explicativas.
El mejor modelo de cada temporada se indica en negrillas.
AIC = criterio de información de akaike, band3 = banda del rojo, band4 = banda del borde rojo, ndvi = índice normalizado de vegetación, evi = índice mejorado de vegetación, Contr = contraste, mean = media de los pixeles, nvar = varianza de los pixeles, Corr = correlación, SE = suma de la entropía, DV = diferencia de la varianza, MOC1 = medida de información de correlación
Temporada | Variables del Modelo | Estimativas de parámetros(Error estándar) | Coeficientes Estandarizados | R2 | R2 Ajustada | Error estándar de la estimación |
---|---|---|---|---|---|---|
Lluvias | Constante | -1.07(1.07) | 0.64 | 0.58 | 0.27 | |
MOC1_evi* | -5.316(1.07) | -0.650 | ||||
band4* | -0.001(0.0001) | -1.145 | ||||
mean_band4* | 0.002(0.0001) | 1.073 | ||||
Corr_band4* | -1.208(0.31) | -0.475 | ||||
Estiaje | Constante | 11.11(1.43) * | 0.69 | 0.63 | 0.22 | |
DV_ndvi* | -541.169(94.35) | -1.022 | ||||
Contr_ndvi* | -0.047 (0.01) | -0.811 | ||||
SE_evi* | -5.210 (0.98) | -0.686 | ||||
Corr_band4* | 0.705(0.19) | 0.429 | ||||
band3* | 0.0001(0.00001) | -0.346 |
* Variables incluidas en el modelo con p < 0.01; band3 = banda del rojo, band4 = banda del borde rojo, ndvi = índice normalizado de vegetación, evi = índice mejorado de vegetación, mean = media de los pixeles, Corr = correlación, Contr = contraste, DV = diferencia de la varianza, SE = suma de la entropiaa, MOC1 = medida de información de correlación
Los modelos de semivariogramas revelaron una estructura espacial en los residuales del IAF para la temporada de lluvias, mientras que para la temporada de estiaje no hubo autocorrelación espacial (Tabla 4, Fig. 3). En la temporada de lluvias se ajustó un modelo esférico a los semivariogramas experimentales, produciendo valores de R2 = 0.89. La varianza estructural, que determina la varianza debida a la dependencia espacial explicada por el modelo, fue de 47.9%. Este resultado sugiere que el modelo incorpora una fracción moderada de variabilidad residual inexplicada (Fig. 3). El rango de influencia mostró valores de 1366.0 m, lo que indica que se esperaría razonablemente una dependencia espacial de los valores residuales entre parcelas separadas por distancias de hasta 1.37 km. En la temporada de estiaje se ajustó un modelo lineal en donde las varianzas nugget y sill tienen el mismo valor y la variación explicada es por el modelo es 0.0, lo cual indica que no existe dependencia espacial de los residuales.
Temporada | Modelo | Varianza nugget | Varianza total | Rango | Varianza estructural relativa (%) | R2 |
---|---|---|---|---|---|---|
Lluvias | Esférico | 0.0360 | 0.0697 | 1366.0 | 47.9 | 0.89 |
Estiaje | Lineal | 0.0330 | 0.0330 | 1267.8 | 0.0 | 0.62 |
La validación cruzada de las predicciones del IAF usando regresión con kriging (Fig. 4) muestra que la autocorrelación espacial en los residuales del IAF mejoró ligeramente la variación explicada por el modelo en la temporada de lluvias en comparación con el análisis que utiliza regresión de manera exclusiva. Los valores de R2 aumentaron de 0.49 a 0.51 y el RMSE disminuyó de 0.30 a 0.28. Por otro lado, los valores de R2 fueron mayores en la temporada de estiaje comparados con los de la temporada de lluvias (0.56 y 0.51 respectivamente). De igual manera los valores de RMSE fueron ligeramente menores en la temporada de estiaje que en la de lluvias (0.25 y 0.28, respectivamente).
Finalmente, los mapas de distribución espacial del IAF para la temporada de lluvias y de estiaje se muestran en las figuras 5, respectivamente. En estos mapas se puede observar que el IAF es más uniforme o con menores variaciones durante la temporada de estiaje.
Discusión
Los resultados obtenidos en este estudio revelan claramente una variación estacional en el IAF; en la temporada de lluvias se obtuvo un valor promedio de 3.37, mientras que en la temporada de estiaje fue ligeramente menor con 2.49. Estas variaciones estacionales han sido registradas en diferentes investigaciones, por ejemplo, Maass et al., (1995) y Kalácska et al. (2005), quienes encontraron variaciones estacionales en el IAF en selvas tropicales secas en varias regiones de Latinoamérica, indicando que esta diferencia se debe a los tiempos de desprendimiento de las hojas según la forma de crecimiento (árboles, arbustos o lianas) y a la edad de sucesión. Esta asincronía en el IAF tiene implicaciones sobre la estimación del IAF a escala regional, en particular cuando se obtienen a partir de datos de sensores remotos. De ahí la necesidad de realizar análisis con imágenes multitemporales como en este estudio.
Los análisis de regresión múltiple mostraron diferencias en el porcentaje de variación explicado por el modelo con R2 = 0.64 y 0.69 (R2 ajustada = 0.58 y 0.63) para la temporada de lluvias y estiaje, respectivamente. En este estudio fue mejor el modelo encontrado en la época de estiaje, incluso la validación del modelo muestra un mayor poder de predicción en la temporada de estiaje con un valor de R2 = 0.56 y un RMSE = 0.25, comparado con R2 = 0.51 y un RMSE = 0.28 en la temporada de lluvias. Estas variaciones se deben a diferentes factores, entre los que se encuentran la apertura del dosel, la presencia o ausencia de tejido foliar, la estructura de la vegetación, entre otras (Colombo, Bellingeri, Fasolini y Marino, 2003). La información de textura extraída de imágenes de satélite ha demostrado mejorar la precisión de las estimaciones del IAF (Gray y Song, 2012; Zhou et al., 2014), dado que la información de la textura de los pixeles y la resolución espacial alta de la imagen se asocian con los atributos estructurales del bosque. Además, se ha demostrado que durante la estación seca puede ser más exitoso en el mapeo de la vegetación debido a que las características de la estructura de la vegetación son más pronunciadas y distintivas durante este periodo, ya que hay pocas hojas (Kalacska et al., 2007). Sin embargo, los sensores remotos que tienen resolución espacial alta para estimar el IAF a escala local carecen de la resolución temporal para mapear la evolución de la vegetación en diferentes estaciones, considerando además que la cobertura de las nubes reduce el número de imágenes que se pueden usar en un año, particularmente en los trópicos. Por lo tanto, se sugiere que, en estudios posteriores, se exploren técnicas de fusión de imágenes de múltiples sensores, para combinar las características deseables de resolución espacial con las de resolución temporal (Ghassemian, 2016).
Por otro lado, los resultados obtenidos en este estudio respaldan trabajos previos, en donde información de textura extraída de imágenes de alta de resolución junto con información espectral de la imagen caracterizan de mejor manera el IAF, en comparación con el uso de información espectral de manera exclusiva (Zhou et al., 2014, Beckschäfer, Fehrmann, Harrison, Xu y Kleinn, 2014, Pu y Cheng, 2015). Por ejemplo, Zhou et al. (2014) encontraron que los modelos de estimación del IAF basados en reflectancia explicaron 68% de la variación observada en el IAF, mientras que los métodos basados en textura de manera exclusiva explicaron hasta 72% de la variación observada. Sin embargo, la combinación de ambas variables explicativas tuvo valores de R2 de 0.84. No obstante, los valores de R2 obtenidos en el presente estudio fueron más bajos que la mayoría de los valores encontrados en estudios similares, posiblemente debido a la poca heterogeneidad entre sitios con altos niveles de biomasa (Tillack, Clasen, Kleinschmit y Förster, 2014; Pu y Cheng 2015). Los bajos valores de R2 obtenidos en este estudio también podrían deberse en parte a diferencias en las condiciones ambientales y, por consiguiente, en la fenología de la vegetación entre la toma de las fotografías hemisféricas en campo y las imágenes satelitales adquiridas. Esto también podría explicar en parte la menor precisión y capacidad de predicción del modelo de la temporada de lluvias comparado con el de sequía, dado que no se pudo adquirir una imagen satelital durante la temporada de lluvias, sino al comienzo de la temporada de sequía (febrero), cuando varios árboles podrían haber perdido sus hojas.
Uno de los principales resultados de este estudio es que el IAF varía no solo debido a la fenología del bosque, sino también debido a la relación entre las mediciones de campo y las estimaciones basadas en imágenes de satélite que cambian entre estaciones, generándose modelos con diferentes variables explicativas. En general, el IAF está fuertemente relacionado con los datos de reflectancia y textura obtenidos con imágenes de alta resolución. Más específicamente, el análisis de regresión sugiere que el IAF está fuertemente relacionado con la reflectancia de las bandas 3 (roja) y 4 (borde rojo) y varios atributos de textura. El coeficiente negativo de la banda 3 indica una respuesta baja en las longitudes de onda rojas visibles debido a la absorción de clorofila y su consiguiente alto potencial para discriminar el IAF en función de las diferencias en sus propiedades de reflectancia. Por otro lado, el coeficiente positivo con el promedio de la banda 4 indica una respuesta alta en las longitudes de onda del borde rojo causada por la dispersión de la luz en el mesófilo de las hojas (Nagendra, 2001; Thenkabail, Enclona, Ashton, Legg y De Dieu, 2004). Además, se encontraron relaciones negativas entre la heterogeneidad espectral y el IAF (Tabla 2), esto indica asociaciones del IAF con texturas homogéneas. Las medidas de textura incluyen información sobre los atributos estructurales del bosque tales como el diámetro de la copa, la densidad de árboles, el tamaño y la distribución de los huecos del dosel, y la iluminación/geometría del mismo. Se espera que, en áreas con una gran cantidad de hojas pero con dosel cerrado, las texturas sean homogéneas (Gray y Song, 2012). Este parece ser el caso del sitio de estudio, en donde existen áreas de bosque muy cerrado, con árboles de muchas hojas. En contraste, las texturas heterogéneas se asociarían con áreas con niveles de tejido foliar intermedio y un dosel más heterogéneo, es decir con aperturas (Gray y Song, 2012). Por otro lado, durante la estación seca se encontró que no existía dependencia espacial en los residuales, lo cual podría deberse a que la caída de las hojas crea áreas más homogéneas en relación con el IAF; dicho de otra manera, existe menos variabilidad o mayor uniformidad en la distribución espacial del IAF durante esta temporada, comparada con la de lluvias, como se observa en la figura 5 (Dale y Fortin, 2002).
Cuando se toma en cuenta la autocorrelación espacial de los puntos de muestreo, se puede mejorar la precisión de las predicciones. En este estudio, la predicción del IAF durante las lluvias mejoró ligeramente al tomar en cuenta la dependencia espacial de los residuales (Fig. 5). Incorporar esta dependencia espacial entre observaciones no solo puede mejorar la exactitud de las estimaciones, sino que también proporciona información sobre la estructura espacial del IAF, más allá de la asociación obtenida entre este índice y las propiedades espectrales y de textura del bosque. (Pfeffer, Pebesma y Burrough, 2003).
Conclusiones
El procedimiento utilizado para el mapeo del IAF en bosques tropicales secos con alto contenido de biomasa para diferentes temporadas es práctico y se complementa con la integración de valores de reflectancia, índices de vegetación y textura de las imágenes de satélite y las mediciones de campo del IAF, además de la estructura espacial de los datos. Los modelos empíricos pueden usarse para predecir el IAF a escala de paisaje en áreas de bosques tropicales secos de la península de Yucatán y de otras regiones con características similares, dado que tuvieron buenos ajustes y predicción. Adicionalmente cuando existe una autocorrelación espacial entre las parcelas de muestreo, la contabilización de esta estructura espacial mejora la predicción del IAF usando modelos que incluyen dicha estructura.
Se confirma que existen variaciones estacionales en el IAF en este tipo de bosques y que estas influyen de manera directa en las estimaciones basadas en imágenes de satélite. También se encontró que información de textura extraída de imágenes de alta de resolución junto con información espectral de la imagen caracterizan adecuadamente el IAF en bosques con alto contenido de biomasa.