# OUTILS MÉTHODOLOGIQUES POUR L'ANALYSE D'IMAGES MSG : ESTIMATION DU MOUVEMENT, SUIVI DE MASSES NUAGEUSES ET DÉTECTION DE FRONTS

Thomas Corpetti¹,², Vincent Dubreuil¹, Etienne Mémin³, Nicolas Papadakis⁴, Olivier Planchon¹, Claire Thomas¹,³

1. CNRS, Université de Rennes 2, COSTEL UMR 6554 LETG / OSUR, Rennes, France

2. CNRS LIAMA / TIPE, Pékin, Chine

3. INRIA, Centre Rennes Bretagne-Atlantique, Rennes, France

4. CNRS, IMB, Bordeaux, France

## Résumé

Cet article présente les travaux que nous avons menés ces dernières années autour de l'analyse d'images Météosat Seconde Génération (MSG). Comparés à la première génération, les données MSG possèdent une résolution spatiale et temporelle plus élevée, autorisant l'accès à un certain nombre d'informations liées aux phénomènes climatiques observés. Cependant, l'étape consistant à remonter à cette information physique à partir des données images s'avère délicate car nous sommes confrontés à des structures soumises à de très fortes déformations, parfois observées en transparence et avec une durée de vie variable. Par conséquent, les outils classiques issus de l'analyse d'images se révèlent souvent limités et il est nécessaire de les adapter à cette spécificité. Nous nous focalisons ici sur trois applications particulières : l'estimation du mouvement, permettant de remonter aux vents atmosphériques, le suivi de masses nuageuses, autorisant par exemple l'analyse de phénomènes convectifs et la détection de fronts, appliquée ici aux brises de mer.

**Mots clés :** Images MSG, estimation du mouvement, cellules convectives, brises de mer.

## Abstract

This article presents some researches we performed the last decade on the analysis of Meteosat Second Generation (MSG) images. Compared to the first generation, spatial and temporal resolution of MSG data is higher and this authorizes now to derive a number of information about observed phenomena. However, extracting a physical quantity from MSG luminance is a tricky task since we have to deal with complex and partly observed structures, submitted to sudden deformations with a limited life cycle. Therefore usual computer vision tools are limited and we need to design adapted approaches to MSG specificities. This is the scope of this paper. We focus in particular on three key applications: motion estimation (in order to recover atmospheric flows), cloud tracking (applied to convective cells) and front detection, applied to sea breeze fronts.

**Keywords :** Meteosat Second Generation, motion estimation, curve tracking, front detection.

## 1. Introduction

La télédétection en général et les satellites géostationnaires en particulier sont des outils cruiaux pour accéder à un certain nombre de paramètres liés à la physique atmosphérique. Par rapport aux outils de mesure locales (sondes, ballons, ...), les images satellites offrent une large couverture spatiale et une très bonne résolution temporelle ce qui est primordial dans certains endroits du globe où ni les sondes locales, ni les modèles physiques usuels et leur simulation numérique ne peuvent fournir une information fiable. A titre d'exemple, dans de nombreux secteurs clés du point de vue climatique (comme la forêt amazonienne ou le désert saharien), il est très difficile d'estimer exactement des paramètres physiques à l'aide des sondes locales en raison de la densité limitée des réseaux de mesure et de la variabilité spatiale et temporelle trop élevée des

paramètres climatiques (précipitations tropicales, température).

Dans ce contexte, les satellites géostationnaires Météosat et en particulier MSG (Météosat Seconde Génération) ont été conçus en réponse aux exigences des utilisateurs pour accéder à des mesures physiques précises et comportent à présent 12 bandes spectrales (au lieu de 3), une meilleure résolution spatiale (3km au lieu de 5km) et temporelle (15 min au lieu de 30 min). Ils préfigurent la nouvelle génération de satellites géostationnaires de GOES (GOES-Next). Avec ces nouvelles données, des études plus fines sur la forme et la dynamique de certaines structures nuageuses sont envisageables et ouvrent un vaste champ applicatif.

Le suivi des formations nuageuses est en effet un enjeu majeur en termes de climatologie et météorologie. Sur le plan global, le couvert nuageux est une des com-
posantes majeures du bilan radiatif terrestre d'une part en raison de son albédo et son effet écran (parasol) dans les courtes longueurs d'ondes et d'autre part en raison de son rôle majeur dans les grandes longueurs d'ondes : absorption du rayonnement infrarouge terrestre et limitation des pertes radiatives du système terre-atmosphère. Le rôle central des nuages dans le mécanisme d'effet de serre justifie les programmes de recherches spécifiques dédiés à leur étude comme ISCCP (International Satellite Cloud Climatology Program) ou AMIP (Atmospheric Model Intercomparison Project). Sur le plan météorologique, le suivi des formations nuageuses est aussi décisif pour la prévision du temps : outre leur importance pour la définition du type de temps (état du ciel, amplitude thermique diurne, etc...), les nuages sont aussi associés à des phénomènes atmosphériques divers comme la pluie, les orages, etc. Sur le plan des mécanismes il convient cependant de distinguer les échelles auxquelles ils se forment :

* à l'échelle synoptique, les systèmes s'organisent en fonction des masses d'air et des centres d'action : Zone de Convergence Inter Tropicale (ZCIT), fronts, cyclones, etc

* à l'échelle régionale (meso-scale), les formations résultent de dynamiques spécifiques où s'imbriquent échelles synoptiques et facteurs locaux comme l'illustrent les fronts de brise de mer ou les cellules convectives tropicales;

* à l'échelle plus locale, certaines formations résultent fréquemment d'interactions avec la surface terrestre : brouillards ou nuages orographiques par exemple.

Analyser de telles structures nuageuses dans les images MSG est donc d'un intérêt majeur en vue d'une meilleur compréhension sur les phénomènes météorologiques ou climatiques. Cependant, à l'inverse de capteurs locaux qui mesurent en général directement une quantité physique désirée (vent, température, ...), l'exploitation d'images satellites requiert le recouvrement des quantités d'intérêt à partir des seules informations de luminance des différentes bandes spectrales. Cette étape d'inversion nécessite le recours à des techniques de traitement d'images où un certain nombre d'outils très évolués sont à l'heure actuelle disponible et utilisés dans les centres opérationnels. Néanmoins, en raison de la complexité des lois qui gouvernement la dynamique atmosphérique, des problématiques de projections 3D/2D, des relations non triviales entre la quantité d'intérêt et la luminance des images, la précision des mesures actuelles par imagerie reste limitée.

Dans cet article, nous présentons des outils d'analyse d'images spécifiquement dédiés aux images MSG. Nous nous focalisons en particulier sur trois tâches primordiales : l'estimation du mouvement, le suivi de courbes et la détection de fronts, présentés respectivement dans les sections 2, 3 et 4. Dans chacun des cas,

nous présentons les fondements scientifiques des outils usuels existants et introduisons nos contributions.

## 2. Estimation du mouvement

### 2.1. Problématique

L'estimation du mouvement est la tache consistant à mesurer le la vitesse apparente à partir d'images satellites. C'est une opération largement effectuée en météorologie opérationnelle. A partir de deux images $I_{t-1}(x) = I(x, t-1)$ au temps $t-1$ et $I_t(x) = I(x, t)$ au temps $t$, définies pour tout point $x \in \Omega$ du domaine image $\Omega$, l'objectif consiste à extraire le champ bidimensionnel de vitesse apparente $v(x) = (u(x), v(x))^T$. Pour cela, il est courant de s'appuyer sur l'Equation de Contrainte du Mouvement Apparent (ECMA) :

$$ \frac{dI}{dt} = \frac{\partial I(x, t)}{\partial t} + v(x, t) \cdot \nabla I(x, t) \sim 0, \quad (1) $$

qui assume que les points $x$ conservent leur intensité le long du déplacement, la luminance $I$ étant vue comme une fonction continue et $\nabla = (\partial/\partial x, \partial/\partial y)^T$ étant l'opérateur de gradient. Appliqué à une paire d'images, cette relation s'écrit :

$$ \begin{gathered} I_t(x + \Delta t v(x)) - I_{t-1}(x) = 0 \Rightarrow \\ \frac{I_t(x) - I_{t-1}(x)}{\Delta t} + v(x) \cdot \nabla I_t(x) = 0, \end{gathered} \quad (2) $$

où nous avons un développement de Taylor au premier ordre de la contrainte de conservation $I_t(x + \Delta t v(x)) - I_{t-1}(x) = 0$ autour de $\Delta t v(x)$, $\Delta t$ étant l'intervalle de temps entre deux images (ici, et par convention, on suppose que $\Delta t = 1$). Cette relation créé un lien entre l'image des différences $I_t(x) - I_{t-1}(x)$, les gradients spatiaux de la seconde image $\nabla I_t(x)$ et la vitesse $v$ à extraire. Les équations (1-2) sont à la base de nombreuses études en estimation du mouvement. Cependant, à ce point, on peut remarquer que :

