Journ´ees Activit´es Universitaires de M´ecanique

La Rochelle, 31 aoˆut et 1er septembre 2006

Suivi de mouvement fluide : approches stochastique et variationelle ´ Nicolas Papadakis, Anne Cuzol et Etienne M´emin IRISA/Univ. Rennes I, Campus Beaulieu, 35042 Rennes Cedex, France [email protected] ;[email protected] ;[email protected]

R´esum´e : Dans ce papier, nous pr´esentons deux m´ethodes diff´erentes permettant de reconstruire le champ des vitesses apparentes sur le plan image d’un e´ coulement fluide a` partir de mesures bruit´ees extraites de la s´equence d’images. La premi`ere m´ethode s’appuie sur le filtrage stochastique tandis que la seconde fait appel a` l’assimilation de donn´ees.

Abstract : In this paper, we present two different methods allowing to reconstruct the apparent motion field of a fluid flow on the image plan. This reconstruction is obtained from noisy measurements extracted from the image sequence. The first method relies on the stochastic filtering, whereas the second one uses data assimilation framework.

Mots-clefs :

1

mesure du mouvement, filtrage stochastique, assimilation de donn´ees

Introduction

L’analyse de s´equences d’images d´ecrivant des ph´enom`enes fluides est de grande importance dans de nombreux domaines. Dans le domaine des sciences g´eophysiques par exemple, on pourra citer la m´et´eorologie et l’oc´eanographie, o`u l’on cherche a` suivre des syst`emes nuageux, a` estimer des courants oc´eaniques ou e´ valuer la d´erive d’entit´es passives telles que des icebergs ou des agents polluants. L’analyse d’images d’´ecoulements fluides est e´ galement cruciale en m´ecanique des fluides exp´erimentale pour un certain nombre d’applications. Pour tous les domaines mentionn´es ci-dessus il est d’un int´erˆet majeur de suivre au cours du temps des structures repr´esentatives de l’´ecoulement, et ceci le plus pr´ecis´ement possible. Lorsque les images sont bruit´ees ou lorsque les d´eplacements sont de grande amplitude et chaotiques (pour des e´ coulements turbulents par exemple), l’estimation de mouvement est difficile et sujette aux erreurs. La plupart des approches n’utilisent en effet qu’un ensemble r´eduit d’images (classiquement deux images cons´ecutives d’une s´equence) et peuvent donc souffrir d’une incoh´erence temporelle d’une image a` l’autre. L’extension des r´egularisations spatiales aux r´egularisations spatio-temporelles [1] ou l’introduction d’une contrainte dynamique simple dans les techniques de segmentation de mouvement reposent sur des hypoth`eses grossi`eres adapt´ees uniquement au mouvement d’objets rigides [2]. De mani`ere a` am´eliorer le suivi de structures dynamiques a` partir de s´equences d’images d’´ecoulements fluides, nous proposons deux m´ethodes de suivi conduisant a` une estimation robuste et coh´erente dans le temps des champs de d´eplacements instantan´es tout au long de la s´equence.

2

Notations

Tout champ de vecteurs bidimensionnel w d´efini sur un ensemble born´e Ω de R 2 et C 2 (Ω), s’annulant a` l’infini, peut eˆ tre d´ecompos´e en somme d’une composante irrotationnelle a` vorticit´e nulle et d’une composante sol´eno¨ıdale a` divergence nulle. Cette d´ecomposition est connue sous le nom de d´ecomposition de Helmholtz. Quand on ne dispose pas de conditions nulles au bord, une nouvelle composante de transport w∞ a` divergence et rotationnel nuls doit eˆ tre ajout´ee a` la d´ecomposition. La d´ecomposition s’´ecrit alors : w = w irr + wsol + w∞. En remplac¸ant les deux composantes w irr et wsol par leurs expressions en termes de fonctions de potentiel et en calculant la divergence et la vorticit´e des champs nous pouvons e´ crire les fonctions de potentiel comme solutions de deux e´ quations de Poisson : ∆φ = div w

et

∆ψ = −curl w,

(1)

avec ∆ l’op´erateur Laplacien. Les champs de vecteurs w irr et wsol s’´ecrivant respectivement comme le gradient et le gradient orthogonal des fonctions de potentiel φ et ψ, on retrouve : 1

Journ´ees Activit´es Universitaires de M´ecanique

wirr = K ∗ div w

La Rochelle, 31 aoˆut et 1er septembre 2006

et

wsol = −K ⊥ ∗ curl w,

