Patrick Héas
email: patrick.heas@inria.fr

Cédric Herzet

Etienne Mémin

Dominique Heitz

Pablo D Mininni

Bayesian estimation of turbulent motion

Keywords: optic flow, turbulence, robust estimation, constrained optimization, Bayesian model selection

Based on physical laws describing the multi-scale structure of turbulent flows, this article proposes a regularizer for fluid motion estimation from an image sequence. Regularization is achieved by imposing some scale invariance property between histograms of motion increments computed at different scales. By reformulating this problem from a Bayesian perspective, an algorithm is proposed to jointly estimate motion, regularization hyper-parameters, and to select the most likely physical prior among a set of models. Hyper-parameter and model inference is conducted by posterior maximization, obtained by marginalizing out non-Gaussian motion variables. The Bayesian estimator is assessed on several image sequences depicting synthetic and real turbulent fluid flows. Results obtained with the proposed approach exceed the state of the art results in fluid flow estimation.

I. INTRODUCTION

The inverse modeling of turbulent motion in images is an important issue in various application areas like meteorology, oceanography, or turbulence studies and experimental fluid mechanics. Motion characterization from experimental images is particularly important since a complete physical theory is still missing for turbulence phenomenology. Optic flow techniques estimating dense velocity fields are interesting in this context since they can describe precisely such complex and continuous flows. For these dense representations, regularization models are commonly used to remove the motion ambiguities and achieve inversion. However, the actual regularizers proposed in the literature are all insufficient since they impose in a small spatial neighborhood a prior smoothness which improperly describes the regularity of turbulent flows [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] [8] [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF] [START_REF] Horn | Determining optical flow[END_REF] [20] [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. Moreover, going further and assuming a set of proper models is given, choosing the most appropriate one and fixing its hyper-parameters is a tricky problem, which has only been studied for Gaussian models in the optic flow literature.

In this paper, we propose a novel motion estimation procedure dedicated to turbulent flows dealing with such limitations. Our work exploits simultaneously two ideas: i) the definition of physically-consistent priors for fluid flows relying on scale invariance properties characterizing turbulent flows; ii) the use of Bayesian modeling for robust optic flow estimation.