1. Dans les zones homogènes, tous les termes s'annulent et il y a une infinité de solutions;

2. En raison de la projection $v \cdot \nabla I$, seule la composante normale aux gradients photométriques peut être extraite (ce problème est connu sous le nom de "problème de l'ouverture").

Par conséquent, les relations (1-2) ne sont pas suffisantes pour extraire un champ de vitesse cohérent et il est nécessaire d'y associer d'autres contraintes. A l'heure actuelle, les techniques utilisées en météorologie opérationnelle s'appuient sur de la corrélation, comme le produit AMV (Atmospheric Motion Vectors, cf. Schmetz et al., 1993, 1995) d'Eumetsat, l'agence européenne en charge des données Météosat. Dans ce contexte, la vitesse $v(x)$ est localement choisie comme étant celle qui maximise un critère de corrélation $C_W$ entre une fenêtre $W(x)$ centrée en $x$ dans la première image et une fenêtre déplacée en $x + v(x)$ dans la seconde :

$$ v = \max_{v = \{-U, \dots, U\} \times \{-V, \dots, V\}} \sum_{x \in W(x)} C_W(I(x + v, t), I(x, t-1)). \quad (3) $$
L'espace des solutions possibles pour $v$ est discret et fixé a priori. En général le critère de corrélation $C_W$ s'appuie sur la conservation de la luminance (relation (2)), ce qui donne:

$$C_W(I(\boldsymbol{x} + \boldsymbol{v}, t), I(\boldsymbol{x}, t - 1)) = I(\boldsymbol{x} + \boldsymbol{v}, t)I(\boldsymbol{x}, t - 1), \quad (4)$$

ou, en version normalisée:

$$C_W(I(\boldsymbol{x} + \boldsymbol{v}, t), I(\boldsymbol{x}, t - 1)) = \frac{(I(\boldsymbol{x} + \boldsymbol{v}, t) - \bar{I}(\boldsymbol{x} + \boldsymbol{v}, t))(I(\boldsymbol{x}, t - 1) - \bar{I}(\boldsymbol{x}, t - 1))}{\sigma_{W,I}(\boldsymbol{x}+\boldsymbol{v},t)\sigma_{W,I}(\boldsymbol{x},t-1)}, \quad (5)$$

où $\bar{I}(\bullet)$ et $\sigma_{W,I}(\bullet,t)$ sont la moyenne et l'écart type de $\bar{I}(\bullet)$ dans la fenêtre $W$. Ces techniques sont à l'heure actuelle très performantes pour recouvrer un champ de vent consistant et sont la base de nombre de produits météorologiques. Néanmoins, deux principaux défauts peuvent être mentionnés :

1. Le choix de la taille de la fenêtre de corrélation est primordial : avec une petite taille on risque de manipuler des blocs homogènes et donc d'être soumis au problème de l'ouverture tandis qu'une grande taille nous privera de l'estimation de mouvements à petite échelle;

2. En raison de la nature discrète de l'ensemble des solutions possibles, une perte de précision en découle. Cela perturbe parfois l'analyse des champs de vitesses, en particulier dans le calcul de ses gradients.

Une manière classique de gérer ces problèmes est d'exploiter des méthodes efficaces de post-traitement, notamment en assimilant les champs de vitesses dans des modèles atmosphériques (cf. par exemple (Menzel et al., 1983; Holmlund, 1998; Schmetz et al., 1993, 1995)). Cependant, la qualité de ce post-traitement par assimilation de "pseudo-observations" (dans la mesure où les données assimilées se sont pas directement observées) est étroitement liée à la qualité des vitesses d'entrées. Il semblerait alors plus pertinent de définir des techniques efficaces d'assimilation de données qui prennent en compte directement les images pour en extraire un champ dense de déplacement, au lieu de calculer ces champs dans un premier temps pour les assimiler dans une seconde étape. C'est ce que nous proposerons dans la section 2.3. Avant cela, la section suivante présente les principaux concepts de l'assimilation de données.

## 2.2. Assimilation de données

Nous introduisons ici l'assimilation de données pour des modèles dynamiques imparfaits. Pour une présentation méthodologique complète ainsi que des applications dans le contexte de fluides géophysique, nous renvoyons le lecteur vers les références suivantes (Bennett, 1992; Le Dimet et Talagrand, 1986; Lions, 1971; Papadakis, 2007; Talagrand et Courtier, 1987; Talagrand, 1997; Vidard et al., 2000). Le problème auquel nous nous attachons ici consiste à recouvrer, à partir d'un état initial,

un espace d'état $X$ observé partiellement et soumis à une loi dynamique imparfaite. Cela se formalise comme la recherche des états $X(\boldsymbol{x}, t)$, en toute position $\boldsymbol{x}$ au temps $t \in [t_0, t_f]$, qui satisfont les conditions:

$$\frac{\partial X}{\partial t}(\boldsymbol{x}, t) + \mathbf{M}(X(\boldsymbol{x}, t)) = \nu_m(\boldsymbol{x}), \quad (6)$$

$$X(\boldsymbol{x}, t_0) = X_0(\boldsymbol{x}) + \nu_n(\boldsymbol{x}), \quad (7)$$

$$\mathcal{Y}(\boldsymbol{x}, t) = \mathbb{H}(X(\boldsymbol{x}, t)) + \nu_o(\boldsymbol{x}, t), \quad (8)$$

où $\mathbb{M}$ est le modèle dynamique non-linéaire, $X_0$ est la condition initiale au temps $t_0$ et $(\nu_n, \nu_m)$ sont des variables de contrôle additives liées au bruit (supposé connu) sur la dynamique et la condition initiale. Par ailleurs, des observations bruitées $\mathcal{Y}$ de l'espace d'état sont disponibles par le biais de l'opérateur d'observation non linéaire $\mathbb{H}$, à un bruit près $\nu_o$. Afin d'estimer $X$, une méthodologie classique consiste à minimiser une fonctionnelle $\mathcal{J}$ du type :

$$\begin{aligned} \mathcal{J}(X) = & \frac{1}{2} \int_{t_0}^{t_f} \|\mathcal{Y} - \mathbb{H}(X(\nu_m, \nu_n))\|_{R^{-1}}^2 dt \\ & + \frac{1}{2} \|X(\boldsymbol{x}, t_0) - X_0(\boldsymbol{x})\|_{B^{-1}}^2 \\ & + \frac{1}{2} \int_{t_0}^{t_f} \|\frac{\partial X}{\partial t}(\boldsymbol{x}, t) + \mathbb{M}(X(\boldsymbol{x}, t))\|_{Q^{-1}}^2 dt, \end{aligned} \quad (9)$$

où nous avons introduit les matrices $R, B, Q$ liées à la covariance des erreurs $(\nu_m, \nu_n, \nu_o)$. La distance de Mahalanobis utilisée ici se lit, pour une matrice d'information $A: \|X\|_{A^{-1}} = X^T A^{-1} X$. L'évaluation de $X$ se fait alors en annulant le gradient $\delta\mathcal{J}_X(\theta) = \lim_{\beta \to 0} \frac{J(X+\beta\theta)-J(X)}{\beta}$ de (9). Malheureusement, l'évaluation de ce gradient est en pratique irréaliste pour un espace d'état de grande dimension (comme un champ dense de déplacement par exemple) car cela requiert d'intégrer le modèle dynamique à chaque fois que nous appliquons une perturbation autour de $X$, ce qui devient impossible du point de vue calculatoire pour les grands espaces d'états tels que ceux rencontrés dans notre cas.

Une manière élégante de résoudre ce problème est de s'appuyer sur une formulation adjointe du problème. Pour cela, des variables adjointes $\lambda$, exprimant les erreurs entre l'état courant et le modèle dynamique, sont introduites:

$$\lambda = Q^{-1} \left( \frac{\partial X}{\partial t} + \mathbb{M}(X) \right). \quad (10)$$

En notant

* $(\frac{\partial \mathbb{M}}{\partial \tilde{X}})$ et $(\frac{\partial \mathbb{H}}{\partial \tilde{X}})$ les opérateurs tangents linéaires de $\mathbb{M}$ et $\mathbb{H}^1$;

$$\left(\frac{\partial \mathbb{A}}{\partial \tilde{X}}\right) (dX) = \lim_{\beta \to 0} \frac{\mathbb{A}(\tilde{X} + \beta dX) - \mathbb{A}(\tilde{X})}{\beta}, \quad (11)$$

<footnote n="1">¹L'opérateur tangent linéaire $(\frac{\partial \mathbb{A}}{\partial \tilde{X}})$ d'un opérateur $\mathbb{A}$ est sa dérivée directionnelle (la dérivée Gâteaux) :</footnote>
• $(\partial_X M)^*$ et $(\partial_X H)^*$ leurs opérateurs adjoints<fnref n="2" />,

