Source: http://mingaonline.uach.cl/scielo.php?script=sci_arttext&pid=S0717-92002017000100003&lng=en&nrm=iso&tlng=es
Timestamp: 2019-05-20 08:54:02+00:00

Document:
Bosque (Valdivia) - Estimation of timber stocks of Pinus radiata stands at the provincial scale using low resolution LiDAR data
Estimación de las existencias maderables de Pinus radiata a escala provincial utilizando datos LiDAR de baja resolución
Estimation of timber stocks of Pinus radiata stands at the provincial scale using low resolution LiDAR data
Saray Martín-García a, Ulises Diéguez-Aranda a, Juan Gabriel Álvarez González a, César Pérez-Cruzado a, Sandra Buján b, Eduardo González-Ferreiro a,c,d*
*Autor de correspondencia: a Universidade de Santiago de Compostela, Departamento de Enxeñería Agroforestal, Unidade de Xestión Forestal Sostible (UXFS), Escola Politécnica Superior, R/ Benigno Ledo, Campus Terra, 27002 Lugo, España, tel.: +34 982 823 222, edu.g.ferreiro@gmail.com
b Universidade de Santiago de Compostela, Departamento de Enxeñería Agroforestal, Laboratorio do Territorio (LaboraTe), Lugo, España.
c Oregon State University, Department of Forest Ecosystems and Society, Corvallis, Oregon, USA.
d USDA Forest Service - Pacific Northwest Research Station, Laboratory of Applications of Remote Sensing in Ecology (LARSE), Corvallis, Oregon, USA.
In this study, the stand volume (V, m3 ha-1) of Pinus radiata plantations of the province of Lugo (NW of Spain) was estimated using LiDAR-derived metrics from the low-density countrywide LiDAR dataset provided by the Spanish National Aerial Photography Program. The estimate was based on a oriented field inventory design to improve the fit of the model and on a inferential method based on the variability of the model parameter estimation and the variability of the LiDAR variable in the target population, which allowed to estimate the population mean and the standard error of V. A total of 25 permanent field inventory plots -that represent adequately the existing range of ages, stand densities and sites of the species in the province- were measured in order to fit the empirical relationship between field-measured V and LiDAR-derived metrics, using linear and (multiplicative) power function models. Regression analysis indicated a strong relationship between V and the 70 percentile of the LiDAR normalized height distribution (h70) using a power function model (R2 = 0.899; RMSE = 63.4 m3 ha-1). The selected model was applied to all timber and high-pole development stands of P. radiata with coverage of this species superior to the 60 % according to the Spanish National Forest Map. Finally, a spatial-explicit map that describes V was produced. In 2009, the estimated mean population for V was 166.3 m3 ha-1 with a relative sampling error of 13.9 %. This methodology will allow generating thematic maps every six years (temporal resolution scheduled for the PNOA LiDAR flights).
Key words: Airborne Laser Scanning (ALS), Forest Inventory (FI), Remote Sensing (RS), estimators, regression models.
Se estimó el volumen de madera de rodal (V, m3 ha-1) de Pinus radiata en la provincia de Lugo (NO España) a partir de métricas de datos LiDAR de baja resolución del PNOA (Plan Nacional de Ortofotografía de España). La estimación se basó en un diseño de muestreo orientado para mejorar el ajuste del modelo y en un método inferencial basado en la variabilidad de la estimación de los parámetros del modelo y en la variabilidad de la variable LiDAR en la población objetivo, lo que permitió estimar la media poblacional y el error en la estimación de V. Se usaron 25 parcelas permanentes, que cubren el rango existente de calidad, edad y densidad en los rodales de la especie en la provincia y se ajustaron modelos lineales y potenciales. El mejor resultado se obtuvo con un modelo potencial entre V y el percentil 70 de la altura normalizada (h70) (R2 = 0,899; REMC = 63,4 m3 ha-1). Dicho modelo se aplicó a todas las teselas de P. radiata de la provincia de Lugo con estados de desarrollo de latizal y fustal y con ocupación de la especie superior al 60 % según el Mapa Forestal de España. Luego, se elaboró un mapa espacialmente explícito que describe V a escala provincial. En 2009 la media poblacional estimada para V fue de 166,3 m3 ha-1, con un error relativo de muestreo del 13,9 %. Esta metodología permitirá generar mapas temáticos cada seis años (periodicidad programada para los vuelos LiDAR del PNOA).
Palabras clave: láser escáner aerotransportado, inventario forestal, teledetección, estimadores, modelos de regresión.
La sociedad actual es cada vez más consciente del papel de los ecosistemas forestales en la protección de los suelos y otros ecosistemas, y en el control de la contaminación ambiental. Sin embargo, el papel ecológico de los montes debe ser considerado conjuntamente con su función económica, con el fin de asegurar la persistencia y viabilidad de todos los recursos forestales. La consecución de estos objetivos pasa por una optimización de la gestión forestal. Para ello, resulta necesario disponer de una información actualizada y de calidad de los recursos existentes, entre los cuales el volumen de madera es indispensable.
Hoy en día en España, la única información a escala nacional, autonómica y/o provincial acerca del volumen de madera disponible es la que proporciona el Inventario Forestal Nacional (IFN) en sus diferentes ediciones (IFN1, IFN2, IFN3 e IFN4). La estimación de las variables de rodal del IFN español se basa exclusivamente en inventarios sistemáticos en campo y en una fase previa de estratificación de las teselas (manchas de vegetación con características homogéneas) mediante fotointerpretación. Este tipo de inventarios no permite la plena cobertura del terreno y, a menudo, tardan varios años en completarse; solamente para la provincia de Lugo (NO de España), las labores de campo del IFN4 se prolongaron durante aproximadamente cinco meses, en los cuales se montaron un total de 2.382 parcelas. El procesado de los datos y la edición y publicación de los mismos demoró aproximadamente dos años y la estratificación previa se basó en ortoimágenes de hasta seis años antes, lo cual revela el elevado coste asociado y además, la información que contiene no refleja adecuadamente las existencias en el momento de su publicación.
En los inventarios forestales, la práctica más común para mejorar la precisión de las estimaciones de los parámetros poblacionales de las diferentes variables dasométricas ha sido aumentar las mediciones de campo, con el coste asociado que supone. Alternativamente, el uso de variables auxiliares cuyas observaciones estén correlacionadas con las observaciones de la variable respuesta (volumen en este caso) puede mejorar la estimación de los parámetros poblacionales en los inventarios forestales. Así, por ejemplo, las estimaciones se pueden hacer empleando modelos que relacionen variables respuesta con una o más variables auxiliares. Los sensores LiDAR aerotransportados (Airborne Light Detection and Ranging) han demostrado ser una interesante fuente de datos auxiliares útiles para inventario forestal (McRoberts et al. 2014), debido a su capacidad para medir de forma directa la estructura tridimensional de la vegetación y otras características de los bosques a diferentes escalas.
Actualmente, España cuenta con una gran cobertura de datos LiDAR, gracias al Plan Nacional de Ortofotografía Aérea (PNOA), que desde el año 2009 ha incorporado a sus productos fotogramétricos convencionales (ortofotografías aéreas digitales y modelos digitales de elevación - MDE) datos LiDAR para el control de la calidad de los MDE. De esta forma, se dispone de datos LiDAR de baja resolución (0,5 primeros retornos m-2 según las especificaciones técnicas del PNOA) para la mayor parte del territorio español, proporcionando de forma continua gran cantidad de datos de elevación georreferenciados, que pueden ser de mucha utilidad para la descripción cuantitativa y cualitativa de los sistemas forestales.
Existen dos niveles principales para realizar inventarios forestales con LiDAR: inventario a nivel de árbol individual con delineación de copas (Individual Tree Crown approach, ITC) e inventario a nivel de rodal (Area-Based Approach, ABA). Este último es el de menor coste y, por tanto, el más apropiado cuando se cuenta con datos LiDAR de baja densidad de pulsos por unidad de superficie, como suele suceder en los vuelos regionales o nacionales. En el estudio de González-Ferreiro et al. (2012) una reducción drástica de la densidad de datos LiDAR (de 8 a 0,5 pulsos m-2) no afectó de forma significativa a la estimación ABA de la mayor parte de las variables forestales de interés para la gestión forestal. Esta metodología establece relaciones empíricas entre las principales variables dasométricas a nivel de rodal y un conjunto de estadísticos y métricas LiDAR, habiendo sido utilizada para estimar el volumen en diferentes tipos de bosque (e.g. Cartus et al. 2012, González-Ferreiro et al. 2012, Guerra-Hernández et al. 2016a).
Los modelos desarrollados mediante ABA pueden utilizarse posteriormente para estimar de forma espacialmente explícita el volumen u otras variables forestales de interés, siendo posible generar mapas temáticos que muestren los resultados del modelo de predicción para cada una de las variables estudiadas (e.g. Cartus et al. 2012, Guerra-Hernández et al. 2016a) y para toda la extensión forestal cubierta por el vuelo. Si las observaciones de las variables auxiliares están disponibles para toda la población, es decir, existen datos LiDAR para toda la zona de estudio, el modelo generado se puede utilizar para predecir la variable de interés para todas las unidades poblacionales. Diversos estimadores inferenciales, entre ellos los estimadores basados en modelos de regresión, han sido desarrollados para este tipo de situaciones (McRoberts et al. 2014).
El objetivo principal de esta investigación es comprobar la validez de los datos LiDAR del PNOA para estimar el volumen de rodales regulares monoespecíficos mediante ABA y el uso de estimadores basados en modelos de regresión y comprobar que el error de estimación del volumen de masa es admisible. Para ello, se trabaja como ejemplo con la especie Pinus radiata D. Don en la provincia de Lugo (NO de España). La hipótesis de este trabajo es que el uso combinado de parcelas de muestreo (con un diseño muestral orientado a la mejora del modelo) y datos LiDAR de baja resolución puede proporcionar estimaciones de volumen por debajo del 20 % en volumen de masa con corteza, que es el error máximo admisible fijado a escala de monte por las instrucciones generales de ordenación y de gestión de montes de Galicia. En caso de verificarse la hipótesis de partida, el estudio tendría importantes implicaciones, ya que la metodología expuesta supondría la reducción del costo de los inventarios forestales a diversas escalas, desde los inventarios de gestión a los inventarios provinciales, regionales o nacionales, de una forma económicamente viable y fácilmente revisable, teniendo en cuenta que la periodicidad programada para los próximos vuelos LiDAR del PNOA es de seis años.
Área de estudio. Para delimitar la zona de estudio se emplearon, como unidades básicas, las teselas del Mapa Forestal Español (MFE), que son manchas de vegetación homogénea con una superficie mínima de 1 ha en el caso de zonas arbolada. La zona de estudio se corresponde con todas aquellas teselas situadas dentro de la provincia de Lugo que contienen P. radiata como especie principal, con una ocupación superior al 60 % y cuyas formaciones se encuentren en estado de desarrollo de latizal y fustal. Dicha zona ocupa 44.161,6 ha, repartidas en 2.712 teselas. Esto representa el 46,9 % de la superficie total de P. radiata de la provincia (figura 1).
Figura 1. Área de estudio: A) Localización de Galicia en España y en Europa. B) Masas de
Pinus radiata en Galicia según el IFN4. C) Masas de Pinus radiata en la provincia de Lugo
según el IFN4 (áreas en color negro y áreas en color gris), teselas de estudio (áreas en
color gris) y centroides de las 25 parcelas de campo establecidas en la provincia de Lugo
(puntos rojos).
Study area: A) Location of Galicia in Spain and Europe. B) Pinus radiata stands in Galicia
according to the 4NFI (4th National Forest Inventory). C) Pinus radiata stands in the province
Lugo according to 4NFI (black and grey areas), analyzed stands (grey areas) and centroids
of the 25 field sample plots in the province of Lugo (red dots).
Datos de entrenamiento. Los datos de campo empleados para desarrollar el modelo de predicción del volumen (datos de entrenamiento para el modelo) se obtuvieron entre marzo de 2009 y abril 2010 de dos fuentes independientes. La primera fuente (A), comprende una red de 10 parcelas rectangulares, de entre 600 y 1.000 m2 de superficie, establecidas por la Universidad de Santiago de Compostela para desarrollar un modelo dinámico de crecimiento para las plantaciones de P. radiata en Galicia. El diseño del inventario tuvo como principal objetivo representar adecuadamente el rango existente de edad, densidad y calidad de estación. La segunda fuente (B), comprende 15 parcelas rectangulares, de 1000 m2 de superficie, establecidas por la Universidad de León para investigar cómo se relacionan las variables de rodal con el riesgo potencial de que exista un incendio de copas. El inventario se diseñó con el fin de representar la variabilidad estacional de los rodales jóvenes de elevada densidad en la región.
En general, las parcelas escogidas se caracterizan por una alta densidad de plantación, una baja intensidad de tratamientos selvícolas y por la presencia de una carga moderada de matorral (ver González-Ferreiro et al. 2014, p. 352 y su bibliografía asociada para una completa descripción de las parcelas). La combinación de ambas fuentes de datos permitió cubrir el rango existente de calidad, edad y densidad en los rodales de la especie en la provincia de Lugo y su distribución espacial engloba la principal área de distribución de P. radiata en Galicia, que se corresponde esencialmente con la provincia de Lugo (figura 1). Todas las parcelas se localizaron en campo de forma precisa, midiendo las coordenadas UTM de sus cuatro esquinas con una estación total y un GPS diferencial.
En todos los árboles de cada parcela se midió el diámetro normal (a 1,3 m sobre el suelo, usando una forcípula con precisión de 1 mm) y la altura total (usando un hipsómetro con precisión de 1 dm). Posteriormente, se estimó el volumen con corteza del tronco de cada árbol utilizando el modelo expresado por la ecuación 1, que es apropiado para la especie y zona geográfica (Diéguez-Aranda et al. 2009):
Donde, v = volumen total con corteza del tronco (m3), d = diámetro normal con corteza (cm) y h = altura total del árbol (m).
El volumen de cada parcela se obtuvo agregando el volumen de los árboles individuales y refiriendo el resultado a la hectárea (V, m3 ha-1). También se calcularon las siguientes variables de rodal a partir de los datos de las parcelas: edad del rodal (t) calculada a partir de la fecha de plantación, número de árboles por hectárea (N), altura dominante (H, m, definida como la altura media de los 100 árboles con mayor diámetro normal por hectárea) e índice de sitio (S, m, definido como la altura dominante del rodal a la edad de referencia de 20 años, estimado con el modelo recogido en Diéguez-Aranda et al. 2009). La figura 2 muestra una matriz de diagramas de dispersión para las variables t, N, H, S y V. Los datos de campo representan el rango de calidades de estación y etapas de desarrollo para los turnos de corta comúnmente empleados para P. radiata en Galicia (como media 25 años en plantaciones privadas, Rodríguez et al. 2002).
Figura 2. Gráfico de dispersión matricial de las variables edad (t), densidad
(N, pies ha4), altura dominante (H, m, definida como la altura media de los 100
árboles con mayor diámetro normal por hectárea) e índice de sitio (S, m, definido
como la altura dominante del rodal a la edad de referencia de 20 años, estimado
con el modelo de Diéguez-Aranda et al. (2005)) y volumen de masa (V, m3 ha-1).
Los círculos sólidos representan los 10 datos de la red de parcelas A y las
circunferencias los 15 datos de la red de parcelas B.
Scatter plot graphic that represent the range of ages (t, ages), densities
(N, stem ha-1), dominant height (H, m, defined as the mean height of the 100
largest diameter trees per hectare) and site index (S, m, defined as the dominant
height at the reference age of 20 years, estimated using the model described in
Diéguez-Aranda et al. (2005)) and stand volume (V, m3 ha-1). Filled circles
represent field data from source A (10 plots). Empty circles represent field data
from source B (15 plots).
Datos LiDAR. Los datos LiDAR para la zona oriental de Galicia (provincias de Lugo y Ourense) se adquirieron en el proyecto PNOA entre el 5 de septiembre y el 29 de octubre de 2009 bajo la dirección del Ministerio de Fomento de la Administración General del Estado (a través de la Dirección General del Instituto Geográfico Nacional -IGN- y del Centro Nacional de Información Geográfica -CNIG) y la Consellería de Medio Ambiente, Territorio e Infraestructuras de la Xunta de Galicia (a través del Instituto de Estudos do Territorio -IET), empleando un sensor RIEGL LMS-Q680, operando a 1064 nm, con una frecuencia de repetición de pulsos de 70 Hz, una frecuencia de escáner de 46 Hz, un ángulo máximo de escaneado de 30° y una altura media de vuelo de 1.300 m sobre el nivel del elipsoide GRS80. Se registró un máximo de cuatro retornos por pulso, con una densidad teórica requerida por el proyecto PNOA de 0,5 primeros retornos m-2. Los estadísticos descriptivos de la densidad total de retornos por metro cuadrado dentro de las parcelas de entrenamiento fueron: media = 0,476, mínimo = 0,194, máximo = 1,098 y desviación estándar = 0,200.
Normalización de la intensidad de los datos LiDAR. Ya que la zona de estudio es extensa y el terreno muestra una alta variabilidad en cuanto a pendientes y altitudes, los valores de intensidad se normalizaron en función de un rango estándar, con el fin de eliminar la dependencia de los valores de intensidad con respecto a la altura de vuelo (García et al. 2010) y al ángulo de incidencia (Kaasalainen et al. 2011). Así, la intensidad normalizada (I') se obtuvo mediante la multiplicación del valor bruto de la intensidad (I) por el cociente entre el rango de cada retorno (rg, m) y el rango estándar (rgs, m) (en este caso 1.000 m) y por el inverso del coseno del ángulo de incidencia del rayo láser (α, radianes). El primer rango se calculó como la diferencia entre la altura media de vuelo (1.300 m sobre el nivel del elipsoide GRS80 en este caso) y la altura elipsoidal de cada retorno (m), es decir, la diferencia entre el sensor y el objeto con el que impacta el láser:
Extracción de las variables LiDAR. Se emplearon varios algoritmos implementados en LTK (LiDARTool Kit) del programa FUSIÓN V. 3.50 (McGaughey 2015) para el filtrado, interpolación y generación de los MDE, para la obtención de la nube de puntos LiDAR con altura normalizada y para el cálculo de los estadísticos relacionados con las distribuciones de altura y de intensidad de los retornos LiDAR en las 25 parcelas de entrenamiento (véanse los pasos descritos en González-Ferreiro et al. 2012 y la descripción de los estadísticos en los cuadros 1 y 2). Todas las métricas se calcularon a partir de los retornos por encima de 1 m para evitar retornos de matorral, arbustos, rocas, troncos, etc. Además, se calcularon un conjunto de métricas LiDAR relacionadas con el cierre de copas, usando varias ratios entre el número de retornos por encima de un umbral de altura de 2 m (cuadro 2).
Cuadro 1. Variables LiDAR, potencialmente explicativas, relacionadas con las distribuciones
de altura e intensidad LiDAR.
Potential explanatory LiDAR variables related with normalized height and normalized intensity
LiDAR distributions.
Nota: Las variables se calcularon usando todos los retornos LiDAR, es decir, primeros,
segundos, terceros y cuartos retornos.
Note: All LiDAR-derived variables were computed from all LiDAR returns in the database,
that is, 1st, 2nd, 3rd and 4th LiDAR returns.
Cuadro 2. Variables LiDAR, potencialmente explicativas, relacionadas el cierre de copas.
Potential explanatory LiDAR variables related with canopy closure.
Elaboración de los modelos. Se establecieron relaciones empíricas entre las variables obtenidas a partir de los datos medidos en campo y las métricas LiDAR, mediante el empleo de modelos de regresión múltiple. Las expresiones generales de los modelos son:
Donde, X1, X2..., Xk = variables explicativas del modelo, extraídas a partir de estadísticos y métricas de las distribuciones de altura e intensidad de las nubes de puntos LiDAR (cuadro 1) o variables LiDAR relacionadas con el cierre de copas (cuadro 2); β0, β1...., βk = parámetros a estimar en el proceso de ajuste; y ε = término aditivo del error en el modelo, que ha de ser normal, independiente e idénticamente distribuido, con media cero.
Todos los análisis estadísticos se realizaron con el programa estadístico R (R Core Team 2014). Con el fin de corregir la no linealidad en las relaciones entre algunas de las variables explicativas y la variable dependiente en el modelo 3, todas las variables potencialmente explicativas fueron transformadas usando una transformación potencial BoxCox. Para ello, se usó la función boxcox implementada en el paquete MASS (R Core Team 2014), que computa y opcionalmente representa el perfil de máxima verosimilitud para el parámetro X. La forma del modelo lineal con las variables independientes transformadas es:
Ajuste y selección del modelo. La selección del mejor conjunto de variables explicativas se realizó mediante una búsqueda con reemplazamiento secuencial aplicando el argumento sequential replacement de la función regsubsets, implementada en el paquete LEAPS (R Core Team 2014). Este método busca los subconjuntos de variables independientes que mejor predicen la variable dependiente por regresión lineal en una muestra dada. Se usó el criterio de información bayesiano (BIC) para ordenar los modelos candidatos de mejor a peor siguiendo un orden de BIC ascendente (Schwarz 1978). En el caso del modelo 4, para poder aplicar la técnica anterior se tomaron logaritmos en ambos términos de la igualdad para linealizarlo y, una vez seleccionado el mejor conjunto de variables explicativas, se ajustó el modelo potencial original con las mismas variables.
Los parámetros se estimaron por mínimos cuadrados ordinarios. Para analizar el comportamiento de los modelos se realizaron análisis numéricos y gráficos de los residuos. Se utilizaron dos estadísticos de bondad de ajuste: el coeficiente de determinación (R2) y la raíz del error cuadrático medio (REMC). Aunque existen limitaciones asociadas con el uso del R2 en regresión no lineal, su uso común como medida global de la adecuación de un modelo anula estos inconvenientes (Ryan 1997, p. 424). El R2 indica la proporción de la varianza total de la variable dependiente que es explicada por el modelo, pero no debe ser utilizado como criterio único en la selección del mejor modelo (Myers 1990, p. 166). El RMSE es útil porque proporciona una idea de la precisión de la estimación en las mismas unidades que la variable dependiente. Además, penaliza a los modelos con más parámetros, de acuerdo con el principio general de simplicidad científica.
Por último, se estudió la multicolinealidad entre variables calculando el índice de condición mediante la función colldiag del paquete PERTURB (R Core Team 2014). De acuerdo con Belsley (1991), se rechazaron los modelos de regresión con índice de condición superior a 30. Finalmente, solo se consideraron los modelos con estimaciones de los parámetros β0, β1...., βk significativas a un nivel del 5 %.
Aplicación del modelo y obtención del volumen. Para extraer las métricas LiDAR seleccionadas como independientes en el modelo de regresión en un formato ráster se emplearon varios algoritmos implementados en LTK del programa FUSION V 3.50 (McGaughey 2015). El primer paso consistió en emplear el algoritmo PolyClipData, para extraer los datos LiDAR contenidos en las teselas de estudio. Para ello se empleó el archivo espacial en formato shapefile procedente del MFE (figura 1), al cual previamente se le aplicó un buffer de 50 m de anchura, con el fin de evitar errores de cálculo en los bordes de las teselas en procesos subsecuentes de normalización de la altura de los datos LiDAR. En segundo lugar, se empleó el algoritmo GroundFilter para extraer de la nube de puntos LiDAR aquellos retornos que pertenecen al suelo y a partir de ellos generar los MDE de las teselas, interpolando con una resolución de 3 m con el algoritmo GridSurfaceCreate. En tercer lugar, se usó el algoritmo GridMetrics para extraer diferentes métricas y estadísticos de las distribuciones normalizadas de altura e intensidad de los datos LiDAR usando un tamaño de píxel de 30 m de lado. Además, GridMetrics permite excluir valores extremos; en este caso se descartaron los valores de altura normalizada fuera del rango de 1 a 40 m, para evitar en la parte inferior retornos pertenecientes al matorral, troncos tirados, piedras, etc. o retornos afectados por el efecto de retrodispersión múltiple, y en la parte superior retornos provenientes de insectos, pájaros, tendidos eléctricos, etc., considerándose 40 m la altura máxima del arbolado en la zona de estudio. En cuarto lugar se utilizó el algoritmo CSV2Grid, que proporciona un ráster espacial-mente explícito del estadístico o estadísticos seleccionados como variable o variables independientes del modelo de estimación; de este modo, se obtuvo un ráster de salida por tesela y estadístico. Por último, se usó el algoritmo MergeRaster para fusionar los archivos correspondientes a las diferentes teselas y obtener un único ráster por cada estadístico seleccionado como explicativo en el modelo.
Una vez seleccionado el modelo de predicción del volumen de rodal, se procedió a su aplicación, usando como variable o variables explicativas el ráster o rásteres obtenidos en el proceso anterior. Para ello, en primer lugar se leyó el archivo o archivos ráster, se aplicó la ecuación del modelo a cada píxel y el volumen obtenido se exportó a formato ráster. Todos estos procesos se realizaron con las funciones scan, ifelse y write.table del paquete BASE (R Core Team 2014). De este modo, se obtuvo un archivo espacialmente explícito con un tamaño de pixel de 30 m de lado, que describe las existencias del volumen de cada rodal de P. radiata seleccionado en la provincia de Lugo. Finalmente, se procedió a estimar los principales estadísticos resumen para la variable volumen de rodal.
Determinación del error de muestreo. Para comprobar la precisión de las estimaciones se recurrió al cálculo del error de las estimaciones basadas en modelos de regresión. La inferencia en los estimadores basados en modelos de regresión no se fundamenta en el diseño muestral con el que han sido elegidas las unidades de muestreo, sino en el modelo empírico que relaciona la variable primaria de interés con una segunda variable (variable suplementaria) en dichas unidades de muestreo (Gregoire 1998). Para el uso de estos estimadores se requiere el conocimiento de los valores de la variable suplementaria para todos los elementos del dominio de la población donde se pretenden hacer estimaciones.
Para la aplicación de procedimientos de estimación basados en modelos se asume que la población objeto de estudio es finita y que se dispone de una medición de la variable auxiliar x para todos los elementos de la población (N). Por su parte, únicamente se dispone de medición de la variable primaria y para los n elementos que componen la muestra. La variable primaria y secundaria se relacionan mediante el modelo ajustado (genéricamente i = f (Xi, )) el cual permite estimar para cada i elemento de la población (i = 1, 2, N) el valor de la variable primaria ( i). El estimador de la media basado en modelos de regresión es (McRoberts 2006):
Donde, i = estimación del valor de la variable primara para cada elemento de la población.
El estimador correspondiente del error estándar es:
Donde, los elementos de la matriz Z son zˑk = ϑf (X., )/ϑβk siendo βk los parámetros del modelo, y σ2{ } es la matriz de varianzas-covarianzas de la estimación de los parámetros k del modelo (McRoberts 2006).
Debido a que la estimación del error estándar mediante la ecuación 7 es computacionalmente muy intensiva, se ha aproximado el valor de S MOD mediante un muestreo aleatorio simple sin reemplazamiento de n = 10.000 elementos (McRoberts et al. 2013). Además, se estimaron el error absoluto (E) y el error relativo (ε):
Donde, t = valor de la t de Student para una probabilidad del 95 % y n-2 grados de libertad.
Análisis de regresión. La transformación potencial Box-Cox permitió corregir la falta de linealidad en la relaciones entre la variable explicativa y la variable dependiente en el modelo 3, aunque se obtuvieron resultados ligeramente mejores empleando el modelo potencial (ecuación 4). El percentil del 70 de la distribución de alturas (h70) fue la única variable explicativa seleccionada en ambos modelos. La variabilidad observada explicada por el modelo potencial fue del 89,9 % con un valor de REMC relativo del 18,0 %. Todos los parámetros del modelo potencial resultaron significativos a un nivel del 5 % y sus estimaciones, estadísticos de bondad del ajuste y covarianza se muestran en el cuadro 3. Los diagramas de dispersión indican que la variable predicha (V) está fuertemente correlacionada con la variable independiente y que no existen tendencias observadas de sobreestimación o subestimación (figura 3).
Cuadro 3. Estimación de los parámetros, estadísticos de bondad del ajuste y covarianza
del modelo seleccionado para estimar volumen de rodal (V, m3 ha-1) a partir de datos LiDAR
(el percentil del 70 de la distribución de alturas, h70).
Parameter estimates, goodness-of-fit statistics and covariance for the selected model used
to estimate stand volume (V, m3 ha-1) from LiDAR data (70 percentile of the LiDAR normalized
height distribution, h70).
Figura 3. Gráfica de dispersión que
representa la variable independiente
seleccionada (h70 m) y la variable
dependiente (V, m3 ha-1). La curva
representa el mejor ajuste potencial
que se describe en el cuadro 3.
Scatter plot representing the chosen
explanatory variable (h70 m) and the
dependent variable (V, m3 ha-1). The
curve represents the best potential
adjustment described in table 3.
El modelo ajustado cumplió todos los supuestos necesarios para realizar regresión paramétrica (independencia de los errores, homocedasticidad y normalidad de la distribución de errores) y a la vista de los resultados permitirá realizar predicciones de V e inferencia en el área de estudio a partir de los datos LiDAR del PNOA.
Volúmenes obtenidos y error de muestreo. La metodología de estimadores basados en modelos de regresión propuesta, en la que se tomó el volumen medido en campo como variable dependiente y la variable h70 calculada a partir de datos LiDAR como independiente, proporcionó una media poblacional estimada ( MOD) para la zona de estudio de 166,3 m3 ha-1. El error estándar de la media poblacional de V (S MOD), el error absoluto (E) y el error relativo (ε) obtenidos a partir de la muestra aleatoria de 10000 elementos fueron 11,7 m3 ha-1, 23,1 m3 ha-1 y 13,9 %, respectivamente.
La aplicación del modelo de estimación del volumen de rodal permitió elaborar un mapa temático espacialmente explícito (figura 4) que representa la distribución espacial de las existencias maderables de P. radiata a escala provincial. La figura 5 representa una muestra del cálculo del volumen obtenido a nivel de tesela.
Figura 4. Representación gráfica de las existencias en volumen sobre píxeles de 30 m
Graphical representation of the stocks of timber volume using pixels of 30 m of side.
Figura 5. Existencias maderables de P. radiata estimadas con LiDAR.
Timber stocks estimation of Pinus radiata stands using LiDAR.
La adquisición de datos de campo para inventario forestal es generalmente costosa tanto en tiempo como en recursos. En España, aproximadamente el 50 % de los costes de inventario provienen de los trabajos de campo, de modo que el coste de un inventario tradicional con fines de gestión forestal es de aproximadamente 20 € ha-1 (Tomé Morán et al. 2013).
Los inventarios forestales basados en datos LiDAR a nivel de rodal (es decir ABA) se están convirtiendo en una realidad rentable en inventarios forestales a gran escala en zonas boreales, pero siguen siendo un reto en las regiones templadas de Europa en las que, a menudo, diferentes especies arbóreas están presentes en pequeñas propiedades forestales que además presentan una gran dispersión geográfica y heterogeneidad. Aun así, Fabra (2012) afirma que, en España, el inventario forestal con datos LiDAR puede abaratar los costes hasta aproximadamente los 5 € ha-1. Además, la posibilidad de contar con datos LiDAR a nivel nacional puede reducir más aún los costes y fomentar el uso múltiple de los datos LiDAR (los costes unitarios disminuyen a medida que aumenta la superficie escaneada y el número de objetivos propuestos para cada vuelo) (González-Ferreiro et al. 2014). En el caso de los datos a nivel nacional del PNOA, las consultas realizadas al Instituto Geográfico Nacional Español (IGN) revelaron que, para grandes adquisiciones (900.000 ha), el coste de los datos LiDAR es de aproximadamente 0,14 € ha-1 en la mayor parte de los casos (duplicándose en el caso de adquisiciones centradas en áreas muy perfiladas e irregulares o pequeñas islas del territorio). De forma adicional, este coste debe ser imputado a los diferentes usos que se pueden dar a los datos LiDAR, reduciendo de este modo el coste global de los inventarios forestales.
Los datos LiDAR del PNOA han demostrado su utilidad para la estimación de un gran conjunto de variables de rodal como el volumen (Guerra-Hernández et al. 2016a), el área basimétrica (Guerra-Hernández et al. 2016a), la altura media de Lorey (González-Ferreiro et al. 2014, Guerra-Hernández et al. 2016a), variables del complejo de combustibles de copa (González-Ferreiro et al. 2014), la severidad del fuego (Montealegre et al. 2014), la biomasa o el carbono (Montealegre et al. 2015, Guerra-Hernández et al. 2016ab), pero no han sido usados aún en España para el inventario forestal de grandes superficies, como la provincia de Lugo.
A continuación se discuten los dos principales procesos metodológicos que se abordan en este estudio: el modelo de regresión y la estimación del volumen a partir del modelo.
Modelo de regresión. En este estudio se han usado datos LiDAR procedentes del PNOA y una metodología basada en el análisis estadístico de las relaciones entre el volumen de rodal medido en campo y una serie de variables regresoras obtenidas de las distribuciones de altura e intensidad de la nube de puntos LiDAR normalizada. La precisión de los MDE y los parámetros de vuelo (densidad de muestreo, altura de vuelo y ángulo de escaneo) son factores muy influyentes en las métricas LiDAR. Sin embargo, y de acuerdo con lo expresado por Guerra-Hernández et al. (2016a), se han ajustado modelos que explican un elevado porcentaje de la variabilidad total del volumen de rodal en masas de Pinus pinaster Ait. (hasta un 98 %), Quercus pyrenaica Willd. (hsta un 88 %), Pinus pinea L. (hasta un 74 %) y masas mixtas de P. pinea, P. pinaster y Q. pyrenayca (hasta un 82 %) utilizando los datos LiDAR del PNOA, con parámetros de vuelo que podrían ser a priori limitantes, como ángulos de escaneo de hasta 50°, densidad de pulsos LiDAR de tan solo 0,5 primeros retornos por m2 y alturas de vuelo de hasta 2.866 m sobre el elipsoide de referencia.
Para discutir el comportamiento del modelo con respecto a otros estudios en masas de P. radiata se han usado los estadísticos R2 y REMC, debido a su uso generalizado como indicadores de la adecuación de un modelo, así como la densidad de pulsos LiDAR, ya que el uso de densidades de pulsos LiDAR elevadas puede compensar o reducir la mayor parte de las fuentes de error provocadas por la elección del resto de los parámetros de vuelo (Hyyppä et al. 2008). El modelo seleccionado para V (R2 = 0,899; REMC = 63,4 m3 ha-1) se ha comportado significativamente mejor que los modelos presentados por González-Ferreiro et al. (2012) para masas de P radiata en el NO de España, con densidades de muestreo LiDAR de entre 0,5 y 8 pulsos por m2, obtenidas a partir de una reducción artificial de la base de datos LiDAR original (R2 = 0,68 - 0,79; REMC = 76,9 - 94,1 m3 ha-1) y ligeramente mejor que los modelos presentados por Cartus et al. (2012) para plantaciones de P. radiata y Eucalyptus globulus Labill. en la zona central de Chile (R2 = 0,70 - 0,85 para densidades de 1 y 3 pulsos m-2, respectivamente).
El modelo seleccionado para V incluye el percentil del 70 de la distribución de alturas (h70) como única variable explicativa, lo que concuerda con las averiguaciones realizadas por Gobakken y Nssset (2007), que constataron que los percentiles intermedios y elevados de la distribución de alturas LiDAR (h50 a h90) permanecen muy estables incluso con bajas densidades de retornos LiDAR, por lo que es de esperar que el modelo presentado tenga un comportamiento robusto, aún con densidades de retornos LiDAR muy bajas.
Por otra parte, teniendo en cuenta la gran extensión de la zona de estudio, la amplia variabilidad en cuanto a estructuras de masa cubiertas por las parcelas de campo usadas para crear el modelo y que estas se distribuyeron a lo largo de la principal área de distribución de P. radiata en Galicia, es de esperar que el modelo empírico aquí propuesto pueda ser aplicado a la mayor parte de los rodales de P radiata en la región. Aun así, se debería ser cauteloso cuando se trate de extrapolar las estimaciones del modelo fuera del dominio de valores de los datos de campo usados en este estudio, de forma particular en masas fuertemente aclaradas o muy jóvenes.
Estimación del volumen. En este estudio se ha caracterizado la localización espacial de las existencias maderables para los rodales de P. radiata en estado de latizal y fustal con una ocupación superior al 60 % y dicha información se ha recogido en un mapa que describe V en cada píxel de 30 m de lado. El tamaño del píxel ha sido escogido según las recomendaciones de Laes et al. (2011), de forma que ocupase una superficie próxima a la de las parcelas usadas en la fase de modelización (entre 600 y 1.000 m2 en este caso). También se consideró la densidad de los datos LiDAR utilizados, con el fin de contar con un suficiente número de retornos por celda que proporcionasen métricas LiDAR no sesgadas y la heterogeneidad de las masas forestales analizadas, con el fin de obtener una buena representación de las mismas.
Las estimaciones poblacionales en el ámbito del inventario forestal tradicionalmente han sido realizadas mediante técnicas de inferencia basadas en el diseño de la selección de las unidades de muestreo (Cochran 1999). Sin embargo, estas técnicas no son demasiado flexibles en cuanto al uso de variables auxiliares que pudieran contribuir a mejorar la precisión de la estimación, puesto que en todos los casos requieren que la muestra tenga un diseño muestral conocido (Gregoire 1998). En este trabajo se han empleado estimadores basados en modelos de regresión, los cuales permiten hacer uso de parcelas ya establecidas sin diseño muestral previo, conjuntamente con la observación de una variable auxiliar sobre todos los elementos de la población objetivo.
Una ventaja añadida de esta metodología es que permite estimar la media y su error asociado para delimitaciones concretas del dominio de estudio ("small area estimation" Chambers y Clark 2012), como pueden ser el rodal, monte o estrato (e.g. Breidenbach et al. 2016). Sin embargo, a pesar de la flexibilidad que otorgan estos métodos, para su uso se debe poner especial énfasis en las propiedades estadísticas de los modelos, ya que el peso de la inferencia recae en la calidad del ajuste y en su capacidad explicativa, junto con la variabilidad de la variable suplementaria en el dominio donde se pretenden hacer estimaciones (Gregoire 1998, McRoberts 2006).
Se debe tener en cuenta que las teselas identificadas en el MFE deben tener una extensión mínima de 1 ha en el caso de zonas arboladas, lo cual hace que en zonas forestales altamente heterogéneas se haga necesario englobar masas que difieren mucho en su estructura y composición florística, situación común en Galicia debido al predominante minifundismo presente en la región, en la que el 64 % de la propiedad forestal es privada y el tamaño medio de superficie forestal por propietario es de 1,78 ha distribuidas en varias parcelas, según el Plan Forestal de Galicia en el año 1992.
Aunque las teselas de la zona de estudio (figura 1) han sido filtradas para contener exclusivamente masas de P. radiata en estado de latizal y fustal con una ocupación superior al 60 %, los píxeles que las constituyen pueden incluir datos de zonas rasas y arbolado joven. Esto explicaría en parte la diferencia entre el valor medio del volumen de la población y el de las parcelas utilizadas para ajustar el modelo.
Finalmente, el modelo propuesto para estimar el volumen a nivel de rodal podría utilizarse regularmente cada seis años, que es la periodicidad programada para los vuelos LiDAR del PNOA. En la actualidad se han licitado futuros vuelos, y es de esperar que los nuevos datos LiDAR estén disponibles en un breve periodo de tiempo. La posibilidad de utilizar los datos LiDAR para estimar las existencias en los próximos IFN supondría reducir sustancialmente su periodicidad actual de 10 años.
Los resultados obtenidos a nivel provincial sugieren que los datos LiDAR a escala nacional del PNOA, con una resolución espacial de aproximadamente 0,5 primeros retornos por metro cuadrado y una resolución temporal programada de seis años, pueden ser útiles para la estimación del volumen de rodal a otras escalas de mayor detalle (rodal, estrato o monte), ya que los errores relativos obtenidos mediante la aplicación de estimadores basados en modelos de regresión para la provincia de Lugo son claramente inferiores a los máximos admisibles según la legislación vigente en la Comunidad Autónoma de Galicia (20 % para cada estrato y 40 % para el monte en todas las masas que previsiblemente entrarán en corta en un decenio), lo cual permitiría reducir el coste de los inventarios forestales haciendo uso de información adquirida mediante sensores remotos en proyectos multipropósito.
Agradecemos la financiación de: 1) el Ministerio Español de Ciencia e Innovación (AGL2008-02259/FOR); 2) el Gobierno de Galicia -Xunta de Galicia, Dirección Xeral de Montes (09MRU022291PR); y 3) la Multinacional energética Norvento (PGIDT09REM023E). También agradecemos la financiación posdoctoral del Dr. Eduardo M. González Ferreiro por parte del Gobierno Gallego y el Fondo Social Europeo (Diario Oficial de Galicia - DOG n° 52, p. 11343, fecha: 17/03/2014), y del Dr. César Pérez-Cruzado por parte de la Secretaría de Estado de Investigación, Desarrollo e Innovación (Juan de la Cierva - Incorporación).
Belsley DA. 1991. A guide to using the collinearity diagnostics. Computer Science in Economics and Management 4(1): 33-50.
Breidenbach J, RE McRoberts, R Astrup. 2016. Empirical coverage of model-based variance estimators for remote sensing assisted estimation of stand-level timber volume. Remote Sensing of Environment 173: 274-281.
Cartus O, J Kellndorfer, M Rombach, W Walker. 2012. Mapping canopy height and growing stock volume using airborne LiDAR, ALOS PALSAR and Landsat ETM+. Remote Sensing 4(11): 3320-3345.
Chambers RL, RG Clark. 2012. An introduction to model-based survey sampling with applications. New York, USA. Oxford University Press. 265 p.
Cochran WG. 1999. Sampling techniques. New York, USA. Willey. 428 p.
Diéguez-Aranda U, A Rojo Alboreca, F Castedo-Dorado, JG Alvarez González, M Barrio-Anta, F Crecente-Campo, JM González, C Pérez-Cruzado, R Rodríguez-Soalleiro, CA López-Sánchez, MA Balboa-Murias, JJ Gorgoso, F Sánchez. 2009. Herramientas selvícolas para la gestión forestal sostenible en Galicia. Santiago de Compostela, España. Xunta de Galicia. 259 p.
Fabra M. 2012. Aplicaciones de la tecnología LiDAR al sector forestal y comparación de costes frente a metodologías tradicionales. Revista Montes 110: 33-37.
García M, D Fiano, E Chuvieco, FM Danson. 2010. Estimating biomass carbon stocks for a Mediterranean forest in central Spain using LiDAR height and intensity data. Remote Sensing of Environment 114(4): 816-830.
Gobakken T, E Nassset. 2007. Assessing effects of laser point density on biophysical stand properties derived from airborne laser scanner data in mature forest. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences 36(3): 150-155.
González-Ferreiro E, U Diéguez-Aranda, F Crecente-Campo, L Barreiro-Fernández, D Miranda, F Castedo-Dorado. 2014. Modelling canopy fuel variables for Pinus radiata D. Don in NW Spain with low-density LiDAR data. International Journal of Wildland Fire 23(3): 350-362.
González-Ferreiro E, U Diéguez-Aranda, D Miranda. 2012. Estimation of stand variables in Pinus radiata D. Don plantations using different LiDAR pulse densities. Forestry 85(2): 281-292.
Gregoire TG. 1998. Design-based and model-based inference in survey sampling: appreciating the difference. Canadian Journal of Forest Research 28(10): 1429-1447.
Guerra-Hernández J, M Tomé, E González-Ferreiro. 2016a. Using low density LiDAR data to map Mediterranean forest characteristics by means of an area-based approach and height threshold analysis. Revista de Teledetección 46: 103-117.
Guerra-Hernández J, E Bastos-Görgens, J García-Gutiérrez, LC Estraviz-Rodriguez, M Tomé, E González-Ferreiro. 2016b. Comparison of ALS based models for estimating above-ground biomass in three types of Mediterranean forest. European Journal of Remote Sensing 49: 185-204.
Hyyppä J, H Hyyppä, D Leckie, F Gougeon, X Yu, M Maltamo. 2008. Review of methods of small-footprint airborne laser scanning for extracting forest inventory data in boreal forests. International Journal of Remote Sensing 29(5): 1339-1366.
Kaasalainen S, U Pyysalo, A Krooks, A Vain, A Kukko, J Hyyppä, M Kaasalainen. 2011. Absolute radiometric calibration of ALS intensity data: Effects on accuracy and target classification. Sensors 11(11): 10586-10602.
Laes D, SE Reutebuch, RJ McGaughey, B Mitchell. 2011. Guidelines to estimate forest inventory parameters from LiDAR and field plot data. Companion document to the Advanced LiDAR Applications-Forest Inventory Modeling class. Salt Lake City, USA. US Forest Service. 22 p.
McGaughey R. 2015. Fusion/LDV version 3.50: Software for LiDAR data analysis and visualization. Portland, USA. US Forest Service. 182 p.
McRoberts RE. 2006. A model-based approach to estimate forest area. Remote Sensing of Environment 103(1): 56-66.
McRoberts RE, E Nassset, T Gobakken. 2013. Inference for lidar-assisted estimation of forest growing stock volume. Remote Sensing of Environment 128: 268-275.
McRoberts RE, HE Andersen, E Nassset. 2014. Using airborne laser scanning data to support forest sample surveys. In Maltamo M, E Nasset, J Vauhkonen eds. Forestry applications of airborne laser scanning: Concepts and case studies. London, UK. Springer. p. 269-292.
Myers RH. 1990. Classical and Modern Regression with Applications. Belmont, USA. Duxbury Press. 488 p.
R Core Team. 2014. R: A language and environment for statistical computing. Vienna, Austria. R foundation for statistical computing. Consultado 27 jul. 2015. Disponible en http://www.R-project.org/
Rodríguez R, F Sánchez, J Gorgoso, F Castedo, C López, KV Gadow. 2002. Evaluating standard treatment options for Pinus radiata plantations in Galicia (north-western Spain). Forestry 75(3): 273-284.
Ryan TP. 1997. Modern Regression Methods. New York, USA. Willey. 515 p.
Schwarz G. 1978. Estimating the dimension of a model. Annals of Statistics 6(2): 461-464.
Tomé Morán JL, P Sanjuanbenito García, A Fernández Landa. 2013. Cartografía de Vegetación en la Comunidad de Madrid utilizando información LiDAR del Plan Nacional de Ortofotografía Aérea (PNOA). In VI Congreso Forestal Español de la Sociedad Española de Ciencias Forestales 6CFE01-421. Vitoria-Gasteiz, España. Sociedad Española de Ciencias Forestales. p. 2-14.
Recibido: 30.03.16
Aceptado: 23.08.16

References: resolución

 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución