la qualité de l'air dans les Pays de la Loire

contexte : modélisation et prévision de l’ozone

Air Pays de la Loire

L’organisme Air Pays de la Loire est agréé par le MEEDDAT pour la surveillance de la qualité de l’air dans les Pays de la Loire. La surveillance est assurée à l’aide d’observations et d’outils de modélisation. Elle couvre l’ensemble de la région des Pays de la Loire. Air Pays de la Loire mène deux missions d'intérêt général : surveiller et informer.

l’ozone O3

L’ozone est un constituant de l’air nécessaire en haute altitude pour protéger les espèces vivantes des rayonnements solaires mais est considéré comme un polluant au niveau du sol (« mauvais ozone »). L'ozone en excès présente des effets indésirables sur la santé et la nature.

L’ozone est un polluant de l’air présentant une élévation des concentrations en période estivale, saison pendant laquelle une vigilance particulière est mise en place. C’est un polluant secondaire issu de réactions chimiques entre des polluants primaires dont le dioxyde d’azote et les Composés Organiques Volatils. Ces réactions sont amplifiées par le rayonnement ultraviolet du soleil. L’ozone est un polluant à grande échelle présentant des concentrations homogènes sur de larges zones.

L’ozone fait l’objet d’une surveillance assurée par des relevés de concentrations en continu effectués sur une vingtaine de sites de mesure de la pollution de fond (typologie : urbaine, périurbaine, rurale) ainsi que par un système régional de modélisation et de prévision des concentrations d’ozone. Des bulletins de recommandation et d’information sont diffusés en cas de dépassement du seuil 180 µg/m³.

La procédure d’information et de recommandation sur la région a été mise en œuvre à 10 reprises en 2003 au cours de l’été caniculaire et 3 fois en 2006.

IRIS

La plate-forme régionale de prévision et de cartographie IRIS fournit des informations sur la qualité de l’air de la région des Pays de la Loire. La modélisation concerne deux principaux polluants : l’ozone (O3) et le dioxyde d’azote (NO2). IRIS couvre deux zones géographiques avec des résolutions différentes : le Grand Ouest (GOUEST – 15km) et les Pays de la Loire (APL - 5 km). La prévision est disponible pour 4 échéances allant de la veille jusqu’au surlendemain avec un pas de temps horaire.

La plate-forme intègre trois modèles déterministes décrivant les émissions (Airemis), la météorologie (MM5), la photochimie (Chimère). Iris est alimenté également par des données grande échelle de météorologie (NCEP) et de pollution (PREV’AIR) afin de prendre en compte les imports de pollution à l’échelle nationale et continentale.

Elle fournit des informations sous forme de cartes de concentrations moyennes horaires. Des cartes de moyennes journalières et de maximums journaliers sont également produites par IRIS à partir des données horaires. Les observations disponibles aux sites de mesure sont utilisées afin de corriger les cartes de prévision pour la veille et le jour même.

La plate-forme Iris répond au besoin d’anticiper les pointes de pollution en donnant une information fiable jusqu’à deux jours avant l’apparition des pointes. Ces prévisions sont réactualisées chaque jour.

OPTIRIS

L’étude se positionne dans le cadre du projet OPTIRIS (OPTimisation de la plate-forme IRIS) qui vise à améliorer le modèle régional IRIS présentant des scores perfectibles sur les niveaux de pointes. Un bilan sur les prévisions au cours de l’été 2006 a mis en avant des scores légèrement inférieurs à ceux du modèle national PREV'AIR ainsi que des prévisions déficientes de dépassement du seuil d’information pour l’ozone, les 3 jours avérés de dépassement du seuil d’information (180 µg/m3) n’ayant été prévus ni par IRIS ni par PREV’AIR.

Un diagnostic a débuté fin juillet 2006 avec la société ACRI-ST afin de comprendre les défauts du modèle pour ensuite l’améliorer. Les premières actions ont consisté en la mise à jour du modèle de photochimie CHIMERE (V200501H) et la modification du module émission. Après recalcul des prévisions, les scores étaient améliorés mais encore légèrement inférieurs à ceux de PREV’AIR.