il peut être démontré que l'annulation du gradient $\delta J_X(\theta)$ peut se faire efficacement à l'aide d'une intégration rétrograde du système dans un modèle d'évolution adjoint permettant d'estimer les variables $\lambda$. Une fois estimées, ces variables permettent de recouvrer l'espace d'état $X$ par la relation (10). Finalement, en présence de modèles non linéaires, l'estimation de $X$ conduit au système incrémental suivant (Bennett, 1992):

1. Partant de $\tilde{X}(x, t_0) = X_0(x)$, réaliser une intégration directe: $\frac{\partial \tilde{X}}{\partial t} + M(\tilde{X}) = 0$

2. Les variables $\tilde{X}(x, t)$ étant à présent disponibles, calculer les variables adjointes $\lambda(x, t)$ avec l'équation rétrograde :

$$ 
\begin{aligned}
\lambda(t_f) &= 0; \\
-\frac{\partial \lambda}{\partial t}(t) + (\partial_X M)^* \lambda(t) &= (13) \\
(\partial_X H)^* R^{-1} (\mathbf{Y} - H(\tilde{X}))(t)
\end{aligned}
$$

3. Mettre à jour la condition initiale : $dX(t_0) = B\lambda(t_0)$;

4. Avec les variables $\lambda$, calculer les variables $dX(t)$ à partir de $dX(t_0)$ avec l'intégration directe :

$$ \frac{\partial dX}{\partial t}(t) + \left( \frac{\partial M}{\partial \tilde{X}} \right) dX(t) = Q\lambda(t) \quad (14) $$

5. Mettre à jour : $\tilde{X} = \tilde{X} + dX$

6. Retour à l'étape (2) jusqu'à convergence.

Intuitivement, les variables adjointes $\lambda$ nous informent sur la différence entre les observations et le modèle dynamique. Elles sont calculées à partir d'une solution courante $\tilde{X}$ par l'équation rétrograde (13) qui prend en compte à la fois le modèle dynamique et les observations. Cette information est ensuite utile pour raffiner la condition initiale (étape (3)) et recouvrer l'espace d'état par un modèle dynamique imparfait comprenant comme erreurs $Q\lambda$ (étape (4)). Il convient de noter que si la dynamique est parfaite, la matrice de covariance associée $Q$ est nulle et l'algorithme revient à un raffinement de la condition initiale. Cela est souvent exploité en météorologie par exemple. Cependant ici, lorsque l'on manipule des images, une dynamique parfaite est difficile à obtenir en raison des problèmes de projection 3D/2D, des conditions d'illumination variables, des conditions aux bordes des images inconnues, ...

Le cadre de l'assimilation de données présenté ci-dessus est donc une solution intéressante pour l'estimation d'un espace d'état de grande dimension comme un champ dense de mouvement ou encore des courbes fermées (comme nous le verrons dans la section

3). Il est nécessaire de définir quatre quantités pour un système complet d'assimilation de données :

1. L'espace d'état $X$;

2. La modèle dynamique $M$, son tangent linéaire $(\frac{\partial M}{\partial \tilde{X}})$ et son adjoint $(\partial_X M)^*$;

3. L'opérateur d'observation $H$ et son adjoint $(\partial_X H)^*$;

4. Les matrices d'erreurs de covariance $B, Q, R$.

En général, pour des applications en sciences atmosphériques comme la météorologie, la quantité d'intérêt est directement observée et les modèles dynamiques sont parfaits. Cependant, quand on manipule des images, en raisons des difficultés mentionnées précédemment, nous devons définir des modèles dynamiques et des opérateurs d'observations adaptés à la spécificité des données, tout comme les matrices de covariances associées. Dans la section suivante, nous présentons une application dans le cadre de l'estimation de mouvements atmosphériques.

## 2.3. Estimation du mouvement par assimilation de données

A partir d'une séquence d'images $I_t(x)$, $t = \{1, ..., N\}$, l'objectif consiste à estimer les champs de déplacements associés en exploitant le formalisme de l'assimilation de données présenté dans la section précédente. Pour cela, nous définissons donc dans les quatre paragraphes qui suivent l'espace d'état, le modèle dynamique, l'opérateur d'observation ainsi que les matrices de covariance d'erreurs.

### 2.3.1. Espace d'état

Une manière naturelle de représenter un champ de vitesses est d'utiliser ses composantes horizontales $u$ et verticales $v$ dans l'espace Euclidien : $X = [u, v]^T$. Cependant, lorsque l'on couple cette représentation avec des modèles dynamiques, nous avons trouvé plus judicieux de s'appuyer sur ses composantes de divergence et de rotationnel (ou "vorticité"). Pour un champ de vitesses $2D v = (u, v)^T$, la divergence et la rotation sont définies par :

$$ 
\left\{ 
\begin{array}{l}
\mathrm{div} v = D = \nabla \cdot v = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \\
\mathrm{curl} v = \zeta = \nabla \wedge v = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}.
\end{array}
\right. 
\qquad (15) $$

Suivant la décomposition de Helmholtz, tout champ de déplacement $v$ peut se décomposer en une partie irrotationnelle $v_{irr}$, solénoïdale $v_{sol}$ et harmonique $v_{har}$. La partie irrotationnelle possède une vorticité nulle tandis que la partie solénoïdale n'a pas de divergence. La composante harmonique (ou laminaire) ne possède ni divergence, ni rotation. Il existe des techniques efficaces basées sur du flot optique variationnel (Horn et Schunck, 1981) pour estimer cette composante harmonique $v_{har}$ et nous la supposerons donc nulle dans nos applications. La figure 1 illustre la décomposition de Helmholtz.

Ainsi, les composantes de rotationnel et de divergence $(\zeta, D)$ permettent de décrire à elles seules le

$$ \forall x_1, x_2 \in D, \langle A x_1, x_2 \rangle = \langle x_1, A^* x_2 \rangle, \quad (12) $$

<footnote n="2">²l'opérateur adjoint $A^*$ d'un opérateur linéaire $A$ sur un espace $D$ se définit par la relation:</footnote>
champ de déplacement et nous choisissons donc cette représentation pour l'espace d'état $X = [\zeta, D]^T$. La connaissance de $(\zeta, D)$ permet de recouvrer les composantes $(u, v)$ du champ de vitesses par la loi de Biot-Savart:

$$v = \nabla^{\perp} G * \zeta + \nabla G * D = \mathbb{H}_G(X), \quad (16)$$

où nous avons introduit l'opérateur $\mathbb{H}_G : [\zeta, D]^T \to v = [u, v]^T$ tel que $\mathbb{H}_G(X) = \nabla^{\perp} G * \zeta + \nabla G * D$, $\nabla^{\perp} = (\partial/\partial y, -\partial/\partial x)^T$ et $G$ représente le noyau de Green associé à l'opérateur de Laplacien ($G = \frac{1}{2\pi} \ln(|x|)$). Le calcul de cet opérateur peut se réaliser de manière efficace dans l'espace de Fourier (Corpetti et al., 2003).

### 2.3.2. Modèle dynamique

Le but de cette partie est de définir un modèle dynamique suffisamment générique pour estimer des champs de vitesses à partir d'images MSG. Nous nous intéressons en premier lieu au modèle sur la vorticité et ensuite sur la divergence.

#### Dynamique de la vorticité

Sous l'hypothèse, réaliste pour les données MSG, d'incompressibilité moyenne (la divergence est faible et le mouvement est surtout influencé par les échelles horonzitales), nous nous appuyons sur la formulation vorticité-vitesse des équations de Navier-Stokes :

$$\frac{\partial \zeta}{\partial t} + v \cdot \nabla \zeta - \nu \Delta \zeta = 0, \quad (17)$$

qui est une équation classique d'advection-diffusion où $\nu$ est un coefficient de dissipation à fixer. Sa valeur peut être obtenue en se basant sur les théories de simulation aux grandes échelles et on peut en particulier s'appuyer soit :

1. Sur le modèle de Smagorinsky : $\nu = (C\delta_x)^2 \sqrt{2(u_x^2 + v_y^2 + (u_y + v_x)^2)}$ où $C$ est le coefficient de Smagorinsky (souvent fixé à 0.17) et $\delta_x$ la taille de la grille;

2. Sur des modèles sous-mailles basés sur l'entropie (Mansour et al., 1978). Cela se base sur les transferts de vorticité entre échelles (Taylor, 1932) et se lit $\nu = (C\delta_x)^2|\zeta|$.

En pratique dans nos applications, les deux approches se sont révélées cohérentes et nous avons privilégié la seconde, qui apparaît plus robuste. Le modèle en (17) est donc celui retenu pour la vorticité.

#### Dynamique de la divergence

Etant donné qu'aux larges échelles, la divergence peut être considérée comme faible, nous supposons simplement que cette quantité est advectée par le flot à un bruit près traduisant l'incertitude faite sur cette hypothèse. Plus précisément, nous modélisons la carte de divergence par le processus stochastique suivant :

