Repte 3. Predicción de O3 a nivel troposférico mediante modelos de Machine Learning
Informe final

Draft Soriano 860606967-image1.png

Índice

Índice 2

Equipo de Trabajo 3

Glosario (resumido) 3

Descripción del Reto 4

Caja azul. Tabla de indicadores 5

Caja amarilla. Resumen 5

Variables de salida de los modelos (outputs) 6

Variables de entrada de los modelos (inputs) 6

Caja verde. Resumen 7

Descripción del motor de cálculo 7

Entorno de programación 8

Caso de uso sobre zonas de aplicación específicas 9

Opción 1 – ¿Qué está pasando? 10

Opción 2 – ¿Qué pasará? 12

Opción 3 – ¿Qué pasaría si…? 13

Opción 4 – ¿Por qué pasa? 14

Resultados globales 15

Resumen de resultados científico-técnicos 15

Otros resultados científico-técnicos asociados al proyecto 17

Escalabilidad de la solución. Combinación con otros Retos y entidades 18

Requerimientos de síntesis 19

Entorno de programación y librerías asociadas 19

Estructura interna de archivos, espacio de almacenamiento y frecuencia de cálculo 19

Conclusiones y flujo de trabajo futuro. Resumen 21

Anexos 23

Glosario (extendido) 24

Términos relativos a modelos de Machine Learning (ML) y ciencias de la computación 24

Términos relacionados con aspectos de calidad del aire 26

Términos relacionados con aspectos meteorológicos 27

Anexo Caja Amarilla 1. Descripción de variables de entrada y salida 30

Variables de salida de los modelos (outputs) 31

Variables de entrada 33

Anexo Caja Amarilla 2. Cálculo de tipos sinópticos 37

Metodología de cálculo para crear el listado de tipos sinópticos o categorías WT (weather types) 37

Resultado de la clasificación sinóptica diaria: 2000-2020. 40

Anexo Caja Amarilla 3. Variables analizadas y descartadas 43

Conclusiones y flujo de trabajo futuro (extendido) 45

Puntos fuertes 45

Puntos débiles 46

Flujo de trabajo a corto plazo (julio 2022 – diciembre 2022) 47

Flujo de trabajo a medio plazo (enero 2023 – diciembre 2023) 47

Referencias 49

Equipo de Trabajo

  • Equipo de CIMNE

  • Equipo de Universidad de Barcelona
  • Marc Lemus – Postdoc Researcher at Andorra Research + Innovation < mlemus@ari.ad>
  • Equipo de Gencat
  • Eva Pérez – Cap de la Secció d'Immissions. Servei de Vigilància i Control de l'Aire. Generalitat de Catalunya <eva.gabucio@gencat.cat>

Glosario (resumido)1

  • Machine Learning (ML) – Disciplina del campo de la inteligencia artificial que, mediante algoritmos, dota a los ordenadores de la capacidad de identificar patrones en datos masivos y elaborar predicciones (análisis predictivo).
  • Variables de entrada “inputs” – Variables independientes de los modelos de ML, a partir de las cuales se trata de predecir valores de la variable dependiente (variable output).
  • Variable de salida “output” – Variable dependiente, o respuesta del modelo, cuyo valor se trata de predecir a partir de los valores de las variables independientes proporcionados a los modelos de ML.
  • Episodios ambientales – Situación en la que las condiciones atmosféricas son desfavorables para la dispersión y la ventilación, lo que hace que la concentración del contaminante en estudio aumente tanto que puede llegar a superar los valores establecidos por la legislación.
  • Ozono (O3) – Molécula triatómica de oxígeno. En la troposfera se forma mediante reacciones fotoquímicas en presencia de contaminantes precursores como el NOx y compuestos orgánicos volátiles (COVs). Su presencia cerca del suelo es peligrosa porque su elevado potencial de oxidación puede causar daños en los tejidos respiratorios de animales y en los tejidos de las plantas.

(1) En la sección de Anexos se recoge un listado completo de términos asociados al Reto 3

Descripción del Reto

En el reto 3 del proyecto Piksel se ha diseñado una herramienta digital para la visualización, predicción y gestión de la concentración de niveles de ozono (O3) troposférico a escala regional, abarcando todo el territorio de Cataluña. Se pretende que esta herramienta sea multi-propósito, estando orientada principalmente a gestores de la contaminación atmosférica. Para ello, se ha desarrollado un motor de cálculo principal basado en modelos de análisis de datos, utilizando técnicas de Machine Learning (ML). El objetivo de estos modelos es la predicción de niveles de concentración de O3 a partir de la introducción de variables relacionadas con este proceso. Las variables a predecir son dos indicadores de gran relevancia en la gestión de este tipo de contaminante según la legislación vigente, los cuales son utilizados para definir diferentes tipos de umbrales de riesgo para la salud pública. Estos indicadores son: valor máximo diario horario de ozono (O3Max) y valor máximo diario de las medias 8-horarias móviles (O38h).

Para predecir estos valores, los modelos de ML se nutrirán de un gran número de series de datos de diferentes variables. La mayoría se ellas provienen de estaciones de monitorización, tanto de equipos de medición de contaminantes atmosféricos como de redes meteorológicas. Además, se han incluido otras variables calculadas a partir de estadísticos básicos o que definen parámetros temporales. Las fuentes de información de las que se han recopilado los datos son oficiales, de carácter regional o europeo.

Además del citado motor de cálculo, a la hora de elaborar esta herramienta se han planteado otra serie de utilidades que permitan mejorar la gestión de este tipo de contaminante desde una perspectiva más global. Así, los principales objetivos a los que se pretende dar respuesta en el Reto 3 son los siguientes:

  • Visualización de datos históricos de concentración de O3, especialmente orientado al seguimiento de la evolución de este contaminante durante los últimos periodos temporales (última semana, mes, año, etc.), así como al registro de episodios de superación de los límites de concentración marcados por la legislación en esta materia.
  • Predicción de niveles futuros de O3 a 24 y 48 horas, considerando dos de los indicadores más habituales en la gestión de este tipo de contaminante: O3Max y O38h. Estas predicciones se realizarán de forma automática, con una frecuencia diaria y ejecutándose durante la noche, disponiendo el usuario de los valores predichos para el día que acontece (24 horas) y el día siguiente (48 horas).
  • Simulación de niveles de O3 sintéticos, siendo el usuario quien pueda reproducir escenarios de contaminación o meteorológicos ficticios, editando los valores de entrada de los modelos a través de la propia interfaz de usuario.
  • Análisis de la relación entre los variables de entrada de los modelos (inputs) y las variables de salida predichas (outputs), evaluando su grado de importancia en dichas relaciones.

Caja azul. Tabla de indicadores

Se puede considerar que el estado actual de la herramienta es operativo, traduciéndose en un nivel de TRL elevado (nivel orientativo 6-7). El grado de madurez, tanto en relación al desarrollo del motor de cálculo como a su integración en la plataforma, se estima haberse completado en un 90% aproximadamente.

El motor de cálculo trabaja en tiempo real, estableciendo las conexiones con los servidores de adquisición de datos de entrada diariamente y de forma completamente automatizada. Todas las utilidades de la interfaz de usuario se han definido en base a los requerimientos propuestos a lo largo del proyecto por parte de los potenciales usuarios finales de la herramienta (equipo de trabajo de Generalitat de Catalunya). Las opciones de cálculo se encuentran actualmente estables e integradas en la plataforma, a excepción de la opción de “¿Qué pasaría sí…?”. De esta utilidad, tanto el código de programación como la interfaz de usuario han sido finalizados, pero falta establecer la conexión entre ambas para que su funcionamiento sea robusto y eficiente en términos computacionales.

Actualmente, la tecnología de este Reto puede ser testada por parte de los usuarios finales. No obstante, se ha consensuado entre los miembros del equipo de trabajo del Reto 3 validar la herramienta durante el verano. Las pruebas serán realizadas por parte del equipo de la Generalitat de Catalunya, centrándose en dos aspectos fundamentales:

  • Robustez de la herramienta – Se evaluará la herramienta desde el punto de vista operativo, detectando posibles fallos relacionados con la selección de opciones a través de la interfaz de usuario y con los datos mostrados en cada una de las utilidades.
  • Validación de los modelos – Se analizarán los resultados arrojados por los modelos de predicción, comparándolos con las observaciones reales de niveles de concentración de O3 medidos posteriormente en las propias estaciones de calidad de aire.

En base a las conclusiones de esta evaluación, se retroalimentará al equipo de desarrollo del Reto 3 enumerando potenciales mejoras. No obstante, en una primera evaluación general, la valoración de la herramienta es muy positiva y será una herramienta de toma de decisiones para los gestores de calidad del aire de la Generalitat de Catalunya.

Caja amarilla. Resumen

En una serie de Anexos recogidos al final del presente documento (ver serie Anexos Caja Amarilla) se puede encontrar información más detallada sobre las variables analizadas en el este estudio:

  • Anexo Caja Amarilla 1. Descripción de variables de entrada y salida – Descripción detallada de la mayor parte de variables de entrada y salida analizadas y utilizadas en el presente proyecto.
  • Anexo Caja Amarilla 2: Cálculo de variable sinóptica – Metodología de cálculo de variable “tipos de clasificación sinóptica”.
  • Anexo Caja Amarilla 3: Variables descartadas – Listado de variables analizadas, pero finalmente descartadas, junto a la justificación que motivó este descarte.

Variables de salida de los modelos (outputs)

Se ha seleccionado como variables a predecir dos indicadores derivados de la concentración de O3 de especial relevancia en la gestión de este tipo de contaminante según la legislación vigente: O3Max y O38h. Ambos indicadores son utilizados en la definición de umbrales de riesgo asociados a la salud pública.

Variables de entrada de los modelos (inputs)

Como resumen, la Tabla 1 incluye un listado completo de variables de entrada introducidas en un modelo de ML (hay que tener en consideración que no todos los modelos incluirán la totalidad de variables).

Tabla 1. Variables de entrada de los modelos ML para un caso genérico
Datos de calidad de aire COMean, NOMean, NO₂-Mean, NOX-Mean, SO₂-Mean, COMax, NOMax, NO₂-Max, NOX-Max y SO₂-Max.

*Asociado a cada de las variables anteriores se ha calculado los valores de media móvil de días anteriores con una ventana de días de 3 y 7 días

Datos meteorológicos Temperatura media (TM), temperatura máxima (TX), temperatura mínima (TN), Precipitación acumulada 24 horas (PPT24), Humedad relativa (HR), Radiación solar (RS24), velocidad del viento a 10 metros de altura tanto escalar como vectorial (VVEM10 y VVVM10 respectivamente), y dirección de viento con módulo unitario o 1 (DVUM10 y DVM10)

Asociado a cada de las variables anteriores se ha calculado los valores de media móvil de días anteriores con una ventana de días de 3, 7 y 30 días

Otros procesos meteorológicos Estabilidad atmosférica* y clasificación sinóptica (1-10)**.

*Variable introducida como variable categórica (A-D). La clasificación en la que se basa de Pasquill-Gifford se basa en categorías de la A-F. Las dos últimas se han descartado por estar asociadas a periodos nocturnos, momento en el cual no se producen los máximos valores de concentración de O3.

**Variable introducida como variable categórica (1-10), ver Anexo Caja Amarilla 2 – Cálculo de tipos de clasificación sinóptica

Variables temporales Número de días (desde inicio de periodo de entrenamiento), día del año, año, mes y día de la semana.


Caja verde. Resumen

Descripción del motor de cálculo

Para alcanzar los objetivos planteados se ha construido un motor de cálculo principal, basado en técnicas de Machine Learning (ML), con el que poder realizar predicciones de niveles de concentración de O3 a escala regional. Alrededor de este motor se han elaborado una serie de utilidades complementarias que ofrecen opciones de visualización y análisis de resultados. La Figura 1 muestra un esquema general de flujo de funcionamiento del motor de cálculo principal de la “Caja Verde” del Reto 3.

Draft Soriano 860606967-image5.png

Figura 1. Esquema metodológico del motor de cálculo principal del Reto 3

