2.1. k Vecinos más cercanos(k-NN)

El método de los k vecinos más cercanos (k-Nearest Neighbors k-NN) es uno de los métodos más utilizados para problemas de clasificación, clustering y regresión (^{Ertu\ugrul and Ta\ugluk, 2017}), la regla de los k vecinos más cercanos indica que la clase asignada a un nuevo caso será la clase más votada entre sus k vecinos más próximos del conjunto de entrenamiento, la letra k indica el número de vecinos a utilizar y cuando se utilizan varios vecinos se aprovecha de forma más eficiente la información que se puede extraer del conjunto de entrenamiento (^{Gallego et al., 2018}).

Debido a que los procesos de generación de datos a menudo producen patrones repetidos de comportamiento es posible aplicar el k-NN para regresión (k-NNR), debido a que estos patrones pueden proporcionar información valiosa para predecir el comportamiento de los datos en el futuro (^{LINDSAY and NORMAN, 1977}). El principio de k-NNR consiste estimar la respuesta del punto de prueba t como un promedio ponderado de las respuestas de los puntos k de entrenamiento más cercanos, x₁, x₂,..., x_k, en el vecindario de x_t (^{Hu et al., 2014}).

2.2. Transformada Wavelet

La transformada Wavelet puede ser utilizada para analizar características en tiempo-frecuencia de cualquier tipo de señal o series temporales. A medida que la Wavelet madre (Función de ventana flexible) se mueve a través de la señal durante el proceso de Transformada Wavelet (Wavelet Transform WT), ésta genera varios coeficientes que representan la similitud entre la señal y la Wavelet madre (en cualquier escala específica) (^{Araghi et al., 2015}). El término madre da a entender que las funciones con diferentes regiones de actuación que se usan en el proceso de transformación provienen de una función principal (^{Martínez and Castro, 2002}).

Hay dos tipos de transformada Wavelet: Continua y discreta. El uso de la transformada continua de Wavelet (Continuos Wavelet Transform CWT) puede generar un gran número de coeficientes, haciendo su uso e interpretación más complicado, la transformada discreta de Wavelet (Discrete Wavelet Transform DWT) simplifica el proceso de transformación al tiempo que proporciona un análisis muy efectivo y preciso (^{Partal and Küçük, 2006}).

Para una serie de tiempo discreta x _i , donde x ocurre en un tiempo discreto i, Los coeficientes de DWT pueden ser calculados por la siguiente ecuación (^{Partal and Küçük, 2006}).

Donde m y n son enteros que controlan, respectivamente la dilatación Wavelet (escala) y la translación (tiempo), W _m,n es el coeficiente Wavelet para la Wavelet discreta, 2 m 2 es la ubicación de la DWT.

Aplicando la DWT a una señal, ésta se descompone en dos componentes llamados Coeficientes de Aproximación (Aproximation Coeficiente ApCo) y Coeficientes de Detalle (Detail Coeficients DeCo). ApCo comprende el componente de la señal a gran escala y baja frecuencia, mientras que el DeCo representa el componente de alta frecuencia y pequeña escala. En general, ApCo muestra las características más importantes de la señal, especialmente en el caso del análisis de variaciones a largo plazo y lo que es más importante para estudios de análisis de tendencias. El proceso de descomposición puede utilizarse como un proceso iterativo, en el cual, ApCo de la primera descomposición se desglosa en nuevos ApCo y DeCo (^{Araghi et al., 2015}).

2.3. Regresión por Vectores de Soporte

La idea básica del algoritmo de Regresión por Vectores de Soporte (Support VectorRegresion SVR) es utilizar un mapeo no lineal para llevar los datos a un espacio de características de alta dimensión (^{Smola and Scholkopf, 2004}). Para un conjunto de entrenamiento dado X=((x₁,y₁ ),(x₂,y₂),...,(x_n,y_n)},, donde y _i es el valor objetivo de un valor de entrada x _i , el objetivo de la SVR es encontrar una función f(x _i ) que tenga como máximo una desviación ε de los valores objetivos reales de y _i y al mismo tiempo sea lo más plana posible. Lo que permite aceptar datos que tengan errores menores que ε y rechazar lo que tengan un error mayor (^{Smola and Scholkopf, 2004}). El caso de una función lineal f puede describirse como:

Donde, w es un vector de coeficientes de pesos y b es un termino de ajuste, ambos pueden ser encontrados resolviendo el problema de optimización (^{Smola and Scholkopf, 2004}):