$$dX_t = v(X_t)dt + \sqrt{2\nu} d\mathbf{B}_t. \quad (18)$$

Dans cette relation, $X_t$ représente un ensemble de "particules de divergence" comprenant chacune une valeur de cette quantité et $B_t$ est un mouvement Brownien standard de $\mathbb{R}^2$. Cette relation indique que chaque particule est transportée par le champ $v$ (advection) et sa position est connue à une incertitude près qui croit linéairement avec le temps. Il peut être démontré, en utilisant la formule d'Ito et l'équation avant de Kolmogorov (nous renvoyons le lecteur à (Oksendal, 1998) pour une introduction au calcul stochastique), que l'espérance $\xi(t, X) = E[\mathrm{div}v(X_t)]$ au temps $t$ d'une telle divergence suit également une équation d'advection-diffusion:

$$\begin{aligned} \xi_t + v \cdot \nabla \xi + \xi \mathrm{div}(v) - \nu \Delta \xi &= 0, \\ \xi(0, X_0) &= \mathrm{div}v(X_0). \end{aligned} \quad (19)$$

En supposant qu'aux larges échelles, la divergence est similaire à son espérance ($D \approx \xi$), nous obtenons le modèle de divergence simplifié :

$$\frac{\partial D}{\partial t} + v \cdot \nabla D + D^2 = \nu \Delta D. \quad (20)$$

Suivant cette relation, l'espérance de la divergence est advectée par le flot et dissipée en raison des incertitudes sous-mailles. Cette hypothèse est assez naturelle en modélisation grandes échelles.

#### Modèle dynamique complet

D'après les relations (17-20), le modèle dynamique complet s'écrit, avec $X = [\zeta, D]^T$ :

$$\frac{\partial}{\partial t} \begin{bmatrix} \zeta \\ D \end{bmatrix} + \underbrace{\begin{bmatrix} v \cdot \nabla \zeta - \nu \Delta \zeta \\ v \cdot \nabla D + D^2 - \nu \Delta D \end{bmatrix}}_{M(X)} = \nu_m \quad (21)$$

Ce modèle est non linéaire et son tangent linéaire autour de $X = \tilde{X} + dX$ se lit :

$$\begin{aligned} \left(\frac{\partial M}{\partial \tilde{X}}\right) (dX) &= \left[ \begin{array}{l} \mathbb{H}_G(\tilde{X}) \cdot \nabla d\zeta + \mathbb{H}_G(dX) \cdot \nabla \tilde{\zeta} - \nu \Delta d\zeta \\ \mathbb{H}_G(\tilde{X}) \cdot \nabla dD + \mathbb{H}_G(dX) \cdot \nabla \tilde{D} + 2\tilde{D}dD - \nu \Delta dD \end{array} \right] \\ &= \left[ \begin{array}{l} \mathbb{H}_G(\tilde{X}) \cdot \nabla + \nabla \tilde{\zeta} \cdot \mathbb{H}_G - \nu \Delta \\ \mathbb{H}_G(\tilde{X}) \cdot \nabla + \nabla \tilde{D} \cdot \mathbb{H}_G + 2\tilde{D} - \nu \Delta \end{array} \right] dX. \end{aligned} \quad (22)$$

La discrétisation de ce modèle requiert une attention particulière. Les termes advectifs de la forme $\bullet \cdot \nabla dX$ (où $\bullet$ est un vecteur) sont en effet susceptibles d'introduire un bruit numérique non négligeable. Afin de garantir une discrétisation stable de ces termes, des schémas conservatifs doivent être utilisés et nous avons eu recours à des schémas TVD ("Total Variation Diminishing") qui empêchent un accroissement des erreurs numériques au cours du temps (Kurganov et Levy, 2000; Kurganov et Tadmor, 2000; Levy et Tadmor, 1997). En ce qui concerne l'intégration temporelle, elle est réalisée par un schéma de Runge-Kutta au troisième ordre, qui respecte également les propriétés TVD (Giles, 2000).
A partir de ce modèle tangent linéaire, il faut dériver son adjoint requis par l'équation (13). Son expression analytique n'est pas triviale car cet opérateur fait intervenir des produits de convolution par le terme $\mathbb{H}_G$. Néanmoins, on peut démontrer dans l'espace de Fourier que l'adjoint de $\mathbb{H}_G$ est $-\mathbb{H}_G$ (voir (Papadakis et Memin, 2008) pour plus de détails). En ce qui concerne les autres opérateurs (advection, Laplacien), une fois créatisés et exprimés sous forme matricielle, nous exploitons le fait que les opérateurs adjoints associés sont obtenus en transposant les matrices associées aux opérateurs tangent linéaires. Plus de détails sur les opérateurs adjoints et leurs construction peuvent être trouvés dans (Talagrand et Courtier, 1987).

Passons à présent à la définition de l'opérateur d'observation.

### 2.3.3. Opérateur d'observation

Le but ici est de définir un opérateur d'observation $\mathbb{H}$ qui établit un lien entre l'espace d'état $X = [D, \zeta]^T$ et les observations images $\mathcal{Y}$. Pour cela, on peut s'appuyer sur l'ECMA vu dans la section précédente :

$$ \frac{\partial I}{\partial t} + \nabla I \cdot v = 0. \qquad (23) $$

Cependant dans le cas de fluides atmosphériques, on peut également de baser sur la loi de conservation de la masse, qui s'exprime, pour un fluide de densité $\rho$ :

$$ \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{V}) = 0, \qquad (24) $$

où $\mathbf{V}$ est la vitesse tridimensionnelle. Cette relation établit un lien entre la perte de densité du fluide et la divergence qu'il subit. En établissant l'hypothèse que l'intensité des image est liée à la densité du fluide (la validité de cette hypothèse est discutée dans (Corpetti et al., 2002, 2009), on obtient comme terme d'observation:

$$ \frac{\partial I}{\partial t} + \nabla I \cdot v + I \text{div} v = 0. \qquad (25) $$

L'ECMA correspond à cette dernière relation lorsque la divergence est nulle. Nous préférons ainsi s'appuyer sur ce modèle qui est plus générique pour définir le système d'observation. En notant que $v = \mathbb{H}_G(X)$ et $\text{div} v = D$, le système $\mathcal{Y} = \mathbb{H}(X)$ est construit à partir des images d'entrées et de la relation (25) de la manière suivante:

$$ \mathcal{Y} = \frac{\partial I}{\partial t} \quad \text{and} \quad \mathbb{H}(X) = (-\nabla I^T \mathbb{H}_G - [0, I]) X. \qquad (26) $$

Son adjoint $(\partial_X \mathbb{H})^*$ s'écrit :

$$ (\partial_X \mathbb{H})^* = \mathbb{H}_G \nabla I - [0, I]^T. \qquad (27) $$

Il convient de noter que cet opérateur d'observation est soumis au problème de l'ouverture, mentionné en introduction. Ainsi dans certaines situations difficiles, comme les zones homogènes, ce modèle n'est pas valide et cela doit être spécifié dans la matrice de covariance d'erreur associée $R$. Dans le cas où de telles situations seraient

fréquentes, il est également possible de supposer que le champ de déplacement est homogène dans un certain voisinage autour de chaque position $x$. Dans ce cas, cela revient à convoluer les relations (26) et (27) par un noyau Gaussien. Passons à présent à la définition des matrices de covariances.

### 2.3.4. Matrices de covariances

#### Observations

Deux manières différentes peuvent être utilisées pour fixer la matrice de covariance d'erreur $R$ liée aux observations. La première s'appuie sur l'importance du gradient de la fonction de luminance $\|\nabla I(x)\|$, dont on sait qu'elle est étroitement liée à la qualité du modèle d'observation (les zones à fort gradient ne sont pas soumises au problème de l'ouverture). Ainsi, $R$ peut se définir par :

$$ R(x) = R_{min} + (R_{max} - R_{min}) \exp - \frac{\|\nabla I(x)\|^2}{\sigma_{obs}^2}, \qquad (28) $$