Este motor de cálculo se puede desagregar en cinco bloques principales que se resumen a continuación.

  • Generación de datos de entrenamiento (Bloque 1) – La primera fase consiste en la creación de bases de datos independientes que permitan entrenar cada modelo de ML. Estas bases de datos contienen información histórica de todas las variables intervinientes en el proceso de entrenamiento, tanto inputs como outputs. La primera acción ejecutada en este bloque es la obtención de las propias series de datos históricas, o bien mediante descarga directa a través de APIs o bien calculadas (por ejemplo, series de clasificación sinópticas). A continuación, se lleva a cabo una serie de acciones de preproceso para filtrar las series de cada modelo de ML: ordenación, eliminación de valores nulos, unión de series de variables de calidad con series meteorológicas correspondientes, etc.
  • Construcción de los modelos de ML (Bloque 2) – La segunda fase del motor de cálculo es la construcción de los propios modelos de ML. Estos se han diseñado para predecir el nivel de O3 en los mismos puntos de medición de este contaminante, es decir, en las propias estaciones. En total, se ha trabajado con 47 puntos. Además, se ha elaborado un conjunto de cuatro modelos para cada estación, correspondientes a los cuatro posibles outputs a predecir: O3Max y O38h, tanto a 24 horas como a 48 horas vista. La construcción de los modelos se puede dividir en tres fases principales: entrenamiento en base a las series históricas de 20 años (2000-2020), validación de los modelos junto a acciones e ajuste de metaparámetros y fase de test, donde se evalúa la precisión de los mismos a la hora de generalizar la acción de predicción. Los algoritmos de ML utilizados han sido dos: Random Forests (RF) y Boosted Regression Trees (BRT).
  • Predicción diaria de valores futuros (Bloque 3) – El tercer gran bloque ejecuta las predicciones de valores futuros de niveles de concentración de O3. Para ello, en primer lugar, se crean las nuevas series de variables de entrada con datos descargados en tiempo real. Esta acción se ejecuta automáticamente con una frecuencia diaria, descargando información de diferentes servidores oficiales (a través de sus servicios APIs correspondientes). Posteriormente, se tratan los datos descargados para introducirlos como nuevas series en los modelos ML desarrollados en el anterior bloque. Por último, cada modelo genera los resultados de predicción de nivel de concentración de O3 futura a corto plazo.
  • Creación de resultados derivados (Bloque 4) – En este bloque se engloban todas las acciones orientadas a la creación de nuevos resultados, entendiendo como tales aquellas que derivan o bien de los datos descargados diariamente o bien de los resultados de predicción arrojados por los modelos de ML. Estos resultados derivados se pueden clasificar en cuatro grupos: (i) resultados asociados a mediciones directas de O3, (ii) agregación de resultados de predicción en una escala espacial correspondiente a las Zonas de Calidad de Aire (ZQA), (iii) extrapolación de resultados de forma distribuida en todo el territorio mediante técnicas de interpolación espacial y (iv) resultados relativos a relación entre variables inputs y outputs.
  • Simulación de Escenarios Sintéticos (SES) (Bloque 5) – La última sección del motor de cálculo es el módulo destinado a poder realizar predicciones de niveles de O3, basados en datos de entrada sintéticos. Este bloque, utiliza los mismos modelos de ML utilizados en la tarea de predicción futura (Bloque 3). La diferencia principal es que los valores de entrada podrán editarse por parte del usuario a través de la interfaz, permitiendo realizar simulaciones de escenarios sintéticos.

Entorno de programación

Todo el motor de cálculo ha sido desarrollado ‘ad-hoc’ para el proyecto Piksel, utilizándose dos entornos de programación diferentes, ambos de libre distribución: Python y R. El primero (Python), es donde se ha desarrollado el código principal de la herramienta o core, integrando la mayor parte de procesos y utilidades del Reto. El segundo, R, se ha empleado para elaborar acciones específicas de análisis. Se han establecido las conexiones necesarias entre ambos entornos, así como entre ellos y el servidor de la plataforma del proyecto Piksel. La Figura 2 muestra un resumen de las acciones principales programadas con cada lenguaje. A la derecha de cada proceso, se puede encontrar el bloque al que corresponde, de los cinco descritos anteriormente (‘B1’ corresponde al Bloque 1, ‘B2’ al Bloque 2, etc.)

Draft Soriano 860606967-image6-c.png

Figura 2. Entornos de programación utilizados y procesos desarrollados en cada uno de ellos

Caso de uso sobre zonas de aplicación específicas

En esta sección se describen dos casos de uso específicos de las 47 estaciones de calidad de aire integradas en la plataforma: estación de (Vic) y estación de (Montseny “La Castanya”). Se han seleccionado estas dos estaciones por tratarse de dos puntos de especial interés, al tratarse de dos localizaciones que recogen algunos de los valores más elevados de concentración de O3. En la Figura 3 se puede ver su localización.

La primera (Vic), se encuentra en una región conocida como la Plana de Vic, caracterizada desde el punto de vista orográfico por una depresión alargada en dirección norte-sur. Debido a su particular configuración, en esta zona se han registrado episodios ambientales con cierta frecuencia, por la escasa ventilación atmosférica que se da en esta depresión. La segunda (Montseny “La Castanya”) se localiza en una zona cercana a la anterior, pero fuera de la citada depresión. También suele registrar niveles de O3 por encima de la media del conjunto de estaciones de la red.

Draft Soriano 860606967-image7.png

Figura 3. Localización de las estaciones utilizadas como casos de uso: estación #1 Vic y estación #2 Montseny “La Castanya”

A continuación, se muestran resultados asociados a cada una de las estaciones, recorriendo algunas opciones de visualización y análisis integradas en la interfaz de usuario.2


(2) El lector puede encontrar más detalle sobre las opciones y el funcionamiento de la interfaz de usuario mostrada en esta sección en otro documento de este proyecto: Manual de Usuario Repte 3

Opción 1 – ¿Qué está pasando?

El módulo “¿Qué está pasando?” está enfocada a la visualización de registros históricos de concentración de O3. El usuario puede analizar datos de cualquier periodo de fechas que desee, tanto de la evolución de O3 (con un paso temporal diario) como del número de episodios ambientales que han ocurrido a lo largo de toda la serie (agregados en un paso mensual).

Para ello, basta con seleccionar sobre el visor la estación correspondiente, apareciendo una ventana emergente con las dos opciones de visualización. En la Figura 4 y la Figura 5 se pueden observar los gráficos de evolución de O3 (basados en los indicadores O3Max y O38h) para las dos estaciones de estudio y distintos periodos de estudio.

Draft Soriano 860606967-image8.png

Figura 4. Evolución de O3Max y O38h en la estación de Vic durante el último mes (izquierda) y últimos 3 años (derecha)

Draft Soriano 860606967-image9.png

Figura 5. Evolución de O3Max y O38h en la estación de Montseny “La Castanya” durante el último mes (izquierda) y últimos 3 años (derecha)

Se puede observar cómo ambas estaciones siguen patrones similares. Se detecta un pico asociado a un episodio ambiental el día 17 de junio de este mismo año (información actualizada 28 de junio de 2022). También, ambas gráficas de evolución se caracterizan por presentar un patrón estacional con incrementos de O3 durante los meses de verano, y valores mínimos durante los meses de invierno.

De esta información se puede extraer también el número de episodios ambientales ocurridos en el periodo de estudio seleccionado. En la Figura 6 se puede observar la información relativa a las dos estaciones analizadas. Se puede ver cómo la información de número de episodios ocurridos es ofrecida en base a tres rangos de concentración diferentes: entre 120 y 180 µg/m3 (azul), entre 180 y 240 µg/m3 (naranja) y mayores de 240 µg/m3 (verde).

Draft Soriano 860606967-image10.png
Figura 6. Episodios ambientales ocurridos en los últimos 3 años en las estaciones de Vic (izquierda) y Montseny “La Castanya” (derecha)

Se observa cómo en las dos estaciones se han registrado episodios de alta concentración de ozono troposférico en los últimos 3 años, sobre todo de primera categoría (entre 120 y 180 µg/m3) siendo el año 2019 especialmente significativo. La información recopilada hasta el momento del presente año (2022) muestra un mayor número de episodios en la estación de Montseny “La Castanya” que en la de Vic.

Opción 2 – ¿Qué pasará?

Esta opción permite al usuario visualizar los resultados de predicción futura de niveles de O3 arrojados por cada uno de los modelos de ML. Para generar estos resultados, el motor de cálculo principal de la herramienta se ejecuta de forma diaria durante la noche, pudiendo disponer el usuario de los valores de predicción desde primera hora del día.

Los resultados mostrados se seleccionan a través de los tres menús desplegables de este menú:

  • Variable output – El usuario podrá seleccionar el output del modelo entre O3Max y O38h.
  • Periodo de predicción – Se pueden seleccionar dos periodos de predicción: 24 horas, con la que se calcula la previsión de niveles máximos alcanzada para el día de hoy y 48 horas, con la que se obtienen la previsión de valores máximos para el día de mañana.
  • Tipo de resultados – Combinando las opciones de los dos desplegables anteriores se puede acceder a los resultados de predicción de cada uno de los cuatro tipos de modelos de ML: ‘24 horas – O3Max’, ‘24 horas – O38h’, ‘48 horas – O3Max’ y ‘48 horas – O38h’. Estos resultados, inicialmente de carácter numérico (µg/m3), pueden transformarse en valores categóricos a través de este desplegable.

Como ejemplo, los resultados arrojados por estos modelos de predicción en la fecha de redacción del presente documento (28 de junio de 2022) se recogen en la Tabla 2.

Tabla 2. Valores de predicción de O3Max y O38h(en µg/m3) arrojados por los modelos de ML la noche del día 28 de junio de 2022, con una ventana de predicción de 24 horas (28 de junio de 2022) y 48 horas (29 de junio de 2022).
24 horas (O3Max) 24 horas (O38h) 48 horas (O3Max) 48 horas (O38h)
Vic 141.91 122.63 130.32 109.77
Montseny 119.59 102.68 117.47 102.83


Estos resultados, junto al del resto de estaciones de calidad de Cataluña, pueden ser agregados en dos escalas de resolución espaciales diferentes (ZQA e interpolación espacial). La Figura 7 muestra dos ejemplos de distintos modelos para cada tipo de resolución.

Draft Soriano 860606967-image11.png

Figura 7. Ejemplos de resultados de predicción de diferentes modelos para el conjunto de estaciones de Cataluña; agregados a nivel de ZQA (izquierda) e interpolación espacial (derecha)

Opción 3 – ¿Qué pasaría si…?

Esta opción tiene como objetivo poder simular escenarios sintéticos, editando los valores de las principales variables tipo input desde la propia interfaz de usuario. Los modelos de cálculo a los que se apunta a través de esta opción de interfaz son de la misma naturaleza que en la opción anterior: modelos predictivos de ML cuyos outputs son los indicadores O3Max y O38h.

En la Figura 8 se muestra un ejemplo de la interfaz de usuario desarrollada para utilizar esta opción junto al formato de resultados generados3. Los valores se muestran en sendas casillas, las cuales automáticamente cambian de color de acuerdo a la escala cromática que corresponde para cada caso.

Draft Soriano 860606967-image12.png
Figura 8. Ejemplo de interfaz de usuario de opción de cálculo asociada a opción “¿Qué pasaría si…?

(3) Los resultados de esta opción de cálculo son una muestra sintética, al tratarse de una opción de cálculo aún en desarrollo

Opción 4 – ¿Por qué pasa?

Por último, la sección “¿Por qué pasa?” muestra resultados, generados durante el entrenamiento y validación de los modelos de ML, relacionados con la relación entre variables tipo input y el output de los mismos. Estos resultados aparecen en pantalla en una ventana emergente al seleccionar una estación individual sobre el visor. Se ofrecen dos tipos de resultados:

  • Importancia de variables – Mediante esta opción se pueden visualizar los resultados relativos a la importancia de las variables. La Figura 9 muestra las variables más influyentes para cada una de las estaciones. Se puede observar cómo, en ambos casos, los indicadores basados en mediciones de O3 (lectura de día anterior, medias móviles de varios días anteriores, etc.) son las variables más determinantes a la hora de predecir los outputs de los modelos. A continuación, se encuentran indicadores basados en lectura de radiación solar acumulada (RS24H) y temperatura.
Draft Soriano 860606967-image13.png
Figura 9. Importancia de las variables ordenados de mayor a menor; resultados de estación de Vic y estación de Montseny “La Castanya”.
  • Gráficos de dependencia parcial – Mediante esta opción se puede observar la relación parcial entre cada uno de los inputs seleccionados y el output de los modelos. La Figura 10 ilustra la relación entre dos de las variables más influyentes en cada estación y la variable a predecir.

Draft Soriano 860606967-image14.png

Figura 10. Resultados de dependencia de parcial de dos de las variables más influyentes para cada caso. Estación de Vic (izquierda) variables: media móvil de últimos 7 días de O38h (O3-78h) y media móvil de últimos 7 días de radiación solar acumulada a 24 horas en estación meteorológica “V3” (RS24H7V3). Estación de Montseny “La Castanya” (derecha) variables: media móvil de último día O38h (O3-18h) y media móvil de últimos 30 días de radiación solar acumulada a 24 horas en estación meteorológica “XK” (RS24H30XK).

