Multi-Δt 3D-PTV based on Reynolds decomposition

A novel approach is investigated to extend the range of measurable velocities by 3D-PTV systems. The method is specifically conceived for robotic volumetric PTV measurements, but it has applications for other similar techniques. The multi-Δt method relies upon combining the information from two or more sets of double-frame images with pulse separation of different time duration. Measurements with a short time separation yield a robust particle velocity field estimation with a higher percentage of valid vectors, yet a low measurement precision. Conversely, measurements with longer time separation potentially offer a higher measurement precision but suffer from an increased probability of spurious particle pairing. Reynolds decomposition is used to combine the two (dual-Δt) sets where a predictor for the mean particle displacement and its statistical dispersion is used to pair particle recordings from a longer time separation. For this reason, this method is aimed at the analysis of turbulent flows where the Reynolds decomposition is meaningful (e.g. turbulent flows with steady/quasi-steady boundary conditions). The extent of the search region is selected dynamically, based on the estimate of the velocity fluctuations from the short time separation evaluation. A more advanced variant of the algorithm contemplates the progressive increase of the pulse separation (multi-Δt) until the expected dispersion of data due to turbulent fluctuations eventually exceeds the distance between neighbouring particles. Flow measurements in the near wake of a truncated cylinder obstacle and of an Ahmed body are carried out to examine the performance of the proposed method. Reference data is taken from time-resolved multi-frame analysis based on the Shake-The-Box (STB) algorithm. The two experiments differ for the measurement principle used: the first one is conducted with a tomographic-like system (large aperture), whereas the latter uses coaxial volumetric velocimetry. The rate of correct pairing as well as the velocity dynamic range dependence upon the choice of the time separation are monitored and discussed. The results compare favourably with the STB analysis, indicating that the measurement of the time-average velocity field can be based on dual-Δt 3D-PTV measurements removing the constraint of time-resolved particle motion recording.


Introduction
The introduction of coaxial volumetric velocimetry (CVV, Schneiders et al 2018), in combination with the use of heliumfilled soap bubbles as flow tracers for large-scale measurements (Bosbach et al 2008), has reduced some requirements of system calibration and optical access for three-dimensional (3D) flow measurements. The coaxial arrangement between the camera's lines of sight (at low tomographic aperture) and the illumination system enable single direction optical access, as opposed to tomographic PIV systems. When combined with robotic manipulation, CVV becomes suited to automated measurements over complex objects, as demonstrated by Jux et al (2018) covering a total measurement volume of about 2 m 3 by means of approximately 400 views around a full-scale cyclist. CVV is based on the time resolved analysis of highspeed recordings of particle tracers motion using the Shake-The-Box algorithm (STB, Schanz et al 2016), an efficient particle tracking method for multi-frame recordings. However, the latter requires high-speed PIV equipment. In the case of CVV measurements, currently available compact CMOS cameras do not exceed 1000 Hz, resulting in a maximum flow speed for measurements up to approximately 10 m s −1 .
Experiments at higher flow velocity are hampered by the above limitations unless based on dual-frame image recording (e.g. by frame-straddling), where image separation ∆t below the microsecond. It is known, however, that the dynamic velocity range (DVR, Adrian 1997) of time-resolved (TR) measurements can be higher than that obtainable with dualframe recordings (Lynch and Scarano 2013). The DVR issue is exacerbated for CVV where the in-depth velocity component is about 10 times less accurate than the other two components due to the low tomographic aperture .
A method is investigated here to perform 3D-PTV analysis in double-frame mode and restore a DVR comparable to that achieved with TR multi-frame techniques. The approaches reported in literature aiming at increasing the DVR of PIV either increase the maximum particle image displacement or decrease the minimum resolvable displacement. To enlarge the maximum resolvable displacement, Fincham and Delerce (2000) developed a multi-∆t approach on three-frame recordings separated by ∆t and 2∆t. Cross-correlation at separation ∆t produces a predictor for the analysis at time separation 2∆t. Multi-∆t acquisitions have been also used to quantify uncertainties (Nogueira et al 2009 andNogueira et al 2011 for the peak locking). Multi-∆t recordings analysis has been exploited by Scharnowski et al (2019) to quantify flow turbulence intensity form PIV measurements.
With the aim of increasing the DVR, Pereira et al (2004), Hain and Kähler (2007) and Persoons and O'Donovan (2011) have developed several multi-frame approaches for TR recordings, where the time separation is locally optimized based on the flow conditions or on the cross-correlation signal to noise ratio. Conversely, as far as the reduction of the minimum resolvable velocity is concerned, techniques of correlation averaging have proven to be effective. The pyramid correlation (Sciacchitano et al 2012) further expands the method through linear combination of correlation maps obtained at different time separation. Non-linear motions were taken into account by Lynch and Scarano (2014) and later by Jeon et al (2014) with a least-squares approach. For particle tracking velocimetry (PTV), Cierpka et al (2013) showed that the use of four or more time steps in combination with a multi-∆t image analysis greatly enhances a reliable particle pairing even with high levels of the seeding concentration. From the above discussion, it emerges that multi-frame approaches (recordings that encompass more than two snapshots) have been most pursued to increase the DVR of PIV and PTV techniques.
In the present work, we investigate the use of multi-step analysis of double-frame recordings making use of a variable time separation between exposures. The work focuses on the potential to increase the DVR of measurements and a specific discussion is made for low-aperture 3D-PTV systems like CVV and astigmatism PTV (Cierpka et al 2010).