où $R_{min}$, $R_{min}$ et $\sigma_{obs}$ sont fixés par l'utilisateur et dépendent des données. L'incertitude associée au modèle d'observation, liée à la valeur de $R(x)$, est ainsi maximale (resp. minimale) quand $\|\nabla I(x)\|$ est faible (resp. élevée). Dans un cadre plus générique où on utilise un modèle d'observation plus évolué (donc moins sensible à l'importance du gradient), la matrice $R^{-1}$ peut également se définir par :

$$ R^{-1}(x) = R_{max} \exp - \frac{[\mathcal{Y} - \mathbb{H}(X)]^2}{\sigma_{obs}^2}. \qquad (29) $$

Comme cela est démontré dans (Corpetti et al., 2009), cette pénalisation revient à considérer une norme robuste dans le premier terme de la fonction de coût en (9). Une telle fonction robuste permet d'accorder une importance moindre aux données "aberrantes" du terme d'observation $[\mathcal{Y} - \mathbb{H}(X)]$ (communément nommées *outliers* en statistique robuste, voir par exemple (Huber, 1981; Geman et Reynolds, 1992; Delanay et Bresler, 1998)). Dans notre cas, cela permet de rejeter les points où le modèle d'observation n'est pas complètement valide.

#### Condition initiale

Nous avons choisi de représenter la matrice de covariance $B$ associée à la condition initiale par $B = \text{Id} - \exp\{-|X(t_0) - X_0|/\sigma_B^2\}$ où $\sigma_B$ est un paramètre à fixer. Elle est liée à l'adéquation entre la solution courante au temps $t_0$ et la condition initiale.

#### Modèle dynamique

Etant donné que nous n'avons pas d'a priori sur les zones où le modèle dynamique présenté dans la relation (21) est valide ou non, la matrice $Q$ a été fixée comme étant constante diagonale, c'est-à-dire qu'on accorde une incertitude identique à tous les points de l'image et qu'on ne suppose pas de corrélations des erreurs entre les différents points.
Maintenant que l'espace d'état, le modèle dynamique, l'opérateur d'observation et les matrices de covariance des erreurs associées sont définis, tous les ingrédients sont réunis pour le système d'assimilation de données. La section suivante présente quelques résultats.

## 2.4. Résultats sur des images MSG

Nous avons appliqué cette technique à deux séquences météorologiques. La première (haut de la figure 2) date du 9 Octobre 2005 et représente un clone. Pour trois instants de la séquence, nous avons représenté le champ de déplacement estimé superposé à l'image correspondante (première ligne). La seconde ligne représente la composante de vorticité associée tandis que la troisième ligne affiche les cartes de vorticités obtenues avec un estimateur de mouvement dédié aux fluides mais dans lequel aucun lien temporel n'est imposé (cet estimateur est présenté dans (Corpetti et al., 2002)). Sur ces deux dernières lignes, il est intéressant d'observer que la structure en rotation apparaît clairement. Cependant, comme l'image contient du bruit, on remarque que même un estimateur de mouvement dédié aux phénomènes fluides fournit des cartes de vorticité parfois bruitées, spatialement et temporellement. Ce bruit est en grande partie retiré à l'aide de la technique d'assimilation présentée dans cet article, grâce au lien temporel imposé par le modèle dynamique. Des résultats plus complets avec évaluations qualitatives et quantitative peuvent être obtenus dans (Papadakis et al., 2007).

La seconde séquence, en bas de la figure 2, a été acquise en 1998 et représente une dépression dans le sud de l'Europe. Nous avons bruité quelques images d'entrée, comme on peut le voir dans la figure 2 (j–m). Les champs de déplacements obtenus par assimilation de données sont représentés sur les figures 2 (n–q) tandis que la dernière ligne (Figure 2 (r–u)) représente les champs de mouvements associés sans cohérence temporelle (technique présentée dans (Corpetti et al., 2002)). Sur ces images, il est évident que l'assimilation permet d'améliorer grandement les résultats.

Ces quelques résultats illustrent l'intérêt de l'utilisation de l'assimilation de données pour estimer, à partir d'une séquence d'images, les champs denses de déplacements associés. D'autres résultats exploitant le formalisme de l'assimilation de données, pour de l'estimation du mouvement (pour des images de pression en particulier) ou encore pour d'autres tâches liées aux géosciences comme la reconstruction de données par peuvent être consultés dans (Ba et al., 2011, 2012; Corpetti et al., 2007, 2009; Papadakis et al., 2007; Papadakis et Memin, 2008). Passons à présent au deuxième outil que nous avons proposé et qui concerne le suivi de courbes.

# 3. Suivi de courbes et application aux cellules convectives

## 3.1. Détection de masses nuageuses : problématique

Suivre des courbes dans une image est un problème ancien qui se décline sous plusieurs formes, selon les situations (une seule ou plusieurs courbe(s), ouverte(s) ou fermée(s), possibilité ou non de fusion, scission, ...). Dans le cas des structures nuageuses, la notion de contours des objets à suivre est plus floue que pour des objets rigides et nous sommes par ailleurs régulièrement confrontés à des disparitions, éclatements, regroupements, ..., de masses nuageuses. Tel est le cas par exemple des cellules convectives. Elles peuvent s'organiser de façon large au sein de la ZCIT ou des secteurs de mousson, parfois sur des systèmes intermédiaires (méso-échelle) mais peuvent aussi résulter de dynamiques diurnes, par exemple au dessus des grandes forêts tropicales ; leur cycle de vie peut alors être assez bref. La contribution globale de ces formations nuageuses est pourtant décisive en terme de pluviométrie dans les régions tropicales (80% des précipitations en Amazonie brésilienne). Elles peuvent également être responsables de violentes pluies, tornades, tempêtes, ... Au sein de la ZCIT notamment et en Afrique centrale, ces cellules génèrent la plupart de la pluie durant la mousson. Elles sont également indirectement liées aux sécheresses et aux inondations et cela a des conséquences parfois désastreuses. Le suivi de telles cellules fait ainsi l'objet de nombreuses recherches (voir par exemple les travaux de (Machado et al., 2004)), généralement en exploitant la haute fréquence d'acquisition d'images par les satellites géostationnaires. Avec de telles études, les phases d'initiation ou de dissipation de ces cellules ont parfois pu être mises en relation avec les hétérogénéités de surface (relief, occupation du sol). Du point de vue fonctionnel, cela a également permis de caractériser les cellules convectives par deux régions, représentant des zones avec des températures caractéristiques : le *coeur*, qui contient les plus fortes pluies, et la *zone d'influence*, caractéristique de l'étendue de la cellule. Afin de systématiser ce genre d'approche, nous avons travaillé sur le développement d'outils automatiques pour extraire et suivre dans le temps les contours des cellules convectives.

En analyse d'images, une méthode efficace pour détecter des contours complexes consiste à faire évoluer d'une manière judicieuse des équations aux dérivées partielles définies sur des courbes (paramétriques ou non, voir par exemple les références (Caselles et al., 1997; Chan et Vese, 2001; Kass et al., 1988)). Les fonctionnelles associées incluent en général un terme d'attaché aux données basé sur les propriétés photométriques des objets à suivre, auquel on ajoute un terme de lissage caractérisant la régularité désirée du contour. En pratique cependant, les équations résultantes sont non linéaires et les contours extraits dépendent fortement de la condition initiale, empêchant de
gérer des cas compliqués tels que la fusion ou scission par exemple. Parmi les différentes extensions proposées, les méthodes Eulériennes d'ensemble de niveaux (communément appelés "level-set", voir (Osher et Sethian, 1988; Sethian, 1999)) ont été spécifiquement proposées pour éviter ces situations difficiles. Dans ce cas, les contours à extraire sont représentés par le niveau zéro d'une surface implicite de dimension supérieure. Le problème revient alors à décrire l'évolution de cette surface implicite de manière à gérer correctement les changements topologiques. Une illustration du fonctionnement de ces techniques est visible sur la figure 3.

De telles approches ont déjà été appliquées dans un contexte météorologique (voir (Cohen et Herlin, 1998; Papin et al., 2000; Yahia et Berroir, 1998)). Cependant elles se concentrent sur la segmentation des masses nuageuses de manière indépendante en chaque image et aucune cohérence temporelle n'est garantie. Cela est parfois préjudiciable lorsque l'on est en présence de mouvements violents engendrant des déformations massives comme c'est le cas des cellules convectives et d'un grand nombre d'autres phénomènes atmosphériques. Ainsi, l'analyse du comportement de ces cellules et leur prévision est donc d'un intérêt majeur.

Nous avons proposé d'exploiter le formalisme de l'assimilation variationnelle de donnée pour effectuer le suivi simultané de deux courbes imbriquées correspondant aux deux régions d'intérêts des cellules convectives (coeur et zone d'influence). Ceci est présenté dans la section suivante.

## 3.2. Suivi de cellules convectives par assimilation de données

Etant donné que le coeur et la zone d'influence des cellules convectives possèdent une température caractéristique, nous utilisons des cartes de température comme données d'entrées, celles-ci étant extraites à partir d'une combinaison de canaux MSG (voir par exemple (Thomas et al., 2010)). Nous devons ainsi définir l'espace d'état, le modèle dynamique, l'opérateur d'observation et les matrices de covariance des erreurs.

### 3.2.1. Espace d'état

En utilisant le formalisme des ensembles de niveaux, une courbe à l'instant $t$, $C_t(p) = C(p, t) : [a, b] \times \mathbb{R} \rightarrow \mathbb{R}^2$ est décrite de manière implicite par le niveau zéro d'une fonction scalaire $\phi(x, t) : \Omega \times \mathbb{R} \rightarrow \mathbb{R}$:

$$C_t(.) = C(., t) = \{x \in \Omega \mid \phi(x, t) = 0\}, \quad (30)$$