Resultados globales

Resumen de resultados científico-técnicos

Los resultados tecnológicos alcanzados en este Reto han sido amplios, materializándose todos ellos en forma de herramienta digital integrada en la plataforma Piksel. Una descripción detallada de este tipo de resultados se recoge a lo largo del presente documento y del documento relativo al “Manual de Usuario”.4

Además, de estos resultados de carácter tecnológico, también se ha podido abordar un estudio preliminar desde el punto de científico. Recordemos que el objetivo principal de los modelos de ML es la predicción de valores futuros de los indicadores O3Max y O38h, tanto a 24 como a 48 horas. En el transcurso del proyecto se han analizado cada uno de los 47 puntos asociados a estaciones de monitorización de nivel de O3, distribuidos por toda la región de Cataluña. En esta línea, se ha llevado a cabo un análisis de los resultados arrojados por cada modelo de ML, focalizando el análisis en tres aspectos principales: precisión general de los modelos, distribución de este error en función del valor de la variable a predecir e importancia de las variables.

  • Precisión de los modelos – Se ha medido la precisión global de los modelos de ML utilizando como métrica la raíz del error cuadrático medio (RMSE – por sus siglas en inglés, Root of Mean Squared Error). Este indicador se ha utilizado para evaluar los subconjuntos de datos del periodo histórico (2000-2020) reservados para la fase de test de los modelos de ML. Se han evaluado los 4 conjuntos de modelos (‘24 horas – O3Max’, ‘24 horas – O38h’, ‘48 horas – O3Max’ y ‘48 horas – O38h’), resumiéndose los resultados obtenidos en la Tabla 3. Se puede observar cómo, para los modelos de 24 horas, son algo mejores los valores de precisión de la variable O3Max que la variable O38h. Mientras que, para la predicción a 48 horas, los modelos ofrecen mejores resultados para la variable O38h. No obstante, se pueden considerar valores similares en los 4 conjuntos de modelos (entre 12.66 y 14.12) de valor promedio.
Tabla 3. Precisión de los modelos de ML desarrollados (valores medios de las 47 estaciones de calidad)
24 horas (O3Max) 24 horas (O38h) 48 horas (O3Max) 48 horas (O38h)
RMSE Promedio 12.66 13.60 14.12 13.25
Máximo 15.42 18.83 18.74 15.62
Mínimo 8.11 8.95 9.43 9.15


  • Relación entre error y valor rango de valores de variable a predecir – Los primeros análisis llevados a cabo sobre la distribución de este error en función del rango de valores a predecir, parecen mostrar un empeoramiento en la capacidad predictiva de valores extremales (tanto en rangos muy bajos como muy altos). En la Figura 11 se muestra la correlación existente entre el error residual obtenido en cada predicción y el valor de la variable a predecir.
Draft Soriano 860606967-image15.png
Figura 11. Resultados, para 4 estaciones (08015021, 08019044, 08019050 y 43148028) de correlación entre error residual y valor de observación de la variable a predecir.

Queda patente en estas figuras, cómo existe un comportamiento generalizado en todas las estaciones en el que, efectivamente, en los rangos de valores extremales de O3 el error cometido por los modelos es mayor. Destaca un gran error sobre todo en los valores más altos registrados. Este hecho es coherente con los procesos algorítmicos de ML utilizados, dado que existe falta de representatividad de dichos rangos en las series de datos de entrenamiento. Desde el punto de vista práctico, esto supone un inconveniente cuando el objetivo principal de la herramienta desarrollada se utilice para tratar de predecir valores asociados a episodios ambientales de superación de límites legislativos.

  • Variables más influyentes – Otro de los aspectos que se ha examinado es la importancia de las variables. Este tipo de análisis aporta información relevante sobre la relación entre las variables de entrada y las de salida, indicando el grado de importancia de cada una de ellas. A pesar de existir algunas variaciones entre los valores obtenidos en cada modelo de ML, de forma generalizada se ha obtenido el mismo tipo de resultados para todas las estaciones. Las variables de entrada asociadas a niveles de concentración de O3 (valores máximos y máximos 8-horarios, y medias móviles de 3 y 7 días) han sido las más significativas en la mayoría de casos. Seguidamente, la radiación solar y medias móviles de ésta suelen ser el segundo gran conjunto de variables significativas. Puede existir, para algunas estaciones, cierta importancia de otras variables (temperatura máxima o media, humedad relativa, estabilidad atmosférica, NO2, etc.). En cualquier caso, el grado de importancia de estas últimas, generalmente, representa valores poco significativos.
Draft Soriano 860606967-image16.png
Draft Soriano 860606967-image17.png
Figura 12. Resultados de importancia de las variables asociados a las estaciones de Berga (derecha) y de Montcada i Reixac (izquierda), donde se aprecia el mayor impacto de aquellas variables relacionadas con niveles de O3 y con la radiación solar acumulada a 24 horas (RS24H).

(4) Para más detalle, el lector puede dirigirse a la sección “Caja Verde” y su serie de “Anexos” asociados. Asimismo, se dispone de un “Manual de Usuario” del presente Reto en el que se muestra el alcance tecnológico de la herramienta.

Otros resultados científico-técnicos asociados al proyecto

Además de la propia consecución de los objetivos planteados en el Reto 3, concretados finalmente con la integración del motor de cálculo en la plataforma Piksel, a lo largo del transcurso del proyecto se han alcanzado otros dos hitos significativos:

  • Trabajo Fin de Máster – Durante la primera fase del proyecto (marzo-julio 2021), se elaboró un Trabajo Fin de Máster (TFM) dentro del programa Euroaquae - MSc. in Hydroinformatics and Water Management” por parte de uno de los miembros del equipo de trabajo: Sergio López Chacón. El título del trabajo fue “Prediction of tropospheric ozone concentration at urban locations using machine learning algorithms. Application to Barcelona, Spain” obteniendo una calificación de sobresaliente: 9.
  • Preparación de propuesta para proyecto de Plan Nacional (proyecto MALECÓN) – Durante la segunda fase del proyecto (enero-febrero 2022), se preparó una propuesta para la convocatoria del plan nacional de I+D+i ‘Colaboración Público-Privada’. Esta propuesta, finalmente no fue presentada, encontrándose actualmente a la espera de próximas convocatorias.

Escalabilidad de la solución. Combinación con otros Retos y entidades

El Reto 3 del proyecto Piksel ha sido desarrollado por equipo de trabajo interdisciplinar, tanto en el campo de la investigación fundamental (CIMNE – grupo ‘Machine Learning in Civil Engineering’ y UB – grupo de Climatología’), como en el área de desarrollo tecnológico (CIMNE – grupos ‘TIC – Information and Communication Tecnologies’ y ‘Machine Learning in Civil Engineering’). Otro de los pilares en los que se fundamenta el proyecto es la colaboración en el mismo de los potenciales usuarios de las herramientas desarrolladas (Gencat – ‘Servei de Vigilància i Control de l’Aire’). Así, el proyecto se considera de carácter transversal y constituye una oportunidad única para el intercambio de conocimiento que es esencial para el desarrollo de metodologías útiles e innovadoras.

A nivel tecnológico el grado de madurez alcanzando se considera óptimo, traduciéndose en un nivel de TRL elevado (nivel orientativo 6-7). Esto se traduce en un motor de cálculo e interfaz de usuario que trabaja en tiempo real, de forma completamente automatizada, y con la mayor parte de opciones de cálculo finalizadas y funcionalmente operativas. Esto ha permitido que la tecnología del Reto 3 puede ser testada por parte de los usuarios finales.

Los resultados preliminares arrojados por los modelos de cálculo son prometedores, aunque requieren un mayor análisis y evaluar la incorporación de cambios a corto-medio plazo que mejoren el rendimiento de los modelos. Es este sentido, destaca la conveniencia de llevar a cabo un análisis profundo comparando los valores de predicción futura a corto plazo arrojados por los modelos de ML y las observaciones reales de O3. Fruto de este análisis, se propondrán los siguientes pasos a acometer, pudiendo integrarse nuevos algoritmos o incorporando otras variables de entrada que aún no se hayan contemplado en esta fase del proyecto.

Alineado con este último aspecto (introducción de nuevas variables), se valora positivamente la posibilidad de colaborar con otros grupos de investigación de CIMNE, especializados en el desarrollo de modelos numéricos para la simulación de contaminantes atmosféricos. Destaca en este sentido el trabajo que se ha realizado en el Reto 4 (liderado por el investigador Ignasi de Pouplana), cuyo objetivo principal ha sido la elaboración de un código de cálculo numérico para simular la concentración y dispersión de NO2 troposférico. La evidente interacción entre este componente atmosférico y el O3 (no en vano el primero se trata de un compuesto químico precursor de la formación del segundo), hacen que sea interesante explorar posibles sinergias.

Requerimientos de síntesis

Entorno de programación y librerías asociadas

La arquitectura del Reto 3 se ha diseñado en dos entornos de programación basados en software libre: Python y R. Para poder ejecutar el código desarrollado, en la plataforma Piksel se han instalado una serie de “dockers” o contenedores, compatibles con las versiones de los entornos de programación en los que se ha desarrollado el código del motor de cálculo. Las versiones mínimas recomendables para ambos lenguajes son las siguientes: Python (versión 3.6), y R (versión 4.2). Además, es necesario instalar una serie de librerías para poder ejecutar algunos de los procesos programados. La Tabla 4 recoge las librerías específicas a instalar en cada uno de los dos entornos de programación.

Tabla 4. Librerías utilizadas en el entorno de programación
asdfsdf
Python (versión mínima recomendable 3.6) datetime, glob, joblib, json, math, numpy, os, pandas, pickle, pyproj, requests, rpy2, sodapy, statistics, suntime
R (versión mínima recomendable 4.2) randomForest, forecast, stringr, DescTools 1

sp, rgdal, rgeos, ggmap, ggplot2, latticeExtra, raster 2

tidyverse, fs, lubridate, raster, curl, metR, rgdal, ncdf4, udunits2, PCICt 3

1librerías relacionadas con la “construcción de modelos de ML”

2librerías relacionadas con proceso de “interpolación espacial”

3librerías relacionadas con cálculo de “clasificación sinóptica”


Estructura interna de archivos, espacio de almacenamiento y frecuencia de cálculo

El Reto 3 conlleva la generación y almacenamiento de un gran número de archivos. Existen 47 conjuntos de modelos de ML para cada estación de calidad de aire (cuatro modelos asociados a cada output a predecir), series de datos históricas para entrenar los modelos, datos generados diariamente para realizar nuevas predicciones, resultados de diferente naturaleza, etc. Para organizar toda esta información, se ha creado un conjunto de carpetas/archivos que se describe a continuación, indicando el volumen en disco ocupado por cada subcarpeta:

  • 00_aux_metainformation (10.4 Mb) – Contiene la información necesaria para caracterizar las estaciones (de calidad y meteorológicas), ZQAs, y parámetros internos de cálculo requeridos por el programa.
  • 01_dataframes_training (3.3 Gb) – Aquí se guardan todas las series de datos históricas utilizadas para entrenar cada modelo. Se almacenan series diarias de todas las variables de entrada y salida (diferentes para cada modelo). Además, se incluyen los archivos originales obtenidos con las series históricas en bruto antes de ser preprocesadas. El periodo de las series abarca 21 años (desde 01-01-2000 hasta 31-12-2020).
  • 02_ML_models (1.2 Gb) – Carpeta destinada al almacenamiento de los modelos de ML ya entrenados (tanto los modelos como los scripts programados en R). Estos modelos, guardados en formato binario “.rds”, posteriormente son cargados desde Python para ejecutar los procesos de predicción.
  • 03_dataframes_vector_inputs (2.1 Mb) – Aquí se alojan los archivos que contiene la información de las variables de entrada que se van descargando y procesando diariamente. Estos archivos se van actualizando, añadiendo una serie cada día de ejecución del programa. El incremento de volumen de archivos no supone un problema en términos de espacio.
  • 04_results (10.6 Mb) – Resultados de las predicciones realizadas diariamente y resto de resultados derivados del resto de utilidades. Esta carpeta va actualizándose diariamente, incrementando el volumen de algunos archivos y reescribiéndose otros archivos. El incremento de volumen de archivos no supone un problema en términos de espacio.
  • 05_WT_prediction (1 Mb) – En esta carpeta se encuentra toda la información relativa al cálculo de la variable “clasificación sinóptica”.