Two-frame particle tracking principles
Particle tracking principles are amply discussed in the literature (Malik andDracos 1993, Pereira et al 2006;amongst others). Here, fundamental definitions and properties are recalled for use in the discussion presented further. Let us consider particle tracers distributed in the physical space of coordinates X, Y, Z. When at uniform concentration C, the average distance λ between neighbouring particles, following Pereira et al (2006), reads as . (1) The nearest neighbour (NN) principle is arguably the simplest approach to pair subsequent images of a particle tracer. Considering a particle displacement ∆X occurring between two subsequent frames with time separation ∆t, the ratio γ = |∆X|/λ between the displacement and the mean particle distance determines the probability of obtaining a correct pairing between the two images of the same particle. A schematic illustration is given in figure 1, where the condition of γ < 1 (left) yields a high probability of successful pairing. Conversely, when γ > 1 (right), the increased search region leads to false pairing when the NN principle is applied.
The NN algorithm is usually coupled with a condition of maximum search distance (Pereira et al 2006), here referred to as search radius R S . The above discussed condition for a high probability of correct detection translates into a relationship between the search radius and the average particle distance, more specifically: Several criteria to optimize the choice of R S are given in the literature. Malik and Dracos (1993) proposed the following: From the above, it can be concluded that choosing a particle displacement significantly smaller than the particle distance Particle images from double-frame recordings. Left: Length travelled by particles between first and second exposure is smaller than the distance λ separating neighbouring particles (γ < 1). Right: Particle displacement (green arrows) exceeds the inter-particle distance (γ > 1) and the search radius R S includes more than one candidate for pairing. is a favourable condition to correct particle pairing. However, the accuracy of the instantaneous velocity measurement directly depends upon the length of the particle displacement. The DVR, defined as the ratio between the maximum and minimum resolvable velocity (viz. displacement) can be written as (Adrian 1997): In the intermediate term, |∆X| max represents the maximum particle displacement in the measurement domain and σ s the minimum resolvable displacement. Following Adrian (1991), the latter can be described as a function of the particle image diameter d τ and the uncertainty c τ of the particle image centroid position.
A theoretical limit for DVR can be formulated considering equations (2) and (4): Equation (5) expresses the trade-off between the DVR and the instantaneous tracers' concentration, through the parameter λ.
As an illustration of the above, reference data on the probability of correct pairing using double-frame recordings with the NN principle is obtained with a Monte-Carlo (MC) simulation. In a volume 150 × 30 × 30 mm 3 , N = 50 particles are randomly distributed, resulting in the average concentration C = 0.4 particles cm −3 . From equation (1) the mean particle distance is λ = 8.6 mm. The relative displacement parameter γ is varied in the interval [0.01-2.75] by changing the particle displacement ∆X. Particles are paired with the NN algorithm and the fraction of correct pairing η p is evaluated. Statistical results are achieved by repeating the random simulation 10 000 times.
For γ < 0.2, the correct pairing is higher than 99% (see grey curve in figure 9, left). When γ > 0.20, false pairings start appearing. If the criterion prescribed by Malik and Dracos (1993) is chosen, γ = 0.3, the ratio of correct pairing reads η p = 0.98, in agreement with the results obtained by the latter authors.
If we now assume σ s = 0.1 mm as the uncertainty of the particle displacement estimation, it follows that for a displacement of 8 mm DVR ≈ 80. However, for the given concentration of 0.4 particles cm −3 , the latter displacement corresponds to γ = 1, leading to a probability of false pairing of approximately 50%. Conversely, imposing a correct pairing probability of 98% leads to a value of γ = 0.3 and a corresponding DVR < 30.
From this discussion, considering a given particle concentration, it becomes clear that there is a fundamental limit in the trade-off between robustness and DVR. When robustness is guaranteed (figure 1, left) a lower DVR is returned, with particle displacement comparable to the particle diameter. A higher DVR (figure 1, right) comes at the cost of pairing reliability, unless the process makes use of a predictor for the particle displacement, as discussed in the next section.