Une dernière action d’amélioration a été la mise en place par Air Pays de la Loire et en collaboration avec la société ACRI-ST d’un post-traitement des sorties « brutes » du modèle IRIS par adaptation statistique.

Il a été montré que cette approche améliorait les scores pour la prévision des concentrations en ozone mais que le modèle conservait le défaut de prévision des pointes de pollution.

objectifs : amélioration de l’adaptation statistique

L’objectif de l’étude est de proposer une nouvelle méthode de correction des prévisions d’ozone par adaptation statistique qui corrigerait les défauts de la méthode actuellement en place.

Un objectif d’opérationnalité de la méthode a guidé toute l’étude. Le modèle retenu à la fin de cette analyse tient compte de cet objectif.

moyens

adaptation statistique

L’adaptation statistique est une méthode de prévision statistique utilisant des sorties de modèles déterministes comme variables explicatives (météo, chimie,…). Cette technique permet notamment de corriger une sortie de modèle déterministe en déterminant une "climatologie" de l’erreur de prévision. Elle est utilisée principalement dans le domaine de la pollution de l’air et en particulier pour prévoir les concentrations d’ozone. Une étude bibliographique a montré l’utilisation des méthodes de régression, CART, de séries temporelles ainsi que de méthodes de type « boîte noire » comme les forêts aléatoires ou les réseaux de neurones.

méthodes statistiques

Dans cette étude, nous nous sommes intéressés principalement aux méthodes de régression paramétriques qui permettent d’expliquer simplement une variable réponse à partir de variables explicatives. D’autres approches ont été testées, notamment la méthode des forêts aléatoires qui est une amélioration de la méthode CART.

Les modèles de régression permettent en pratique de corriger la prévision déterministe aux niveaux des sites de mesure. Une procédure de krigeage des innovations permet d’étendre la correction à l’ensemble de la grille de prévision en tenant compte de la structure spatiale de l’erreur de modélisation déterministe.

données

L’historique de mesure et prévision utilisé correspond à la période de l’été 2006 (15 juin au 15 septembre) ainsi que la période allant du 15 juin 2007 au 15 juin 2008. Ces dates correspondent à la période de fonctionnement de la plate-forme IRIS avec la version V200501H du modèle CHIMERE. Les données sont disponibles pour 23 stations mesurant la pollution de fond.

Plus précisément, l’inventaire des variables, donné en annexe, est constitué de l’ozone mesuré aux stations, des prévisions déterministes de l’ozone et du dioxyde d’azote issues de la plate-forme IRIS, de l’heure TU ainsi que de 13 variables correspondant à des prévisions météorologiques issues du modèle MM5 (température, hauteur de couche limite, précipitation,…).

résultats

variable réponse

Une première étape a consisté à choisir la variable à expliquer dans le modèle de régression : soit la concentration d’ozone observée, soit l’erreur de modélisation déterministe (erreur modèle = [O3]prévue - [O3]mesurée). Concernant les hypothèses nécessaires à l’utilisation de la régression, l’erreur modèle apparaît comme le meilleur candidat. La variable erreur modèle possède une distribution plus gaussienne (« en cloche ») que l’ozone mesuré. L’utilisation de l’erreur modèle comme variable à expliquer permet également de réduire les phénomènes de colinéarité entre les variables explicatives. La prévision déterministe est en effet très corrélée à la température ou à la hauteur de couche limite.

étude de l’ensemble des stations

Dans cette partie nous présentons les résultats de l’étude sur l’ensemble des stations. Nous disposons d’un jeu de données pour chacune des 23 stations, chaque modèle (APL et GOUEST) et chaque échéance (J-1, J, J+1, J+2). Nous montrons l’homogénéité spatiale du comportement de l’ozone et les similitudes du modèle déterministe de prévision quel que soit le domaine et la profondeur de prévision.

La figure ci-dessous montre l’homogénéité de l’erreur modèle d’un point de vue géographique.