La Figura 13 muestra un esquema relacional de la estructura de carpetas.

Draft Soriano 860606967-image18.png
Figura 13. Entornos de programación utilizados y procesos desarrollados en cada uno de ellos

El tiempo de cálculo y la frecuencia con la que se ejecutan los scripts desarrollados son:

  • Procesos que se ejecutan una sola vez (posible actualización periódica anual) – Acciones de preproceso (duración: 20 minutos) y construcción de modelos ML (duración: 2 días).
  • Procesos de ejecución automática diaria – Creación de vectores de predicción y generación de resultados (duración: 5 minutos aprox.).
  • Procesos de cálculo opcional – Simulación de escenarios sintéticos (duración: 30 seg.).

Conclusiones y flujo de trabajo futuro. Resumen

En términos generales, el desarrollo del Reto 3 ha dado como resultado la generación de conocimiento científico relevante para entender y manejar mejor el problema de la contaminación por ozono. También se han alcanzado avances significativos en cuanto a desarrollo e integración de tecnologías. No obstante, se han encontrado puntos débiles que deben ser analizados y mejorados en próximas fases de desarrollo (flujo de trabajo futuro).

Se puede encontrar una versión extendida de esta sección en el anexo correspondiente (ver Anexo Conclusiones y flujo de trabajo futuro).

Puntos fuertes

  • Creación de una amplia base de datos con datos provenientes de redes de monitorización de calidad de aire y redes meteorológicas, a escala regional.
  • Estudio sobre el estado del arte relacionado con procesos complejos asociados al O3 y crear nuevas variables e indicadores.
  • Desarrollo de modelos de ML para la predicción de niveles de concentración de O3 (extrapolable a otros contaminantes).
  • Creación de una herramienta con gran capacidad de automatización y estructura modular que facilita la integración de nuevo modelos.
  • Desarrollo de herramientas que ofrecen distintos servicios útiles para gestores de políticas de calidad de aire.

Puntos débiles

* Falta de información sobre variables relevantes.
* Distribución espacial de la red de monitorización poco homogénea.
* Falta de series de datos sobre predicciones futuras.
* Problemas de descarga a través de APIs: falta de datos puntuales.
* Error medio global de los modelos relativamente bajo.
* Menor precisión en la predicción de valores extremales.

Flujo de trabajo futuro a corto plazo (julio 2022 – diciembre 2022)

  • Buscar otras fuentes de información, que puedan complementar las series de datos utilizadas en esta fase del proyecto.
  • Aplicar otros algoritmos y metodologías de ML, que puedan mejorar los resultados de precisión obtenidos, sobre todo para valores extremales.
  • Realizar tareas de testeo de la herramienta, detectando y corrigiendo posibles errores, y aportando mayor robustez al código de programación.

Flujo de trabajo futuro a medio plazo (enero 2023 – diciembre 2023)

  • Continuar con el desarrollo de la utilidad asociada a la actual opción de cálculo “¿Qué pasaría si?”.
  • Implementar otro enfoque a la hora de construir los modelos de ML, utilizando series de datos de las variables de entrada basadas en predicciones futuras.
  • Valorar la incorporación, como variables de entrada, de resultados generados mediante modelos numéricos.
Repte 3. Predicció de O3 a nivell troposfèric mitjançant models de Machine Learning

Anexos

Glosario (extendido)

Términos relativos a modelos de Machine Learning (ML) y ciencias de la computación

  • Application Programming Interface (API) – Conjunto de definiciones y protocolos para que dos aplicaciones, software o servicios se comuniquen entre sí.
  • Aprendizaje supervisado – Rama del ML, donde se utilizan algoritmos en los que se entrena con series de datos etiquetados. Se entiende como serie de "datos etiquetados" aquella que incluye los valores de las variables independientes y el valor de la variable dependiente a predecir.
  • Bagging – Método algorítmico asociado a técnicas de ML, en el que se selecciona una muestra aleatoria de datos de un conjunto de entrenamiento con sustitución, lo que significa que los datos individuales pueden elegirse más de una vez.
  • Boosting – Método algorítmico asociado a técnicas de ML, mediante el cual los modelos aprenden secuencialmente. Esto significa que se construyen una serie de modelos y con cada iteración del nuevo modelo, se incrementan los pesos de los datos mal clasificados en el modelo anterior.
  • Boosted Regression Trees (BRT) – Son técnicas de ML de aprendizaje supervisado, cuya estructura algorítmica se basa en la generación de un conjunto (ensemble) de árboles de decisión (DT), sobre los que se aplica técnicas de “boosting”.
  • Calibración de modelos – Proceso iterativo de modificación de los meta-parámetros numéricos de los modelos de ML para minimizar la diferencia entre los valores medidos y los valores calculados.
  • Código abierto – Modelo de desarrollo de software basado en la colaboración abierta.
  • Datos de entrenamiento – Conjunto de datos utilizados para la construcción de los modelos de ML. Este conjunto generalmente incluye un amplio número de observaciones, proporcionando valores tanto de las variables independientes como de las dependientes.
  • Datos de validación – Conjunto de datos, independientes de los datos de entrenamiento, utilizados para seleccionar el mejor modelo de ML. Para él, se pueden efectuar acciones como la calibración de los modelos (ajuste de metaparámetros) o la reducción de variables input (seleccionando las variables más influyentes).
  • Datos del test – Conjunto de datos, independientes de los datos de entrenamiento y validación, utilizados para medir la precisión real del modelo de ML finalmente seleccionado. En esta fase, se proporcionan las entradas del conjunto de datos, se realizan las predicciones asociadas a cada serie y se comparan los valores predichos con los reales correspondientes a cada observación. Esto permite cuantificar su capacidad de generalización antes de la introducción de nuevas series de entrada.
  • Decision Trees (DT) o árboles de decisión – Son técnicas de ML de aprendizaje supervisado, cuya estructura algorítmica está formada por nodos internos (representan cada una de las variables de entrada a considerar para tomar una decisión) y ramas (representan la decisión en función de una determinada condición), para llegar a determinar mediante procesos iterativos los valores de una serie de nodos finales (representan el resultado de la decisión).
  • Interpolación espacial de datos – Parte de la geoestadística basada en el cálculo de los valores desconocidos de una variable espacial a partir de otros puntos cuyo valor y localización son conocidos.
  • K-Fold Cross Validation (CV) – Técnica utilizada para evaluar los resultados de los modelos de ML y garantizar que sean independientes de la partición entre datos de entrenamiento y validación. Consiste en dividir el conjunto de datos utilizados para el entrenamiento y la validación en K particiones y realizar un proceso de entrenamiento/validación en K iteraciones, siendo diferente cada subconjunto a cada iteración. Por último, la precisión de los modelos se calcula como la media aritmética obtenida de las medidas de evaluación sobre diferentes particiones.
  • Machine Learning (ML) – Disciplina del campo de la inteligencia artificial que, mediante algoritmos, dota a los ordenadores de la capacidad de identificar patrones en datos masivos y elaborar predicciones (análisis predictivo).
  • Metaparámetros – Valores de los parámetros de un algoritmo de ML utilizados durante el proceso de entrenamiento.
  • Modelos de ML específicos ( ) – Conjunto de modelos predictivos de ML, generados específicamente para el presente proyecto, cuyos outputs difieren en cuanto al período (Period) y al parámetro de predicción (Par.): (i) (output del modelo: O3Max a 24-horas), (ii) (output del modelo: O3Max a 48-horas), (iii) (output del modelo: O38h a 24-horas) y (iv) (output del modelo: O38h a 48-horas).
  • Python - Lenguaje de programación, interpretado y orientado a objetos.
  • R – Lenguaje de programación, interpretado y orientado a objetos.
  • Random Forest (RF) o bosques aleatorios – Son técnicas de ML de aprendizaje supervisado, cuya estructura algorítmica se basa en la generación de un conjunto (ensemble) de árboles de decisión (DT), sobre los que se aplica técnicas de “bagging” para que cada árbol difiera entre sí. El resultado final del modelo se calcula a partir de la media de los resultados obtenidos en cada uno de los árboles generados.
  • Regresión – El análisis de regresión es un subcampo del aprendizaje supervisado que tiene como objetivo establecer un método para la relación entre cierto número de características y una variable objetivo numérica continua.
  • Resolución espacial – Tamaño característico de la discretización utilizada para resolver el dominio de cálculo.
  • Simulador de Escenarios Sintéticos (SES) – Herramienta para simular escenarios sintéticos de O3, mediante la definición de los valores de las variables más relevantes por parte del usuario.
  • Technology Readiness Levels (TRL) – Los niveles de madurez tecnológica son los bloques constitutivos de un método para estimar la madurez de tecnologías.
  • Variables de entrada “inputs” – Variables independientes de los modelos de ML, a partir de las cuales se trata de predecir valores de la variable dependiente (variable output).
  • Variable de salida “output” – Variable dependiente, o respuesta del modelo, cuyo valor se trata de predecir a partir de los valores de las variables independientes proporcionados a los modelos de ML.

Términos relacionados con aspectos de calidad del aire

  • Concentración – Proporción de la cantidad de sustancia disuelta en otra sustancia. Por ejemplo, podemos hablar de la concentración de 150 µg de O3 por cada m3 de aire (150 µg/m3).
  • Compuestos orgánicos volátiles (COVs) – Los compuestos orgánicos volátiles son un tipo de sustancias (hidrocarburos) que se evaporan a la temperatura ambiente. Participan junto con los óxidos de nitrógeno en la formación del O3, siendo por tanto precursores de este contaminante.
  • Dióxido de nitrógeno (NO2) – Compuesto químico formado por dos elementos de oxígeno y uno nitrógeno. Es un gas tóxico, irritante y precursor de la formación de partículas de nitrato.
  • Dióxido de azufre (SO2) – Compuesto químico formado por dos elementos de oxígeno y uno de azufre.
  • Episodios ambientales – Situación en la que las condiciones atmosféricas son desfavorables para la dispersión y la ventilación, lo que hace que la concentración del contaminante en estudio aumente tanto que puede llegar a superar los valores límite establecidos por la legislación.
  • Monóxido de carbono (CO) – Compuesto químico formado por un elemento de oxígeno y uno de carbono.
  • Monóxido de nitrógeno (NO) – Compuesto químico formado por un elemento de oxígeno y uno nitrógeno.
  • Óxidos de nitrógeno (NOx) – Término genérico que hace referencia a un grupo de gases que contienen nitrógeno y oxígeno en diversas proporciones.
  • Ozono (O3) – Molécula triatómica de oxígeno. En la troposfera se forma mediante reacciones fotoquímicas en presencia de contaminantes precursores como el NOx y compuestos orgánicos volátiles (COVs). Su presencia cerca del suelo es peligrosa porque su elevado potencial de oxidación puede causar daños en los tejidos respiratorios de animales y en los tejidos de las plantas.
  • Ozono: Valor máximo diario horario (O3Max) – Variable utilizado para establecer umbrales de información y alerta, asociados a la legislación vigente de contaminación atmosférica.
  • Ozono: Valor máximo diario de las medias 8-horarias móviles (O38h) – Variable utilizada para establecer indicadores asociados a Valor Objetivo de Protección de la Salud humana (VOPS).
  • PM2.5 – Partículas sólidas o líquidas de polvo, cenizas, hollín, partículas metálicas, cemento o polen, dispersas a la atmósfera, y cuyo diámetro es menor a 2,5 µm.
  • PM10 – Partículas sólidas o líquidas de polvo, cenizas, hollín, partículas metálicas, cemento o polen, dispersas a la atmósfera, y cuyo diámetro varía entre 2,5 y 10 µm.
  • Predicción de ozono a 24 horas – Valor de ozono que, siguiendo una metodología determinada y a través de los resultados de los modelos de predicción, van dirigidas a definir el valor más probable de los indicadores de este contaminante en las próximas 24 horas.
  • Predicción de ozono a 48 horas – Valor de ozono que, siguiendo una metodología determinada y a través de los resultados de los modelos de predicción, van dirigidas a definir el valor más probable de los indicadores de este contaminante en un periodo futuro entre 24 y 48 horas.
  • Puntos de medición y equipamiento de la XVPCA – Cada uno de los puntos donde se localiza un medidor de niveles de contaminación atmosférica.
  • Red de Vigilancia y Previsión de la Contaminación Atmosférica (XVPCA) – Sistema de detección de los niveles de inmisión de los principales contaminantes en Cataluña. Red creada por la Ley 22/1983, de 21 de noviembre, definida por la Orden de 20 de junio de 1986 y actualmente está adscrita administrativamente al Departamento de Acción Climática, Alimentación y Agenda Rural de la Generalidad de Cataluña.
  • Umbral de información – Umbral basado en nivel de variable O3Max, establecido en valor de 180 µg/m3.
  • Umbral de alerta – Umbral basado en nivel de variable O3Max, establecido en valor de 240 µg/m3.
  • Valor Objetivo de Protección de la Salud humana (VOPS) – Umbral basado en nivel de variable O38h, que establece que no podrá superarse el valor de 120 μg/m3 en más de 25 ocasiones por año y con una media de 3 años consecutivos.
  • Zonas de Calidad de Aire (ZQA) – Regiones espaciales donde se ha distribuido el territorio de Cataluña, que tienen como objetivo que las medidas que se realizan en una zona sean representativas de la calidad del aire de toda el área que la comprende. Con el objetivo de cada ZQA sea homogénea, los criterios para definir su superficie se basan en aspectos como la orografía, la climatología, la densidad de población y el volumen de emisiones industriales y de tráfico.

