INTRODUCCIÓN
El calentamiento global es el aumento de la temperatura media de la superficie terrestre como resultado del aumento de concentraciones de gases de efecto invernadero (GEI). Es bien conocido que la liberación de CO2 a través de diversas fuentes es uno de los causales de este calentamiento. A través del proceso de fotosíntesis los bosques capturan CO2, lo fijan en sus estructuras vivas, lo acumulan en su biomasa y lo transfieren al mantillo en descomposición, y así se constituyen en reservas de carbono [1]. Esta importancia mundial de los ecosistemas forestales insiste en la necesidad de determinar con precisión la cantidad de carbono almacenado en los bosques.
Existen diferentes métodos convencionales para cuantificar la captura de CO2, muchos de los cuales son complejos, costosos y de cobertura limitada [2]. A este respecto, la posibilidad de estimar la biomasa y carbono a través de técnicas de teledetección representa una alternativa importante para sortear estas limitaciones, y ha sido ampliamente reportada en estudios a diversas escalas y ambientes especiales [3-5]. Para estimaciones de almacenamiento de CO2, IPCC [6] sugiere aplicar una fracción del 50 % de la biomasa aérea.
La estimación de la biomasa es una tarea difícil, especialmente en zonas con estructuras de condiciones ambientales complejas [7]. Los modelos alométricos han sido desarrollados para aplicaciones en bosques naturales y plantaciones forestales, pero su uso reviste limitaciones debido a condiciones naturales locales -tales como forma del árbol y arquitectura del dosel-, lo que sugiere desarrollar modelos locales de estimación [8-9]. A este respecto, los métodos de estimación de biomasa y carbono a partir de sensores remotos necesitan la integración de información de campo, basada en un muestreo espacial suficiente, y el uso de las ecuaciones alométricas apropiadas para mejorar su estima y asegurar su robustez [10-11]
Los vehículos aéreos no tripulados (unmanned aerial vehicles, UAV) equipados con sensores espectrales ligeros facilitan el análisis de vegetación con procedimientos no destructivos [12]. En los últimos años, los sistemas UAV han llamado la atención de las geociencias debido a la posibilidad de capturar datos de forma rentable a altas resoluciones espaciales y temporales [13]. En los trabajos con UAV donde se busque discriminar entre especies individuales de vegetación, Torres-Sánchez et al. [14] recomiendan resoluciones espaciales inferiores a cuatro centímetros que se logran con alturas de vuelo inferiores a 100 metros. A fin de corregir el posicionamiento geográfico de los elementos capturados en la escena, es necesario establecer puntos de control terrestre (GCP) distribuidos de manera uniforme en la zona de trabajo [15-16].
La atmósfera es la fuente primaria de ruido para la medición precisa de la reflectancia superficial con sensores remotos ópticos [17-18]. Los efectos de la atmósfera sobre los índices de vegetación son relevantes: Agapiou et al. [19] encontraron diferencias significativas de hasta 15 % entre hacer o no correcciones atmosféricas en el cálculo del NDVI. Homem et al. [20], a su turno, evaluaron los efectos atmosféricos en el NDVI y el SAVI utilizando el algoritmo de corrección 6S: encontraron que para la obtención de valores ajustados de estos índices es necesario corregir atmosféricamente las imágenes; así mismo, concluyen que estos efectos tienen alta incidencia cuando se busca relacionar los índices de vegetación con variables biofísicas como LAI y biomasa, por lo cual recomiendan efectuar mediciones de campo para ajustar los modelos.
El dosel forestal es la frontera entre la atmosfera y la superficie terrestre; la comprensión de su estructura es importante para modelar la captura de gases atmosféricos como el CO2 [21]. El LAI es un atributo estructural de los ecosistemas forestales a través del cual se han realizado investigaciones en estimativos de producción forestal y capacidad de captura de carbono [22], su correcta estimación es esencial para impactos de ecosistemas sobre el clima, el ciclo biogeoquímico y el balance energético [23]. Para medir de manera precisa y eficiente el LAI y Fapar se utiliza el instrumento óptico que recibe el nombre de Tracing Radiation and Architecture of Canopies (Trac): a través de la interceptación de la luz que pasa por el dosel, logra definir su estructura [24].
Con el trabajo presentado aquí se pretende calcular la capacidad de captura de CO2 del relicto boscoso "Jardín Botánico Cedro Rosado" de la Universidad del Quindío mediante la integración de técnicas de teledetección con datos de campo. Para lograr este objetivo, en la fase inicial del trabajo se parte de la compactación del ortomosaico de la zona de estudio con los valores de reflectancia necesarios para calcular los índices de vegetación NDVI, EVI y SAVI. Será necesario discriminar las diversas especies vegetales a través de la clasificación de imagen basada en objetos (OBIA), lo que llevará a definir las áreas de mediciones in situ de variables arbóreas. En la segunda fase, para el cálculo de las variables biofísicas LAI y Fapar, se utilizará el instrumento óptico Trac, el cual registra la transmisión de luz directa del dosel a través de los transectos diseñados. Finalmente, para encontrar las constantes del modelo exponencial de regresión propias para calcular biomasa, se medirá el grado de correlación de las variables obtenidas con los datos de campo, de modo que se pueda hacer la interpolación de la variable dependiente para la totalidad de la zona de estudio.
1. MATERIALES Y MÉTODOS
1.1 Caracterización biológica de la zona de estudio
El Jardín Botánico Cedro Rosado es un bosque secundario ubicado en el campus principal de la Universidad del Quindío (Municipio de Armenia, Depto. del Quindío, Colombia), con posición 4°33'N , 75°39'W; comprende, principalmente, vegetación natural intervenida con manchas de bosque secundario y guaduales. De acuerdo con la clasificación ecológica de zonas de vida expuestas por Holdridge [25], el espacio que nos ocupa corresponde a bosque muy húmedo premontano (BMH - PM), cuya superficie aproximada es de 13 Ha y se encuentra a una altitud que oscila entre 1490 y 1530 msnm. Su temperatura promedio asciende a 19,5 °C; su precipitación, a 2436 mm/año; y su humedad relativa varía entre 65-75 %. Posee dos microcuencas que, en conjunto, forman la Quebrada del Cedro Rosado [26]. La diversidad biológica encontrada incluye 64 especies de aves; 53 de hongos macroscópicos; 3 de anfibios; 2 de reptiles; 7 de mamíferos; y 247 de plantas (230 nativas y 17 cultivadas)
1.2 Diseño del plan de vuelo
Esta labor definió la ruta a seguir por la aeronave desde su despegue hasta su aterrizaje. La misión se planificó con el software Map Pilot® instalado en la estación de cómputo terrestre. La línea de vuelo, a su turno, se diseñó sobre las imágenes del área de trabajo importadas desde la aplicación Google Earth®, y para planificarla fueron necesarios tres parámetros: zona de estudio, especificaciones de la cámara y tareas propias del UAV. La información del área de vuelo incluye el ancho y la longitud de la zona de estudio, el ángulo de dirección del lado principal y el traslapo deseado en las imágenes, capturadas con 60 % de recubrimiento longitudinal y 30 % transversal. Se empleó la cámara Survey 2 fue empleada para registrar longitudes de onda del infrarrojo cercano en la vegetación, con una resolución espectral de 760 a 840 nm, y una resolución radiométrica de 24 bits; y con la cámara Phantom 3 Advanced se registraron los valores de RGB con una resolución radiométrica de 8 - 10 bits [27]. Las tareas del UAV son las acciones que el sensor tiene que realizar una vez llega a cada punto de captura de la imagen; aquí se incluyen el número de fotos y el tiempo de permanencia en cada punto. Ingresadas estas variables y la altura del vuelo, el aplicativo estableció de forma automática la ruta autónoma de vuelo y estimó su duración de acuerdo con el número total de imágenes planificadas.
1.3 Vuelo UAV y adquisición de imágenes
Una vez diseñada la misión de vuelo, se cargó la ruta en la plataforma Map Pilot. Luego, se realizó el montaje de las cámaras al vehículo que registraron información de manera simultánea: una registraría radiación visible en el sensor RGB; mientras que la otra poseía un filtro que reflejaba la luz del infrarrojo cercano en el canal azul del sensor RGB. Mediante el transmisor de radiocontrol se ejecutó el lanzamiento de la aeronave activando la ruta de vuelo trazada para que el vehículo llegara hasta el punto de inicio de la primera fotografía y se desplazara por las líneas de vuelo hasta cubrir por completo la zona estudiada.
1.4 Procesamiento de restitución para compactación del ortomosaico
Para iniciar el procedimiento de restitución por coincidencia espacial entre elementos de cada imagen se utilizó el software especializado Agisoft PhotoScan Pro®. A través de esa herramienta se alienaron imágenes para identificar puntos coincidentes entre ellas, de tal suerte que se compuso una nube de puntos para definir el modelo; y luego se definieron los puntos de control terrestre pertenecientes a la red topográfica de alta precisión que se encuentra en el campus de la universidad, a fin de georreferenciar el ortomosaico.
1.5 Preprocesamiento digital de imágenes
El modelo de procesamiento digital de la imagen exige dos tareas fundamentales para obtener reflectancias corregidas, esto es, las calibraciones radiométrica y atmosférica. Los procedimientos de calibración radiométrica aplicados a la imagen tienen en cuenta los errores que afectan el valor del brillo de cada pixel de la escena debido a dos fenómenos fundamentales: el error del detector del sistema sensor y el de atenuación ambiental [28]. El sistema MAPIR contiene un aditamento especial con tres objetivos, de los cuales se conoce su curva de reflectancia. Antes del vuelo se toman fotografías con cada cámara al conjunto de blancos que se usarán para calibrar los valores de los pixeles en el software de posprocesamiento; ello permite mejorar el nivel de normalización de las mediciones y capturar imágenes en días de poca claridad [29]. La calibración aplica los modelos señalados en las ecuaciones (1) y (2).
Donde:
Lx: radiancia espectral en la apertura del sensor.
ND: valor crudo del pixel.
Gain y Bias son las constantes de calibración de la cámara.
Donde: (2)
RT: factor de reflectancia de un objetivo desconocido.
θz: ángulo cenital solar en cualquier tiempo t medido con el radiómetro.
DNT(t): matriz digital producida por el UAV cuando el instrumento visa el objetivo.
DNR(t): matriz digital del panel de referencia blanco cuando la cámara está visando al panel de referencia entre los tiempos t0 y t1.
RBaSo4 [0°/θZ]: factor de reflectancia del panel BaSO4 con ángulo 0°y ángulo cenital solar de θZ.
La interacción de la radiación electromagnética con las partículas de agua y polvo que se pueden encontrar entre el sensor de la cámara y el dosel causan perturbaciones en las señales captadas, debido a que actúan como medio absorbente y dispersante de la señal. La corrección atmosférica, que es la eliminación de estas alteraciones, constituye uno de los factores de más influencia al estimar con precisión los índices de vegetación.
Luego de realizar las calibraciones radiométrica y atmosférica se obtienen los valores de reflectancias necesarios para calcular los índices de vegetación NDVI, EVI y SAVI, señalados en las ecuaciones (3), (4) y (5).
Donde:
ρNIR: radiancia en unidades de reflectancia en la banda del infrarrojo cercano.
ρred: radiancia en unidades de reflectancia en la banda del rojo.
Donde:
ρ: reflectancia atmosféricamente corregida en el NIR, regiones espectrales rojo y azul. C1=6,0 y C2=7,5: coeficientes de resistencia atmosférica.
L=1,0: corrección al efecto del fondo del follaje.
G=2,5: factor de ganancia.
Donde:
ρ NIR: reflectancia atmosféricamente corregida en la banda del infrarrojo cercano.
ρ red: reflectancia atmosféricamente corregida en la banda del rojo.
L=0,5: constante de compensación promedio.
El cálculo de estos índices de vegetación mediante sensores instalados en UAV se ha convertido en una alternativa costo-efectiva más flexible y precisa que el procesamiento de imágenes satelitales [12] [30], además de facilitar la realización de un monitoreo periódico del área de estudio [31].
1.6 Clasificación de imagen basada en objetos (object-based image analysis, OBIA)
A través de este modelo se pretende superar los problemas de clasificación basada en pixeles mediante la definición de segmentos, permitiendo que la variabilidad de la reflectancia espectral se use como atributo para discriminar las características de la segmentación. Se realiza esta clasificación OBIA a fin de obtener el reconocimiento de patrones que permitan discriminar las distintas especies vegetales en la zona de estudio.
OBIA se ha usado durante la última década para definir patrones específicos en una imagen. Esta tecnología se orienta hacia la identificación de objetos a través de la información espectral, las formas geométricas, el análisis de vecindad, la textura y el tono de la imagen; y supera la clasificación basada en pixeles debido a que esta última no tiene en cuenta la forma, textura y estructura, clasificadores de baja eficiencia en el reconocimiento de patrones [32].
1.7 Mediciones in situ de variables arbóreas
Una vez hecha la clasificación, se harán mediciones in situ de especies para cuantificar los almacenamientos de carbono en las parcelas seleccionadas a través de la ecuación (6), mediante la medición de variables arbóreas como altura y diámetro a la altura del pecho (DAP). La densidad de la madera (ρ) se tomará de la base de datos reportada por Zanne et al. [33], en la cual se encuentran las especies arbóreas identificadas en la zona de estudio, como se muestra en la tabla 1.
Familia | Nombre científico | Densidad madera (g/c3) masa secada en horno/ volumen fresco | Región |
---|---|---|---|
Fabaceae | Bauhinia variegata | 0,606 | Sureste de Asia (tropical) |
Malvaceae | Heliocarpus americanus | 0,200 | América del Sur (tropical) |
Urticaceae | Cecropia peltata | 0,310 | América del Sur (tropical) |
Juglandaceae | Juglans neotropica | 0,490 | América del Sur (tropical) |
Myrtaceae | Syzygium jambos | 0,700 | América del Sur (tropical) |
Myrtaceae | Psidium guajava | 0,629 | América del Sur (tropical) |
Sapindaceae | Cupania Americana | 0,730 | América del Sur (tropical) |
Fuente: Adaptación de [33].
Para aplicar el modelo alométrico referido en la ecuación 6 se usarán las constantes B1 reportadas en el Protocolo para la estimación nacional y subnacional de biomasa - Carbono en Colombia [34]: allí se establece que para la zona de vida "bosque muy húmedo premontano" (BMH-PM) se puede utilizar la ecuación (6) asociada a ella: B1 = 0,932 y a = -2,289.
Donde:
BA: biomasa aérea en kg.
D: diámetro a la altura del pecho medido a 1,30 m de altura sobre el suelo (cm).
H: altura total del árbol (cm).
ρ: densidad de la madera (g/cm3).
Para las especies como el Otobo (Otobo lehmannii) y Guama (Inga densiflora), de las cuales no se encuentra información de la densidad de la madera en las bases de datos, se usará el promedio del nivel taxonómico superior (familia) [34].
1.8 Obtención de variables biofísicas
Una vez diseñados los transectos, se realizará el recorrido con el instrumento óptico Trac, que registrará la transmisión de luz directa a una frecuencia alta (32 Hz) a través del dosel como insumo fundamental para calcular la fracción de hueco (gap fraction P0(Θ)), lo que definirá las variables biofísicas LAI [35] y Fapar [36] a través de las ecuaciones (7), (8) y (9).
Donde:
Donde:
Existen diversos instrumentos que permiten capturar estas variables en campo, entre los cuales se encuentran el LAI-2200C Plant Canopy Analyzer, el cleptómetro, Trac y cámaras equipadas con lentes de tipo "ojo de pescado" (fisheye) que capturan fotografías hemisféricas. Tanto estas últimas como el LAI-2200C dependen de un cielo uniforme y cubierto ya que, cuando las hojas se encuentran muy iluminadas por el Sol, se producen efectos difíciles de medir. El cleptómetro trabaja bajo condiciones de cielo despejado y utiliza sensores para estimar la luz promedio bajo el dosel [37]; se ha encontrado que estos instrumentos ópticos tienden a subestimar el LAI de bosques donde la distribución de los elementos del follaje no es aleatoria.
Para minimizar la incidencia de estos efectos en la medición de las variables de campo, en la segunda fase de este trabajo se utilizará el mencionado Trac, instrumento óptico desarrollado por [24] y empleado para medir el LAI y la Fapar absorbida por la copa de los árboles. Esta tecnología utiliza la información de la "fracción de hueco" (gap size) del dosel para proporcionar un índice de aglomeración del follaje que cuantifica el efecto de distribución espacial no aleatoria del follaje. Tales elementos conducen a la definición de la arquitectura del dosel a través del modelo de Poisson [38].
1.9 Modelación estadística
Con esta se pretende medir el grado de relación entre las variables biofísicas relacionadas, los índices de vegetación y los datos híbridos de campo. Así mismo, se busca detectar el comportamiento de conglomerados de las variables a través del índice de Moran -ecuación (10)- estadístico, que mide la autocorrelación espacial a través de la identificación de clases de distinta homogeneidad de la muestra; su objetivo es evaluar si el patrón se encuentra agrupado, disperso, o aleatorio.
Donde:
Al medir el grado de relación entre las variables se tendrá la representación gráfica de la correlación a través de diagramas de dispersión, así como la relación que presenta cada una de las variables calculando su correlación cruzada, la correlación de Pearson y la bondad del ajuste a partir del coeficiente de determinación R2. Esto definirá los parámetros necesarios para hallar las constantes del modelo exponencial de regresión, de modo que se pueda hacer la interpolación de la variable dependiente (biomasa) para la totalidad de la imagen.
En este trabajo se aplicó la primera parte de la metodología: esto es, captura y procesamiento de las imágenes, cálculo de los índices de vegetación y OBIA, como se muestra en la figura 1.
2. RESULTADOS
2.1 Vuelo UAV y adquisición de imágenes
Se realizaron tres vuelos para garantizar cobertura y autonomía de vuelo, a partir de los cuales se capturaron 326 fotografías a una altura de 90 m, con traslapes mínimos longitudinal de 80 % y transversal de 60 %. La duración de los vuelos osciló entre 8 y 14 minutos, y se cubrieron cerca de 38 Ha en total (tabla 2 y figura 3).
Parámetro | Itinerario 1 | Itinerario 2 | Itinerario 3 |
---|---|---|---|
Área de cobertura | 14.47 hectáreas | 15.46 hectáreas | 8.2 hectáreas |
Altura de vuelo | 90 m | 90 m | 90 m |
Distancia recorrida | 5,09 km | 4,96 km | 3,14 km |
Velocidad | 7 m/s | 7 m/s | 7 m/s |
Cantidad de baterías | 1 | 1 | 1 |
Duración del recorrido | 13 m, 30s | 13 m, 13s | 8 m, 56 s |
Numero de imágenes | 121 | 130 | 75 |
Espacio necesario en memoria SD | 0,60 GB | 0,65 GB | 0,37 GB |
Resolución | 3,9 cm/ pixel | 3,9 cm/ pixel | 3,9 cm/ pixel |
Fuente: elaboración propia.
Antes de efectuarse los vuelos se realizaron tomas con la cámara NIR al target de calibración suministrado por el fabricante. Estas proporcionan parámetros de reflectancia conocidos y estandarizados para los colores y materiales en los que se encuentra fabricado dicho target, y son necesarios en el proceso de calibración de la ortofoto NIR.
2.2 Procesamiento de restitución para compactación del ortomosaico
Para la construcción de los ortomosaicos se empleó el software Agisoft Photoscan®, aplicación especializada que extrae información tridimensional de las fotografías aéreas, al tiempo que permite la elaboración de ortofotos georreferenciadas mediante la identificación y relacionamiento de características similares (homólogas) entre imágenes en áreas comunes o de traslape [39]. Para la orientación exterior del ortomosaico se utilizaron los puntos de control terrestre pertenecientes al campus de la Universidad del Quindío, de los cuales se estableció su posicionamiento con instrumental topográfico de alta precisión como marco de referencia para trabajos cartográficos y fotogramétricos. Este proceso se ilustra en la figura 3: el color indica el error en Z; el tamaño y la forma de la elipse representan el error en XY; y las posiciones estimadas de los puntos de control se indican con puntos negros.
Sumado a lo anterior, se realizó el modelo digital de elevaciones (MDE) presentado en la figura 4 con la información de la nube densa de puntos obtenida de las imágenes. Este tiene una resolución de 3,89 cm/pixel y una densidad de puntos de 662 puntos/m2. Con la nube de puntos y el MDE se generó el ortomosaico RGB con una resolución de 3 cm/pixel; y se exportó en formatos TIFF y JPG para su posterior procesamiento digital (figura 5).
El ortomosaico NIR, representado en la figura 6, contiene valores de reflectancia en dos longitudes de onda -infrarrojo cercano a 850 nm y luz roja a 660 nm-, los cuales fueron calibrados con la curva de reflectancia de los objetivos de calibración propios del sistema Mapir.
2.3 Estimación de los índices de vegetación
Una vez calibradas las reflectancias de las imágenes, se realizó el cálculo de los índices de vegetación NDVI, EVI y SAVI con el software QGIS 2.18. Los valores obtenidos en el NDVI se presentan en la figura 7. Se pudo observar que permite la diferenciación entre diferentes tipos de coberturas del suelo; no obstante, los valores obtenidos dentro del relicto boscoso son muy similares entre especies, lo que no permite realizar una diferenciación basada en estos índices.
2.4 Clasificación de imagen basada en objetos (OBIA)
El primer paso del proceso de clasificación correspondió a la segmentación de la imagen. En el desarrollo de dicha labor se dio un peso superior a la reflectancia banda infrarroja, con lo cual se buscaba un mayor nivel de diferenciación debido a los atributos espectrales de la vegetación. Para llevar a cabo la clasificación, y con la colaboración del personal del Jardín Botánico, se llevó a cabo el entrenamiento que dio lugar a las siguientes clases: agua, arbustos, asbesto, asfalto, bosque mixto, caracolí, carbonero, casco de buey, caucho, cedro negro, cedro rosado, construcciones, guadual, gualanday, guayabo, mestizo, otobo, pasto, pavimento, suelo desnudo, vehículos y yarumo. Una vez seleccionadas las muestras de cada una de las clases -para lo cual se empleó el software eCognition® 9.0.1-, se realizó la OBIA de la imagen total, como se observa en la figura 8.
3. DISCUSIÓN DE RESULTADOS
En esta fase del trabajo se obtuvieron dos imágenes, una en el espectro visible y otra en el infrarrojo cercano empelando un VANT, con las que se calcularon los índices de vegetación NDVI, EVI y SAVI. Estos no posibilitaron una disgregación entre tipos de vegetación, dada su similitud en coberturas boscosas; sin embargo, mediante la técnica OBIA, y a través del algoritmo multirresolución para definir la segmentación, se dio peso doble a los valores de reflectancia del infrarrojo cercano sobre las reflectancias de las otras bandas, con lo que se logró detectar diferentes especies. De acuerdo con los resultados obtenidos se concluye que, mediante el procesamiento digital de las imágenes obtenidas a través de sensores instalados en VANT, se pueden calcular diferentes índices de vegetación. Así mismo, la alta resolución obtenida (3 cm/pixel) hizo posible realizar una clasificación OBIA de relictos boscosos urbanos, con la que se consiguió una diferenciación ente tipos de vegetación y especies arbóreas; esto último no sería posible si se emplearan imágenes satelitales de resoluciones gruesas, dado que no bastan para discriminar entre especies.
4. CONCLUSIONES
Los resultados obtenidos constituyen un punto de partida importante para caracterizar relictos boscosos mediante el uso de sensores remotos, hecho que propone una metodología novedosa de estimación de variables forestales. El uso de técnicas de observación de la tierra en la estimación de biomasa arbórea disminuye costos en tiempo y dinero sin implicaciones destructivas del ecosistema; sin embargo, las mediciones de campo son necesarias para calibrar los modelos generados. Las variables derivadas NDVI, EVI, SAVI obtuvieron valores con alto grado de similitud y no se logró discriminación de especies a través de estos índices. Como alternativa efectiva se utilizó la clasificación por objetos OBIA, dentro de la cual se asignó peso doble a la banda infrarroja y, en combinación con el algoritmo de segmentación multiresolución, se obtuvo la disgregación de especies o agrupaciones de especies con características similares. Los índices NDVI, EVI, y SAVI derivados a partir de la combinación de bandas multiespectrales, capturadas con sistemas VANT, deberán ser objeto de análisis estadístico para hallar su correlación con las características estructurales del ecosistema.