Les produits de Niveau 3A

=>

Parmi les produits que nous préparons pour les traitements de données Sentinel-2 du pôle THEIA, les produits de Niveau 3A n'ont pas encore été abordés dans ce blog. Ceux-ci fournissent une synthèse mensuelle des produits de niveau 2A. Les produits de Niveau 3A devraient être bien utiles pour les raisons suivantes :

  • Le produit de Niveau 3A, fourni une fois par mois, représente un volume de données jusqu'à 6 fois inférieur à celui des niveaux 2A acquis pendant un mois.
  • Le produit de niveau 3 fournit un échantillonnage temporel régulier de l'évolution des réflectances, alors que l'échantillonnage du Niveau 2A est dépendant de la couverture nuageuse
  • Beaucoup de méthodes de traitement sont perturbées par la présence de trous dans les images dus à la couverture nuageuse. Le Niveau 3A vise à minimiser les pixels manquants.

 

Grâce au jeu de données SPOT4(Take5), nous avons pu expérimenter différentes méthodes de création des produits de Niveau 3A sur différents types de sites. C'est Mohamed Kadiri, au CESBIO, qui a pris en charge ce travail, financé par le budget CNES du projet THEIA, et avec le soutien de Mirelle Huc. La méthode consiste à calculer, pour chaque pixel, la moyenne des réflectances de surface des observations non nuageuses, obtenues pendant une demi période de N jours autour de la date T0 du produit de niveau 3A. Par exemple, dans l'exemple ci-dessous obtenu avec N=21, pour le produit du 15 mars (T0), les données utilisées couvrent la période du 24 février au 4 avril.

 

 

La moyenne calculée est une moyenne pondérée, qui attribue plus de poids :

  • aux images sans nuages
  • aux pixels situés loin de nuages
  • aux images avec peu d'aérosols
  • aux images proches de la date du produit de Niveau 3

 

Les poids et surtout la demi-période N ont une forte influence sur la qualité des résultats. Pour déterminer leurs valeurs, nous avons mis au point 3 critères de qualité :

  • le pourcentage de pixels dont la réflectance est manquante en raison des nuages
  • la différence du produit de Niveau 3A avec un produit de Niveau 2A faiblement nuageux acquis à quelques jours de la date centrale T0.
  • une mesure des artefacts. Les artefacts apparaissent en bordure des trous (dus aux nuages ou aux ombres) présents sur chacune des images.

 

Voici par exemple les performances obtenues sur le site de Versailles, site fortement nuageux lors du printemps 2013. Cette courbe confirme que malgré le mauvais  temps, Sentinel-2 devrait pouvoir fournir des produits quasiment sans nuages chaque mois sur ce genre de sites. Sur ce site, l'optimum de la durée de synthèse est compris entre 2*21 et 2*28 jours :

Performances obtenues sur le site de Versailles, pour plusieurs valeurs de la longueur de période N. En rouge, les trous résiduels (échelle à gauche), en jaune et vert, l'écart maximal à l'image centrale exprimé en réflectance, pour les 70% et 90% de pixels les meilleurs, et en bleu l'écart-type moyen des artefacts

 

Les produits de niveau 3A de Sentinel-2 devront aussi inclure une correction des effets directionnels, afin de pouvoir inclure dans un même produit de Niveau 3A des données issues de traces orbitales différentes, et donc acquises avec des angles de prise de vue différents. Enfin, en option, nous pourrons proposer une opération de bouchage des trous résiduels par interpolation temporelle ou en utilisant le comportement moyen de pixels similaires. Bref, il nous reste du pain sur la planche. Une comparaison avec la méthode classique de maximum de NDVI est fournie dans cet article.

The product level names, how they work ?

Simulation of Sentinel-2 products from Formosat-2 data (CESBIO)

 

Many users are confused by the earth observation product names. Maybe detailing the logic behind the names will help recall them. Here is how the THEIA Land Data Center product names were chosen.

  • we decided to use Sentinel-2 product names 1C, 2A, 3A, since we are sure Sentinel-2 will become a reference in high resolution earth observation.
    • Level 1C is a monodate ortho-rectified image expressed in TOA reflectance
    • Level 2A is a monodate ortho-rectified image expressed in surface reflectance, provided with a cloud/cloud shadow/snow/water mask
    • Level 3A is a monthly composite of Level2A Cloud/Cloud shadows free pixels
  • this nomenclature, defined by ESA and CNES, complies with the norms of the Committee on Earth Observation Satellites (CEOS)

 

CEOS naming rules are quite difficult to find, and I had searched them several times unsuccessfully. But recently, I found the list of members of the CEOS product harmonization committee, and two of its members (Frédéric Baret (INRA, France), and Kenneth McDonald (NOAA, USA)) replied very quickly to my questions.

 

The CEOS norm is based on a nomenclature defined by NASA in 1996, which is available on  Wikipedia :

 

Data Level NASA-EOSDIS Definition
Level 0 Reconstructed, unprocessed instrument and payload data at full resolution, with any and all communications artifacts (e.g., synchronization frames, communications headers, duplicate data) removed.
Level 1A Reconstructed, unprocessed instrument data at full resolution, time-referenced, and annotated with ancillary information, including radiometric and geometric calibration coefficients and georeferencing parameters (e.g., platform ephemeris) computed and appended but not applied to Level 0 data.
Level 1B Level 1A data that have been processed to sensor units (not all instruments have Level 1B source data).
Level 2 Derived geophysical variables at the same resolution and location as Level 1 source data.
Level 3 Variables mapped on uniform space-time grid scales, usually with some completeness and consistency.
Level 4 Model output or results from analyses of lower-level data (e.g., variables derived from multiple measurements).

 

I do not know what you think of it, but my sense is that it is quite vague in some aspects (what is a  "sensor unit") and too directive in some other aspects : a resampling of data on a cartographic grid is only allowed at level 3. If a uses does not want to handle the always complex reprocessing of data, he has to use the level 3 data.

 

The CEOS norm provided below was clearly inspired by NASA's norm, but it allows a data resampling  starting at level 1, and the data can be expressed in Physical units and not only "Sensor units". The CEOS norm does not detail the sub-levels (1A, 2B...). However, the distinction between Level 1 and Level 2 is still a bit fuzzy, as it is not always easy to tell a physical unit from a geophysical unit. We often consider a top of atmosphere reflectance as a Level 1 product and a surface reflectance after atmospheric correction a Level 2 product. Is a surface reflectance a physical unit or a geophysical unit?

 

Data Level CEOS Definition
Level 0 Reconstructed unprocessed instrument data at full space time resolution with all available supplemental information to be used in subsequent processing (e.g., ephemeris, health and safety) appended.
Level 1 Unpacked, reformatted level 0 data, with all supplemental information to be used in subsequent processing appended. Optional radiometric and geometric correction applied to produce parameters in physical units. Data generally presented as full time/space resolution. A wide variety of sub level products are possible.
Level 2 Retrieved environmental variables (e.g., ocean wave height, soil moisture, ice concentration) at the same resolution and location as the level 1 source data.
Level 3 Data or retrieved environmental variables which have been spatially and/or temporally re-sampled (i.e., derived from level 1 or 2 products). Such re-sampling may include averaging and compositing.
Level 4 Model output or results from analyses of lower level data (i.e., variables that are not directly measured by the instruments, but are derived from these measurements).

 

Having had some difficulties finding the famous CEOS norm, I had build my own idea of what the CEOS norm should be, probably, from dscussions with colleagues during the POLDER.project preparation. So here is my own version of the product level names norm (but I know a personal norm is not a norm anymore)

 

