Introducción

El inventario del volumen maderable o la biomasa contenidos en los bosques es considerado un proceso básico en su gestión, permite medir y monitorear a escala local los cambios que suceden en este tipo de ecosistemas (^{Urbazaev et al., 2018}). Su realización implica la medición de árboles ubicados en unidades de muestreo de dimensión fija, comúnmente distribuidas en una red sistemática sobre un área determinada (^{Grafström y Ringvall, 2013}). En las últimas décadas se ha impulsado el uso combinado de datos medidos en campo y datos colectados mediante sensores remotos para la realización de este tipo de inventarios (^{Grafström y Ringvall, 2013}; ^{Crespo-Peremarch et al., 2018}).

La tecnología LIDAR (Light Detection and Ranging) aerotransportada ha demostrado utilidad para estimar variables biométricas del bosque con precisión aceptable a partir del registro de retornos (rebotes) de un haz de luz previamente emitido por un equipo especializado (^{Ortiz-Reyes et al., 2015}; ^{Marchi et al., 2018}). Desafortunadamente, la adquisición de datos de LIDAR es costosa en comparación con datos provenientes de sensores ópticos o de radar. El costo está influenciado por el tiempo de vuelo y por la densidad de retornos por metro cuadrado que se desea registrar, los cuales a su vez varían con la velocidad, la altitud y el ancho de barrido del vuelo (^{Hansen Hofstad et al., 2015}).

En años recientes se han estudiado estrategias para abaratar el costo de usar datos LIDAR en la estimación de inventarios forestales; particularmente la reducción de la densidad de retornos por metro cuadrado y el uso de transectos en lugar de datos de cobertura completa (^{Almeida et al., 2019}; ^{Ortiz-Reyes et al., 2019}). Evaluaciones del efecto de la densidad de retornos en la precisión de la estimación de variables como la altura, diámetro cuadrático, área basal, volumen y biomasa, indican que las densidades óptimas varían de 0.5 m-2 a 15 m-2 en plantaciones, bosques tropicales (^{Hansen Hofstad et al., 2015}; ^{Silva et al., 2017}; Almeida et al., 2019) y en bosques templados puros o mixtos (^{Treitz et al., 2012}; ^{Singh et al., 2015}). Respecto al uso de datos LIDAR para solo una fracción del área a inventariar (a manera de una muestra), y a partir de esta expandir las estimaciones mediante modelos matemáticos alimentados con datos de otras imágenes de costo menor, se han publicado resultados prometedores (^{Ortiz-Reyes et al., 2019}). Aunque el método tiene limitaciones debido a la trayectoria lineal de la plataforma que porta el instrumento (avión) (^{Wulder et al., 2012}), éste es actualmente un procedimiento común en la selección de las muestras de campo, además de optimista para la estimación de variables de interés en el inventario (^{Mahoney et al., 2018}; ^{Ver et al., 2018}). De acuerdo con ^{Chen et al. (2012)} y ^{Wulder et al. (2012)} la adquisición de una muestra de datos LIDAR en transectos, en combinación con un marco de muestreo apropiado, genera estadísticos precisos y de bajo costo para áreas extensas.

Los transectos LIDAR pueden incorporarse a un diseño de muestreo en dos etapas en el que las unidades de muestreo primarias (UMP) son los transectos y las secundarias (UMS) son parcelas medidas en campo (^{Gregoire et al., 2011}; ^{Wulder et al., 2012}; ^{Strîmbu et al., 2017}). Mediante un estimador complejo en dos etapas, se puede estimar la media general de variables de interés. Para este propósito regularmente se han utilizado métodos paramétricos como la regresión lineal (^{Strîmbu et al., 2017}) y no lineal (^{Saarela et al., 2015}), además de modelos semiparamétricos o no paramétricos (^{McRoberts et al., 2011}; ^{Kangas et al., 2016}).

Aunque el uso de estimadores complejos representa una solución para la estimación del inventario de variables de interés forestal a través de transectos LIDAR, paralelamente se ha propuesto la extracción de métricas derivadas de sensores ópticos para combinarlos con datos LIDAR y con ello estimar parámetros forestales a escala regional (^{Chen et al., 2012}) y nacional. En México se han realizado estudios a escala nacional, combinando imágenes satelitales con datos LIDAR colectados en transectos para la estimación de biomasa aérea (^{Nelson et al., 2017}; ^{Urbazaev et al., 2018}). A escala regional este enfoque se ha explorado para la estimación de biomasa aérea, carbono, riqueza y diversidad de especies en la península de Yucatán (^{Hernández-Stefanoni et al., 2014}; ^{Ochoa-Franco et al., 2019}; ^{George-Chacon et al., 2019}; ^{Ortiz-Reyes et al., 2019}). Si bien esto demuestra el interés nacional por fomentar el uso de estas herramientas de teledetección y enfoques de muestreo para la estimación de variables forestales, su uso es todavía limitado, particularmente para la estimación de inventarios operativos.

Este estudio tuvo como propósito analizar de manera conjunta el efecto del diseño de muestreo LIDAR y la densidad de retornos por metro cuadrado sobre la precisión del inventario de área basal (AB), volumen de madera (V) y biomasa (B) en un bosque templado bajo manejo ubicado en Zacualtipán, Hidalgo.

Objetivos

Los objetivos de este estudio fueron analizar el efecto de diferentes diseños de muestreo y densidades de retornos (puntos) por unidad de área (m2) sobre la estimación del inventario de AB, V y B en áreas con datos LIDAR y combinar los resultados del mejor diseño encontrado para estimar y mapear el inventario total de AB, V y B en un predio forestal mediante estimadores de muestreo complejos asistidos por modelos aditivos generalizados (GAM).

Materiales y métodos

Área de estudio

El estudio se realizó en el área forestal que conforma el Sitio de Monitoreo Intensivo de Carbono (SMIC) en Atopixco, localizado en el municipio de Zacualtipán, Hidalgo, entre las coordenadas 20°40’17” y 20°34’51” latitud norte, 98°40’07” y 98°34’22” longitud oeste (Fig. 1). Su fisiografía comprende parte de la Sierra Alta Hidalguense y del Eje Neovolcánico. La precipitación promedio anual es de 1601 mm, con una temperatura media anual de 17 °C y clima templado húmedo con lluvias todo el año (Instituto Nacional de Estadística y Geografía [Inegi], 2017). El SMIC tiene una superficie de 900 ha cubiertas mayoritariamente de bosque coetáneo. Lo conforman rodales dominados por Pinus patula con edad variable desde 0 años hasta 31 años y algunas áreas sin intervención silvícola en los últimos 80 años (naturales) (^{Torres-Vivar et al., 2017}).

a) transectos en dirección este-oeste; b) transectos en dirección norte-sur; y c) estratificado (estratos definidos con base en el año (periodo) de cosecha del rodal).

FIGURA 1 Mapa del Sitio de Monitoreo Intensivo de Carbono Atopixco y diseños de muestreo LIDAR. 

Datos del sensor LIDAR

Los datos LIDAR fueron obtenidos en mayo del 2013 con el sensor Riegl-VQ480 a una altitud de 397 m, una frecuencia de 200 kHz, un ángulo de barrido de ±15º y una intensidad de retornos de 5 m-2 a 10 m-2. Los retornos registrados permitieron generar una “nube de puntos” tridimensional, la cual se constituyó como la base de datos LIDAR (^{Ortiz-Reyes et al., 2015}). Además, se generó un modelo digital de elevación (MDE) de 1 m de resolución para normalizar las alturas de los árboles sobre el suelo. La Figura 2 muestra los pasos metodológicos generales implementados en el presente estudio. En los párrafos siguientes se detalla el procedimiento.

FIGURA 2 Flujo de trabajo utilizado para la modelización espacial de variables de inventario forestal mediante diferentes densidades de retornos de datos LIDAR y espectrales. 