Particle pairing aided with a displacement predictor
The use of multi-step analysis of particle motion has been demonstrated to effectively improve the probability of correct pairing even at high concentration of particles (Bastiaans et al 2002;among others). The super-resolution method proposed by Keane et al (1995) makes use of cross-correlation analysis to produce a predictor for the displacement of individual particles inside the interrogation window. Particle pairing is then obtained based on a NN search in the second exposure at a position given by the predictor. For low image-density recordings typical of 3D-PTV measurements, however, the crosscorrelation approach becomes unsuited due to two main reasons: (i) the particle field is often represented in the physical space by their positions and not by voxel intensities; (ii) with a large inter-particle distance, a low signal-to-noise ratio is expected.
Here, similarly to the super-resolution method, an estimator of the tracer velocity is considered based on a first level analysis that yields an estimation of particles time-average velocity and expected level of the fluctuations. This is based on the Reynolds decomposition of the local flow velocity: The time average velocity V is determined from the previous analysis performed at short time separation as detailed in the next section. Such time-average velocity is used to offset the search region by an amount corresponding to the local mean predicted displacement ∆X pred = V · ∆t. In the turbulent flow regime, the actual position of an individual tracer will not coincide with the position predicted with the time average. Let us define such discrepancy by the average positional disparity vector, which reads as where σ V is the velocity standard deviation and ∆t the pulse time separation. This way, the choice of the search radius R S needs to account for the expected fluctuations only, since the displacement due to the mean velocity is considered with the predictor. The search radius reads as follows the relation: It is proposed more specifically that R S = 2 ∆X ′ , corresponding to a confidence level of 95% when the fluctuations follow a Gaussian distribution. Since the second part of equation (2) is still valid also when a predictor is available, consequently, the use of a mean velocity predictor turns the restriction posed in equation (2) into As a result, for a given velocity field, the value of the time separation ∆t can be increased of a factor V max /|σ V | when a predictor for the mean displacement is available. Due to the velocity prediction, the occurrence of correct pairings is no longer directly related to the particle displacement, but it is proportional to the ratio between the radius of search R S and the mean particle distance λ. The dynamic range of a velocity measurement making use of a predictor reads, therefore, as It can be concluded that the DVR making use of a predictor is extended with respect to the case of a single-step particle tracking according to As an illustration, if the method is used to measure a turbulent flow with fluctuations of the order of 10% of the maximum velocity, equation (11) indicates a potential order of magnitude increase of the velocity dynamic range. We should keep in mind, however, that the above analysis relies on a number of hypotheses: (1) the increase of time separation shall remain limited to the range where truncation errors are negligible with respect to random errors (Boillot and Prasad 1996); (2) the operations that determine the mean velocity predictor (binning process discussed in chapter 3) are performed at a sufficient spatial resolution and with statistical convergence of the velocity and its fluctuations to reliably apply Reynolds decomposition.
As shown by Hain and Kähler (2007), the truncation error appears in the presence of acceleration in the flow. The truncation error scales with the square of the pulse separation time ∆t when the velocity is evaluated with a central-difference scheme. This appears to be problematic along curved streamlines (radial acceleration) and when the flow rapidly decelerates or accelerates (tangential acceleration). Let us consider the former case within the core of a steady vortex, where the flow rotates like a rigid body. A relative error on the velocity magnitude due to truncation lower than 10% corresponds to a pulse separation time of 1/4 T, where T is the core turnover time. The flow vorticity ω is often monitored with PIV measurements (T = 4π/ω); imposing a time separation one order of magnitude smaller than the reciprocal of the local flow vorticity can be seen as a conservative criterion to prevent that truncation errors affect the measurement accuracy.

