Vélocimétrie par images de particules

Jérôme Lemelin

La vidéo suivante a été produite lors du confinement dû à la pandémie de COVID-19 au printemps 2020 afin d’offrir un supplément de matière pour le cours de Méthodes d’analyse en océanographie physique (OCE77016-T3 et OCE72616-TA) qui a été forcé de se terminer précipitamment.

Pour réaliser cette vidéo, j’ai [enfin] tiré avantage de la fonctionnalité tactile de l’écran de mon ordinateur portable avec les logiciels open-source et gratuits Xournal et Vokoscreen, respectivement pour prendre des notes manuscrites à l’écran et pour la capture vidéo de l’écran.

Ici, je présente d’une part la capsule réalisée expliquant la méthode numérique qu’est la vélocimétrie par images de particules, et d’une autre part un complément théorique de la méthode qui suit le schéma explicatif, en y apportant ponctuellement des précisions. Le complément reprend également l’exemple de la capsule, soit celui d’un intrigant vortex dans le fjord du Saguenay.

Vélocimétrie numérique par images de particules

La vélocimétrie par images de particules, ou la particle image velocimetry (PIV) en anglais, est une méthode de vélocimétrie qui cherche à décrire l’écoulement d’un fluide en suivant un traceur passif prenant la forme d’un ensemble de particules sur une séquence d’instantanés ou d’images du fluide en question. La méthode se base sur le déplacement des particules d’une image à l’autre pour en tirer une mesure de vélocité.

La méthode peut être divisée en cinq grandes étapes présentées ci-dessous. À l’exception de l’acquisition des images, la méthode peut être effectuée entièrement numériquement et on parle alors de vélocimétrie numérique par images de particules, ou de digital particle image velocimetry (DPIV) en anglais.

Acquisition des images

Pour être en mesure d’effectuer une analyse par PIV de l’écoulement d’un fluide, il faut s’assurer que celui-ci contient un ensemble de particules discernables dont la trajectoire épouse fidèlement l’écoulement du fluide. Traditionnellement, il est possible d’«ensemencer» un fluide de particules pour s’en assurer. En océanographie physique, la méthode est plutôt appliquée sur un plan d’eau qui répond déjà à la contrainte de présence et de discernabilité des particules, pouvant se manifester par des débris ou des floes, par exemple.

Une particule idéale est une particule dont le nombre de Stokes $\mathrm{Stk}$ est très petit, à savoir

$\mathrm{Stk}=\dfrac{t_0\, u_0}{l_0}\ll 1$ ,

où $t_0$ , $u_0$ et $l_0$ sont respectivement le temps de relaxation de la particule, la vitesse caractéristiques de l’écoulement et la longueur caractéristique de la particule. Cette condition est généralement satisfaite lorsque la particule est petite et que sa masse volumique est sensiblement la même que celle du fluide.

En océanographie physique, les images sont usuellement celles de plans d’eau, prises à partir de la côte avec des caméras traditionnelles ou du dessus avec un drone ou un satellite, par exemple. Dans le présent exemple, les images ont été prises par drone et utilisent des floes comme particules.

Pré-traitement des images

Une séquence d’images d’un écoulement est rarement analysable pour la PIV immédiatement après l’acquisition. Il faut donc traiter les images acquises pour être en mesure d’en évaluer la vélocimétrie. Ces traitements peuvent être d’une part tout ce qui peut faire en sorte que les particules sont plus facilement discernables sur l’image, ou distinguable du fond (du fluide). D’une autre part, ces traitements peuvent aussi être une série de déformations sur l’image pour rendre la PIV applicable.

Numériquement, une image peut être représentée par une matrice (tableau) dont les éléments correspondent à une valeur d’intensité lumineuse. Les traitements peuvent donc être des ajustements de différents paramètres associés à la luminosité comme le contraste, par exemple. En terme de déformations, les traitements effectués cherchent plutôt à régulariser le domaine spatial, ou éviter des projections irrégulières de l’image.