Our approach for the definition of physically-consistent priors has been driven by the important advances made in the field of statistical modeling of turbulence since the precursor work of Kolmogorov in 1941 [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds number[END_REF]. Kolmogorov hypothesized, using exact properties followed by the fundamental Navier-Stokes equations, that one can expect some universal scaling properties on the empirical variance of motion increments at different scales. Although limitations to his universality assumption are now known, his theory agrees up to a good degree with observations from experiments and numerical simulations. Therefore, this point of view of turbulence has been the starting point of many recent studies [START_REF] Frisch | Turbulence : the legacy of A.N. Kolmogorov[END_REF]. In this paper, we exploit these results to design a family of physicallysound priors for our motion estimation problem.

We place the problem of model selection within a Bayesian framework. Bayesian analysis has been intensively studied in the past for hyper-parameter estimation and for model selection [START_REF] Jaynes | Bayesian methods: General background[END_REF] [28] [START_REF] Molina | Bayesian and regularization methods for hyperparameter estimation in image restoration[END_REF] [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. However, in optic flow estimation, state-of-the-art inference techniques remained very limited. Indeed, neither proper model selection, nor hyper-parameter estimation [START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF] with model deviations from Gaussianity, have been considered, except in the very recent work on Bayesian inference for optic flow [START_REF] Heas | Bayesian inference of models and hyper-parameters for robust optic-flow estimation[END_REF]. Such non-Gaussian models are nevertheless very common in computer vision, for instance to cope with observation outliers due to noise or varying lighting conditions. In this paper, we propose a Bayesian methodology for the selection of proper models in the context of fluid flow estimation from image sequences. The proposed inference is a generalization of preliminary work [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF], [START_REF] Heas | Power laws and inverse motion modeling: application to turbulence measurements from satellite images[END_REF] to non-Gaussian model of unknown hyper-parameters.

The paper is organized as follows. In section III, we present the problem of fluid motion estimation from an image sequence and highlight the paper's contributions. In section IV, we introduce the scale invariance properties derived from the physics of fluids. Then, in section V, we define and solve a constrained motion estimation problem where the motion estimate is forced to possess these scale invariance properties. In section VI, by reformulating the constrained estimation problem from a Bayesian perspective, we show how we can achieve model selection and hyper-parameter estimation. Finally, the experimental evaluation of the Bayesian method is presented in section VIII.

II. NOTATIONS

Except if otherwise specified, the notational conventions adopted in this paper are as follows. Italic lowercase indicates a scalar quantity, as in a; boldface lowercase indicates a vector quantity, as in a; the kth element of vector a is denoted a(k); capital boldface letters indicate matrices, as in A; the element corresponding to the ith row and jth column of A is denoted as A(i, j); the trace (resp. determinant) of A writes trA (resp. detA); R(A) is the range of the columns of A; we will use the notation Λ a to define a diagonal matrix whose elements are those of vector a; calligraphic letters, e.g., A, represent the set of values that a variable or a vector can take on; capital normal letters, as A, denote random variables. The distribution of a random variable A evaluated at a is denoted p A (a); when clear from the context, we will usually use the short-hand notation: p A (a)

p(a). The gradient and the Laplacian operators are respectively defined by ∇ s A(s , t )

[ ∂A ∂s(1) (s , t ) . . . ∂A ∂s(n) (s , t )] T and by ∆ s A(s , t )

∂ 2 A ∂s(1) 2 (s , t ) + . . . + ∂ 2 A ∂s(n) 2 (s , t ).
The symbol ∝ will denote equality up to a multiplicative factor.

III. PROBLEM AND CONTRIBUTIONS

A. Problem formulation

Estimating the motion of a physical system from the observation of a sequence of images is a well-known problem in the computer vision community. At the heart of state-ofthe-art motion estimation algorithms is the definition of an observation and a prior model on the sought velocity field. The observation model relates the motion of the physical system to the spatial and temporal variations of the image intensity. The prior model defines the spatio-temporal constraints that the motion has to satisfy. In a nutshell, state-of-the-art algorithms intend to find the motion field as the "best" compromise between these two models. We provide a formal description of this approach hereafter.

Let I : (s, t) ∈ S × T → I(s, t) ∈ R be an image intensity function where S ⊆ R 2 (resp. T ⊆ R) is the image spatial (resp. temporal) domain. Let the motion of the physical system be defined by a function v : (s, t) ∈ S × T → v(s, t) ∈ R 2 which associates a two-dimensional motion vector to every spatio-temporal position. Many motion estimation algorithms can be regarded as procedures intending to solve a constrained optimization problem of the form:

v = arg min v f (v, I) subject to h(v) = 0. (1)
For example, one standard way to relate v(s) to I(s, t) is by means of the well-known "Optical Flow Constraint" (OFC),

g(v, I, s, t) = ∂I ∂t (s, t) + ∇ T s I(s, t) v(s, t) = 0, ∀ ∈ S × T (2)
which is valid under intensity conservation along motion trajectories. Quadratic penalization of OFC discrepancies yields [START_REF] Horn | Determining optical flow[END_REF]:

f (v, I) = T S g(v, I, s, t) 2 ds dt. (3)
For other configurations, many other models have been proposed in the literature to relate the image intensity function to the sought motion fields [START_REF] Corpetti | Dense estimation of fluid flows[END_REF], [START_REF] Liu | Fluid flow and optical flow[END_REF]. h(v) defines a set of constraints that have to be satisfied by the sought velocity field.

In the standard computer vision literature, dealing with the motion analysis of natural scenes, these constraints are usually chosen to enforce some spatial coherence of the solution, see e.g., [START_REF] Horn | Determining optical flow[END_REF], [START_REF] Lucas | An iterative image registration technique with an application to stereovision[END_REF]. In the problem considered in this paper, we have the advantage of a perfect knowledge of the equations governing the fluid motion, i.e., the Navier-Stokes equations. Hence, a perfect characterization of the velocity field for incompressible flows should, in principle, solve a problem of the form of (1) with the following constraints:

∇ T s v = 0 ∂(ρv) ∂t + (v T ∇ s )v = -1 ρ ∇ s p + ν∆ s v + F (4)
where v denotes here an unknown 3-dimensional velocity field, p and ρ are scalar variables representing the 3-dimensional pressure and density fields, F and ν represent an external force acting on the fluid flow and its cinematic viscosity. Unfortunately, this choice of constraints leads to an optimization problem which is computationally very demanding and often under-determined since p, ρ, F and the motion component of v perpendicular to the image plane are usually unknown [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. Some recent works propose to alleviate this problem by using simplified constraints of the Navier-Stokes equations for global [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF] or sequential [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF] [35] regularization. Alternatively, in this paper we propose a novel physicallysound prior model on the velocity field and derive an efficient algorithm to solve the corresponding optimization problem. Let us remark that the proposed prior can in principle be combined with the former simplified dynamical constraints. In practice, the estimation of the motion field is performed over a finite set S r × T r of S × T . Hereafter, except if otherwise specified, we will consider the scenario where the subset T r only contains one elements and, accordingly, we will drop the velocity field temporal index. Moreover, we will assume that S r collects the spatial samples on a two-dimensional cartesian grid. We will then use the notation v to denote the m-dimensional vector made up of the concatenation of v(s)'s ∀ s ∈ S r , row by row of the cartesian grid. Hence, in the sequel we will consider an optimization problem with N constraints of the form:

v = arg min v f (v, I) subject to h i (v) = 0, i ∈ [1, ..., N ]. (5)
The data term, obtained by discretizing (2) on the grid S r , has the following quadratic form:

f (v, I) = Av + b 2 2 (6)
where A and b represent respectively a rectangle matrix and a vector depending on the image intensity function I.

B. Overview on the proposed solution

Instead of using the original Navier-Stokes equations describing a fluid flow trajectory, we will rather rely on a simple set of constraints {h i (v), i ∈ [1, ..., N ]} which only involve the instantaneous flow v. These constraints rely on theoretical and experimental results obtained as the continuation of Kolmogorov's works at the beginning of the 20 th century. They describe some scale invariance properties of turbulent flows.

We will see in section IV, that these constraints can be defined in the following quadratic form:

h i (v) = v T A i v -c i , i ∈ [1, ..., N ], (7)
where A i denotes a semi-positive definite matrix and where c i ∈ R + denotes a scalar parameter. In the sequel, we will consider problem [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF] where f (v, I) and h(v) are defined in ( 6) and [START_REF] Carlier | Report 1 on production and diffusion of fluid mechanics images and data[END_REF]. This leads to a particular instance of quadratic programming under quadratic constraints (QPQC) optimization problem. In section V, we will then consider a dual optimization algorithm to solve this problem. Since constraints [START_REF] Carlier | Report 1 on production and diffusion of fluid mechanics images and data[END_REF]

( , s, θ) ∈ R + × S × [0, 2π] → δv ( , s, θ) ∈ R, which is given by: δv ( , s, θ) [v(s + n θ ) -v(s)] T n θ , (8)
where the scalar represents a spatial scale and where n θ denotes the unitary vector in the direction θ. The increment function ( 8) is thus defined as a velocity difference vector between two points distant of in the direction of vector n θ , which is then projected onto vector n θ . Let us give an example to illustrate the definition. We denote the two components of the discrete velocity field by v = (v T 1 , v T 2 ) T , where each component is a m-dimensional vector v i = (v i (1), ..., v i (m)) T . Then, if we consider that an angle θ = 0 corresponds to a vector n θ=0 collinear to the velocity field component v 1 , we have:

δv (1, s, 0) = v 1 (s(1) + 1, s(2)) -v 1 (s) (9) δv (1, s, π/2) = v 2 (s(1), s(2) + 1) -v 2 (s). (10)
The spatial average of the velocity increment function raised to some power p, so-called the p th longitudinal structure function, is given by:

δv ( ) p = 1 2π|S| S [0,2π] δv ( , s, θ) p dθds, (11)
where |S| represents the surface of spatial domain S. Since a finite grid S r with |S r | nodes is considered, we need to approach the spatial integral by a sum over the finite subset S r , while the integral with respect to θ can be approximated by a sum over the finite angular subset A r ⊂ [0, 2π] of |A r | elements:

δv ( ) p ≈ 1 |A r ||S r | s∈Sr θ∈Ar |δv ( , s, θ)| p . ( 12
)
For the particular case p=2, the longitudinal structure function is proportional to the two-point cross-correlation function of the velocity, and therefore related (through the convolution theorem) with the power spectrum of the velocity field. Both the structure functions and the power spectrum are often used in experimental and numerical studies to quantify scaling properties of turbulent flows. Moreover, in fluid mechanics, several exact relations for the third-order longitudinal structure function (p = 3) have been demonstrated for turbulent flows [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]. In particular, for three-dimensional isotropic and homogeneous flows, Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds number[END_REF] demonstrated from the Navier-Stokes equations that the third-order moment is linear with respect to scale, i.e., δv ( ) 3 ∝ . This linear relation theoretically arises in the so-called inertial scale range of the flow, that is to say from the smallest scale of turbulence, also known as Kolmogorov scale, to the scale of the energycontaining vortices of the flow. Analogously, for bidimensional turbulence, Kraichnan [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF] showed that there exists at small scales a phenomenology implying that δv ( ) 3 ∝ 3 . More generally, for any turbulent flow, the p th -order structure function can be described in some scale range I ⊂ R + , so-called inertial range, by a power-law:

δv ( ) p = γ p ζp , ∈ I, (13)
where γ p and ζ p denote the prefactor and the exponent of the power-law.

B. Constraints definition

We now propose constraints for the motion estimation problem [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF], derived from [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF]. We choose to consider the secondorder structure function δv ( ) 2 , i.e., our set of constraints will be derived from the relation:

δv ( ) 2 -γ 2 ζ2 = 0, ∈ I, (14)
for sufficiently small. In practice, only a finite number of longitudinal velocity increments can be computed from v since the latter vector collects the values of the velocity field v(s) on a discrete spatial grid S r . It can then easily be observed that ( 14) can be rewritten in the following matrix form:

v T D T D v -γ 2 ζ2 = 0, ∈ I, (15)
for some compatible with the considered spatial grid. For example, in the case of = 1, a 2 × 2 spatial grid S r and periodic boundary conditions, using the subset A r = {0, π 2 }, i.e., motion increments defined in ( 9)-( 10), we obtain:

δv (1) 2 -γ 2 1 ζ2 ≈ 1 |A r ||S r | Sr θ∈Ar |δv (1, s, θ)| 2 (16) =v T D T 1 D 1 v -γ 2 (17)
with

D 1 = 1 √ 8     -1 1 0 0 0 -1 1 0 0 0 0 -1 1 1 0 0 -1 -1 0 1 0 0 -1 0 1 0 1 0 -1 0 0 1 0 -1     .
For experimental evaluation, in section VIII, we will choose the following discretization subsets:

A r = {-π, -π 2 , 0, π 2 } for ∈ N + ∩ I and A r = {-3π 4 , -π 4 , π 4 , 3π 4 } for ∈ √ 2N + ∩ I.
Obviously, choosing larger subsets S r and A r would result in more accurate approximations of integral [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF].

Before concluding this section, let us point out that parameters γ 2 , ζ 2 are usually unknown. In practice, these parameters should therefore be estimated. We will however postpone this problem to section VI-C. In the next section, we will assume that γ 2 , ζ 2 are perfectly known and will address the problem of motion estimation under constraints [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF].

V. FLUID FLOW ESTIMATION UNDER PHYSICAL CONSTRAINTS

A. A dual approach for constrained motion estimation

In this section, we address the constrained motion estimation problem:

v = arg min v f (v, I) (18)
subject to v T D T D v -γ 2 ζ2 = 0, ∀ ∈ I,
where I denotes the set of constrained scales. In this section, we will particularize our development to the case of a quadratic objective function, i.e., f (v,

I) = Av + b 2 2 . ( 19
)
Problem ( 18) thus belong to the general family of QPQC optimization problems. The extension to the more general case where f (v, I) has the structure of an "M-estimator" will be discussed in section V-B. We consider a dual optimization algorithm, where the solution of ( 18) is searched as the minimum of a Lagrangian function:

L(v, λ) = Av + b 2 2 + λ v T D T D v -γ 2 ζ2 .
The equivalence between the points minimizing the Lagrangian, i.e., vλ = arg min

v L(v, λ), (20)
and the solution of ( 18) is not straightforward and requires additional conditions [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF]. In particular, vλ is a solution of (18

) if vT λ D T D vλ -γ 2 ζ2 = 0, ∀ ∈ I, (21)
λ = arg max λ∈D q(λ), (22)
where q(λ)

inf v L(v, λ)
is the dual function of ( 18) and D {λ|q(λ) > -∞} is its domain.

The optimization procedure considered in this paper is then as follows:

1) Compute λ as the solution of ( 22), 2) Evaluate vλ from (20), 3) Check that vλ satisfies conditions [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds number[END_REF]. We discuss hereafter each step of the procedure:

1) Resolution of the dual problem [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF]: The dual function q(λ) is concave and its domain is convex [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF]. Moreover, it can readily be shown that q(λ) is differentiable in the interior of its domain. We can thus access the global maximum of q(λ) with a gradient ascent method1 :