où $\Omega$ est le domaine image. Par convention, cette surface est généralement choisie de telle manière à ce que les points à l'intérieur (resp. à l'extérieur) d'un objet aient des valeurs positives (resp. négatives). Un choix consiste à prendre la distance signée par rapport aux contours d'un objet mais toute autre surface dont le niveau zéro correspond aux contours d'intérêts est possible. Dans cette application, afin de représenter deux courbes emboîtées,

l'espace d'état $X$ utilisé est composé de deux surfaces implicites :

* La première, $\phi^e$, est attachée aux contours externes (la zone d'influence) des cellules convectives, $C^e(t)$. Ces courbes sont associées à un niveau de température $T_e$;

* La seconde, $\phi^c$, représente le coeur de la région de chaque cellule, $C^c(t)$. Elle est reliée à une température plus froide $T_c$.

La première surface est définie sur l'ensemble du domaine image $\Omega$: $\phi^e(x, t) : \Omega \times \mathbb{R} \rightarrow \mathbb{R}$ tandis que la seconde est définie à l'intérieur du domaine encapsulé par le contour des cellules $\mathcal{I} = \{x \mid \phi^e(x, t) > 0\}$ (i.e. le domaine délimité par le niveau 0 de $\phi^e$), $\phi^c(x, t) : \mathcal{I} \times \mathbb{R} \rightarrow \mathbb{R}$. Ainsi, l'espace d'état composé des deux surfaces s'écrit : $X = [\phi^e, \phi^c]^T$.

### 3.2.2. Modèle dynamique

Afin de se doter d'un modèle dynamique pour $X$, nous supposons que chaque courbe est soumise à un champ de vitesse et propagée par celui-ci, à une incertitude près le long de sa normale et le long de sa tangente. Pour chacune des deux courbes, cela peut s'écrire :

$$dC_t^i = (w^i \cdot n^i)n^i dt + \sigma_n n^i dB_t^n + \sigma_\tau n^{i\perp} dB_t^\tau. \quad (31)$$

Ici, les termes $n^i dB^n$ et $n^{i\perp} dB^\tau$ sont deux mouvements Browniens indépendants dirigés le long de la tangente et de la normale à chaque courbe et $\sigma_\tau$, $\sigma_n$ sont les incertitudes associées. Dans ce modèle, la courbe est transportée par un mouvement déterministe (premier terme de la relation (31)) et cela est vrai à un bruit Gaussien isotrope près (le long de la normale et de la tangente) dont les covariances croient linéairement avec le temps (deuxièmes et troisièmes termes de (31)). Comme mentionné précédemment, la position de la courbe est définie par l'évolution de la surface implicite $\phi^i(x, t)$. Par conséquent, cette surface dépend du processus stochastique $C_t$ et sa différentielle est obtenue en utilisant des règles de calcul stochastique et en particulier la formule d'Ito, mentionnée dans la section précédente. Nous obtenons ainsi l'équation différentielle stochastique suivante (où nous avons retiré les indices de courbes $i = \{e, c\}$ pour plus de lisibilité):

$$\left( \partial_t \phi + w \cdot \nabla \phi + \frac{1}{2} \sum_{\{x_i, x_j\} = \{x, y\}} \frac{\partial^2 \phi}{\partial x_i \partial x_j} d\langle C_t^{x_i}, C_t^{x_j} \rangle \right) dt + \sigma_n |\nabla \phi| d\hat{B}_t^n = 0, \quad (32)$$

où $C_t^x$ et $C_t^y$ correspondent aux coordonnées des points de la courbe dans le plan Euclidien au temps $t$. Nous avons ici introduit les déplacements $w(x)$ et les mouvements Browniens $\hat{B}_t^n(x)$, qui sont des extensions dans le plan image des bruits associés aux courbes. En ce qui concerne la vitesse, nous avons considéré un champ de déplacement unique pour les deux courbes, ce champ
pouvant être en pratique mesuré sur les images à l'aide d'un estimateur de mouvement externe (par exemple issu des techniques de corrélation) ou encore issu de la technique présentée dans la section précédente, ce qui est le cas des exemples de cet article. Afin de dériver un modèle dynamique que suit chaque surface $\phi$ à partir de la relation (32), nous développons les variations quadra-tiques de la forme $\langle C_t^{x_i}, C_t^{x_j} \rangle$, ce qui donne :

$$ \underbrace{\left( \frac{\partial \phi}{\partial t} + \mathbf{M}(\phi) \right)} \text{Modèle dynamique} \quad dt + \sigma_n |\nabla \phi| d\hat{B} = 0, \quad (33) $$

avec

$$ \mathbf{M}(\phi) = w \cdot \nabla \phi + \frac{\sigma_{\tau}^2}{2} |\nabla \phi| + \frac{\sigma_n^2}{2|\nabla \phi|^2} \nabla \phi^T \nabla^2 \phi \nabla \phi, \quad (34) $$

où $\nabla^2 \phi$ représente la matrice Hessienne et $\kappa$ est la courbure moyenne : $\kappa = \text{curv}(\phi) = \frac{1}{\|\nabla\phi\|} (\Delta\phi - \nabla\phi^T \nabla^2\phi \nabla\phi)$. L'espérance temporelle de la relation (33) est sa partie déterministe et nous la considérons comme notre modèle d'évolution pour les courbes $\phi^i$, $i = \{e, c\}$. Tous les détails concernant ce modèle dynamique, avec en particulier la manière d'estimer les certitudes ($\sigma_n, \sigma_\tau$), les modèles adjoints et les détails de l'implémentation peuvent se trouver dans (Thomas et al., 2009, 2010). Passons à présent à l'opérateur d'observation.

### 3.2.3. Opérateur d'observation

Dans le cas des cellules convectives, il n'est pas possible de construire un opérateur d'observation qui s'appuie sur un histogramme de référence, comme cela peut se faire en suivi d'objets plus conventionnels. En effet, l'organisation spatiale de la température le long du cycle de vie des cellules possède une forte variabilité interne, ce qui n'assure pas de cohérence temporelle des histogrammes associés. En revanche, on peut raisonnablement supposer, comme cela se fait pour l'étude de phénomènes convectifs, que la zone du coeur (resp. d'influence) vit pour des températures en dessous d'un certain seuil $T_c$ (resp. $T_c$). Ainsi, la fonction $\phi$, dont le niveau zéro correspond au contour à extraire, doit être telle que :

$$ H(I(x) - T_{\bullet})(1 - H(\phi^{\bullet}(\mathbf{x}))) + H(T_{\bullet} - I(x))H(\phi^{\bullet}(\mathbf{x})) = \epsilon, \quad (35) $$

où $\bullet = \{e, c\}$ et $H$ est la fonction de Heaviside ($H(x) = 0$ si $x < 0$ et $H(x) = 1$ sinon). Cette relation est nulle si $I(x) > T_{\bullet}$ et $\phi(\mathbf{x}) > 0$ ou si $I(x) < T_{\bullet}$ et $\phi(\mathbf{x}) < 0$. Cette relation encourage ainsi l'intérieur de la courbe à être sous un seuil donné et l'extérieur de la courbe au dessus de ce seuil. Dans les autres cas cette fonction est toujours positive. Ainsi, à partir de (35), on peut définir un modèle d'observation $\mathcal{Y} = \mathbb{H}(\phi)$ par :

$$ \left\{ \begin{array}{l} \mathcal{Y} = 0 \\ \mathbb{H}(\phi) = H(I(x) - T_{\bullet})(1 - H(\phi(\mathbf{x}))) \quad (36) \\ \qquad + H(T_{\bullet} - I(x))H(\phi(\mathbf{x})) \end{array} \right. $$

L'opérateur $\mathbb{H}$ est non linéaire par rapport à la variable d'état $\phi$. Son tangent linéaire, l'opérateur gradient associé ainsi que la définition des matrices de covariances peuvent se trouver dans (Thomas et al., 2009, 2010).

### 3.2.4. Quelques résultats

Nous avons testé cette méthode sur une première séquence d'images qui met en jeu l'évolution d'une importante cellule convective. Les résultats sur visibles sur la figure 4. La première ligne représente quatre images de la séquence où la propagation du contours initial par le modèle dynamique est superposée (le contour en bleu représente la zone d'influence de la cellule et le contours en rouge représente son coeur). La seconde ligne représente les résultats après le processus d'assimilation de données. Sur la première ligne de cette figure, on peut observer que l'intégration du contours initial dans le modèle dynamique fournit déjà une bonne approximation de la solution finale pour la zone d'influence. Dans cet exemple particulier, le développement de la cellule est progressif et assez linéaire. Ainsi le modèle dynamique linéarisé fournit des contours assez proches de la solution finale. A l'inverse, le coeur de la cellule se dissipe lentement au cours de la séquence et possède un comportement plus chaotique. L'intérêt du processus d'assimilation apparaît clairement sur la seconde ligne de la figure 4 puisque les courbes estimées, liées au coeur et à la zone d'influence, semblent cohérentes. La seconde expérience représente une scène plus complexe qui contient plusieurs cellules à différents degrés d'évolution. Par ailleurs, étant donné qu'aucun phénomène convectif n'est présent dans la première image, il n'y a aucune courbe pour la condition initiale et par conséquent, la première étape de la méthode (intégration du modèle dynamique) ne fournit aucune courbe pour toutes les images de la séquence. Malgré cette initialisation difficile, nous pouvons observer sur la figure 5 que le processus d'assimilation permet d'estimer correctement les concours associés à la zone d'influence et à la partie active de chaque cellule. Ces deux expériences démontrent la capacité de la technique d'assimilation de données proposée à recouvrer les contours de masses nuageuses, déformables et à durée de vie variable. Passons à présent à la troisième application qui concerne le suivi de fronts de brise de mer.