Figure 1 : boxplots de l’erreur modèle par station pour le domaine APL et l’échéance J+1

Nous étudions maintenant l’influence du domaine de calcul et de l’échéance de prévision sur la distribution de l’erreur de modélisation déterministe. Intéressons nous au cas des stations de Rezé (BALIO3) et de Pornichet (EMILO3). Ces exemples sont représentatifs des situations aux autres stations.

Figure 2 : boxplots des erreurs modèle sur deux stations (BALIO3 à gauche et EMILO3 à droite) selon l’échéance de prévision et le domaine de calcul

Au regard de ces graphiques, nous pouvons conclure à une influence non significative de l’échéance de prévision dans la répartition de l’erreur modèle pour un même domaine de calcul. Le domaine de calcul influe en termes de moyenne sur la l’erreur modèle. Un biais est visible entre les distributions des erreurs modèles des domaines APL et GOUEST pour la station BALIO3. Ce biais se retrouve pour d’autres stations. En revanche, la dispersion de l’erreur ne dépend ni de l’échéance ni du domaine.

Une étude des relations inter-variables a également montré une homogénéité des liens sur l’ensemble des stations et ceci pour les deux domaines et les quatre échéances. En particulier, l’erreur modèle est expliquée par les mêmes variables selon les situations.

étude de modèles sur un cas test

Après avoir montré les faibles différences entre les sites, modèles et échéances, nous nous intéressons uniquement au cas de la station rurale la Tardière (FOLIO3) à l’échéance J+1 et pour le domaine de calcul APL. L’objectif de cette partie est de tester différents modèles prédictifs sur ce cas particulier afin d’améliorer la méthode d’adaptation en place.

choix de l’échantillon test

Dans un premier temps, les données ont été séparées en deux échantillons distincts, un échantillon de calibration servant à construire les modèles et un échantillon test servant à les comparer entre eux. L’échantillon test est composé de 3 périodes d’une durée d’un mois soit environ un quart des données disponibles. Ces mois ont été choisis de manière à être représentatifs de l’ensemble des situations rencontrées en termes de concentrations d’ozone. Nous les décrivons dans le tableau
ci-dessous :

La figure ci-dessous montre que le choix effectué permet d’avoir un échantillon test représentatif des situations rencontrées en termes d’erreur de prévision déterministe.

Figure 3 : résumé graphique de l’erreur modèle

sélection des variables explicatives

Une étude des corrélations entre variables a montré que les variables pbl_hgt (hauteur de couche limite, 26%), t2 (température, 15%), vnorm (norme de la vitesse du vent, 13%), NO2 (dioxyde d’azote prévu, 12%) et v10 (composante Nord-Sud de la vitesse du vent, 12%) sont significativement corrélées avec l’erreur modèle. Les significations des variables ainsi que les valeurs des corrélations avec l’erreur sont données entre parenthèses.

Les autres variables explicatives se sont révélées peu ou pas corrélées avec la variable réponse et ont donc été rejetées de l’étude. L’analyse des corrélations a montré également des phénomènes de colinéarité (corrélations significatives entre variables explicatives), notamment entre les variables température et taux de mélange ainsi qu’entre les variables épaisseur optique et atténuation photolytique qui sont liées par une relation déterministe.

L’heure de prévision a également fait l’objet d’une étude afin d’analyser son impact sur l’erreur modèle. La figure ci-dessous représente la distribution de l’erreur sous forme de boîtes à moustache (boxplots) en fonction de l’heure TU.

Figure 4 : boxplots de l’erreur modèle en fonction de l’heure de la journée

L’erreur modèle évolue en moyenne en fonction de l’heure de la journée. Elle est la plus élevée le matin et la plus faible en fin d’après-midi. Une analyse de la variance a permis de prouver la pertinence de ce facteur heure. Dans un but de réduire le nombre d’effets à estimer concernant ce facteur, des comparaisons multiples ont été effectuées entre les moyennes aux différentes heures. Ces comparaisons ont permis d’agglomérer les heures dans des groupes (modalités) homogènes.