λ (n+1) = λ (n) + α (n) ∇ λ q(λ (n) ), (23)
where α (n) is a "properly-chosen" step factor and the components of the gradient write

∂q(λ) ∂λ = -v T λ D T D vλ + γ 2 ζ2 . (24)
Note that a global maximizer of q(λ) may not exist, i.e., there may be no λ such that q(λ ) = sup λ q(λ).

This scenario may for example occur when D is an open or unbounded set2 . We have therefore to distinguish between several cases. If a global maximizer λ exists and is an interior point of D, the (unique) fixed point of the gradient-ascent procedure corresponds to the global maximizer λ and satisfies

∇ λ q(λ ) = 0. ( 26
)
If λ exists but is not an interior point of D, a gradient procedure can still access the maximum but the latter does not necessarily satisfy [START_REF] Liu | Fluid flow and optical flow[END_REF].

Finally, if a global maximizer does not exist, the dual optimization procedure [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] fails in solving problem ( 18) since [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF] can not be satisfied. We can however find an approximation of the solution. To do so, we proceed exactly in the same way as if a global maximizer exists but return instead the value of λ (n) achieved when ∇ λ q(λ (n) ) drops below a given threshold. With a slight abuse of notation, we also denote the latter value as λ in the sequel. We will discuss the quality of the resulting approximation below.

Note that resorting to an approximation of the global maximizer does not lead to any particular problem for the rest of the optimization procedure. In particular, as will be shown shortly, a solution to [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] exists for any λ (n) provided that λ (0) ∈ D.

2) Resolution of the dual problem (20): First, let us emphasize that any feasible λ (i.e., λ ∈ D) has to satisfy

A T A + λ D T D 0, Ab ∈ R(A T A + λ D T D ), (27)
where denotes the fact the matrix is semi-positive definite. Indeed, if [START_REF] Lucas | An iterative image registration technique with an application to stereovision[END_REF] does not hold for some λ, it can be seen that

q(λ) inf v L(v, λ) = -∞, (28)
so that λ / ∈ D. From this observation, it follows that ∀λ ∈ D, problem [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] has an analytic solution, namely:

vλ = -A T A + λ D T D † Ab, (29)
where † denotes the Moore-Penrose pseudo inverse of the matrix.

Coming back to the gradient procedure [START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF], we can note that, if λ (0) ∈ D, then λ (n) ∈ D ∀n since the gradient algorithm [START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF] can only increase the value of q(λ) throughout the iterations. The solution of ( 20) is therefore well-defined for any value of λ generated by the gradient algorithm [START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF].

3) Constraint satisfaction (21): Let us first consider the case where a global maximizer of q(λ) exists and is an interior point of D. Then, the unique fixed point of the gradient algorithm must satisfy [START_REF] Liu | Fluid flow and optical flow[END_REF]. Hence, taking the expression of the gradient elements [START_REF] Lavoie | Effects of initial conditions in decaying turbulence generated by passive grids[END_REF] into account, we have that ( 21) is necessarily satisfied.