Dans l’exemple, les images sont obliques et ont dû être d’abord géorectifiées. Ce traitement (g_rect) consiste à géoréférencer une image en tous points et la projeter sur un repère cartographique régulier. Les images ont ensuite été binarisées. Ce traitement (binarize) fixe un seuil d’intensité en-dessous duquel toute intensité est retournée vers le noir et au-dessus duquel toute intensité est retournée vers le blanc.

fig_preproc — Exemple de pré-traitement pour rendre l’analyse PIV exécutable.

Évaluation et analyse de vélocimétrie

L’analyse est le cœur de la méthode de la PIV (ou DPIV). La méthode se veut essentiel-lement statistique puisqu’elle doit suivre un ensemble de particules. Tel que mentionné précédemment, les images étudiées sont considérées comme des champs d’intensité lumineuse représentés numériquement par des matrices. Comme les particules sont suivies d’images en images, l’analyse est effectuée par paires d’images successives. Pour une paire d’images, la première va être représentée par la matrice $I$ prise au temps $t$ et la suivante par la matrice $I'$ prise au temps $t'=t+\Delta t$ telles que

$I(\mathbf{X}; \mathbf{x},t) \qquad \longrightarrow\quad \boxed{\quad\mathbf{v}(\mathbf{x},t)\quad}\quad \longrightarrow \qquad I'(\mathbf{X};\mathbf{x'},t'=t+\Delta t)$ ,

où $\mathbf{v}$ correspond à la vélocité des particules qui agit comme une déformation de $I\mapsto I'$ en déplaçant les particules de leurs positions initiales $\mathbf{x}$ à leurs positions finales $\mathbf{x'}$ sur le domaine $\mathbf{X}$ . Cette vélocité peut être exprimée telle que