Multi-step algorithm
The approach presented here relies on the acquisition of two or more sets of double-frame images with varying (increasing) pulse separation time. In this section, the case in which two sets are acquired is considered, with respectively, pulse separation time ∆t 0 and ∆t 2 , with ∆t 0 < ∆t 2 . The analysis of the data at ∆t 0 features high robustness but low precision and is used for a first estimate of time-averaged and fluctuating velocity with criteria defined in equations (1)-(5). These estimates are then used to aid the analysis at separation ∆t 2 .
The analysis of the dataset ∆t 0 is described first: the 3D particle detection is based on the iterative particle detection algorithm (IPR, Wieneke 2013). Particle pairs are determined selecting the closest (in 3D space) particle between the two frames (NN approach, Pereira et al 2006). The search radius R S here needs to account for the maximum expected particle displacement (equation 2).
The result of this evaluation yields the instantaneous flow velocity. For each particle, the velocity vector is placed at the midpoint between the two positions. The time-averaged velocity field is reconstructed with the binning procedure as described by Agüera et al (2016): (1) all the velocity vectors pertaining to the series of recordings are combined into a single ensemble, which increases the spatial density of the velocity information; (2) the measurement volume is divided in sub-volumes (bins) with dimension L bin ∼ O 10 −2 m arranged on a Cartesian grid. Similar to what is done in PIV image processing, overlap between adjacent bins (e.g. by 75%) decreases the distance between neighbouring vectors. The data captured inside a single bin features a cloud of  velocity samples as a result of local turbulence and the measurement uncertainties (figure 2, left). Performing the Reynolds decomposition (equation 6), we obtain the average displacement ∆X 0 and its fluctuations ∆X ′ 0 . In order to decrease the error due to unresolved velocity gradient, the velocity samples are weighted according to their distance with respect to the centroid of the bin. A Gaussian weighting function is applied, following the approach proposed by Agüera et al (2016). The Gaussian is then centred in the bin center and has a standard deviation equal to half of the bin size.
In the second stage, recordings acquired at a time separation ∆t 2 > ∆t 0 are interrogated, making use of the above results. Also in this case, particle detection is performed using IPR, resulting in a cloud of particles for both the exposures of all the recordings. Considering the particles triangulated in the first exposure, the time-averaged velocity field measured in the previous stage is interpolated at particle positions and the predicted displacement is calculated trough a linear scaling (homothety) and reads as Similarly, the choice of the search radius is locally determined based on the estimated level of velocity fluctuations: In synthesis, both average displacement and the radius of search are obtained through homothety with the coefficient given by the ratio of time separation, as shown in figure 3.
Then, the NN criterion is applied between the predicted arrival position and the particles detected within the spherical After pairing, the binning procedure yields again data on a Cartesian grid.
The logics of the entire algorithm are illustrated in figure 4 and comprise the following operations.
• Acquisition of multiple double-frame datasets with increasing time separation; • 3D particle detection by IPR at ∆t 0 ; • particle pairing with NN principle; • ensemble-average of sparse velocity vectors within bins; • displacement predictor and search radius are built upon the Reynolds decomposition of the velocity inside the bin; • step 2 is repeated at ∆t 2 ; • the predictor for displacement and fluctuations is applied at the position of the tracers; • step 3 is repeated based on the predicted position and the search radius; • step 4 yields data on a Cartesian grid, with time separation ∆t 2 .