## 4. Détection de fronts

### 4.1. Problématique

La détection de fronts est un point critique dans de nombreux domaines liés aux géosciences. En effet, un front est en général le marqueur d'une séparation entre deux zones bien spécifiques : par exemple, le front marin représente une transition entre courants chauds/froids, le front anticyclonique localise deux masses d'aires entrant en contact ou encore le front de brise de mer matérialise la limite de pénétration de la brise dans les terres. Localiser un front dans une image est un problème délicat: il est représenté par une courbe ouverte, ce qui d'une
part empêche l'emploi de techniques efficaces comme les ensembles de niveaux vus précédemment et d'autre part cela pose un certain nombre de difficultés liées aux conditions aux bords.

Ici, nous nous focalisons sur la détection du front de brise de mer à partir d'images MSG. La brise de mer est un vent prenant naissance lorsqu'un gradient important de température entre l'océan et la terre apparaît. Une circulation allant de la mer vers le continent et prenant fin au front de brise apparaît. Au niveau du front, l'air se condense et cela provoque le développement de nuages cumuliformes parallèles. Cette ligne de nuages matérialise le front de brise de mer, comme on peut le voir sur la figure 7 (a). L'identification de cette formation des nuages dans les images MSG permet de déterminer l'importance de la pénétration de la brise dans les terres et par ailleurs, sa localisation en fonction des saisons ou encore les liens entre l'importance de la pénétration du front de brise avec le climat local sont des phénomènes qui ne sont à l'heure actuelle pas entièrement maîtrisés.

Du point de vue de l'analyse d'images, le front n'apparaît pas de manière nette mais il est plutôt caractérisé par la frontière entre deux motifs texturés composés d'un côté d'une zone sans nuages et de l'autre de nuages cumuliformes. Par conséquent, les techniques usuelles de détection/caractérisation de fronts, comme par exemple les techniques de type contrast-invariant image representations (cf. (Alvarez et al., 1993; Monasse et Guichard, 2000)) ou les approches basées sur les coefficients de Lyapunov, sont inefficaces ici. La principale contribution que nous avons proposée se résume en un détecteur de front P original à partir duquel une technique de segmentation a ensuite été appliquée. Celle-ci s'appuie sur les contours actifs introduits dans (Kass et al., 1988). Dans le cas particulier où les points extrêmes de la ligne de front sont connus (cela est parfois le cas dans certains endroits géographiques), on peut également utiliser la méthode des chemins minimaux présentée dans (Cohen et Kimmel, 1997). La carte de potentiel P (i.e. le détecteur de front) proposée en tout point x de l'image s'écrit :

$$P(x) = \max_{\theta} \left\{ \begin{array}{l} \kappa(R_{\theta_1}, R_f) \overline{\kappa(R_{\theta_2}, R_f)} + \\ \overline{\kappa(R_{\theta_1}, R_f)} \kappa(R_{\theta_2}, R_f) \end{array} \right\}, \text{où} \quad (37)$$

* $\theta$ est un angle;

* $R_f$ est une région d'apprentissage qui caractérise les textures formées par les nuages cumuliformes;

