ASSIMAGE

Compte-rendu de réunion

Etienne Huot

31 Mars - 02 Avril 2004

Participants

• Etienne Mémin,
• Etienne Huot,
• Thomas Corpetti.

Ordre du jour

1. Présentation des programmes d'estimation de mouvement (Thomas Corpetti),
2. Discussion sur les différentes approches d'estimation et leur implémentation (Etienne Huot et Etienne Mémin),
3. Prise en main des programmes (Etienne Huot),
4. Exécution sur les données océanographiques du LODYC (Etienne Huot).

Estimation du mouvement

Problèmes d'instabilité

Des problèmes d'instabilité numérique peuvent apparaître à plusieurs niveaux du processus d'estimation du mouvement. En ce qui concerne les approches par EDP, les instabilités apparaissent:

1. lorsque les hypothèses de conservation et de régularité qui sont faites ne sont pas respectées par les données,
2. lorsque les échelles de temps sont trop grandes par rapport aux vitesses (grands déplacements),
3. lors du calcul des fonctions ${\mathrm {div\,}}$et ${\mathrm {curl\,}}$,
4. lors de l'utilisation de la régularité ${\mathrm {div\,}}-{\mathrm {curl\,}}$d'ordre deux.

A chacun de ces problèmes, Etienne Mémin propose une solution.

Modèle de données erroné

Lorsque les hypothèses de conservation ou de régularité ne sont pas respectées par les données, on peut proposer de remplacer les normes au sens $L_2$par des estimateurs robustes. Ainsi l'hypothèse de conservation de la luminance qui permet de déduire la fonctionnelle suivante:

$$E_{obs} = \int_\Omega \vert\frac{\displaystyle \partial I}{\displaystyle \partial t} + \mathbf {w}.\nabla I\vert^2 dx dy,$$

s'écrit:

$$E_{obs} = \int_\Omega \Psi_1(\frac{\displaystyle \partial I}{\displaystyle \partial t} + \mathbf {w}.\nabla I) dxdy.$$

où $\Psi_1$est un estimateur robuste, c'est-à-dire une fonction semi-quadratique qui possède les propriétés suivantes:

1. $\Psi_1$est croissante sur ${\rm I\!R}^+$,
2. $\lim_{t\rightarrow \infty} \Psi_1'(t) < \infty$et
3. $\Psi_1(\sqrt(t))$est concave.

C'est-à-dire que nous obtenons:

$$\forall t, \Psi_1(t) = \min_z \{zt^2 + \psi_1(z) \} \mbox{~ ... ...} = \arg \min_z \{zt^2 + \psi_1(z) \} = \frac{\Psi_1(t)}{2t},$$

où $\psi(t)$est une fonction strictement convexe définie explicitement à partir de $\Psi_1$. L'utilisation d'un tel estimateur permet de proposer une fonction de coût robuste, en conservant des propriétés intéressantes du point de vue de la minimisation puisqu'elle peut se faire de manière alternée. En effet, $\min_w(\Psi_1(g(w))$est remplacé par le problème auxiliaire$\min_{z,w} zg^2(w) + \psi(w)$, qui peut-être résolu par une technique de moindres carrés pondérés-itérés. La même technique est utilisée pour le terme de régularité, où:

$$E_{reg} = \int_\Omega (\vert\nabla u\vert^2 + \vert\nabla v\vert^2) dxdy$$

est remplacé par

$$E_{reg} = \int_\Omega [\Psi_2(\vert\nabla u\vert) + \Psi_2(\vert\nabla v\vert)] dxdy$$

dans le cas d'une régularisation au premier ordre. Les fonctions $\Psi_1$et $\Psi_2$utilisées, sont données par l'estimateur de Leclerc:

$$\Psi_i(w) = 1 - e^{-w^2/\sigma^2_i}, i = 1,2,$$

où les paramètres $\sigma_1$et $\sigma_2$sont à définir.

Grands déplacements

Lorsque les déplacements sont trop grands, il est intéressant d'utiliser une version intégrée de l'équation de conservation:

$$\frac{dI}{dt} = 0,$$

qui devient en intégrant:

$$I_1(\mathbf {x},t_1) = I_2(\mathbf {x}+\mathbf {d}, t_2),$$

qui fait explicitement intervenir le déplacement $d$entre l'image au temps $t_1$et au temps $t_2$. De même lors de l'utilisation de l'équation de continuité:

$$\frac{dI}{dt} + I{\mathrm {div\,}}\mathbf {w}= 0,$$

qui devient:

$$I_1(\mathbf {x},t_1) = I_2(\mathbf {x}+\mathbf {d}, t_2) \exp({\mathrm {div\,}}\mathbf {d}).$$