(2)

o`u ∗ d´esigne le produit de convolution et K le gradient du noyau de Green associ´e au Laplacien : x K(x) = . (3) 2π|x|2 La deuxi`eme e´ quation de (2) est connue sous le nom d’int e´ grale de Biot-Savart. Les deux e´ quations indiquent que les composantes sol´eno¨ıdale et irrotationnelle du champ (et donc le champ complet) peuvent eˆ tre reconstruites si la divergence et la vorticit´e du champ de d´eplacements sont connues.

3

Suivi stochastique

Dans cette partie, nous d´ecrivons un algorithme de filtrage stochastique non lin´eaire adapt´e au suivi de mouvement fluide et bas´e sur une repr´esentation de faible dimension du champ de d´eplacements.

3.1 Particules de vortex et de source L’id´ee des m´ethodes de particules de vortex [3, 4] consiste a` approcher la vorticit´e d’un champ w par une somme discr`ete de Diracs r´egularis´ees et situ´ees en des vortex ponctuels x i : X curl w(x) ≈ γi f ∗ δ(x − xi ), (4) i

o`u le coefficient r´eel γi repr´esente la force associ´ee a` la particule i, f est une fonction de lissage radiale de domaine d’influence  et δ est la fonction de Dirac. Cette discr´etisation de la vorticit´e en un nombre limit´e d’´el´ements permet d’´evaluer le champ de d´eplacements directement a` partir de l’int´egrale de BiotSavart . Une repr´esentation analogue peut eˆ tre e´ crite pour la divergence du champ, a` partir de particules de source.

3.2 Estimation a` partir d’une s´equence d’images La discr´etisation de la vorticit´e a` l’aide de p particules de vortex conduit, grˆace a` l’int´egrale de BiotSavart, a` la repr´esentation suivante pour la composante sol´eno¨ıdale du champ de d´eplacement : p X wsol (x) ≈ (5) γisol K⊥sol (xsol i − x), i

i=0

o`u xi et i indiquent respectivement le centre et le domaine d’influence de K ⊥sol , le gradient orthogoi nal r´egularis´e du noyau de Green. Nous consid´erons ici que tous ces param`etres sont libres de varier d’une fonction a` l’autre. Une repr´esentation similaire peut eˆ tre obtenue pour la composante irrotationnelle a` l’aide de q particules de source. En conclusion, nous obtenons une approximation du champ de d´eplacements par une somme pond´er´ee de fonctions de base d´efinies par la position de leur centre et leur domaine d’influence. Cette repr´esentation param´etrique peut eˆ tre li´ee a` un mod`ele de variation spatio-temporelle de la luminance, de mani`ere a` estimer le mouvement fluide a` partir de deux images cons´ecutives de la s´equence d’images [5].

3.3 Probl`eme de filtrage Avant de d´ecrire en d´etails la technique de suivi que nous proposons, nous rappelons d’abord les principes du filtrage stochastique discret. Le filtrage stochastique a pour but de fournir des estimations sur une s´equence de variables al´eatoires d´ecrivant une cible d’int´erˆet a` diff´erents instants. La s´equence est suppos´ee eˆ tre une chaˆıne de Markov cach´ee x 0:n = {x0 , ..., xn } a` instants discrets, de distribution initiale ` un instant donn´e, l’inf´erence est effectu´ee sur la base des estimations p(x0 ) et transition p(xk |xk−1 ). A pr´ec´edentes et des observations bruit´ees et incompl`etes. Ces mesures sont not´ees z 1:n = {z1 , ..., zn } et sont suppos´ees ind´ependantes conditionnellement a` l’´etat, de distribution de probabilit´e p(z k |xk ). Les filtres stochastiques constituent des proc´edures permettant d’estimer la distribution de probabilit´e a posteriori p(xk |z1:k ) de l’´etat sachant toutes les mesures jusqu’`a l’instant k. L’inf´erence est r´ealis´ee en deux e´ tapes : 2

Journ´ees Activit´es Universitaires de M´ecanique

La Rochelle, 31 aoˆut et 1er septembre 2006