Resultados y discusión

Los valores promedio de AB, V y B (ha-1) más altos, estimados con datos de campo, se obtuvieron en el diseño de muestreo dirección Este-Oeste seguido del diseño estratificado (Tabla 1). Del total de métricas LIDAR calculadas, la “altura media del dosel” (Elev_mean) y la “proporción de todos los retornos por encima de la media entre el total de primeros retornos” (Allretabmean _Tot1stret100) mostraron una correlación alta con el AB, V y B; alcanzando coeficientes de correlación de hasta 0.92 (Tabla 1). Estas dos métricas se utilizaron para predecir el AB, V y B en los transectos con datos LIDAR adicionando como variables las coordenadas de longitud (x) y latitud (y), las cuales contribuyeron a los modelos GAM con un efecto espacial de suavizado bivariado.

La evaluación de los modelos GAM mediante validación cruzada (LOOCV) demostró que la reducción de la densidad de retornos a escala de sitio afectó la bondad de ajuste calculada, aumentando los valores de RMSE y disminuyendo el R2 en los tres diseños de muestreo implementados en densidades menores o iguales a 1 retorno por metro cuadrado. Además, los promedios de R2 calculados para los modelos de V, B y AB en las diferentes densidades de retornos oscilaron entre 0.86 y 0.90, 0.82 y 0.92, 0.76 y 0.83, respectivamente. El diseño estratificado (DE) presentó los valores promedio más bajos de RMSE (volumen = 37.01 m3 ha-1, biomasa = 18.26 Mg ha-1 y área basal = 4.97 m2 ha-1), seguido del diseño por transectos este-oeste (DH). Estos resultados demuestran que, si se utiliza el diseño estratificado, se captura la mayor variabilidad del atributo de interés. Por otro lado, de acuerdo con ^{Wulder et al. (2012)}, los transectos LIDAR obtenidos de forma sistemática, operativamente resultan menos costosos en áreas grandes, sin embargo, se debe identificar el gradiente que logre capturar la mayor variación de la masa forestal dentro de cada línea de vuelo. Lo anterior puede explicar el hecho de que la variabilidad en los valores de RMSE y R2 se mantuvo constantes y simétricos en estos diseños aun en densidades bajas de retornos por metro cuadrado (Fig. 3). Los resultados obtenidos concuerdan con lo señalado por ^{Hansen Hofstad et al. (2015)}, quienes indicaron que las métricas derivadas de datos LIDAR aerotransportado son confiables cuando se construyen a partir de densidades mayores a 0.5 retornos por metro cuadrado. En el mismo sentido, ^{Treitz et al. (2012)} utilizaron una densidad de 0.5 retornos por metro cuadrado para el ajuste de modelos lineales múltiples y demostraron que fue suficiente para la estimación del inventario forestal a escala de sitio en Ontario, Canadá. Conviene mencionar que el ajuste de los modelos puede aumentar si se normaliza la nube de puntos (retornos) respecto al suelo a partir de un MDE de buena calidad; este proceso fue evaluado por ^{Silva et al. (2017)}, quienes encontraron que las estimaciones de biomasa a escala de paisaje se subestimaron fuertemente a densidades por abajo de 0.8 retornos por metro cuadrado, principalmente en áreas con pendientes pronunciadas donde un MDE de baja calidad es más sensible en la normalización de las alturas. Al respecto, ^{Hansen Hofstad et al. (2015)} observaron que a partir de una densidad mayor a 1 retorno por metro cuadrado los efectos sobre la generación del MDE no fueron significativos. Este estudio demostró que una densidad promedio de 1 retorno por metro cuadrado puede considerarse como valor óptimo para escaneos LIDAR posteriores en el SMIC.

R2 y RMSE calculados mediante la validación cruzada repetida para; V (a-a1); B (b- b1) y AB (c-c1); los puntos negros representan puntos aberrantes.