Chain-variant of the multi-step algorithm
Inferring a displacement predictor and search radius from measurements at the shortest pulse separation time requires a robust and accurate estimate. The value of R S may, however, be affected by a high relative uncertainty. When amplified by the homothety the uncertainty may lead to overestimating the value of R S , in turn increasing the false detection probability at time separation ∆t 2 . This effect is mitigated if one or multiple additional steps are included between ∆t 0 and ∆t 2 , as shown in figure 5. As presented in the previous section, R S is built from ∆X ′ , more specifically from the standard deviation σ V , that can be decomposed as (Sciacchitano and Wieneke 2016): While the first term under the square root, representing the physical flow fluctuations (∆X' fluc ), scales linearly with ∆t, the fluctuations associated with measurement noise (∆X' err ) can be considered independent of the particle displacement. Given the above, any measurement with time separation larger than ∆t 0 yields an estimation of the velocity fluctuations that is less affected by σ V,err . By this method, the choice of the search radius R S for the final and largest time separation becomes significantly less affected by noise, reducing the search area and the probability of erroneous pairing.

Application to turbulent wake flows
Two experiments dealing with the wake of bluff objects have been conducted at the Aerodynamics Laboratory of the Aerospace Engineering Department of TU Delft. In the first experiment a large aperture tomographic setup is used, whereas the second experiment makes use of robotic volumetric PTV (called robotic volumetric PIV in the original work from Jux et al 2018), which makes use of a CVV system .