Ces groupes définissent un nouveau facteur heure_groupe à 4 modalités :

• modalité 0 : 22 heures à 4 heures ;
• modalité 1 : 5 heures à 9 heures ;
• modalité 2 : 10 heures à 13 heures ;
• modalité 3 : 14 heures à 21 heures.

L’étude des facteurs jours de la semaine et mois de l’année n’a pas permis de retenir ces variables.

dépendance temporelle de l’erreur

Une analyse de la dépendance de l’erreur avec son passé proche a montré un phénomène marqué d’autocorrélation périodique selon une période de 24 heures. Ce phénomène décroît globalement au fur et à mesure que le décalage augmente mais reste significatif pour des décalages très grands.

Figure 5 : autocorrélations de la série des erreurs modèles en fonction du décalage (lag)

Afin de prendre en compte cette dépendance temporelle de l’erreur modèle, nous avons introduit les variables décrivant l’erreur retardée par pas de 24 heures. La variable erreur_Jm1 correspond ainsi à l’erreur modèle retardée de 24 heures. Une régression linéaire multiple de l’erreur sur son passé a permis de ne retenir que les variables erreur_Jm1, erreur_Jm2, erreur_Jm3 et erreur_Jm7 correspondant respectivement aux erreurs retardées de 24, 48, 72 et 168 heures.

estimation du modèle de régression retenu

L’ensemble des variables explicatives retenues a été utilisé afin d’expliquer les variations de l’erreur modèle au moyen d’un modèle de régression linéaire multiple. L’équation du modèle est ainsi :

êrreur~erreur_Jm1+erreur_Jm2+erreur_Jm3+erreur_Jm7+NO2+pbl_hgt+t2+v10+vnorm+heure_groupe

Afin de retenir uniquement les variables les plus pertinentes et également de réduire les phénomènes de colinéarité, une sélection des variables par procédure stepwise utilisant le critère BIC a été réalisée sur ce modèle. La variable erreur_Jm3 a ainsi été rejetée. Le modèle final explique près de 25% de l’erreur modèle. Cette information vient en complément de celle apportée par la prévision déterministe.

Une analyse graphique des résidus a permis de valider l’ensemble des hypothèses de la régression exceptées celles concernant l’indépendance des résidus. Un phénomène d’autocorrélation subsiste en effet pour les décalages inférieurs à 24 heures.

Ce modèle a ensuite été testé avec succès sur l’échantillon test, les scores obtenus étant meilleurs après adaptation pour l’ensemble des concentrations ainsi que pour les pointes.

scores globaux

Nous comparons l’adaptation statistique actuellement en place à la méthode retenue après étude. Chacun des modèles de régression a été ajusté station par station et sur chacune des situations (2 domaines x 4 échéances). Cet ajustement a été effectué sur la base de l’échantillon de calibration défini dans le cadre de l’étude de la station la Tardière, en utilisant une procédure de sélection des variables pas à pas (stepwise). Cette procédure consiste à sélectionner le modèle qui minimise un critère d’information (AIC, BIC,…) en retirant et ajoutant des variables de manière itérative.

Les deux approches sont détaillées dans le tableau ci-dessous :

Les prévisions, dites adaptées, issues de ces deux méthodes, ainsi que les prévisions déterministes sont comparées aux concentrations observées à l’aide d’indicateurs statistiques horaires et journaliers (scores).

Ces scores sont déterminés sur les données de l’échantillon test précédemment défini. L’échantillon est filtré selon un seuil (seules les concentrations supérieures à ce seuil sont conservées) afin de permettre de se focaliser sur les niveaux élevés (seuil 120) ou de s’affranchir des erreurs commises sur les faibles concentrations qui peuvent biaiser les scores (seuil 10).

comparaison avec la méthode en place

Nous présentons ici les scores Biais Normalisé et RMSE pour les seuils 10 et 120 à l’échéance de prévision J+1. Ces filtrages conduisent à considérer des échantillons de taille 30000 environ pour le seuil 10 et 1500 pour le seuil 120.