In the general case (i.e., if a global maximizer of q(λ) does not exist or is not an interior point), v λ is not necessarily ensured to satisfy constraints [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds number[END_REF]. If needed, it is possible to compute a feasible solution and assess its quality as follows. We first project v λ onto the set of feasible points and obtain a feasible solution v ⊥ λ . We then obtain the following bounds on the optimal solution f (v ),

f (v ⊥ λ ) (a) ≥ f (v ) (b) ≥ sup λ q(λ) (c) ≥ q(λ ), (30)
by using the following observations:

(a) follows from the fact that v is a global minimizer of f over all feasible points, (b) is a consequence of the Weak Duality Theorem [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF], (c) follows from the definition of the supremum. The left-hand (resp. right-hand) side of ( 30) provides an upper (resp. lower) bound on the value of f (v ). In other words, knowing v ⊥ λ and λ , we can derive an objective measure of the quality of the proposed solution. Note, in particular, that if λ is an interior point of D, all the inequalities in (30) become equalities.

B. Extension to "M-estimators"

We now extend the procedure described in section V-A to the case where the objective function f is an "M-estimator":

f (v, I) = ρ(Av + b). (31)
This type of function is broadly used in the computer vision literature and enjoys several desirable properties. In particular, it can be shown [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF] that any M-estimator ρ can be expressed as the solution of a convex optimization problem, i.e., ρ(y) = min

z y T Λ z y + ψ(z), (32)
where ψ(z) is the Fenchel-Legendre dual function of ρ. Thus, problem (18) can be rewritten as

(v , z ) = arg min v,z (Av + b) T Λ z (Av + b) + ψ(z) (33) subject to v T D T D v -γ 2 ζ2 = 0, ∀ ∈ I,
and a block-coordinate minimization procedure can be used to access to the solution, i.e.,

z (n+1) = arg min z (Av (n) + b) T Λ z (Av (n) + b) + ψ(z) (34)
v (n+1) = arg min v (Av + b) T Λ z (n) (Av + b) (35) subject to v T D T D v -γ 2 ζ2 = 0, ∀ ∈ I,
Problem (34) has the following analytic solution [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF] From an initial point z (i=0) , iterate until convergence: • Reach the fixed point λ solution of (22) by gradient ascent; i.e., from an intitial point λ (n=0) ∈ D, iterate until convergence:

gradient computation [START_REF] Lavoie | Effects of initial conditions in decaying turbulence generated by passive grids[END_REF] with v (i)

λ (n)
given by ( 29) at z (i) gradient ascent step [START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF]:

λ (n+1) = λ (n) + α (n) ∇ λ q(λ (n) ) • Compute v (i) λ using (29) • Compute z (i+1) using v (i) λ in (36) TABLE I FLOW CHART OF THE OPTIMIZATION ALGORITHM z (n+1) = 1 2 Λ -1 Av (n) +b ∇ρ(Av (n) + b), (36)
which is typically easy to evaluate. Problem ( 35) is a quadratic optimization problem with quadratic equality constraints. Hence, a solution can be searched by applying the procedure described in section V-A.

A flow chart of the overall optimization procedure is summarized in Table I VI. A BAYESIAN FRAMEWORK FOR SELECTING PHYSICAL

CONSTRAINTS

In this section, we address the problem of estimating the parameters γ 2 , ζ 2 of model ( 14). Estimation will be addressed by exploiting Bayesian model selection theory.

A. Bayesian estimation theory

Let us first recall the basics of Bayesian estimation theory [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. The latter provides nice answers to the general problematic of inference of the value of some vector θ 0 from some observations b, which may be noisy and incomplete. To that aim, this theory makes the fundamental assumption that the observations, but also the unknown parameter θ 0 , have a probabilistic nature, i.e., are described by inter-dependent random variables. The Bayesian approach then builds an estimator exploiting the knowledge of the conditional law p(θ 0 |b), so-called a posteriori probability distribution. The estimator noted θ0 is designed according to the optimization of some criteria. Common Bayesian estimators are for example the conditional-mean estimator: θ0 = θ 0 p(θ 0 |b)dθ 0 ; or the well-known Maximum A Posteriori (MAP) estimator: θ0 = arg max θ0 p(θ 0 |b). However, in practice, the observations b rarely only depend on θ 0 but often also on a set of other variables {θ 1 , θ 2 , ...}. Let us restrict ourselves to the set of extra variables {θ 1 , θ 2 }. By making analogous assumptions on their probabilistic nature, it is possible to design the joint posterior probability distribution p(θ 0 , θ 1 , θ 2 |b). Thus, the construction of the posterior p(θ 0 |b) now requires to integrate the joint probability over θ 1 and θ 2 :

p(θ 0 |b) = p(θ 0 , θ 1 , θ 2 |b)dθ 1 dθ 2 . ( 37
)
The random variables θ 1 and θ 2 often live in high dimensional spaces. Consequently, the computation of the latter integral may be quite cumbersome. Trackable alternatives consist in relying on: 1) the joint MAP solution of the problem:

arg max θ0,θ1,θ2 p(θ 0 , θ 1 , θ 2 |b), (38)
2) the evidence framework [START_REF] Mackay | Bayesian interpolation[END_REF] which solves the following maximization problem:

arg max θ0 p(θ 0 |θ map 1 , θ map 2 , b), (39)
where

   θ map 1 = arg max θ1 p(θ 1 |θ map 2 , b) θ map 2 = arg max θ2 p(θ 2 |b). (40)