Near-wake of truncated cylinder
The turbulent flow developing behind a truncated cylinder interacting with a flat plate turbulent boundary layer was described in the study of Schneiders et al (2016). The experiments were performed in a low-speed wind tunnel with a crosssection of 0.4 × 0.4 m 2 at free-stream velocity of 5 m s −1 . A flat plate produces a turbulent boundary layer of approximately 25 mm thickness 1 m downstream of its leading edge. A truncated cylinder of 100 mm diameter and height is placed in the symmetry plane of the plate. The Reynolds number based on the cylinder diameter is Re D = 3.5 × 10 4 . The measurement volume was 30 × 15 × 20 cm 3 and is schematically represented in figure 6 left. The use of helium-filled soap bubbles as flow tracers was necessary to produce sufficient light scattering over such volume.
The data consist of three sequences of 2000 frames acquired at 2 kHz (∆t 0 = 0.5 ms) with four high-speed CMOS cameras (Photron FastCAM SA1, 1024 × 1024 pixels) subtending a solid angle of approximately 40 × 40 square degrees.
The particle image recordings were evaluated with the algorithm STB (Schanz et al 2016). The velocity field obtained with STB is considered as reference to evaluate the performances of the multi-step analysis based on doubleframe recordings.
The STB analysis makes use of long tracks produced by a particle imaged at multiple time instants. In this case, a particle has been considered valid if it was tracked for at least six time steps. Given the track length and the least-squares fit used to model the trajectory, a high measurement precision is achieved (estimated in the order of 10 -3 m s −1 ), as discussed in Schanz et al (2016).
To apply the proposed method, a double-frame dataset has been constructed from the original multi-frame recordings. Particles detected at the time step t by STB have been considered for the first frame, assigning to the second frame the particles found at the time step t + α∆t 0 , with α = [1, 2, 3]. The algorithm STB assigns a unique track ID number to each particle tracked across the domain. This number permits us to evaluate if the particle pairing performed by the proposed method is correct.
The results obtained by STB and with the dual-frame analysis are subject to the same binning process to yield the velocity distribution on a Cartesian grid, as described in section 3. Volumes 2 × 2 × 2 cm 3 have been considered, with an overlap factor of 75%. The final vector pitch is then 0.5 mm.
The flow field around the cylinder exhibits large vortices and separated regions (figure 7, left), making it well suited to analyse the accuracy and robustness of the tracking algorithm under varying flow properties. The near-wake is characterized by regions of high fluctuations due to the interaction between the shear layers created at the sides, shown in figure 7 (right), and the two counter-rotating vortices that originate from the top of the object. double-frame image analysis yields a similar value when the multi-step algorithm is used. The single-step analysis at time separation ∆t 0 exhibits a significant velocity bias (mean velocity of 5.5 m s −1 ) and a dispersion one order of magnitude larger than the STB measurement (1.49 m s −1 and 0.16 m s −1 , respectively). The single-step analysis at larger time separation is directly compromised by a large number of spurious pairs (71.5%), leading to very large bias and random errors.
Two kinds of analysis are conducted based on the Reynolds averaged predictor. In both analyses, the velocity predictor is built from the time-average velocity and the velocity . Left: Correct pairing probability ηp versus the ratio γ between the particle displacement and the mean particle distance. Results obtained with single-step analysis and comparison to MC simulation of free-stream conditions. Right: ηp variation with the ratio between search radius and mean particle distance (legend same as for left figure). Data obtained by multi-step analysis with ∆t 2 = 3 • ∆t 0.
fluctuations evaluated with the ∆t 0 single-step analysis. In the first case, indicated with ∆t 0 multi-step in table 1, the time separation ∆t 2 = ∆t 0 , whereas in the second analysis, indicated with 3 · ∆t 0 multi-step, ∆t 2 = 3 · ∆t 0 . Both multi-step methods yield a major reduction of the number of spurious pairs, leading to a probability of correct pairing of 99%. Furthermore, the mean velocity and the velocity fluctuations evaluated with these analyses agrees well with the reference data. Nevertheless, the ∆t 0 multi-step analysis exhibits slightly larger fluctuations than the reference and 3 · ∆t 0 multi-step analysis due to the larger relative uncertainty of the measured displacement, thus confirming the enhanced precision achievable with a larger time separation. Figure 9 (left) illustrates the correct pairing probability η p versus the relative displacement γ evaluated in the regions A, B and C shown in figure 7. The MC simulation of the free-stream flow is taken as the reference behaviour. Although reproducing a similar trend, the correct pairing probability η p obtained with MC simulation slightly overpredicts the results obtained by the single-step approach in the free-stream domain. The most evident behaviour observed by this analysis is that the regions with an increased level of turbulent fluctuations exhibit a more rapid drop of correct pairing probability when the single-step time separation is increased.
The introduction of the displacement predictor based on Reynolds average increases the percentage of correct pairing in all the considered region of the flow: in the free-stream (region A) ∆t can be extended up to 8 times with the probability of correct pairing remaining above 98%. In the turbulent regions, such as the recirculation region in the wake (region B), and the lateral shear layer (region C), the use of the predictor yields benefits up to ∆t 2 = 3∆t 0 . The latter behaviour is consistent with equation (11), given the higher level of velocity fluctuations in the object wake. Figure 9 (right) shows the probability of correct pairing ratio η p with respect to the ratio R S /λ for the multi-step analyses. The curves collapse approximately onto the same behaviour, indicating a universal relation between R S /λ and the probability of successful pairing. Considering η p = 0.95 as acceptance criterion, the corresponding search radius becomes approximately R S = 0.5λ. The latter may be proposed as design criterion to choose the upper limit of ∆t 2 for a given experiment comprising recordings at different values of the time separation.
Given the spatial variability of the flow properties a single optimum value for ∆t 2 cannot be identified. An illustration figure 10 shows the spatial distribution of η p at the plane Y = 0 mm, while increasing ∆t. The single-step analysis rapidly degrades in regions of large displacement. For instance, when ∆t = 3•∆t 0 , only the low-velocity region in the near wake exhibits a high percentage of correct pairing (top-right). The use of the predictor for the case ∆t 2 = 3•∆t 0 (bottom-left) leads to η p > 0.9 in most of the flow field, except for the wake with high fluctuations (see figure 10 top-right), where η p > 0.7.
Further extension of the pulse separation time (figure 10, bottom-right) results in frequent false pairing, mostly in the wake, which is due to the condition given by equation (2) not being respected.
The results shown until now are obtained by applying the pairing strategy on particles previously detected by the STB algorithm, which is also considered unaffected by the phenomenon of ghost-particles (Schanz et al 2016). Therefore, by knowledge of the particle tracks, it has been possible to distinguish correct and incorrect particle pairings for the doubleframe analysis.
A more realistic situation has been simulated using the Iterative Particle Reconstruction algorithm proposed by Wieneke (2013), thus following the steps illustrated in the flow chart of figure 4.
The standard deviation of the streamwise velocity component at X = 100 mm, illustrated in figure 11 (bottom-left), shows that the single-step algorithms yield spurious fluctuations due to a significant percentage of incorrect pairings. The level of fluctuations is clearly unacceptable for ∆t = 3 • ∆t 0 , where the velocity standard deviation in the outer region exceeds 50% of the free-stream value.
Conversely, the adoption of the multi-step ∆t methodology suppresses the spurious velocity fluctuations due to incorrect pairing, yielding measured fluctuations of the same order as the reference ones obtained by STB. It must be noticed that in the regions of highest flow fluctuations, namely the two free shear layers at the sides of the models Y = ± 60 mm, the difference between the multi-step analysis and the reference data is the largest. As indicated by equation (11), the increase of the local turbulence intensity reduces the maximum achievable extension of ∆t.
The effects of the chain-variant method have been assessed considering the same final pulse separation time ∆t = 3 • ∆t 0 . The chain algorithm permits us to decrease R S in most of the flow field, with the maximum decrease that occurs in the free-stream, where R S = 1.7 mm for the dual-step algorithm and R S = 0.8 mm for its the chain-variant. For what concerns the standard deviation of the velocity, due to the limited time increase, no substantial differences is noticed.