Data Level Product nomenclature according to... myself
Level 0 Reconstructed unprocessed instrument data at full space time resolution with all available supplemental information to be used in subsequent processing (e.g., ephemeris, health and safety) appended.
Level 1 All pixels were acquired at the same time (within a few instants, during one satellite overpass), and their processing does not make assumptions on the nature of the observed pixel. Each pixel is processed in the same way, whatever it is made of (cloud, forest, sea...). The values are expressed in physical units or the product provides all the necessary information to convert the values to physical units. The product may be resampled onto a cartographic grid, or may just provide the necessary information to resample it.
Level 2 All pixels were acquired all at the same time (within a few instants, during one satellite overpass), but here, the processing may include assumptions on the nature of the pixel, for instance concerning the atmosphere, vegetation of sea state. The pixels may be processed differently according to their nature.
Level 3 The product is obtained from data acquired at different dates, often with different footprints. As for Level 2, processing may differ according to the pixel nature, and assumptions on this nature are allowed. 

Level 3 products are often composite products based on the level 2 data acquired during a certain period of time (10 days, one month...)

 

The CEOS norm does not define sub-levels, and for that aspect, NASA's norm has some influence. For instance, with Sentinel-2, Level 1A and Level 1B exist as internal products and are quite similar to what is in NASA norm, while the Level 1C is ortho-rectified and expressed in top of atmosphere reflectance. The level 1C only will be distributed to standard users.

 

Sentinel-2 Mission requirement document also defines a level 2A, expressed in surface reflectance, and a Level 2B that provides biophysical variables such as LAI, fAPAR.

 

Finally, let's recall that several Earth observing missions do not respect the CEOS norms, often because their nomenclature was defined before the norm existed. It is the case of SPOT and followers (Pleiades, Rapid Eye...)

 

 

Les noms des niveaux de produits, comment ça marche ?

=>

Simulations des produits Sentinel-2 à partir d'images Formosat-2

 

Vous êtes nombreux à ne pas retenir les noms des produits d'observation de la terre, et si je n'étais pas tombé dedans quand j'étais petit, je ne ferais pas mieux.  Peut-être qu'en détaillant la logique qui a conduit à leur définition, il sera plus facile de les retenir. Voici donc comment les noms des produits du pôle THEIA ont été définis :

 

  • nous avons décidé d'utiliser la nomenclature de Sentinel-2 : 1C, 2A,  3A, car nous sommes convaincus que Sentinel-2 deviendra rapidement la référence en imagerie spatiale
    • Le produit de Niveau 1C est une image monodate ortho-rectifiée, exprimée en réflectance au sommet de l'atmosphère.
    • Le produit de Niveau 2A est une image monodate ortho-rectifiée, exprimée en réflectance de surface, accompagnée d'un masque de nuages/ombres de nuages/neige/eau
    • le produit de Niveau 3A est une synthèse mensuelle ortho-rectifiée des pixels non nuageux provenant des Niveau 2A.
  • cette nomenclature, définie par l'ESA et le CNES respecte les normes du Comittee on Earth Observation Satellites (CEOS)

 

Les normes du CEOS sont assez difficiles à trouver, je les avais cherchées à plusieurs reprises. Heureusement, j'ai fini par trouver quelques membres du comité travaillant sur l'interopérabilité des systèmes, et deux d'entre eux (Frédéric Baret (INRA), et Kenneth McDonald (NOAA)) ont répondu immédiatement à mes questions.

 

La norme du CEOS reprend en fait les termes d'une norme de la NASA, datant de 1986, qui elle est disponible sur Wikipedia :

 

Data Level NASA-EOSDIS Definition
Level 0 Reconstructed, unprocessed instrument and payload data at full resolution, with any and all communications artifacts (e.g., synchronization frames, communications headers, duplicate data) removed. (In most cases, the EOS Data and Operations System (EDOS) provides these data to the data centers as production data sets for processing by the Science Data Processing Segment (SDPS) or by a SIPS to produce higher-level products.)
Level 1A Reconstructed, unprocessed instrument data at full resolution, time-referenced, and annotated with ancillary information, including radiometric and geometric calibration coefficients and georeferencing parameters (e.g., platform ephemeris) computed and appended but not applied to Level 0 data.
Level 1B Level 1A data that have been processed to sensor units (not all instruments have Level 1B source data).
Level 2 Derived geophysical variables at the same resolution and location as Level 1 source data.
Level 3 Variables mapped on uniform space-time grid scales, usually with some completeness and consistency.
Level 4 Model output or results from analyses of lower-level data (e.g., variables derived from multiple measurements).

 

Je ne sais pas ce que vous en pensez, mais je trouve cette description assez floue, la notion de "sensor units" n'est pas claire, et surtout, il faut attendre le produit de niveau 3 pour que les données soient rééchantillonnées sur une grille cartographique. Si l'utilisateur ne veut pas avoir à se soucier de rééchantillonnage géométrique, il doit utiliser le produit de Niveau 3.

 

La norme du CEOS (ci dessous), clairement inspirée de celle de la NASA, permet de rééchantillonner les données dès le Niveau 1 et de les exprimer en grandeur physiques et non pas pas en "sensor units". La norme du CEOS ne détaille pas les sous niveaux (1A, 1B). Il me semble que cette norme ne détaille cependant pas bien la différence entre les grandeurs physiques et variables géophysiques (par exemple, pourquoi une réflectance au sommet de l'atmosphère serait -elle une grandeur physique et une réflectance de surface une variable biophysique), et laisse beaucoup de marge d’interprétation.

 

Data Level CEOS Definition
Level 0 Reconstructed unprocessed instrument data at full space time resolution with all available supplemental information to be used in subsequent processing (e.g., ephemeris, health and safety) appended.
Level 1 Unpacked, reformatted level 0 data, with all supplemental information to be used in subsequent processing appended. Optional radiometric and geometric correction applied to produce parameters in physical units. Data generally presented as full time/space resolution. A wide variety of sub level products are possible.
Level 2 Retrieved environmental variables (e.g., ocean wave height, soil moisture, ice concentration) at the same resolution and location as the level 1 source data.
Level 3 Data or retrieved environmental variables which have been spatially and/or temporally re-sampled (i.e., derived from level 1 or 2 products). Such re-sampling may include averaging and compositing.
Level 4 Model output or results from analyses of lower level data (i.e., variables that are not directly measured by the instruments, but are derived from these measurements).

Ayant eu bien du mal à trouver ces fameuses normes du CEOS, je m'en étais fait ma propre idée, probablement à partir de discussions avec les collègues qui travaillaient avec moi sur le projet POLDER. Voici ma version personnelle de la norme (sachant qu'une norme personnelle ne sert à rien bien sûr...)

 