FIGURA 3 Efecto de la densidad de retornos LIDAR en la modelización de área basal (AB), volumen (V) y biomasa (B) para los sitios muestreados bajo los diferentes diseños (n=96); DH (este-oeste), DV (norte-sur); DE (estratificado). 

Las estimaciones de AB y V para los transectos del diseño DH fueron mapeados a nivel pixel con fines comparativos, mediante la aplicación de modelos GAM ajustados utilizando tres densidades de retornos: 0.8 m-2, 1 m-2 y 10 m-2 (Fig. 4). Se ejemplifica este diseño por considerarlo una manera de uso operativo actual de un sensor LIDAR, además de que presentó bondades de ajuste adecuadas. Las estimaciones a partir de 0.8 m-2 y 1 m-2 retornos (densidades bajas) se contrastaron con las estimaciones a partir de 10 m-2 (densidad alta) mediante la prueba estadística no paramétrica Mann-Whitney-Wilcoxon (Fig. 4), encontrándose que las estimaciones de AB fueron diferentes (P < 0.5), sin embargo, las estimaciones de V no fueron significativamente diferentes en términos de distribución y media (P > 0.5). Además, se encontró que las estimaciones derivadas de la densidad de 1 retorno por metro cuadrado fueron más conservadoras en comparación con las de 10 retornos por metro cuadrado.

FIGURA 4 Estimacion de área basal y volumen en transectos este -oeste a partir de variables LIDAR calculados a densidades de retornos diferentes (izquierda). a y d) 0.8 pulsos por metro cuadrado; b y e) 1 pulso por metro cuadrado; c y f) 10 pulsos por metro cuadrado. Prueba Mann-Whitney-Wilcoxon (0.6-10 y 1-10) (derecha). 

Las estimaciones de AB, V y B generadas para el área de los transectos a partir de un modelo GAM y una densidad de 1 retorno por metro cuadrado (considerado el más adecuado), sirvieron como variables de respuesta para ajustar otros modelos GAM y RF para estimar y mapear las mismas variables para todo el SMIC. Se utilizaron como predictores los valores de las bandas espectrales, el NDVI y métricas de textura derivadas del sensor RapidEye que registraron un coeficiente de correlación de Pearson alto con las variables de interés (AB, V y B): texturas construidas a partir de la banda roja (B3) (media y varianza), azul (media) (B1_mean), verde (media) (B2_mean) y el NDVI (media y varianza) (ventana móvil de 15 × 15) (Tabla 2). El modelo GAM presentó menor valor del porcentaje de varianza explicada (R2) y mayor RMSE en la estimación de AB, V y B, comparado con RF (Tabla 3). Sin embargo, todos los parámetros de suavizado que corresponden a las texturas del RGB y NDVI fueron altamente significativos (P < 0.0001).

El modelo RF resultante incluyó cinco bandas del sensor Rapideye y 42 texturas, además de las coordenadas longitud (x) y latitud (y). Este algoritmo logró explicar las variables dasométricas de interés por arriba de 70%, teniendo mejor rendimiento en la estimación de la B (79%, RMSE = 23.67 Mg ha-1) seguido del AB (73%, RMSE = 5.69 m2 ha-1) y V (71%, RMSE = 46.95 m3 ha-1).

RF mostró mejores bondades de ajuste al utilizar texturas derivadas del sensor RapidEye, lo que coincide con lo observado en otros estudios (^{Matasci et al., 2018}; ^{Zhu et al., 2020}). Una ventaja de RF es que implementa una validación cruzada a partir de conjuntos de datos de entrenamiento y validación para evitar la inflación de la varianza causada por el uso de variables redundantes en la explicación de la variable de interés (^{Fassnacht et al., 2014}). Además, se ha observado una mejoría en las estimaciones con el uso de texturas usando el método GLCM, con un tamaño de ventana de 15 × 15 (^{Zhu et al., 2020}) o mayor dependiendo la variabilidad del dosel. A diferencia de las métricas comunes, las métricas de textura espectral pueden reducir los impactos de sombra y extraer más eficazmente la información de bosques maduros, con ello pueden aumentar las correlaciones entre la textura y la masa forestal madura (^{Lu y Batistella, 2005}), como en el caso del presente estudio.