Coaxial velocimetry in the near-wake of Ahmed Body
Experiments are performed in the Open Jet Facility (OJF) of TU Delft Aerospace Engineering Laboratories. The near-wake of the Ahmed body (Ahmed et al 1984) at a free-stream velocity of 12 m s −1 and a turbulence intensity of 0.5% (Lignarolo et al 2015) is investigated by robotic volumetric PTV . The Reynolds number based on the height H of the model is Re H = 115 000 and the selected slant angle is 25 • . Table 2 summarises the characteristics of the robotic volumetric PTV system.
The considered volume is 200 × 200 × 450 mm 3 , obtained from a single view of the CVV system (figure 12). Both TR and double-frame acquisitions have been performed. For the former, the acquisition frequency is f acq = 700 Hz. In doubleframe mode, sets of image-pairs are acquired at a rate of 340 Hz. The TR dataset is analysed with the STB algorithm from LaVision DaVis 8.1 software.
Multiple datasets with ∆t 0 = 61 µs and larger separation ∆t = [2,4,6,8,10] • ∆t 0 were acquired. The minimum pulse separation time is selected for a conservative value of γ = 0.07, guaranteeing a high probability of correct pairing at ∆t 0 . For the binning process, 2 × 2 × 2 cm 3 volumes have been considered, with an overlap factor of 75%. The final vector pitch is then 0.5 mm. The wake of the Ahmed body features the so-called Cpillar vortices: a set of counter-rotating large-scale streamwise vortices emanating from the upstream edge of the slant. As they develop downstream, the C-pillar vortices interact with the recirculation region at the back of the object, creating a complex 3D flow field. The organization of the velocity field is illustrated at X obj = 0.5H in figure 13. The presence of the two vortical structures is confirmed by the vectors in the The robotic system is characterized by the ability of measuring multiple portion of the flow that can be stitched together to obtain the final global average velocity field. For this reason, two relevant coordinate systems are mentioned: the global (object) and the intrinsic (CVV) one, shown schematically in figure 12. The following analysis is performed using the CVV coordinate system in order to analyse the properties of in-plane and the coaxial velocity components separately.
A theoretical estimate of the DVR is given, based on the characteristics of the coaxial velocimeter. Considering a c τ = 0.2, the resulting error on particle position is ε = 0.13 mm along the x-and y-directions, which becomes ε = 2.2 mm in the depth (viz. coaxial) direction . The above translate in terms of velocity uncertainty relative to free-stream value with ε u = ε v = 0.26 for a single-step double-frame measurement with ∆t = ∆t 0 = 61 µs. For a TR measurement with particle tracks comprising five samples of the particle position, the uncertainty reduces to ε u = 0.0025 (0.25%). From the above, a DVR of 4 and 400 can be inferred for the double-frame single-step and STB analysis (with five recordings separated by 1.43 ms each) respectively. This large difference is the result of two factors: the STB analysis encompasses a significantly longer time separation (approximately 23 times larger than ∆t 0 for double-frame); second, the velocity measurement is the result of a least-squares polynomial fit that reduces random errors. The multi-step analysis based on Reynolds average predictor allows increase the time separation, therefore increase of DVR, but only from the former of the two factors. In the present case the recordings with the longest pulse separation time, ∆t = 10·∆t 0 , potentially lead to a DVR = 40. However, along the coaxial direction, the DVR remains fairly limited (DVR ∼ 3-5), considering the small angular aperture of the coaxial velocimeter.
The probability density of x-and coaxial-component in the outer flow region are shown in figure 14. Synthesis of the results in terms of mean and standard deviation are presented in table 3 for the different methods. The reference is assumed to be the TR analysis from STB, which also exhibits the lowest dispersion of the velocity data.
The single-step analysis of double-frame recordings at time separation ∆t 0 = 61 µs exhibits the widest dispersion of the data, with σ u being approximately four times larger than that given by STB. The multi-step analysis progressively reduces the data dispersion by increasing ∆t 2 . A standard deviation 14% higher than the reference is obtained when ∆t 2 = 10 • ∆t 0 .
Along the coaxial direction, a much wider data dispersion is observed and the single-step analysis with ∆t = ∆t 0 returns almost a flat distribution. Increasing the time separation by the multi-step analysis, although the overall uncertainty remains large: at ∆t 2 = 10·∆t 0 , the coaxial velocity component is underestimated by approximately 30%.
The analysis until now has been carried out in the outer region, where the maximum displacement is expected. The high particle displacement corresponds to high values of γ, leading to an increase of false pairing appearance.
The amplitude of velocity fluctuations plays a crucial role for determining the success rate of correct pairing, even when a predictor velocity is available. For this reason, a region of strong spatial and temporal and fluctuations has been considered in the shear region of the C-pillar vortices. The probability distribution of u and w is analysed and shown in figure  15. The smaller displacement and the relatively low particle    concentration yield γ = 0.04 for ∆t = ∆t 0 . In this region, given the wider dispersion of the value due to the physical fluctuations exhibited by the flow, the in-plane and coaxial component show similar behaviour and the results are more closely comparable to those obtained with STB. A tenfold increase of the pulse separation time returns a velocity distribution not affected by false pairing (R S /λ = 0.4 for ∆t 2 = 10∆t 0 ). A final analysis is made to investigate the measurements of the 3D vorticity field, often inspected to understand the topology of vortices emanating from complex bluff bodies. Figure 16 shows the 3D distribution of the time-average streamwise vorticityω x by two iso-surfaces selected at ±250 Hz. The comparison is made between STB, singlestep (∆t = ∆t 0 ) and two multi-step analyses: ∆t 2 = ∆t 0 and ∆t 2 = 10 • ∆t 0 , respectively. The C-pillars vortices visualisation using the single-step analysis suffers from random fluctuations appearing in the entire measurement domain.
These fluctuations are mostly associated with the large uncertainty on the coaxial component (and its spatial derivative) that takes part in the formulation of the streamwise vorticity. The multi-step analysis at shortest time separation exhibits some noise reduction, ascribed to the reduction of incorrect pairings when a displacement predictor and a smaller search radius are used. When the pulse separation is extended, with ∆t 2 = 10·∆t 0 , noisy fluctuations are considerably attenuated and a more regular vorticity iso-surface is obtained, in better agreement with the STB analysis.

Conclusions
A novel method for the analysis of 3D-PTV experiments based on double-frame recordings has been proposed, which is based on Reynolds decomposition. The method yields the time-average velocity field from the analysis of the instantaneous particle velocity (compared to other multi-frame methods proposed in literature, like from Hassan andCanaan 1991, Schanz et al 2016).
In the multi-∆t method two or more sets of recordings are necessary to produce first a robust displacement predictor, based on a short-time separation, and afterwards extend the displacement with a dataset obtained at larger time separation. A theoretical analysis shows that the DVR of the multi-step method can be significantly higher than single step analysis in flows with low to moderate turbulence. A chain-like variant of the multi-step method has the additional benefit of reducing the bias that overestimates the amplitude of turbulent fluctuations.
Two experimental databases are used to assess the multi-∆t method and compare it to TR analysis made with the STB algorithm. The multi-step analysis clearly extends the DVR of single-step analysis. The fundamental limit to extending the time separation for the multi-step method lies in the condition where the displacement dispersion caused by the turbulent velocity fluctuations become comparable to the interparticle distance. A posteriori analysis suggests R s < 0.5 λ as an experiment design criterion for the optimal extension of the time separation, which corresponds to the condition ∆t 2max = 0.25λ/|σ V |.