Data Level La nomenclature des produits selon OH
Level 0 Données reconstruites mais non traitées, fournies à pleine résolution, accompagnées de toutes les informations nécessaires pour les traitements de niveaux supérieurs.
Level 1 Le produit de Niveau 1 est un produit dont tous les pixels ont été acquis à la même date (en quelques instants, au cours d'un passage du satellite), et dont le traitement ne fait pas d'hypothèses sur la nature du paysage observé. Chaque pixel de l'image est traité de la même manière, quelle que soit sa nature (nuage, forêt, mer...). Le produit exprime les données en grandeur physique ou fournit toutes les informations nécessaires pour les convertir. Il peut être rééchantillonné sur une grille cartographique ou doit fournir toutes les informations nécessaires à cette conversion.
Level 2 Le produit de Niveau 2 doit être lui aussi acquis en quelques instants, mais cette fois, des traitements différenciés par type de pixel sont autorisés, et l'on peut faire des hypothèses sur l'atmosphère, l'état de la végétation ou de la mer, et traiter les différentes classes de pixels d'une manière différente.
Level 3 Le produit de Niveau 3 est constitué de données acquises à des dates différentes, souvent avec des emprises différentes, et, comme sur le Niveau 2, les traitements peuvent faire des hypothèses sur ce qui est observé. Les produits de niveau 3 sont en général des synthèses périodiques (hebdomadaires, décadaires, mensuelles) des produits de Niveau 2.

 

 

La norme du CEOS ne définit les sous-niveaux de produits, mais dans ce domaine, la norme de la NASA a une certaine influence. Pour Sentinel-2, les produits de niveau 1A et 1B sont des produits internes assez proche de la norme de la NASA. mais le produit de niveau 1C a fournit des données ortho-rectifiées dans une projection cartographique. Seul le niveau 1C devrait être distribué aux utilisateurs normaux.

 

 

La spécification de mission de Sentinel-2 définit également un produit de Niveau 2A, exprimé en réflectance de surface après correction atmosphérique, et un produit de niveau 2B fournissant des variables biophysiques (LAI, FAPAR...). Aucun de ces produits ne sera cependant fourni par le segment sol de Sentinel-2. Ces produits existeront en revanche dans le centre de traitement de THEIA, et nous les avons déja produits pour SPOT4 (Take5).

 

Enfin, notons que de nombreux produits d'observation de la terre ne respectent pas ces normes, souvent parce que la nomenclature a été définie avant l'adoption de la norme : c'est le cas de SPOT et de ses successeurs, y compris les plus récents (Pleiades, Rapid Eye).

 

The terrain effect correction : how it works

=>


Caution, this post contains formulas.

 

The topographic (or terrain) effects on the observed reflectances are due to several phenomena, illustrated below :

  • the closer the surface is perpendicular to sun direction, the more energy it receives per surface unit (we talk about irradiance). If the surface is parallel to sun direction, it does not receive direct sunlight. We can model it this way :
    • For an horizontal surface :  E_h= E_0.T_{dir}^\downarrow.cos(\theta_s)

      Definition of sun zenith angle and sun incidence angle

    • For a sloped surface  E_i=E_0.T_{dir}^\downarrow.cos(\theta_i)
    •  E_0 is the Top of Atmosphere irradiance, and  T_{dir}^\downarrow is the downward direct transmission, i.e. the proportion of the light that reaches directly the surface without being scattered by the atmosphere.
    • Assuming that all the irradiance is direct, the measured reflectance if the surface was horizontal is calculated from the following formula:  \rho_h=\rho_i \frac{cos(\theta_s)}{cos(\theta_i)} . However, the above assumption is not true and this formula tends to over correct terrain effects
  • The surfaces also receive a diffuse sun irradiance scattered by the atmosphere. If the surface is not horizontal, a part of the sky is obscured by the slope reducing the diffuse irradiance. Moreover, the diffuse irradiance depends on the amount of aerosols (and clouds) in the atmosphere. In addition, the surrounding terrain can also hide a part of the sky, but we do not take this effect into account here in our modelling. We use the following approximation, which is equivalent to assuming that the slope is alone in a horizontal region.
    • If surface is horizontal, the visible sky fraction is 1, if it is vertical, this fraction is 1/2
    •  \displaystyle F_{sky}= \frac{1+cos(slope)}{2}
  • Finally, the slope can receive light from surrounding surfaces, which become directly visible. In our simplified model, we always assume the entire environment of our slope is flat and,for instance, we do not take the effect of the opposite side in a valley into account :
    • If surface is horizontal, the visible ground fraction is 0, if it is vertical, it is 1/2.
    •  \displaystyle F_{fround}= \frac{1-cos(slope)}{2}

 

Finally, we use the following formula  to compute the reflectance that would be observed if the surface was horizontal   \rho_{h} , as a function of the slope (inclined) reflectance   \rho_{i} :

  \rho_{h}=\displaystyle \rho_{i}.\frac{T^{\downarrow}}{T_{dir}^{\downarrow}.\frac{cos(\theta_i)}{cos(\theta_s)} + T_{dif}^{\downarrow} F_{sky} + T^{\downarrow} F_{ground} \rho_{env}}

 T^{\downarrow} is the downward transmission, sum of direct and diffuse irradiances :  T^{\downarrow}= T_{dir}^{\downarrow}+ T_{dif}^{\downarrow} , and \rho_{env} is the average reflectance of the neighbourhood.

 

Finally, we can also account for bidirectional reflectance effects, but this correction is tricky since directional effects depend on the surface cover type. See for instance : Dymond, J.R.; Shepherd, J.D. 1999: Correction of the topographic effect in remote sensing. IEEE Trans. Geosci. Remote Sens. 37(5): 2618-2620.

 

It is very difficult to validate the correction of directional effects : , we could compare the correction results for satellite overpasses at different times in the day. But all the satellite optical imagers have nearly the same overpass time. A qualitative way of estimating the accuracy is to check that similar land covers on opposite slopes in a valley ( a meadow, a forest) have a similar reflectance after correction. The most suitable points are North-South valleys. Here are some examples of terrain effects correction results.

Formosat-2 image in the Alps, before (left) and after (right) terrain effect correction. The image on the right seems to have been flattened, and the opposite slope reflectances seem much more alike after correction.

Finally, an essential part of the method's accuracy is the availability of a highly accurate digital  elevation model (DEM),  up to now, only the SRTM DEM is available globally, and it only has a 90 meter resolution. Its accuracy is somewhat inadequate and sometimes leaves artefacts if the slope changes are poorly located.

 

 

La correction des variations d'éclairement dues au relief

=>


Attention, cet article contient des formules.

 

Les effets du relief sur les réflectances observées proviennent de plusieurs phénomènes, schématisés dans l'illustration ci-dessous.

  • plus la surface est perpendiculaire à la direction solaire, plus elle va recevoir d'énergie par unité de surface (on parle d'éclairement). Si la surface est parallèle à la direction solaire, elle ne reçoit pas de lumière directe du soleil, l'éclairement direct reçu est donc nul.
    • Pour une surface horizontale :  E_h= E_0.T_{dir}^\downarrow.cos(\theta_s)

      Définition de l'angle zenithal solaire et de l'angle d'incidence

    • Pour une surface inclinée  E_i=E_0.T_{dir}^\downarrow.cos(\theta_i)
    •  E_0 est l'éclairement solaire en haut de l'atmosphère, et  T_{dir}^\downarrow est la transmission directe descendante, c'est à dire la proportion de cet éclairement qui rejoint directement la surface sans être diffusé par l'atmosphère.
    • Si on suppose que tout l'éclairement est direct, la réflectance mesurée si la surface avait été horizontale se calcule à partir de la formule suivante :  \rho_h=\rho_i \frac{cos(\theta_s)}{cos(\theta_i)} . Cette formulation a tendance à surcorriger l'effet des pentes.
  • Les surfaces reçoivent aussi un éclairement solaire diffusé par l'atmosphère. Si la surface n'est pas horizontale, une partie du ciel est masquée par la pente: l'éclairement diffusé par l'atmosphère et qui parvient à la surface est donc réduit. Cet éclairement va varier en fonction de la quantité d'aérosols (et de nuages) présents dans l'atmosphère. Par ailleurs, les reliefs environnants peuvent aussi masquer une partie du ciel, mais nous ne prenons pas cet effet là en compte dans notre modélisation. Nous utilisons l'approximation suivante, qui revient à supposer que la pente est le seul endroit incliné dans un environnement horizontal.
    • Si la surface est horizontale, la fraction de ciel Fciel vue est 1, si la surface est verticale, la fraction de ciel vue est 1/2
    •  \displaystyle Fciel= \frac{1+cos(pente)}{2}
  • enfin, la pente peut recevoir de la lumière des surfaces environnantes, qui deviennent directement visibles. Dans notre modélisation simplifiée, nous supposons toujours que tout l'environnement de notre pente est plat et nous ne prenons pas en compte l'effet du versant opposé dans une vallée.
    • Si la surface est horizontale, la fraction de sol vue est 0, si la surface est verticale, la fraction de sol vue est 1/2.
    •  \displaystyle Fsol= \frac{1-cos(pente)}{2}

 

