Le programme lidar_ozone_actris_v01.py
Les grandes étapes du traitement lidar
Généralités
Tout le programme a été concu autours de grandes étapes de traitement, elles mêmes sous classé en deux grandes parties. En effet, la donnée brute lidar est une grandeur physique compté par de détecteurs. Cette donnée est grandement impacté par l’instrument qui a réalisée la mesure. Le but de la première étape, appelée prétraitement consiste à supprimer tous les effets dus à l’instrument. On va donc retrouver la suppression du fond de ciel, la correction de la saturation des détecteurs et la correction du facteur géométrique. Globalement, la grandeur physique en sortie de ce premier traitement est similaire à celle d’entrée.
La deuxième étape, appelée courament traitement consiste à transformer la grandeur physique mesurée en grandeur géophysique recherché. Par exemple, dans le cas d’un lidar ozone, on va transformer un nombre de photon reçu en concentration de molécule d’ozone.
Le prétraitement est unique en fonction de l’instrument, en revanche, le traitement doit être identique en foncion de la donnée géophysique recherchée.
Fond de ciel, bruit du détecteur et signal induit
On appelle “fond de ciel” 3 phénomènes spécifiques :
la lumière issue du ciel : cela peut être la lumière du soleil, de la lune, des étoiles, de la lumière artificielle des villes, et du batiment le plus proche. Globalement, c’est un bruit blanc, dépendant de la longueur d’onde, mais non dépendant de l’altitude. Les variations en fonction du temps sont lentes, mais non nulles.
le bruit de fond de détection : de base, le détecteur et la carte d’acquisition ne sont pas parfait, et, même plongés dans le noor total, le détecteur ‘compte’ quelques photons de temps en temps. C’est en général un bruit blanc, dépendant du détecteur, mais non dépendant de l’altitude. En théorie, ce bruit n’évolue pas dans le temps. Il est cependant intéressant à vérifier de temps en temps, car le vieillissement de celui-ci peut avoir un impact sur cette valeur.
le signal induit : c’est un signal qui est dépendant de la lumière reçue. Globalement, c’est photons rentrant dans le détecteur sont captés par celui-ci et relachés plus tard. Ce phénomène peut être important lorsque le signal reçue sur les détecteurs est important. Cette fois ci, le signal n’est pas un bruit blanc car dépendant de l’altitude
Dans la partie “fond de ciel”, il s’agit de supprimer ces 3 phénomènes en même temps. En théorie, une simple moyenne est suffisant, mais du fait su signal induit, on va souvent utiliser une régression linéaire sur le fond du signal, afin de corriger la valeur de signal induit. Cette correction linéaire doit être faible, une forte pente de la régression indique un problème important de signal induit, soit un autre problème.
Example of background correction
Désaturation
Le détecteur n’est pas parfait, et il peut être éblouit s’il recoit trop de signal. Cela veut dire qu’il a une valeur seuil de comptage. Cependant, avant d’être totalement éblouie, il y a une zone où il compte “mal” le nombre de photons, et plus particulièrement, où le nombre de photon comptés est infénieur au nombre de photons reçus. Il existe plusieurs équations mathématiques permettant de simuler cette saturation du détecteur, souvent paramétrisables, c’est-à-dire où il faut choisir un paramètre, qui peut évoluer avec le veillissement du détecteur.
Etude de la désaturation des signaux par la méthode de Muller
Facteur géométrique
Le facteur géométrique correspond au taux de visualisation du laser par le télescope : lorsqu’il est égal à 0, le laser n’est pas dans le champ de vue du télescope, lorsqu’il est égale à 1, il est totalement, et entre les deux lorsque le laser est partiellement dans le champ de vue du télescope.
Exemple de facteur géométrique.
Gluing des données
Cette technique permet de relier des signaux entre eux. En effet, les détecteurs n’était pas parfait, ils ne peuvent travailler que sur une gamme d’altitude. Pour couvrir une large gamme d’altitude, il faut donc utiliser plusieurs détecteurs, et donc relier les signaux par la suite.
Exemple de gluing des signaux, entre une voie haute (en orange) et une voie basse (en bleu). A gauche, on peut voir les signaux bruts, au milieu les signaux glués, et à droite le type de raccordement effectué (ici, un système en cosinus et sinus carrés. En plus du raccordement, il y a eu un décalage de la voie haute sur la voie basse.
Filtrage des signaux
Le filtrage des signaux lidars va permettre de réduire le bruit, et de pouvoir aller plus haut en altitude. Comme le bruit du issue du détecteur est dépendant de la puissance du signal, et donc de l’altitude, il est courant d’utiliser un filtrage dont le nombre de point du filtre (et donc sa résolution) dépend de l’altitude. Dans les basses altitudes, lorsque le signal est fort, on filtre peu, mais dans les hautes altitudes, lorsque le signal est faible, on filtre beaucoup plus.
Le fait de filtrer va permettre de réduire le bruit au détriment de la résolution. En effet, pour filtrer, il faut utiliser un certain nombre de points, et donc le nombre de points indépendant va diminuer.
Exemple du filtrage d’un signal. A gauche, le signal filtré (en orange) est moins bruité sur le signal brut (en bleu). Au niveau de l’erreur, on voit que le signal filtrée à une erreur plus faible. Enfin, à droite, on voit la résolution, qui passe de 50 mètres à 4 km d’altitude à presque 400 mètres à 40 km d’altitude.
Le traitement aérosols
Il s’agit de récupérer un profil de rétrodiffusion par la méthode Klett. Il s’applique seulement sur les voies Rayleigh Mie, à des longueurs d’onde non perturbées par d’autres molécules (comme la voie ON de l’ozone).
On va souvent visualiser le signal en backscatter ratio, c’est-à-dire le rapport entre le signal diffusé par les aérosols et le signal rétrodiffusé total (aérosols et molécules de l’atmosphère). Lorsque ce signal est à 1, il n’y a pas d’aéorosls, et lorsque celui-ci est supérieur à 1, c’est qu’il y a présence d’aérosols. En dessous de 1, il y a un problème dans la mesure.
Exemple de backscater ratio
Le traitement ozone
On part de l’équation lidar pour un signal Rayleigh :
avec :
\(P\) : le signal reçu par le détecteur
\(\lambda_e\) : longueur d’onde d’emission du laser
\(\lambda_r\) : longueur d’onde de reception du détecteur
\(z\) : altitude
\(K\) : constante instrumentale
\(A\) : angle solide
\(beta\) : coefficient de diffusion
\(\tau_{i}\) : transmission de la moléculue i, les molécules principales de l’atmosphère (N2 et O2), les aérosols et les autres constituants minoritaires de l’atmosphère (ozone, SO2, etc)
La transmission est l’intégration de l’extinction en fonction de l’altitude, et s’écrit donc ainsi :
Soit le signal corrigé de l’altitude \(Pz2\) :
On dérive puis on prend le logarithme pour obtenir l’équation suivante :
En prenant le signal à deux longueurs d’ondes \(\lambda_{on}\) et \(\lambda_{off}\), on a alors deux équations :
Afin de simplifier l’équation suivante, nous posons :
On soustrayant les deux équations (1), on obtient alors :
On peut également écrire \(\alpha\) de la façon suivante :
avec :
\(\sigma_i\) : section efficace du constituant atmosphérique i
\(n_i\) : concentration du constituant atmosphérique i
On peut réécrire une partie de la somme de la façon suivante :
Si on pose :
L’équation devient alors :
L’équation principale devient donc :
D’où la concentration d’ozone peut s’écrire de la façon suivante :
Exemple de profil d’ozone
Le traitement température
Fonctionnement de la fonction lidar_actris_v01
Lancement de la fonction
La fonction lidar_ozone_actris_v01 accepte deux arguments : le fichier paramètre et la date à traiter (optionnel). La date est optionnelle, car elle peut se trouver dans le fichier de paramètre également. On le verra plus tard, mais cette fonctionnalité permet de faire tourner le programme avec un fichier paramètre par station.
Le fichier paramètre est une chaine de caractère, correspondant au chemin complet du fichier paramètre à utiliser. On peut également mettre en entrée un objet de classe ParametersFromXml (plus d’information concernant cette classe ici : Parameter file in xml).
def lidar_ozone_actris_v01(path_parameters_file: str | ParametersFromXml | None = None,
dt: str | date | None = None) -> Lidar | None:
"""
Data process for DIAL ozone lidar, ACTRIS.
:param path_parameters_file: path of the parameters file (default: None)
:type path_parameters_file: str or ParametersFromXml or None
:param dt: date to process (default: None)
:type dt: str or datetime.date or None
:return: the Lidar object
:rtype: Lidar or None
"""
Initialisation
Cette partie du code permet de retourner la valeur None si la varaible parameter en entrée est égale à None. Si le programme renvoie un None et non une erreure, cela permet de pouvoir relancer le programme de façon automatique, comme on le verra plus tard.
########
# Init #
########
if path_parameters_file is None:
return None
Creation de l’objet Lidar
Tout le programme repose sur la classe Lidar. La première étape est donc de l’initialiser. Elle a en entrée les deux paramètres d’entrée de la fonction lidar_actris_v01. C’est à partir de cet objet que sera stockées les différentes données (lidars, auxilières, traitements, etc), ainsi que toutes les fonctions de traitement.
A noter que les données sont stockées selon le formalisme NetCDF.
#####################
# 0. Parameter File #
#####################
lidar = Lidar(path_parameters_file, dt)
Si la date est l’un des paramètres, le programme va enregistrer cette date dans le fichier xml, à l’endroit suivant :
<Dial>
<Read>
<date_to_process>
<date type="date" format_date="iso"></date>
</date_to_process>
</Read>
</Dial>
Le format iso est le suivant : YYYY-MM-DD. Exemple : 2023-03-28 pour le 28 mars 2023.
Importation des données lidars
Première étape, il faut importer les données lidars dans l’objet Lidar. Ce sont ces données qui vont créées les coordonnées du fichier NetCDF : altitude et temps.
#############################
# 1. Reading Raw lidar Data #
#############################
lidar.import_lidar_data()
Pour paraméter cette partie, il faut remplir les champs suivants dans le fichier xml :
<Dial>
<Read>
<Lidar_files>
<filenames type="relativepath">pathname_lidar</filenames>
<format>Teslas</format>
<Signal_type type="int">
102, 103, 100, 101, 104, 105
</Signal_type>
</Lidar_files>
</Read>
</Dial>
Dans ce cas, le programme va chercher à lire le fichier pathname_lidar. Il est nécessaire de rajouter le format du fichier (ici, Teslas), ainsi que les signaux types des différentes voies. Les valeurs de ces différents Signal_Type sont présenté ici : Signal_Type.
La section filenames peut également contenir plusieurs fichiers, s’il y a plusieurs fichiers à lire. Il suffit de séparer les différents noms de fichier par un délimiteur (par defaut : “,”, mais peut être remplacé par une option dans le fichier xml, comme ceci :
<Dial>
<Read>
<Lidar_files>
<filenames type="relativepath" delimiter=";">
pathname_lidar_1;pathname_lidar_2
</filenames>
</Lidar_files>
</Read>
</Dial>
Dans le cas d’un fichier dépendent d’une hiérarchie de classement dépendant de la date, le champs filenames peut être décomposé en 4 sous-parties : le répertoire de base, un formalisme de répertoire avec la date à l’intérnieur (au format datetime de Python), le prefixe du nom de fichier, le suffixe du nom de fichier, et au milieu, le format de la date dans le nom du fichier. Dans l’exemple suivant, le programme va chercher tous les fichiers dont la structure est la suivante : Configuration_files/Main_files/Lidar/sho240328*.s*. En effet, dans notre cas, la partie date dans le répertoire est vide, donc non utilisé, et le format de la date utilisé dans le fichier est l’année sur 2 digits, le mois sur 2 digits et le jour sur 2 digits.
<Dial>
<Read>
<date_to_process>
<date type="date" format_date="iso">2024-03-28</date>
</date_to_process>
<Lidar_files>
<filenames type="relativepath">
<repertory check="repertory" type="relativepath">
Configuration_files/Main_files/Lidar
</repertory>
<repertory_date_format></repertory_date_format>
<data_name_prefixe>sho</data_name_prefixe>
<data_name_suffixe>*.s*</data_name_suffixe>
<data_name_date_format>%y%m%d</data_name_date_format>
</filenames>
</Lidar_files>
</Read>
</Dial>
On peut ajouter l’unité de la donnée dans le fichier d’entrée : il y a deux choix possibles, “photocouting” (par défaut) ou “MHz”. On le rajoute de la façon suivante :
<Dial>
<Read>
<unit_photocounting>MHz</unit_photocounting>
</Read>
</Dial>
On peut également décider de ne pas sauvegarder toutes les données issues du fichier d’entrée. Pour cela, on va utiliser l’option “Remove_Channel” : il suffit d’ajouter les Signal_Type qu’il ne faudra pas prendre en compte. Exemple :
<Dial>
<Read>
<Remove_channel type="int">101</Remove_channel>
</Read>
</Dial>
Lorsque le fichier brut contient des données avec des résolutions spatiales différente, il est nécessaire de réaliser une régression afin de mettre toutes les données dans la même résolution. Pour cela, on peut ajouter la ligne suivante. Dans ce cas, on force la résolution à 15 mètres.
<Dial>
<Read>
<ForceResolution type="float" units="m">15</ForceResolution>
</Read>
</Dial>
On peut ajouter un délai au niveau de l’électronique, pour le calcul du vecteur altitude. Pour cela, on utilise le champ electronic_delay. Ce délai peut être différent en fonction de la voie. Dans l’exmple suivant, c’est un délai de 100ns pour les 6 premières voies.
<Dial>
<Read>
<electronic_delay type="float" units="ns">100, 100, 100, 100, 100, 100</electronic_delay>
</Read>
</Dial>
Importation des données auxilières
Pour fonctionner, le programme a, au minimum, besoin d’un profil de pression et de température, en fonction de l’altitude (composition atmosphérique de l’atmosphère). Ajouter à cela, le programme a besoin pour de l’affichage des composants atmosphériques liées à l’ozone. Cela donne le code suivant :
############################
# 2. Import auxiliary data #
############################
lidar.import_atmospheric_component()
lidar.import_other_atmospheric_component("ozone")
Pour la partie climatologie de l’atmosphère, il faut remplir les champs Atmosphere_file du fichier xml. Le programme accepte plusieurs modèles possibles, qui seront joint en fonction de l’altitude. En effet, en régle général, on va utiliser un radiosondage dans les premières couches de l’atmosphère (jusqu’à 30km) qu’on va compléter avec un ou des modèles climatologiques. Il faut donc renseigner les informations de chacun de ces modèles. Pour cela, on va créer des objets modele_X, avec X un numéro (ce numéro n’est pas important, mais pour des soucis de lisibilités, il est préréfable de mettre le numéro d’ordre).
Pour chaque modèle, il faudra remplir :
le ou les nom des fichiers où se trouvent le ou les fichiers (en utilisant le même modèle que pour le nom du fichier lidar)
le format du ou des fichiers (radiosondage, ncep, arletty)
l’ordre d’importante pour la concaténation
la valeur par défaut
des variables liées au modèle (exemple : le nom de la station pour les fichiers Arletty
A noter que, pour le nom des fichiers, une nouvelle variable est apparut : time_shift. Il permet de décaler en temps la recherche du fichier. Par exmple, si on souhaite récupérer le radiosondage de la veille, on peut ajouter -1 comme valeur pour ce champ.
Voici un exemple de code xml :
<Dial>
<Read>
<Atmosphere_file>
<modele_0>
<filenames type="relativepath">
<repertory check="repertory" type="relativepath">
Configuration_files/Main_files/Atmosphere
</repertory>
<repertory_date_format></repertory_date_format>
<data_name_prefixe delimiter=",">ht, te</data_name_prefixe>
<data_name_suffixe>.nmc</data_name_suffixe>
<data_name_date_format>%y%m%d</data_name_date_format>
<time_shift type="float">0</time_shift>
</filenames>
<format>ncep</format>
<order type="int">1</order>
<fill_value>nan</fill_value>
<station_name>OHP</station_name>
</modele_0>
<modele_1>
<filenames type="relativepath">
Configuration_files/Main_files/Atmosphere/MAP85_PRE.DAT,
Configuration_files/Main_files/Atmosphere/MAP85_TEMn.DAT
</filenames>
<format>MAP85</format>
<order type="int">1</order>
<fill_value>nan</fill_value>
</modele_1>
</Atmosphere_file>
</Read>
</Dial>
Le programme a besoin des sections efficaces de l’ozone pour fonctionner. Pour les ajouter, il faut remplir le champ Cross_section_file du fichier xml :
<Dial>
<Read>
<Cross_section_file>
<filenames type="relativepath">
Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_218K.dat,
Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_228K.dat,
Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_243K.dat,
Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_273K.dat,
Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_295K.dat
</filenames>
<format>Bass And Paur files</format>
</Cross_section_file>
</Read>
</Dial>
Le programme a besoin d’un fichier climatologie de l’ozone (pour afficher des courbes). Il faut donc remplir le champ Climatology_ozone_file du fichier xml :
<Dial>
<Read>
<Climatology_ozone_file>
<filenames type="relativepath">
Configuration_files/Climatology_files/climatology_ozone_ohp.nc
</filenames>
<format>netcdf</format>
</Climatology_ozone_file>
</Read>
</Dial>
Création d’un signal lidar simulé
Pour les figures, et valider les données, on va comparer les données en fonction d’un signal lidar théorique, à chaque étape de traitement. Pour cela, il est nécessaire de le calculer. On utilise alors la commande suivante :
###########################
# Lidar Signal Simulation #
###########################
lidar.create_simulated_lidar()
Pré-traitement des données
Le programme est fait pour réaliser 3 prétraitements : calcul du fond de ciel, désaturation et correction du facteur géométrique. Pour appliquer ces différents paramètres, il faut utiliser la fonction apply_process de la classe Lidar. Cela donne le code suivant :
###############
# Pre-Process #
###############
lidar.background()
lidar.saturation()
lidar.geometrical_factor()
lidar.apply_preprocess()
Exemple de programme xml pour le paramètrage de ces fonctions :
<Dial>
<!-- DESATURATION -->
<Desaturation>
<Method>pelon, pelon, pelon, pelon, pelon, pelon</Method>
<Parameters>
<tau type="float" units="ns">3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9</tau>
<N_cm type="float">77.895,79.84,80.7487,70.1884,60,60</N_cm>
</Parameters>
</Desaturation>
<!-- SKY BACKGROUND -->
<SkyBackground>
<Method>Weighted linear regression, Weighted linear regression, Average, Average, Average, Average</Method>
<Parameters>
<alt_min units="m" separator=",">100220,92270,87620,92270,65870,60020</alt_min>
<alt_max units="m" separator=",">149870,138000,149870,149870,149870,149870</alt_max>
<order type="int">1,1,0,0,0,0</order>
<log_regression type="bool">False,False,False,False,False,False</log_regression>
</Parameters>
</SkyBackground>
<!-- GEOMETRIC FACTOR -->
<GeometricFactor>
<filenmes relativepath="True"></filenmes>
<format>mrd</format>
</GeometricFactor>
</Dial>
Traitement aérosols
Afin de pouvoir visualiser la présence ou non d’aérosols, et de pouvoir corriger des aérosols le profil ozone, le programme calcule la rétrodiffusion due aux aérosols.
####################
# Aerosols Process #
####################
lidar.aerosols_process()
Exemple de programme xml pour le paramètrage de ces fonctions :
<Dial>
<Aerosols>
<altitude_max type="float" units="km">40</altitude_max>
<altitude_ref type="float" units="km">25</altitude_ref>
<beta_aero_ref type="float" units="m**-1">0</beta_aero_ref>
<method>Klett</method>
</Aerosols>
</Dial>
Calcul du profil ozone par la méthode DIAL
Afin de calculer l’ozone, il est nécessaire de filtrer le signal (dans le cas de la méthode DIAL, en le dérivant), de liées les mesures hautes et basses, de calculer l’ozone à partir des signaux dérivés, puis de lier les voies Raymeigh et les voies Raman.
#################
# Ozone Process #
#################
lidar.filtering()
lidar.glue(name_param_tag="MergingSlope")
lidar.ozone_process()
lidar.glue()
Exemple de programme xml pour le paramètrage de ces fonctions :
<Dial>
<Filtering>
<parameters>
<filter_name>Savitzky-Golay</filter_name>
<mode>dial poly</mode>
<ncan1>3,3,4,4,5,5</ncan1>
<ican1>120,120,120,120,60,60</ican1>
<ncan2>35,35,30,30,40,40</ncan2>
<ican2>300,300,300,300,200,200</ican2>
<alt_unit>km</alt_unit>
</parameters>
<resolution_methodology options="df,ir,origin" multiple_option="True">df, ir, origin</resolution_methodology>
</Filtering>
<MergingSlope>
<channel_1 type='int' separator=",">100, 101</channel_1>
<channel_2 type='int' separator=",">102, 103</channel_2>
<method separator=",">linear, linear</method>
<shift type='int' separator=",">0, 0</shift>
<alt_min type="float" units="km" separator=",">17.12,17.12</alt_min>
<alt_max type="float" units="km" separator=",">18.02,18.02</alt_max>
</MergingSlope>
<MergingOzone>
<channel_1>132</channel_1>
<channel_2>136</channel_2>
<method>linear</method>
<shift>0</shift>
<alt_min units="km">9.02</alt_min>
<alt_max units="km">10.07</alt_max>
</MergingOzone>
<Ozone>
<binomial_filter type="int">2</binomial_filter>
</Ozone>
</Dial>
Domaine de validité
Calcul des colonnes partielles d’ozone
Ajout d’un facteur géométrique
Export données
Lancement du programme en bash
Afin de pouvoir lancer le programme avec un script bash, on programme a été crée : run_dial.py (Run Dial software)