El problema de optimización convexa presentado anteriormente sólo es viable cuando se tiene un f existente que aproxime todos los pares (x _i , y _i ) con ε de precisión. Encontrar un conjunto de datos que permitan esto es realmente difícil, por esto se permiten algunos errores, al agregar variables de holgura, ξ _i , ξ _i ^* haciendo que el problema de optimización sea replanteado como (^{Kazem et al., 2013}):

Donde C > 0 es la constante de costo, cuya función es encontrar un punto de equilibrio entre la complejidad del modelo y el valor más alto de (^{Smola and Scholkopf, 2004}). Usando la formulación dual a través de multiplicadores de Lagrange, la SVR puede extenderse para funciones no lineales, en este sentido el problema de optimización se define sólo en términos del multiplicador de Lagrange α _i , α _i ^*. Esto es posible porque la función Kernel φ(x _i , x _j ) devuelve el producto punto entre los pares en un espacio dimensional de orden superior, sin un mapa explícito de los datos (^{Kazem et al., 2013}).

La función de predicción se puede formular en términos de los multiplicadores de Lagrange y la función del kernel como:

En este artículo, se emplea el Kernel de Base Radial (Radial Basis Kernel RBK), ya que el estado de la técnica muestra buenos resultados a lo largo de una gran variedad de aplicaciones, la función RBK es 𝜑 𝑥 𝑖 , 𝑥 𝑗 = 𝑒 −𝛾 𝑥 𝑖 − 𝑥 𝑗 2 , donde y es un parámetro libre que controla la amplitud de la función.

2.4. Optimización por Enjambre de Partículas

Optimización por Enjambre de Partículas (Particle Swarm Optimization PSO), es una técnica de optimización metaheurística inspirada en el comportamiento colectivo de los animales sociales (^{Marini and Walczak, 2015}). Muchas partículas se distribuyen aleatoriamente en el espacio de búsqueda de dimensión n. El estado actual de cada partícula se describe con posición y velocidad. La partícula actualiza constantemente su posición y velocidad de acuerdo con dos indicadores en cada proceso de iteración. Un indicador es la solución óptima individual de partículas y la otra es la solución óptima global. Esas dos soluciones óptimas se actualizan nuevamente después de cada proceso de iteración. Todas las partículas se agregan a la solución óptima global. Cuando la solución óptima global tiende a ser estable finalmente, el resultado es la solución óptima de la partícula en el espacio variable dimensional N (^{Hu, Dong and Yu, 2016}).

La posición x _i y velocidad v _i de cada partícula es actualizada de acuerdo con las siguientes fórmulas respectivamente:

Donde ϵ ₁ y ϵ ₂ son vectores aleatorios, con valores entre 0 y 1. Los términos α y β son parámetros de aprendizaje, donde ambos pueden tomar valores iguales a 2, finalmente x_i ^*(t) es la mejor partícula actual i y g * ≈ min f(x _i )(∀i=1,2,⋯,n) es el mejor global actual.

2.5. Base de datos

Para este estudio, se utilizaron los datos suministrados por el Sistema de Alerta Temprana de Medellín y el Valle de Aburrá (SIATA). Cada estación realiza mediciones horarias de diferentes variables meteorológicas (velocidad y dirección del viento, temperatura y humedad del aire y radiación solar) y contaminantes del aire (PM₁₀, PM_2,5, NO, NO₂ y O₃). Los datos usados comprenden las mediciones obtenidas entre el 1 de enero de 2013 y el 31 de diciembre del 2016 de dos estaciones de monitoreo ITA-CJUS (Itagüí - Casa de Justicia) y BEL-USBV (Bello - Universidad de San Buenaventura). Se seleccionaron estas dos estaciones, ya que estas realizan mediciones de los contaminantes PM₁₀ y PM₂₅ con menor cantidad de datos faltantes, permitiendo obtener espacios de representación compuestos por 35040 filas para ITA-CJUS y BEL-USBV por 9 columnas para ITA-CJUS y 10 columnas para BEL-USBV.

3.1. Caracterización

Con el fin de obtener un espacio de representación para realizar una predicción con 24 horas de anticipación de la concentración de PM₁₀ y PM_{2 5} utilizando una SVR-PSO (Figura 1), se llevó a cabo una caracterización en tiempo y tiempo-frecuencia de la concentración de los contaminantes y de las variables meteorológicas medidas por el SIATA, la caracterización en tiempo, ha sido con base al trabajo desarrollado previamente en (^{Murillo-Escobar et al., 2019}).

