Xu Chen

Pascal Zillé
email: pascal.zille@gmail.com

Liang Shao
email: liang.shao@ec-lyon.fr

Thomas Corpetti
email: thomas.corpetti@univ-rennes2.fr

Optical flow for incompressible turbulence motion estimation

Keywords: Optical flow, Large Eddy Simulation, Turbulence, Mixed model, incompressible flows

We propose in this paper a new formulation of optical flow dedicated to 2D incompressible turbulent flows. It consists in minimizing an objective function constituted by an observation term and a regularization one. The observation term is based on the transport equation of the passive scalar field. For non-fully resolved scalar images, we propose to use the Mixed model in Large Eddy Simulation to determine the interaction between large scales and unresolved ones. The regularization term is based on the continuity equation of 2D incompressible flows. Compared to prototypical method, this regularizer preserves more vortex structures by eliminating constraints over the vorticity field. The evaluation of the proposed formulation is done over synthetic and experimental images, the improvements in term of estimation are discussed.

Introduction

For decades, the study of fluid flow from image sequence is of great interest to fluid mechanics and computer vision community. Retrieving an accurate velocity field from image pairs is one of the main task in various scientific domains. For environmental science such as meteorology and oceanography, the satellite images provide information for the estimation of cloud motions and ocean current. In the field of experimental fluid dynamics, the imaging and measurement techniques constitute an unique basis for the study of unsteady flow. For example, Particle Image Velocimetry (PIV) technique [START_REF] Adrian | Particle-imaging techniques for experimental fluid mechanics[END_REF]) is a well established method, which consists in computing the instantaneous velocity field of flows based on the local spatial correlation.

Despite of the efficiency of PIV, it is less applicable when the seeding density is low and inhomogeneous, since the correlation peaks are hard to detect. In addition, the assumption that the variation of velocity field within the interrogation windows is negligible may lead to no reliable estimation in some particular cases. Furthermore, PIV technique is not able to cope with scalar images, though this latter constitute the main source of data in many environmental domains such as meteorology (satellite images).

Compared to correlation techniques used in PIV, optical flow performs more efficiently in dense estimation of velocity field from image sequence. This family of methods, originally proposed in computer vision community with the aim of extracting the apparent rigid motions from image pair [START_REF] Horn | Determining optical flow[END_REF], has become during the last decade more and more popular in fluid mechanic community for its ability of dense velocity estimation with continuous patterns. Besides, in comparison with PIV technique, optical flow is more adaptive to various physical constraint [START_REF] Barron | Performance of optical flow techniques[END_REF]Fleet and Weiss 2006).

The usual assumption behind optical flow is the wellknown brightness conservation [START_REF] Horn | Determining optical flow[END_REF], namely the classical Optical Flow Constraint Equation (OFCE):

∂ I ∂t (x,t) + u(x,t) • ∇I(x,t) = 0 (1)
where ∇ is the spatial gradient. This term assumes that the gray-level I of the particular point x = (x, y) T in two successive images is constant. Here, u represents the velocity field. However, OFCE can, by no means, resolve the two components of velocity field in iso-brightness regions as in homogeneous areas, since all terms in (1) vanish and the problem is ill-posed. This is known as the "aperture" problem in the computer vision community. To cope with this, a regularization term has to be introduced as a constraint over the spatial distribution of the velocity field. [START_REF] Horn | Determining optical flow[END_REF] have suggested to embed relation (1) in an optimization problem, which consists in minimizing an energy defined by an observation term based on OFCE and a regularization one represented by the magnitude of velocity gradient (the idea behind this assumption being that all neighbor points have similar motions):

J (u(x,t), I(x,t)) = J O (u(x,t), I(x,t)) + αJ R (u(x,t))

= Ω ∂ I ∂t (x,t) + u(x,t) • ∇I(x,t) 2 dx + α Ω ∇u(x,t) 2 dx (2)
where Ω denotes the image domain and α is a weighting coefficient. We refer to this equation as H&S for abbreviation in this paper. The minimization of this energy allows an unique solution since the function is convex.

Based on H&S, considerable work has been done for the application of optical flow in a context of fluid mechanics. [START_REF] Ruhnau | Variational optical flow estimation for particle image velocimetry[END_REF] proposed a continuous variational formulation, along with a multi-scale strategy for particle images; [START_REF] Liu | Fluid flow and optical flow[END_REF] used the projected-motion equations as the optical flow constraint; [START_REF] Héas | Dynamic consistent correlation-variational approach for robust optical flow estimation[END_REF] reinforced optical flow method by adding the robustness of correlation techniques; [START_REF] Papadakis | Variational assimilation of fluid motion from image sequence[END_REF], based on the optimal control theory, proposed a variational technique by imposing a temporal consistency. A detailed overview over different estimation techniques for fluid motion measurement is presented in [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF].