Los mapas de AB y el V correspondientes al área boscosa del SMIC (815.03 ha) muestran que ambos modelos (GAM y RF) estiman con tendencia similar (Fig. 5). Los datos indican que el modelo GAM estimó valores totales de AB y V de menor magnitud que RF (Tabla 4), pero valores absolutos más altos en ambas variables (Fig. 6). Conforme a lo esperado, los valores más bajos de AB y V corresponden a rodales jóvenes, cosechados en 2003-2013 y los valores máximos corresponden a los rodales maduros y a áreas sin intervención silvícola (Tabla 5).

FIGURA 5 Variabilidad espacial del área basal y volumen de madera estimada a partir de modelos aditivos generalizados y el algoritmo random forest utilizando métricas derivadas del sensor RapidEye. 

FIGURA 6 Diagrama de caja de las estimaciones de volumen de madera generadas a partir de modelos aditivos generalizados y el algoritmo random forest ajustados utilizando variables del sensor RapidEye. 

La muestra de datos LIDAR (transectos) que representó 30% del total de la superficie del área de estudio, fue suficiente para estimar con precisión el AB, V y B, y para expandir las estimaciones a toda el área mediante la incorporación de datos espectrales RapidEye. Lo anterior sugiere que la estimación del tamaño óptimo de la muestra requiere de estudios más detallados al respecto. ^{Chen y Hay (2011)} observaron que una muestra de 8.8% (un solo transecto) de la superficie total aumentó el rendimiento (Volumen: R2= 0.64 y RMSE= 156.9 m3 ha-1) del algoritmo de aprendizaje utilizado [regresión de soporte de vectores (SVR Support Vector Regression)] para predecir variables de interés. Estos autores diseñaron cuatro muestras de transectos LIDAR (8.8%, 17.6%, 26.4% y 35.2% del área total) y las combinaron con imágenes QuickBird para estimar la altura del dosel, biomasa y volumen en Vancouver Island, Canadá. En otro estudio, ^{Chen et al. (2012)} evaluaron la proporción y dirección (Norte-Sur y Este-Oeste) de transectos LIDAR para estimar altura, biomasa y volumen en Quebec, Canadá y encontraron que un solo transecto que representó 7.6% del área total presentó mejores bondades de ajuste (volumen: R2 = 0.72; RMSE = 52.59 m3 ha-1) en el algoritmo SVR.

Conclusiones

Los modelos GAM ajustados bajo tres diseños de muestreo para inventario forestal asistido por datos de un sensor LIDAR aerotransportado permitió identificar al diseño estratificado (DE) como el mejor (RMSE: V = 37.01 m3 ha-1, B = 18.26 Mg ha-1 y AB= 4.97 m2 ha-1), aun en densidades de retornos bajas. Aunque este diseño de muestreo no se utilizó para la estimación del inventario total del SMIC, los resultados alientan la posibilidad de su aplicación en inventarios forestales operativos mediante la captura de datos LIDAR a través de sensores montados en drones.

Se encontró que una densidad de datos LIDAR de 1 retorno por metro cuadrado fue suficiente para obtener estimaciones precisas de AB, V y B en bosques templados bajo manejo coetáneo, lo que permitirá considerar en futuros planes de vuelo para colecta de datos LIDAR una mayor altura de vuelo, con la consecuente reducción de tiempo y costo de escaneo. Los modelos de regresión GAM y al algoritmo Random Forest demostraron propiedades importantes en la estimación de las variables de interés forestal tanto a escala de transectos como en áreas donde la información LIDAR no está presente, pero se cuenta con otro tipo de imágenes ópticas de resolución espacial fina.

Los estimadores de muestreo en dos fases y en dos etapas demostraron ser estimadores compatibles con las características de vuelo LIDAR mediante transectos, permitiendo alcanzar precisiones de muestreo menores o iguales a 6%.