Figura 1 Metodología Propuesta 

Las características temporales utilizadas se dividieron en dos grupos, estos son: promedios, máximos y mínimos de las últimas 24 horas de las variables meteorológicas y de las concentraciones de contaminantes, este grupo de características fue denominado Características Temporales (CT).

Con el objetivo de convertir series temporales de la concentración de PM₁₀ y PM₂₅ en señales estacionarias y regulares con menores fluctuaciones, se utilizó la transformada Wavelet estacionaria, lo que permite conservar la misma cantidad de puntos iniciales. Las Wavelet madres más comúnmente utilizadas son Daubechies (db) y Morlet (^{Chen et al., 2013}; ^{Kalteh, 2015}), para este artículo se utilizó Daubechies 5 con 4 niveles de descomposición. Debido a que este tipo de caracterización requiere un mínimo de muestras, antes de aplicar la transformada Wavelet se realizó una interpolación con la función Spline Cubic para obtener una mayor resolución temporal y así dar más argumentos a la transformada Wavelet

Se utilizaron diferentes escalas temporales para caracterizar la señal de Material Particulado y obtener el mejor espacio de representación, información necesaria para poder predecir la concentración de este en determinado momento del día. En este contexto, se analizaron dos escalas temporales, la primera permite el análisis de la tendencia previa del contaminante durante más tiempo, con el fin de obtener la dinámica subyacente de este, mientras la segunda, pretende obtener información en una menor escala de tiempo, con el fin de detectar la influencia de posibles fenómenos climáticos o sociales.

Escala temporal 1 (E1)

Esta escala temporal reúne información sobre el comportamiento del contaminante durante las últimas 7 semanas en la hora y día que se realizará la estimación, de este modo teniendo como objetivo predecir la concentración de PM₁₀ o PM₂₅ a la hora n (e.g. lunes 14:00), se toma el valor de la concentración del contaminante 30 horas previas de la semana anterior n-30 (e.g. domingo 8:00), además, se obtiene una muestra de las 30 horas previas a la hora objetivo durante las 7 semanas anteriores, posterior a ello, se deriva una serie de tiempo promedio a partir de las 7 series temporales, para luego aplicar la TW con 4 niveles de descomposición. Finalmente se calcula la energía y la media de cada uno de los DeCo y del ultimo ApCo de cada una de las bandas, estas características fueron almacenadas en dos grupos denominados Energías Wavelet (WaEn) y Medias Wavelet (WaMe) (Figura 2).

Figura 2 Caracterización y Escala Temporal 1 

Escala temporal 2 (E2)

Esta escala temporal pretende obtener información más reciente sobre el comportamiento del contaminante, por ello sólo se utilizaron las 30 horas previas de la semana anterior a la hora objetivo en la que se realiza la predicción, de esta escala se obtuvieron de igual manera la energía y media de los coeficientes Wavelet después de 4 niveles de descomposición.

Adicionalmente, se creó un tercer grupo de características, la cual consistió en combinar las variables obtenidas en la Escala 1 y en la Escala 2 (E1+E2).

3.2. Selección de características y escala temporal

Con el fin de determinar la escala temporal y el grupo de características que ofrece las mejores capacidades predictivas del Material Particulado, se realizaron todas las combinaciones posibles de los tres grupos de características propuestos (CT, WaEn y WaMe), dando como resultado un total de cuatro combinaciones.

En la construcción del modelo predictivo de PM₁₀ y PM₂₅ con un horizonte predictivo de 24 horas utilizando SVR-PSO se tomó como conjunto de entrenamiento los datos del año 2013 los cuales fueron divididos aleatoriamente en 70/30 y se realizó la validación en los años 2014, 2015 y 2016. El PSO buscó maximizar el R² y como parámetros se determinaron 100 iteraciones y un total de 20 partículas.

Para determinar la mejor escala temporal se aplica inicialmente el test de normalidad Shapiro Wilk, posterior a ello se aplica el test Friedman para determinar si existe una diferencia estadísticamente significativa al realizar la predicción del Material Particulado al usar cada una de las escalas temporales propuestas. Luego de encontrar en cuál subconjunto de las escalas temporales se logra el mejor desempeño se procede a realizar la misma estrategia, pero con el fin de determinar cuál de las cuatro combinaciones de características presenta el mejor desempeño.