• Connaissant p(xk−1 |z1:k−1 ), l’´etape de pr´ediction utilise la loi de transition p(x k |xk−1 ) pour effectuer une premi`ere approximation de l’´etat suivant. • La vraisemblance p(zk |xk ) fournie par la nouvelle observation z k est utilis´ee pour calculer la loi a posteriori a` l’instant k. Dans le cas g´en´eral o`u le mod`ele d´ecrivant l’´evolution de l’´etat ainsi que le mod`ele de mesures sont non lin´eaires, les e´ quations de filtrage peuvent eˆ tre approch´ees. La m´ethode du filtrage particulaire [6] est une m´ethode de Monte-carlo s´equentielle qui permet d’estimer cette loi a posteriori, encore appel´ee loi de filtrage. Dans notre cas, l’ensemble z1:n = {z1 , ..., zn } des mesures est extrait de la s´equence d’images a` partir ` chaque e´ tape k, l’´equation de mesure est d’un mod`ele de variation temporelle de l’intensit´e lumineuse. A k construite a` partir du champ estim´e w (obtenu a` partir des positions pr´edites des particules de vortex) et de la paire d’images (Ik , Ik+1 ).

3.4 Application au suivi de vortex Dans cette partie nous montrons comment un tel sch´ema peut eˆ tre adapt´e pour suivre la composante sol´eno¨ıdale (ie. les particules de vortex) du champ de d´eplacements. Par souci de simplicit´e nous supposons ici que le champ de d´eplacements complet est a` divergence nulle et se r´eduit donc a` une composante sol´eno¨ıdale. Nous rappelons la formulation de l’´equation de Navier-Stokes en termes de vitesse et de vorticit´e : ∂ξ + (w.∇)ξ = ν4ξ, (6) ∂t ` chaque pas de o`u ξ repr´esente la vorticit´e, d´efinie par ξ = curl w et ν est le coefficient de viscosit´e. A temps, l’´evolution de la vorticit´e est d´ecrite par une e´ quation de convection-diffusion. Cette e´ quation (6) peut eˆ tre r´esolue par deux e´ tapes distinctes. Cette mani`ere de traiter l’´equation, –connue sous le nom de m´ethode de viscous splitting [3] – permet de traiter successivement la partie non visqueuse de convection et la partie visqueuse de diffusion de l’´equation. Le transport total de vorticit´e dˆu a` la convection et a` la diffusion est finalement donn´e par la combinaison de ces deux solutions. L’´etape de convection est r´ealis´ee par une int´egration avant. Pour des forces et domaines d’influence fix´es les positions successives des particules d´ecrivent le champ de d´eplacements sous-jacent a` partir de l’´equation (5) section (3.2). La partie diffusion peut eˆ tre simul´ee par une m´ethode de marche al´eatoire [3]. Cette technique repose sur l’interpr´etation probabiliste de la solution de l’´equation de la chaleur [7]. L’avantage de cette approche est qu’elle permet d’int´egrer naturellement un terme d’incertitude (brownien) sur l’erreur du mod`ele. La m´ethode consiste finalement a` d´eplacer les particules de vortex par leur propre vitesse (´etape de convection) et a` leur ajouter une perturbation gaussienne appropri´ee pour simuler la diffusion.

3.5 R´esultats Nous pr´esentons dans cette partie les r´esultats obtenus sur une s´equence r´eelle de l’ONERA de 80 images montrant l’´evolution de vortex g´en´er´es au bout d’une aile d’avion. Le suivi est initialis´e a` l’aide de la m´ethode d’estimation d´ecrite dans la partie (3.2). L’´etat initial est pr´esent´e figure 1. Sur cette figure sont repr´esent´es la premi`ere image de la s´equence, le champ de vecteurs initial et la carte de vorticit´e correspondante. Le r´esultat a e´ t´e obtenu avec un ensemble de 15 particules de vortex. La figure 2 montre le r´esultat obtenu par la m´ethode de filtrage non lin´eaire d´ecrite dans la section pr´ec´edente. On peut constater sur la ligne (a) que la solution obtenue par la m´ethode de suivi est toute a` fait satisfaisante. Le mouvement du vortex est bien reconstruit a` chaque instant. Ceci est e´ galement illustr´e par les cartes de vorticit´e de la ligne (b). La d´eformation du vortex est bien reconstruite et suit les contours photom´etriques.

3

Journ´ees Activit´es Universitaires de M´ecanique

(a)

La Rochelle, 31 aoˆut et 1er septembre 2006

(b)

(c)

F IG . 1 – (a) Premi`ere image de la s´equence ; (b) champ de d´eplacements estim´e au temps k = 0 avec la m´ethode pr´esent´ee en section 3.2 ; (c) distribution de vorticit´e correspondant a` l’ensemble des particules de vortex utilis´ees pour initialiser le suivi.

(a)

(b) k = 30

k = 50

k = 70

´ F IG . 2 – Evolution du champ de d´eplacements (a) et des cartes de vorticit´e (b) obtenues par filtrage particulaire pour diff´erents instants du suivi.

4

Suivi variationnel

Nous proposons a` pr´esent une technique de suivi variationnel s’appuyant sur des concepts de contrˆole optimal d´evelopp´es en assimilation de donn´ees [8, 9, 10]. L’approche propos´ee consiste a` formaliser l’estimation d’une fonction d’´etat, comme la minimisation d’une fonction de coˆut spatio-temporelle globale et continue d´eduite d’une estimation bay´esienne au sens du Maximum A Posteriori. La technique a e´ t´e appliqu´ee au suivi de champs de vecteurs, au travers de leur vorticit´e et divergence. Comme le d´emontrent quelques exp´erimentations, cette m´ethode permet de g´erer naturellement le bruit ainsi que l’occlusion partielle ou compl`ete des donn´ees de la s´equence d’images sur une longue dur´ee.