$\begin{aligned}\int_{t}^{t'} \mathbf{v}\,\mathrm{d}t = \mathbf{D}(\mathbf{x};t,t') \iff \mathbf{v}=\dfrac{\mathbf{D}}{\Delta t} \quad \text{avec}\quad \mathbf{D}=\mathbf{x'}-\mathbf{x}\end{aligned}$ ,

de sorte qu’uniquement à partir du déplacement $\mathbf{D}$ des particules et du pas de temps $\Delta t$ entre les deux images, il est possible de déterminer la vélocité. Ce déplacement $\mathbf{D}$ des particules d’une image à l’autre peut réexprimer les matrices $I(\mathbf{x'}-\mathbf{D})$ et $I'(\mathbf{x'})$ de sorte à introduire une intuition de décalage entre les deux images.

La DPIV se définit alors un décalage $\mathbf{s}$ de positions entre les deux images tel que si $\mathbf{s}=\mathbf{D}$ , alors $\mathbf{x}+\mathbf{s}=\mathbf{x'}$ , pour comparer les images $I$ et $I'$ en sommant sur les particules par

$\begin{aligned}C(\mathbf{s})=(I\star I')(\mathbf{s})=\sum_{i} I(\mathbf{x'}_{i}-\mathbf{s})\,I'(\mathbf{x'}_{i})\end{aligned}$ ,

où $C(\mathbf{s})$ est la matrice de corrélation croisée entre les images telle qu’une est décalée de $\,\mathbf{s}\,$ sur l’autre. Cette opération calcule la corrélation entre les deux images en tous points pour tout décalage $\mathbf{s}$ .

piv_shift — Schématisation de la matrice de corrélation croisée. La valeur de chaque élément s de C est la corrélation entre I et I’ tel que I est décalé de s par rapport à I’.

La position du maximum de corrélation correspond alors au décalage d’une image sur l’autre qui maximise la corrélation entre les deux, donc correspondant le plus exactement au déplacement des particules d’une image à l’autre, soit

$\mathrm{max}\Big\{(I\star I)(\mathbf{s})\Big\} = (I\star I)(\mathbf{0}) \quad \text{et}\quad \mathrm{max}\Big\{(I\star I')(\mathbf{s})\Big\}= (I\star I')(\mathbf{D})$ .

Il est donc possible d’estimer le déplacement $\mathbf{D}$ des particules d’une image à l’autre en calculant la position par rapport au centre du maximum de corrélation croisée entre les deux images.

Pour une description plus globale et précise de l’écoulement, les images sont divisées en sous-figures par une fenêtre d’interrogation mobile qui évalue la PIV localement partout sur la paire d’images. La corrélation croisée est alors calculée pour chaque paire de sous-figures d’une image à l’autre et la vélocité peut être estimée localement partout.

piv_comp — Exemple de corrélations croisées pour des sous-figures générées par une fenêtre d’interrogation de dimensions de 128 px × 128 px. À gauche, une paire de sous-figures issues de la même image; le maximum de corrélation est au centre. À droite, une paire de sous-figures issues d’images successives; la position du maximum de corrélation par rapport au centre est le déplacement des particules.

La taille de la fenêtre d’interrogation doit être déterminée entre autres en termes de la densité de particules qui y est comprise et du déplacement maximal des particules qu’elle peut admettre, qui peut être mis en relation avec la fréquence d’échantillonnage. La taille doit être suffisamment grande pour conserver une corrélation entre deux images successives et doit être choisie de sorte à satisfaire localement le théorème d’échantillonnage de Nyquist. Ces conditions sont typiquement satisfaites par une taille de fenêtre d’interrogation qui isole en moyenne six particules.

Par sa définition et puisqu’on l’applique une multitude de fois, le calcul de la corrélation croisée sur une paire d’images par le biais d’une fenêtre d’interrogation mobile se veut une tâche numérique laborieuse. Ceci dit, par le théorème d’inversion (ou théorème intégral) de Fourier, il est possible de simplifier la tâche computationnelle pour obtenir le théorème de corrélation croisée.

Démonstration: En reprenant la définition de la corrélation croisée $C=A\star B$ avec $\mathbf{s}=(\imath,\jmath)$ telle que $A$ et $B$ forment une paire de sous-figures respectives de $I$ et $I'$ , on a

$\begin{aligned}(A\star B)(\imath,\jmath)=\sum_{m,n} A(m-\imath,n-\jmath)B(m,n)\end{aligned}$ .

Le théorème d’inversion de Fourier stipule que tout signal absolument intégrable et continu peut être réexprimé par la transformée inverse de sa transformée de Fourier. Ainsi, avec $\mathcal{A}$ et $\mathcal{B}$ respectivement les transformées de Fourier discrétisées de $A$ et $B$ , on a

$\begin{aligned}(A\star B)(\imath,\jmath)=\sum_{m,n}\left[\sum_{\mu,\nu} \mathcal{A}(\mu,\nu)e^{-2\pi i(m-\imath,n-\jmath)(\mu,\nu)}\sum_{\mu',\nu'}\mathcal{B}(\mu',\nu')e^{-2\pi i(m,n)(\mu',\nu')}\right]\end{aligned}$ .

Parce que $A$ est réelle, $\mathcal{A}$ peut être réexprimée en termes de sa conjuguée complexe $\mathcal{A}^{\ast}$ de sorte que

$\begin{aligned}\mathcal{A}(\mu,\nu)e^{-2\pi i(m-\imath,n-\jmath)(\mu,\nu)}=\mathcal{A}^{\ast}(\mu,\nu)e^{-2\pi i(\imath-m,\jmath-m)(\mu,\nu)}=\mathcal{A}^{\ast}(\mu,\nu)e^{-2\pi i(\imath,\jmath)(\mu,\nu)+2\pi i(m,n)(\mu,\nu)}\end{aligned}$ .

En distribuant le signe dans l’argument de l’exponentielle de la transformée inverse de $\mathcal{B}$ et en isolant les différentes dépendances des sommations par indépendance linéaire, on a

$\begin{aligned}~\hspace{6pt}(A\star B)(\imath,\jmath)&=\sum_{\mu,\nu} \sum_{\mu',\nu'}\mathcal{A}^{\ast}(\mu,\nu)\mathcal{B}(\mu',\nu')e^{-2\pi i(\imath,\jmath)(\mu,\nu)}\left(\sum_{m,n}e^{-2\pi i(m,n)(\mu'-\mu,\nu'-\nu)}\right)\\ ~ &=\sum_{\mu,\nu}\sum_{\mu',\nu'}\mathcal{A}^{\ast}(\mu,\nu)\mathcal{B}(\mu',\nu')e^{-2\pi i(\imath,\jmath)(\mu,\nu)} \delta_{\mu'\mu}\delta_{\nu'\nu}\end{aligned}$

où $\delta$ est un delta de Kronecker (équivalent à la fonction delta de Dirac dans le cas discret occurrent). On trouve finalement

$\begin{aligned}(A\star B)(\imath,\jmath) = \sum_{\mu,\nu}\mathcal{A}^{\ast}(\mu,\nu)\mathcal{B}(\mu,\nu) e^{-2\pi i (\imath,\jmath)(\mu,\nu)}\end{aligned}$ ,

qui correspond, finalement, à la transformée inverse du produit des transformées de Fourier de $A$ et $B$ , soit

$\begin{aligned}A\star B = \mathcal{F}^{-1}\big\{\mathcal{A}^{\ast}\cdot\mathcal{B}\big\}\end{aligned}$ .

Ce qu’il fallait démontrer. Ainsi, seulement à partir des transformées de Fourier des images, numériquement représentées par des Fast Fourier Transform (FFT), le calcul de corrélation au cœur de la méthode PIV peut être rendu plus efficace et simplifié.

Finalement, pour déterminer le déplacement des particules d’une sous-figure à l’autre, il ne reste qu’à trouver le maximum de corrélation croisée. Ce maximum se doit idéalement d’être confiné le mieux possible en un seul pic maximal de corrélation. Pour affiner ainsi le pic de corrélation croisée d’une paire de sous-figure, il est possible de calculer la corrélation croisée de sous-figures légèrement décalées dans le voisinage de la sous-figure initiale.

shifts_cut — Exemples de sous-figures légèrement décalées dans le voisinage diagonal de la sous-figure initiale, au centre.

La corrélation finale sera simplement le produit de toutes les corrélations croisées.

corr_comp — Exemple d’affinement de corrélation croisée. À gauche, la corrélation croisée d’une simple paire de sous-figures. À droite, le produit des corrélations croisées des paires de sous-figures décalées dans le voisinage des sous-figures initiales.

Pour trouver le pic, il ne suffit alors qu’à trouver la position du maximum de corrélation. Pour une mesure du déplacement pour exacte, il est possible d’interpoler le pic maximal de corrélation à l’échelle du sous-pixel en assumant un ajustement gaussien des points avoisinants le maximum. Un ajustement sur trois points peut être suffisant. Ayant calculé la corrélation croisée en tous points et ayant trouvé la position par rapport au centre à l’échelle du sous-pixel, la vélocité peut être calculée.

Validation et post-traitement

Ici, les vélocités calculées sont encore en px/s. L’analyse étant pratiquement exclusivement statistique, certaines mesures retournées de vélocité risquent d’être irréalistes physiquement. Seulement en se basant sur la distribution des composantes de vélocité, il est possible de discriminer une partie des valeurs irréalistes. Les valeurs conservées sont généralement contenues dans

$\begin{aligned}\mu_{u,v}-n\cdot\sigma_{u,v}\ \leq\ u,v\ \leq\ \mu_{u,v}+n\cdot\sigma_{u,v}\end{aligned}$ ,

avec $\mu$ et $\sigma$ respectivement la moyenne et la déviation standard des composantes $u$ ou $v$ de la vélocité et $n$ un nombre, généralement 3.

Mesures vélocimétriques

Ici, les vélocités calculées sont encore en px/s. La première étape − qui pourrait entrer dans le post-traitement − est de convertir les mesures de vélocité de px/s à m/s. Cette opération est plutôt aisément réalisée lorsque la résolution de l’échelle de l’image est connue, c.-à-d. la correspondance px : m.

L’analyse de la PIV (ou DPIV) est désormais terminée et différentes mesures dérivées du champ de vélocité peuvent alors être calculées, telle la grandeur $U$ de la vitesse, la vorticité $\omega$ , la divergence $\nabla\cdot\mathbf{v}$ , la cisaille $S$ , etc. telles que, par exemple

$\begin{aligned}U=\sqrt{u^2 + v^2}\ , \qquad \omega = \nabla\times \mathbf{v}\ ,\qquad S=\dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x}\ , \qquad \ldots\end{aligned}$