3.3. Imputación

Las bases de datos que se utilizaron contienen una gran cantidad de mediciones ausentes debido a que las diferentes estaciones de monitoreo pueden sufrir daños o ser sometidas a mantenimiento y calibración, desencadenando que sean registrados datos faltantes en las bases de datos, afectando el desempeño del sistema de predicción (^{Chen et al., 2015}), a la vez que limita su aplicación en tiempo real.

Por esta razón se utilizó una imputación de datos faltantes con base a la técnica de los k-NNR (^{Ahmat Zainuri, Aziz Jemain and Muda, 2015}; ^{Zhang et al., 2017}), con el fin de evaluar el efecto de los nuevos valores generados utilizando k-NNR al momento de la predicción, cabe resaltar que sólo se imputaron valores perdidos en las muestras que no superaban 5 variables ausentes, con el fin de no agregar demasiada incertidumbre al sistema. Se evaluó el desempeño usando la mejor escala temporal y el mejor grupo de características para diferentes valores de k (k=3, 6 y 9). Con el fin de encontrar el mejor k para la imputación se utilizó el test de Friedman y un test de comparaciones múltiples. Finalmente, se realiza la comparación entre SVR-PSO con y sin imputación por medio del test de Mann-Whitney.

4.1. Caracterización

Luego de la descomposición Wavelet estacionaria de 4 niveles, se obtuvieron 4 componentes de detalle y de aproximación, de los cuales se generaron 2 grupos de características provenientes de dicha descomposición, la energía de los coeficientes de detalle, el último de aproximación (WaEn) y las medias de las mismas(WaMe), obteniendo así un total de 10 características en la matriz proveniente de la caracterización Wavelet, por otro lado, se obtuvieron las características temporales hora a hora, en donde se agruparon la concentración de contaminantes (PM₁₀, NO, NO₂ y O₃ para BEL-USBV y PM_{2 5}, NO y NO₂ para ITA-CJUS) y las características meteorológicas (temperatura, dirección y velocidad del viento, humedad del aire y radiación solar) de ambas estaciones.

4.2. Selección de características y escala temporal

Para seleccionar el mejor grupo de características, así como la mejor escala temporal, se realizó la predicción de PM₂₅ y PM₁₀ durante el periodo 2014-2016 usando SVR-PSO, donde obtuvo el RMSE por mes.

En la Figura 3 y Figura 4 se pueden observar los errores de predicción de PM10 en BEL-USBV y PM2.5 en ITA-CJUS respectivamente usando SVR-PSO a partir de los 3 grupos de características con las diferentes escalas temporales. El test Shapiro Wilk arrojó un P-Valor < 0,05, indicando que los datos se comportan de manera no normal.

Figura 3 Desempeños predicción SVR-PSO BEL-USBV 

Figura 4 Desempeños predicción SVR-PSO ITA-CJUS 

Posteriormente, el test no paramétrico de Friedman arrojó un P-Valor < 0,05, comprobando que existe una diferencia estadísticamente significativa entre las escalas temporales analizadas E1, E2, E1+E2, tanto en BEL-USBV como ITA-CJUS. El ranking de medias (Tabla 1) del test de Friedman evidencia un mejor desempeño en E2, tanto para BEL-USBV como para ITA-CJUS, además, un test de comparaciones múltiples (Tabla 2) muestra que existe una diferencia estadística entre E2 y las demás escalas para BEL-USBV. En ITA-CJUS no existe una diferencia estadística entre E1 y E2, teniendo en cuenta que el ranking de medias para ITA-CJUS E2 presenta el menor error y para BEL-USBV E2 supera estadísticamente el desempeño de los demás, este fue seleccionado como la mejor alternativa de trabajo.

Una vez determinada la mejor escala temporal, definimos cuál de los 3 grupos de características es el que mejor desempeño brinda, obteniendo un P-Valor $<$ 0,05 en Friedman para BEL-USBV e ITA-CJUS. En la Tabla 3 se observa el ranking de medias y destaca el grupo WaEn+CT con respecto a los demás. Para comprobar que WaEn+CT presenta una diferencia con respecto a los demás se realizó un test de comparaciones múltiples (Tabla 4) en donde se observa que en todas sus comparaciones se obtiene un P-Valor < 0,05 en BEL-USBV. Para ITA-CJUS, sólo se obtuvo significancia estadística respecto a CT+WaMe, no obstante, el ranking de medias en esta estación lo muestra como el mejor grupo de características.