Au final, pour calculer la réflectance mesurée si la surface avait été horizontale   \rho_{h} , en fonction de la réflectance mesurée sur la surface inclinée   \rho_{i} , on obtient donc la formule suivante :

  \rho_{h}=\displaystyle \rho_{i}.\frac{T^{\downarrow}}{T_{dir}^{\downarrow}.\frac{cos(\theta_i)}{cos(\theta_s)} + T_{dif}^{\downarrow} F_{ciel} + T^{\downarrow} F_{sol} \rho_{env}}

 T^{\downarrow} est la transmission descendante composée d'une partie directe et diffuse :  T^{\downarrow}= T_{dir}^{\downarrow}+ T_{dif}^{\downarrow} , et \rho_{env} est la réflectance moyenne autour du pixel à corriger.

 

Enfin, on peut aussi prendre en compte un facteur de correction lié aux effets directionnels, mais cet aspect est délicat, car les effets directionnels varient en fonction du type de surface observé. Voir par exemple  : Dymond, J.R.; Shepherd, J.D. 1999: Correction of the topographic effect in remote sensing. IEEE Trans. Geosci. Remote Sens. 37(5): 2618-2620.

 

Il est très difficile de valider cette correction : pour bien le faire, il faudrait comparer des images acquises à différentes heures de la journée. Malheureusement, tous les satellites à haute résolution passent à peu près à la même heure. On peut essayer cependant de vérifier que des surfaces ayant un couvert semblable (une prairie, une forêt) sur des versants opposés, présentent bien la même réflectance après correction. Les points les plus propices à ce genre de vérifications sont des vallées Nord Sud. Voici quelques exemples de résultats après correction.

Image Formosat-2 acquise dans les alpes, avant (à gauche) et après (à droite) correction des effets du relief. L'image de droite donne l'impression d'avoir été aplatie, les versants de part et d'autre des vallées présentent des réflectances bien plus semblables après correction.

Enfin, un élément essentiel de la précision de cette méthode est la disponibilité d'un modèle numérique de terrain très précis, pour le moment, le seul MNT disponible globalement est SRTM à 90 mètres de résolution. Sa précision est un peu insuffisante et laisse parfois des artefacts, si les ruptures de pentes sont mal localisées.

 

The adjacency effects, how they work.

As explained in the post about atmospheric effects, the scattering of light by molecules and aerosols in the atmosphere brings about several effects : scattering adds some haze on the images (the atmospheric reflectance), lessens the signal from the surface (the atmospheric transmission), and blurs the images (the adjacency effects). This post is about the adjacency effects, the other aspects have already been quickly explained in the above post.

 

The figure on the right shows the types of paths that light can follow before getting to the satellite. Path 1 corresponds to the atmospheric reflectance, path 2 is path that interacts with the target, it is the one which is useful to determine the surface reflectance, paths 3 and 4 contribute to the total reflectance but interact with the surface away from the target. These paths are thus the cause of adjacency effects and they blur the images.

 

 

If not corrected, adjacency effects may cause large errors. Let's take the case of a fully developed irrigated field surrounded by bare soil. For such a case, the second figure on the right shows the relative percentage of errors for reflectances and NDVI as a function of aerosol optical thickness, if adjacency effect is not corrected.

 

 

 

An approximate correction can be applied, but it thus requires to know the aerosol optical thickness. In our MACCS processor, here is how it works :

 

  1. We first correct the images under the assumption that the Landscape is uniform. We obtain a surface reflectance under uniform absorption which is noted  \rho_{s,unif} .
  2. We compute the neighbourhood reflectance (  \rho_{s,adj} ) using a convolution filter with a 2km radius, that computes the average neighborhood reflectance weighted by the distance to the target. To be fully rigorous, this filter should depend on the optical thickness and on the viewing and sun angle (The less aerosols, the larger radius), but as we did not work on an accurate model, we used a constant radius.
  3. We correct for the contribution of paths 3 and 4 using :

 \rho_{s}=\frac{\rho_{s,unif}.T^{\uparrow}.\frac{1-\rho_{s,unif}.s}{1-\rho_{s,adj}.s}-\rho_{s,adj}. T_{dif}^{\uparrow}}{T_{dir}^{\uparrow}}

  • where  T^{\uparrow}=T_{dif}^{\uparrow}+T_{dir}^{\uparrow} is the total upward transmission, sum of diffuse and direct upward transmissions, and s is the atmosphere spheric albedo. These quantities depend on the wavelength, on the aerosol model and on the AOT. They are computed using Look up Tables based on radiative transfer calculations.

 

As this processing uses convolution with a large radius, it takes quite a large part of the atmospheric processing time.

 

Result Exemples

The images below show 3 stages of the atmospheric processing, for 2 Formosat-2 images obtained over Montreal (Canada) with a 2 days interval. The first image was acquired on a hazy day (aerosol optical thickness (AOT) of 0.47 according to MACCS estimate); and the second one on a clear day (AOT=0.1).

  • The first line corresponds to the Top Of Atmosphere images, without atmospheric correction. The left image is obviously blurred compared to the right image.
  • The second line corresponds to the atmospheric correction under uniform landscape assumption (as in step 1). The left image is still obviously blurred compared to the right image.
  • the third line show the same images after adjacency effect correction. In that case, the left image is not blurred any more, it is even maybe a little over corrected as it seems somewhat sharper that the right image.

TOA Images (On the left, the hazy image)


Surface reflectance under uniform landscape assumption (on the left, the hazy image)

 

Surface reflectance after adjacency effect correction (on the left, the hazy image)

 

The pixel wise comparison of reflectances is also a way to show the enhancement due to the adjacency effect correction. The plot below compares the images of both dates corrected under the uniform landscape assumption (on the left), and after adjacency effect correction (on the right). You may observe that the dots are closer the the black diagonal on the right. On the hazy image (May 27th), the high reflectances are a little too low, while the low reflectances are a little too high, which is the symptom of a loss of contrast.

Les effets d'environnement, comment ça marche ?

=>

Comme expliqué dans l'article sur les effets atmosphériques, la diffusion de la lumière par les molécules et les aérosols présents dans l'atmosphère provoque plusieurs effets. La diffusion ajoute un voile aux données (la réflectance atmosphérique), atténue le signal en provenance de la surface (la transmission atmosphérique), et rend les images floues (les effets d'environnement). Cet article s'intéresse aux effets d'environnement, les autres aspects ont été abordés dans le lien fourni ci-dessus.

 

Le schéma ci-joint montre les différents types de trajets que peut suivre la lumière avant d'arriver au capteur. Le Trajet 1 correspond à la réflectance atmosphérique, le trajet 2 est proportionnel à la réflectance de la cible observée atténué par sa traversée de l'atmosphère, c'est celui qui nous intéresse et nous permet de retrouver la réflectance de surface. Les trajets 3 et 4 apportent au capteur une part de signal qui ne provient pas directement de la surface que le satellite observe mais de son voisinage (d'où le nom d'"effets d'environnement"). Ce sont ces trajets qui apportent du flou sur l'image.

 

Les effets d'environnement peuvent engendrer de fortes erreurs lorsqu'on observe une parcelle de végétation entourée de sols nus ou de végétation senescente. Pour un tel cas, la figure ci-dessous présente les erreurs en pourcentage de réflectance et de NDVI, si on ne prend pas en compte les effets d'environnement, en fonction de l'épaisseur optique des aérosols.

 

Il est possible de corriger ces effets de manière approchée, à condition de connaître la quantité d'aérosols. Dans les traitements de la chaîne MACCS, nous procédons de la manière suivante :

 

  1. Nous procédons à la correction atmosphérique en supposant que le paysage est uniforme. Nous obtenons une réflectance de surface sous hypothèse uniforme que nous notons  \rho_{s,unif} .
  2. Nous calculons la réflectance de l'environnement du pixel (  \rho_{s,env} ) en utilisant un filtre de convolution gaussien de 2 km de diamètre, qui calcule une moyenne pondérée de la réflectance environnante. En toute rigueur, ce filtre devrait dépendre de la quantité d'aérosols présents dans l'atmosphère (moins il y a d'aérosols, plus le rayon devrait être grand), et des angles de prise de vue, mais nous n'avons pas encore travaillé sur cet aspect, nous avons donc utilisé un filtre constant.
  3. Nous corrigeons finalement la réflectance des trajets 3 et 4 par la formule suivante :

 \rho_{s}=\frac{\rho_{s,unif}.T^{\uparrow}.\frac{1-\rho_{s,unif}.s}{1-\rho_{s,env}.s}-\rho_{s,env}. T_{dif}^{\uparrow}}{T_{dir}^{\uparrow}}

  •  T^{\uparrow}=T_{dif}^{\uparrow}+T_{dir}^{\uparrow} est la transmission atmosphérique montante totale, somme de la transmission atmosphérique diffuse et directe. s est l'albedo atmosphérique. Toutes ces grandeurs sont déduites de calculs de transfert radiatif et dépendent de la quantité et du type d'aérosols.

 