It is well known that the joint MAP criterion (38) yields an estimate of θ 0 overfitting the observations [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. The evidence framework usually provides a better approximation of the marginalized MAP estimate of θ 0 , i.e., the maximizer of [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF].

In particular, for Gaussian variables, solving (39) provides the exact marginalized MAP estimate [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. However, the two maximization problems in (40) may remain potentially difficult due to high-dimensional non-Gaussian integrals. Nevertheless, as we shall see, there exist trackable algorithms yielding relevant approximated solutions for the evidence framework. Thus, in the following, we will stick to this estimation framework.

B. Probabilistic model and inference problem

Now, coming back to our specific problem, we consider a joint probability distribution on the following random variables: image observations b, the motion field variable v and parameters z, β, γ 2 and ζ 2 . The purpose of the introduction of a new parameter β will be clarified latter. We propose to use the evidence framework (39) in order to solve a problem of the form:

arg max v,z p(v, z|β map , γ map 2 , ζ map 2 , b), (41)
where β map , γ map are the marginalized MAP estimates given by: To this aim, we propose the following factorization:

β map = arg max β p(v, z, β|γ map 2 , ζ map 2 , b)dvdz, (42)
(γ map 2 , ζ map 2 ) = arg max γ2,ζ2 p(v, z, β, γ 2 , ζ 2 |b)dvdzdβ. (43
p(b, v, z, β, γ 2 , ζ 2 ) = p(b|v, z, β)p(v|γ 2 , ζ 2 )p(z|β)p(γ 2 , ζ 2 , β) (44)
and we define the following probabilistic models:

p(b|z, v, β) ∝ exp - (b + Av) T Λ z (b + Av) 2β -1 , (45)
p(z|β) ∝ exp - ψ(z) 2β -1 det(βΛ z ) -1 , (46)
p(v|γ 2 , ζ 2 ) ∝ exp - v T Γ(γ 2 , ζ 2 ) -1 v 2 , ( 47
) p(β, γ 2 , ζ 2 ) = p(β)p(γ 2 )p(ζ 2 ). (48)
In order to obtain factorization (44), we have exploited obvious conditional independences: p(b|v, z, β) = p(b|v, z, β, γ

Γ(γ 2 , ζ 2 ) (β λ (γ 2 , ζ 2 )D T D ) -1 , (49)
where λ (γ 2 , ζ 2 ) denotes the solution of ( 22) when the value of the parameters appearing in the constraints in [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF] are equal to γ 2 , ζ 2 . To alleviate notations, we will omit in the following the dependence of the covariance to γ 2 , ζ 2 . Now, we need to clarify the choice of models ( 44)-(48). To this aim, let us emphasize the equivalence between the dual formulation [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] of optimization problem [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF] and the MAP problem (41). We consider the general case, where the cost function f is an Mestimator. The Lagrangian function corresponding to (33) then writes:

L(v, z, λ ) =(Av + b) T Λ z (Av + b) + ψ(z) + λ (v T D T D v -γ 2 ζ2 ). (50)
Then, using the probabilistic models (45)-( 47), it is easy to see (by direct substitution) that

(v λ , ẑλ ) arg min v,z L(v, z, λ ), (51)
= arg max v,z log p(b, v, z|β, γ 2 , ζ 2 ). (52)
Hence, relying on model ( 45)-(47), we obtain that the solution of the Lagrangian optimization problem (51) and of problem (52), i.e., of the MAP estimation problem (41), are identical.

The MAP estimation problem is conditioned by the parameters β, γ 2 , ζ 2 . Note however that the variable β does not affect the equivalence between (51) and (52) since it represents a multiplicative constant. Therefore, the problem of estimating (v λ , ẑλ ) does not depend on β. However, we will see in the following that this variable is important for model selection.

C. Estimation method

Considering the MAP interpretation of the Lagrangian optimization problem, we now propose to infer γ 2 and ζ 2 by solving the marginalized MAP problem (43). Nevertheless, we avoid the complex calculation of the non-Gaussian integral with respect to z by substituting (43) with the following problem:

(γ map 2 , ζ map 2 ) = arg max γ2,ζ2 log p(γ 2 , ζ 2 |ẑ λ , b) (53)
conditioned from (36) by: ẑλ

1 2 Λ -1 Av λ +b ∇ρ(Av λ + b). (54)
The probability to be maximized in (53) is defined through two successive marginalizations of the joint distribution:

p(γ 2 , ζ 2 |ẑ λ , b) p(β, γ 2 , ζ 2 |ẑ λ , b)dβ, (55)
p(β, γ 2 , ζ 2 |ẑ λ , b) p(β, γ 2 , ζ 2 , v|ẑ λ , b)dv. ( 56
)
One can provide an exact analytical expression of the integral (56) over v since the integrand is a m-dimensional Gaussian. The Laplace's method yields an approximation to the monodimensional integral (55) over variable β. The Laplace's method approximates the integral of a function by fitting a Gaussian at its maximum and computing the volume under the Gaussian. For a uni-dimensional variable x and a function f (x), the Laplace's approximation reads:

f (x)dx f (x) √ 2π - ∂ 2 ∂x 2 log f (x) -1/2 |x , (57)
where x = arg max x f (x). Hence, in order to evaluate (55) we consider the marginal approximation:

p(γ 2 , ζ 2 |ẑ λ , b) = p(b|ẑ λ , β, γ 2 , ζ 2 ) p(β, γ 2 , ζ 2 ) dβ, p(b|ẑ λ , βλ , γ 2 , ζ 2 ) √ 2πσ β ( βλ )p( βλ , γ 2 , ζ 2 ), (58)
where βλ = arg max

β log p(b|ẑ λ , β, γ 2 , ζ 2 ). ( 59
)
σ β ( βλ ) = - ∂ 2 ∂β 2 log p(b|ẑ λ , β, γ 2 , ζ 2 ) -1/2 | βλ (60)
Note that, contrarily to problem (52), estimate βλ is no longer a multiplicative constant and appears in the evaluation of the marginal probabilities. In summary, in order to solve problem (53), we search the element of a finite set Υ ⊂ R 2 of parameters (γ 2 , ζ 2 ) maximizing model probability p(γ 2 , ζ 2 |ẑ λ , b), using the following procedure:

(γ map 2 , ζ map 2 ) = arg max (γ2,ζ2)∈Υ log p(b|ẑ λ , γ 2 , ζ 2 )p(γ 2 , ζ 2 ) (61)
where, in order to evaluate p(b|ẑ λ , γ 2 , ζ 2 ) for some (γ 2 , ζ 2 ) ∈ Υ, we solve:

(v λ , ẑλ ) = arg max v,z log p(b, v, z|β, γ 2 , ζ 2 ), (62)
βλ = arg max β log p(b|ẑ λ , β, γ 2 , ζ 2 ). (63)
We detail hereafter the strategy we considered to implement each step of the proposed algorithm:

1) Robust constraint motion estimation (62): Given some power-law parameters (γ 2 , ζ 2 ), the problem is equivalent to the constraint optic flow estimation problem (18) discussed and solved in section V. It can be solved by the proposed dual optimization algorithm, which uses iterative conditional maximizations to estimate (v λ , ẑλ ), and a gradient algorithm to converge towards λ related to (γ 2 , ζ 2 ).

2) Parameter estimation (63): Given the power-law parameters (γ 2 , ζ 2 ), the estimate βλ is obtained by canceling out the gradient of the marginalized logarithm likelihood log p(b|ẑ λ , β, γ 2 , ζ 2 ). Let • to denote the expectation with respect to p(v|b, ẑλ , β, γ 2 , ζ 2 ). By noticing that:

∇ β log p(b|ẑ λ , β, γ 2 , ζ 2 ) = ∇ β log p(b, v|ẑ λ , β, γ 2 , ζ 2 ) p(v|b, ẑλ , β, γ 2 , ζ 2 )∇ β log p(b, v|ẑ λ , β, γ 2 , ζ 2 )dv (64)
and using the definition of the joint distribution in (45)-(47), the gradient of the marginalized log likelihood can be efficiently computed as:

∇ β log p(b|ẑ λ , β, γ 2 , ζ 2 ) = 1 2 β -1 |S r | -(b -Av) T Λ ẑλ (b -Av) . (65)
Then, after some manipulation of (65) described in Appendix I, one can show that the gradient vanishes for:

βλ = |S r | -tr( A T Λ ẑλ A + λ D T D -1 A T Λ ẑλ A) (b -Av λ ) T Λ ẑλ (b -Av λ ) . ( 66
)
βλ is the unique global maximum since it can be shown that the function log p(b|ẑ λ , β, γ 2 , ζ 2 ) is concave. The computation of (66) only requires tractable linear operations and the MAP estimates (v λ , ẑλ ), which are the solutions of (62). Details are given in section VII.

3) Prior model inference (61): The decision on the prior model parameter (γ 2 , ζ 2 ) is made by maximizing log p(b|ẑ λ , γ 2 , ζ 2 ). Since we assume that (γ 2 , ζ 2 ) takes its values in a finite set Υ, solving (61) only requires the evaluation of log p(b|ẑ λ , γ 2 , ζ 2 ) for these values. Since p(γ 2 ) and p(ζ 2 ) are assumed to be non-informative priors, the logarithm of the Laplace approximation of this probability (58) reads:

log p(b|ẑ λ , γ 2 , ζ 2 ) ∝ log σ β ( βλ )+log p(b|ẑ λ , βλ , γ 2 , ζ 2 ), (67)
where, after some calculations detailed in appendix II, the logarithm of the marginalized likelihood can be rewritten as:

log p(b|ẑ λ , βλ , γ 2 , ζ 2 ) = 1 2 log det λ D T D 1 2 log det βλ Λ ẑλ - 1 2 (b -Av λ ) T βλ Λ ẑλ (b -Av λ ) - 1 2 log det A T Λ ẑλ A + λ D T D ,
and the logarithm of its Hessian evaluated at βλ can be expressed as:

log σ β ( βλ ) = 1 2 log 2 βλ (b -Av λ ) T Λ ẑλ (b -Av λ )
.

(68)

VII. ALGORITHM COMPLEXITY

The flow chart of the algorithm for the selection of the power-law prior and the estimation of the instantaneous velocity field from an image couple is presented in Table II.

Let κ denote the conditioning number of the posterior Hessian matrices, which have been preconditioned by an incomplete LU decomposition. The complexity of the overall algorithm presented in Table II is of O(κm), where we recall that m denotes the dimension of the motion vector.

More precisely, the complexity of two computation steps dominates the algorithm. First, the complexity of calculating

• For all power-laws (γ 2 , ζ 2 ) ∈ Υ:

-Compute MAP estimate (v λ , ẑλ ) using algorithm of [START_REF] Mininni | Direct simulations of helical hall-mhd turbulence and dynamo action[END_REF], i.e. solving a linear system using for instance a conjugate gradient squared (CGS) algorithm, is of O(κm). The additional use of the incomplete LU preconditioner enables: 1) to lower the conditioning number κ and therefore to reduce the complexity of ( 29), 2) to access with a negligeable complexity to the determinants appearing in (67). Indeed the determinant is given by the product of the diagonal elements u ii of the upper triangular matrix U:

det(LU) = i u ii , ( 69
)
where the incomplete LU decomposition is performed on the posterior and prior Hessian matrix. Second, a straightforward calculation of the trace appearing in (66) can be very demanding. Instead, we rely on a trace randomization technique [START_REF] Skilling | The eigenvalues of mega-dimensional matrices[END_REF] which is efficiently achieved using the CGS algorithm. Details are given in Appendix I. The method complexity for the calculation of the trace is also O(κm)

VIII. EXPERIMENTAL EVALUATION A. Simulated 2D turbulence

The Bayesian method is first assessed on a synthetic image sequence database of numerically generated 2D homogeneous, isotropic and incompressible turbulent flow fields [START_REF] Carlier | Report 1 on production and diffusion of fluid mechanics images and data[END_REF]. This database has proven its relevance in previous studies in the field of fluid flow estimation [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF], [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF], [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF], [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. It relies on a Direct Numerical Simulation (DNS) of Navier-Stokes equations at Reynolds number Re = 3000. The output of this simulation is then used to compute the Lagrangian trajectories of particle tracers transported by the flow or, in the case of scalar images, to solve an advection-diffusion equation. The latter equation reads:

∂I ∂t (s, t) + ∇ T s I(s, t) v(s) = κ∆ s I(s, t), ∀s ∈ S r , ∀ t ∈ T r , (70)
where κ ∈ R + denotes a diffusion coefficient. Motion and scalar transport equations are solved in Fourier space using dealiased Fourier expansions in two directions with periodic boundary conditions. In order to sustain turbulence, the flow was stirred by an external random force in Fourier space acting at low wavenumbers. The time integration is done using a third-order Runge-Kutta scheme. The simulation resulted in two synthetic sequences of 100 images of 256 × 256 pixels depicting scalar or particles, together with the true underlying motion fields. A precise description of the simulation can be found in [START_REF] Carlier | Report 1 on production and diffusion of fluid mechanics images and data[END_REF]. Fig. 1 presents three samples of the synthetic scalar or particle image sequences and a color representation of the underlying velocity fields. In these visualizations color Synthetic data set:

Motion estimates:

t = 30∆t t = 60∆t t = 90∆t
Fig. 1. 2D turbulence: synthetic data set and motion estimates. From top to bottom at time t = 30∆t (left), t = 60∆t (middle) and t = 90∆t (right): simulated velocity fields visualized using the color system of [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF], scalar and particle images, motion estimates in the case of scalar and particle images . and intensity code vector orientations and magnitudes [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]. In the case of scalar images (resp. particle images), the observation model was the advection-diffusion equation (70) (resp. the optic flow constraint equation ( 2)) penalized by the robust M-estimators approaching the L 1 norm (see e.g., [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF]). A multi-resolution approach was used to cope with large displacements [START_REF] Bergen | A 3-frame algorithm for estimating two-component image motion[END_REF]. We chose to perform model inference relying on motion field increment estimated at the last level of this pyramidal approach.

The inference of the prior power-law model (γ map 2 , ζ map 2 ) is performed for the two experiences, in the scale range of 1 to 4 pixels and at the arbitrary chosen time of t = 20∆t. Fig. 6 shows, for both experiments, the shape (up to an irrelevant constant) of the a posteriori probability distribution of the power-law model with respect to its parameters (γ 2 , ζ 2 ), that is p(γ 2 , ζ 2 |b, ẑλ ). The global maximum of the probability is located for both the cases of scalar or particle imagery, at the coordinates (γ map 2 , ζ map 2 ) =(2.5e-3,1.9). It is important to note that, as shown in Fig. 6, these coordinates also correspond to an "optimal" choice of parameters in the sense of the minimum of MBA error. Moreover, we can observe in Fig. 2 that the inferred power-law describes perfectly the true structure function computed with the ground truth motion field. Finally, the inferred scaling exponent is consistent with the theoretical value of 2 predicted for 2D turbulence in the inertial or dissipative range [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF].

Assuming constant second order statistics along time, the selected power-law model (γ map 2 , ζ map 2 ) is used to process the entire image sequence. The evolutions of the resulting mean barron angular (MBA) error and root mean square (RMS) error are plotted in Fig. 7 and compared to recent state of the art approaches 3 . Note that MBA and RMS errors were computed with divergence-free motion components obtained by Helmoltz decomposition [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF]. Indeed, this decomposition was used here since a subset of state-of-the-art methods code implicitly an additional divergence-free constraint [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] [37], leading possibly to a bias in the regularization scheme evaluation.

The proposed estimation method reaches the level of the best state-of-the-art algorithms in the case of the particle image sequence, while for scalar imagery it outperforms them all. In 3 Let us mention that, for consistency of the evaluation, we have compared the proposed method to algorithms implying neither temporal regularization schemes (see e.g., [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF] [35] [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF]) nor specific configurations of fluid flows (see e.g., [START_REF] Ruhnau | Optical stokes flow estimation: An imagingbased control approach[END_REF] for the particular case of stokes flows). Note that, if necessary, the proposed Bayesian method could in principle be extended without any major difficulties to include an additional temporal regularization term. particular, although second order regularizations [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF] [20] or higher-order regularization schemes [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF] yield some improvements in comparison to standard first order regularization [START_REF] Horn | Determining optical flow[END_REF], there is a significant gap in accuracy with the proposed method in the case of scalar imagery. This gap is due to the weak spatial gradients of scalar images (i.e., poorly conditioned data terms), which make crucial the information brought by an appropriate prior. For particle imagery, the data term is better conditioned and penalization of high-order derivatives is sufficient to obtain relevant turbulent flow approximations. Let us also mention that classical correlation-based techniques or the adaptive correlation method proposed by [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] perform reasonably well in the case of particle images, while such local schemes usually fail in the case of scalar imagery. The method is then evaluated on two synthetic image couples obtained by deformation of instances of the scalar and particle image sequences presented in the last section. The velocity field used for deformation is the plane components of 3D motion fields computed by numerical simulation of 3D homogeneous, isotropic and incompressible turbulence [START_REF] Mininni | Direct simulations of helical hall-mhd turbulence and dynamo action[END_REF] at Reynolds number Re = 667. On the contrary to the incompressible 2D case, the apparent motion field is here characterized by divergence-free and irrotational components of comparable magnitude. Fig. 8 presents the scalar or particle synthetic image couples and a color representation of the underlying velocity field. For both the scalar and particle cases, the data model used in these experiments is the optic flow constraint equation ( 2) penalized by the robust M-estimators approaching the L 1 norm (see e.g., [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF]). A multi-resolution approach is also used to cope with large displacements [START_REF] Bergen | A 3-frame algorithm for estimating two-component image motion[END_REF].

B. Plane components of simulated 3D turbulence

The inference of the prior power-law model (γ map 2 , ζ map 2 ) is performed for the two experiences, in the scale range of 1 to 4 pixels. Fig. 9 depicts, for both experiments, the shape of the a posteriori probability distribution of the powerlaw model parameters (γ 2 , ζ 2 ). The global maximum of the probability is located for the case of particle imagery at (γ map 2 , ζ map 2 ) =(1.5e-3,2.0) and for the case of scalar imagery at (γ map 2 , ζ map 2 ) =(1.5e-3,1.9). It is important to note that, as shown in this figure, these coordinates also correspond to a nearly "optimal" choice of parameters in the sense of the minimum of MBA error for both experiments. Moreover, we can observe in Fig. 3 that the inferred power-law describes almost perfectly the true structure function computed with the ground truth motion field.

A comparison with recent state of the art approaches in terms of MBA and RMS errors is given in the table of Fig. 10. The proposed estimation method reaches the level of the best state of the art algorithms in the scalar case or outperforms them in the case of particle images. Let us mention that, conversely to the case of 2D incompressible turbulence, the adaptive correlation method of [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] performs relatively well in the case of particle images.

C. Experimental 3D turbulence

The Bayesian approach is also assessed on a real image sequence recorded in a wind tunnel 4 by particle image velocimetry. This image sequence shows three-dimensional homogeneous and isotropic turbulent motion at Reynolds number Re = 2300 in the wake of a biplane grid (see [START_REF] Lavoie | Effects of initial conditions in decaying turbulence generated by passive grids[END_REF]). The size of images is 512 × 512 pixels corresponding to 34 × 34 mm 2 . A spatial subset of an image couple is depicted in Fig. 4. A tiny time delay between frames was chosen in order to limit the loss of particles due to the "out of plane" component.

The data model used in these experiments was the OFC equation ( 2) penalized by the robust M-estimators approaching the L 1 norm (see e.g., [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF]). The inference of the power-law parameters (γ map 2 , ζ map 2 ) was performed on an image couple of the sequence in the scale range of 1 to 4 pixels. The global maximum of the model probability is located near the coordinates (γ map 2 , ζ map 2 ) =(1e-3,2.0). The inferred scaling exponent seems to be consistent with physics: the scales of the power-law mainly belongs to the dissipative range which is known to be ruled theoretically by the exponent γ 2 = 2.0. Indeed, the dissipative scales was estimated below 3.3 pixels with hot wire anemometry measurements.

The selected prior power-law yielded highly structured motion fields as shown in the zoom in Fig. 4. The strength of multi-scale regularization is particularly striking in this case: tiny turbulent structures are revealed. Some vortex structures are smaller than 10 pixels diameter. It can be noted that for this flow with energy level sustained across the small scales, classical optical flow methods and correlation approaches estimated only inaccurate vector fields. At larger scales, differences are Fig. 4. 3D experimental turbulence. From top to bottom, spatial subset of the two input images, estimated velocity in color representation and divergencefree component obtained by a Helmoltz decomposition [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] for the proposed method, divergence-free component given by a first-order regularizer [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF] (below). The inferred power-laws parameters are (γ less prominent and the proposed regularization seems coherent with a first-order scheme [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF].

D. Atmospheric turbulence

We finally test the proposed modeling on a satellite image sequence. The analysis was done using METEOSAT Second Generation meteorological images acquired above the north Atlantic Ocean at a rate of an image every 15 min. The spatial resolution is 3 × 3 km 2 at the center of the Earth image disk. Following the methodology in [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF], a set of sparse pressure difference images of 256 × 256 pixels associated to a stack of layers at various altitude was derived. As detailed in the latter work, each image of the stack corresponds to a given altitude and can be related to a vertically averaged horizontal wind field via an observation model relying on layer mass conservation ∀s ∈ S r , ∀ t ∈ T r :

∂I ∂t (s, t) + ∇ T s I(s, t) v(s) + I(s, t) ∇ s • v(s) = 0. (71)
A spatial subset of the very sparse image observations used in these experiments is displayed in Fig. 5. These images represent incomplete pressure difference maps of an atmospheric layer at intermediate altitude. We used in this case the data model proposed in [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF], which describes discrepancies from the mass conservation equation (71) only in observed areas. Discrepancies are given by the Leclerc's robust M-estimators [START_REF] Corpetti | Dense estimation of fluid flows[END_REF].

The inference of the power-law parameters (γ map 2 , ζ map 2 ) was performed in the scale range of 1 to 4 pixels. A multiresolution approach was used again to cope with large displacements. A clear maximum is located at (γ map 2 , ζ map 2 )=(2.5e-5,2.0). According to physical consideration, the inferred scaling law exponent suggests two-dimensional turbulence, i.e., scaling exponent γ 2 = 2.0. Although some controversy remain on the spectrum of atmospheric turbulence at intermediate scales [START_REF] Lindborg | Horizontal velocity structure functions in the upper troposphere and lower stratosphere 2. theoretical considerations[END_REF], this result seems reasonable since the atmosphere can be modeled as a quasi two-dimensional flow in a wide range of scales.

As illustrated in a zoom displayed in Fig. 5, the proposed modeling produced motion fields visually consistent with the sparse observations and realistic in non-observed regions. The velocity field at this resolution appears to be very smooth and structured with vortices of only a few pixels size. Improvements brought by the method is in this case particular striking when compared to the second order regularization scheme proposed in [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF].

IX. CONCLUSION

We have presented a method for fluid flow estimation exploiting physical prior knowledge on the scale invariant properties of turbulent flows. Multi-scale regularization is achieved by imposing some power-law behavior to the second order structure function of motion. Bayesian modeling treating motion fields, prior power-law models together with their hyper-parameters as interdependent random variables allows us to simultaneously provide solutions for three problems: robust motion estimation, hyper-parameter estimation, and power-law model selection. Optimal model and hyperparameter estimates are derived in analytical forms in a marginalized a posteriori sense. Consequently, motion estimation does not involve neither tuning regularization hyperparameter nor choosing an arbitrary power-law model. Experiments on synthetic and real image sequences prove that the proposed Bayesian inference scheme succeeds to select adapted model and hyper-parameters and enhances state of the art in fluid flow estimation. Extending the Bayesian method to the modeling of non-homogeneous motion fields constitutes a promising research direction. Nevertheless, this extension is likely to face complexity bottlenecks, necessitating probably the introduction of variational Bayesian inference schemes [START_REF] Bernardo | The variational bayesian em algorithm for incomplete data: with application to scoring graphical model structures[END_REF]. Horizontal wind at intermediate altitude. From top to bottom, spatial subset of the two input incomplete images, estimated wind field in color representation and divergent-free component obtained by Helmoltz decomposition [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] for the proposed method, divergence-free component given by a second-order regularizer [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF] (below). The inferred power-laws parameters are (γ Kadri12 Proposed Fig. 7. 2D turbulence: evolution with respect to time index, in the case of scalar (above) or particle (below) imagery, of MBA (left) and RMS (right) errors of the divergence-free motion component relative to motion estimate with adaptive correlation [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF], with first order regularization [START_REF] Horn | Determining optical flow[END_REF], with second order regularizations [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF] [20], with higher-order regularization [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF], or with the proposed method. Divergence-free motion components are obtained by Helmoltz decomposition [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] Fig. 8.

3D turbulence: synthetic data set and motion estimates. From left to right, scalar (above) and particle (below) image couple, simulated velocity field (above) and motion estimates for scalar (above) and particle (below) images visualized using the color system of [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]. Becker11 [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] Horn81 [START_REF] Horn | Determining optical flow[END_REF] Derian11 [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF] Kadri12 [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] with first order regularization [START_REF] Horn | Determining optical flow[END_REF] [9], with second order regularizations [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF], or with the proposed method.

2 and ζ map 2

 2

  ) Clearly, system (41)-(43) corresponds to the evidence framework defined in (39)-(40) with θ 0 = {v, z}, θ 1 = {β} and θ 2 = {γ 2 , ζ 2 }. The practical implementation (41)-(43) requires the definition of an explicit joint probability distribution.

Fig. 2 .

 2 Fig.2. 2D turbulence: second order structure function estimation δv ( ) 2 at t=20∆t. Above: inferred power-law model (straight line), true (continuous curves) and estimated (stars and crosses) structure function in horizontalvertical and diagonal direction for particle images. Below: identical legend for scalar imagery.

Fig. 3 .

 3 Fig.3. 3D turbulence: second order structure function estimation δv ( ) 2 . Above: inferred power-law model (straight line), true (continuous curves) and estimated (stars and crosses) structure function in horizontal-vertical and diagonal direction for particle images. Below: identical legend for scalar imagery.

Fig. 5 .

 5 Fig.5.Horizontal wind at intermediate altitude. From top to bottom, spatial subset of the two input incomplete images, estimated wind field in color representation and divergent-free component obtained by Helmoltz decomposition[START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF] for the proposed method, divergence-free component given by a second-order regularizer[START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF] (below). The inferred power-laws parameters are (γ

map 2 , ζ map 2 )Fig. 6 .

 226 Fig. 6. 2D turbulence: comparison at at t=20∆t between minus of the logarithm of unnormalized power-law model probability (left) and MBA error (right) with respect to power-law exponent ζ 2 (y-axis) and prefactor γ 2 (x-axis) for particle (above) or scalar (below) imagery. Iso-contours of decreasing values are plotted around the minimum (γ map 2 , ζ map 2 ) =(2.5e-3,1.9).

Fig. 9 .

 9 Fig. 9. 3D turbulence: comparison between minus of the logarithm of unnormalized power-law model probability (left) and MBA error (right) with respect to power-law exponent ζ 2 (y-axis) and prefactor γ 2 (x-axis) for scalar (above) or particle (below) imagery. Iso-contours of decreasing values are plotted around the MBA error minimum and the probability minimum, that is to say (γ map 2 , ζ map 2 ) =(1.5e-3,2.0) for particle imagery or (γ map 2 , ζ map 2 ) =(1.5e-3,1.9) for scalar imagery.

  2 , ζ 2 ), p(v|γ 2 , ζ 2 ) = p(v|z, β, γ 2 , ζ 2 ) and p(z|β) = p(z|β, γ 2 , ζ 2 ).We consider non-informative probability distributions for p(β), p(γ 2 ) and p(ζ 2 ) in (48). Models (45) and (47) are Gaussian, with β denoting a positive parameter and the covariance matrix:

Table I -

 I Access βλ with (66) -Evaluate the probability p(b|ẑ λ , γ 2 , ζ 2 ) of the power-law (γ 2 , ζ 2 ) with (67)

	map 2 )
	satisfying (61)

