1. INTRODUCCIÓN
La contaminación atmosférica se ha convertido en uno de los principales problemas en Colombia debido al deterioro continuo en la calidad del aire y el impacto sobre la salud humana y el medio ambiente [1]. Informes de diferentes instituciones públicas y privadas de vigilancia de la calidad del aire confirman el incremento en las concentraciones de contaminantes ambientales que superan lo permitido en la norma nacional de calidad del aire [2]. La concentración de contaminantes se ve afectada por factores topográficos y meteorológicos como: depresiones por altas montañas que impidan que circule el aire, precipitaciones, formación de neblina, reducción de la radiación solar y alteración de la temperatura entre otros [3]. Este contexto muestra la necesidad de usar herramientas de análisis espacial de contaminación ambiental como apoyo a las decisiones administrativas de Colombia y el mundo.
El análisis de la distribución espacial de la contaminación atmosférica se ha abordado mediante redes de vigilancia [4] y sensores remotos. Las metodologías basadas en sensores remotos incluyen las siguientes etapas: preprocesamiento, en la cual se disminuyen los efectos de ruido y se corrigen radiométrica y geométricamente las imágenes; una etapa de extracción de características, en la que se estiman las variables o índices de interés a partir de las imágenes con el fin de discriminar áreas con características similares [5]-[7]; en la extracción de características se obtienen diferentes índices que son usados, junto a información extra de la calidad el aire [8]-[10] como variables para el análisis de la calidad ambiental urbana. Tal como presenta [11] al estimar la concentración de material particulado PM 10 en zonas con pocas estaciones usando un modelo de regresión empírico del uso del suelo LUR e imágenes de los sensores Landsat 7 ETM+, Landsat 8 OLI y MODIS Aqua-Terra para extraer los índices NDVI, NDSI, SAVI, NDWI y TS como variables de entrada.
Entre los métodos más usados para el análisis de la calidad ambiental se encuentra el Análisis por Componentes Principales PCA. Esta técnica estudia la relación entre variables correlacionadas permitiendo reducir su número a un conjunto de variable sin correlacionar [12]. Tal es el caso de [13], que evalúa la calidad del medio ambiente urbano por medio de un modelo lineal construido a partir de componentes principales extraídos de variables socioeconómicas, datos ecológicos e indicadores ambientales obtenidos de imágenes Landsat TM. Algo similar realiza [14] por medio del análisis de componentes principales al estimar el índice de calidad urbano usando imágenes Landsat ETM+ para el cálculo de los indicadores ambientales. Por su parte [15], mediante el análisis de componentes principales, calcula un índice de calidad ambiental EQI considerando variables agrupadas por dominios (socioeconómico, construcciones, calidad del aire, calidad del agua y uso del suelo) y aplicando PCA sobre cada uno de los grupos de variables.
Este documento presenta una metodología para el análisis de la contaminación ambiental en Medellín a partir de datos de estaciones meteorológicas e imágenes de los satélites Landsat 7 y Landsat 8. El artículo se ordena de la siguiente manera: la sección 2, muestra los materiales y métodos; sección 3, los resultados obtenidos; y, la sección 4, las conclusiones.
2. MATERIALES Y MÉTODOS
La metodología propuesta para el análisis de la calidad ambiental aplicando técnicas de teledetección y análisis de componentes principales se desarrolla en el Software ENVI 5.0 y permite el cálculo de los índices TS, NDVI, TSAVI, NDWI y NSI; la interpolación de variables de calidad del aire PM 10 , PM 2.5 , NO 2 y O 3 ; y el Análisis de Componentes Principales, (verFigura 1). Se utilizan imágenes entre el año 2016 al 2019 obtenidas de la USGS (United States Geological Survey) de los satélites Landsat 7, que contiene 3 bandas del visible (B1 azul, B2 verde y B3 rojo), 2 bandas del infrarrojo cercano (B4 y B5), 1 banda térmica (B6), 1 banda infrarrojo medio (B7) y 1 banda pancromática (B8); Landsat 8, contiene 3 bandas del visible (B2 azul, B3 verde y B4 rojo), 1 banda de aerosol (B1), 1 banda de infrarrojo cercano (B5), 2 bandas SWIR (B6 y B7), 1 banda pancromática (B8), 1 banda cirrus (B9) y 2 bandas TIRS (B10 y B11) y mediciones meteorológicas de las variables de calidad del aire NO 2 , O 3 , PM 2.5 , y PM 10 de las mismas fechas. La información relacionada con estas variables es obtenida del SIATA (Sistema de Alerta Temprana de Medellín y el Valle de Aburrá) (ver Tabla 1).
2.1 Preprocesamiento de las imágenes
Las imágenes requieren de un preprocesamiento para disminuir los efectos del ruido geométrico y radiométrico presente, y ser el insumo en las etapas de análisis visual o digital.
La corrección atmosférica disminuye las distorsiones producidas por la atmosfera que se introducen en los valores de radiancia en las bandas de las imágenes, esta se calcula por medio de (1) [16].
Donde ρ es la reflectividad, E 0 es la irradiancia exoatmosférica, τ 1 es el coeficiente de transmisión atmosférica sol-tierra para cada una de las bandas B1τ 1=0.70, B2τ 1=0.78, B3τ 1=0.85, B4τ 1=0.91, B5τ 1=0.95 y B7τ 1=0.97, τ 2 =1 es el coeficiente de transmisión atmosférica tierra-sensor, L es la radiancia del pixel a corregir, L a es el valor mínimo de niveles digitales por cada banda, θ es el ángulo cenital solar y d=1.49598* 108 Km es la distancia tierra-sol [5].
El preprocesamiento se realiza en el Software ENVI. Por medio del toolbox Radiometric Calibration se convierten los niveles digitales a valores de radiancia. Por otra parte, se realiza la corrección radiométrica del bandeado que corrige el efecto de desplazamiento del histograma de la imagen con el fin de obtener el mismo valor promedio y una misma desviación típica para todas las bandas [17]. Esto se realiza para las imágenes con bandeado del sensor Landsat 7, por medio de la herramienta Landsat_gapfill de ENVI.
2.2 Procesamiento de las imágenes
En esta etapa se usan las imágenes corregidas radiométrica y geométricamente como insumos para estimar distintas variables que son indicadores de calidad del aire.
A continuación, se describe cómo estimar cada una de estas variables.
2.2.1 Temperatura de la superficie TS
La temperatura de la superficie se calcula por medio de (2) [6].
Donde T l es la temperatura de brillo superficial obtenida en función de las constantes de calibración del sensor, λ=0.00115 μm es la longitud de onda media de la banda térmica, ρ=1.438*10-2 mK y ε es la emisividad de la superficie, la cual se calcula teniendo en cuenta la proporción de vegetación [6].
2.2.2 Índice de vegetación normalizado NDVI
Los índices de vegetación son transformaciones matemáticas entre distintas bandas espectrales de la misma escena. Estos permiten evaluar la proporción de reflectancia en la vegetación verde para las longitudes de onda asociadas a la luz roja e infrarroja, pues la clorofila absorbe la luz roja en el proceso de fotosíntesis. El índice de vegetación normalizado NDVI, es uno de los índices más usados en la literatura científica, ya que facilita la interpretación de los parámetros biofísicos de la vegetación fotosintéticamente activa [18], como en (3).
Donde BIRC es la banda del infrarrojo cercano y BR la banda rojo del visible [6].
2.2.3 Índice de vegetación ajustado al suelo transformado TSAVI
El índice TSAVI es un ajuste al NDVI en el que se busca corregir el efecto del brillo del suelo en zonas con escasa cobertura vegetal [19], (4).
En (4) el numerador aparece con la diferencia entre el valor de la banda BIRC y el valor que tendría esta banda si toda la luz en el rojo procediese del suelo, y el denominador se define de tal forma que se corrija los efectos debido a la visión del suelo entre la vegetación, s y a son la pendiente y el intercepto con el suelo [20].
2.3 Métodos de interpolación
La interpolación es usada para predecir valores en una superficie en lugares que no tienen información [23]. Esta es usada para la determinación de las superficies de las variables de calidad del aire NO2, O3, PM2.5, y PM10. Entre los diferentes métodos de interpolación se encuentran los métodos determinísticos y geoestadísticos. La ponderación de distancia inversa IDW es uno de los métodos determinísticos más comúnmente empleados. Por su parte, los interpoladores geoestadísticos cuantifican la estructura espacial de los datos por medio de variogramas. Estos están en el grupo de los kriging [24].
El kriging simple es usado si los fenómenos son estacionarios, con varianza y esperanzas conocidas y constantes [23]. Como en (7):
Donde z(x0) es el valor estimado de la variable en la posición x0, λi son los pesos calculados por medio de la matriz de covarianza, m es la media del valor de la variable, Z(xi) son las mediciones de las variables de interés en los puntos i=1 , 2,…, n [23].
El IDW asigna los pesos a los datos muestreados entorno a una función inversa de la distancia que los separa, por lo cual los puntos más cercanos tienen un peso mayor en el cálculo [25], como en (8):
Donde es el valor estimado para el punto, j.n es el número de puntos utilizados para la interpolación, z j es el valor en el punto i- ésimo, k ij es el peso asociado al dato i en el cálculo del nodo j. Los pesos k varían entre 0 y 1 según el dato y la suma total de los pesos es la unidad [25].
2.4 Análisis de componentes principales
El análisis de componentes principales es una técnica multivariante que estudia la relación que se presenta entre variables correlacionadas y reduce el número de dimensiones mediante la creación de un conjunto de variables sin correlación. Las componentes principales son combinaciones lineales de las variables correlacionadas y se construyen en orden de importancia en cuanto a la variabilidad total de la muestra [26].
El análisis de componentes principales permite el paso de un espacio vectorial Rn a un subespacio Rm (n>m) sin pérdida de información relevante, maximizando su varianza. El PCA calcula una matriz de autocorrelación de los datos con sus vectores propios y los ordena de acuerdo con su valor propio para posteriormente normalizarlos [27]. La componente se calcula mediante un conjunto de variables originales (x 1,x,…,x p ), siendo x un vector con p variables aleatorias [28], como en (9):
Donde ∝' x es la funciona lineal de los elementos de x de máxima varianza, y ∝ 1 es un vector de p constantes y denota la transpuesta ∝ 11 ,∝ 12,…,x 1p . El segundo componente se calcula con ∝' 2 x incorrelacionada con ∝' 1 x . Así se eligen los componentes no correlacionados entre sí, de modo que las variables aleatorias tengan cada vez menor varianza [28]. Por medio de PCA se estima el mapa de calidad ambiental a partir de las variables de calidad del aire NO 2 , O 3 , PM 2.5 , y PM 10 y los índices (TS, NDVI, TSAVI, NDWI y NSI).
3. RESULTADOS Y DISCUSIÓN
Con las imágenes Landsat 7 y 8 se realizó un preprocesamiento que consistió en la corrección de bandeado para el caso del Landsat 7, la corrección atmosférica para todas las imágenes y el recorte para ajustar a la zona de estudio la ciudad de Medellín. Seguidamente, se calculó cada uno de los índices (TS, NDVI, TSAVI, NDWI y NSI) para cada una de las fechas. También, con las mediciones meteorológicas se estimaron las superficies de las variables NO 2 , O 3 , PM 2.5 , y PM 10 usando los métodos de interpolación kriging simple e IDW.
Las superficies asociadas a todas las variables se obtienen con un tamaño de píxel de 30 metros. Después, mediante el análisis de componentes principales se generó un mapa de calidad ambiental asociado a cada fecha estudiada y se identificaron las áreas con calidad de aire más deficiente. Esta información fue cruzada con información socioeconómica descargada del Portal Geográfico del Municipio de Medellín, la cual contenía la localización de los barrios y el estrato socioeconómico.
A continuación, se describe cada superficie encontrada con el fin de identificar los sectores con la temperatura más alta, las zonas con mayor cobertura de vegetación, las áreas con mayor contenido de agua y las zonas con mayor suelo construido.
3.1 Temperatura de la Superficie TS
Se calculó la temperatura de la superficie para las imágenes entre 2016 y 2019. Al igual que en el estudio de [29] las zonas con mayores temperaturas están ubicadas alrededor del corredor vial del río Medellín y la zona centro de la ciudad como se observa en la Figura 2.
En estos sectores se presentan temperaturas por encima de los 30°C y son sectores con construcciones grandes y de baja altura.
3.2 Índice de vegetación NDVI
El índice de vegetación NDVI permitió identificar las zonas con cobertura de vegetación. Se evidenció que en los extremos del casco urbano se presenta una vegetación dispersa. El sector de la ciudad con mayor vegetación dispersa se encuentra en el sureste, siendo un sector residencial de estrato socioeconómico 5 y 6, con una gran cobertura de árboles. Por otro lado, el centro de Medellín es la zona donde se presenta menos cobertura de vegetación al igual que en los barrios aledaños al corredor vial del río Medellín. La Figura 3 muestra el índice NDVI obtenido en las fechas estudiadas.
3.3 Índice de vegetación TSAVI
El índice de vegetación TSAVI puede verse afectado por las variaciones en el brillo del suelo, por lo cual los valores de la cubierta vegetal son más independientes al reflejo del suelo [30]. Este índice permitió determinar las coberturas de vegetación en el casco urbano de Medellín. A diferencia del NDVI, TSAVI identificó menos áreas con vegetación en la zona céntrica de la ciudad, lo cual se asemeja a lo identificado por [29], en su estudio de islas de calor urbana en el Valle De Aburrá. La Figura 4 ilustra el índice calculado para cada una de las fechas de estudio.
3.4 Índice de vegetación NDWI
Este índice determinó las zonas con mayor concentración de humedad, donde al igual que el NDVI y el TSAVI en la zona suroriental y noroccidental, es donde se presenta mayor vegetación y por ende mayor contenido de humedad. En la Figura 5 se presenta el índice obtenido para cada una de las fechas.
3.5 Índice de vegetación NSI
Este índice determinó las áreas construidas con el fin de conocer cómo ha cambiado el uso del suelo y su expansión a lo largo del tiempo y el espacio. Para el año 2019 se presentaron cambios en el uso del suelo (ver Tabla 2), siendo ese año el que mostró la mayor variación. Por otro lado, entre los años 2017 y 2018 no se presentó variación en cuanto al área construida y por lo general las fechas de estudio no presentaron una variación significativa. La Figura 6 ilustra el índice de vegetación obtenido.
3.6 Análisis de relaciones entre los índices
Por medio de matrices de correlación se analizaron las relaciones entre los diferentes índices (NDVI, TSAVI, NDWI, NSI y TS) para cada una de las fechas estudiadas (ver Tabla 3). Se observó que se presenta una correlación alta entre el NDWI, TSAVI y NDVI, debido a que el primero analiza la cantidad de humedad en el suelo mientras el TSAVI y NDVI identifican las zonas con cobertura de vegetación.
Por otro lado, se presentó una correlación negativa o inversa entre el NDVI y el NSI, dado que el NSI busca identificar las zonas con mayores construcciones y el NDVI las coberturas de vegetación.
La temperatura TS presentó una correlación positiva con el NSI, dado que a mayor área construida mayor aumento de temperatura, es por ello por lo que en la zona centro de la ciudad donde se presenta una alta densidad de construcciones, la temperatura es mayor.
Lo contrario ocurre con los índices NDVI, TSAVI y NDWI, ya que la correlación es negativa debido que a mayor temperatura, menor cobertura de vegetación y menor humedad del suelo lo cual concuerda con el estudio de [14], por lo tanto en la zonas con mayor cobertura de vegetación la temperatura es menor, lo cual se ve reflejado en la zona sureste de la ciudad.
El TSAVI y el NDVI presentaron una alta correlación debido a que ambos índices identifican cobertura de vegetación.
3.7 Interpolación de variables de calidad del aire
Las variables de calidad del aire descargadas del SIATA fueron adquiridas de la misma fecha de la toma de las imágenes. La red de monitoreo del SIATA cuenta con 43 estaciones a lo largo del Valle de Aburrá, de las que se tomó la información de 30 (ver Figura 7).
Las variables PM 10 , PM 2.5 , NO 2 y O 3 fueron interpoladas por el método kriging simple para cada una de las fechas de estudio. Este método generó una superficie homogénea que representaba la tendencia de los datos y la distribución de los mismos.
Para las variables PM 2.5 y PM 10 del año 2016, se utilizó el método de interpolación IDW debido a que los datos poseían un distanciamiento grande y el Kriging muestra problemas en la generación de los variogramas, ocasionando que los datos no interpolen o no representen la realidad, por lo cual el IDW se perfiló mejor, para este tipo de distribución de datos.
3.8 Mapa de calidad ambiental
Al igual que en [14], se determinó que las variables ambientales (NDVI, TSAVI, NDWI, NSI y TS) presentaban una correlación alta, positiva o negativa, y por ello eso se decidió usar el PCA. Se calcularon los componentes principales con variables normalizadas entre 0 y 1, dado que la información debe estar en un mismo rango. El primer componente principal contenía el 90 % de la variabilidad en la información y se tomó este como el mapa de calidad de ambiental. Este mapa fue categorizado mediante la desviación típica de cada uno de los datos para cada fecha de estudio. La Figura 8 ilustra los mapas de calidad ambiental para cada una de las fechas.
En el año 2016, 16.47 hectáreas presentaban calidad ambiental muy deficiente y para el año 2017 ese valor aumentó a 40.32 hectáreas, ya que zonas con calidad ambiental deficiente pasaron a la categoría muy deficiente. En el año 2018 el área disminuyó a 36.36 Ha y esta tendencia se mantuvo hasta el año 2019 con 31.32 Ha.
La calidad ambiental buena tendió a aumentar, ya que para el año 2016 era de 1884.24 Ha y para el año 2019 paso a ser 1984.05 Ha, debido a que zonas que presentaban una calidad ambiental moderada pasaron a tener una buena calidad ambiental. Por el contrario, zonas categorizadas con muy buena calidad pasaron a tener calidad buena, como en la zona sureste de la ciudad.
Se evidenció que las zonas con calidad muy deficiente son sectores con una gran cantidad de construcciones y poca cobertura de vegetación, como la zona centro de la ciudad, donde se integran rutas de buses, el metro y se presenta una gran concurrencia de personas.
Las zonas con una calidad del aire muy buena y buena presentan una mayor cobertura de vegetación. Por lo general están ubicadas en los extremos de la ciudad, zonas de estrato socioeconómico 4, 5, y 6, lo cual evidencia que la calidad ambiental está muy asociada a la cobertura de vegetación que es un factor clave de balance y al estrato socioeconómico, debido al área construida y el diseño urbanístico. Los focos de contaminación se ubican en el centro de la ciudad y en algunos sectores aledaños al Aeropuerto Olaya Herrera, debido a su diseño urbanístico con poca cobertura de vegetación, una gran cantidad de construcciones y un gran flujo vehicular.
En este estudio se analizó la calidad ambiental por medio de PCA usando variables ambientales y variables de calidad de aire, a diferencia [12] que usó en su estudio en los Estados Unidos variables de calidad de aire, calidad de agua, índices ambientales, variables de construcción y sociodemográficas. Por otro lado [13] realizó su estudio usando variables socioeconómicas, ecológicas e índices ambientales.
La metodología propuesta para el análisis de la contaminación ambiental se usó para analizar un periodo de 4 años con datos anuales, pero puede adaptarse para utilizar datos con periodicidad trimestral para analizar efectos estacionales. Asimismo, incluyendo variables de calidad del agua, socioeconómicas y sociodemográficas que permitan la creación de índice de calidad ambiental que analice la contaminación ambiental desde varios componentes.
4. CONCLUSIONES
Se propuso una metodología para analizar la distribución espacial de la contaminación ambiental en el casco urbano de Medellín utilizando técnicas de percepción remota y análisis estadísticos. Con la metodología propuesta se obtiene un mapa de calidad ambiental que se enfoca en la contaminación atmosférica dado el uso de variables de calidad del aire como PM 10 , PM 2,5 , NO 2 y O 3 . Este mapa se obtuvo a partir del primer componente, el cual contenía el 90 % de variación de la información y permitió ubicar los sectores con la calidad ambiental más deficiente, siendo insumo para la planificación urbana y para dar atención prioritaria a las áreas más críticas.
De acuerdo con la metodología propuesta, los principales focos de contaminación del casco urbano de Medellín se presentaron en la zona centro de la ciudad y en sectores aledaños al Aeropuerto Olaya Herrera. Estos focos corresponden a zonas con poca cobertura de vegetación, un área construida extensa y elevado flujo vehicular.
El mapa de calidad ambiental puede ser usado como insumo en la oportuna toma de decisiones en cuanto a la planificación urbana por parte de entidades como la Alcaldía y las secretarías de Medio Ambiente y de Planeación, ya que posibilita la pronta intervención en las zonas donde la calidad ambiental es muy deficiente. Proporciona también que las entidades de salud puedan identificar los sectores más vulnerables por la Muy Deficiente calidad ambiental y hacer seguimiento a los residentes de estos sectores, dado los problemas de salud que esto conlleva.
La metodología propuesta puede ser aplicada en diferentes zonas y se recomienda incluir variables socioeconómicas como uso del suelo urbano, densidad de población y variables de calidad de agua, entre otras, debido a que la incorporación de otros indicadores robustece el análisis y permite una mejor comprensión de los factores que influyen en la contaminación ambiental.