* $R_{\theta_1}$ et $R_{\theta_2}$ sont deux régions adjacentes autour du point x et délimitées par un angle $\theta$ (voir l'illustration dans la figure 6);

* $\kappa(R_{\theta_\ell}, R_f)$, $\ell = \{1, 2\}$ est un critère de similarité entre deux régions $R_{\theta_\ell}$ et $R_f$ et

* $\overline{\kappa(R_{\theta_\ell}, R_f)} = 1 - \kappa(R_{\theta_\ell}, R_f)$ est une croyance sur le fait que la région $R_{\theta_\ell}$ est différente de $R_f$.

En tout point $x$, le potentiel est maximal lorsque seulement une des deux régions $R_{\theta_1}$ et $R_{\theta_2}$ adjacentes à $x$ appartient à $R_f$. Il est donc lié à la présence/absence du front de brise. En pratique, le critère $\kappa$ de similarité utilisé est un critère de comparaison de textures s'appuyant sur les ondelettes et la théorie des évidences défini dans (Corpetti et Planchon, 2011).

## 4.2. Quelques résultats

La figure 7 (a) est une image issue du canal visible de MSG le long de la côte brésilienne où nous avons superposé le contour initial sur la ligne de côte (tracé à la main). On peut observer le front de brise au nord de l'image et en figure 7 (b) un détail de la figure 7 (a) pour observer précisément le front texturé. Les fronts sont détectés en intégrant la fonction P définie dans la section précédente dans une technique de contours actifs. Afin d'évaluer l'apport du détecteur de front présenté dans la relation (37), nous avons comparé les résultats avec la même technique mais qui utilise un détecteur de contours classique basé sur la norme des gradients de la fonction de luminance de la forme $P(x) = \lambda \|\nabla I(x)\|$.

Les figures 7 (c) et 7 (d) représentent le front obtenu quand les points extrêmes sont respectivement inconnus et connus tandis que la figure 7 (e) correspond au potentiel obtenu avec la relation $P(x) = \lambda \|\nabla I(x)\|$. En utilisant la fonction de potentiel P de la relation (37), les fronts obtenus sont lorsque les points extrêmes sont inconnus ou non sont respectivement visibles dans les figures 7(f,g) tandis que la fonction de potentiel est représentée en 7(h). Sur ces exemples, l'intérêt de la fonction de potentiel P est bien mis en évidence car il identifie nettement la ligne de front, ce qui n'est pas le cas du potentiel issu de la figure 7 (e) qui met en évidence les zones de forts gradients mais qui ne correspond pas à un front réel. Ces résultats démontrent l'intérêt du détecteur de fronts qui s'appuie sur la relation (37). Pour de détails sont disponibles dans (Corpetti et Planchon, 2011).

# 5. Conclusion

Dans cet article, nous avons passé en revue trois méthodes développées pour analyser les images MSG : estimation du mouvement, suivi de courbes fermées et ouvertes. Les deux premières techniques s'appuient sur l'assimilation variationnelle de données qui est un cadre particulièrement adapté à la manipulation d'espaces d'états de grandes dimension soumis à des dynamiques imparfaitement modélisées. L'approche pour détecter les fronts s'appuie quant à elle sur un détecteur original basé sur la comparaison de zones texturées. Ces développements permettent de fournir des outils adaptés à la spécificité des images MSG : données en grandes dimension avec une résolution spatiale et temporelle améliorée par rapport à la première génération. Cela autorise la visualisation plus fine et l'analyse plus précise de phénomènes atmosphériques et météorologiques.
Références

Alvarez, L., Guichard, F., Lions, P.-L., Morel, J.-M., 1993. Axioms and fundamental equations of image processing. Archive for Rational Mechanics and Analysis 123 (3), 199–257.

Ba, S., Corpetti, T., Fablet, R., September 2011. Multi-resolution missing data interpolation in SST image series. Dans: IEEE International Conference on Image Processing, ICIP '11. Brussels, Belgium, pp. 1517–1520.

Ba, S., Fablet, R., Corpetti, T., Chapron, B., March 2012. Reconstruction of velocities and SST with the variational method. Dans: Ifremer-ESA Globcurrent workshop. Brest, France.

Bennett, A., 1992. Inverse Methods in Physical Oceanography. Cambridge University Press, U.K.

Caselles, V., Kimmel, R., Sapiro, G., 1997. Geodesic Active Contours. International Journal of Computer Vision 22 (1).

Chan, T., Vese, L., 2001. Active contours whitout edges. IEEE Transactions on Image Processing 10 (2), 266–277.

Cohen, I., Herlin, I., 1998. Tracking meteorological structures through curve(s) matching using geodesic paths. Dans: Proc of 6th IEEE International Conference on Computer Vision. Bombay, India, pp. 396–401.

Cohen, L., Kimmel, R., 1997. Global Minimum for Active Contour Models: A Minimal Path Approach. International Journal of Computer Vision 24 (1), 57–78.

Corpetti, T., Héas, P., Mémin, E., Papadakis, N., 2009. Pressure image assimilation for atmospheric motion estimation. Tellus Series A: Dynamic Meteorology and Oceanography 61 (1), 160–178.

Corpetti, T., Mémin, E., Pérez, P., 2002. Dense Estimation of Fluid Flows. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (3), 365–380.

Corpetti, T., Mémin, E., Pérez, P., 2003. Extraction of Singular Points from Dense Motion Fields: an Analytic Approach. Journal of Mathematical Imaging and Vision 19 (3), 175–198.

Corpetti, T., Papadakis, N., Mémin, E., 2007. Dense estimation of motion fields on meteosat second generation images using a dynamical consistency. Dans: IEEE International Geoscience and Remote Sensing Symposium. Barcelona, Spain, pp. 4749–4752.

Corpetti, T., Planchon, O., 2011. Front detection on satellite images based on wavelet and evidence theory: application to the sea breeze fronts. Remote Sensing of Environment 115 (2), 306–324.

Delanay, A., Bresler, Y., 1998. Globally convergent edge-preserving regularized reconstruction: an application to limited-angle tomography. IEEE Transactions on Image Processing 7 (2), 204–221.

Geman, D., Reynolds, G., 1992. Constrained Restoration and The Recovery Of Discontinuities. IEEE Transactions on Pattern Analysis and Machine Intelligence 14 (3), 367–383.

Giles, M., 2000. On the use of Runge-Kutta time-marching and multigrid for the solution of steady adjoint equations. Rapport technique 00/10, Oxford University Computing Laboratory.

Holmlund, K., 1998. The utilization of statistical properties of satellite-derived atmospheric motion vectors to derive quality indicators. Weather and Forecasting 13 (4), 1093–1104.

Horn, B., Schunck, B., 1981. Determining optical flow. Artificial Intelligence 17 (1-3), 185–203.

Huber, P., 1981. Robust Statistics. John Wiley & Sons.

Kass, M., Witkin, A., Terzopoulos, D., 1988. Snakes : Active Contour Models. International Journal of Computer Vision 1 (4), 321–331.

Kurganov, A., Levy, D., 2000. A Third-Order Semidiscrete Central Scheme for Conservation Laws and Convection-Diffusion Equations. SIAM Journal on Scientific Computing 22 (4), 1461–1488.

Kurganov, A., Tadmor, E., 2000. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. Journal of Computational Physics 160 (1), 241–282.

Le Dimet, F.-X., Talagrand, O., 1986. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus, 97–110.

Levy, D., Tadmor, E., 1997. Non-oscillatory Central Schemes for the Incompressible 2-D Euler Equations. Mathematical Research Letters 4, 321–340.

Lions, J.-L., 1971. Optimal control of systems governed by PDEs. Springer-Verlag.

Machado, L., Laurent, H., Dessay, N., Miranda, I., 2004. Seasonal and diurnal variability of convection over the Amazonia: a comparison of different vegetation types and large scale forcing. Theoretical Applied Climatology 78, 61–77.

Mansour, N., Ferziger, J., Reynolds, W., 1978. Large-eddy simulation of a turbulent mixing layer. Rapport technique, Report TF-11, Thermosciences Div., Dept. of Mech. Eng., Standford University, USA.

Menzel, W. P., Smith, W. L., Stewart, T., 1983. Improved cloud motion wind vector and altitude assignment using vas. Journal of Applied Meteorology 22, 377–384.

Monasse, P., Guichard, F., 2000. Fast computation of a contrast-invariant image representation. IEEE Transactions on Image Processing 9 (5), 860–872.

Oksendal, B., 1998. Stochastic differential equations. Springer-Verlag.

Osher, S., Sethian, J., 1988. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics 79 (1), 12–49.

Papadakis, N., 2007. Assimilation de données images : application au suivi de courbes et de champs de vecteurs. Thèse de doctorat, University Rennes 1, France.

Papadakis, N., Corpetti, T., Mémin, E., 2007. Dynamically consistent optical flow estimation. Dans: IEEE International Conference on Computer Vision. Rio de Janeiro, Brasil, pp. 1–7.

Papadakis, N., Mémin, E., 2008. Variational Assimilation of Fluid Motion from Image Sequence. SIAM Journal on Imaging Sciences 1 (4), 343.

Papin, C., Bouthemy, P., Mémin, E., Rochard, G., 2000. Tracking and Characterization of Highly Deformable Cloud Structures. Dans: European Conference on Computer Vision. Vol. II. Dublin, Ireland, pp. 428–442.

Schmetz, J., Geijo, C., Menzel, W. P., Strabala, K., Van De Berg, L., Holmlund, K., Tjemkes, S., 1995. Satellite observations of upper tropospheric relative humidity, clouds and wind field divergence. Contributions to Atmospheric Physics 68 (4), 345–357.

Schmetz, J., Holmlund, K., Hoffman, J., Strauss, B., Mason, B., Gaertner, V., Koch, A., Berg, L. V. D., 1993. Operational cloud-motion winds from meteosat infrared images. Journal of Applied Meteorology 32 (7), 1206–1225.

Sethian, J., 1999. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science. Cambridge University Press.

Talagrand, O., 1997. Assimilation of observations, an introduction. Journal of the Meteorological Society of Japan 75, 191–209.

Talagrand, O., Courtier, P., 1987. Variational assimilation of meteorological observations with the adjoint vorticity. I: Theory. Quarterly Journal of the Royal Meteorological Society 113 (478), 1311–1328.

Taylor, G., 1932. The transport of vorticity and heat through fluids in turbulent motion. Dans: Proc London Math Soc. Ser A. pp. 151–421.

Thomas, C., Corpetti, T., Mémin, E., 2009. Data Assimilation for Convective Cells Tracking in MSG images. Dans: IEEE International Geoscience and Remote Sensing Symposium.
Vol. 2. Cape Town, South Africa, pp. 813–816.

Thomas, C., Corpetti, T., Mémin, E., 2010. Data Assimilation for Convective Cells Tracking on Meteorological Image Sequences. IEEE Transactions on Geosciences and Remote Sensing 48 (8), 3162–3177.

Vidard, P., Blayo, E., Le Dimet, F.-X., Piacentini, A., 2000. 4D Variational Data Analysis with Imperfect Model. Flow, Turbulence and Combustion 65 (3-4), 489–504.

Yahia, H., Berroir, J., 1998. Segmentation of deformable templates with level sets characterized by particle systems. Dans: International Conference on Pattern Recognition. Brisbane, Australia, pp. 1421–1423.
<figure><img src="image_1.png" /><figcaption>**Figure 1: Décomposition de Helmholtz** (a): un champ de vecteurs; (b) sa composante solénoïdale; (c) sa composante irrotationnelle; (d) sa composante harmonique</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 2: Résultats sur des données météorologiques** HAUT: (a–c): trois images représentant un cyclone et les champs de déplacements associés; (d–f): la vorticité estimée correspondante; (g–i): La vorticité obtenue avec un estimateur de mouvement dédié aux fluides mais sans cohérence temporelle. BAS: séquence représentant une dépression: (j–m) 4 images de la séquence où certaines ont été bruitées; (n–q) champs de déplacements obtenus avec assimilation de données; (r–u) champs de déplacements obtenus sans assimilation de données.</figcaption></figure>
<figure><img src="image_3.png" /></figure>

**Figure 3: Principe des méthodes à ensemble de niveaux.** Le niveau zéro de la surface correspond aux contours à extraire. Par cette représentation il est possible de gérer les changements de topologie, comme la fusion des deux objets dans (a) en un seul objet dans (b).

<figure><img src="image_4.png" /></figure>

**Figure 4: Suivi de cellules convectives : séquence #1.** 4 images de la séquence. La première ligne correspond aux contours obtenus après la propagation du contour initial par le modèle dynamique ; la seconde ligne présente le résultat final après l'assimilation de données. La zone d'influence est délimitée par les courbes bleues tandis que la partie active est représentée par le contours en rouge.

<figure><img src="image_5.png" /></figure>

**Figure 5: Suivi de cellules convectives : séquence #2.** 4 images de la séquence où la zone d'influence est délimitée par les courbes bleues tandis que la partie active est représentée par le contours en rouge.

<figure><img src="image_6.png" /></figure>

**Figure 6: Regions $R_{\theta_1}$ et $R_{\theta_2}$ en fonction de l'angle $\theta$ en chaque point $x$ de l'image.**
<figure><img src="image_7.png" /><figcaption>**Figure 7: Détection du front de brise de mer:** (a): image satellite avec le contours initial en blanc et le front de brise réel en rouge; (b): un détail de l'image (a); (c,d): contours final obtenu avec une fonction de potentiel s'appuyant sur la norme des gradients lorsque les points extrêmes sont connus ou non; (f,g): idem que (c,d) avec une fonction de potentiel s'appuyant sur (37); (e,h): les cartes de potentiel correspondants.</figcaption></figure>