Términos relacionados con aspectos meteorológicos

  • Análisis de Componentes Principales (ACP) – Método estadístico que permite simplificar un conjunto de valores de variables mediante combinaciones lineales de nuevas variables, o componentes, no correlacionadas entre sí, sin pérdida significativa de la información. Es útil para analizar grandes conjuntos de datos de muchas variables al expresarlos mediante varios componentes.
  • Clasificación sinóptica – Catálogo de situaciones sinópticas para una región (y un determinado objetivo).
  • Estabilidad atmosférica (EA) – Término relacionado con la mayor o menor facilidad de dispersión (mezcla) que experimenta un paquete de aire en la atmósfera. La estabilidad atmosférica es clave en la dispersión de un contaminante en la atmósfera.
  • Estabilidad atmosférica: Clasificación de Pasquill y Gifford – Esquema utilizado para cuantificar de forma sencilla la estabilidad atmosférica en un punto. Se calcula a partir de los datos de velocidad del viento e irradiación solar y el resultado es una clasificación en 6 categorías: de la A (la más inestable) a la F (la más estable), pasando por la D, que corresponde a una atmósfera neutra (en la que el gradiente real coincide con el adiabática).
  • European Centre for Medium-Range Weather Forecasts (ECMWF) – Centro Europeo de Predicción a Medio Plazo. Organización intergubernamental europea de investigación y servicios operativos de predicción numérica del tiempo.
  • Humedad relativa media (HR) – Humedad relativa media diaria (%).
  • Irradiación solar global (RS24) – Irradiación solar global diaria (MJ/m2).
  • Isobara – Línea que une puntos con la misma presión atmosférica reducida al nivel del mar.
  • Precipitación acumulada (PPT24) – Precipitación acumulada diaria (mm).
  • Presión atmosférica – Peso del aire por unidad de superficie, expresada normalmente en hectopascales (hPa).
  • Reanálisis ERA-5 – Quinta generación del reanálisis atmosférico del ECMWF que cubre el planeta y el período de 1950 hasta la actualidad.
  • Red de Estaciones Meteorológicas Automáticas (XEMA) – Sistema de medición de variables meteorológicas en Cataluña, que gestiona el Servicio Meteorológico de Cataluña (SMC), y que está integrada en la Red de Equipamientos Meteorológicos de la Generalidad de Cataluña (Xemec), creada por la Ley 15/2001, de 14 de noviembre, de meteorología. La XEMA se compone actualmente (diciembre 2017) por un total de 186 Estaciones Meteorológicas Automáticas (EMA).
  • Situación sinóptica – Conjunto típico de configuraciones de isobaras (y, en altura, de isohipsas) de una región.
  • Temperatura media (TM) – Temperatura media diaria obtenida como media de valores de temperatura horarios (ºC).
  • Temperatura máxima (TX) – Temperatura máxima absoluta diaria obtenida como valor máximo de la temperatura horario (ºC).
  • Temperatura mínima (TN) – Temperatura mínima absoluta diaria obtenida como valor mínimo de la temperatura horario (ºC).
  • Viento: Velocidad del viento escalar (VVEM10) – Velocidad media diaria del viento a 10 metros de altura (escalar) (m/s).
  • Viento: Velocidad del viento vectorial (VVVM10) – Velocidad media diaria del viento a 10 metros (vectorial) (m/s).
  • Viento: Dirección del viento m. u. (DVUM10) – Dirección media diaria del viento a 10 metros (módulo unitario).
  • Viento: Dirección del viento m. 1 (DVM10) – Dirección media diaria del viento a 10 metros (módulo 1).

Anexo Caja Amarilla 1. Descripción de variables de entrada y salida

Para la elaboración del Reto 3 ha sido necesario utilizar un gran número de variables con las que poder nutrir o entrenar los modelos predictivos de ML. Inicialmente, es difícil determinar las variables predictoras, o inputs de los modelos, que más influencia tienen en la concentración de O3. No en vano, este es uno de los objetivos del estudio. No obstante, se han seleccionado aquellas que, a priori, se considera que guardan una relación estrecha con este tipo de procesos. La naturaleza de estas variables es heterogénea: contaminantes precursores de la generación de O3 (Ghazali et al. 2010), procesos meteorológicos con influencia en su generación/dispersión o indicadores temporales que puedan aportar información sobre actividades antropogénicas que puedan tener relación. Las variables de salida de los modelos, outputs, son los valores de predicción a 24 y 48 horas de los indicadores O3Max y O38h. La Figura 14 muestra un esquema general de las variables asociadas al motor de cálculo principal con las que se ha trabajado.

Draft Soriano 860606967-image19-c.png
Figura 14. Esquema de “caja roja” de Reto3; relación entre variables de entrada y salida

Las series se datos han sido recopiladas de fuentes de información oficiales, dos de ellas a escala regional (red de monitorización de calidad de aire y red meteorológica) y otra internacional (modelos de situación sinóptica). Se resumen a continuación:

  • Servicio ‘Datos abiertos’ de la Generalitat de Cataluña (Departamento de Medio Ambiente y Sostenibilidad) De esta fuente se han recopilado los datos asociados a variables de contaminación atmosférica, provenientes de equipos de medición automática de la red de ‘Red de Vigilancia y Previsión de Contaminación Atmosférica (XVPCA).5
  • Servicio meteorológico de Cataluña (Meteocat) – De esta fuente se han recopilado los datos asociados a variables meteorológicas, a escala regional. Esta información proviene de equipos de medición automáticos de la ‘Red de Estaciones Meteorológicas Automáticas (XEMA)’.6
  • Centro Europeo de Predicción Meteorológica a Medio Plazo (ECMWF)7 De esta fuente se han recopilado los datos asociados a la variable ‘clasificación sinóptica’ a escala europea. Esta información, proveniente de modelos europeos, ha sido tratada y escalada al ámbito regional por parte del equipo de Climatología de la Universidad de Barcelona (UB). Alternativamente, en el caso de ocurrir algún fallo de conexión con este servidor, también se han establecido mecanismos de conexión con el servidor norteamericano Global Forecast System (GFS), perteneciente a su servicio nacional meteorológico: National Weather Service (NWS)8

(5) https://mediambient.gencat.cat/ca/05_ambits_dactuacio/atmosfera/qualitat_de_laire/vols-saber-que-respires/

(6) https://www.meteo.cat/wpweb/divulgacio/equipaments-meteorologics/estacions-meteorologiques-automatiques/xarxa-destacions-meteorologiques-automatiques-xema/

(7) https://www.ecmwf.int/

(8) https://nomads.ncep.noaa.gov/txt_descriptions/GFS_doc.shtml

Variables de salida de los modelos (outputs)

Se ha seleccionado como variables a predecir dos indicadores derivados de la concentración de O3 de especial relevancia en la gestión de este tipo de contaminante según la legislación vigente: O3Max y O38h. Ambos indicadores son utilizados en la definición de umbrales de riesgo asociados a la salud pública. En el caso del O3Max se han establecido dos límites significativos: umbral de información (180 µg/m3) y umbral de alerta (240 µg/m3). En el caso del O38h el valor límite a considerar es de 120 µg/m3, el cual define el umbral para el indicador de Valor Objectiu de Protecció de la Salut humana (VOPS).

La información asociada al O3 se ha tomado íntegramente de la red XVPCA. Esta red, está compuesta por un total de 105 estaciones (última fecha de actualización 21/01/2022)9, de las cuales 49 realizan mediciones de O3. De este listado de estaciones, dos de ellas fueron descartadas en distintas fases del proyecto, o bien por no disponer de un periodo suficiente de datos para la fase de entrenamiento de los modelos de ML (estación de Vila-Seca) o bien por dejar de proporcionar datos de O3 en el transcurso del proyecto (estación de Rubí). Así, finalmente se dispone de 47 estaciones de medida de O3 (el listado completo se recoge en la Tabla 5). Para cada una de ellas, se ha construido un modelo de ML independiente, por lo que se dispone de 47 conjuntos de modelos de ML, uno por cada una de estas localizaciones. En la Figura 15 se puede observar la distribución espacial de las estaciones de calidad de aire utilizadas.

Draft Soriano 860606967-image20.png
Figura 15. Distribución espacial de las estaciones de calidad de aire utilizadas.
Tabla 5. Listado de estaciones de calidad de aire utilizadas y ZQA asociado
ZQA 1 – Àrea de Barcelona Badalona, Barcelona (Ciutadella), Barcelona (Eixample), Barcelona (Gràcia - Sant Gervasi), Barcelona (Observatori Fabra), Barcelona (Palau Reial), Barcelona (Parc Vall Hebron), El Prat de Llobregat (Sagnier), Gavà, Sant Adrià de Besòs, Sant Vicenç dels Horts (Ribot), Viladecans – Atrium.
ZQA 2 – Vallès – Baix Llobregat Granollers, Montcada i Reixac, Sabadell, Sant Cugat del Vallès, Terrassa
ZQA 3 – Pendès – Garraf Vilafranca del Penedès, Vilanova i la Geltrú
ZQA 4 – Camp de Tarragona Alcover, Constantí, Reus, Tarragona (Parc de la Ciutat)
ZQA 5 – Catalunya Central Berga, Igualada, Manresa
ZQA 6 – Plana de Vic Manlleu, Tona (Zona Esportiva), Vic
ZQA 7 – Maresme Mataró
ZQA 8 – Comarques de Girona Agullana, Montseny (La Castanya), Sant Celoni, Santa Maria de Palautordera, Santa Pau
ZQA 9 – Empordà Begur
ZQA 11 – Pirineu Oriental Bellver de Cerdanya, Pardines
ZQA 12 – Pirineu Occidental Sort
ZQA 13 – Prepirineu Montsec, Ponts
ZQA 14 – Terres de Ponent Juneda (Pla del Molí), Lleida
ZQA 15 – Terres de l’Ebre Amposta, Els Guiamets, Gandesa, La Sénia



(9) https://mediambient.gencat.cat/web/.content/home/ambits_dactuacio/atmosfera/qualitat_de_laire/avaluacio/xarxa_de_vigilancia_i_previsio_de_la_contaminacio_atmosferica__xvpca/Equipament.pdf

Variables de entrada

La precisión de los modelos de ML varía mucho en función del número de observaciones utilizadas y la fiabilidad de los valores, pero también es fundamental que la información introducida sea lo más completa posible y comprenda la mayor parte de fenómenos que intervienen en la variable objetivo a predecir. En ocasiones, esto entraña una gran dificultad debido a la imposibilidad de disponer de datos que monitoricen de forma directa todos los procesos que intervienen. La formación de O3 es un buen ejemplo de este tipo de casos, dado que interviene un elevado número de variables y de mecanismos físico-químicos asociados que en muchas ocasiones son muy difíciles o imposibles de monitorizar. A pesar de que se ha demostrado la gran dependencia de algunas variables en la formación de O3(por ejemplo, la temperatura o la radiación solar) no existe pleno consenso sobre el grado de importancia de todos los mecanismos que intervienen. Esto queda patente en el gran número de trabajos previos existentes cuyo objetivo es similar (estimar o predecir niveles de concentración de O3) pero que en su mayoría se enfocan en aspectos específicos, tratando el problema de forma parcial: influencia de variables meteorológicas, afección de otros componentes de calidad de aire, influencia de algún tipo de actividad antropogénica (por ejemplo, el tráfico o la actividad industrial), etc.

