INTRODUCCIÓN
Las plantaciones forestales ocupan un área aproximada de 53.4 millones de hectáreas, de las cuales el 46.4 % se distribuyen en América Latina () y representan la principal fuente de celulosa del mundo, lo que las convierte en determinantes para las industrias del papel, cartón, energía eléctrica, conglomerados y de mobiliario, entre otras (). Adicional a esto, estas plantaciones juegan un papel importante en la conservación de los bosque naturales al reducir su presión, mejorar la conectividad del paisaje y la biodiversidad, así como secuestrar carbono y restaurar las tierras degradadas y servicios ecosistémicos ().
La biomasa aérea (AGB, por su siglas en inglés) es un indicador de la productividad de las plantaciones forestales que permite conocer su estado de desarrollo y el potencial de rendimiento en volumen (m3.ha-1; ). Los métodos convencionales para estimar AGB de plantaciones forestales se basan en el levantamiento de datos dasométricos de inventarios forestales y el subsiguiente ajuste de modelos alométricos (), los cuales tienen la desventaja de poderse poner en marcha en grandes áreas por sus costos económicos, de tiempo y de mano de obra, además de la complejidad en la interpretación de los datos para áreas superiores a las inventariadas (). Por otra parte, los métodos basados en percepción remota son potencialmente adecuados para suministrar datos de buena exactitud sobre algunas variables biofísicas, a partir de los cuales se pueden obtener estimaciones confiables, eficientes y oportunas tanto de la estructura del bosque como de la biomasa aérea (). Ese potencial se debe, entre otros aspectos, a la capacidad de proporcionar observaciones completas y frecuentes en escalas locales y globales (), a la disponibilidad de datos gratuitos y al uso de una amplia gama de índices espectrales. Este es el caso de las imágenes obtenidas mediante los satélites Sentinel de la Administración Espacial Europea (ESA), los cuales tienen una resolución espacial, espectral y temporal apropiada para muchas aplicaciones de evaluación de los recursos naturales ().
Los métodos basados en percepción remota para estimar la biomasa aérea han sido aplicados con éxito, encontrándose que al combinar datos ópticos y datos radar de apertura sintética ( SAR) junto con algoritmos de aprendizaje automático, como Random Forest (RF) o Support Vector Machine (SVM), se obtiene una mejor correlación entre la biomasa aérea y las variables obtenidas de los datos espectrales (; ; ). Random Forest es un método estadístico no paramétrico considerado como uno de los mejores de clasificación y regresión debido la alta precisión en la estimación, su gran velocidad de cálculo, la robustez y la capacidad para predecir las variables importantes; por esto, ha sido ampliamente usado para estimar valores dendrométricos como AGB o volumen de un bosque (; ; ).
En plantaciones forestales en algunas regiones de Latinoamérica, y en Argentina y Perú respectivamente, se hallaron que existe una correlación significativa entre la biomasa total de Eucalyptus grandis y Pinus patula con algunos índices de vegetación extraídos de imágenes Landsat 5 y Sentinel-2A. En México encontró relaciones significativas entre las métricas de LiDAR y datos de campo para variables como altura, biomasa aérea, biomasa total y cobertura arbórea para especies de pino y encino. Sin embargo, en Colombia el uso de información de sensores remotos para estimar biomasa en plantaciones forestales no ha sido estudiado de manera sistemática y la mayoría de investigaciones que incorporan el uso de estos datos se han enfocado en bosques naturales (; ; ; ), dejando de lado estudios orientados al sector forestal comercial.
El vacío de conocimiento en el tema hace necesario realizar investigaciones que incorporen el uso de nuevas tecnologías satelitales dentro del sector forestal del país, las cuales podrán optimizar la estimación de la biomasa aérea mediante el diseño y prueba de métodos económicos y rápidos que presenten resultados confiables. El objetivo de esta investigación fue estimar la biomasa aérea de eucalipto (Eucalyptus grandis) y pino (Pinus spp.) en plantaciones forestales comerciales localizadas en el departamento del Cauca, a partir de imágenes ópticas Sentinel-2A, datos SAR Sentinel-1A y datos de inventarios forestales y mediante el uso del algoritmo de Random Forest (RF) para generar los mejores modelos de regresión asociados a AGB.
MATERIALES Y MÉTODOS
Área de estudio
El estudio se realizó en los municipios de Popayán, Cajibío, Sotará y Timbío, en el departamento del Cauca, Colombia, entre las coordenadas geográficas 2°38’44.92’’-2°16’10.49’’ norte y 76°56’0.68’’-76°32’40.57’’ oeste (figura 1). La temperatura promedio anual varía entre 18°C y 20°C, el régimen de precipitación es bimodal que anualmente oscila entre 1 900 mm y 2 400 mm, la altitud varía entre los 1674 y 3 116 m y las pendientes medias son del 24 % . Gran parte de la actividad económica de la zona está asociada a plantaciones forestales comerciales pertenecientes a la empresa papelera Smurfit Kappa Colombia (SKC), cuyo objetivo es la elaboración de pulpa () y fueron establecidas desde el año 2009 para las plantaciones de E. grandis, y desde el año 1999 para las plantaciones de Pinus spp. (P. oocarpa, P. patula y P. tecunumanii).
Datos de campo
Los datos fueron tomados en el mes de diciembre de 2015, mediante levantamiento de inventarios forestales en parcelas circulares de 250 m2 donde se registraron las variables altura (h) y diámetro a la altura del pecho (DAP), número de individuos y coordenadas del centro de las parcelas, censando todos los árboles vivos con DAP mayores a 5 cm . En total la zona tiene plantada un área de 39 km2 de E. grandis y 61 km2 de P. spp y se inventariaron 23 lotes, 13 de E. grandis y 10 de P. spp, los cuales tienen una extensión variable que oscila entre 0.15 km2 a 18 km2 para E. grandis y de 8 km2 a 27 km2 para P. spp (). Mediante modelos alométricos construidos por la empresa Smurfit Kappa Colombia , se estimó la biomasa aérea forestal para las dos especies en t.ha-1 usando DAP como variable explicativa. La tabla 1 presenta la distribución de la biomasa aérea en los lotes muestreados de acuerdo con la especie inventariada. Para el caso de datos del género Pinus (P. oocarpa, patula y tecunumanii) estos se agruparon en P. spp.
Datos ópticos Sentinel-2A y radar Sentinel-1ª
Las imágenes usadas en el presente estudio fueron obtenidas por la Agencia Espacial Europea (ESA) el 24 de diciembre de 2015 (Sentinel-2A_MSIL1C) y 11 de junio de 2016 (Sentinel-1A IW GRD). La imagen Sentinel-2A cuenta con 13 bandas espectrales con tamaño de píxel 10 metros (B2 azul, B3 verde, B4 roja y B8 infrarrojo cercano), 20 metros (B5, B6 y B7 de borde rojo, B8a de infrarrojo cercano, B11 y B12 de infrarrojo de onda corta) y 60 metros (B1 de aerosoles, B9 de vapor de agua y B10 de cirrus). Tres bandas corresponden a los rangos espectrales de “borde rojo”, las cuales son importantes para derivar información asociada al estado de la vegetación . La imagen Sentinel-1A de radar de apertura sintética (SAR) opera en la banda C y cuenta con polarización dual VV (Transmisión vertical y recepción vertical) / VH (Transmisión vertical y recepción horizontal), fue obtenida mediante una órbita ascendente con un ancho de franja de 250 km, tiene una resolución espacial de 5 × 5 m y un ángulo de incidencia entre 31.09 y 46.23 grados ().
Preprocesamiento de imágenes
Los datos Sentinel-2A fueron corregidos atmosféricamente transformando los datos de reflectancia en el tope de la atmósfera (TOA) (L1C) a reflectancia a nivel de la superficie (BOA) (L2A), a través de algoritmos aplicados en el módulo Sen2Cor del programa SNAP (). Dado que los datos de AGB en campo fueron capturados en un área circular de 250 m2, se re-muestrearon las bandas 5, 6, 7 y 8a, a un tamaño de píxel de 10 m empleando el método de convolución cúbica.
La imagen SAR fue calibrada radiométricamente usando la herramienta Terrain calibrate del programa SNAP, para representar los valores del píxel original a valores reales de retrodispersión de la superficie reflectante. Posteriormente, se aplicó un filtro Gamma MAP de 3 x 3 píxeles para disminuir el efecto moteado o “speckle” preservando los detalles más finos de la imagen . La imagen filtrada fue georrectificada usando el modelo de elevación digital (DEM) de Alos Palsar correspondiente a octubre de 2010 cuya resolución espacial es de 12.5 metros. Como resultado de este preprocesamiento se obtuvieron dos coeficientes de retrodispersión: σ°VV y σ°VH en decibeles (dB).
Índices espectrales a partir de datos Sentinel-2 y texturas a partir de datos Sentinel-1
Los índices espectrales fueron seleccionados de acuerdo con su poder predictivo en la estimación de AGB. El índice de vegetación de diferencia normalizada (NDVI) relaciona las bandas espectrales correspondientes al espectro rojo e infrarrojo y ha demostrado ser una variable de predicción efectiva para modelar la biomasa aérea en estudios previos (), aun cuando en algunos casos este índice se satura en coberturas densas (). El índice de vegetación de diferencia normalizada verde (GNDVI) es más sensible a la variación en el contenido de clorofila que el NDVI (). El índice de vegetación ajustado al suelo (Savi) presenta un adecuado resultado cuando la cobertura vegetal es densa, ya que tienen rangos dinámicos amplios y menor susceptibilidad a las perturbaciones (). Los índices se calcularon utilizando las expresiones presentadas en la tabla 2.
Los coeficientes de retrodispersión se utilizaron como entrada para el cálculo de las variables texturales para cada polarización (σ°VH y 𝜎°VV), usando la matriz de coocurrencia de niveles de gris (GLCM) (), con una ventana de 7 por 7 píxeles y para todas las direcciones (). En la tabla 3 se encuentran las ecuaciones de cada textura generada con la matriz GLCM.
Posterior a esto, los valores de los píxeles se extrajeron utilizando la herramienta de estadística focal usando un kernel de 3 x 3 píxeles, calculando la media y la desviación para disminuir el error de los datos obtenidos con el dispositivo GPS. La utilidad de este método fue recomendada por bajo el supuesto de que los valores de píxel representan un valor promedio de un área particular en el suelo.
Estimación de AGB con la técnica Random Forest (RF)
Mediante la librería Random Forest del programa estadístico R () se usó el algoritmo de regresión Random Forest (RF), basado en el aprendizaje automático no paramétrico a partir de una serie de árboles de decisión para identificar las variables predictoras importantes en la estimación de la AGB de cada especie (). Una variable predictora es importante en el modelo de regresión si al omitirla de la lista de variables predictoras aumenta el error “out of bag” (OOB) ). La importancia de las variables se cuantificó utilizando el porcentaje de aumento en el error cuadrático medio (IncMSE %), eliminando las variables menos importantes para luego construir el modelo final de cada especie (). En este estudio se utilizaron inicialmente 500 árboles (ntree) para cada especie y para tres subconjuntos de datos (Sentinel-2A, Sentinel-1AA y la combinación de ambos). Para el parámetro mtry, es decir, el número de variables a probar en cada nodo, se utilizó el número total de variables predictoras dividido sobre 3 (). Los parámetros de RF (mtry y ntree) se optimizaron con el fin de obtener el mejor poder predictivo en la estimación de AGB. La evaluación del rendimiento de los modelos se efectuó con el 70 % de los datos de campo para ajustes (n = 50 y n = 57 en E. grandis y P. spp, respectivamente) y el 30 % para validación (n = 22 y n = 24 E. grandis y P. spp, respectivamente), aplicando la técnica de validación cruzada 10 veces para probar de manera sólida el rendimiento del algoritmo, tomando como criterio el menor valor de EMC (). Finalmente, se mapeó la AGB forestal sobre cada lote usando el modelo AGB óptimo para cada especie con la interpolación Kriging (). Los mapas resultantes fueron generados en el programa RStudio (versión 1.1.442; ).
RESULTADOS
Importancia de las variables predictoras
Para E. grandis las variables derivadas de los datos Sentinel-2A más determinantes en la estimación de AGB fueron las bandas SWIR (B11 y B12), mientras que para Sentinel-1A fueron DisimilitudVH, CorrelaciónVV y σ°VV (figura 2). Al combinar los datos ópticos con SAR se encontró que las variables más importantes fueron las bandas SWIR y las texturas media, varianza y correlación de la polarización VV. Para P. spp GNDVI, bandas B2 y B3 y σ° VV, media y correlación de la polarización VV fueron las variables más importantes en la estimación de AGB (figuras 3a y 3b) y al combinar datos ópticos y SAR, CorrelaciónVV, σ°VH, σ°VV, GNDVI y B2 (figura 3c) están más asociadas a AGB mientras que ContrasteVV, B6 y DisimilitudVH no son influyentes.
Modelo y mapa de AGB forestal
De acuerdo con los R2 generados para cada grupo de datos, los resultados indican para E. grandis una relación negativa con la estimación de AGB para las imágenes de Sentinel-1A (tabla 4) a través del método RF y, por ende, se interpreta como nulo (). Sin embargo, al combinar los datos ópticos con SAR se mejoró la capacidad predictiva (R2 = 0.27) generándose el siguiente modelo con el método RF:
AGB E. grandis = 711.18-3 247.65*B12 - 283.2*B11- 428.36*(VV_MEA) + 218.92*(VV_VAR) + 120.84* (VV_COR) +5.784*σ° (VV) + 2.03*(VV_DIS)-346.81*B5 - 2 469.49*B3 - 4.36*σ° (VH) + 6 820.57*B4
*De acuerdo con el método RF las relaciones negativas se interpretan como un 0 %.
Para P. spp se encontró que los datos ópticos generaron EMC más bajos que los datos SAR (146.12 t.ha-1; tabla 4). Sin embargo, su combinación aumentó significativamente el coeficiente de determinación del modelo (R2 = 0.36). El modelo generado con el método RF para estimar AGB de P. spp fue:
Finalmente, la figura 4 presenta la distribución espacial de la AGB estimada obtenida a partir del método kriging usando la función Wave, en la que el valor medio estimado de AGB fue de 185.32 t.ha-1 y 499.06 t.ha-1 para E. grandis y P. spp respectivamente, superior al observado en campo.
DISCUSIÓN
Ranking de variables predictoras importantes con RF
En general, las variables derivadas de Sentinel mostraron respuestas diferentes para las especies forestales incluidas en el estudio; del mismo modo, la combinación de datos ópticos y SAR incrementó el potencial para estimar AGB, tal como lo han reportado los estudios de , para plantaciones en el trópico.
A diferencia de las plantaciones de P. spp, en las que todas las bandas e índices espectrales mostraron ser importantes para la estimación de AGB, posiblemente por la forma de la copa y las hojas, las bandas del infrarrojo de onda corta SWIR (B11 y B12) fueron las variables más importantes en E. grandis y esto se asocia a características de especies latifoliadas, relacionadas con contenido de humedad de las hojas, una mayor exposición del suelo de la plantación y una mayor eficiencia en la captura de energía por parte del dosel, lo que incide en su respuesta frente a la AGB
Si bien la importancia de las bandas de borde rojo en la estimación de AGB se ha abordado en investigaciones anteriores (), en este estudio se encontró que para el caso del índice GNDVI construido con la tercera banda de borde rojo (B7) la inclusión de este rango del espectro es importante para la estimación de AGB en E. grandis, dada su efectividad para monitorear la información del estado de la vegetación (; ; ). En P. spp GNDVI mostró ser la variable predictora más importante contradiciendo lo hallado por , pese a que en nuestro estudio la cantidad de AGB observado fue alta (489, 49 t.ha-1) y no se evidenció un efecto de saturación, posiblemente por el reemplazo de la banda roja por la banda verde (GNDVI en lugar de NDVI) (). Esto puede ser similar a lo presentado en E. grandis, donde NDVI no fue una variable importante a diferencia de lo hallado por en Brasil y que GNDVI puede ser más sensible a la concentración de clorofila-a que el NDVI (Askar, Nuthammachot, Phairuang, Wicaksono y Sayektiningsih, 2018).
La variable retrodispersión σ°VV (VV) derivada de los datos SAR mostró mayor importancia al momento de estimar biomasa aérea en plantaciones E. grandis y P. spp, comparado con la retrodispersión VH. Esto se contradice con el trabajo de , en el cual encontraron que la polarización VH obtuvo las correlaciones más fuertes (r > 0.917), en relación con las dos polarizaciones simples (HH y VV), y con , donde usaron datos Alos-Palsar para estimar AGB en T. grandis (r2VH= 0.46 y r2VV= 0.44). No obstante, coincide con lo hallado para S. robusta (r2VH= 0.52y r2VV= 0.67) y puede deberse a los valores que surgen de la dispersión del volumen, relacionada directamente con la AGB, y a la influencia de la humedad de la vegetación (). De igual manera, la retrodispersión en VV mostró mayor importancia en P. spp, debido a su relación con la irregularidad del dosel que se ha demostrado que es positivamente correlacionado con la biomasa de bosques tropicales más allá del punto de saturación de la relación de retrodispersión-biomasa (; ).
Por otro lado, las variables texturales derivadas de SAR mostraron bajas relaciones con los datos de AGB tomados en campo para ambas especies. Sin embargo, las texturas CorrelaciónVH y CorrelaciónVV para E. grandis y CorrelaciónVH y EntropíaVV para P. spp, derivadas de la matriz GLCM arrojaron mejores resultados que los coeficientes de retrodispersión VH y VV. Estos resultados son similares a los encontrados por , en el cual las medidas de textura media, correlación o ambas estuvieron involucradas en casi todos los modelos de AGB en bosques de coníferas, lo que implica que estas texturas tuvieron contribuciones significativas para mejorar las predicciones de AGB. encontraron que la adición de texturas derivadas de datos Alos-Palsar a los datos multiespectrales Landsat TM mejoró la estimación de AGB en bosques de Fagus orientalis. La razón por la cual algunas características texturales se correlacionaron mejor con la AGB que la retrodispersión SAR, se puede deber a la homogeneidad de los tipos de vegetación estudiados dado que los parámetros de textura son más aplicables en condiciones de alta variación local ).
Estimación de AGB de E. grandis y P. spp usando el método Random Forest
Los R2 encontrados en los modelos generados con RF son bajos comparados con otros estudios en plantaciones de coníferas y latifoliadas usando datos Landsat-5 y GLAS a través del método RF con R2 entre 0.73 y 0.96 (; ), y con el uso de datos Alos2-Palsar-2 y Sentinel-2A en un bosque tropical con R2 cercano a 0.70 ). Nuestros resultados pudieron verse afectados por la relación entre el tamaño del píxel de los datos Sentinel y el tamaño de la parcela usada para estimar la AGB en campo, tal como lo reporta y por la distribución espacial de las muestras dentro del área de estudio.
El EMC obtenido para E. grandis fue menor que el obtenido para P. spp; esto puede deberse a que se usaron datos de tres especies del género Pinus, las cuales logran tener respuestas espectrales diferentes asociadas a la variación en las tasas de crecimiento de cada especie y, por consiguiente, en la predicción de su biomasa, lo que aumentaría la variabilidad en los resultados.
Tanto para E. grandis como para P. spp el mejor modelo para estimar AGB forestal se obtuvo combinando los datos derivados de Sentinel-2A y los derivados de Sentinel-S1A. Este resultado coincide con y con los estudios de , y , quienes sugieren que el uso sinérgico de múltiples sensores ópticos y SAR tiene un mejor potencial para estimar la AGB.
También se demostró que al usar algoritmos de aprendizaje de máquina como Random Forest, en combinación con índices espectrales y datos SAR como variables predictoras, se generan mejores resultados para la estimación de AGB (; ; ; ).
CONCLUSIONES
Los resultados indicaron que la combinación de datos Sentinel-2A y Sentinel-1A pueden ser usados para generar mejores estimaciones de AGB forestal para las especies E. grandis y P. spp, usando el método Random Forest y apoyando el mapeo y monitoreo de plantaciones forestales a un bajo costo y de manera más rápida que las mediciones tradicionales en campo. La retrodispersión con polarización VV mostró ser una variable importante para la estimación de AGB en el área de estudio, así como el GNDVI y lo anterior está relacionado con la rugosidad de la vegetación y el contenido de humedad. Sin embargo, los resultados indican bajos valores de R2 en los modelos, por lo que se recomienda considerar un mayor tamaño de muestra y una mayor área de las parcelas en campo que incluyan más de un píxel para reducir el nivel de incertidumbre respecto a la asociación espacial entre los datos de las imágenes y los datos levantados en campo. De igual manera, para Pinus spp. es necesario discriminar las especies para la estimación de la biomasa usando el método Random Forest con el fin de reducir los EMC y obtener predicciones más confiables, aumentando el número de la muestra.
Finalmente, uno de los beneficios clave que ofrece el uso de imágenes de Sentinel es que tanto los datos de radar SAR como las imágenes multiespectrales están disponibles de forma gratuita que pueden procesarse en el programa gratuito SNAP. Por último, en países donde los recursos para la adquisición de imágenes y programas especializados son limitados el potencial de aplicación de esta metodología es alta.