Censo forestal de la parcela

El censo de la parcela de estudio fue realizado entre los años 2002 y 2010. La identificación de árboles siguió una metodología estandarizada utilizada en una red mundial de grandes parcelas forestales (^{Condit, 1998}). Durante estos trabajos de campo fueron etiquetadas, identificadas y localizadas todas las plantas leñosas ≥ 1 cm de diámetro de tallo. ^{Valencia et al. (2004)} identificaron todos los árboles con diámetro ≥ 1cm en 25 hectáreas de la parcela de estudio. Hasta el año 2010 se realizaron tres censos de individuos de diámetro ≥ 1 cm en la mitad occidental de la parcela (25 ha) mientras que en el resto de la parcela se realizó un censo de individuos de diámetro <10 cm y dos censos de individuos de diámetro ≥ 10 cm. En el censo de 25 hectáreas se identificaron 1 104 morfoespecies entre 152 353 individuos mientras que en las 50 hectáreas se censaron 341 339 individuos.

Se dispuso de mediciones directas del diámetro del tronco a la altura de pecho de todos los individuos de la parcela de 50 hectáreas, y datos de densidad de madera (gravedad específica) tomados de literatura especializada y de mediciones directas en el entorno de la parcela o áreas del PNY (^{Valencia et al., 2009}). En estos censos no se realizaron mediciones de alturas de los árboles debido al dosel denso y cerrado de la zona.

La elevación media del dosel en la parcela es de 18-25 metros. También se encuentran árboles emergentes de hasta 40 m de altura y en escasas ocasiones árboles de 50 m de altura. El mayor número de especies de árboles tenían un dosel cuyo tamaño medio tenía una altura de 10 a 20 m, con presencia de árboles de menor tamaño en el sotobosque con alturas máximas de 5-10 m. Las especies de árboles pequeños dominantes en el sotobosque eran de los géneros Matisia y Rinorea, mientras que las especies de dosel importantes fueron las Iriartea deltoidea y Eschweilera coriácea (^{Valencia et al., 2004}b).

Modelo general para estimación de biomasa con variables medidas en campo

Las estimaciones de biomasa de los bosques tropicales están limitadas por el conocimiento de la alometría de árboles tropicales. Su extrema diversidad de especies de árboles en general impide que se desarrollen alometrías de especies específicas, por lo que generalmente se aplican relaciones generales. La precisión para el cálculo del stock de carbono depende de la disponibilidad de modelos alométricos fiables que permitan inferir la biomasa aérea de los árboles en bosque tropical a partir de datos de inventario forestal (^{Chave et al., 2005}, 2014). Por la imposibilidad de desarrollar relaciones alométricas para un área específica de interés, es común utilizar ecuaciones alométricas ya existentes. Debido a la gran cantidad de estudios que documentan tales relaciones (^{Feldpausch et al., 2012}), es importante identificar correctamente las ecuaciones más representativas.

Para la estimación alométrica de la biomasa aérea, se aplicó el enfoque general de ^{Chave et al. (2014)}. Este modelo (Ecuación 1) se usa en aquellos casos en los que no se cuenta con mediciones de la altura total de los árboles, e incorpora variables como diámetro del tronco, densidad de la madera, y la variable E. Este modelo mejora las predicciones de biomasa propuestos en Chave et al. (2005).

Donde D representa el diámetro del tronco a la altura de pecho en cm, y ρ es la densidad de madera en g cm^-3 de cada árbol. La relación genérica entre el valor de biomasa total sobre el suelo y las características forestales del área de estudio, depende linealmente de una variable de estrés bioclimático (E), que a su vez depende de índices de variabilidad de temperatura, del índice de precipitaciones, y de la intensidad de la sequía.

Con los datos recogidos en el censo y el valor del parámetro E en el Parque Nacional Yasuní (-0.0228121), se calculó la cantidad de biomasa sobre el suelo (en unidades de Mg C ha-1), de todos los individuos existentes en la parcela. La sumatoria de las estimaciones en cada árbol, nos dio el valor total de biomasa sobre el suelo para toda la parcela.

Determinación del tamaño de la parcela para el análisis de estimación de carbono