Esta etapa arroja como resultado que la escala temporal E2 usando los grupos de características WaEn+CT ofrece las mejores capacidades predictivas. De la Tabla 3 se puede notar que usar sólo las características Wavelet (WaEn+WaMe) genera el desempeño más bajo, no obstante, al combinar CT con los grupos de características Wavelet da lugar a una disminución significativa en el error, lo cual concuerda con (^{Feng et al., 2015}), en donde obtuvieron los mejores resultados al realizar un modelo que combina características temporales y características provenientes de la transformada Wavelet luego de descomponer la señal en 5 niveles y utilizar una red neuronal para cada señal de descomposición, lo cual supone un método altamente complejo y costoso a nivel computacional, ya que, a pesar de usar redes neuronales, no se define claramente un método de optimización como sí se hizo en este estudio. Además de esto, utilizaron características temporales con una baja resolución, ya que sólo se estimó la concentración promedio diaria, a diferencia del estudio actual que realiza la estimación hora a hora.

4.3. Imputación

Se realizó la predicción con SVR-PSO imputando la base de datos con k-NNR usando diferentes valores de k (k=3, 6 y 9), posteriormente se siguió con la estructura de predicción propuesta anteriormente utilizando el mejor grupo de características. Finalmente, al comparar el desempeño obtenido con los diferentes k (Figura 5) se obtuvo un P-valor < 0,05 en el test de Friedman para BEL-USBV e ITA-CJUS, lo cual indica que el valor de k sí afecta el rendimiento al momento de realizar la predicción, lo cual hace necesario utilizar un test de comparaciones múltiples para identificar el mejor k para cada una de las bases de datos.

Figura 5 Desempeño diferentes valores de k 

En la Tabla 5 se observa que para BEL-USBV el mejor k, en términos de desempeño es k =6, por el contrario, para ITA-CJUS, el mejor es k =9, teniendo en cuenta que no hay diferencia estadísticamente significativa entre k=6 y k=9 para BEL-USBV (Tabla 6), se concluye que k=9 es el más adecuado para imputar ambas bases de datos.

El porcentaje de datos perdidos antes de la imputación era de 8,78% para BEL-USBV y 7,22% para ITA-CJUS, luego de llevar a cabo la imputación con k=9 se redujo a 3,32% y 2,08% respectivamente.

En la Figura 6 se pueden observar los resultados obtenidos luego de la imputación, en donde se evidencia el espacio perdido que hay entre determinadas horas, además se observa que luego de realizar la imputación con k-NNR se logran recuperar estos datos perdidos, generando cierta confiabilidad ya que la figura muestra cómo estos datos siguen la tendencia de la señal y no muestran un comportamiento errático a simple vista, además, se pueden observar las líneas de predicción, las cuales continúan comportándose de manera similar luego de imputar el sistema, convirtiendo la imputación en un elemento de gran importancia porque permite que la predicción siga funcionando ante la presencia de datos perdidos, lo cual es de vital importancia en la implementación real de este tipo de metodologías.

Figura 6 Resultado de la imputación ITA-CJUS 

En la Figura 7 se observa la comparación entre los desempeños obtenidos luego de la predicción SVR-PSO y SVR-PSO-k-NNR. Se aplicó el test Mann-Whitney, el cuál arrojó un P-valor de 0,9416 para BEL-USBV y 0.0990 para ITA-CJUS, lo cual indica que no existe una diferencia estadísticamente significativa luego de realizar la imputación, lo que sugiere que este proceso no afectó el desempeño de la predicción como se constata en la Tabla 7 en donde se muestra que los errores en ambos sistemas SVR-PSO y SVR-PSO-k-NNR se conservan.

Figura 7 Comparación SVR-PSO vs SVR-PSO-k-NNR 

Estos hallazgos son importantes ya que a diferencia de otros estudios, el trabajo actual permito predecir brechas de hasta cuatro a diferencia de otros trabajos como (^{Qin et al., 2014}) y (^{Feng et al., 2015}), Además el manejo valores no tuvo un impacto negativo sobre el desempeño del predictor como en (Shen, ^{Huang and Yan, 2016}).