Yin Yang
email: yin.yang@irstea.fr

Dominique Heitz

Etienne Mémin

Lagrangian Particle Image Velocimetry

published or not. The documents may come

Introduction

Accurate and efficient reconstruction of volumetric Eulerian velocity field of fluid flow from particle images is an important means to study the turbulence. To obtain a higher spatial resolution, higher seeding densities of particle must be dealt with. For many years, the dominant approach for volumetric reconstruction is Tomo-PIV (Elsinga et al., 2006) that is able to deal with particle images of density up to 0.05 particles in pixel (ppp). In this approach, first the distribution of 3D voxel intensity discretized from the original particle fields is reconstructed, then the velocity is obtained through the cross-correlation analysis over interrogation volumes (IV). In 2D planar PIV, the cross-correlation algorithm is proved to be robust as well as efficient. This motion analysis strategy contributes to the robustness of Tomo-PIV since the crosscorrelation algorithm is a direct extension from 2D pixel space to 3D voxel space. However, for volumetric reconstruction, the 3D cross-correlation algorithm is computational demanding. In addition, the quality of voxel intensity reconstruction is likely to be tempered by the presence of ghost particles. Ghost particles appear naturally for particle image of high ppp due to the ambiguities associated with the triangulation process. Consequently, the velocity result is less convincing.

Another way of reconstructing 3D Eulerian velocity is first to obtain particle's velocity through Lagrangian particle tracking (LPT), followed by some interpolation algorithm that interpolates the Lagrangian velocity field to Eulerian grid in an optimal manner. Current approaches include FlowFit [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF] and VIC+ [START_REF] Schneiders | Dense velocity reconstruction from tomographic PTV with material derivatives[END_REF]. These approaches have gained growing interests since the introduction of the so called Shake-The-Box (STB) method, a 4DPTV scheme proposed by [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF]. STB is able to reconstruct flow scenes with high particle seeding density compared to standard 3DPTV where particle positions are sought only through triangulation. STB method, rooted in Iterative Particle Reconstruction (IPR) approach [START_REF] Wieneke | Iterative reconstruction of volumetric particle distribution[END_REF], features an active image pattern matching scheme. This scheme searches for the best particle's parameters (position, intensity) iteratively so that the resulting projected image through an explicit image model is most similar to the particle image captured by cameras. Besides, STB introduces a predicting phase for the particle position. Such usage of temporal information reduces the difficulties of the particle tracking procedure and can produce nearly ghost-free particle fields.

After the Lagrangian flow field is obtained, FlowFit or VIC+ post-processing procedures can be used to provide Eulerian vector fields. The optimal interpolation processes of FlowFit and VIC+ minimize the difference between the reconstructed field and the velocity and acceleration measurements, respectively. Furthermore, FlowFit introduces additional physical constraints into its cost function. Whereas VIC+ employes the velocity-vorticity formulation and adoptes an adjoint approach to compute the gradient of the cost function.

More recently, Lasinger et al. ( 2018) proposed a method that jointly reconstruct individual particles as well as volumetric velocity fields directly from particle images. This approach alternates between a particle triangulation procedure (using simple triangulation or IPR) and an optimization procedure that finds the optimal parameter field (particle position, intensity and Eulerian velocity) explaining the variations of particle fields between two consecutive snapshots. This method is a direct 3D extension of variational flow estimation.