En el proyecto Piksel-Reto3 se han introducido un gran número de variables provenientes de redes de monitorización (otros contaminantes y procesos meteorológicas), así como otras variables derivadas que se considerada puedan tener relevancia a la hora de predecir la concentración de O3.

  • Parámetros de otros contaminantes monitorizados en estaciones de calidad de aire – Los parámetros de calidad de aire seleccionados para cada modelo de ML, han sido recabados de las propias estaciones medidoras de la red XVPCA de las que se han extraídos los valores de O3. No obstante, no en todas las estaciones se mide el mismo número de parámetros, por lo que los inputs difieren en cada modelo. Los contaminantes utilizados han sido los siguientes: monóxido de carbono (CO), monóxido de nitrógeno (NO), dióxido de nitrógeno (NO2), otros óxidos de nitrógeno (NOX), partículas PM10 y partículas PM2,5. Estos, junto a los COVs, son los que se ha considerado que tienen mayor relación con la formación de O3. Precisamente los COVs, se ha descartado al no disponer de casi ninguna estación automática que los monitorice. Otros parámetros que sí miden algunas estaciones, pero se han descartado (por considerarse poco representativo o de menor relevancia en la formación de O3) han sido: benceno (C6H6), ácido sulfhídrico (H2S) y mercurio (Hg).
  • Variables meteorológicas – Se han seleccionado una serie de variables asociadas a fenómenos meteorológicos que pueden influir tanto en la formación como en la dispersión del O3. Las series de datos se han recopilado de la Red de Estaciones Meteorológicas Automáticas (XEMA). Las variables en cuestión han sido: temperatura media (TM), temperatura máxima (TX), temperatura mínima (TN), precipitación acumulada a 24 horas (PPT24), humedad relativa (HR), radiación solar global a 24 horas (RS24), velocidad del viento a 10 metros de altura tanto escalar como vectorial (VVEM10 y VVVM10 respectivamente), y dirección de viento con módulo unitario o 1 (DVUM10 y DVM10). Del total de estaciones meteorológicas que compone la red, para cada modelo de ML se han tomado las series de datos de las estaciones más cercanas geográficamente (ver Figura 16) que miden al menos dos variables: temperatura y radiación solar. Este criterio se ha introducido por considerarse dos de las variables de mayor influencia en la formación de O3. El resto de parámetros puede ser variable en cada modelo, en función de si la estación seleccionada dispone o no de dichas mediciones.
Draft Soriano 860606967-image21.png
Figura 16. Proceso de selección de estación meteorológica más cercana a estación de calidad de aire para la que se construye cada modelo de ML (obligatoriedad de medición de variables T y RS24).
  • Indicadores relacionados con procesos meteorológicos complejos – Existen fenómenos meteorológicos que debido a su gran complejidad son difícilmente monitorizables, y sin embargo pueden afectar directamente a la generación de O3. Un claro ejemplo son los procesos relacionados con la estabilidad atmosférica, los cuales pueden afectar significativamente a la concentración de O3 y a su dispersión, en direcciones tanto horizontal como vertical. Para considerar este tipo de fenómenos se han calculado dos variables meteorológicas además de las anteriormente citadas: estabilidad atmosférica y clasificación sinóptica. El primer término está asociado a la altura de la capa de mezcla y se ha calculado a partir de valores de radiación solar y velocidad de viento, de acuerdo a la clasificación de Pasquill-Gifford (Pasquill 1961). Las clasificaciones sinópticas son agrupaciones de mapas sinópticos que ayudan a entender la variabilidad y complejidad de los patrones meteorológicos en la escala sinóptica. El resultado de la clasificación es un catálogo de los tipos de tiempo que se dan en la zona estudiada a la que se asocian unos patrones de circulación atmosférica. Esta variable ha sido calculada a partir de datos de reanálisis por el grupo de Climatología de la Universidad de Barcelona (UB). Una explicación más detallada sobre su cálculo la podemos encontrar en el anexo siguiente (ver Anexo Caja Amarilla 2 – Cálculo de tipos de clasificación sinóptica).
  • Estadísticos básicos asociados a las variables monitorizadas – A partir de los valores directamente medidos por las estaciones de monitorización (tanto de calidad de aire como meteorológicas), se han generado nuevas variables basadas en parámetros estadísticos sencillos que aporten información adicional a los modelos de ML. Dependiendo de la naturaleza de la variable y del paso temporal con el que se haya medido cada variable, se han considerado estadísticos como: valores medios o máximos diarios, o medias móviles con distinto ancho de ventana (3, 7 y 30 días). Este último tipo de estadísticos puede aportar información relevante si existe una inercia temporal de determinados procesos que pueda afectar a la concentración de O3. Cabe destacar que los valores de medias móviles del propio O3 (3 y 7 días) también han sido incluido como variable de entrada de los modelos.
  • Indicadores relacionados con actividad antropogénica – Por último, la actividad humana es otro factor que debe considerarse en el estudio. Las emisiones asociadas a distintos tipos de actividad, como el tráfico o la industria, pueden afectar de forma directa a los niveles de concentración de ozono. Si bien muchas de estas emisiones están monitorizadas directamente por las estaciones de la red XVPCA, hay consideraciones relevantes asociadas al factor humano a tener en cuenta si se pretende predecir los niveles de O3 cuando aún no se dispone de datos medidos. En este sentido, variables de naturaleza temporal (día del año, mes, día de la semana, festivo/laboral) están asociadas estrechamente con la intensidad de diferentes tipos de actividad antropogénicas. En el presente estudio se han considerado las siguientes variables: número de días desde el inicio del estudio, día del año, día de la semana, mes y año.

Como resumen, la Tabla 6 incluye un listado completo de variables de entrada introducidas en un modelo de ML (hay que tener en consideración que no todos los modelos incluirán la totalidad de variables).

Tabla 6. Variables de entrada de los modelos ML para un caso genérico
Datos de calidad de aire COMean, NOMean, NO₂-Mean, NOX-Mean, SO₂-Mean, COMax, NOMax, NO₂-Max, NOX-Max y SO₂-Max.

*Asociado a cada de las variables anteriores se ha calculado los valores de media móvil de días anteriores con una ventana de días de 3 y 7 días

Datos meteorológicos Temperatura media (TM), temperatura máxima (TX), temperatura mínima (TN), Precipitación acumulada 24 horas (PPT24), Humedad relativa (HR), Radiación solar (RS24), velocidad del viento a 10 metros de altura tanto escalar como vectorial (VVEM10 y VVVM10 respectivamente), y dirección de viento con módulo unitario o 1 (DVUM10 y DVM10)

Asociado a cada de las variables anteriores se ha calculado los valores de media móvil de días anteriores con una ventana de días de 3, 7 y 30 días

Otros procesos meteorológicos Estabilidad atmosférica* y clasificación sinóptica (1-10)**.

*Variable introducida como variable categórica (A-D). La clasificación en la que se basa de Pasquill-Gifford se basa en categorías de la A-F. Las dos últimas se han descartado por estar asociadas a periodos nocturnos, momento en el cual no se producen los máximos valores de concentración de O3.

**Variable introducida como variable categórica (1-10), ver Anexo Caja Amarilla 2 – Cálculo de tipos de clasificación sinóptica

Variables temporales Número de días (desde inicio de periodo de entrenamiento), día del año, año, mes y día de la semana.


Anexo Caja Amarilla 2. Cálculo de tipos sinópticos

Los mapas del tiempo de superficie, que muestran los valores de la presión atmosférica en superficie mediante el trazado de isobaras, expresan la dinámica atmosférica para un momento determinado, determinando la dirección y la velocidad del viento, las áreas de estabilidad e inestabilidad, etc. El reconocimiento en los mapas del tiempo de las formas de isobaras típicas (anticiclones, borrascas, dorsales, vaguadas, etc.) es el primer paso para deducir el tipo de tiempo, es decir, la temperatura, la humedad del aire, la precipitación, el viento, etc., para cada área. Desde un punto de vista climático, tiene interés definir las situaciones sinópticas, o patrones sinópticos, que son los conjuntos típicos de configuraciones de isobaras sobre una determinada región. Su determinación resulta del máximo interés, porque a ellas van asociadas unos tipos de tiempo característicos. Una clasificación sinóptica es un catálogo de situaciones sinópticas para una determinada región compuesto por las situaciones más frecuentes, que expresan la variabilidad atmosférica.

Los mapas del tiempo previstos a 24, 48, 72 h, u otros intervalos, construidos hoy por los modelos numéricos de predicción meteorológica, constituyen la mejor herramienta para la predicción del tiempo a una escala sinóptica (de unos miles de kilómetros) y de mayor detalle, para esos plazos temporales.

Metodología de cálculo para crear el listado de tipos sinópticos o categorías WT (weather types)

Los datos meteorológicos usados para el cómputo de la clasificación sinóptica provienen del modelo de reanálisis atmosférico ERA-5 (Hersbach et al. 2020), a escala diaria, y a una resolución espacial nativa de 0.25º, remuestreada a 1º, para el período 2000-2020. El dominio geográfico usado es la ventana 20ºW-20ºE, 30ºN–60ºN. La variable usada ha sido la presión atmosférica reducida a nivel del mar (SLP, del inglés sea level pressure).

El primer paso para la obtención de una clasificación de tipos de circulación atmosférica a partir de la aproximación estadística multivariante del Análisis de Componentes Principales (ACP) es la creación de una matriz que recoja los datos atmosféricos a clasificar, en nuestro caso obtenidos a partir de datos de reanálisis.

Es el caso de la matriz en modo S, con una disposición tiempo (t) x espacio (g) en las columnas y filas, respectivamente. Una representación esquemática de esta matriz y su disposición se presenta en la Figura 17.

Draft Soriano 860606967-image22-c.png
g1 (lon1,lat1) g2 (lon2,lat2) g3 (lon3,lat3) gn (lonn,latn)
t1 SLP(t1,g1) SLP(t1,g2) SLP(t1,g3) SLP(t1,gn)
t2 SLP(t2,g1) SLP(t2,g2) SLP(t2,g3) SLP(t2,gn)
t3 SLP(t3,g1) SLP(t3,g2) SLP(t3,g3) SLP(t3,gn)
tn



Figura 17. En la izquierda representación esquemática de la disposición de un conjunto multidimensional de datos climáticos a la hora de aplicar un ACP en modo S. Las matrices de color verde representan la presión atmosférica en superficie (SLP). La variable tiene n pasos temporales (tn) y puntos de retícula (gn). A la derecha encontramos la disposición de la matriz de datos utilizada a la hora de realizar este ACP en modo S, que se concreta en que las variables, o cabezas de columna, son los puntos de retícula, y las filas son las observaciones de todos los casos de SLP.

Al aplicar el análisis de componentes principales (ACP) sobre esta matriz de datos original (t)x(g), se obtiene un conjunto de componentes principales (CPs) que básicamente son la combinación lineal –es decir, una suma ponderada– de todos los pares de coordenadas de la retícula de valores de SLP. Los pesos en esta suma ponderada son correlaciones de Pearson –a las que llamaremos loadings, a partir de ahora–, que se emplean a la hora de ajustar cada uno de los ejes de estas CPs, con el fin de obtener el mejor resumen en los datos originales.

Hay que tener en cuenta que estos ejes derivados de la aplicación del ACP son ortogonales entre sí. Debido a la ortogonalidad de estos componentes, la segunda CP será independiente de la primera, mostrando una correlación 0 entre ambas CPs, a la vez que captando la máxima varianza posible. La tercera CP también será ortogonal a las anteriores, por tanto, no estará correlacionada con las dos CPs anteriores. Este proceso sigue hasta que se han calculado el total de n componentes principales, es decir, el número de componentes principales es el mismo que el número original de variables, que en este caso son los puntos de retícula de SLP, pero optimizando la distribución de la varianza explicada del total de los datos.

Sin embargo, y como se ha comentado anteriormente, el ACP supone la rotación del conjunto original de datos, lo que implica que estas observaciones se encuentran en posiciones distintas a las originales por el hecho de quedar explicadas por el nuevo sistema de coordenadas que representan las CPs. Estas coordenadas derivadas de la aplicación del ACP se denominan scores y son calculados como combinaciones lineales de las variables originales (serie temporal de campos de SLP) y los pesos (loadings). Dicho de otra forma, estas nuevas coordenadas multivariantes permiten, precisamente, conocer el grado de representatividad de cada uno de nuestros casos originales en cada una de las componentes principales, a la vez que permiten buscar grupos dentro de estos casos por proximidad (estadística) en el espacio multivariante.

Para lograr, pues, la deseada clasificación sinóptica será necesario trabajar con los scores, por los motivos ya mencionados anteriormente. Sin embargo, previamente es necesario decidir qué número de CPs a retener para seguir con el tratamiento posterior de los resultados manteniendo el máximo de información (varianza). En esta investigación se ha optado por hacer uso del Scree test (Catell, 1966), que se basa en el cambio de pendiente en la línea que une los puntos que representan porcentajes de varianza de cada componente principal. El punto previo a la caída de la varianza se convierte en el punto de corte (Figura 18). En este caso, se decidió usar los 5 primeros CPs, cuya varianza explicada acumulada está cerca del 75%.