Es necesario establecer parcelas de gran tamaño para realizar estudios de dinámica del bosque porque abarcan poblaciones importantes de la mayoría de las especies de árboles en la comunidad, especialmente en bosques tropicales (^{Condit et al., 2000}). Hay una tendencia a que los errores en las estimaciones de biomasa disminuyan con el aumento de tamaño de la parcela, ya que parcelas grandes reducen la probabilidad de los efectos de borde, que se producen cuando el dosel de los árboles se encuentra a lo largo de los límites de la parcela. También hay una disminución de errores en general, porque las parcelas grandes capturan una cantidad adecuada de la variabilidad estructural en el campo (^{Chave y Condit, 2004}). Una compilación de investigaciones publicadas sobre estimaciones de biomasa sobre el suelo usando LiDAR (^{Zolkos et al., 2013}) concluye que la incertidumbre se acerca al 10% cuando el tamaño de la parcela se aproxima a 1 ha. Las estimaciones de biomasa sobre el suelo a escalas espaciales por debajo de 1 ha están dominadas por la incertidumbre, por los efectos de corte de parcela y por las diferencias en las características de los sensores (^{Mascaro et al., 2011}). Se admite que escalas espaciales, típicamente del orden de 1 ha, dan resultados más precisos y son más adecuadas en estudios para escalas más finas (^{Meyer et al., 2013}).

Tomando en cuenta los hallazgos de las investigaciones mencionadas, el diseño de muestreo para el análisis fue un cuadriculado regular de subparcelas de 1 ha de superficie cada una de ellas (Figura 2), que permitió capturar el rango de variación de la estructura forestal en la zona de estudio. Las coordenadas X, Y de la esquina SW obtenidas en el levantamiento planimétrico realizado en la parcela de estudio (Figura 3) se usaron como punto de partida; de este modo las subparcelas quedaron georreferenciadas y numeradas de 1 a 50 comenzando en la esquina inferior izquierda.

Figura 2 Sistema de cuadriculado regular de subparcelas de 1 ha de superficie, utilizado para el cálculo de la biomasa de la zona de estudio e identificación de las parcelas para el ajuste (trama horizontal roja) y validación (trama vertical verde). 

Figura 3 Ubicación de los 4 puntos GNSS y las 4 esquinas de la parcela; Mapa base: C.T Rio Tivacuno, 1:50.000, IGM. 

Geoposicionamiento y levantamiento planimétrico de la parcela de estudio

Para determinar con precisión submétrica los vértices de la parcela de estudio, y relacionar con la máxima exactitud la información del censo con los datos obtenidos mediante el sistema LiDAR aerotransportado, se realizó un levantamiento sobre la parcela de estudio. El Sistema de Referencia Horizontal utilizado en este levantamiento fue el Sistema de Referencia Geocéntrico para las Américas (SIRGAS, Ecuador), compatible con el sistema de posicionamiento por satélite GPS.

Para el levantamiento planimétrico se implantaron 4 vértices geodésicos (NW-1, NW-2, NE-1, NE-2), que sirvieron de línea base (Figura 3), se utilizó la tecnología GNSS para determinar su ubicación. En la observación se utilizó el vértice geodésico GNSS “Y-NPF”, que forma parte de la Red Nacional GPS del Ecuador (RENAGE). Esta red se encuentra enlazada a SIRGAS (Sistema de Referencia Geocéntrico para las Américas) - ITRF94 (International Terrestrial Reference Frame), época de referencia 1995.4, que como se ha indicado es un marco de referencia compatible con el sistema de posicionamiento satelital GNSS.

La determinación de las coordenadas de los vértices geodésicos se realizó por el método de posicionamiento GPS estático diferencial de fase, utilizando receptores doble frecuencia marca TRIMBLE, modelo R8, que permiten la medida de líneas base hasta 200 kilómetros, con una precisión de ± 0.005m + 1 ppm. Para el procesamiento de los datos GPS se utilizó el software de procesamiento TRIMBLE BUSINESS CENTER (TBC) Versión 2.5. Para el procesamiento de cada vector GPS, se fijó la precisión horizontal de < 0.050 m + 1 ppm, precisión vertical < 0.100m + 1 ppm.

Posteriormente, se realizaron las observaciones necesarias para calcular las coordenadas de las 4 esquinas de la parcela (Tabla 1), mediante una estación total TRIMBLE modelo S3, con una precisión de 2”.

Toma de datos con sensor LiDAR aerotransportado

Los datos LiDAR fueron colectados en mayo de 2014, la plataforma aérea utilizada fue una avioneta CESSNA 172 modelo Skyhawk que permite cubrir grandes extensiones de terreno a bajas altitudes y velocidades bajas. El equipo LiDAR utilizado fue un Optech ALT Gemini multiretorno. Los datos y configuración del vuelo y equipo se detallan en la Tabla 2.