4.1 Formulation du suivi variationnel Dans ce paragraphe, nous d´ecrivons le syst`eme g´en´eral mis en place pour notre probl`eme de suivi. Nous supposons ainsi que notre fonction d’´etat, repr´esentant la vorticit´e et la divergence d’un champ de T vecteurs et not´ee X = [ξ, ζ] , suit un mod`ele d’´evolution temporel M, que la valeur initiale X 0 de cette fonction est connue, et que l’on poss`ede des observations de X not´ees y n a` certains temps discrets tn . On obtient alors le probl`eme mod`ele suivant :  ∂X  ∂t + M(X) = η(x, t) (7) X(x, t0 ) = X0 (x) + ν(x)  Y (x, t) = H(X) + ε(x, t) 4

Journ´ees Activit´es Universitaires de M´ecanique

La Rochelle, 31 aoˆut et 1er septembre 2006

η, ν et ε d´esignent des bruits gaussiens d´efinis sur le domaine image Ω, de moyennes nulles et de covariances respectives Q(x, t, x 0 , t0 ), B(x, x0 ) et R(x, t, x0 , t0 ). Ces bruits repr´esentent les diff´erentes erreurs associ´ees aux diff´erentes composantes de notre syst`eme, ils sont suppos´es non correll´es en espace et en temps. Le mod`ele d’´evolution temporel est donn´e par la formulation vorticit´e-vitesse de l’´equation de Navier-Stokes et une diffusion Brownienne de la divergence [11] : ∂ζ − νζ ∆ζ = ηζ . ∂t

∂ξ + w · ∇ξ − νξ ∆ξ = ηξ ∂t

(8)

Les observations discr`etes yn sont quant a` elles r´eunies dans une fonction d’observation continue Y . L’op´erateur H permet alors de passer de l’espace d’´etat dans celui des observations. On consid`ere ici que ces deux espaces sont les mˆemes, i.e H = Id. L’estimation de l’esp´erance conditionnelle E[X|X 0 , Y ] du syst`eme (8) conduit a` minimiser la fonctionnelle suivante : J(X) =

1 2 Z

Z Z Ω

T

(X(x, t0 ) − X 0 (x)) B −1 (x, x0 ) (X(x0 , t0 ) − X 0 (x0 )) dx0 dx Ω

 T  ∂X ∂X 1 −1 0 0 (x, t)Q (x, t, x , t ) + M(X) + M(X) (x0 , t0 )dt0 dx0 dtdx + 2 Ω,t Ω,t ∂t ∂t Z Z T 1 + (Y − H(X)) (x, t)R−1 (x, t, x0 , t0 ) (Y − H(X)) (x0 , t0 )dt0 dx0 dtdx. 2 Ω,t Ω,t Z



(9)

La solution obtenue correspond a` un lissage des observations y n suivant le mod`ele dynamique M.

4.2 Technique de r´esolution En calculant les e´ quations d’Euler-Lagrange et en introduisant une variable adjointe λ : λ=

Z Z Ω

τ

Q−1 0



 ∂X + M(X) dt0 dx0 , ∂t

(10)

ainsi que les op´erateurs tangent lin´eaires ∂X M et ∂X H d´efinis par dM(X + βθ) = ∂X M(θ), β→0 dβ lim

(11)