In the same vein, we propose in the present study a direct estimation of the Eulerian velocities from Lagrangian particle positions. However, the novelty of our method lies in the following points: First, we measured the data discrepancy term on local patches rather than for individual particle and we also assumed a local homogeneity of the velocity; Then, we solved the optimization problem using a strategy inspired from Ensemble-based Data Assimilation approach (EnDA) [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. This strategy avoids the calculation of gradient of the cost function with respect to the original control variables. The efficiency and accuracy of this method supersedes STB as shown in [START_REF] Yang | An ensemble filter estimation scheme for lagrangian trajectory reconstruction[END_REF] in the context of particle tracking; Finally, we employed a transport equation that propagated the particles to their positions at next time level according to a given velocity candidate on the fixed grid. This particle transport equation played a central role in the formulation because it related the particle position to the Eulerian velocities.

Note that our method can be applied to either time-resolved sequential data or double-frame data. In either case, the triangulation or IPR technique is still needed to provide the prior particle field to initialize the velocity estimation procedure. We evaluate in this paper a double-frame strategy where the IPR was required to build the initial particle field for the first frame as well as to improve this particle field along with the velocity estimation procedure.

Methods

In the following, we formulate the Eulerian velocity estimation task as an optimization problem and propose relevant solution methods. In principe, at frame k -1, we search for the 3D Eulerian velocity field within the flow domain that propagates the particles with known coordinates to the next frame k by minimizing the image discrepancy between the one captured by camera and the one active-generated by OTF. We use sum-square-difference (SSD) measuring the image discrepancy. This strategy translates into following cost function:

J(u) = ∑ i ||I i rec,k -I i proj,k (u)|| 2 , ( 1
)
where i is the camera index and I rec is the image recorded by camera. I proj is an active-projection image model and a function of u. This image model can be developed into a superposition of total Q particles' projection respectively:

I i proj,k = Q ∑ p I i part,k = Q ∑ p OTF i {M i [X k p ]}, (2)
where the OTF is modeled as Gaussian blob [START_REF] Wieneke | Iterative reconstruction of volumetric particle distribution[END_REF][START_REF] Schanz | Non-uniform optical transfer functions in particle imaging: calibration and application to tomographic reconstruction[END_REF], M is the mapping function linking world coordinate X p to camera coordinate x i p . The particle position X k p at frame k is related to the velocity u k-1 through a transport model X k p = T k (u k-1 , X k-1 p ). This transport model computes the pathline of particle p according to the particle velocity v k-1 p that depends on the space and time. Still, we must establish the operator L p such that v p = L p (u). Lasinger et al. ( 2018) used a trilinear interpolation scheme to determine the particle velocity v p from u on nearby grids. We employ a different strategy inspired by the idea of Lucas-Kanade algorithm and we elaborate this strategy in the next section.

Local approach

Our approach follows from an discretization of the whole domain of interest into small 3D local patches similar to the interrogation volumes used in Tomo-PIV. Suppose q particles are bounded in a local patch l, the cost function at frame k reads,

J loc (u l ) = ∑ i ||I i rec,k - q ∑ p=0 I i part,k [T k,p (u l )]|| 2 . ( 3
)
This is a general formulation because it does not show how the particle velocity is computed. We also assume that each patch has independent but homogeneous velocity. This suggests all particles within the patch should have the same velocity. The homogeneous assumption is not required as we can indeed solve above cost function with any interpolation scheme. Also this approach seems less accurate compared to the interpolation. However we adopt this strategy because first, for any realistic flow analysis problem, the resolution of the reconstructed velocity field is limited by the particle density. For example, we have found that our approach is relatively more robust under noise with an adequate size of local patch containing 5 to 10 particles. This empirical parameter coincides with the number of particles contained in an IV on average used in Tomo-PIV. Additional overlapping window techniques can be used to increase the spatial resolution.

Second, the ensemble technique used to solve the optimization problem performs better with smaller control space. E.g. for a single particle, the size of control vector of our strategy is eight times smaller compared to the size of control vector if adopting the trilinear interpolation.

Adding regularizations

One drawback of local approach is that the velocity can not be inferred when no particle is found in the local patch. This can happen when the particle density is not high enough or the particle distribution is not homogeneous enough. One strategy to handle this issue, without the loss of resolution on the reconstructed velocity field, is to consider some regularization on the velocity field. The common choices for the regularization terms include a term penalizing the gradient of velocity component as well as a term penalizing the divergence for incompressible flows.

After optimizing u l for each patch l in Eq.(3), we obtain u L (L denotes the local analysis). We formulate another cost function in terms of u G (G denotes global analysis) taking account of the regularizations on velocity fields:

J tot (u G ) = ||u G -u L || 2 + Ω (α 3 ∑ j=1 ||∇u G j || 2 + β||∇ • u G ||)∂ω (4)
where α and β control the relative weight of corresponding regularization term compared to the first term.

Solution

The system (3) constitutes a non-linear optimization problem that can be solved using the iterative Gauss-Newton algorithm. At iteration m, we can linearize the image model around an estimation of the velocity field ûl m-1 (either deduced from previous iteration or u 0 = u init to initialize the first iteration m = 0). After linearization, we have the following linearized cost function (5):

J(δu l m ) = ∑ i ||I i res,k - q ∑ p=0 ∂ X I i part,k ∂ u T k,p δu l m || 2 , ( 5
)
where ∂ z F denotes the tangent linear model of non-linear model F with respect to its variable z linearized around the prior z m . We also note that,

I i res,k = I i rec,k,U l -∑ q p=0 I i part,k [T k ( Xk-1,p , u l m )],
is the difference, between the recorded image and the active-projected image of particles driven by the prior velocity u m . Consequently, we search for the increment δu l m at each iteration. Inspired by the work of [START_REF] Yang | An ensemble filter estimation scheme for lagrangian trajectory reconstruction[END_REF], in which the authors proposed to solve the Lagrangian particle field estimation problem within the framework of EnDA [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF], we also employ the ensemble technique to solve the optimization problem. The basic idea is to employ an ensemble of velocity candidates that produce an ensemble of particle image patterns through the transport equation, the mapping function and the OTF. Then the optimal velocity can be sought, in form of a linear combination of all candidates, by minimizing a modified empirical cost function. Eventually the original non-linear optimization problem is transformed into a linear one that is easy to solve.

The system (4) is relatively easy to solve since the discretization of gradient or divergence operators on Eulerian grids constitutes eventually a linear system. We use conjugate-gradient algorithm to solve for the final velocity. The details of the whole approach is listed in section 2.5.

Discussion

In Tomo-PIV, we assign E k-1 (X V ) as local 3D voxel intensity field, then the task is to find u minimizing

∑ X V ||E k (W (X V , u)) -E k-1 (X V )||
2 where E k is the 3D voxel intensity field at k and W is a 3D warp mapping the voxel coordinates X V to a new set. Tomo-PIV solves this problem by the cross-correlation technique. It is clear that our approach intends to solve a similar template matching task as TomoPIV without building voxel intensity field. Our solution method uses a random sampling strategy on control variables yielding samples of projected image patches. This strategy is more efficient than the cross-correlation technique regarding computation in 3D space.

It is also interesting to discuss the relationship between our approach and the scene flow estimation since the scene flow can be seen as the extension of 2D optical flow to 3D. In our approach, we search for the velocity vector centered on a local 3D patch that transports the particles within the patch and produces the most accurate particle distribution for the next frame. The accuracy is thus examined on 2D image space. This can also be interpreted as we modulate the 3D velocity vector so that its projection on 2D image best fits the corresponding image flow vector in the sense of optical flow. Note that we can not compare the 3D velocity vector directly to the 2D image velocity vector calculated by an optical flow algorithm because the 2D optical flow in corresponding patch is a superposition of the whole 3D vector along the line of sight.

Summary of algorithm

Here we briefly state the algorithm of estimating velocity fields using an ensemble:

1. We need to do triangulation or IPR to obtain initial particles' positions Xk-1,p . The initial velocity field u init is also needed and can be set to zeros.

2. First we need to decompose the domain into small local patches. For different particle densities, we make sure that each local patch contains 5 to 10 particles in average. This give us directly the first reconstruction resolution depending on the particle density.

3. We perturb the velocity u using some errors statistics,

u j = u init + ξ j , yielding the ensemble E u ∈ R n×N .
4. We start iterating from m = 0: 5. For the lth patch, we need to determine the q particles located in the patch at time index k -1.

6. Still in lth patch, we try to minimize the linearized cost function (5).

Assuming δu

l m = A u γ = 1 √ N-1 (E u -Ēu )γ
, a linear combination of the column space of ensemble anomaly matrix A u .

Denoting

I i, j sample,k,p = 1 √ N-1 ∑ q p=0 I i part,k {T k,p [ Xk-1,p , E j u ]} -I i part,k {T k,p [ Xk-1,p , Ēu ]}
, where j is the sample number index. This term can be computed efficiently because for each particle, it has N possible locations at time k. The term I i part,k {T k,p [ Xk-1,p , Ēu ]} is common to all. 9. Finally, we have the following linear system to solve:

∑ i q ∑ p=0 (I i sample,k,p ) T I i sample,k,p γ = ∑ i q ∑ p=0 (I i sample,k,p ) T I i res,k
where I i sample,k,p = [I i,0 sample,k,p , . . . , I i, j sample,k,p , . . . , I i,N sample,k,p ] has one sample j per column.

Recover ûl

m = A u γ.
11. Solve (4) for u G and set ûl m = ûG . 12. Do reduced IPR on image residual at frame k -1 to produce more particles.

13. Update I i res,k using ûl m and iterate back to step 5 with m = m + 1 and u l m+1 = ûl m .

To summarize, our Lagrangian PIV approach is a hybrid approach that reconstructs the 3D Eulerian velocity field immediately through tracking particles.

Results of Re3900 wake flow behind a cylinder

For the synthetic particle image setup, we followed the same procedures as in [START_REF] Yang | An ensemble filter estimation scheme for lagrangian trajectory reconstruction[END_REF]. The synthetic particle images were constructed by projecting virtual particles onto 4 virtual cameras. Those particles were transported by velocity fields obtained with Large Eddy Simulation [START_REF] Jeon | 4D flow field reconstruction from particle tracks by VIC+ with additional constraints and multigrid approximation[END_REF] characterizing the turbulent wake-flow behind a cylinder at Reynolds number equal to 3900. We used a uniform Gaussian form to model the OTF that led to an averaged particle diameter of 2 px.

We chose one snapshot as the source of velocity field and sampled two datasets at different locations with different size of domain and different size of images. We started with a small test case with a volume size of 1.5D × 1.5D × D/3, where D is the cylinder diameter. The image size was 256 × 256 pixels for this case. Then we proceed to a full scale case extending over a volume of size 6D × 1.9D × D with image sizes of 1280 × 800 pixels. For both cases, the volume was located such that the inlet face was 3D behind the center of the cylinder, i.e. in the transitional region of the wake flow where large three-dimensional motions and high shear take place [START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF]. The simulated ppp was at the level of 0.01, and the cameras were pre-calibrated accordingly. The time separation was set to a level so that the mean particle displacement between two consecutive snapshots was around 1 pixel. As discussed before, a prior guess on the velocity field was needed to initialize our approach. Here we initialized the velocity field as zero field to assess the robustness of our method. We also employed a window over-lap of 75%. The regularization coefficient α and β are set as 1 for all cases.

In the following, we assessed the quality of the reconstructed flow quantities by comparing it to the reference. Note that the source quantity was first interpolated to the reconstructed grid to provide a suitable ground truth.

A small test case

Figure 1 shows the reconstructed streamwise velocity field, the reference streamwise velocity field as well as the error distribution in the center plane Z = D/6. The reconstructed velocity field was in a good correspondence with the LES reference. We want to emphasize that although initialized by zero velocities, our method was able to converge. The final average analysis velocity error was only 2% of the initial error. In figure 2, we also plot the reconstructed normal vorticity component, the reference normal vorticity component as well as the error distribution in center plane Z = D/6. Based on both the streamwise velocity error plot and the vorticity error plot, we observed that the locations that had large vorticity were likely to have relatively large errors. This was expected because, at locations with large vorticity, our assumption of the homogeneity on the velocity in local patch was not valid anymore. A direct solution is to add more particles to the flow so that we can work on finer local path where the homogeneity assumption can be justified. However, high ppp introduces much more ambiguities for triangulation which in turn deteriorate the accuracy of the estimated velocity field as well. There is no clear solution to this problem, yet we argue that one possible solution is to consider Computational Fluid Dynamic models since they can provide finer and physical-induced perturbations to our ensemble-based data assimilation scheme.

A full scale test case

In this case the spatial resolution of the reference velocity was largely increased (nearly 4 times finer in x direction and 3 times finer in z direction). As a result, it is more difficult to obtain accurate estimation on finer resolution bounded by the particle density that stays the same (ppp=0.1) as in the previous small case. This is a direct consequence of our algorithm that requires each local patch must contain a minimum number of particles.

We compared the performance of our method against VIC+ and TomoPIV implemented in Davis10. Note that this version of VIC+ in Davis10 is labeled as VIC#, with additional constraints and multigrid approximations according to [START_REF] Jeon | 4D flow field reconstruction from particle tracks by VIC+ with additional constraints and multigrid approximation[END_REF]. To provide a fair comparison, we fed both the synthetic data sets as well as the mapping function to Davis 10. Note that we provided time-resolved sequential data for VIC+ method since it relies on the particle field reconstructed by the STB method. We also chose the reconstructed velocity field at the converged phase of STB when almost all particles were tracked successfully.

Figure 3 shows the results of isosurfaces of velocity magnitude obtained by our method (3b) compared to the reference (3a) as well as VIC+ (3c) and TomoPIV (3d). Among the three methods, TomoPIV diverged the most from the reference and it failed to provide a good estimation at regions (3D < x < 4D, 0 < y < 1D) where the velocity gradient is high. Both VIC+ and our method yield results that were in good accordance with the reference. Our method provided more local coherent structures compared to VIC+, which as a global method, tended to smooth out those small scale variations.

The isosurfaces of low and high vorticity magnitude of different approaches are visualized in figure 4. Both VIC+ and our method could recover the regions with strong vorticity. VIC+ filtered out the high frequency variations that were nevertheless maintained with our method. TomoPIV failed to reconstruct the small scale vortex with high vorticity. At regions with low vorticity, only our method is able to provide a good estimation. Both VIC+ and TomoPIV tended to overestimate the vorticity in regions where the rotation was trivial. This reflected that these two methods were not sensitive to small scale variations. Our method provided a visually more consistent vorticity field compared with the reference. This is because our method, acting as a local-based method, captured small scale variations compared to VIC+. At the meantime, our method, concentrating on a particle-based image matching scheme, avoided the averaging effect introduced in TomoPIV.

Conclusion

In this work, we proposed a new method to infer the Eulerian velocity directly from particle images. Combined with IPR, the proposed method was able to reconstruct the underlying velocity field through a single optimization procedure. The proposed scheme was assessed with synthetic images. It was proved that our approach was quite accurate and in addition robust because it could converge although initialized with zero velocities. We are currently working on real datasets in order to demonstrate the potential of the proposed technique. A viable perspective is to incorporate CFD models under current ensemble framework. It is expected that the CFD model will bring finer flow resolution and physical balanced/coherent perturbations into the ensemble approached.

Figure 1 :

 1 Figure 1: The streamwise velocity field of center plane at Z = D/6: ground truth (Left), reconstructed (Middle) and their absolute difference (Right).

Figure 2 :

 2 Figure 2: The vorticity field of center plane at Z = D/6: ground truth (Left), reconstructed (Right) and their absolute difference (Right).

Figure 3 :

 3 Figure 3: Isosurfaces of velocity magnitude |u|/U c = 0.6 (silver) and |u|/U c = 1.1 (red) of reference field (a) and reconstructed field using LaPIV (b), VIC+ (c) and TomoPIV (d) in Davis10.

Figure 4 :

 4 Figure 4: Isosurfaces of vorticity magnitude |ω| = 0.05U c /D (silver) and |ω| = 0.5U c /D (red) of reference field (a) and reconstructed field using LaPIV (b), VIC+ (c),and TomoPIV (d) in Davis10.