Cette correction qui implique l'utilisation de convolutions est assez lourde et prend près d'un quart du temps de correction atmosphérique.

 

Exemples de résultats

Les images ci dessous présentent 3 stades de la correction atmosphérique pour deux images Formosat-2 acquises au dessus du Canada, à deux jours d'intervalle, la première image est acquise un jour il y a beaucoup d'aérosols (épaisseur optique de 0.47 d'après nos calculs), alors que la seconde est acquise un jour très clair (épaisseur optique de 0.1 selon nos calculs).

 

  • La première ligne correspond aux images au sommet de l'atmosphère, sans correction atmosphérique. On voit bien que l'image de gaucheest plus floue.
  • La deuxième ligne correspond aux images corrigées en supposant le paysage uniforme. Il s'agit de l'image obtenue à l'issue de l'étape 1 dans la méthode décrite ci-dessus. L'image de gauche est toujours plus floue.
  • La troisième ligne présente ces mêmes images après la correction d'environnement. Dans ce cas, l'image de gauche n'est plus floue, elles est même légèrement trop nette (un peu de sur correction).

Images TOA (à gauche, l'image avec fort contenu en aérosols)

 

Images en réflectance de surface, en supposant le paysage uniforme (à gauche, l'image avec fort contenu en aérosols)

 

Images en réflectance de surface après correction des effets d'environnement (à gauche, l'image avec fort contenu en aérosols)

 

On peut aussi comparer point à point les réflectances pour juger de l'amélioration après correction des effets d'environnement. La courbe ci-dessous compare les images corrigées en supposant le paysage uniforme, et les images corrigées en tenant compte des effets d'environnement. On constate que les points se rapprochent de la diagonale après correction des effets d'environnement. Sur l'image du 27 mai, pour laquelle l'épaisseur optique est la plus forte, on note que les fortes réflectances sont un peu trop faibles, alors que les faibles réflectances sont un peu trop fortes, ce qui correspond bien à une perte de contraste.

How to estimate Aerosol Optical Thickness

=>

Caution ! This post contains formulas !


Aerosols play a great role in the atmospheric effects. Aerosols are particles suspended in the atmosphere, which can be of several types: sand or dust, soot from combustion, sulfates or sea salt, surrounded by water... Their size ranges between 0.1 micron and a few microns, depending on the type of aerosol or on the air moisture. Their quantity is also extremely variable : rain can suddenly reduce their abundance (known as "aerosol optical thickness"). The abundance variations result in great variations of observable reflectances from one day to the next, and it is therefore necessary to know the quantity and type of aerosols, in order to correct their effects.

 

Unfortunately, to correct the effects of aerosols, there is no global aerosol observation network, and the only available data are local observations from the few hundred points of Aeronet network. Therefore, this network can not be used operationally to correct the satellite images over large areas.

Weather forecast models just start predicting the amounts of aerosols, based on satellite observations and modeling of sources and sinks and of the transport of aerosols by the winds, but these data do not seem to have sufficient accuracy yet to be used for the atmospheric correction of images.

 

Our atmospheric correction method, named MACCS, is therefore based on an estimate of aerosol optical depth from the images themselves. To understand how this method works, one must already understand the effects of aerosols on radiation. We have seen in this post, that the effects of diffusion can be modelled as follows (assuming the corrected gas absorption):

ρTOA = ρatm +Td ρsurf

The reflectance at the top of the atmosphere ρTOA (Top of Atmosphere) is the sum of the atmospheric reflectance  ρatm and of the surface reflectance ρsurf transmitted by the atmosphere. We seek to know the surface reflectance, but for each measurement made at the top of the atmosphere, there are three unknowns to be determined. To separate the effects of the atmosphere and surface effects, we must use other information.

 

Dark pixel method

When the image includes a surface whose surface reflectance is nearly zero, the reflectance observed at the top of the atmosphere becomes ρTOA = ρatm. We can therefore deduce the atmospheric reflectance and using a radiative transfer model, the aerosols optical thickness (AOT). Finally, knowing the AOT, we can compute the diffuse transmission, and finally calculate ρsurf. An even simpler and more approximate version of this method consists in subtracting directly the reflectance of the dark pixel (or ρatm) to the entire image (neglecting the transmission) [Chavez, 1988].

 

However, this method assumes that there is a very dark area in the image (which is not always the case), and that the reflectance of the dark surface is known. The method also assumes that the amount of aerosols is constant over the image and it neglects the effect of terrain. The results obtained by this method can be quite inaccurate. In our method (MACCS), however, we use the method of black pixel determine the maximum value of the optical thickness in the area.

 

Multi Spectral Method, called "DDV"

If you know the type of aerosols in the atmosphere, it is possible to deduce the properties of aerosols in a spectral band from the optical properties in another spectral band.

 

If there are two spectral bands, there are two measures ρsurf and three unknowns (both surface reflectance in these bands, and the amount of aerosols). An additional equation can be obtained if we know the relationship between the surface reflectance of the two bands.

 

The method named "Dark Dense Vegetation" (DDV) is based on assumptions about relationships between surface reflectances of the dense vegetation exploiting the fact that the spectrum of dense green vegetation is quite constant. The most famous version of this method is that used by NASA for MODIS project [Remer 2005]. It connects the surface reflectance in the blue and red with those in the SWIR. This provides two equations for estimating the type of aerosol optical thickness. This method works well in temperate and boreal zones, but not in arid areas where it is difficult to find the dense vegetation. Early versions used the following equations:

 

ρBlue = 0.5 ∗ ρSWIR

ρRed = 0.25 ∗ ρSWIR

 

The following versions of the MODIS DDV algorithm are a bit more complicated but follow the same principle. Our work has shown that using the equation below allows a more accurate determination of the optical thickness, for less dense vegetation cover (NDVI to a 0.2) because bare soil brown also respect this relationship.

 

ρBlue = 0.5 * ρRed

(the exact value of the coefficient is adjusted according to the spectral bands of the instrument)


This version of the  method, however, does not allow to determine the aerosol model. In the case of SPOT4 (Take5), the absence of a blue band does not allow us to use this equation, resulting in a slight loss in accuracy.

This diagram shows that the correlation between surface reflectance above vegetation is much better for the (blue, red) couple of spectral bands than for couples including using (SWIR).

 

 

 

Multi Temporal Method

In most cases, the reflectance of the land surface changes slowly over time, while the aerosol optical properties vary rapidly from one day to another. We can therefore consider what changes from one image to another (apart from special cases often linked to human intervention) is associated with aerosols, and deduce the properties of aerosols and then correct for atmospheric effects. This method is too complex to be explained in detail here, interested readers can refer to [Hagolle 2008].

 

So that surface reflectance be nearly constant from one image to another, however, it is required that images be acquired at a constant angle. Indeed, the reflectance depend on the viewing angles: this is what we call directional effects. This method therefore applies only to satellite observations obtained with constant angle. It does not apply to standard SPOT data, but this condition is true for SPOT4 (Take5) data. It also applies to Landsat Venμs and Sentinel-2.

 

Finally :

 

Validation of aerosol optical thickness (AOT) from time-series of FORMOSAT-2 images, depending on the method (multi-spectral, multi-temporal, combined), compared with the measurements provided by the Aeronet network of in-situ measurements. The multi-spectral method works best on sites covered with vegetation and is much less accurate on arid sites, while the multi-temporal method performs a little worse on green sites, but much better on dry sites. The combination of the two methods retains the best of the two basic methods.

The MACCS/MAJA method, used for SPOT4 (Take5) experiment, and also for LANDSAT, VENμS and Sentinel-2 data, combines the three methods described above to obtain robust estimates of aerosol optical thickness. These methods work in many cases, but sometimes fail when the assumptions on which they are based prove to be incorrect. They generally tend to work better on vegetated areas rather than in arid areas. for now, they assume the model known aerosol and in the coming years, we will look for reliable ways to identify the type of aerosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENμS and Sentinel-2 Images. Remote Sens. 2015, 7, 2668-2691.

Les aérosols jouent un rôle prépondérant dans les effets atmosphériques. Les aérosols sont des particules en suspension dans l'atmosphère, qui peuvent être de plusieurs types : grains de sable ou poussières, suies issues de combustion, sulfates ou sels marins entourés d'eau... Leur taille peut varier de 0.1 µm à quelques microns, en fonction du type d'aérosols ou de l'humidité de l'air. Quant à leur quantité, elle est extrêmement variable, une pluie pouvant réduire brutalement leur abondance (on parle d'"épaisseur optique d'aérosols"). Ils peuvent faire varier fortement d'un jour à l'autre les réflectances observables depuis le sommet de l'atmosphère et il est donc nécessaire de connaître leur quantité et leur type afin de pouvoir corriger leurs effets.

 

Malheureusement, pour corriger les effets des aérosols, on ne dispose pas de réseau global d'observation des aérosols, seulement d'observations locales, sur les quelques centaines de points du réseau Aeronet. Ce réseau ne peut donc pas être utilisé pour corriger opérationnellement les images de satellites sur de grandes étendues.
Des modèles météorologiques commencent à prédire les quantités d'aérosols, en se basant sur les observations de satellites et la modélisation des sources et du transport des aérosols par les vents, mais ces données ne semblent pas encore avoir une précision suffisante pour être utilisées pour la correction atmosphérique des images.

 

Notre méthode de correction atmosphérique (MACCS) repose donc sur une estimation de l'épaisseur optique des aérosols à partir des images elles-mêmes. Pour bien comprendre le fonctionnement de cette méthode, il faut déjà comprendre les effets des aérosols sur le rayonnement. On a vu, dans ce billet, que les effets de la diffusion peuvent être modélisés ainsi (on suppose l'absorption gazeuse corrigée) :

ρTOA = ρatm +Td ρsurf

La réflectance au sommet de l'atmosphère ρTOA (Top of Atmosphere) est la somme de la réflectance atmosphérique ρatm et de la réflectance de surface ρsurf transmise par l'atmosphère. On cherche à connaître la réflectance de surface, mais à chaque mesure réalisée au sommet de l'atmosphère, on a trois inconnues à déterminer. Pour séparer les effets de l'atmosphère et les effets de la surface, il faut donc utiliser d'autres informations.

 

Méthode du pixel noir

Lorsque l'image contient une surface dont la réflectance est quasi nulle, la réflectance observée au sommet de l'atmosphère devient ρTOA= ρatm. On peut donc en déduire la réflectance atmosphérique, et en utilisant un modèle de transfert radiatif, l'épaisseur optique des d'aérosols. On peut enfin en déduire la transmission diffuse, et finalement calculer ρsurf. Une version encore plus simple et plus approximative consiste à soustraire directement la réflectance du pixel sombre (soit ρatm) à toute l'image. [Chavez, 1988]

 

Cependant, cette méthode revient à supposer qu'il existe bien une surface très sombre dans l'image (ce qui n'est pas toujours le cas), et que la réflectance de cette surface sombre est connue. La méthode suppose aussi que la quantité d'aérosols est constante dans l'image et elle néglige les effets du relief. Les résultats obtenus par cette méthode peuvent donc être assez imprécis. Dans notre méthode (MACCS), nous utilisons cependant la méthode du pixel noir déterminer la valeur maximale de l'épaisseur optique dans la zone.

 

Méthode Multi Spectrale, dite "DDV"

Si on connaît le type d'aérosols présent dans l'atmosphère, il est possible de déduire les  propriétés des aérosols dans une bande spectrale, à partir des propriétés optiques dans une autre bande spectrale.

 

Si on dispose de deux bandes spectrales, on dispose de deux mesures ρsurf et de trois inconnues( les deux réflectances de surface dans ces bandes, et la quantité d'aérosols). Une équation supplémentaire peut être obtenue si on connaît la relation entre les réflectances de surface des deux bandes.

 

La méthode  méthode "Dark Dense Vegetation" (DDV ) est basée sur des hypothèses de relations entre réflectances de surface sur la végétation dense exploitant le fait que le spectre de la végétation dense et verte est un peu toujours le même. La version la plus connue de cette méthode est celle utilisée par la NASA pour le projet MODIS [Remer 2005]. Elle relie les réflectances de surface dans le bleu et dans le rouge avec celles dans le moyen infra-rouge. On dispose ainsi de deux équations qui permettent d’estimer le type d’aérosols et l’épaisseur optique. Cette méthode fonctionne bien en zones tempérées et boréales, mais pas en zones arides, où il est difficile de trouver de la végétation dense. Les premières versions utilisaient les équations suivante

ρBleu = 0.5 ∗ ρSWIR

ρRouge = 0.25 ∗ ρSWIR

Les versions suivantes ont un peu compliqué ces équations, sans en modifier le principe. Nos travaux ont montré que l’utilisation de l'équation ci dessous  (la valeur exacte du coefficient est à ajuster en fonction des bandes spectrales de l'instrument):

ρBleu = 0.5 ∗ ρRouge

permet une détermination plus précise de l’épaisseur optique, pour des couverts végétaux moins denses (jusqu’à un NDVI de 0.2), car les sols nus de couleur marron respectent aussi cette relation. La méthode ne permet pas, par contre, de déterminer le modèle d’aérosols. Dans le cas de SPOT4 (Take5) l'absence d'une bande bleue ne nous permet pas d'utiliser cette dernière équation, d’où une légère perte en précision.

Ce diagramme montre que la corrélation entre réflectances de surface au dessus de la végétation est bien meilleure pour le couple de bandes spectrales (bleu, rouge) que pour les couples incluant le moyen infra rouge. (SWIR)

 

Méthode Multi Temporelle

On observe dans la plupart des cas que les réflectances de la surface terrestre évoluent lentement avec le temps, alors que le propriétés optiques des aérosols varient très rapidement, d'un jour à l'autre. On peut donc considérer que ce qui change d'une image à l'autre (en dehors de cas particuliers souvent liées à des interventions humaines) est lié aux aérosols, et donc en déduire les propriétés des aérosols pour ensuite corriger les effets atmosphériques. Cette méthode est un peu trop complexe pour être expliquée en détails ici, les lecteurs intéressés pourront se reporter à [Hagolle 2008].

 

Pour que les réflectances de surface soient quasi constantes d'une image à l'autre, il faut cependant que les images soient acquises sous un angle de vue constant. Les changements d'angles d'observation font en effet varier les réflectances (ce phénomène sera prochainement expliqué dans un autre article). Cette méthode ne s'applique donc qu'aux seuls satellites permettant des observations à angle constant.  Elle ne s'applique donc pas aux données SPOT normales mais par contre convient parfaitement aux données SPOT4 (Take5). Elle s'appliquera aussi à Landsat, Venµs et Sentinel-2.

En résumé :

Performance de l'estimation de l'épaisseur optique des aérosols sur des séries temporelles d'images Formosat-2,, en fonction de la méthode (multi-spectrale, multi-temporelle, combinée), par comparaison avec les mesures fournies par le réseau de mesures in-situ Aeronet. La méthode multi spectrale fonctionne mieux sur des sites couverts de végétation et moins bien sur des sites arides, la méthode multi-temporelle marche un peu moins bien sur les sites verts, mais beaucoup mieux sur les sites arides. La combinaison des deux méthodes garde le meilleur des deux méthodes élémentaires.

 

Notre méthode MACCS, utilisée pour l'expérience SPOT4 (Take5), et pour les données LANDSAT, VENµS et Sentinel-2, combine les trois méthodes présentées ci-dessus pour obtenir des estimations robustes des épaisseurs optiques d'aérosols. Ces méthodes fonctionnent dans un grand nombre de cas, mais peuvent parfois échouer quand les hypothèses sur lesquelles elles reposent s'avèrent fausses. Elles ont en général tendance à mieux fonctionner sur des zones couvertes de végétation plutôt que dans des zones arides. pour le moment, elles supposent le modèle d'aérosol connu, et dans les prochaines années, nous chercherons des manières fiables d'identifier le type d'aérosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

L'estimation du contenu atmosphérique en aérosols

=>

Attention, cet article contient des formules !

 

Les aérosols jouent un rôle prépondérant dans les effets atmosphériques. Les aérosols sont des particules en suspension dans l'atmosphère, qui peuvent être de plusieurs types : grains de sable ou poussières, suies issues de combustion, sulfates ou sels marins entourés d'eau... Leur taille peut varier de 0.1 µm à quelques microns, en fonction du type d'aérosols ou de l'humidité de l'air. Quant à leur quantité, elle est extrêmement variable, une pluie pouvant réduire brutalement leur abondance (on parle d'"épaisseur optique d'aérosols"). Ils peuvent faire varier fortement d'un jour à l'autre les réflectances observables depuis le sommet de l'atmosphère et il est donc nécessaire de connaître leur quantité et leur type afin de pouvoir corriger leurs effets.

 