et ∂X M∗ et ∂X H∗ leur adjoints d´efinis au sens du produit scalaire < ∂ X MX 1 , X 2 >=< X 1 , ∂X M∗ X 2 >, on obtient une expression num´erique du gradient de la fonctionnelle, au moyen de l’algorithme it´eratif suivant [12] : λ(x, τ ) = 0 Z Z τ ∂λ − + (∂X M∗ λ = ∂X H∗ R−1 (Y − H(X))dtdx ∂t Ω 0 Z
B −1 (X(x0 , 0) − X 0 (x0 ) dx0 . λ(x, 0) =

(12) (13) (14)

Ω

La m´ethode consiste en une premi`ere int´egration ”avant” de la condition initiale X 0 a` l’aide de la loi dynamique M. La solution est ensuire corrig´ee par une int´egration r´etrograde de la variable adjointe selon le mod`ele adjoint (12, 13). L’´evolution de λ est guid´ee par l’´ecart entre les mesures et l’estimation : Y − H(X). La condition initiale est alors mise a` jour par l’´equation (14). Le processus est r´ep´et´e iterativement jusqu’`a convergence.

4.3 R´esultats Nous pr´esentons ici des r´esultats obtenus sur une s´equence m´et´eorologique montrant l’´evolution d’un cyclone sur la ligne (a) de la figure 3. La ligne (b) pr´esente les vorticit´es des champ de vecteurs utilis´ees comme mesures d’observation. Ces vecteurs de d´eplacements sont obtenus par une version rapide et grossi`ere d’une m´ethode d’estimation de champs dense [13]. On peut observer qu’il s’agit de mesures tr`es bruit´ees de la vorticit´e que l’on souhaite suivre. La ligne (c) illustre le r´esultat de l’assimilation de la vorticit´e. Le lissage a ainsi bien e´ t´e r´ealis´e dans la direction donn´ee par l’´equation dynamique, tout en respectant l’orientation des observations. 5

Journ´ees Activit´es Universitaires de M´ecanique

La Rochelle, 31 aoˆut et 1er septembre 2006

(a)

(b)

(c) k=1

k=7

k=14

F IG . 3 – (a) s´equence cyclone. (b) Vorticit´e observ´ee . (c) R´esultats de la vorticit´e assimil´ee.

5

Conclusion

Nous avons d´ecrit deux m´ethodes d´edi´ees au suivi de mouvement fluide. La premi`ere repose sur un filtrage stochastique non lin´eaire adapt´e, tandis que la seconde consiste a minimiser une fonction de coˆut spatio-temporelle. De telles m´ethodes permettent de reconstituer un ensemble de champs de d´eplacements a` partir de donn´ees extraites d’une s´equence d’images.

R´ef´erences [1] J. Weickert and C. Schn¨orr. Variational optic-flow computation with a spatio-temporal smoothness constraint. J. Mathematical Imaging and Vision, 14(3) :245–255, 2001. [2] G. Fanerback. Very high accuracy velocity estimation using orientation tensors parametric motion, and segmentation of the motion field. Proc. Int. Conf. Computer Vision, pages 5–26, 1999. [3] A. Chorin. Numerical study of slightly viscous flow. J. Fluid Mech., 57,, 57 :785–796, 1973. [4] A. Leonard. Vortex methods for flow simulation. J. Comp. Phys., 37, 1980. [5] A. Cuzol and E. M´emin. Vortex and source particles for fluid motion estimation. In Scale-Space’05. [6] A. Doucet, S. Godsill, and C. Andrieu. On sequential monte carlo sampling methods for bayesian filtering. Statistics and Computing, 10(3) :197–208, 2000. [7] B. Oksendal. Stochastic differential equations. Spinger-Verlag, 1998. [8] A.F. Bennet. Inverse Methods in Physical Oceanography. Cambridge University Press, 1992. [9] F.-X. Le Dimet and O. Talagrand. Variational algorithms for analysis and assimilation of meteorological observations : theoretical aspects. Tellus, pages 97–110, 1986. [10] O. Talagrand and P. Courtier. Variational assimilation of meteorological observations with the adjoint vorticity equation. I : Theory. Journ. of Roy. Meteo. soc., 113 :1311–1328, 1987. [11] A.J. Chorin and O.H. Hald. Stochastic Tools in mathematics and Science. Springer-Verlag, 2000. [12] N. Papadakis, E. M´emin, and F. Cao. A variational approach for object contour tracking. In ICCV’05 VLSM. [13] T. Corpetti, E. M´emin, and P. P´erez. Dense estimation of fluid flows. IEEE Trans. Pattern Anal. Machine Intell., 24(3) :365–380, 2002.

6