Área de estudio

El APFFLP comprende un área de 30 500 ha y está ubicada en la región central del estado de Jalisco, al oeste de Guadalajara, la tercera ciudad más importante de México (Figura 1). Ubicada en la sierra del mismo nombre, La Primavera es uno de los relieves volcánicos más diversos de México, donde cúpulas anulares, mesetas, colinas y montañas siguen las líneas de fractura de la caldera volcánica y colinas irregulares modeladas por la erosión; estas formaciones también están influenciadas por las fuerzas fluviales y tectónicas. El pico más alto de La Primavera es el Cerro Planillas (2250 m), seguido por el Cerro San Miguel (2170 m), ambos ubicados en la parte sur del área protegida. Los valles ubicados alrededor son llanuras originadas por deposiciones de espumas aportadas por la formación La Primavera (^{Semarnat, 2000}).

Fuente: elaboración propia.

Figura 1 Localización del Área de Protección de Flora y Fauna La Primavera (APFFLP). 

Según la clasificación Rzedowski, en La Primavera hay cuatro tipos de vegetación: bosque de encino (Quercus), bosque de encino-pino (Quercus-Pinus), bosque de pino-encino (Pinus-Quercus), bosque tropical caducifolio, así como tres comunidades vegetales riparia, rupí y ruderal, que se desarrollan dentro de los diferentes tipos de vegetación antes mencionados. La vegetación representativa es el bosque de encino-pino, el cual está presente desde los 1800 hasta los 2225 msnm, en los puntos más altos del área (cerro San Miguel y cerro Planillas) (^{Semarnat, 2000}).

Procesamiento de imágenes

Se utilizaron imágenes Landsat-7 del path/row 29/46, las cuales son gratuitas y proveen la mejor resolución espacial (30 m) y temporal (16 días) para el monitoreo del APFFLP durante el periodo estudiado. Estos datos fueron generados a partir del Landsat Ecosystem Disturbance Adaptative Processing System (LEDAPS) (^{Vermote et al., 1997}; ^{Masek et al., 2006}) con un nivel de procesamiento de L1T de reflectancia de la superficie y enmascaramiento de nubes a través de FMASK (^{Zhu y Woodcock, 2012}). Estas imágenes se descargaron del repositorio del servicio geológico de EE. UU. (https://earthexplorer.usgs.gov/).

A pesar de su resolución temporal (16 días), de enero de 2003 a diciembre de 2016, la misión Landsat-7 solo generó 238 imágenes, y en este estudio todas estas imágenes fueron usadas. Visto como un cubo de información de tres dimensiones (longitud, latitud, tiempo), nuestro análisis exploratorio reveló la existencia de zonas en donde el porcentaje de datos faltantes por pixel es al menos del 54% (véase la Figura en el Apéndice I). Basándonos en estas imágenes, generamos series de tiempo de NDVI y NBR. Recordamos que:

donde NIR, RED y SWIR refieren a los valores de las bandas de reflectancia del infra rojo cercano, reflectancia del rojo, y reflectancia del infra rojo de onda corta, respectivamente. En el sensor Landsat-7 RED, NIR y SWIR corresponden a las bandas 3, 4, y 7, respectivamente.

El llenado de huecos en estas series de tiempo, a nivel pixel, se atendió por medio de interpolación lineal y splines. Específicamente, en el lenguaje de programación R hemos empleado la rutina na_interpolation del paquete imputeTS de ^{Moritz y Bartz-Beielstein (2017)} con las opciones linear y spline. La opción linear permite aplicar el clásico método de interpolación lineal mientras que la opción spline ajusta un spline cúbico natural, siguiendo la rutina desarrollada en ^{Forsythe (1977)}, los valores de este spline son usados para rellenar los huecos.

Nuestro estudio evalúa el desempeño del método de estimación de área quemada en función de la calidad de los datos disponibles (porcentaje de datos faltante por píxel). Consideramos que esta evaluación es necesaria para aportar evidencia sobre la pertinencia del uso del método propuesto. Por tal motivo, y adicionalmente a las series de tiempo rellenadas vía interpolación lineal y splines, generamos otras series de tiempo de NDVI y NBR a las cuales aplicamos suavizado espacial a cada una de las 238 imágenes del acervo Landsat-7, vía el filtro de paso bajo de ArcMap 10.3.1 con un kernel espacial de 3×3 y el algoritmo gdal_fillnodata de la biblioteca GDAL (^{GDAL, 2019}). Los resultados de aplicar nuestro algoritmo de área quemada a todas estas series de tiempo se discute en la sección Resultados.

Datos de referencia

Para validar el mapa de área quemada obtenido con nuestra metodología empleamos un mapa de área quemada para 2012 hecho con imágenes RapidEye; la alta resolución de una imagen RapidEye (5 m) se obtiene a cambio de un alto costo económico. Es un procedimiento estándar validar un producto de mediana resolución con otro cuya resolución es más alta (véase ^{Justice et al., 2000}). Para la elaboración del mapa de área quemada basado en imágenes RapidEye se aplicó una corrección radiométrica estándar para obtener la reflectancia del techo de la atmósfera (TOA, por sus siglas en inglés), se generaron los índices NDVI y el índice de área quemada propuesto por ^{Martín y Chuvieco (2001)}, y se establecieron umbrales para determinar el área quemada; con base en interpretación visual se eliminaron los píxeles erróneos. Esta interpretación visual es también un procedimiento estándar para identificar áreas quemadas (véase ^{Boschetti et al., 2015}).

Determinación de cambios abruptos en series de tiempo de NDVI

Para estimar estadísticamente cambios abruptos en la tendencia de la serie de tiempo NDVI, se aplicó BFAST de ^{Verbesselt et al. (2010a)}, un método implementado en R. BFAST determina si existen cambios abruptos significativos en la tendencia de una serie de tiempo. Posteriormente, la serie de tiempo se segmenta en tantas partes como cambios abruptos se hayan detectado. BFAST se ha utilizado en múltiples aplicaciones para evaluar diferentes tipos de cambios en la vegetación de varios ecosistemas (^{Watts y Laffan, 2014}; ^{DeVries et al., 2015}; ^{Schultz et al., 2016}; ^{Zewdie et al., 2017} y ^{Wu et al., 2020}, entre muchos otros). En lo que sigue se usa el término puntos de cambio en lugar de cambios abruptos.

Técnicamente, BFAST puede explicarse de la siguiente manera, Sea y _t el valor del NDVI al tiempo y supongamos que la siguiente relación se cumple:

Los tres componentes que forman esta representación son coloquialmente conocidos como tendencia (T _t ), estacionalidad o temporalidad (S _t ) y los errores de medición (ε _t ) los cuales se suponen aleatorios, independientes, con media cero y varianza común σ². La letra N denota el número de observaciones.

BFAST supone que las observaciones pertenecientes al intervalo determinado por cada dos puntos de cambio consecutivos poseen una tendencia lineal. En matemáticas este tipo de tendencia se denomina lineal a pedazos y su representación es la siguiente. Sea K el número total de puntos de cambio en nuestras N observaciones. Usemos t_j para denotar el j-ésimo punto de cambio en la tendencia; j puede tomar valores en el conjunto de números {1, …, K}. Cuando la observación número t pertenece al intervalo (t_j, t_j+1) entonces:

donde α_j y β_j son denominados como el intercepto y la pendiente de esta tendencia lineal. Nótese que tanto el número de puntos de cambio (K), como la ubicación de los puntos de cambio (t_js) como los interceptos (α_js)y las pendientes (β_js) de las tendencias lineales entre puntos de cambio son todos parámetros desconocidos a priori. Para estimar K se minimiza una cantidad conocida como el criterio de información bayesiana (BIC, en inglés). Desde su introducción en estadística en ^{Schwarz (1978)}, el uso de este criterio ha permeado en muchas ramas científicas, (véase ^{Neath y Cavanaugh, 2012}). El BIC, también conocido como criterio de Schwarz, permite seleccionar un modelo dentro de un conjunto finito de opciones (modelos). Este criterio incluye una penalización por el número de parámetros a estimar, lo cual garantiza evitar una sobre estimación.

Una vez estimado K, se deben estimar los K puntos de cambio (t_js) o , o equivalentemente los K-1 segmentos en los que se dividirá la tendencia. Para estimar los t_js se aplica una estrategia de programación dinámica a K-1 problemas de mínimos cuadrados. Concluido este paso es estiman los interceptos (α_js) y las pendientes (β_js) de manera robusta (mediante el método conocido como M -estimation). Al final, obtenemos T^t, el estimador de la tendencia de las observaciones y _t s .

El método de estimación de S _t es similar al descrito arriba para la tendencia. El primer paso consiste en remover la tendencia estimada (T^t) de las observaciones y _t , es decir, se calcula la serie de tiempo yt-T^t. Posteriormente, con base en estas observaciones sin tendencia el componente estacional es estimado mediante una regresión armónica con cambios abruptos; para una representación de esta componente véase la sección 2.1 de ^{Verbesselt et al. (2010b)}. Denotamos por S^t al estimador de la estacionalidad de las observaciones y _t . En BFAST la estimación de T^t y S^t se realiza de forma iterativa y alternada hasta que se alcanza un criterio de convergencia.

A continuación se explica el procedimiento empleado por BFAST para identificar estadísticamente los puntos de cambio en la tendencia de NDVI. Supongamos que en la i-ésima iteración del algoritmo para estimar la tendencia tenemos T^t(t) (estimador de la tendencia en la iteración) y S^t(t) (estimador del componente estacional obtenido en alguna iteración previa por lo cual no requiere un súper índice para propósitos de la siguiente explicación). A partir de estos dos estimadores, el residual al tiempo t (y en la iteración i) es, por definición,

En lo que sigue suprimiremos el súper índice de la notación para simplificar la presentación. Siguiendo la notación dada en ^{Chu et al. (1995)}, dado un ancho de banda h (el cual, por definición, es un número en el intervalo (0,1)) el l-ésimo residual suavizado por una suma móvil (MOSUM, en inglés) es, por definición, igual a:

donde l toma valores en el conjunto de números 0,1,…,N-κ-N-κh, σ^ es un estimador de la varianza de los errores de medición, y N-κh denota el número entero menor o igual a N-κh; por ejemplo, 4.3=4, 10.7=10, etc. Arriba k=2(K-1) y este factor aparece debido a que al estimar la tendencia lineal a pedazos en realidad estamos estimando 2 parámetros (intercepto y pendiente) por cada segmento (en total hay K-1 segmentos). Obsérvese que el número de sumandos en el numerador de r _l,h es siempre igual a N-κh. De modo que, en realidad, r _l,h es un tipo de media móvil con coeficientes constantes e iguales a N-κh1/2/σ^.

El estadístico que BFAST utiliza para determinar la existencia de puntos de cambio en la tendencia está basado en los residuales r _l,h . Específicamente, BFAST emplea el estadístico OLS-MOSUM (OLS, en inglés, son las siglas de mínimos cuadrados ordinarios)

Bajo la hipótesis nula de no existencia de puntos de cambio se sabe que la distribución asintótica de MS _N,h se aproxima a la distribución de la variable aleatoria max_{0≤t≤1/h-1}|W(t+1)-W(t)|, donde W denota el proceso de Wiener (también conocido como el movimiento browniano estándar) (véase el Teorema 2.1 en ^{Chu et al., 1995}, para el resultado de aproximación mencionado arriba y consúltese ^{Karatzas y Shreve, 2012}, entre muchos otros, como referencia de las propiedades del proceso de Wiener). Una consecuencia de este resultado es que el p-valor asociado a esta prueba puede aproximarse utilizando la probabilidad de que |W(t+1)-W(t)|>q_α en algún intervalo 0≤t≤1/h-1; q_α denota el cuantil (1-α)%. Observe que al adaptar los elementos de arriba al caso de S _t , el estadístico OLS-MOSUM sirve también para determinar puntos de cambio en la componente estacional.

La presentación anterior muestra que el estadístico usado por BFAST depende del parámetro ancho de banda h. En consecuencia, la identificación de los puntos de cambio en la tendencia depende del parámetro h. En términos algorítmicos, h controla la porción de residuales que forman parte de la suma móvil. En términos prácticos, este parámetro, básicamente, controla la separación admisible entre dos puntos de cambio estadísticamente significativos. Observe que la separación entre puntos de cambio es relativa al número de observaciones disponibles.

En este trabajo consideramos un acervo de imágenes Landsat-7 durante un periodo de 14 años. De acuerdo con la resolución temporal de este acervo (16 días), por año tendríamos 23 observaciones, y en total cada serie de tiempo de NDVI tendría N=322 observaciones. Con esta N como referencia, identificamos tres regímenes para el valor de h: pequeño, medio y grande. En el régimen pequeño encontramos valores tales como h=23/322(0.075. En este caso, los residuales r _l,h son poco suavizados y tendrán una apariencia similar a los residuales originales (sin suavizar). Este caso es análogo a tener una media móvil con un ancho de banda muy pequeño. Esta característica puede producir que el estadístico MS _N,h rechaze la hipótesis nula cuando esta es verdadera en la población, es decir, se identifiquen más puntos de cambio de los que en realidad existen. En el régimen con h medio consideramos, por ejemplo, h=0.15 y h=0.23. Utilizar un ancho de banda h=0.15, equivale a requerir que existan al menos 48(≈322×0.1) observaciones entre cualesquiera dos cambios abruptos estadísticamente significativos; con base en la resolución temporal del acervo Landsat-7, 48 observaciones equivalen a tener imágenes para cubrir un periodo de 2 años y un mes. Con h=0.23, BFAST puede detectar cambios abruptos siempre que estos tengan una separación de al menos 3 años y 5 semanas (aproximadamente 74 observaciones). En la sección 2.1 de ^{Verbesselt et al. (2010a)} se hace mención a la pertinencia de usar valores de h en el régimen medio para series de tiempo de NDVI con longitud de al menos 10 años. Con esta referencia en mente, y debido a la gran falta de datos en algunas zonas, la evaluación de nuestro método de estimación de área quemada también depende del valor de h y en la sección de Resultados reportamos el análisis con h=0.15 y h=0.23. Para el régimen h grande, se condicionó a BFAST a identificar cambios abruptos siempre que estos pertenezcan a años distintos y estén considerablemente separados entre sí. Por ejemplo, con h=0.45 y N=322, dos puntos de cambio significativos deberán estar separados por al menos 6 años para que BFAST pueda identificarlos apropiadamente. A partir de los registros históricos de incendios en La Primavera sabemos que estos ocurren en lapsos menores a 6 años, por esta razón, en este texto no consideramos valores grandes de h. Nuestro código, disponible en https://github.com/inder-tg/burnSeverity, permite que los usuarios interesados puedan generar mapas de áreas quemadas con BFAST utilizando el valor de apropiado para sus propósitos.

Determinación de área quemada

Se denota con NBR _t el valor del NBR al tiempo t. A continuación suponemos que BFAST ha sido aplicado a una serie de tiempo de NDVI utilizando un ancho de banda h. Denotemos con τ a un punto de cambio en la tendencia de NDVI estimado mediante BFAST. Usamos el símbolo τ para generalizar las ideas presentadas en esta sección. Nótese que τ puede tomar valores en el conjunto de números Nh,Nh+1,…,N-Nh. Obsérvese que τ puede interpretarse como una fecha plausible para la presencia de un cambio abrupto en la tendencia. Dado que la vegetación dominante en La Primavera es el bosque de encino-pino (^{INEGI, 2016}), calculamos:

Note que en la resolución temporal de la misión Landsat-7, la fecha τ-23 corresponde a una fecha con anterioridad de un año a la fecha τ. Similarmente, τ+1 es una fecha justo después del cambio abrupto estimado. Como argumentan ^{Escuin et al. (2008)}, considerando estas dos fechas se minimizan las diferencias específicamente ligadas a cambios fenológicos o condiciones de iluminación. Esta circunstancia, sin embargo, condiciona el valor del ancho de banda h. En efecto, primero observemos que, de acuerdo a la definición de dNBR(τ), necesariamente τ debe ser mayor que 24 para poder comparar los valores de NBR de fechas con cambios abruptos plausibles en NDVI detectadas a partir del segundo año con aquellos valores nominales de NBR del año anterior. Además, por lo mencionado arriba, el mínimo valor que τ puede tomar es Nh. Combinando estas dos condiciones, deducimos que h<24/N es una condición necesaria para que nuestra metodología pueda ligar cambios abruptos en la tendencia del NDVI con afectaciones en la vegetación debidas a quemas.

Habiendo calculado dNBR(τ), empleamos los valores de la Tabla 1, adaptada de ^{Key y Benson (2006)} (véase también ^{Lutes et al., 2006}), para clasificar el tipo de afectación de la vegetación asociada a la fecha τ. Con base en esta tabla, podemos caracterizar áreas con crecimiento de vegetación (rebrote) y áreas quemadas sufriendo algún tipo de severidad. Para la determinación del área quemada consideramos dos categorías únicamente: no quemado (dNBR <0.1) y quemado (dNBR ≥ 0.1). De acuerdo con la Sección 4.2 de ^{Verbesselt et al. (2010a)}, agrupar cambios abruptos anualmente facilita el análisis y la comparación de los resultados. Siguiendo este enfoque, calculamos mapas anualizados de áreas quemadas en La Primavera y determinamos sus niveles de severidad.

Validación del mapa de área quemada del 2012

Usamos el porcentaje de datos faltantes por pixel como indicador de la calidad de los datos empleados para elaborar nuestro mapa de área quemada de 2012; esta calidad de datos ha sido considerada al calcular la exactitud general de estos mapas. Concretamente, se definió un píxel de baja calidad si su porcentaje de datos faltantes oscila entre el 50 y el 53%. De forma similar, se definió un pixel de calidad moderada si este presenta entre un 47% y un 49% de datos faltantes. Con estas definiciones, calculamos la exactitud general de cada mapa mostrado en la Figura 3 de acuerdo a ^{Congalton y Green (2019)}; el mapa de referencia contra el que cada uno de nuestros mapas fue comparado es el mapa de área quemada descrito en la sección Datos de referencia. Los métodos usados para rellenar las series de tiempo que dan lugar a los mapas de la Figura 3 se resumen en la Tabla 3.

Figura 3 Mapas del área quemada en el 2012. Interpolación lineal: (A) h=0.15 (B) h=0.23; Interpolación spline (C) h=0.15 (D) h=0.23; (E) Interpolación con gdal_fillnodata, h=0.15 (F) Interpolación con ArcMap, h=0.23. 

El Mapa A muestra una clara discontinuidad espacial en el área quemada estimada. Hasta cierto punto las características exhibidas en este mapa también están presentes en el Mapa F. Esta aparente discontinuidad es atenuada cuando aumenta el ancho de banda de BFAST (h=0.23) y cuando se usa spline como método de interpolación (Mapas C y D). Los mapas en los que muestran la mayor exactitud general independientemente de la calidad de los datos, lo que h=0.15 quizás se deba al régimen del fuego característico de La Primavera, al presentar incendios anualmente de baja severidad, lo que hace que exista cierta recuperación, que al analizar los cambios con (h=0.23) donde se requiere más de tres años de observación, las condiciones de las áreas quemadas cambien.

Finalmente, la combinación de un enfoque de llenado de huecos espacial y temporal (Mapa E y cuarto renglón en la Tabla 3) produce los niveles de exactitud más altos. Comparado con el enfoque basado únicamente en llenado temporal (Mapas A y B), aparentemente los recursos computacionales invertidos en el llenado de datos espacial y temporalmente propicia un incremento (marginal) en la exactitud global de los mapas.

En promedio, al enfocarse solo en áreas con pixeles de baja calidad, la exactitud general de los mapas de la Figura 3 es de 69%. Esta exactitud general alcanza 91.32% (en promedio) cuando consideramos pixeles de calidad moderada. Esta mejora en la exactitud general sugiere que nuestro método determinará áreas quemadas apropiadamente cuando los datos (series de tiempo de NDVI y NBR) tengan una buena calidad.