En la parcela de 50 hectáreas se registraron 9,603.411 pulsos en total, si bien el número de retornos fueron 4, el 89.17 % de los datos LiDAR registrados (8,563.795) pertenecieron al primer retorno. Este hecho puede ser explicado por la altura y densidad de la vegetación en la zona de estudio.

Procesamiento de los datos LiDAR

La cuantificación de la biomasa sobre el suelo (AGB) depende de las relaciones entre métricas LiDAR y estimaciones de la biomasa a nivel de parcela a partir de mediciones de campo y ecuaciones alométricas (^{Zolkos et al., 2013}). La estimación de la altura del dosel forestal utilizando datos LiDAR depende de una representación exacta de la superficie del suelo en los modelos digitales del terreno (MDT). Para los estudios forestales en particular, la tecnología LiDAR es capaz de caracterizar tanto terreno y la estructura de la vegetación de manera efectiva. Sin embargo, cualquier error en el MDT se propaga y afecta a la precisión de las métricas derivadas de la vegetación (^{Tinkham et al., 2012}) y modelos de altura de dosel.

Los datos LiDAR obtenidos se normalizaron a la superficie del suelo a fin de expresar los retornos en términos de alturas sobre el suelo y se recortaron en forma de malla de 100 m x 100 m usando el software LAStools. A continuación, utilizando el programa FUSION, se calculó un conjunto de datos estadísticos a partir de la nube de puntos LiDAR de cada una de las parcelas de 1 ha y para todo el grupo de parcelas de muestreo (50 has). Estos datos son métricas descriptivas de alturas obtenidas con LiDAR como la mediana, la moda, valores máximos y mínimos, el rango intercuartílico, los distintos percentiles y cuartiles entre otros.

Estas métricas fueron utilizadas para relacionar la nube de puntos con la información alométrica del inventario de parcela de campo (diámetro del tronco, área basal, densidad de madera) y para el cálculo de biomasa sobre el suelo calculado con fórmula propuesta por ^{Chave et al. (2014)}.

Metodología de calibración y ajuste del modelo general

Las mediciones con LiDAR producen una cobertura detallada de la altura y métricas tridimensionales del dosel que se pueden correlacionar con el área basal, dada la naturaleza volumétrica de las mediciones en 3-D (^{Drake et al., 2003}). La mayoría de los estudios se han centrado hasta el momento en bosques de coníferas y plantaciones forestales en zonas templadas. En esta investigación se analizó esta correlación en los bosques tropicales Amazónicos occidentales que contienen mucha diversidad florística (^{Bass et al., 2010}).

Modelo general de estimación de densidad de carbono sobre el suelo (ACD)

Al igual que Asner et al. (²⁰¹¹), consideramos que existen tres factores estructurales básicos que describen la cantidad de carbono almacenado en los bosques tropicales: área basimétrica o basal (AB), la altura de los árboles (H) y la densidad específica de la madera (WD). La información obtenida a través de la tecnología LiDAR aerotransportada, está relacionada con las principales variables de masa forestal mediante el ajuste de modelos estadísticos, pudiendo estudiarse variables tales como la altura dominante, el área basimétrica, y la biomasa.

Para eso es necesario tener dos fuentes de datos: 1) Mediciones dasométricas en parcelas de campo convencionales geoposicionadas con muy alta precisión (submétricas) para garantizar que todas las fuentes de datos disponibles para el estudio coincidan espacialmente; 2) nubes de puntos obtenidos mediante tecnología LiDAR aerotransportada de la zona de estudio, referenciándose de este modo la posición geométrica de la nube de puntos con la información dasométrica de las parcelas. Errores en la identificación y la baja precisión de las variables estimadas, son a menudo atribuidas a discrepancias entre el co-registro de las parcelas de referencia en tierra y la información LiDAR (^{Næsset, 2002}).

Alometría conjunta de parcela

El enfoque de alometría conjunta de parcela para calibración con LiDAR propuesta por Asner et al. (2012), sigue la teoría general de alometría de árbol (^{Chave et al., 2005}), y refleja las propiedades de la estructura forestal de la parcela completa en conjunto, en lugar de las propiedades de cada organismo en particular. Este enfoque radica en un modelo universal simplificado que demuestra que la biomasa seca de árboles y por lo tanto el contenido de carbono, que es ~ 48% de la biomasa seca en peso (^{Martin y Thomas, 2011}), puede estimarse a partir de mediciones de alturas obtenidas con datos LiDAR y otras variables medidas en campo (Asner et al., 2012).

El modelo universal utilizado en este estudio (Ecuación 2) ha sido calibrado y validado experimentalmente en cuatro regiones tropicales: Perú, Panamá, Madagascar y Hawai. El interés de esta investigación es continuar con la validación de este modelo, evaluando si es posible su implementación en otras regiones geográficas o queda relegado a la zona de estudio de la que surgió. El modelo universal se define de la siguiente:

Donde la densidad de carbono sobre el suelo (ACD) corresponde al valor de la biomasa sobre el suelo obtenido con la Ecuación 1 y multiplicado por el factor 0.48. MCH es la altura media al tope del dosel obtenida con LiDAR, AB es el área basal medida en campo, y WD _BA es la densidad de la madera promedio (promedio ponderado con respecto al área basal para cada parcela de 1 ha).

Selección de las subparcelas para el ajuste y validación del modelo general

Los resultados obtenidos por ^{Réjou-Méchain et al. (2014)} refuerzan la importancia de la topografía como un factor que debe ser tomado en cuenta en el diseño de los inventarios forestales y las calibraciones cuando utilizamos sensores remotos. El estudio sugiere que el muestreo en general, debe ser estratificada por hábitats topográficos (por ejemplo, crestas o cimas, valles y laderas), sobre todo si las estimaciones consideran un enfoque basado en el terreno. Por el contrario, cuando la finalidad del muestreo de campo es calibrar productos de detección remota, las zonas topográficamente complejas se deben evitar para reducir al mínimo los errores de muestreo asociados a la variabilidad espacial local. En los inventarios forestales y los sistemas de calibración, una importante fuente de error es la distribución desigual y no aleatoria de parcelas en escalas espaciales amplias, especialmente en zonas tropicales.

Para la selección de las subparcelas para el ajuste y validación del modelo en nuestro estudio, se consideró la distribución espacial y hábitats topográficos dentro de la parcela de 50 ha (valle, ladera y cima). Se utilizaron 32 subparcelas para el ajuste y 18 subparcelas para la validación del modelo (Figura 2). La subparcela No. 2, contiene un remanente de bosque secundario (0.48 ha) resultado de la tala de árboles para la construcción de un antiguo helipuerto abierto en épocas de prospección petrolera en la zona (alrededor del año 1987). Esta subparcela fue considerada en la validación del modelo. En el análisis exploratorio de datos se identificaron dos subparcelas atípicas (la subparcela 36 y la subparcela 50) que contienen los árboles más grandes de la zona. Ambas subparcelas fueron incluidas en la validación del modelo.

Para el ajuste del modelo, se consideraron sólo los árboles ≥10 cm de diámetro a la altura de pecho, ya que los arboles <10 cm generalmente contribuyen con menos del 5% de la biomasa aérea total (AGB) en los bosques tropicales maduros (^{Chave et al., 2003}), y es el valor de corte usual en las medidas de inventario forestal y estimaciones de biomasa usando LiDAR (Feldpausch et al., 2011). Para cada una de las subparcelas, se elaboró una tabla con datos del número de árboles, código de especie, diámetro del árbol, ubicación espacial (X, Y), área basal, valor de biomasa estimada (Mg_ha-1) a partir de la fórmula propuesta por Chave et al., (2014), densidad de madera por especie y la densidad de madera ponderado por el área basal (gr_cm-3) para cada uno de los individuos en las subparcelas.

Ajuste del modelo local a partir del modelo general de cálculo de densidad de carbono sobre el suelo (ACD) para la parcela de estudio

El modelo universal (Ecuación 2) se ajustó mediante regresión lineal múltiple por el método de mínimos cuadrados, transformando logarítmicamente las variables ACD, TCH, AB, y WD _BA en la forma (Ecuación 3):

Se utilizaron los datos LiDAR y la información dasométrica de 32 subparcelas. El ajuste y diagnosis del modelo se lo realizó con el software R Commander (^{Fox, 2005}).

Modelo final de cálculo de ACD

Luego del análisis y ajuste del modelo general propuesto, el resultado del modelo matemático local que mejor se ajusta es el siguiente:

Los estadísticos obtenidos para cuantificar la bondad de las estimaciones son: el error cuadrático medio (RMSE) igual 0.03182 y el coeficiente de determinación R²= 0.966 que determina que porcentaje de la varianza de la variable dependiente es explicado por el modelo de regresión.

A continuación, se revirtió el modelo final ya que nos interesa el parámetro ACD por ha y no su ln y se multiplicó por el factor de corrección (FC) para eliminar el sesgo producido (subestimación) por la transformación logarítmica en el modelo (^{Baskerville, 1972}). El factor de corrección está dado por CF= e MSE/ 2, donde MSE es el error cuadrático medio del modelo de regresión. Para este caso es igual a 1.00044. La ecuación quedó determinada de la siguiente manera: (Ecuación 6)