Draft Soriano 860606967-image23-c.png

Figura 18. Scree test. Los puntos rojos muestran la varianza explicada por cada CP, mientras que las barras indican la varianza acumulada al incluir un nuevo CP.

Una vez retenido el número de componentes principales deseado se ha optado por utilizar la rotación ortogonal VARIMAX, un algoritmo que precisamente tiene por objetivo forzar la máxima varianza para cada componente principal y que ha sido ampliamente usado (Esteban et al. 2006; Martin-Vide et al. 2008).

En este punto, se inicia la explotación de los resultados rotados, es decir, sacar provecho a estas coordenadas que han sido optimizadas inicialmente por la ACP y, posteriormente, reajustadas con la rotación VARIMAX. Así pues, ahora ya tenemos nuestros casos originales simplificados y estandarizados, y dispersos en el espacio multivariante descrito por un número reducido de ejes de coordenadas (CP retenidas); ya es momento de encontrar una metodología que nos permita agruparlos.

Esteban et al. (2005) propone el método de las puntuaciones extremas, el cual se basa en identificar los scores positivos y negativos por encima y por debajo de +2 y -2, respectivamente, para cada componente principal retenida, con el objetivo de establecerlos como centroides de un futuro algoritmo de agrupaciones o clustering. Esta metodología puede establecer hasta dos centroides por cada CP. Por ejemplo, en el caso del CP1, en primer lugar, se buscará si hay días con un resultado > 2 en la primera coordenada (CP1+); en caso afirmativo, se recogerán todos los días que superen este valor y se calculará la media. Esta media no sólo se realiza con los casos del CP1, sino que se hace extensible al resto de CPs (Tabla 7). En caso de que no exista ningún día con un resultado superior a 2, este grupo queda directamente excluido. Así pues, n CPs pueden convertirse en hasta 2n tipos de circulación, en función de si se puede establecer un centro de coordenadas por exceso de umbral. Dicho de otra manera; si no existe algún caso real que se ajuste mínimamente a la estructura espacial explicada por esa CP (en su fase positiva o negativa), esta quedará descartada al considerarse un artefacto estadístico derivado del ACP.

Tabla 7. Ejemplo ficticio sobre el establecimiento de los centroides de cada grupo. En este caso se buscan los resultados superiores a 2 en el PC1, y se encuentran 4 fechas distintas. Se extraen los resultados de todas las componentes (rotadas) para estos días, y posteriormente se realiza la media para obtener el centro de clúster del primer grupo. Este proceso se realiza de forma iterativa para el extremo positivo y negativo de cada componente principal.
Días CP1 CP2 CP3 CP4 CP5 CP6
19/02/1993 2.50 1.49 -1.46 -1.00 -0.74 -0.43
06/11/2007 2.49 1.11 -0.14 -0.16 0.67 1.59
21/11/2008 2.70 0.47 -1.91 0.57 0.04 0.59
09/12/2010 2.59 0.87 -1.17 -1.82 0.75 -0.27
X (centroide) 2.57 0.99 -1.17 -0.60 0.18 0.37


Calculados los centros de cada clúster, se hace uso del método de agrupamiento K-means para clasificar todos los días de estudio en cada uno de los distintos grupos. El método de K-means se aplica sin iteraciones, puesto que se valora que los centros de clúster establecidos con el método de las puntuaciones extremas son altamente representativos de cada uno de los grupos (Esteban et al. 2012). Sin embargo, al aplicar el método, estos centros de clúster iniciales se recalculan para agrupar todos los días, y no dejarlos sin clasificar. Los nuevos centros mantendrán unos valores cercanos a los iniciales, pero ligeramente modificados, para así clasificar la totalidad de la muestra y obtener el catálogo sinóptico deseado.

Resultado de la clasificación sinóptica diaria: 2000-2020.

La implementación de la metodología detallada en el apartado anterior termina con el resultado presentado en la Figura 19, donde aparecen los principales patrones sinópticos que explican el 75 % de la variabilidad (5 PCs) del conjunto de días que conforman la serie diaria del período 2000-2020.

Draft Soriano 860606967-image24-c.png

Figura 19. Catálogo o clasificación sinóptica hallada. Los mapas muestran los valores promedios de presión en superficie (SLP, en hPa) para cada tipo sinóptico.

En él se muestran desde situaciones anticiclónicas (tipo 1) o de bajo gradiente bárico (tipo 7 i 10), favorecedoras de valores altos de O3, pero también situaciones que dispersan tales concentraciones como advecciones marcadas acompañadas de situaciones ciclónicas (tipos 2 y 5). Precisamente, estas situaciones ciclónicas son menos frecuentes que los casos con situaciones sin un gradiente marcado o anticiclónicas, durante este período 2000-2020 (Figura 20 y Figura 21).

Draft Soriano 860606967-image25-c.png
Figura 20. Frecuencias absolutas (en días) para cada patrón sinóptico (1-10), para la serie 2000-2020.
Draft Soriano 860606967-image26-c.png
Figura 21. Frecuencias relativas mensuales por patrón sinóptico (WT).

Anexo Caja Amarilla 3. Variables analizadas y descartadas

En el transcurso del proyecto se han analizado una serie de variables que se consideraron inicialmente de interés, pero que finalmente fueron descartadas por diversos motivos. A continuación, se describen las mismas y la motivación de su descarte.

  • Orografía – La orografía es un factor que se considera muy relevante en la concentración de O3. No en vano, existen áreas de Cataluña con comportamientos singulares debido a su orografía. Un caso paradigmático es la zona de “la Plana de Vic”, región caracterizada desde el punto de vista orográfico por una depresión alargada en dirección norte-sur (ver Figura 22). Debido a su particular configuración, en esta región se han registrado episodios ambientales con cierta frecuencia debido a la escasa ventilación atmosférica que tiene lugar en la misma. No obstante, a pesar de un factor de gran importancia, la naturaleza con la que se han configurado el motor de cálculo de este Reto, hacen que no sea conveniente introducir este tipo de parámetro como variable tipo input. Anteriormente se ha explicado que se ha elaborado un modelo de ML para cada estación de calidad de aire. Esto se traduce en que los parámetros que caracterizan cada estación sean invariantes a lo largo del tiempo, y por tanto se trate de valores fijos en lugar de variables. Este tipo de parámetros no aportarán información relevante en el proceso algorítmico de entrenamiento de los modelos. En el caso de haberse desarrollado otro enfoque completamente diferente en el que se generara un modelo de ML único para todas las estaciones, sí que daría cabida a introducir parámetros fijos de caracterización de cada estación como variables a considerar en el entrenamiento: altura de la estación, tipo de orografía, localización, etc. En conclusión, se determinó que no debía considerarse esta variable como variable input en este estudio.
Draft Soriano 860606967-image27.png
Figura 22. Localización geográfica de “plana de Vic” (derecha) y vista de detalle de la depresión orográfica que caracteriza dicha zona (izquierda)

* Ventilación atmosférica – La ventilación atmosférica es un proceso meteorológico el cual parece mantener una correlación elevada con la concentración de la mayor parte de contaminantes atmosféricos. De forma simplificada se puede citar que, a mayor ventilación atmosféricas, mayor dispersión de contaminantes y por tanto menor nivel de concentración. Existe un gran número de estudios que profundizan en estos procesos. Algunos estudios proponen indicadores mediante los cuales se trata de relacionar la ventilación atmosférica con la concentración de determinados compuestos químicos aéreos (Haeger-Eugensson et al. 2005, Genc et al. 2010, UNTEC 2014, Yim et al. 2015). No obstante, existe un gran número de factores que pueden afectar a este proceso, haciendo que su cálculo adquiera un carácter de gran complejidad. Esta dificultad se acentúa en zonas urbanas, donde la interacción entre la circulación del viento con grandes infraestructuras puede derivar en “efecto pared” (Yim et al. 2009). A pesar de que sí existen estudios que tratan de modelizar procesos de ventilación atmosférica, estos suelen basarse en técnicas de modelación numérica (Leelőssy et al. 2014), fuera del alcance de este Reto del proyecto Piksel. Por todo lo descrito, fundamentalmente asociado a la complejidad de cálculo de este tipo de procesos, se decidió no incluir esta variable como parte del presente estudio. No obstante, se considera que este efecto puede estar suficientemente representado con la inclusión de otras variables como son las relacionadas con la variable viento (tanto velocidad como dirección), así como el cálculo e inclusión del concepto de estabilidad atmosférica.

* Tráfico – El tráfico es otra de las variables cuya relación con la concentración de O3 se presupone de gran relevancia. La intensidad de tráfico, al igual que un gran número de actividades industriales de origen antropogénico, se consideran perjudiciales desde el punto de vista de calidad medioambiental al tratarse de factores generadores de emisiones contaminantes, algunas de ellas precursoras de la formación de O3. Debido a esta preocupación, desde hace varias décadas se trata de correlacionar intensidad de tráfico con concentración de contaminantes atmosféricos como el NO2 o el propio O3. (Reis et al. 2000, Ibarra-Berastegi et al. 2003, Munir et al. 2012, Munir et al. 2014, Nawahda 2016, Rodriguez-Rey et al. 2021). Tras llevar a cabo el estudio del arte del arte en este campo, se concluye que su análisis entraña una gran complejidad, al intervenir un gran número de parámetros y no poder aplicarse de forma evidente sobre un proyecto a escala regional, por lo que no se ha podido incluir este parámetro como input de los modelos de ML. No obstante, a pesar de no haberse considerado directamente esta variable como parte del estudio, se han incluido variables asociadas al tráfico, como son las mediciones de contaminantes precursores emitidos por este tipo de actividad, y variables temporales, las cuales pueden aportar información relativa a la variación de intensidad del tráfico a lo largo del tiempo.

* COVs – Tal y como se ha descrito en la introducción de este documento (ver sección ‘Problemática y contextualización del Reto’), la relación entre COVs y formación de O3 ha quedado demostrada en diferentes estudios, al originarse a partir de reacciones químicas entre COVs y óxidos de nitrógeno, en presencia de radiación solar (Feng et al. 2019, Wang et al. 2019). Si bien se considera una variable que puede estar altamente correlacionada con el output a predecir, y por tanto de gran interés para introducirse como input en los modelos de ML, no se ha podido considerar debido a la insuficiencia de series de observaciones de en la red XVPCA.

Conclusiones y flujo de trabajo futuro (extendido)

Puntos fuertes

  • Creación de una amplia base de datos con datos provenientes de redes de monitorización de calidad de aire y redes meteorológicas, a escala regional – Se dispone de una amplia base de datos compuesta por series de datos de variables relacionadas con los procesos de formación y concentración de O3 la cual podrá ser utilizada, no sólo con el objetivo de seguir mejorando el desarrollo del Reto 3, sino para otros posibles proyectos futuros. Esta base de datos está pretratada y corresponde a una casuística de casos de estudio amplia (47 estaciones de calidad de aire) con diferentes características. Los datos recopilados, basados en fuentes de información oficiales públicas, pondrán ponerse a disposición de la comunidad científica.
  • Estudio sobre el estado del arte relacionado con procesos complejos asociados al O3 y crear nuevas variables e indicadores – Otro de los logros del proyecto ha sido alcanzar un profundo grado de conocimiento sobre los procesos pueden afectar a la formación y concentración de O3: procesos fotoquímicos o físicoquímicos, fenómenos meteorológicos o actividades relacionadas con la acción antropogénica. También se han examinado los parámetros e indicadores más relevantes en este campo.
  • Desarrollo de modelos de ML para la predicción de niveles de concentración de O3 (extrapolable a otros contaminantes) – Se han construido modelos de ML para predecir niveles de O3 a nivel de estación de calidad de aire individual. La metodología llevada cabo puede ser generalizable para otras estaciones de monitorización que midan niveles de concentración de este contaminante, en otras zonas que disponga de los datos necesarios para su aplicación. En la construcción de estos modelos se han analizado variables poco frecuentes en este tipo de modelos, y se han testado varios algoritmos comparando su rendimiento, y optimizando sus metaparámetros. La metodología llevada a cabo puede ser extrapolable para la construcción de modelos de ML con el objetivo de predecir otras variables atmosféricas.
  • Creación de una herramienta con gran capacidad de automatización y estructura modular que facilita la integración de nuevo modelos – La mayor parte del código desarrollado se ha realizado con una alta capacidad de automatización de los procesos de cálculo y con una estructura modular. A cada sección principal del programa (preproceso, modelos de ML, predicción diaria, generación de resultados o simulación de escenarios sintéticos) se ha incorporado una serie de parámetros internos de cálculo que pueden definirse de forma flexible y rápida. Esta estructura permitirá experimentar con nuevos parámetros en fases posteriores del proyecto sin interferir en el resto del código y manteniéndose la robustez del motor de cálculo. Este aspecto cobra especial relevancia para el testeo de nuevos algoritmos de ML, los cuales podrán ser sustituidos por los actuales de una forma rápida y cuasi-automática.
  • Desarrollo de herramientas que ofrecen distintos servicios útiles para gestores de políticas de calidad de aire – Se ha desarrollado una herramienta digital de apoyo a los encargados de la administración en la gestión de este tipo de procesos. Está compuesta por un gran número de utilidades (visualización de datos históricos y episodios ambientales, predicción futura a corto plazo, simulación de escenarios sintéticos, etc.) que permitirán mejorar las acciones tomadas por los gestores.