Despite the fact that great efforts have been done to improve fluid flow estimation, most of the approaches applied a spatial regularization term by assuming the smoothness of the velocity field over a small zone [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF][START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF][START_REF] Kadri-Harouna | Divergence-free wavelets and high order regularization[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. However, as [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF] demonstrated, the minimization of this term is equivalent to the penalization of the magnitude of divergence and vorticity field, which is critical in a context of fluid mechanics. They proposed a new formulation based on continuity equation with a secondorder div-curl regularizer able to keep the divergence and the vorticity of the flow for turbulent flows. For incompressible flows, the former should be zero, but the latter is related to the dissipation, which is not of the same magnitude with the first one. Interestingly, [START_REF] Guichard | Accurate estimation of discontinuous optical flow by minimizing divergence related functionals[END_REF] proposed to minimize a L 1 norm of divergence instead of the gradient but one should notice that the aim of this strategy is for discontinuity detection of 2D rigid motions, which is not based on the principles of fluid dynamics. In this work, we propose to develop a physics-based regularizer for 2D incompressible turbulent flows.

In addition, although numerous researches have been done over the variational optical flow formulations for flow estimation, the problem of turbulence has not been given enough attention. The turbulent flow is known as highly chaotic, and different physical processes may occur on various scales [START_REF] Pope | Turbulent flows[END_REF]. In optical flow, the image sequence is limited by the spatial resolution. More precisely, the concentration field characterized by image intensity is represented till the scale of pixel. The missing information under this scale may be negligible for rigid motion. However, when it comes to turbulent flow, the interaction between subpixel scales (unresolved scales) and resolved ones plays an important role in energy transport, which should be taken into account. Recently, [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF] and [START_REF] Zille | Observation models based on scale interactions for optical flow estimation[END_REF] introduced the concept of Large Eddy Simulation (LES), by considering the motions under pixel as subgrid turbulent dissipation. However, the coefficient of turbulent dissipation is computed a priori from velocity spectrum, which is not practically accessible since the velocity field is unknown neither in physical space nor in spectral one. To cope with this issue, we suggest in this study to apply physicalspace subgrid models and to compute systematically the subgrid scale term embedded in OFCE.

The objective of this paper thus is to propose a new formulation for 2D turbulent flows estimation. It is based on the prototypical method proposed by [START_REF] Horn | Determining optical flow[END_REF], which consists of the minimization of an objective function where we modified the observation and regularization terms: the observation models the SubGrid Scale (SGS) fluxes so as to predict them as accurately as possible; the regularization one eliminates constraints on the vorticity and is dedicated to incompressible flows. The evaluation of this formulation is done over synthetic and experimental images, and improvement of estimation is shown. The rest of the paper is organized as follows. In section 2, we present our physicsbased improved regularizer for 2D incompressible turbulent flows. Then in section 3, we propose a novel observation term based on the filtered transport equation with a physicalspace subgrid model. We describe the optimization strategy based on multiresolution and discretization scheme in section 4. The experimental evaluation of our method is described in section 5 and we compare the result with other ap-proaches in the literature. Finally, a conclusion is proposed in section 6.

A regularization term dedicated to turbulence

To deal with the aperture problem, the seminal work of [START_REF] Horn | Determining optical flow[END_REF] proposed to introduce an additional constraint based on the assumption on the smoothness of velocity field. The smoothness constraint, which is called the regularization term in optical flow community, is expressed as the module of the velocity gradient:

Ω ∇u(x,t) 2 dx (3)
This constraint is particularly suited for stable object movement or rigid motion. However, it is less appropriate when it comes to the fluid motion estimation. It can be indeed demonstrated that the penalization of the magnitude of the gradient of the velocity is equivalent to the penalization of the magnitude of divergence and vorticity [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF]:

Ω ∇ • u(x,t) 2 + ∇ × u(x,t) 2 dx (4)
Therefore such a penalization term does not seem adapted because the vorticity, directly related to the dissipation rate of fluid, plays an important role in the energy cascade. Although other higher-order regularizers have been proposed, most of them minimize more or less the magnitude of vorticity or its gradient [START_REF] Su | Scalar imaging velocimetry measurements of the velocity gradient tensor field in turbulent flows[END_REF][START_REF] Corpetti | Dense estimation of fluid flows. Pattern Analysis and Machine Intelligence[END_REF].

The key information contained in the vorticity field is then erased and the underestimation of vortices may become a source of error for velocity estimation. A new regularizer dedicated to flow estimation is then required in order to estimate accurately turbulent flows. For 2D incompressible flows, the material density in the fluid remains constant and the divergence of velocity is zero:

∇ • u(x,t) = 0 (5)
In addition, if we compare the OFCE in (1) to the transport equation of a passive scalar θ :

∂ θ ∂t (x,t) + ∇ • (u(x,t)θ (x,t)) = κ∆ θ (x,t) (6)
(∆ being the Laplace operator), we notice that two equations are equivalent when the divergence free constraint ( 5) is satisfied and the molecular diffusivity is neglected (κ ≪ 1). This means that the classical OFCE has already implicitly taken into account the divergence free condition. Therefore, we propose to use a simpler regularization term by penalizing the divergence of L 2 norm: Similar to the classical first-order regularization term, the proposed one constrains the global distribution of velocity field as well, since it is a gradient-based regularizer that enables estimation in regions where the spatial brightness gradient is weak. Besides, the new regularizer is much more physically meaningful since the key information in vorticity field is preserved.

J ′ R (u(x,t)) = Ω ∇ • u(x,t) 2 dx (7)

Observation term based on subgrid modeling

The classical observation term, namely OFCE, links the velocity field u(x,t) to the observed image sequence I(x,t).

For scalar images, the gray-level I is related to the passive scalar concentration field θ . [START_REF] Corpetti | Dense estimation of fluid flows. Pattern Analysis and Machine Intelligence[END_REF] proposed that I ∝ θ dz, where z denotes the observation depth. Fig. 1 shows the scalar variance spectrum along with the kinetic energy spectrum for small Schmidt number (Sc < 1), from which one can notice that there exists two inertial ranges for scalar spectrum: the inertial-convective range (k ≪ k θ ), in which the scalar fluctuation is driven by that of velocity; and the inertial-diffusive range (k θ ≪ k ≪ k η ), where the scalar fluctuation is mainly influenced by the molecular diffusivity. The smallest scales (k ≥ k η , diffusive range) are governed by the molecular viscosity, characterized by a exponential decay. In order to extract the total velocity field, the scalar concentration field must contain information till the Kolmogorov scale k η [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF].

In optical flow, however, due to the limitation of acquisition period and the spatial resolution of the image, it is rarely possible to extract the scalar concentration field till the diffusive range. Information represented by the gray-level I contains only the scalar concentration θ till the scale of pixel, which corresponds to the large scale motions (inertialconvective range) in most cases. The smaller ones (subpixel) are eliminated by a filter during image acquisition step. Nevertheless, for fluid flows, these small scale motions, contributing to the energy transfer in the wavenumber space, have an effect on the large scales. The interaction between Fig. 2 Decomposition of energy spectrum in LES [START_REF] Pope | Turbulent flows[END_REF]: the exact solution splits into the resolved scale and the subgrid scale (should be modeled) in spectral space. k c is the cutoff wave number associated with grid size in physical space.

the unresolved scales and resolved ones cannot be neglected, especially in the case of turbulent flows. Unfortunately the classical OFCE in (1) neglects such interactions. It is then of prime importance to rely on alternative models.

As a matter of fact, these sub-pixel motions correspond in essence to the unresolved scale motions in LES. This latter consists in calculating directly the low-frequency modes in space, while modeling the subgrid scale motions. A symbolic representation of the concept of LES is shown in Fig. 2. In optical flow, the cutoff wave number k c is related to the pixel scale ∆ of image sequence:

k c = π ∆ (8)
Applying a cutoff filter to the transport equation of θ and rewriting in dimensionless form, we obtain:

∂ θ ∂t + ∇ • uθ = 1 ReSc ∆ θ (9)
Here, θ and u denote respectively the scalar and velocity field in pixel grid. Re and Sc are Reynolds number and Schmidt number. We divide the filtered non-linear term uθ into the resolved part and residual part, then we have:

∂ θ ∂t + ∇ • uθ + ∇ • τ θ = 1 ReSc ∆ θ (10)
where τ θ is the subgrid scalar flux [START_REF] Sagaut | On the relationship between the mean flow and subgrid stresses in large eddy simulation of turbulent shear flows[END_REF], defined as:

τ θ = uθ -uθ (11)
Cassisa et al (2011) proposed that this subgrid term τ θ is related to the turbulent viscosity which can be expressed as:

τ θ = -D T ∇θ (12)
with D T a turbulent diffusion coefficient, which is a statistical constant calculated according to the formulation of Deardorff (1970):

D T = ∆ × ( 3 2C sgs ∞ k c E uu (K)dK) τ θ = c r u θ -u θ - ν sgs Pr sgs ∇θ (14)
with c r a global constant and Pr sgs subgrid Prandtl number, ranging from 0.1 to 1 [START_REF] Pope | Turbulent flows[END_REF]. In fact, the Mixed model ( 14) is a linear combination of a structural model and a viscosity model in LES.

The structural model is based on the hypothesis of scale similarity that the statistical structures of the flow at different scales are locally auto-similar. [START_REF] Bardina | Improved subgrid-scale models for large-eddy simulation[END_REF] assumed that the residual stress tensor is similar to that evaluated on the basis of the smallest resolved scales, since the subgrid motion shares a common local scaling law with the smallest resolved one. Later on, in [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible turbulent flows[END_REF], they proposed to apply the analytical filter twice for the simulation of subgrid tensor. It is known that the original Bardina model, i.e., the structural model without an eddy viscosity, is not purely dissipative, but the backscatter effect has been taken into account [START_REF] Sagaut | On the relationship between the mean flow and subgrid stresses in large eddy simulation of turbulent shear flows[END_REF]. In practice, to avoid the time integration instability, the Mixed model, i.e., the structural model with an eddy viscosity, is widely used in computational fluid community [START_REF] Sagaut | On the relationship between the mean flow and subgrid stresses in large eddy simulation of turbulent shear flows[END_REF].

The viscosity model is based on the hypothesis that the interaction between subgrid and resolved scales is essentially an energy cascade from large scales to smaller ones, and the effect of which is analogous to molecule diffusion and can be considered as a turbulent dissipation. By introducing the concept of turbulent subgrid eddy viscosity ν sgs , [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment*[END_REF] developed the first SGS model:

ν sgs = c s ∆ 2 2S i j S i j 1 2 (15)
where ∆ denotes grid size (pixel) and c s a universal constant between 0.1 and 0.2. S i j is called the filtered rate-of strain tensor, defined as:

S i j = 1 2 ∂ u i ∂ x j + ∂ u j ∂ x i (16)
Therefore, the Mixed model ( 14) is simply a combination of a subgrid viscosity model for the representation of the energy cascade mechanism, with a scale similarity model for the prediction of the subgrid tensor structure. It should be mentioned that the turbulence in the third direction is supposed to be zero since we focus on 2D incompressible turbulence within the scope of this paper.

Based on this modeling, we may end up with the definition of our novel observation term by representing all physical mechanism of full field:

J ′ O (u, θ ) = Ω ∂ θ ∂t + u • ∇θ + ∇ • τ θ - 1 ReSc ∆ θ 2 dx ( 17
)
with the assumption that θ is proportional to the image brightness. In practice, for particle or scalar images of which the concentration field is fully resolved (till the diffusive range), no subgrid modeling is needed. In these cases, the observation term will be written as:

J ′ O (u, θ ) = Ω ∂ θ ∂t + u • ∇θ - 1 ReSc ∆ θ 2 dx ( 18
)
In previous sections, we have defined our novel regularizer and the observation term for 2D incompressible fluid flow. The estimation of fluid velocity field is conducted by minimizing the sum of the two terms:

min u J ′ (u, θ ) = J ′ O (u, θ ) + β J ′ R (u) (19)
Here, β is a weighting coefficient and J ′ O (u, θ ) is expressed in (17) for non-fully-resolved scalar image, and in (18) for particle image and fully-resolved scalar image. The assumption that θ ∝ I has been made in this paper. Let us now turn to some practical issues related to images, optimization and discretization.

Discretization and Multiresolution technique

It should be mentioned that the input images have to be preprocessed before computation since the artifact and noise may lead to a misestimation of fluid velocity. A recursive Gaussian filter [START_REF] Deriche | Recursively implementating the gaussian and its derivatives Fleet D, Weiss Y[END_REF] is used to remove these noises. Moreover, in order to reserve more information in images, it is preferred to use the convolution property for spatial gradients calculation:

∂ ∂ x i (G * θ) = ∂ G ∂ x i * θ (20)
Thus, the spatial derivative is firstly computed on the filter and then convolved with the scalar concentration field.

Discretization scheme

The discretization process of the minimization of ( 19) could be conducted by resolving the associated Euler-Lagrange equations:

∂ J ′ ∂ u i - ∂ ∂ x j ∂ J ′ ∂ u i, j = 0 ( 21
)
where

u i, j = ∂ u i ∂ x j .
It should be noted that the SGS flux term τ θ in the observation term (17) could be considered practically homogeneous at pixel level, so that its dependency on local velocity u i and its gradient u i, j is weak. For the simplicity of minimization scheme, we suppose that:

∂ τ θ ∂ u i = 0 (22) ∂ τ θ ∂ u i, j = 0 (23)
Based on these assumptions, we obtain:

∂ θ ∂ x i ( ∂ θ ∂t + u j ∂ θ ∂ x j + ∂ τ θ j ∂ x j - 1 ReSc ∂ 2 θ ∂ x j ∂ x j ) = β ∂ ∂ x i u j, j (24)
We denote u by (u, v), then by neglecting higher order terms, ∂ ∂ x i u j, j can approximatively be discretized as ( 25) and ( 26), where h represents the pixel length and (m, n) is the spatial coordinate. ( 24) can be solved iteratively by using ( 27), in which we use the value of τ θ at iteration p to compute the velocity field at iteration p + 1.

∂ ∂

x ∂ u ∂ x + ∂ v ∂ y ≈ (u m+1,n + u m-1,n -2u m,n ) /h 2 + (v m+1,n+1 -v m-1,n+1 -v m+1,n-1 + v m-1,n-1 ) /(4h 2 ) (25) ∂ ∂ y ∂ u ∂ x + ∂ v ∂ y ≈ (v m,n+1 + v m,n-1 -2v m,n ) /h 2 + (u m+1,n+1 -u m-1,n+1 -u m+1,n-1 + u m-1,n-1 ) /(4h 2 ) (26) u p+1 v p+1 = β    ∂ θ ∂ x 2 + 2β h 2 ∂ θ ∂ x ∂ θ ∂ y ∂ θ ∂ x ∂ θ ∂ y ∂ θ ∂ y 2 + 2β h 2    -1   2 h 2 u p + ∂ ∂ x ∂ u p ∂ x + ∂ v p ∂ y -1 β ∂ θ ∂ x ∂ θ ∂t -1 ReSc ∆ θ + ∇ • τ p θ 2 h 2 v p + ∂ ∂ y ∂ u p ∂ x + ∂ v p ∂ y -1 β ∂ θ ∂ y ∂ θ ∂t -1 ReSc ∆ θ + ∇ • τ p θ   (27)

Multiresolution

Generally, optical flow has proven to be high-performance for estimation of apparent motions of small magnitude. However, for large displacements between two images, the discordance between the temporal derivation and spatial gradient in (1) may lead to a poor estimation. To get rid of this issue, multiresolution technique is introduced (Bergen et al 1992), which consists of a incremental coarse-to-fine strategy. The image pyramid structure is built by a series of low-pass filters, such that image resolution is reduced successively by subsampling process (Fig. 3(a)). More precisely, let k denote of the current pyramid level, thus k + 1 is the upper level with a coarser pixel grid. At Level k, we firstly warp the second image to the first one with the estimation of Level k + 1 by using the bicubic interpolation:

θ w k (x,t + 1) = θ k (x + u k+1 ,t + 1) (28)
Then, we compute the increment velocity field δ u k based on θ k (x,t) and θ w k (x,t + 1) at Level k with (24). The estimation at Level k is finally updated by:

u k = u k+1 + δ u k (29)
Fig. 3(b) illustrates briefly this coarse-to-fine strategy. One may refer to [START_REF] Burt | Smart sensing within a pyramid vision machine[END_REF][START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]) for more details. Let us now present our experiments.

Experimental validation

Synthetic images

Data description

The sequence is downloadable from the internet computation. The homogeneous and isotropic turbulent flow represented in the image sequence, with Reynolds number Re = 3000 and Schmidt number Sc = 0.7, has an average velocity around 3.5 pixels for each time interval. Both of scalar and particle sequence are composed by 100 successive images, as shown in Fig. 4(a) and in Fig. 4(b), respectively.

Parameters setting -Scalar images

To remove the noise and artifacts, we firstly filtered the images with a gaussian filter. At the same time, some small scale structures have also been filtered during this process, as shown in Fig. 5. From the point of view of LES, these small scales correspond to the filter width ∆ f ilter that is larger than the pixel size ∆ pixel . A subgrid model based on the filter width ∆ = ∆ f ilter (σ ) is therefore needed for the estimation of motion within these scales. To ensure the pertinence of the subgrid model, the variance σ is need to be chosen in such a way that the corresponding filter width ∆ lays in the inertial subrange. In practice, we use a Gaussian filter of variance σ = 0.9. ∆ can then be computed with ∆ = 2 √ 3σ ≈ 3.2 [START_REF] Pope | Turbulent flows[END_REF], and is set to be 4 approximatively. This choice leads to a small ∆

L θ ≈ 4 62.4 < 0.1 (L θ = π 2θ 2 ∞ 0 E θ (k)
k dk, is the integral length scale), indicating that the energy-containing large scales is not altered by the gaussian filtering process.

The multiresolution technique is not used in this case since the magnitude of displacements is small. As for the subgrid model, the global constant c r and Pr sgs are set to be 1.1 and 0.6 [START_REF] Bardina | Improved turbulence models based on large eddy simulation of homogeneous, incompressible turbulent flows[END_REF][START_REF] Sagaut | On the relationship between the mean flow and subgrid stresses in large eddy simulation of turbulent shear flows[END_REF], respectively. In addition the value c s = 0.19 is assumed and the weighting coefficient β in ( 19) is set to be 150 in this case. The velocity field estimated is filtered with a variance σ = 1 to remove the noise. Lastly the subgrid viscosity ν sgs is averaged within the whole image for a better estimation, since the flow is homogeneous and isotropic. In order to introduce a temporal consistency, at each time step, the velocity field is initialized by the previous one.

-Particle images

We use a 2-level pyramidal decomposition for coarse-to-fine estimation. PIV images are filtered by a recursive Gaussian filter of variance σ = 0.6 at each level and the velocity field estimated is equally filtered with a variance σ = 1.25 in order to remove the noise. The weighting coefficient β in ( 19) is set to be 4000. Compared to the scalar images, the value of β is greater in this case in order to compensate the differences in term of gray-level between these two sequences. Therefore in both cases the observation term and regularization term have the same order of magnitude.

Validation criteria

In order to evaluate statistically the performance of our proposed approach, we calculate the Average Angle Error (AAE) and Root Mean Square Error (RMSE) as evaluation criteria, which are commonly used in optical field community [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]:

AAE(t) = 1 N ∑ x∈Ω arccos u c (x,t)u e (x,t) u c (x,t) u e (x,t) (30)
RMSE(t) = 1 N ∑ x∈Ω u c (x,t) -u e (x,t) 2 (31)
Here, N represents the total number of pixels and u c (x,t) and u e (x,t) are respectively the exact and the estimated velocity field.

Results

-Scalar images

We present in Fig. 6 the vorticity map for a given instantaneous velocity field (at time t = 50) obtained by different methods. The associated velocity field is superimposed to each map and the results share a common color legend. We can clearly observe that the proposed formulation enables to improve the estimation. The vortex structures of different scales, which are comparatively blurred with H&S and are slightly detected by [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF], are much more accurately estimated with the subgrid modeling, such that the vortices have a similar behavior to that of DNS .

To go into more details, Fig. 7 presents the the spectra of velocity field and vorticity field in log-log and semi-log representation. Compared to the classical H&S (green), the formulation of [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF] (blue) and the proposed one (red) lead to a better recovery of the field of small scales, since both of the two approaches have introduced subgrid models that manage to predict small scale motion effects with resolved scales information. Furthermore, the proposed method exhibits a better performance for large scale motion (see semi-log representation in Fig. 7(b,d)). It is due to the new regularization term that maintains the rotational structure, and to the Mixed modeling that predicts the behavior of small scales (Shao et al 1999).

In terms of quantitative values, the evaluation criteria AAE and RSME are shown in Fig. 8. Results of the proposed formulation (solid line) are compared to H&S (dotted line) and Cassisa et al (2011) (dashed line). It is visible that the introduction of sub-pixel motion modeling [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF] has already improved the precision of estimation by 40% for AAE and 34% for RMSE. However, even compared to this high-performance method, our new formulation leads to a further improvement by 50% approximatively for both criteria. This significant improvement may be due to the performance of the Mixed modeling ( 14), that better predicts the interaction between resolved scale motions and those under pixel. Whereas the proposed regularizer that avoids penalizing vortex structures of the flow, helps additionally to improve the precision of the estimation. This first series of experiments in synthetic scalar images demonstrates the ability of the proposed method to properly recover turbulent velocity fields. Let us now turn to particle data.

-Particle images

In Fig. 9, we depicted our velocity vector maps over a zoom area (10 × 10 pixels) of PIV image (Fig. 4(b)), obtained by our methods (red), H&S (blue) and compared to the ground truth (black). It is obvious that the velocity field from the proposed formulation is more closer to the exact one than the classical optical flow method, which highlights the performance of the novel regularizer for 2D incompressible flow estimation. The first order regularization term of the H&S method prevents from a proper recovery of the vorticity which is clearly underestimated. This improvement can be illustrated statistically in Fig. 10, where AAE and RSME are calculated for particle sequence. Compared to H&S, the proposed formulation allows an improvement of 22% for AAE and 24% for RSME. This improvement is not as great as in scalar case, since the spatial gradient of the images, based on which the velocity is calculated, is better formulated for particle sequence. And the classical optical flow method, along with the multiresolution technique, has already allowed an estimation with a satisfactory precision. The novel regularization term, even in this case, manages to increase the precision by more than 20% by amending the spatial distribution of velocity field, which highlights the performance of the proposed regularizer for 2D incompressible flows.

-Comparison with other methods

In order to show the performance of the proposed formulation, we compare our results with other state-of-the-art methods in Fig. 11. We plot AAE and RSME errors for both particle (above) and scalar (below) sequences. Error values are taken from [START_REF] Héas | Dynamic consistent correlation-variational approach for robust optical flow estimation[END_REF].

For scalar images, the proposed estimation outperforms other state-of-the-art methods in optical flow community. This gap is mainly due to the introduction of subgrid model into the observation term, which takes into account the role of sub-pixel motions of turbulent flows.

As for particle images, it is shown that the proposed regularizer outperforms some of the high-accuracy algorithms in the literature (compared to [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] and [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF], the estimation improved by 43.5% and by 9.6% respectively). The proposed formulation is comparable to the regularizer based on the scale invariance property [START_REF] Héas | Dynamic consistent correlation-variational approach for robust optical flow estimation[END_REF], while higher-order regularizers [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF][START_REF] Kadri-Harouna | Divergence-free wavelets and high order regularization[END_REF] lead to estimations with best accuracy (with an improvement of 13% in comparison with our formulation). A further research dedicated to particle images is needed to find a better formulated observation term Fig. 12 A symbolic representation of Batchelor regime and regularizer. From these experiments, one concludes that our technique is vey well adapted to the estimation of turbulence motions. In the next section we show some experiments on real data.

Experimental scalar image

Data description

An experimental study of the dispersion of passive scalar in a 2D turbulence at high Péclet number is reported in [START_REF] Jullien | Experimental observation of batchelor dispersion of passive tracers[END_REF]. The Batchelor regime, characterized by k -1 spectrum, is observed. In this regime, only one velocity scale is active in dispersion process. [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid[END_REF] derived a diffusive cutoff wave number, namely the Batchelor wave number, for scalar field with a very-high Schmidt number (Sc ≫ 1):

k B = ε νκ 2 1/4 (32)
where ε, ν and κ represent respectively turbulent dissipation, kinematic viscosity and molecular diffusivity. Fig. 12 is a symbolic representation of the scalar variance spectrum with the kinetic energy spectrum. Compared to Fig. 1, a viscousconvective range [START_REF] Sagaut | On the relationship between the mean flow and subgrid stresses in large eddy simulation of turbulent shear flows[END_REF] is observed, where the scalar fluctuations are influenced mainly by the viscous effect of velocity field.

The scalar concentration field is visualized by using a 512 × 512 CCD camera with 0.2mm overall spatial resolution (Fig. 14(a) and Fig. 14(b)). The Péclet number (Pe = UL/κ = ReSc) is approximatively equal to 10 7 . One may refer to [START_REF] Paret | Intermittency in the two-dimensional inverse cascade of energy: Experimental observations[END_REF] for details of experiment setup.

Parameters setting

Fig. 13 shows the scalar variance spectrum of the experimental image. The k -1 law is fairly observed, and the spectrum drops rapidly beyond these scales, indicating an effective dissipation region [START_REF] Jullien | Experimental observation of batchelor dispersion of passive tracers[END_REF]. Since the wave number associated with pixel is within the diffusive range, the subgrid-scale model is unnecessary for flow estimation. But this case is still of interest to test the new regularization term. We use a 2-Level multiresolution pyramid and β is set to be 100. For the purpose of noise removing, we filter the image with σ = 4 at finer grid level and with σ = 2 at coarser grid level. Test with H&S has also been conducted for comparison. .14(c) and Fig.14(d) show respectively the vorticity field estimated based on H&S and the proposed formulation. It is visible that the novel formulation extracts more precisely the vortex structure of various scales. The gradient-based regularizer (3) constrains the velocity field to be oversmoothed in order to satisfy the continuity condition. The proposed regularization scheme, as expected, enables to estimate vortex structures more clearly and more distinguishably, while allowing a sufficient smoothness for the velocity field. This is to our opinion a very good behavior and once again, this proves the great efficiency of the proposed technique.

Results

Fig

Conclusion

In this paper, we have proposed a new formulation in the framework of optical flow for turbulence motion estimation. The formulation is based on the minimization of an objective function constituted by original observation and regularization terms: the observation is based on the transport equation of the passive scalar field, where we proposed to determine the value of SGS flux term τ θ by using the Mixed model for non-fully resolved scalar images; the regularization term is based on the continuity equation of 2D incompressible flows, which penalizes exclusively the divergence of the velocity field without any constraint on the vorticity field.

Tests on particle images and experimental images have proven the performance of the novel regularizer, results have shown that it allows a better prediction of the vortex structures compared to [START_REF] Horn | Determining optical flow[END_REF]. As for the synthetic scalar sequence, the Mixed model along with the novel regularizer, estimates precisely the behavior of subgrid scale motion and its effect to the resolved one (as shown in velocity spectra). The proposed formulation leads to increase the precision by 50% in comparison with [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF] and outperforms other state-of-the-art optical flow methods in the literature.

Fig. 1

 1 Fig.1A symbolic illustration of the kinetic energy spectrum (dashed line) and scalar variance spectrum (solide line) for Sc < 1[START_REF] Pope | Turbulent flows[END_REF]. k θ and k η represent respectively the wave number of Obukhov-Corrin and Kolmogorov.

Fig. 3 AFig. 4

 34 Fig. 3 A symbolic representation of the multiresolution pyramid (a) and the minimization process (b) with coarse-to-fine strategy (Bergen et al 1992).

Fig. 5

 5 Fig. 5 Scalar variance spectra of the original image and the filtered one with variance σ = 0.9.

Fig. 6

 6 Fig. 6 Instantaneous vorticity field on scalar sequence at t = 50 : exact DNS vorticity field (a), result of H&S (b), result of Cassisa et al (2011) (c), estimation with proposed formulation (d).

Fig. 7

 7 Fig. 7 Spectra of velocity and vorticity with log-log and semi-log representation: results of DNS (black), H&S (green), Cassisa et al (2011) (blue) and proposed formulation (red).

Fig. 8

 8 Fig. 8 AAE (a) and RMSE (b) errors of estimation for synthetic scalar images.

Fig. 9

 9 Fig.9The velocity vector map over a zoomed area in Fig.4(b): DNS result (black), estimation with proposed formulation (red) and H&S (blue).

  Fig. 10 AAE (a) and RMSE (b) errors of estimation for synthetic particle images.

  Fig. 11 AAE and RMSE errors for particle (above) and scalar (below) images of different methods. The data are taken form Héas et al (2012).

Fig. 13

 13 Fig. 13 Scalar variance spectrum of experimental image (t = 94).

Fig. 14

 14 Fig. 14 Experimental scalar images (a) and (b), with the estimated vorticity map from Horn and Schunck (1981) (c) and proposed formulation (d).

(13) Despite of the performance of this SGS modeling (Cassisa et al 2011), the main drawback lies in the fact that the power spectrum E uu cannot be obtained explicitly if the velocity field is unknown. Consequently, the value of D T remains indeterminable. Besides, the SGS viscous models is purely dissipative, the backscatter effect in energy cascade has not been taken into account[START_REF] Bertoglio | A stochastic subgrid model for sheared turbulence[END_REF]. To avoid this ambiguity, we propose to determine the value of SGS flux term τ θ by using the Mixed model proposed by[START_REF] Bardina | Improved subgrid-scale models for large-eddy simulation[END_REF]:

http://fluid.irisa.fr/data-eng.htm