Le biais Normalisé (en %) décrit la tendance du modèle à sous ou sur-estimer les concentrations relativement à celles-ci tandis que la RMSE décrit l’erreur moyenne. La RMSE est plus sensible aux forts écarts que l’erreur absolue.

Légende des graphiques ci-dessus

Une amélioration significative est apportée par cette nouvelle méthode qui fait intervenir dans les modèles l’erreur passée, la prévision de NO2, le facteur heure_groupe ainsi que les variables météorologiques les plus pertinentes (t2, pbl_hgt, v10, vnorm). Les scores horaires sont améliorés avec cette nouvelle méthode pour les 2 seuils de filtrage. En particulier, la nouvelle méthode améliore les prévisions brutes sur les fortes concentrations (seuil 120) alors que la méthode en place les détériore.

La nouvelle approche donne des résultats proches de ceux de la méthode en place en termes de prévision de dépassements journaliers du seuil 170. Cette comparaison est néanmoins à relativiser du fait du faible nombre de cas d’étude (15). Les comparaisons entre les méthodes sur le seuil 150 (70 cas) permettent d’apprécier le gain apporté par la nouvelle approche. Nous donnons les résultats ci-dessous pour l’échéance de prévision J+1.

Alors que l’adaptation en place dégrade les prévisions de dépassements du seuil 150, la nouvelle méthode les améliore significativement (+8% pour APL et +6% pour GOUEST). Les conclusions sont les mêmes pour l’échéance J. En revanche, les adaptations ne parviennent pas à améliorer les prévisions brutes de dépassements pour l’échéance J+2 pour le seuil 150.

conclusions

Dans le cadre du projet OPTIRIS et du stage, la nouvelle approche consistant à régresser l’erreur modèle en utilisant les variables erreur passée, météorologiques et prévision de NO2 se révèle comme la meilleure méthode d’adaptation qui ait été testée. Cette approche corrige les prévisions déterministes d’ozone sur l’ensemble de la plage de concentration ainsi que sur les pointes. Elle améliore ainsi la méthode d’adaptation en place. Son implémentation ne nécessite que peu de modifications en raison de sa proximité avec la méthode en place. Le principal défaut de cette méthode concerne l’hypothèse d’indépendance des résidus du modèle de régression. Malgré cette restriction, les modèles obtenus possèdent de bonnes propriétés prédictives. Le développement de méthodes plus complexes comme les forêts aléatoires ou la régression GAM ne semble pas pertinent au vu des résultats obtenus sur le cas test.

perspectives

L’amélioration des prévisions météorologiques ainsi que des prévisions de NO2 est une voie d’amélioration de la prévision d’ozone.

L’utilisation de nouvelles variables explicatives est également une piste d’amélioration. En ce sens, il avait été envisagé en début de stage d’utiliser une approche « multimodèles » consistant à considérer plusieurs sources de prévisions de l’ozone (PREV’AIR, IRIS, MOCAGE) dans le but de compenser les défauts propres à chaque modèle. La prise en compte de l’import d’ozone en tant que variable explicative est également une piste qui pourrait être explorée.

Une autre piste d’amélioration reste l’accroissement de l’historique afin de prendre en compte davantage d’épisodes de pollution. Pour réaliser ceci, un recalcul des prévisions avec la version actualisée du modèle CHIMERE sur les années antérieures (2005) est à envisager.

L’utilisation de modèles spécifiques à chaque heure est une approche qui pourrait apporter un gain en termes de prévision en prenant en compte les variations au cours de la journée de manière plus fine. La seule contrainte liée à cette approche concerne la « lourdeur » de sa mise en place au sein de la plate-forme qui pourrait rendre cette approche non implémentable.

Dans une optique de prévision des seules pointes de pollution par l’ozone, la segmentation des données selon la concentration observée semble apporter une réponse à la fois simple à mettre en place et efficace. La prévision des maximums journaliers pourrait également être une approche à exploiter.