Validación del modelo local para la estimación de densidad de carbono sobre el suelo (ACD)

Un procedimiento riguroso de validación del modelo consiste en verificar los resultados con una muestra independiente de la utilizada para construirlo. Para este estudio se realizó la validación del modelo aplicándola a 18 subparcelas que fueron seleccionadas para tal fin. Al aplicar el modelo resultante (Ecuación 5) a las subparcelas de validación, se obtuvo el resultado que se muestra en la Figura 4.

Figura 4 Validación del modelo general de ACD en 18 subparcelas independientes. 

Los resultados de la validación directa ponen de manifiesto que las estimaciones de ACD usando el modelo general son precisas (RMSE= 5.81 Mg C ha-1), y presentan una alta consistencia con las medidas de campo (R2= 0.94), sin apenas desviaciones sistemáticas (sesgo próximo a 0), recogiendo la variabilidad del ecosistema en la zona de estudio.

Validación del modelo por posiciones topográficas

En la Figura 5 se presenta el desempeño del modelo ajustado para cada una de las posiciones topográficas. Los resultados muestran que el modelo es espacialmente sensible a las variaciones de ACD en las 3 posiciones topográficas. La selección de subparcelas para el ajuste del modelo considerando estas posiciones, incide en la precisión de las estimaciones de ACD en la zona de estudio. El sesgo en cada posición topográfica proporciona un conocimiento general de cómo las parcelas de campo representan el paisaje. La disminución del sesgo es independiente del número de parcelas de calibración; sólo depende del error de muestreo asociado a estas parcelas, que está determinado en gran medida por el tamaño de la parcela (^{Réjou-Méchain et al., 2014}).

Figura 5 Validación del modelo general de ACD por posiciones topográficas en: a) valle; b) ladera; c) cima incluyendo el remanente de bosque secundario; d) cima sin bosque secundario. 

La mayor diferencia entre las observaciones de campo y las estimadas por el modelo se presenta en la cima que contiene el remanente de bosque secundario (Figura 5c). Esta diferencia se debe a dos situaciones muy concretas: a) el mayor crecimiento de los árboles del bosque secundario en comparación con los árboles del bosque maduro en la parcela; y b) el tiempo transcurridos entre el censo de árboles y la toma de datos LiDAR. La Figura 5d muestra la validación del modelo sin la parcela 2, donde se puede comprobar la mejora en la precisión del modelo.

Geografía del carbono en la zona de estudio

Los mapas de biomasa y contenido de carbono sobre el suelo de los bosques tropicales son importantes para la conservación de los bosques, ya que constituyen un material básico para la elaboración de políticas de mitigación del cambio climático, y la comprensión del papel que desempeñan los bosques tropicales en el ciclo global del carbono.

Para obtener el mapa de alta resolución de las reservas de carbono en el área de estudio, se aplicó la Ecuación 5 a toda la superficie de cobertura LiDAR utilizando algebra de mapas. La zona total de cobertura LiDAR se dividió en unidades de análisis de 1 ha (malla regular de celdas cuadradas -fishnet), y se calculó el TCH (Figura 6) a partir del modelo de digital de superficie normalizado.

Figura 6 Mapa de altura media al tope del dosel (TCH) para la zona de estudio 

El mapa resultante de la densidad de carbono sobre el suelo (Figura 7), a pesar de cubrir una pequeña extensión de territorio, es una muestra representativa de la diversidad forestal del Parque Nacional Yasuní, y consecuentemente de una gran parte del bosque tropical amazónico del Ecuador.

Figura 7 Mapa de densidad de carbono sobre el suelo en la zona de estudio. 

Los valores más altos de ACD (168 Mg C ha-1) se localizan en los bosques ubicados en los sitios más elevados (cimas), mientras que las menores densidades de ACD se encuentran en los bosques en los valles (82 Mg C ha-1), y en zonas inundadas, o propensas a inundaciones (47 Mg C ha-1). Las zonas con almacenamiento bajo de carbono se deben probablemente a las condiciones hidrológicas como la anoxia, y posiblemente a cambios en la fertilidad del suelo asociados con el substrato geológico subyacente (^{Asner et al., 2010}). Los resultados obtenidos concuerdan con el estudio realizado por el Ministerio del Ambiente de Perú el año 2014 en la Amazonía Norte del Perú, donde encontraron que las mayores densidades de carbono se localizan al norte del río Napo, en el límite fronterizo con el Ecuador.