• Find among the elements of Υ the MAP estimate (γ map 2 , ζ

TABLE II FLOW

 II CHART OF THE OVERALL ESTIMATION ALGORITHM

When λ (n) is a boundary point of D, the subgradient of q(λ) should in principle be considered. However, we avoid such details hereafter for the sake of keeping the discussion as simple as possible.

It is however not possible to predict the existence of a solution for general QPCP optimization problems.

Irstea research institute, center of Rennes, France

APPENDIX I GRADIENT OF THE MARGINALIZED LIKELIHOOD

We want to compute the gradient given in (65) which writes ∇ β log p(b|ẑ λ , β, γ 2 , ζ 2 ):

where we recall that • is the expectation with respect to p(v|b, ẑλ , β, γ 2 , ζ 2 ). Note that in this expression, v only appears in linear and quadratic forms. As a consequence, the latter derivative is only a function of the mean and covariance of p(v|b, ẑλ , β, γ 2 , ζ 2 ). Now, it is easy to see that this distribution is a Gaussian distribution with mean and covariance defined respectively as:

Let us write the second term of (72) as follows:

The first term on the right hand side of (75) can be rewritten as:

where U and Ω are defined by the eigen-value decomposition of

Now,

Plugging this expression in (76), we obtain

where we use the fact that tr(AB) = tr(BA) for square invertible matrices. Combining (65), ( 75) and (81), we finally get:

The trace operator in (82) applies on a large and sparse matrix. However, it can easily be evaluated by randomization. Let I n represent a n-dimensional identity matrix. More precisely, for J

1 random samples of a n-dimensional normalized and centered Gaussian distribution r j ∼ N (0, I n ), by the weak law of large numbers one gets the the approximation that:

where • denotes here the expectation with respect to the ndimensional Gaussian and where the vector Γr j is the solution provided by the CGS algorithm of the problem Γ -1 X = r j , where X is the unknown.

APPENDIX II EXPRESSION OF THE MARGINALIZED LIKELIHOOD

Introducing the definitions (45)-(47) in the integral (56), we obtain:

Let us rewrite the integrant as

where the posterior covariance Γ has been defined in (74). Now, the last two factors in (86) do not depend on v, and the integration of the first term alone gives:

Therefore, we have:

Since we have [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds number[END_REF], using the definition (73)-(74), this can also be rewritten as:

The prior used here is a degenerated distribution and its Hessian is not a full rank matrix i.e., its determinant is equal to zero. We need to impose some Dirichlet boundary conditions before evaluating the logarithm of its determinant appearing in the first term on the right hand side of (89).

Note that the precise value on the boundaries do not need to be specified. Indeed, considering some Dirichlet boundary conditions results in slightly modifying the form of the Hessian matrix but its content remains independent of the boundary function.