Puntos débiles

A continuación, se exponen los puntos más débiles e impedimentos que se han detectado a lo largo del proyecto, con el fin de que puedan mejorarse en una proyección futura a corto-medio plazo.

Aspectos relacionados con los datos

* Falta de información sobre variables relevantes Un ejemplo claro es la concentración de COVs. Se descartó su integración en esta fase del proyecto por falta de representatividad, al detectar que esta variable sólo era monitorizada en unas pocas estaciones.
* Distribución espacial de la red de monitorización poco homogénea – De las 47 estaciones de calidad tratadas en el Reto, muchas de ellas están concentradas en zonas urbanas (la mayoría en la zona de Barcelona), mientras que otras áreas tienen muy pocas estaciones. En términos de ZQA, se puede observar esta irregularidad en la distribución espacial. Por ejemplo, el ‘ZQA 1 – Área de Barcelona’, dispone de 12 estaciones de medición de O3, mientras que las ZQAs del interior de Cataluña o cercanas a la zona pirenaica, tan sólo tienen dos (‘ZQA 11 – Pirineu Oriental’, ‘ZQA 13 – Prepirineu’ o ‘ZQA 14 – Terres de Ponent’) o incluso una estación (‘ZQA 12 – Pirineu Occidental’). Esto se traduce en una falta de representatividad en dichas zonas, que puede afectar sobre todo en el cálculo de resultados con resolución espacial distribuida.
* Falta de series de datos sobre predicciones futuras – En el caso de haber dispuesto de series de datos de entrada basadas en predicciones futuras, se podría haber abordado el problema de entrenar los modelos de ML con valores de entrada basados en dichas predicciones. Si bien, en el caso de variables meteorológicas sí se dispone de predicciones de este tipo (hasta 72 horas vista), la resolución espacial con la que se ofrecen estos datos es a nivel de municipio, dificultando la integración directa en la predicción de valores de O3 a nivel de estación individual. Por otro lado, no se han encontrado valores relativos a predicciones futuras relativas a las variables relacionadas con otros contaminantes atmosféricos, imposibilitando la aplicación de esta metodología.
* Problemas de descarga a través de APIs: falta de datos puntuales En ocasiones, se han detectado fallos puntuales de descarga de datos a través de alguna de las APIs oficiales con las que conecta la herramienta para efectuar las predicciones diarias. Ante estos problemas se han introducido mecanismos para que la herramienta siga siendo robusta, mediante acciones de “completado de datos”.

Aspectos relacionados con la precisión de los modelos de ML

* Error medio global de los modelos – Atendiendo a los resultados de error ofrecidos por los modelos de ML (considerando como métrica el RSME), encontramos valores de entre 8 y 18 (µg/m3), siendo el error medio global de en torno a 13.4.
* Menor precisión en la predicción de valores extremales – Los primeros análisis llevados a cabo sobre la distribución de este error, en función del rango de valores a predecir, parecen indicar un empeoramiento en la capacidad predictiva de los modelos para valores extremales (tanto muy bajos como sobre todo muy altos). Esto es coherente con los procesos algorítmicos de ML utilizados, dado que existe falta de representatividad de dichos rangos en las series de datos de entrenamiento. Desde el punto de vista práctico, esto supone un inconveniente cuando el objetivo principal de la herramienta desarrollada se utilice para tratar de predecir valores asociados a episodios ambientales, los cuales, en general tendrán un mayor error que los valores medios de concentración de O3.

Flujo de trabajo a corto plazo (julio 2022 – diciembre 2022)

  • Buscar otras fuentes de información, que puedan complementar las series de datos utilizadas en esta fase del proyecto. Se incidirá especialmente en variables que se consideran relevantes en los procesos de formación y concentración de O3, pero de las cuales no se disponga información suficiente o fiable en la actualidad (véase, series de datos sobre COVs).
  • Aplicar otros algoritmos y metodologías de ML, que puedan mejorar los resultados de precisión obtenidos, sobre todo para valores extremales. En primera instancia, se probarán nuevos algoritmos cuya algoritmia interna no se fundamente en conjuntos de árboles. Posteriormente, se evaluará la conveniencia de integrar nuevos procesos metodológicos enfocados específicamente en mejorar este tipo de predicción extremal (como puede ser la aplicación de métodos relacionados con “Imbalanced Regression”).
  • Realizar tareas de testeo de la herramienta, detectando y corrigiendo posibles errores, y aportando mayor robustez al código de programación.

Flujo de trabajo a medio plazo (enero 2023 – diciembre 2023)

  • Continuar con el desarrollo de la utilidad asociada a la actual opción de cálculo “¿Qué pasaría si?”. Redefinir aquellos aspectos requeridos por parte del usuario final, para que esta opción de respuesta de la mejor forma posible a los requerimientos por parte de los gestores.
  • Implementar otro enfoque a la hora de construir los modelos de ML, utilizando series de datos de las variables de entrada basadas en predicciones futuras. Este enfoque se basa en la predicción de valores de O3 futuros (correspondientes al día 'd+1' o 'd+2') en base a valores de predicciones (día ‘d+1’ o ‘d+2’, respectivamente) de la mayor parte de variables de entrada. También se analizará la posibilidad de combinar series de datos basadas en predicciones (día ‘d+1’ o ‘d+2’) con otras basadas en mediciones de día ‘d’ o ‘d-1’, como es el actual enfoque. Convendría un análisis cuidadoso sobre el comportamiento de los modelos de ML en el caso de "mezclar" valores de distintas variables asociados a diferentes periodos temporales.
  • Valorar la incorporación, como variables de entrada, de resultados generados mediante modelos numéricos. Este enfoque se abordará en el caso de detectar la necesidad de incorporar variables más complejas, de las cuales no se disponga datos monitorizados, o de los que no se dispongan series asociadas a predicciones futuras.

Referencias

  • BOE – Boletín Oficial del Estado (2008). Directiva 2008/50/CE del Parlamento Europeo y del Consejo. Referencia: DOUE-L-2008-81053.
  • Feng, R., Zheng, H. J., Zhang, A. R., Huang, C., Gao, H., & Ma, Y. C. (2019). Unveiling tropospheric ozone by the traditional atmospheric model and machine learning, and their comparison: A case study in Hangzhou, China. Environmental pollution, 252, 366-378.
  • Esteban, P., Martin-Vide, J. and Mases, M. (2006): Daily atmospheric circulation catalogue for western Europe using multivariate techniques. Int. J. Climatol., 26: 1501-1515. https://doi.org/10.1002/joc.1391
  • Esteban, P., Jones, P.D., Martín-Vide, J. and Mases, M. (2005): Atmospheric circulation patterns related to heavy snowfall days in Andorra, Pyrenees. Int. J. Climatol., 25: 319-329. https://doi.org/10.1002/joc.1103.
  • Ghazali, N. A., Ramli, N. A., Yahaya, A. S., Yusof, N. F. F. M., Sansuddin, N., & Al Madhoun, W. A. (2010). Transformation of nitrogen dioxide into ozone and prediction of ozone concentrations using multiple linear regression techniques. Environmental monitoring and assessment, 165(1), 475-489.
  • Genc, D. D., Yesilyurt, C., & Tuncel, G. (2010). Air pollution forecasting in Ankara, Turkey using air pollution index and its relation to assimilative capacity of the atmosphere. Environmental monitoring and assessment, 166(1), 11-27.
  • Haeger-Eugensson, Borne, Chen & Persson (2005). Development of a New Meteorological Ventilation Index for Urban Air Quality Studies. Technical Report, IVL Swedish Environmental Research Ltd.
  • Ibarra-Berastegi, G., Madariaga, I., Agirre, E., & Uria, J. (2003). Short-term forecasting of ozone and NO2 levels using traffic data in Bilbao (Spain). WIT Transactions on The Built Environment, 64.
  • Leelőssy, Á., Molnár, F., Izsák, F., Havasi, Á., Lagzi, I., & Mészáros, R. (2014). Dispersion modeling of air pollutants in the atmosphere: a review. Central European Journal of Geosciences, 6(3), 257-278.
  • Martin-Vide, J., Sanchez-Lorenzo, A., Lopez-Bustins, J. A., Cordobilla, M. J., Garcia-Manuel, A., and Raso, J. M. (2008): Torrential rainfall in northeast of the Iberian Peninsula: synoptic patterns and WeMO influence, Adv. Sci. Res., 2, 99–105, https://doi.org/10.5194/asr-2-99-2008.
  • MITERD - Ministerio para la Transición Ecológica y el Reto Demográfico (2021). “Evaluación de la Calidad del Aire en España. Informe Anual 2020”. Ministerio para la Transición Ecológica y el Reto Demográfico, Secretaría General Técnica. Centro de Publicaciones. NIPO: 665-21-045-X.
  • Montenegro, E. & Luján-Pérez, M. (2017). Análisis de la variación estacional de la contaminación atmosférica y su relación con variables climáticas en el valle central de Cochabamba, Bolivia. Acta Nova, 8(3), 451-466. ISSN: 1683-0768.
  • Munir, S., Chen, H., & Ropkins, K. (2012). Modelling the impact of road traffic on ground level ozone concentration using a quantile regression approach. Atmospheric environment, 60, 283-291.
  • Munir, S., Chen, H., & Ropkins, K. (2014). Characterising the temporal variations of ground-level ozone and its relationship with traffic-related air pollutants in the United Kingdom: A quantile regression approach. International Journal of Sustainable Development and Planning, 9(1), 29-41.
  • Nawahda, A. (2016). An assessment of adding value of traffic information and other attributes as part of its classifiers in a data mining tool set for predicting surface ozone levels. Process Safety and Environmental Protection, 99, 149-158.
  • OMS – Organización Mundial de la Salud (2005). Air quality guidelines: global update 2005: particulate matter, ozone, nitrogen dioxide, and sulfur dioxide. Organización Mundial de la Salud.
  • Pasquill, F. (1961). The estimation of the dispersion of windborne material. Met. Mag., 90, 33.
  • Reis, S., Simpson, D., Friedrich, R., Jonson, J. E., Unger, S., & Obermeier, A. (2000). Road traffic emissions–predictions of future contributions to regional ozone levels in Europe. Atmospheric Environment, 34(27), 4701-4710.
  • Rodriguez-Rey, D., Guevara, M., … & García-Pando, C. P. (2021). A coupled macroscopic traffic and pollutant emission modelling system for Barcelona. Transportation Research Part D: Transport and Environment, 92, 102725.
  • UNTEC (2014). “Determinacion de un índice de remoción de contaminantes atmosféricos para el territorio nacional”. Informe técnico. Universidad y Tecnologías, Fundación para la transferencia tecnológica.
  • Wang, N., Lyu, X., Deng, X., Huang, X., Jiang, F., & Ding, A. (2019). Aggravating O3 pollution due to NOx emission control in eastern China. Science of the Total Environment, 677, 732-74
  • Yim, S. H., Fung, J. C. H., Lau, A. K. H., & Kot, S. C. (2009). Air ventilation impacts of the “wall effect” resulting from the alignment of high-rise buildings. Atmospheric Environment, 43(32), 4982-4994.
  • Yim, S. H., Fung, J. C. H., & Ng, E. Y. (2014). An assessment indicator for air ventilation and pollutant dispersion potential in an urban canopy with complex natural terrain and significant wind variations. Atmospheric Environment, 94, 297-306.
Back to Top

Informació del document

Publicat a 06/07/22
Presentat el 05/07/22

llicència: CC BY-NC-SA license

Categories

Ozó troposfèric

Puntuació document

0

Visites 0
Recomanacions 0