Malheureusement, pour corriger les effets des aérosols, on ne dispose pas de réseau global d'observation des aérosols, seulement d'observations locales, sur les quelques centaines de points du réseau Aeronet. Ce réseau ne peut donc pas être utilisé pour corriger opérationnellement les images de satellites sur de grandes étendues.
Des modèles météorologiques commencent à prédire les quantités d'aérosols, en se basant sur les observations de satellites et la modélisation des sources et du transport des aérosols par les vents, mais ces données ne semblent pas encore avoir une précision suffisante pour être utilisées pour la correction atmosphérique des images.

 

Notre méthode de correction atmosphérique (MACCS) repose donc sur une estimation de l'épaisseur optique des aérosols à partir des images elles-mêmes. Pour bien comprendre le fonctionnement de cette méthode, il faut déjà comprendre les effets des aérosols sur le rayonnement. On a vu, dans ce billet, que les effets de la diffusion peuvent être modélisés ainsi (on suppose l'absorption gazeuse corrigée) :

ρTOA = ρatm +Td ρsurf

La réflectance au sommet de l'atmosphère ρTOA (Top of Atmosphere) est la somme de la réflectance atmosphérique ρatm et de la réflectance de surface ρsurf transmise par l'atmosphère. On cherche à connaître la réflectance de surface, mais à chaque mesure réalisée au sommet de l'atmosphère, on a trois inconnues à déterminer. Pour séparer les effets de l'atmosphère et les effets de la surface, il faut donc utiliser d'autres informations.

 

Méthode du pixel noir

Lorsque l'image contient une surface dont la réflectance est quasi nulle, la réflectance observée au sommet de l'atmosphère devient ρTOA= ρatm. On peut donc en déduire la réflectance atmosphérique, et en utilisant un modèle de transfert radiatif, l'épaisseur optique des d'aérosols. On peut enfin en déduire la transmission diffuse, et finalement calculer ρsurf. Une version encore plus simple et plus approximative consiste à soustraire directement la réflectance du pixel sombre (soit ρatm) à toute l'image. [Chavez, 1988]

 