Le fait d'utiliser une telle équation intégrée, apporte de nouvelles instabilités puisque le problème devient fortement non-linéaire. Ce problème est géré en utilisant une implémentation multiéchelle et multigrille.

Terme centré pour le calcul des fonctions ${\mathrm {div\,}}$et ${\mathrm {curl\,}}$

Lors du calcul des fonctions ${\mathrm {div\,}}$et ${\mathrm {curl\,}}$, le masque généralement utilisé ne prend pas en compte le point central ce qui peut générer des instabilité lorsque le rapport signal sur bruit est faible. L'utilisation d'un masque centré, permet en général d'être plus robuste, on remplace:

$$\left\{ \begin{array}{rcl} {\mathrm {div\,}}\mathbf {d}(i,j)... ...1} - v_{i,j-1}) -(u_{i+1,j} - u_{i-1,j})] \end{array} \right.$$

par:

$$\left\{ \begin{array}{rcl} {\mathrm {div\,}}\mathbf {d}(i,j)... ...} + u_{i-2,j} - 6u_{i-1,j} + 2u_{i+1,j} ]. \end{array}\right.$$

Régularisation ${\mathrm {div\,}}-{\mathrm {curl\,}}$d'ordre deux

Lorsqu'on veut utiliser une régularisation ${\mathrm {div\,}}-{\mathrm {curl\,}}$d'ordre deux, on s'aperçoit que l'écriture des équations d'Euler-Lagrange associées à la fonctionnelle:

$$E_{reg} = \int_\Omega (\vert\nabla {\mathrm {div\,}}\mathbf ... ...rt^2 + \vert\nabla {\mathrm {curl\,}} \mathbf {w}\vert^2) dxdy$$

font apparaître des dérivées d'ordre quatre, difficiles à estimer de manière robuste. On peut, introduire des fonctions scalaires intermédiaires $\zeta$et $\xi$correspondant à des estimations des fonctions de divergence et de rotationnel de la vitesse $\mathbf {w}$à estimer:

$$\begin{array}{rcl} E_{reg}(\mathbf {w}) & = & \displaystyle ... ...ert^2 + \lambda\vert\nabla \zeta\vert^2\right)dxdy. \end{array}$$

où $\lambda$étant un paramètre positif. Cette fonctionnelle peut se minimiser de manière alternée par rapport aux variables $\xi$, $\zeta$et $\mathbf {w}$. Dans un premier temps, $\xi$et $\zeta$sont fixés et le champ $\mathbf {w}$est estimé. Une fois ce champ estimé, il est lui même fixé tandis que $\xi$et $\zeta$sont estimées successivement. Cette opération est itérée jusqu'à convergence.

Décomposition du champ de vitesse

On peut écrire que le champ de vitesse se décompose en trois:

$$\mathbf {w} = \mathbf {w}_{lam} + \mathbf {w}_{irr} + \mathbf {w}_{sol}.$$

La composante laminaire $\mathbf {w}_{lam} = \mathbf {w}_{glb}$correspond au transport global dans l'image. Tandis que les composantes irrotationelle $\mathbf {w}_{irr}$et solénoïdales $\mathbf {w}_{sol}$décrivent le mouvement local $\mathbf {w}_{irr} + \mathbf {w}_{sol} = \mathbf {w}_{loc}$. Ces composantes sont telles que:

$$\left\{ \begin{array}{ll} {\mathrm {div\,}}(\mathbf {w}_{irr... ...; & {\mathrm {div\,}}(\mathbf {w}_{sol})=0. \end{array}\right.$$

soit:

$$\left\{ \begin{array}{ll} \nabla.\mathbf {w}_{irr}=\nabla.\m... ....\mathbf {w}; & \nabla.\mathbf {w}_{sol}=0, \end{array}\right.$$

où $\nabla^\bot = (-\partial/\partial y, \partial / \partial x)$. Ce système est résolu dans le domaine de Fourier puis la transformée de Fourier inverse permet d'obtenir $\mathbf {w}_{irr}$et $\mathbf {w}_{sol}$. Thomas Corpetti doit m'envoyer le code Matlab qui effectue cette décomposition.

Notons que ces deux composantes peuvent également être exprimées sous la forme de deux fonctions potentielles $\phi$et $\psi$, telles que:

$$\left\{ \begin{array}{l} \Delta \phi = {\mathrm {div\,}}{\ma... ...lta \psi = {\mathrm {curl\,}}{\mathbf {w}}. \end{array}\right.$$

Une fois estimés, ces potentiels permettent de déterminer les composantes $\mathbf {w}_{irr}$et $\mathbf {w}_{sol}$par:

$$\begin{array}{l} \mathbf {w}_{irr} = \nabla \phi\\ \mathbf {w}_{sol} = \nabla^\bot \psi, \end{array}$$

Il est possible d'obtenir $\phi$et $\psi$directement à partir de la séquence d'image, mais ce travail est en cours pour le moment.

Travail réalisé

Les données

Nous disposons de trois types de données image:

• données de température de surface SST,
• données de couleur de l'océan CHL (chlorophylle),
• données de hauteur de surface GPS (altimétrie).

Ces images étant des simulations générées par le modèle de circulation océanographique du LODYC: OPA, nous disposons également des vitesses instantanées théoriques calculées par OPA:

• composante du vecteur vitesse $u_{{modèle}}$,
• composante du vecteur vitesse $v_{{modèle}}$.

Bilan des méthodes testées

Le tableau ci-dessous reprend les différentes méthodes d'estimation de mouvement testées:

conservation$\rightarrow$ régularisation $\downarrow$	luminance	volume
$L_2^*$	SST/CHL	SST/CHL
robuste	SST/CHL	SST/CHL
${\mathrm {div\,}}/ {\mathrm {curl\,}}$	SST/CHL	SST/CHL

$^*$: La régularisation notée $L_2$ici a été ``simulée'' en utilisant l'estimateur robuste $\Psi_2$mais en forçant le paramètre $\sigma$, plus ce paramètre est fort, plus on se rapproche d'une fonction quadratique, il a été mis à 100.0 pour les tests.

La composante laminaire a également été calculée pour les deux types d'image: SST et CHL. En utilisant l'hypothèse de conservation de la luminance et une régularisation $L_2$très fortement pondérée afin d'obtenir le mouvement global à divergence et rotationnel nuls. Pour les tests, le coefficient de pondération $\alpha$a été placé à 10000.

Dépouillement des résultats

Les trajectoires calculables par intégration des champs de vitesses n'ont pas été calculées sur place. Pour la comparaison entre les vecteurs $\mathbf {w}_{{img}}$et $\mathbf {u}_{{modèle}}$, on propose d'utiliser le critère de Barron. Ce critère mesure la moyenne et l'écart type d'une déviation angulaire 3D entre le mouvement réel et $\mathbf {u}_{{modèle}}$et le mouvement estimé $\mathbf {w}_{{img}}$:

$$\alpha_s = \arccos\left( \frac{\mathbf {w}_{{img}}}{\vert\ve... ...dèle}}}{\vert\vert\mathbf {u}_{{modèle}}\vert\vert} \right).$$

Vers un modèle dédié à l'océanographie

Xavier Vigan a développé lors de sa thèse au LODYC (manuscrit commandé en prêt au près de l'université de Paris 6, reçu en microfiches aujourd'hui), une méthode d'estimation des champs de vitesse à partir d'images de température de surface et d'une équation basée sur l'équation de continuité dite d'advection-diffusion pour la chaleur qui dans le cas général s'écrit pour un fluide quelconque:

$$\frac{d(\rho T)}{dt} + \rho T \, {\mathrm {div\,}}\mathbf {u... ...o}{\partial t}\frac{dp}{dt} = {\mathrm {div\,}}(K_T \nabla T),$$

où $\mathbf {u}$est le vecteur vitesse instantanée tridimensionnel. Dans le cas d'un fluide incompressible et sous les hypothèses d'approximation de Boussinesq on a:

$$\rho \frac{dT}{dt} = K_T \Delta T.$$

où $K_T$est le coefficient de diffusion de la chaleur supposé constant. Si on considère une couche comprise entre la surface et une profondeur $h$, on peut intégrer verticalement sur cette profondeur, on obtient quelque chose du type:

$$\frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} = K_T \Delta T + Src.$$

où $Src$est un terme de source qui dépend de la position (x,y) et où $\Delta T$représente ici le laplacien horizontal. D'après Vigan la diffusion horizontale est négligeable aux échelles de temps qui nous intéresse, donc $\Delta T$est négligeable. Il obtient donc:

$$\frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} = Src.$$

D'après Etienne Mémin, on devrait pouvoir considérer que le terme $Src$représente un ``écart'' par rapport au modèle de données considéré, et le gérer ainsi par l'utilisation d'un estimateur robuste $\Psi_1$adapté. D'autre part, n'oublions pas que le teme ${\mathrm {div\,}}\mathbf {u}=0$valable pour les fluides incompressibles n'est valable que lorsque l'on considère le mouvement fluide dans ses trois dimensions. Or on ne peut observer, dans le domaine image, que les composantes horizontales de ce vecteur, le terme en $I{\mathrm {div\,}}\mathbf {w}$devient donc intéressant pour expliquer les mouvements convectifs. Ces termes en $I{\mathrm {div\,}}\mathbf {w}$et $Src$sont donc concurrents.