Cependant, cette méthode revient à supposer qu'il existe bien une surface très sombre dans l'image (ce qui n'est pas toujours le cas), et que la réflectance de cette surface sombre est connue. La méthode suppose aussi que la quantité d'aérosols est constante dans l'image et elle néglige les effets du relief. Les résultats obtenus par cette méthode peuvent donc être assez imprécis. Dans notre méthode (MACCS), nous utilisons cependant la méthode du pixel noir déterminer la valeur maximale de l'épaisseur optique dans la zone.

 

Méthode Multi Spectrale, dite "DDV"

Si on connaît le type d'aérosols présent dans l'atmosphère, il est possible de déduire les  propriétés des aérosols dans une bande spectrale, à partir des propriétés optiques dans une autre bande spectrale.

 

Si on dispose de deux bandes spectrales, on dispose de deux mesures ρsurf et de trois inconnues( les deux réflectances de surface dans ces bandes, et la quantité d'aérosols). Une équation supplémentaire peut être obtenue si on connaît la relation entre les réflectances de surface des deux bandes.

 

La méthode  méthode "Dark Dense Vegetation" (DDV ) est basée sur des hypothèses de relations entre réflectances de surface sur la végétation dense exploitant le fait que le spectre de la végétation dense et verte est un peu toujours le même. La version la plus connue de cette méthode est celle utilisée par la NASA pour le projet MODIS [Remer 2005]. Elle relie les réflectances de surface dans le bleu et dans le rouge avec celles dans le moyen infra-rouge. On dispose ainsi de deux équations qui permettent d’estimer le type d’aérosols et l’épaisseur optique. Cette méthode fonctionne bien en zones tempérées et boréales, mais pas en zones arides, où il est difficile de trouver de la végétation dense. Les premières versions utilisaient les équations suivante :

 

ρBleu = 0.5 ∗ ρSWIR

ρRouge = 0.25 ∗ ρSWIR

 

Les versions suivantes ont un peu compliqué ces équations, sans en modifier le principe. Nos travaux ont montré que l’utilisation de l'équation ci dessous  (la valeur exacte du coefficient est à ajuster en fonction des bandes spectrales de l'instrument):

ρBleu = 0.5 ∗ ρRouge

 

permet une détermination plus précise de l’épaisseur optique, pour des couverts végétaux moins denses (jusqu’à un NDVI de 0.2), car les sols nus de couleur marron respectent aussi cette relation. La méthode ne permet pas, par contre, de déterminer le modèle d’aérosols. Dans le cas de SPOT4 (Take5) l'absence d'une bande bleue ne nous permet pas d'utiliser cette dernière équation, d’où une légère perte en précision.

Ce diagramme montre que la corrélation entre réflectances de surface au dessus de la végétation est bien meilleure pour le couple de bandes spectrales (bleu, rouge) que pour les couples incluant le moyen infra rouge. (SWIR)

 

Méthode Multi Temporelle

On observe dans la plupart des cas que les réflectances de la surface terrestre évoluent lentement avec le temps, alors que le propriétés optiques des aérosols varient très rapidement, d'un jour à l'autre. On peut donc considérer que ce qui change d'une image à l'autre (en dehors de cas particuliers souvent liées à des interventions humaines) est lié aux aérosols, et donc en déduire les propriétés des aérosols pour ensuite corriger les effets atmosphériques. Cette méthode est un peu trop complexe pour être expliquée en détails ici, les lecteurs intéressés pourront se reporter à [Hagolle 2008].

 

Pour que les réflectances de surface soient quasi constantes d'une image à l'autre, il faut cependant que les images soient acquises sous un angle de vue constant. Les changements d'angles d'observation font en effet varier les réflectances : c'est ce qu'on appelle les effets directionnels. Cette méthode ne s'applique donc qu'aux seuls satellites permettant des observations à angle constant.  Elle ne s'applique donc pas aux données SPOT normales mais par contre convient parfaitement aux données SPOT4 (Take5). Elle s'appliquera aussi à Landsat, Venµs et Sentinel-2.

 

En résumé :
Performance de l'estimation de l'épaisseur optique des aérosols sur des séries temporelles d'images Formosat-2,, en fonction de la méthode (multi-spectrale, multi-temporelle, combinée), par comparaison avec les mesures fournies par le réseau de mesures in-situ Aeronet. La méthode multi spectrale fonctionne mieux sur des sites couverts de végétation et moins bien sur des sites arides, la méthode multi-temporelle marche un peu moins bien sur les sites verts, mais beaucoup mieux sur les sites arides. La combinaison des deux méthodes garde le meilleur des deux méthodes élémentaires.

 

Notre méthode MACCS, utilisée pour l'expérience SPOT4 (Take5), et pour les données LANDSAT, VENµS et Sentinel-2, combine les trois méthodes présentées ci-dessus pour obtenir des estimations robustes des épaisseurs optiques d'aérosols. Ces méthodes fonctionnent dans un grand nombre de cas, mais peuvent parfois échouer quand les hypothèses sur lesquelles elles reposent s'avèrent fausses. Elles ont en général tendance à mieux fonctionner sur des zones couvertes de végétation plutôt que dans des zones arides. pour le moment, elles supposent le modèle d'aérosol connu, et dans les prochaines années, nous chercherons des manières fiables d'identifier le type d'aérosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENμS and Sentinel-2 Images. Remote Sens. 2015, 7, 2668-2691.

Land cover map production: how it works

=>

Land cover and land use maps

Although different, the terms land use and land cover are often used as synonymous. From Wikipedia Land cover is the physical material at the surface of the earth. Land covers include grass, asphalt, trees, bare ground, water, etc. There are two primary methods for capturing information on land cover: field survey and analysis of remotely sensed imagery. and Land use is the human use of land. Land use involves the management and modification of natural environment or wilderness into built environment such as fields, pastures, and settlements. It also has been defined as "the arrangements, activities and inputs people undertake in a certain land cover type to produce, change or maintain it" (FAO, 1997a; FAO/UNEP, 1999).


A precise knowledge of land use and land cover is crucial for many scientific studies and for many operational applications. This accurate knowledge needs frequent information updates, but may also need to be able to go back in time in order to perform trend analysis and to suggest evolution scenarios.

 

Satellite remote sensing offers the possibility to have a global point of view over large regions with frequent updates, and therefore it is a very valuable tool for land cover map production.

 

However, for those maps to be available in a timely manner and with a good quality, robust, reliable and automatic methods are needed for the exploitation of the available data.

 

 

 

Classical production approaches

The automatic approaches to land cover map production using remote sensing imagery are often based on image classification methods.

 

This classification can be:

  • supervised: areas for which the land cover is known are used as learning examples;
  • unsupervised: the image pixels are grouped by similarity and the classes are identified afterwards.

Supervised classification often yields better results, but it needs reference data which are difficult or costly to obtain (field campaigns, photo-interpretation, etc.).

 

 

 

What time series bring

Until recently, fine scale land cover maps have been nearly exclusively produced using a small number of acquisition dates due to the fact that dense image time series were not available.

 

The focus was therefore on the use of spectral richness in order to distinguish the different land cover classes. However, this approach is not able to differentiate classes which may have a similar spectral signature at the acquisition time, but that would have a different spectral behaviour at another point in time (bare soils which will become different crops, for instance). In order to overcome this problem, several acquisition dates can be used, but this needs a specific date selection depending on the map nomenclature.

 

For instance, in the left image, which is acquired in May, it is very difficult to tell where the rapeseed fields are since they are very similar to the wheat ones. On the right image, acquired in April, blooming rapeseed fields are very easy to spot.

 

May image. Light green fields are winter crops, mainly wheat and rapeseed. But which are the rapeseed ones?

April image. Blooming rapeseed fields are easily distinguished in yellow while wheat is in dark green.

 

If one wants to build generic (independent from the geographic sites and therefore also from the target nomenclatures) and operational systems, regular and frequent image acquisitions have to be ensured. This will soon be made possible by the Sentinel-2 mission, and it is right now already the case with demonstration data provided by Formosat-2 and SPOT4 (Take 5). Furthermore, it can be shown that having a high temporal resolution is more interesting than having a high spectral diversity. For instance, the following figure shows the classification performance results (in terms of  \kappa index, the higher the better) as a function of the number of images used. Formosat-2 images (4 spectral bands) and simulated Sentinel-2 (13 bands) and Venµs (12 bands) data have been used. It can be seen that, once enough acquisitions are available, the spectral richness is caught up by a fine description of the temporal evolution.

kappaVFS.png

 

 

What we can expect from Sentinel-2

Sentinel-2 has unique capabilities in the Earth observation systems landscape:

  • 290 km. swath;
  • 10 to 60 m. spatial resolution depending on the bands;
  • 5-day revisit cycle with 2 satellites;
  • 13 spectral bands.

Systems with similar spatial resolution (SPOT or Landsat) have longer revisit periods and fewer and larger spectral bands. Systems with similar temporal revisit have either a lower spatial resolution (MODIS) or narrower swaths (Formosat-2).

 

The kind of data provided by Sentinel-2 allows to foresee the development of land cover map production systems which should be able to update the information monthly at a global scale. The temporal dimension will allow to distinguish classes whose spectral signatures are very similar during long periods of the year. The increased spatial resolution will make possible to work with smaller minimum mapping units.

 

However, the operational implementation of such systems will require a particular attention to the validation procedures of the produced maps and also to the huge data volumes. Indeed, the land cover maps will have to be validated at the regional or even at the global scale. Also, since the reference data (i.e. ground truth) will be only available in limited amounts, supervised methods will have to be avoided as much as possible. One possibility consists of integrating prior knowledge (about the physics of the observed processes, or via expert rules) into the processing chains.

 

Last but not least, even if the acquisition capabilities of these new systems will be increased, there will always be temporal and spatial data holes (clouds, for instance). Processing chains will have to be robust to this kind of artefacts.

 

 

Ongoing work at CESBIO

 

Danielle Ducrot, Antoine Masse and a few CESBIO interns have recently produced a a large land cover map over the Pyrenees using 30 m. resolution multi-temporal Landsat images. This map, which is real craftsmanship, contains 70 different classes. It is made of 3 different parts using nearly cloud-free images acquired in 2010.

 

70-class land cover map obtained from multi-temporal Landsat data.

In his PhD work, Antoine works on methods allowing to select the best dates in order to perform a classification. At the same time, Isabel Rodes is looking into techniques enabling the use of all available acquisitions over very large areas by dealing with both missing data (clouds, shadows) and the fact that all pixels are not acquired at the same dates.

 

These 2 approaches are complementary: one allows to target very detailed nomenclatures, but needs some human intervention, and the other is fully automatic, but less ambitious in terms of nomenclature.

 

A third approach is being investigated at CESBIO in the PhD work of Julien Osman: the use of prior knowledge both quantitative (from historical records) and qualitative (expert knowledge) in order to guide the automatic classification systems.

 

We will give you more detailed information about all those approaches in coming posts on this blog.